text:
# Convert Julia expressions to K parse trees
const OPCODE1 = Dict(
    :+     => 1,
    :-     => 2,
    :flip  => 1,
    :neg   => 2,
    :first => 3,        # *:
    :reciprocal =>  4,  # %:
    :where      =>  5,  # &:
    :reverse    =>  6,  # |:
    :null       =>  7,  # ^:
    :group      =>  8,  # =:
    :hopen      =>  9,  # <:
    :hclose     => 10,  # >:
    :string     => 11,  # $:
    #              12   #:,:
    :count      => 13,  # #:
    :floor      => 14,  # _:
    :not        => 15,  # ~:
    :key        => 16,  # !:
    :distinct   => 17,  # ?:
    :type_      => 18,  # @:
    :value      => 19,  # .:
    :read0      => 20,  # 0::
    :read1      => 21,  # 1::
    #              22   # 2::
    :avg        => 23,
    :last       => 24,
    :sum        => 25,
    :prd        => 26,
    :min        => 27,
    :max        => 28,
    :exit       => 29,
    :getenv     => 30,
    :abs        => 31,
    :sqrt       => 32,
    :log        => 33,
    :exp        => 34,
    :sin        => 35,
    :asin       => 36,
    :cos        => 37,
    :acos       => 38,
    :tan        => 39,
    :atan       => 40,
    :enlist     => 41,
    :var        => 42,
    :dev        => 43,
)
const OPCODE2 = Dict(
    :(=)=>0,      # :
    :+ => 1,
    :- => 2,
    :* => 3,
    :/ => 4,  # %
    :& => 5,
    :| => 6,
    :^ => 7,  # should this be 32 (xexp)?
    :(==) => 8,
    :< => 9,
    :> => 10,
    :$ => 11,
)

# q)asc string (key .q)where{$[100=type x;1=count(value x)1;0b]}each get .q
const FUNC1 = Set([
    :asc,
    :avgs,
    :cols,
    :desc,
    :fkeys,
    :gtime,
    :iasc,
    :idesc,
    :keys,
    :lower,
    :ltrim,
    :med,
    :meta,
    :next,
    :parse,
    :rand,
    :rank,
    :rtrim,
    :sdev,
    :show,
    :signum,
    :svar,
    :tables,
    :til,
    :trim,
    :ungroup,
    :upper,
    :view,
    :views,
])

# q)asc string (key .q)where {$[100=type x;2=count(value x)1;0b]}each get .q
const FUNC2 = Set([
    :asof,
    :cross,
    :cut,
    :each,
    :ema,
    :except,
    :fby,
    :ij,
    :ijf,
    :inter,
    :lj,
    :ljf,
    :mavg,
    :mcount,
    :mdev,
    :mmax,
    :mmin,
    :mod,
    :msum,
    :over,
    :peach,
    :pj,
    :prior,
    :rotate,
    :scan,
    :scov,
    :set,
    :sublist,
    :sv,
    :uj,
    :ujf,
    :vs,
    :xasc,
    :xbar,
    :xcol,
    :xcols,
    :xdesc,
    :xgroup,
    :xkey,
    :xlog,
    :xprev,
    :xrank,
])

# q)asc string (key .q)where {$[100=type x;3=count(value x)1;0b]}each get .q
const FUNC3 = Set([
    :aj,
    :aj0,
    :ej,
    :ssr,
])

func(::Type{Val{false}}, x::Symbol) = ks(string(".q.", x))
func(::Type{Val{true}}, x::Symbol) = k(0, string(".q.", x))

function op1(x::Symbol)
    x in FUNC1 && return func(Val{GOT_Q}, x)
    opcode = get(OPCODE1, x, -1)
    opcode == -1 ? ks(x) : ktj(101, opcode)
end

function op2(x::Symbol)
    x in FUNC2 && return func(Val{GOT_Q}, x)
    opcode = get(OPCODE2, x, -1)
    opcode == -1 ? ks(x) : ktj(102, opcode)
end

call(f) = knk(2, op(f, 1), ktj(101, 0))

function call(f, a)
    g = K_Ref(op1(f))  # guarded reference
    knk(2, r1(g.x), K_new(a))
end

function call(f, a, b)
    g = K_Ref(op2(f))
    x = K_Ref(K_new(a))
    y = K_Ref(K_new(b))
    knk(3, r1(g.x), r1(x.x), r1(y.x))
end

function call(f, a, b, c)
    g = K_Ref(op2(f))
    x = K_Ref(K_new(a))
    y = K_Ref(K_new(b))
    z = K_Ref(K_new(c))
    f in (:+, :*) ?
        knk(3, r1(g.x), r1(x.x), knk(3, r1(g.x), r1(y.x), r1(z.x))) :
        knk(4, r1(g.x), r1(x.x), r1(y.x), r1(z.x))
end

function call(f, args...)  # 4 or more args
    args = map(a->K_Ref(K_new(a)), args)
    if f in (:+, :*)
        g = K_Ref(op2(f))
        r = foldr((x,y)->K_Ref(knk(3, r1(g.x), r1(x.x), r1(y.x))), args)
        return r1(r.x)
    end
    g = f === :enlist ? K_Ref(op1(f)) : K_Ref(ks(f))
    knk(length(args) + 1, r1(g.x), map(a->r1(a.x), args)...)
end

const ASSIGN = K_Ref(GOT_Q ? k(0, ":") : ktj(102, 0))

function set(x::Symbol, y)
    v = K_Ref(K_new(y))
    knk(3, r1(ASSIGN.x), ks(x), r1(v.x))
end

function block(args)
    args = Iterators.filter(a->(a.head !== :line), args)
    args = map(a->K_Ref(K_new(a)), args)
    knk(length(args) + 1, kc(';'), map(a->r1(a.x), args)...)
end

unmangle(x::String) = x[1] == '_' ? "."*x[2:end] : x
dotted(ex) = ks(unmangle(dotted_string(ex)))
dotted_string(ex::QuoteNode) = string(ex.value)
dotted_string(ex::Symbol) = string(ex)
dotted_string(ex::Expr) = join(map(dotted_string, ex.args), ".")

function tree(ex::Expr)
    if ex.head === :call
        call(ex.args[1], ex.args[2:end]...)
    elseif ex.head === :(=)
        set(ex.args...)
    elseif ex.head === :block
        block(ex.args)
    elseif ex.head === :.
        dotted(ex)
    else
        error("expression $ex is too complex")
    end
end

K_new(ex::Expr) = tree(ex)


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

using Markdown
using InteractiveUtils

# ╔═╡ 71a70a62-28fb-11eb-2092-7d3640b6ddf5
using Pkg, DrWatson

# ╔═╡ 71d2cdaa-28fb-11eb-3826-cf35b99cde1f
begin
	#@quickactivate "StatisticalRethinkingStan"
	using Distributions
	using StatisticalRethinking
end

# ╔═╡ d31428e6-85d2-11eb-396b-c9d6da91e9ae
md" ### Clip-10-05-13s.jl"

# ╔═╡ a53d874e-8cd2-11eb-3b71-211732d6f4de
md" ### Snippet 10.5"

# ╔═╡ ed43435a-85d2-11eb-0aea-2ffa1acc2cc5
p = [
	1/4 1/4 1/4 1/4;
	2/6 1/6 1/6 2/6;
	1/6 2/6 2/6 1/6;
	1/8 4/8 2/8 1/8
]

# ╔═╡ af15068e-8cd2-11eb-1389-af62bfdfaab8
md" ### Snippet 10.6"

# ╔═╡ 5480e430-876a-11eb-0c63-aff9c5cf17d2
begin
	r, c = size(p)
	h = zeros(c)
	for i in 1:r
		s = [p[i, j] == 0 ? 0 : p[i, j] * log(p[i, j]) for j in 1:c]
		h[i] = -sum(s)
	end
	h
end

# ╔═╡ b5d131fa-8cd2-11eb-2302-456dde8b349f
md" ### Snippet 10.7"

# ╔═╡ c5dbfbae-8c12-11eb-2174-c9794a1c39d0
begin
	pr = 0.7
	a = [(1-pr)^2, pr*(1-pr), (1-pr)*pr, pr^2]
end

# ╔═╡ 99f30866-8cd2-11eb-30c8-4f54788620c2
md" ### Snippet 10.8"

# ╔═╡ 0c763214-8c13-11eb-3acb-35a54497de6c
-sum(a .* log.(a))

# ╔═╡ 87fa5676-8cd2-11eb-3a6b-b563dae2ae8d
md" ### Snippet 10.9"

# ╔═╡ 35ef080a-8c13-11eb-24ed-0d9e65e13069
function sim_p(g = 1.4)
	x123 = rand(Uniform(), 3)
	x4 = (g .* sum(x123) .- x123[2] .- x123[3]) ./ (2 - g)
	z = sum(vcat(x123, x4))
	p = vcat(x123, x4) ./ z
	(h = -sum(p .* log.(p)), p = p,)
end

# ╔═╡ 7499df52-8cd2-11eb-3490-6dc28108d33d
md" ### Snippet 10.11"

# ╔═╡ d4beaeba-8c13-11eb-2801-23af8049ca4a
begin
	h2 = [sim_p() for i in 1:100000]
	df = sort(DataFrame(NamedTuple{keys(h2[1])}.(h2), copycols=false), :h)
end

# ╔═╡ cc3bb1ea-8cd2-11eb-20f6-8176bc9783a9
md" ### Snippet 10.10"

# ╔═╡ f39459c6-8cad-11eb-242a-e3f110056f13
begin
	res = df[[91000, 12000, 700, 5], :]
	fig2 = Vector{Plots.Plot{Plots.GRBackend}}(undef, size(res,1))
	xlabels = ["ww", "wb", "bw", "bb"]
	titles = [:A, :B, :C, :D]
	xs = [xlabels[i] for i = 1:4]
	for i in 1:size(res,1)
		fig2[i] = plot(leg=false)
		plot!(xs, res[i, :p], title=titles[i])
	end
end

# ╔═╡ 13ac6c8e-8cb7-11eb-37f7-e7221773fb74
begin
	fig1 = density(df.h, xlab="Entropy", ylab="Density", leg=false)
	for i in 1:size(res, 1)
		vline!([res[i, :h]], lab=String(titles[i]))
		annotate!([(res[i, :h]-0.01, 2,
			Plots.text(titles[i], 10, :red, :right))])
	end
	fig3 = plot(fig2..., layout=(2, 2))
	plot(fig1, fig3, layout=(1, 2))
end

# ╔═╡ 4df0844e-8cd2-11eb-0c08-519ceb986561
md" ### Snippet 10.12"

# ╔═╡ 8b76cf64-8caa-11eb-10e9-3962fa03da68
findmax(df.h)

# ╔═╡ 29e01bac-8cd2-11eb-1613-11110d0bb1c8
md" ### Snippet 10.13"

# ╔═╡ f4dc21e2-8cd1-11eb-1f5f-5950f5c2c812
df[end, :p]

# ╔═╡ 22d43054-8661-11eb-0dd7-fb846edb81c4
md" ### End clip-10-05-13s.jl"

# ╔═╡ Cell order:
# ╟─d31428e6-85d2-11eb-396b-c9d6da91e9ae
# ╠═71a70a62-28fb-11eb-2092-7d3640b6ddf5
# ╠═71d2cdaa-28fb-11eb-3826-cf35b99cde1f
# ╟─a53d874e-8cd2-11eb-3b71-211732d6f4de
# ╠═ed43435a-85d2-11eb-0aea-2ffa1acc2cc5
# ╟─af15068e-8cd2-11eb-1389-af62bfdfaab8
# ╠═5480e430-876a-11eb-0c63-aff9c5cf17d2
# ╟─b5d131fa-8cd2-11eb-2302-456dde8b349f
# ╠═c5dbfbae-8c12-11eb-2174-c9794a1c39d0
# ╟─99f30866-8cd2-11eb-30c8-4f54788620c2
# ╠═0c763214-8c13-11eb-3acb-35a54497de6c
# ╟─87fa5676-8cd2-11eb-3a6b-b563dae2ae8d
# ╠═35ef080a-8c13-11eb-24ed-0d9e65e13069
# ╟─7499df52-8cd2-11eb-3490-6dc28108d33d
# ╠═d4beaeba-8c13-11eb-2801-23af8049ca4a
# ╟─cc3bb1ea-8cd2-11eb-20f6-8176bc9783a9
# ╠═f39459c6-8cad-11eb-242a-e3f110056f13
# ╠═13ac6c8e-8cb7-11eb-37f7-e7221773fb74
# ╟─4df0844e-8cd2-11eb-0c08-519ceb986561
# ╠═8b76cf64-8caa-11eb-10e9-3962fa03da68
# ╟─29e01bac-8cd2-11eb-1613-11110d0bb1c8
# ╠═f4dc21e2-8cd1-11eb-1f5f-5950f5c2c812
# ╟─22d43054-8661-11eb-0dd7-fb846edb81c4

ext: jl


text:
function read_obj(filepath::String)
    if splitext(filepath)[2] == ".obj"
        reader = vtk.vtkOBJReader()
        reader.SetFileName(filepath)
        reader.Update()

        polydata = reader.GetOutput()
        return extract_simple_data(polydata)
    else
        throw("The file input is not an obj file.")
    end
end

function read_ply(filepath::String)
    if splitext(filepath)[2] == ".ply"
        reader = vtk.vtkPLYReader()
        reader.SetFileName(filepath)
        reader.Update()

        polydata = reader.GetOutput()
        return extract_simple_data(polydata)
    else
        throw("The file input is not a ply file.")
    end
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

"""
```julia
abstract type FiniteField end
```
"""
abstract type FiniteField end

"""
```julia
abstract type AbstractCode end
```
"""
abstract type AbstractCode end

"""
```julia
NonStaticAbstractWord{N, T} = Union{NTuple{N, T}, AbstractVector{T}, AbstractString} where {N, T}
```
"""
NonStaticAbstractWord{N, T} = Union{NTuple{N, T}, AbstractVector{T}, AbstractString} where {N, T}

"""
```julia
mutable struct Word{N, T}
Word(w::NTuple{N, T})
Word(w::AbstractVector{T})
Word(w::AbstractString)
Word(i::T...)
```

A Word is an `StaticArrays.MVector` which is efficient (like tuple) yet mutable.
"""
mutable struct Word{N, T}
    w::MVector{N, T}

    Word(w::NTuple{N, T}) where {N, T} = new{N, T}(MVector{N, T}(w))
    Word(w::AbstractVector{T}) where {T} =
        (len = length(w); new{len, T}(MVector{len, T}(w)))
    Word(w::AbstractString) =
        MVector{length(w), eltype(w)}(collect(w))
    Word(i::T...) where {T} =
        (len = length(i); new{len, T}(MVector{len, T}(i)))
end

# Indexing Interface
Base.getindex(W::Word{N, T}, i::Int) where {N, T} = getindex(W.w, i)
Base.setindex!(W::Word{N, T}, v, i::Int) where {N, T} = setindex!(W.w, v, i)
Base.firstindex(W::Word{N, T}) where {N, T} = firstindex(W.w)
Base.lastindex(W::Word{N, T}) where {N, T} = lastindex(W.w)

# Abstract Array Interface
Base.size(::Word{N, T}) where {N, T} = tuple(N)
Base.length(::Word{N, T}) where {N, T} = N

"""
```julia
AbstractWord{N} = Union{NonStaticAbstractWord{N, T}, MVector{T}} where {T}
```
"""
AbstractWord{N, T} = Union{Word{N, T}, NonStaticAbstractWord{N, T}, MVector{N, T}} where {N, T}

isword(w) = w isa Word
isword(i...) = isword(i)
isabstractword(w) = w isa AbstractWord
isabstractword(i...) = isabstractword(i)

"""
```julia
Codewords{N} <: AbstractCode
```

Simply a wrapper type for a vector of abstract words of length `N`.
"""
Codewords{N} = Vector{AbstractWord{N, T}} where T
# Codewords{N} = Vector{Word{N, T}} where T

"""
```julia
struct Alphabet <: AbstractCode
```

Has the parameter Σ, which is the alphabet; a collection of strings, characters, symbols, or Ints.

---

```julia
Alphabet(Σ::AbstractArray)
```

A constructor method for the struct Alphabet.  Takes an array of letters in the alphabet, and attempts to parse them as 64-bit Ints.

---

```julia
Alphabet{N}(Σ::AbstractArray)
Alphabet{N}(Σ::AbstractString)
```

A constructor method for the struct Alphabet.  Takes in a symbols and splits it into constituent characters.  Those symbols are the letters in the alphabet.  `N` is the number of letters in the alphabet
"""
struct Alphabet{N} <: AbstractVector{Symbol} where N
    Σ::AbstractVector{Symbol}

    Alphabet(Σ::Union{Vector{T}, String}) where {T} =
        (Σ_unique = Set(Σ); new{length(Σ_unique)}(ensure_symbolic(Σ_unique)))
end # end struct

# Indexing Interface
Base.getindex(A::Alphabet, i::Int) = getindex(A.Σ, i)
Base.setindex!(A::Alphabet, v, i::Int) = setindex!(A.Σ, v, i)
Base.firstindex(A::Alphabet) = firstindex(A.Σ)
Base.lastindex(A::Alphabet) = lastindex(A.Σ)

# Abstract Array Interface
Base.size(A::Alphabet{N}) where {N} = tuple(N)
Base.length(A::Alphabet{N}) where {N} = N
Base.rand(A::Alphabet{N}) where {N} = rand(A.Σ)

"""
```julia
gensym(q::Int) -> Vector{Symbol}
```

Generates a vector of unique symbols of length q.
"""
Base.gensym(q::Int) = Symbol[gensym() for _ in 1:q]

"""
```julia
genalphabet(q::Int)
```

Generates an alphabet of q unique symbols.
"""
genalphabet(q::Int) = Alphabet(gensym(q))

"""
```julia
struct UniverseParameters <: AbstractCode
```

Defines a structure for the messages in the code.  Parameters are the alphabet `Σ`, size of alphabet `q`, and block length `n`

```julia
UniverseParameters(Σ::Alphabet, n::Int)
UniverseParameters(Σ::AbstractArray, n::Int)
UniverseParameters(Σ::Alphabet, q::Int, n::Int)
UniverseParameters(Σ::AbstractArray, q::Int, n::Int)
UniverseParameters(q::Int, n::Int)
```

An inner constructor function on the structure `UniverseParameters`.
"""
struct UniverseParameters <: AbstractCode
    Σ::Alphabet{N} where N
    q::Int # size of alphabet
    n::Int # block length

    UniverseParameters(Σ::Alphabet{N}, n::Int) where {N} = new(Σ, N, n)
    UniverseParameters(Σ::Alphabet{N}, q::Int, n::Int) where {N} = new(Σ, q, n)
    UniverseParameters(Σ::AbstractVector{T}, n::Int) where {T} = (Σ = Alphabet(Σ); new(Σ, length(Σ), n))
    UniverseParameters(Σ::AbstractVector{T}, q::Int, n::Int) where {T} = new(Σ, q, n)
    UniverseParameters(q::Int, n::Int) = new(genalphabet(q), q, n)
end

# Iteration inteface functions
Base.length(iter::UniverseParameters) = big(iter.q)^iter.n
Base.eltype(iter::UniverseParameters) = NTuple{iter.n, Symbol}

# Other Base interface functions for UniverseParameters
Base.rand(𝒰::UniverseParameters) = ntuple(_ -> rand(𝒰.Σ), 𝒰.n)

"""
```julia
rand(𝒰::UniverseParameters, C::AbstractArray) -> Tuple
```

Given universe parameters 𝒰 and a code C, return a tuple including
  - A random word in C;
  - A random letter in the alphabet; and
  - A random index in the block length.
"""
Base.rand(𝒰::UniverseParameters, C::AbstractArray) = tuple(Word(rand(C)), rand(𝒰.Σ), rand(1:𝒰.n))

"""
```julia
struct CodeUniverseIterator <: AbstractCode
```

A structure used to iterate through a code universe, with specified universe parameters.  E.g.,

```julia
for c in CodeUniverseIterator(["a", "b", "c"], 3)
    println(c)
end
```

Fields:
  - `𝒰::UniverseParameters`

Methods:

```julia
CodeUniverseIterator(𝒰::UniverseParameters)
CodeUniverseIterator(Σ::Union{Alphabet, AbstractArray}, q::Int, n::Int)
CodeUniverseIterator(Σ::Union{Alphabet, AbstractArray}, n::Int)
CodeUniverseIterator(q::Int, n::Int)
```
"""
struct CodeUniverseIterator <: AbstractCode
    𝒰::UniverseParameters

    function CodeUniverseIterator(𝒰::UniverseParameters)
        Σ = 𝒰.Σ
        return Iterators.product(Vector{eltype(Σ)}[Σ for _ in 1:𝒰.n]...)
    end

    CodeUniverseIterator(Σ::Union{Alphabet{N}, AbstractVector{T}}, q::Int, n::Int) where {N, T} =
        CodeUniverseIterator(UniverseParameters(Σ, q, n))
    CodeUniverseIterator(Σ::Union{Alphabet{N}, AbstractVector{T}}, n::Int) where {N, T} =
        CodeUniverseIterator(UniverseParameters(Σ, n))
    CodeUniverseIterator(q::Int, n::Int) =
        CodeUniverseIterator(UniverseParameters(q, n))
end

"""
```julia
struct CodeUniverse <: AbstractCode
```

Defines a structure for the messages in the code.  Parameters are the abstract array of messages `𝒰`, and the length of the messages `n`.

```julia
CodeUniverse(𝒰::AbstractArray, Σ::Alphabet)
```

An inner constructor function on the structure `CodeUniverse`.
"""
struct CodeUniverse <: AbstractCode
    𝒰::Codewords{M} where {M}
    Σ::Alphabet{N} where N
    q::Int # alphabet size
    n::Int # block length

    function CodeUniverse(𝒰::Codewords{M}, Σ::Alphabet{N}) where {N, M}
        message_length_error = "We have fixed the block length of each message.  Please ensure all messages are of equal length."
        _allequal_length_(𝒰) || throw(error(message_length_error))

        new(𝒰, Σ, length(Σ), length(first(𝒰)))
    end # end constructor function
end

"""
```julia
struct Rounding end
```

An abstract structure used as a parameter in a non-rounding method of `hamming_bound`.  Use `no_round` for this; a predefiend variable of the structure `Rounding`.
"""
struct Rounding end
no_round = Rounding()


text:
using Quon
cz = tait_cz()
contract!(cz, tait_cz())
Quon.plot_tait(cz)
cz.locations[59] = (2.0, 1.0)
cz.locations[45] = (1.0, 1.0)
cz.locations[60] = (3.0, 1.0)
start_hes = Dict(
    1 => (3, -π/2),
    29 => (6, π),
    24 => (63, -π/2),
    45 => (26, 0),
    59 => (53, 0),
    60 => (54, π),
    57 => (131, π/4),
    32 => (87, π/2),
    56 => (153, π/2)
)
p = plot(cz; start_hes = start_hes, radius = 0.6, show_tait = true)
p |> SVG("1.svg")

m = match(Rule{:perm_rz}(), cz)
Quon.rewrite!(cz, m[3])
Quon.rewrite!(cz, m[9])
start_hes[45] = (4, π/2)
start_hes[60] = (64, π/2)
start_hes[57] = (103, π)
delete!(start_hes, 59)
p = plot(cz; radius = 0.5, start_hes = start_hes, show_tait = true)
p |> SVG("2.svg")

m_string_genus = match(Quon.Rule{:string_genus}(), cz)
Quon.rewrite!(cz, m_string_genus[1])
p = plot(cz; start_hes = start_hes, show_tait = true)
p |> SVG("3.svg")

m_fusion = match(Quon.Rule{:z_fusion}(), cz)
Quon.rewrite!(cz, m_fusion[1])
Quon.rewrite!(cz, m_fusion[2])
Quon.rewrite!(cz, m_fusion[3])
start_hes[57] = (156, -π/4)
p = plot(cz; start_hes = start_hes, show_tait = true)
p |> SVG("4.svg")

m_c_rm = match(Quon.Rule{:charge_rm_f}(), cz)
rewrite!(cz, m_c_rm[1])
start_hes[57] = (156, 0)
cz.locations[32] = (1.0, 3.0)
cz.locations[56] = (3.0, 3.0)
p = plot(cz; start_hes = start_hes, show_tait = true)
p |> SVG("5.svg")

m_gf = match(Rule(:genus_fusion), cz)
rewrite!(cz, m_gf[1])
start_hes[29] = (86, -3π/4)
cz.locations[29] = (2.0, 1.0)
p = plot(cz; start_hes = start_hes, show_tait = true)
p |> SVG("6.svg")


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

"""
    LogGrowthSDE <: AbstractModel

dX_t = r X_t (1- X_t/K) dt + β X_t dB_t
r := exponential growth
K := saturation parameter
β := multiplicative noise factor
a1 := 2*r*r/(β*K) factor often used
a2 := a1/K  factor often used
b1,b2, tt := needed for lambda ratio
"""
struct LogGrowthSDE <: AbstractModel
    r::Float64
    K::Float64
    β::Float64
    a1::Float64 #precomputing some useful factors
    a2::Float64 #precomputing some useful factors
    b1::Vector{Float64} #free vectors
    b2::Vector{Float64} #free vectors
    tt::Vector{Float64} #free vectors
    function LogGrowthSDE(r, K, β, L)
        new(r, K, β, 2*r*r/(β*K) , 2*r*r/(β*K*K), fill(0.0, 2<<L -1), fill(0.0, 2<<L -1), fill(0.0, 2<<L -1))
    end
end

# sampling scheme for each stochastic differential equation
sampling_scheme(::LogGrowthSDE) = SubSampling()
# dependence structure for each stochastic differential equation
dependence_strucute(::LogGrowthSDE) = PartialIndependence()

"""
    waitgaus(ξ, θ, u)

Poisson rate  λ(s) = max(0, θ(ξ + θs))
invert function Λ(t) = int_0^t λ(s) ds + ln(u)
"""
function waitgaus( ξ::Float64, θ::Float64, u::Float64)
    b = ξ*θ
    if b>0
        return Sol2E(0.5, b, log(u))
    else
        return Sol2E(0.5, b, log(u) + ξ*ξ*0.5)
    end
end
"""
    waitexp(a, b, u)

invert function  Λ(t) = a/b (e^(bt) - 1) + ln(u)
"""
function waitexp(a, b, u)
    if (log(u)*b/a > 1)
        return error("exponential gave infinity")
    else
        return log(1-log(u)b/a)/b
    end
end
"""
    λbar(n, S::System, X::LogGrowthSDE , u, v)

Poisson time (upper bound) for the coefficient `n`
of model `LogGrowthSDE` starting at `u` and ending at `v`
"""
function λbar(n, S::System, X::LogGrowthSDE , u::Float64, v::Float64, t::Float64)
    w1 = waitgaus(S.ξ[n], S.θ[n], rand())
    X.b1[n] = minimum(fs_expansion(S.ϕ[n], S.ξ, u, v, S.L, S.T))
    X.b2[n] = minimum(fs_expansion(S.ϕ[n], S.θ, u, v, S.L, S.T))
    X.tt[n] = t
    if S.θ[n] > 0
        w2 = waitexp(0.5*S.ϕ[n].δ*X.a1*exp(-X.β*X.b1[n]), -X.β*X.b2[n], rand())
        return min(w1,w2)
    else
        w2 = waitexp(0.5*S.ϕ[n].δ*X.a2*exp(-2X.β*X.b1[n]), -2X.β*X.b2[n], rand())
        return min(w1, w2)
    end
end

"""
    λratio(n::Int64, S::System, X::LogGrowthSDE, u::Float64, v::Float64)

accept reject time drwan from upper bound λbar relative to the coefficient `n`
of model `LogGrowthSDE` starting at `u` and ending at `v`
"""
function λratio(n::Int64, S::System, X::LogGrowthSDE, u::Float64, v::Float64, t::Float64)
    time = t - X.tt[n]
    ω = MCintegration(S.ϕ[n])
    XX = fs_expansion(S, ω, u, v)
    ϕ = Λ(ω, S.ϕ[n].i, S.ϕ[n].j, S.T)   #need to be changed
    num = max(0, S.θ[n]*(0.5*S.ϕ[n].range*ϕ*(X.a1*exp(-X.β*XX) - X.a2*exp(-2X.β*XX)) + S.ξ[n]))
    den = max(0, S.θ[n]*S.ξ[n]) + max(0, 0.5*S.θ[n]*S.ϕ[n].δ*X.a1*exp(-X.β*X.b1[n])*exp(-X.β*X.b2[n]*time)) + max(0, -0.5*S.θ[n]*S.ϕ[n].δ*X.a2*exp(-2X.β*X.b1[n])*exp(-2X.β*X.b2[n]*time))
    return num/den
end

function sol_exp_growth(u, r, k, T)
    C = u*k/(k - u)
    time = 0:0.001:T
    y = zeros(length(time))
    j = 1
    for i in time
        y[j] = C*exp(r*i)/(1 + C/k*exp(r*i))
        j = j + 1
    end
    return (time, y)
end

function runall(SHORT = false)
    Random.seed!(1)
    T = 200.0
    clock = 200.0
    L = 6
    K = 2000
    r = 0.08
    β = 0.1
    u = -log(50)/β     # end points in the lamperti transform
    v= -log(1000)/β
    X = LogGrowthSDE(r, K, β, L)   #end points in the lamperti tranform
    XX = zz_sampler(X, T, L, u, v, clock)
    if SHORT == false
        burning = 10.0    #burning
        f = clock - 1.0; n = 100
        db = (f-burning)/n
        b =  burning:db:f
        p = plotmixing(XX, b, T, L, u, v, x -> exp(- x*β))
        plot!(p, sol_exp_growth(50, r, K, T), color = :blue)
        xaxis!(p, "t")
        yaxis!(p, "X_t")
        display(p)
    end
    return XX
end

error("STOP HERE")

x = runall()
savefig("../../output/exp_growth_01_01_2000.pdf")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
using GenieWebApp
using Documenter

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

makedocs(;
    modules=[GenieWebApp],
    authors="Jaan Tollander de Balsch",
    repo="https://github.com/csc-training/GenieWebApp.jl/blob/{commit}{path}#{line}",
    sitename="GenieWebApp.jl",
    format=Documenter.HTML(;
        prettyurls = true,
        canonical="https://csc-training.github.io/GenieWebApp.jl",
        assets=String[],
        highlights=String["dockerfile", "nginx"],
    ),
    pages=[
        "About Web Technologies and Cloud Computing" => "index.md",
        "Developing a Genie Application (v3.0)" => [
            joinpath("genie", "development.md"),
            joinpath("genie", "testing.md"),
        ],
        "Deploying to Virtual Machine using OpenStack (v5.5)" => [
            joinpath("openstack", "web-user-interface.md"),
            joinpath("openstack", "command-line-interface.md"),
            joinpath("openstack", "deployment.md"),
            joinpath("openstack", "deletion.md"),
        ],
        "Deploying to OpenShift Container Platform (v3.11)" => [
            joinpath("openshift", "container.md"),
            joinpath("openshift", "web-user-interface.md"),
            joinpath("openshift", "command-line-interface.md"),
            joinpath("openshift", "deletion.md"),
        ],
    ],
)

deploydocs(;
    repo="github.com/csc-training/GenieWebApp.jl",
)

max_stars_repo_name: jaantollander/genie-webapp-db

max_issues_repo_name: jaantollander/genie-webapp-db

max_forks_repo_name: jaantollander/genie-webapp-db


text:
export Block

type Block{T}
    length::Integer
end

Block(s::Symbol, n) = Block{s}(n)

@doc """bootstrap_index! produces an array of integers that produces a draw from
the bootstrap distribution when used to index the original data. The first
argument is a vector of length n or an n by m matrix that will store the index.
The integer n is the number of observations and m is the number of bootstrap
replications. The second argument is a `Block` object that specifies the
blocking method and the block length.

Only the circular and moving blocks bootstrap are supported, which can be
invoked by using `Block(:circular, n)` and `Block(:moving, n)` respectively.

Examples
--------
i = Array(Int, 90) ## initialization

bootstrap_index!(i, Block(:circular, 10))
bootstrap_index!(i, Block(:moving, 9)) ## Overwrites previous line's results
""" ->
function bootstrap_index!{T<:Integer}(index::AbstractVector{T}, block::Block{:circular})
    n = length(index)
    nblocks = cld(n, block.length)
    for i in 1:nblocks
        blockstart_i = rand(1:n)
        for t in ((i-1) * block.length + 1):min(i*block.length, n)
            ## We're being sort of sloppy with the offsets here. The
            ## key issue is that there is a single startpoint for each
            ## block, and that starting point is uniformly distributed
            ## on 1:n. Even though t increases as i increases,
            ## randomizing blockstart_i ensures this uniform
            ## distribution.
            index[t] = mod(blockstart_i + t, n) + 1
        end
    end
end

function bootstrap_index!{T<:Integer}(index::AbstractVector{T}, block::Block{:moving})
    n = length(index)
    nboot = 0
    while nboot < n
        blockstart_i = rand(0:n-1)
        ## Shorten the block if it is larger than the remaining spots
        ## in index or if it runs over the end of the dataset.
        blocklength_i = min(n - blockstart_i, n - nboot, block.length)
        index[nboot + (1:blocklength_i)] = blockstart_i + (1:blocklength_i)
        nboot += blocklength_i
    end
end

function bootstrap_index!{T<:Integer}(index::AbstractMatrix{T}, block::Block)
    for i in 1:size(index, 2)
        bootstrap_index!(sub(index,:,i), block)
    end
end


text:
"""
"""
function is_travis_ci(a::AbstractDict = ENV)
    result = ( lowercase(strip(get(a, "CI", ""))) == "true" ) &&
        ( lowercase(strip(get(a, "TRAVIS", ""))) == "true" ) &&
        ( lowercase(strip(get(a, "CONTINUOUS_INTEGRATION", "")) ) == "true")
    return result
end

"""
"""
function is_travis_ci_on_linux(a::AbstractDict = ENV)
    result = is_travis_ci(a) && Sys.islinux()
    return result
end

"""
"""
# function is_travis_ci_on_apple(a::AbstractDict = ENV)
#     result = is_travis_ci(a) && Sys.isapple()
#     return result
# end

"""
"""
# function is_appveyor_ci(a::AbstractDict = ENV)
#     result = ( lowercase(strip(get(a, "CI", ""))) == "true" ) &&
#         ( lowercase(strip(get(a, "APPVEYOR", ""))) == "true" )
#     return result
# end

"""
"""
is_ci(a::AbstractDict = ENV) = is_travis_ci(a)
# is_ci(a::AbstractDict = ENV) = is_travis_ci(a) || is_appveyor_ci(a)

"""
"""
is_ci_or_runtests(a::AbstractDict = ENV) = is_ci(a) || is_runtests(a)

"""
"""
is_ci_or_runtests_or_docs_or_examples(a::AbstractDict = ENV) =
    is_ci_or_runtests(a) || is_docs_or_examples(a)


text:
module Oceananigans

if VERSION < v"1.1"
    @error "Oceananigans requires Julia v1.1 or newer."
end

export
    # Architectures
    CPU, GPU,

    # Logging
    ModelLogger, Diagnostic, Setup, Simulation,

    # Grids
    RegularCartesianGrid, VerticallyStretchedCartesianGrid,

    # Fields and field manipulation
    Field, CellField, FaceFieldX, FaceFieldY, FaceFieldZ,
    interior, set!,

    # Forcing functions
    ModelForcing, SimpleForcing,

    # Coriolis forces
    FPlane, BetaPlane,

    # Buoyancy and equations of state
    BuoyancyTracer, SeawaterBuoyancy,
    LinearEquationOfState, RoquetIdealizedNonlinearEquationOfState,

    # Surface waves via Craik-Leibovich equations
    SurfaceWaves,

    # Boundary conditions
    BoundaryCondition,
    Periodic, Flux, Gradient, Value, Dirchlet, Neumann,
    CoordinateBoundaryConditions, FieldBoundaryConditions, HorizontallyPeriodicBCs, ChannelBCs,
    BoundaryConditions, SolutionBoundaryConditions, HorizontallyPeriodicSolutionBCs, ChannelSolutionBCs,
    BoundaryFunction, getbc, setbc!,

    # Time stepping
    time_step!,
    TimeStepWizard, update_Δt!,

    # Models
    Model, ChannelModel, NonDimensionalModel,

    # Utilities
    prettytime, pretty_filesize,

    # Turbulence closures
    ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity,
    AnisotropicBiharmonicDiffusivity,
    ConstantSmagorinsky, AnisotropicMinimumDissipation

# Standard library modules
using Printf
using Logging
using Statistics
using LinearAlgebra

# Third-party modules
using Adapt
using OffsetArrays
using FFTW
using JLD2
using NCDatasets

import CUDAapi
import GPUifyLoops

using Base: @propagate_inbounds
using Statistics: mean
using GPUifyLoops: @launch, @loop, @unroll

import Base:
    +, -, *, /,
    size, length, eltype,
    iterate, similar, show,
    getindex, lastindex, setindex!,
    push!

#####
##### Abstract types
#####

"""
    AbstractGrid{T}

Abstract supertype for grids with elements of type `T`.
"""
abstract type AbstractGrid{T} end

"""
    AbstractPoissonSolver

Abstract supertype for solvers for Poisson's equation.
"""
abstract type AbstractPoissonSolver end

"""
    AbstractDiagnostic

Abstract supertype for types that compute diagnostic information from the current model
state.
"""
abstract type AbstractDiagnostic end

"""
    AbstractOutputWriter

Abstract supertype for types that perform input and output.
"""
abstract type AbstractOutputWriter end

#####
##### Place-holder functions
#####

function TimeStepper end
function run_diagnostic end
function write_output end

#####
##### Include all the submodules
#####

include("Architectures.jl")

using Oceananigans.Architectures: @hascuda
@hascuda begin
    # Import CUDA utilities if it's detected.
    using CUDAdrv
    using CUDAnative
    using CuArrays

    println("CUDA-enabled GPU(s) detected:")
    for (gpu, dev) in enumerate(CUDAnative.devices())
        println(dev)
    end
end

include("Utils/Utils.jl")
include("Logger.jl")
include("Grids/Grids.jl")
include("Fields/Fields.jl")
include("Operators/Operators.jl")
include("Coriolis/Coriolis.jl")
include("Buoyancy/Buoyancy.jl")
include("SurfaceWaves.jl")
include("TurbulenceClosures/TurbulenceClosures.jl")
include("BoundaryConditions/BoundaryConditions.jl")
include("Solvers/Solvers.jl")
include("Forcing/Forcing.jl")
include("Models/Models.jl")
include("Diagnostics/Diagnostics.jl")
include("OutputWriters/OutputWriters.jl")
include("TimeSteppers/TimeSteppers.jl")
include("AbstractOperations/AbstractOperations.jl")

#####
##### Re-export stuff from submodules
#####

using .Architectures
using .Utils
using .Grids
using .Fields
using .Coriolis
using .Buoyancy
using .SurfaceWaves
using .TurbulenceClosures
using .BoundaryConditions
using .Solvers
using .Forcing
using .Models
using .TimeSteppers

end # module


text:
"""
    GammaExpert(k, θ)

PDF:

```math
f(x; k, \\theta) = \\frac{x^{k-1} e^{-x/\\theta}}{\\Gamma(k) \\theta^k},
\\quad x > 0
```

See also: [Gamma Distribution](https://en.wikipedia.org/wiki/Gamma_distribution) (Wikipedia), shape-scale parameterization

"""
struct GammaExpert{T<:Real} <: NonZIContinuousExpert
    k::T
    θ::T
    GammaExpert{T}(k::T, θ::T) where {T<:Real} = new{T}(k, θ)
end

function GammaExpert(k::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(GammaExpert, k >= zero(k) && θ > zero(θ))
    return GammaExpert{T}(k, θ)
end

## Outer constructors
GammaExpert(k::Real, θ::Real) = GammaExpert(promote(k, θ)...)
GammaExpert(k::Integer, θ::Integer) = GammaExpert(float(k), float(θ))
GammaExpert() = GammaExpert(1.0, 1.0)

## Conversion
function convert(::Type{GammaExpert{T}}, k::S, θ::S) where {T <: Real, S <: Real}
    GammaExpert(T(k), T(θ))
end
function convert(::Type{GammaExpert{T}}, d::GammaExpert{S}) where {T <: Real, S <: Real}
    GammaExpert(T(d.k), T(d.θ), check_args=false)
end
copy(d::GammaExpert) = GammaExpert(d.k, d.θ, check_args=false)

## Loglikelihood of Expert
logpdf(d::GammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? -Inf : Distributions.logpdf.(Distributions.Gamma(d.k, d.θ), x...)
pdf(d::GammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? 0.0 : Distributions.pdf.(Distributions.Gamma(d.k, d.θ), x...)
logcdf(d::GammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? -Inf : Distributions.logcdf.(Distributions.Gamma(d.k, d.θ), x...)
cdf(d::GammaExpert, x...) = (d.k < 1 && x... <= 0.0) ? 0.0 : Distributions.cdf.(Distributions.Gamma(d.k, d.θ), x...)

## expert_ll, etc
expert_ll_exact(d::GammaExpert, x::Real) = LRMoE.logpdf(d, x)
function expert_ll(d::GammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
    expert_ll = (tu == 0.) ? -Inf : expert_ll
    return expert_ll
end
function expert_tn(d::GammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu))
    expert_tn = (tu == 0.) ? -Inf : expert_tn
    return expert_tn
end
function expert_tn_bar(d::GammaExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::GammaExpert; exposure = 1) = d

## Parameters
params(d::GammaExpert) = (d.k, d.θ)
function params_init(y, d::GammaExpert)
    pos_idx = (y .> 0.0)
    μ, σ2 = mean(y[pos_idx]), var(y[pos_idx])
    θ_init = σ2/μ
    k_init = μ/θ_init
    if isnan(θ_init) || isnan(k_init)
        return GammaExpert()
    else
        return GammaExpert(k_init, θ_init)
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::GammaExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-0.0), (1-0.0)*HypothesisTests.ksstats(y[y .> 0.0], Distributions.Gamma(d.k, d.θ))[2])
end

## Simululation
sim_expert(d::GammaExpert) = Distributions.rand(Distributions.Gamma(d.k, d.θ), 1)[1]

## penalty
penalty_init(d::GammaExpert) = [2.0 10.0 2.0 10.0]
no_penalty_init(d::GammaExpert) = [1.0 Inf 1.0 Inf]
penalize(d::GammaExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.θ) - d.θ/p[4]

## statistics
mean(d::GammaExpert) = mean(Distributions.Gamma(d.k, d.θ))
var(d::GammaExpert) = var(Distributions.Gamma(d.k, d.θ))
quantile(d::GammaExpert, p) = quantile(Distributions.Gamma(d.k, d.θ), p)
lev(d::GammaExpert, u) = d.θ*d.k*gamma_inc(float(d.k+1), u/d.θ, 0)[1] + u*(1-gamma_inc(float(d.k), u/d.θ,0)[1])
excess(d::GammaExpert, u) = mean(d) - lev(d, u)

## Misc functions for E-Step
function _diff_dist_series(d::GammaExpert, yl, yu)
    return cdf.(Distributions.Gamma(d.k+1, d.θ), yu) - cdf.(Distributions.Gamma(d.k+1, d.θ), yl)
end

function _int_obs_Y(d::GammaExpert, yl, yu, expert_ll)
    if yl == yu
        return yl
    else
        return exp(-expert_ll) * ( _diff_dist_series(d, yl, yu) * d.k * d.θ )
    end
end

function _int_lat_Y(d::GammaExpert, tl, tu, expert_tn_bar)
    if tl == tu
        return tl
    else
        return exp(-expert_tn_bar) * ( (1 - _diff_dist_series(d, tl, tu)) * d.k * d.θ )
    end
end

function _int_logy_func(d::GammaExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : log(x) * pdf.(d, x)
end

function _int_obs_logY_raw(d::GammaExpert, yl, yu)
    if yl == yu
        return log(yl) * pdf.(d, yl)
    else
        return quadgk.(x -> _int_logy_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_logY_raw(d::GammaExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_logy_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_logy_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

function _gamma_k_to_θ(k, sum_term_zkz, sum_term_zkzy;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    if penalty
        hyper1, hyper2 = pen_pararms_jk[3], pen_pararms_jk[4]
        if hyper1 == 1.0 && hyper2 == Inf
            # return sum(term_zkzy)[1] / (sum(term_zkz)[1] * k)
            return sum_term_zkzy / (sum_term_zkz * k)
        end
        quad1, quad2 = ( k * sum_term_zkz - (hyper1 - 1) )^2, (4/hyper2) * sum_term_zkzy
        return (hyper2/2) * ( (hyper1-1) - k*sum_term_zkz + sqrt(quad1 + quad2) )
    else
        return sum_term_zkzy / (sum_term_zkz * k)
    end
end

function _gamma_optim_k(logk,
                        sum_term_zkz, sum_term_zkzy, sum_term_zkzlogy;
                        penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])
    # Optimization in log scale for unconstrained computation
    k_tmp = exp(logk)

    θ_tmp = _gamma_k_to_θ(k_tmp, sum_term_zkz, sum_term_zkzy, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = (k_tmp-1)*sum_term_zkzlogy - (1/θ_tmp)*sum_term_zkzy - (k_tmp*log(θ_tmp) + loggamma(k_tmp))*sum_term_zkz
    p = penalty ? (pen_pararms_jk[1]-1)*log(k_tmp) - k_tmp/pen_pararms_jk[2] + (pen_pararms_jk[3]-1)*log(θ_tmp) - θ_tmp/pen_pararms_jk[4] : 0.0
    return (obj + p) * (-1.0)
end

## EM: M-Step
function EM_M_expert(d::GammaExpert,
                     tl, yl, yu, tu,
                     exposure,
                    #  expert_ll_pos,
                    #  expert_tn_pos,
                    #  expert_tn_bar_pos,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Further E-Step
    expert_ll_pos = expert_ll.(d, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(d, tl, yl, yu, tu)

    Y_e_obs = vec(_int_obs_Y.(d, yl, yu, expert_ll_pos))
    Y_e_lat = vec(_int_lat_Y.(d, tl, tu, expert_tn_bar_pos))
    nan2num(Y_e_obs, 0.0) # get rid of NaN
    nan2num(Y_e_lat, 0.0) # get rid of NaN

    yl_yu_unique = unique_bounds(yl, yu)
    int_obs_logY_tmp = _int_obs_logY_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    logY_e_obs = exp.(-expert_ll_pos) .* int_obs_logY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    nan2num(logY_e_obs, 0.0) # get rid of NaN

    tl_tu_unique = unique_bounds(tl, tu)
    int_lat_logY_tmp = _int_lat_logY_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    logY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_logY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    nan2num(logY_e_lat, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (yu .!= 0.0)
    term_zkz = z_e_obs[pos_idx] .+ (z_e_lat[pos_idx] .* k_e[pos_idx])
    term_zkz_Y = (z_e_obs[pos_idx] .* Y_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* Y_e_lat[pos_idx])
    term_zkz_logY = (z_e_obs[pos_idx] .* logY_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* logY_e_lat[pos_idx])

    logk_new = Optim.minimizer( Optim.optimize(x -> _gamma_optim_k(x, sum(term_zkz)[1], sum(term_zkz_Y)[1], sum(term_zkz_logY)[1],
                                               penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                               max(log(d.k)-2.0, 0.0), log(d.k)+2.0 )) # ,
                                               # log(d.k)-2.0, log(d.k)+2.0 )) # ,
                                               # GoldenSection() ))
                                               # , rel_tol = 1e-8) )
    k_new = exp(logk_new)
    θ_new = _gamma_k_to_θ(k_new, sum(term_zkz)[1], sum(term_zkz_Y)[1], penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    # println("$k_new, $θ_new")
    return GammaExpert(k_new, θ_new)
end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::GammaExpert,
                            ye, exposure,
                            z_e_obs;
                            penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Further E-Step
    Y_e_obs = ye
    logY_e_obs = log.(ye)
    # nan2num(logY_e_obs, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (ye .!= 0.0)
    term_zkz = z_e_obs[pos_idx]
    term_zkz_Y = z_e_obs[pos_idx] .* Y_e_obs[pos_idx]
    term_zkz_logY = z_e_obs[pos_idx] .* logY_e_obs[pos_idx]

    logk_new = Optim.minimizer( Optim.optimize(x -> _gamma_optim_k(x, sum(term_zkz)[1], sum(term_zkz_Y)[1], sum(term_zkz_logY)[1],
                                               penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                               max(log(d.k)-2.0, 0.0), log(d.k)+2.0 ))

    k_new = exp(logk_new)
    θ_new = _gamma_k_to_θ(k_new, sum(term_zkz)[1], sum(term_zkz_Y)[1], penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    # println("$k_new, $θ_new")
    return GammaExpert(k_new, θ_new)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
using SimpleValueGraphs.AbstractTuples

# TODO we might want to replace this with a proper wrapper around Adjlist
const Adjlist{T} = Vector{Vector{T}}

# ==========

function Adjlist{T}(n::Integer) where {T}
    result = Adjlist{T}(undef, n)
    @inbounds for i in 1:n
        result[i] = T[]
    end
    return result
end

function deepcopy_adjlist(T::Type, adjlist::Adjlist)
    n = length(adjlist)
    result = Vector{Vector{T}}(undef, n)
    @inbounds for i in OneTo(n)
        list = adjlist[i]
        n_list = length(list)
        result_list = Vector{T}(undef, n_list)
        for j in OneTo(n_list)
            # TODO it seems to be wasteful to always convert & deepcopy
            result_list[j] = deepcopy(convert(T, list[j]))
        end
        result[i] = result_list
    end
    return result
end

# TODO do we have a nicer way to specify the type?
function copy_vertexvals(vertexvals::AbstractTuple)

    return map(copy, vertexvals)
end

# TODO do we have a nicer way to specify the type?
function copy_edgevals(edgevals::AbstractTuple)

    return map(vv -> map(v -> copy(v), vv), edgevals)
end

function reverse_adjlist(adjlist::Adjlist{T}) where {T <: Integer}

    n = length(adjlist)
    result = Vector{Vector{T}}(undef, n)

    # calculate the lengths of the resulting vectors, so that
    # we can allocate the right amount of space from the beginning
    lengths = zeros(T, n)
    @inbounds for i ∈ OneTo(n)
        for v ∈ adjlist[i]
            lengths[v] += one(T)
        end
    end

    @inbounds for i ∈ OneTo(n)
        result[i] = Vector(undef, lengths[i])
    end

    insert_after = fill!(lengths, zero(T)) # reuse space
    @inbounds for i ∈ OneTo(n)
        for v ∈ adjlist[i]
            index = (insert_after[v] += one(T))
            result[v][index] = i
        end
    end

    return result
end

# TODO check if still correct
@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:Tuple}
    R = Tuple{( Adjlist{T} for T in E_VAL.types )...}
    return :($R)
end

# TODO check if still correct
@generated function edgevals_container_type(::Val{E_VAL}) where {E_VAL <:NamedTuple}
    R = NamedTuple{ Tuple(E_VAL.names), Tuple{( Adjlist{T} for T in E_VAL.types )...}}
    return :($R)
end

@generated function set_values_for_index!(edgevals::AbstractTuple, s::Integer, index::Integer, values::T) where {T <: AbstractTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds edgevals[$i][s][index] = values[$i])
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

@generated function insert_values_for_index!(tup_adjlist::AbstractTuple,
                                 s::Integer,
                                 index::Integer,
                                 values::T) where {T <: AbstractTuple}
    len = length(T.types)
    exprs = Expr[]
    for i in OneTo(len)
        e = :(@inbounds insert!(tup_adjlist[$i][s], index, values[$i]))
        push!(exprs, e)
    end
    return Expr(:block, exprs...)
end

function delete_values_for_index!(tup_adjlist::AbstractTuple,
                                 s::Integer,
                                 index::Integer)
    @inbounds for adjlist in tup_adjlist
        deleteat!(adjlist[s], index)
    end
end

# TODO E_VAL could probably be extracted without being explicitly passed to the function
function values_for_index(tup_adjlist::AbstractTuple, E_VAL::Type, s::Integer, index::Integer)
    @inbounds return E_VAL( adjlist[s][index] for adjlist in tup_adjlist )
end


text:
module TestGPA
using GaussianProcesses, Calculus, Distributions
using Test, Random

Random.seed!(1)

@testset "GPA" begin
    d, n = 3, 20
    ll = rand(d)
    X = 2π * rand(d, n)
    f = [sum(sin, view(X, :, i)) / d for i in 1:n]

    # Test Bernoulli, binomial, exponential, Gaussian, Poisson, and Student-t likelihood
    liks = (BernLik(), BinLik(n), ExpLik(), GaussLik(-1.0), PoisLik(), StuTLik(3, 0.1))
    ys = ([Bool(rand(Distributions.Bernoulli(abs(f[i])))) for i in 1:n],
          [rand(Distributions.Binomial(n,exp(f[i]) / (1 + exp(f[i])))) for i in 1:n],
          [rand(Distributions.Exponential(f[i]^2)) for i in 1:n],
          f,
          [rand(Distributions.Poisson(exp(f[i]))) for i in 1:n],
          f .+ rand(Distributions.TDist(3), n))

    # Mean and kernel function
    mZero = MeanZero()
    kern = SE(0.0, 0.0)

    @testset "Likelihood $(typeof(lik))" for (lik, y) in zip(liks, ys)
        println("\tTesting ", nameof(typeof(lik)), "...")

        # Fit GP
        gp = GP(X, y, mZero, kern, lik)

        # Sample random parameters
        params = 0.5 * randn(GaussianProcesses.num_params(gp))
        set_params!(gp, params)

        # Exact gradient
        GaussianProcesses.update_target_and_dtarget!(gp)
        exact_grad = gp.dtarget

        # Numerical approximation
        num_grad = Calculus.gradient(params) do params
            set_params!(gp, params)
            update_target!(gp)
            gp.target
        end

        @test num_grad ≈ exact_grad
    end

end
end


text:
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# corr_sdp.jl
#
# Given three random variables A,B,C and given bounds on two of the three
# correlation coefficients:
# -0.2 <= ρ_AB <= -0.1
#  0.4 <= ρ_BC <=  0.5
# We can use the following property of the correlations:
# [  1    ρ_AB  ρ_AC ]
# [ ρ_AB   1    ρ_BC ]  ≽ 0
# [ ρ_AC  ρ_BC   1   ]
# To determine bounds on ρ_AC by solving a SDP
#############################################################################

using JuMP

m = Model()

@variable(m, X[1:3,1:3], SDP)

# Diagonal is 1s
@constraint(m, X[1,1] == 1)
@constraint(m, X[2,2] == 1)
@constraint(m, X[3,3] == 1)

# Bounds on the known correlations
@constraint(m, X[1,2] >= -0.2)
@constraint(m, X[1,2] <= -0.1)
@constraint(m, X[2,3] >=  0.4)
@constraint(m, X[2,3] <=  0.5)

# Find upper bound
@objective(m, Max, X[1,3])
solve(m)
println("Maximum value is ", getvalue(X)[1,3])
@assert +0.8719 <= getvalue(X)[1,3] <= +0.8720

# Find lower bound
@objective(m, Min, X[1,3])
solve(m)
println("Minimum value is ", getvalue(X)[1,3])
@assert -0.9779 >= getvalue(X)[1,3] >= -0.9799


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - https://github.com/strickek/ReadWriteDlm2.jl

"""
## ReadWriteDlm2
`ReadWriteDlm2` functions `readdlm2()`, `writedlm2()`, `readcsv2()` and
`writecsv2()` are similar to those of stdlib.DelimitedFiles, but with additional
support for `Date`, `DateTime`, `Time`, `Complex`, `Rational`, `Missing` types
and special decimal marks. `ReadWriteDlm2` supports the `Tables` interface.

### `readcsv2(), writecsv2()`:
For "decimal dot" users the functions `readcsv2()` and `writecsv2()` have the
respective defaults: Delimiter is `','` (fixed) and `decimal='.'`.

### `readdlm2(), writedlm2()`:
The basic idea of these functions is to support the "decimal comma countries".
They use `';'` as default delimiter and `','` as default decimal mark.
"Decimal dot" users of these functions need to define `decimal='.'`

### Detailed Documentation:
For more information about functionality and (keyword) arguments see `?help` for
`readdlm2()`, `writedlm2()`, `readcsv2()` and `writecsv2()`.
"""
module ReadWriteDlm2

using Dates
using DelimitedFiles
using DelimitedFiles: readdlm_string, val_opts
import Tables

export readdlm2, writedlm2, readcsv2, writecsv2

include("rwd2_dateregex.jl")
include("rwd2_readutils.jl")
include("rwd2_read.jl")
include("rwd2_writeutils.jl")
include("rwd2_write.jl")
include("rwd2_csv.jl")
include("rwd2_tables.jl")

end # End module ReadWriteDlm2


text:
using StatsBase, StatsFuns, JLD,  Hwloc

counts = Hwloc.histmap(Hwloc.topology_load())
ncores = counts[:Core]
if nworkers()<ncores
    addprocs(ncores   - nprocs())
end

using DynNets
@everywhere using DynNets
@everywhere using StatsFuns
N_est = 900
Tvals = [50 500 5000]#500#
N_ind = length(Tvals)[1]

N = 10 ;
GBAdgp = 1
GBAest =GBAdgp;  useStartVal = false

groupsIndsDgp = DynNets.distributeAinVecN(Array{Int}(1:GBAdgp),N)
groupsIndsEst = groupsIndsDgp

est_times = Array{Float64}(zeros(N_est))
est_conv_flag = falses(N_est)
est_flag = falses(N_est)

sizesPar = [N,GBAest,GBAest]
for i =1:3
    ex1 = "estPar_$(i)  = SharedArray{Float64}(zeros(N_ind,N_est, $(sizesPar[i]) ))"
    ex2 = "simPar_$(i) = SharedArray{Float64}(zeros(N_ind,N_est, $(sizesPar[i])))"

    eval(parse(ex1))
        eval(parse(ex2))
 end

ObsAllT = Array{Array{Int16,3},1}(N_ind)
ParTvDgpAllT = Array{Array{Float64,3},1}(N_ind)

#define save path and list variables to save
save_fold = "./data/estimatesTest/"
    file_name = "EstTest_N_$(N)_T_$(Tvals)_$(N_est)_GBAdgp_$(GBAdgp)_GBAest_$(GBAest)_$(useStartVal)_target.jld"
    save_path = save_fold*file_name#

est_times = SharedArray{Float64}(zeros(N_ind,N_est))
est_conv_flag = SharedArray{Bool,2}(N_ind,N_est)
est_flag = SharedArray{Bool,2}(N_ind,N_est)

simUncMeans = DynNets.linSpacedFitnesses(N,N)
sim_B = Array(linspace(0.4,0.9,GBAdgp+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GBAdgp+1)[2:end]) #
sim_W = simUncMeans.*(1-sim_B[groupsIndsDgp])
@progress "T"  for ind = 1:N_ind
#   ind=1

    T = Tvals[ind]
    ObsTmp = SharedArray{Int16,3}(Tvals[ind],N,N_est)
    ParTvDgpTmp = SharedArray{Float64,3}(Tvals[ind],N,N_est)

    @sync @parallel   for n = 1:N_est
#     n=1
        simPar_3[ind,n,1:GBAdgp] = sim_A
        simPar_2[ind,n,1:GBAdgp] = sim_B
        simPar_1[ind,n,:] = sim_W
        Model2Sample = DynNets.SdErgm1(ones(T,N),[sim_W,sim_B,sim_A],
                                                    groupsIndsDgp)

        Y_T,fit_T  = DynNets.sample(Model2Sample)
        sampleDegsT = squeeze(sum(Y_T,2),2)

        ParTvDgpTmp[:,:,n] = fit_T
        ObsTmp[:,:,n] = sampleDegsT

        Model2Est =  DynNets.SdErgm1(sampleDegsT,[zeros(N),zeros(GBAest),zeros(GBAest)],
                                                    groupsIndsEst)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[ind,n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimateTargeting(Model2Est )

        estPar_1[ind,n,:] = hat_pars[1]
        estPar_2[ind,n,:] = hat_pars[2]
        estPar_3[ind,n,:] = hat_pars[3]
        est_conv_flag[ind,n] = conv_flag
        est_flag[ind,n] = true

   end

   ObsAllT[ind] = ObsTmp
   ParTvDgpAllT[ind] = ParTvDgpTmp

    # save

    @save(save_path,groupsIndsEst,ObsAllT,ParTvDgpAllT,simPar_1,simPar_2,simPar_3,
            estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)
    println([N,T,sum(est_times[ind,:]),sum(est_conv_flag[ind,:])/sum(est_flag[ind,:])])
end

##

 ##


text:
interc = readdlm("/home/alkan/c/nese/jtb/Adj_matrix_yeast.dat",' ',Int64);

n = length(interc[:,1])                 # number of nodes
selfact = diag(interc)                  # yeast specific
interc -= diagm(selfact)                # interaction matrix

# --- generate the symmetric adjacency matrix with elements 0,1 ---
function symadj(m)
    n = length(m[:,1])
    ad = zeros(Int64,n,n)
    for (i=1:n) for (j=i+1:n)       # adjacency matrix - without diagonals
        ad[i,j] = m[i,j]==0 && m[j,i]==0 ? 0 : 1
        ad[j,i] = ad[i,j]
    end end
    return ad
end

## ad = symadj(interc)

## --- check connectivity -----------
function isconnected (adj)
    n = length(adj[:,1])                # number of nodes
    nsq = n*n
    dum = adj
    (sumprev,sumdum) = (0,0)
    for (i=1:n)
        dum = sign(dum + dum*adj)
        sumdum = sum(dum)
        if (sumdum == sumprev)
            break
        end
        sumprev = sumdum
    end
    return (sumdum==nsq)
end

## --- node update function (integrator) -----------
function upd(s,intr,slf,n)
    sss = similar(s)
    inp = intr*s
    for (i=1:n)
        sss[i] = (inp[i]>0)?1:(inp[i]<0)?0:(slf[i]==0)?0:s[i]
    end
    return sss
end

## --- check fixed point ---
function isfixed (s,intr,slf,n)
    snew = upd(s,intr,slf,n)
    return (snew==s)
end

## ---- find point attractors of a small network ---
function getpointattr(intr,slf)
    n = length(slf)
    attrlist = zeros(Int64,n)                 # s=0 is always an attractor
    for (i=1:2^n-1)                           # visit each state except s=0
        s = digits(parseint(bits(i)))
        lens = length(s)
        if (lens < n)
            s = [s;zeros(Int64,n-lens)]
        end
        if (isfixed(s,intr,slf,n))
            attrlist = [attrlist s]
        end
    end
    return attrlist
end

## ---- find ALL attractors of a small network ---
## ---- by exhaustive search ---
## ---- output is the attractor list + periods ---
function getcycles(intr,slf)
    n = length(slf)
    nm1 = n-1
    attrlist = zeros(Int64,n)           # s=0 is always an attractor
    keyzeros = join(zeros(Int64,n))
    visdict = [keyzeros=>keyzeros] # add the zero fixed point to visited list
    attdict = [keyzeros=>1]    # values are the basin sizes
    for (i=1:2^n-1)            # visit each state except s=0
        skey = bits(i)[end-nm1:end]
        s = [int(skey[x])-48 for x = 1:n]
        if (!haskey(visdict,skey))      # if not seen this state before
            pathdict = Dict()
            while (!haskey(visdict,skey))   # generate a path starting here
                pathdict[skey] = "-"
                visdict[skey] = "-"
                s = upd(s,intr,slf,n)
                skey=join(s)         # join(s) == bits(i)[end-nm1:end]
            end
            attkey=visdict[skey]
            if (haskey(attdict,attkey)) # reached a known cycle or fp
                for (k in keys(pathdict))
                    visdict[k]=attkey
                    attdict[attkey] += length(keys(pathdict)) # basin
                end
            else                        # this is a new cycle of fp
                for (k in keys(pathdict)) visdict[k] = skey end
                attdict[skey]=length(keys(pathdict)) # basin
                attrlist = [attrlist s]
            end
        end
    end
    basinlist = [attdict[join(attrlist[:,i])]/(2^n-1) for i in 1:length(attrlist[1,:])]
    return attrlist,basinlist
end

## ---- output is the attractor list + basin sizes ---
## ---- calculated from a subset of initial conditions ---
function getcyclessample(intr,slf,ninit)

    n = length(slf)
    nm1 = n-1
    attrlist = zeros(Int64,n)           # s=0 is always an attractor
    keyzeros = join(zeros(Int64,n))
    attdict = [keyzeros=>0]    # values are the basin sizes
    cyclen = [keyzeros=>1]    # values are the cycle lengths
    for (i=1:ninit)       # calculate basins using ninit initial states
        s = rand(0:1,n)
        skey = join(s)
        pathdict = Dict()
        pathlen = 1
        while (!haskey(attdict,skey) && !haskey(pathdict,skey))
            pathdict[skey] = pathlen
            pathlen += 1
            s = upd(s,intr,slf,n)
            skey=join(s)         # join(s) == bits(i)[end-nm1:end]
        end
        if (haskey(attdict,skey)) attdict[skey] += 1
        else                        # this is a new cycle of fp
            attdict[skey] = 1
            cyclen[skey] = pathlen - pathdict[skey]
            attrlist = [attrlist s]
        end
    end
    basinlist = [attdict[join(attrlist[:,i])]/ninit for i in 1:length(attrlist[1,:])]
    cyclelengthlist = [cyclen[join(attrlist[:,i])] for i in 1:length(attrlist[1,:])]
    return attrlist,basinlist,cyclelengthlist
end

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
using CoexistenceHoles
using Test
import LinearAlgebra: norm

@testset "random_growthvector.jl" begin
    N = 8; μ = 0.1; σ = 0.1; # create growth vectors from these specs

    # create some vectors to play with
    r1 = random_growthvector(N, μ, σ, seed=1234);
    r2 = random_growthvector(N, μ, σ, seed=1234);
    r3 = random_growthvector(N, μ, σ, seed=1235);

    # check that the seed is reproducible
    @test r1 == r2
    @test r1 != r3

end

@testset "ransomize_growthvector.jl" begin
    r = [0.1, -0.3, 2.1, -1, 0.5];
    # create randomizations of r
    r_rand1 = randomize_growthvector(r, method="preserve_norm", seed=1234);
    r_rand2 = randomize_growthvector(r, method="preserve_norm", seed=1234);
    r_rand3 = randomize_growthvector(r, method="preserve_norm", seed=1235);

    @test r_rand1 == r_rand2 # reproducibility
    @test r_rand1 != r_rand3 # reproducibility
    @test norm(r_rand1) - norm(r) < 1e-10 # norm_preserved

    r_psign_samp1 = randomize_growthvector(r, method="preserve_sign_sample", seed=1234);
    r_psign_samp2 = randomize_growthvector(r, method="preserve_sign_sample", seed=1234);
    r_psign_samp3 = randomize_growthvector(r, method="preserve_sign_sample", seed=1235);

    @test r_psign_samp1 == r_psign_samp2 # reproducibility
    @test r_psign_samp1 != r_psign_samp3 # reproducibility
    @test sign.(r) == sign.(r_psign_samp1) # same signs
    @test isempty(filter(x -> !in(x, abs.(r)), abs.(r_psign_samp1))) # contains same values (perhaps different order)

    r_psign_shuf1 = randomize_growthvector(r, method="preserve_sign_shuffle", seed=1234);
    r_psign_shuf2 = randomize_growthvector(r, method="preserve_sign_shuffle", seed=1234);
    r_psign_shuf3 = randomize_growthvector(r, method="preserve_sign_shuffle", seed=1235);

    @test r_psign_shuf1 == r_psign_shuf2 # reproducibility
    @test r_psign_shuf1 != r_psign_shuf3 # reproducibility
    @test sign.(r) == sign.(r_psign_shuf1) # same signs
    @test isempty(filter(x -> !in(x, abs.(r)), abs.(r_psign_shuf1))) && isempty(filter(x -> !in(x, abs.(r_psign_shuf1)), abs.(r))) # both contain the values of the other (i.e. no new value introduced in the shuffle)

    r_samp1 = randomize_growthvector(r, method="sample", seed=1234);
    r_samp2 = randomize_growthvector(r, method="sample", seed=1234);
    r_samp3 = randomize_growthvector(r, method="sample", seed=1235);

    @test r_samp1 == r_samp2 # reproducibility
    @test r_samp1 != r_samp3 # reproducibility
    @test isempty(filter(x -> !in(x, abs.(r)), abs.(r_samp1))) # contains same values (perhaps different order)

    r_shuf1 = randomize_growthvector(r, method="shuffle", seed=1234);
    r_shuf2 = randomize_growthvector(r, method="shuffle", seed=1234);
    r_shuf3 = randomize_growthvector(r, method="shuffle", seed=1235);

    @test r_shuf1 == r_shuf2 # reproducibility
    @test r_shuf1 != r_shuf3 # reproducibility
    @test isempty(filter(x -> !in(x, abs.(r)), abs.(r_shuf1))) && isempty(filter(x -> !in(x, abs.(r_shuf1)), abs.(r))) # both contain the values of the other (i.e. no new value introduced in the shuffle)

end

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes


text:
export demo

using Gadfly

function demo()
# Setup
    x = collect(1.0:1.0:20)
    u = x.*0.0
    s = x
    v = x.*0.0
    t = x
    w = x.*0.0

    for i in 6:10
        u[i] = 1.0
    end
    for i in 11:14
        u[i] = -0.2 * i + 3.0
    end
    for i in 11:14
        v[i] = -0.2
    end

    p = collect(1.0:0.2:20)
    r = p.*0.0
    tol = 1e-7
    i = 0
    for pi in p
        i += 1
        if pi >= (6.0 - tol) && pi <= (10.0 + tol)
            r[i] = 1.0
        end
        if pi > (10.0 + tol) && pi < (15.0 - tol)
            r[i] = -0.2 * p[i] + 3.0
        end
    end
####
    interpolate(x, u, RK_W3())

    cond = get_cond()
    ε = get_epsilon()
    println("cond = $cond, ε = $ε")

    σ = evaluate(p)

    δ = σ - r
    rms = round(norm(δ)/sqrt(length(δ)); digits=2)
    println("rms = $rms")

    set_default_plot_size(13cm, 13cm)
    plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
          layer(x = p, y = r, Geom.line, Theme(default_color=colorant"red")),
          layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
          Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
          Guide.manual_color_key("Legend", ["Points", "True", "Spline"], ["orange", "red", "blue"]),
          Guide.xlabel("knots, points"), Guide.ylabel("f, σ"),
          Guide.title("Fig.1a"))
    draw(SVG("plot_W3.svg", 13cm, 13cm), plt)

    σ = evaluate(p, 1)
    plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
          layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
          Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
          Guide.manual_color_key("Legend", ["Points", "Spline 1st derivative"], ["orange", "blue"]),
          Guide.xlabel("knots, points"), Guide.ylabel(raw"σ'"),
          Guide.title("Fig.2a"))
    draw(SVG("plot_W3_d1.svg", 13cm, 13cm), plt)

    σ = evaluate(p, 2)
    plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
          layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
          Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
          Guide.manual_color_key("Legend", ["Points", "Spline 2nd derivative"], ["orange", "blue"]),
          Guide.xlabel("knots, points"), Guide.ylabel(raw"σ''"),
          Guide.title("Fig.3a"))
    draw(SVG("plot_W3_d2.svg", 13cm, 13cm), plt)

###
      interpolate(x, u, s, v, t, w, RK_W3())

      cond = get_cond()
      ε = get_epsilon()
      println("cond = $cond, ε = $ε")

      σ = evaluate(p)

      δ = σ - r
      rms = round(norm(δ)/sqrt(length(δ)); digits=2)
      println("rms = $rms")

      set_default_plot_size(13cm, 13cm)
      plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
            layer(x = p, y = r, Geom.line, Theme(default_color=colorant"red")),
            layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
            Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
            Guide.manual_color_key("Legend", ["Points", "True", "Spline"], ["orange", "red", "blue"]),
            Guide.xlabel("knots, points"), Guide.ylabel("f, σ"),
            Guide.title("Fig.1b"))
      draw(SVG("plot_d12_W3.svg", 13cm, 13cm), plt)

      σ = evaluate(p, 1)
      plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
            layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
            Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
            Guide.manual_color_key("Legend", ["Points", "Spline 1st derivative"], ["orange", "blue"]),
            Guide.xlabel("knots, points"), Guide.ylabel(raw"σ'"),
            Guide.title("Fig.2b"))
      draw(SVG("plot_d12_W3_d1.svg", 13cm, 13cm), plt)

      σ = evaluate(p, 2)
      plt = plot(layer(x = x, y = u, Geom.point, Theme(default_color=colorant"orange")),
            layer(x = p, y = σ, Geom.line, Theme(default_color=colorant"blue")),
            Scale.y_continuous(minvalue=-0.5, maxvalue=1.5),
            Guide.manual_color_key("Legend", ["Points", "Spline 2nd derivative"], ["orange", "blue"]),
            Guide.xlabel("knots, points"), Guide.ylabel(raw"σ''"),
            Guide.title("Fig.3b"))
      draw(SVG("plot_d12_W3_d2.svg", 13cm, 13cm), plt)

end


text:
## Lambidfy an expression
## https://github.com/jverzani/SymPy.jl/issues/60
## Hack until SymPy -> Julia converter written
_func(x) = (x.x)[:func]
_funcname(x) = _func(x)[:__name__]
_args(x) = (x.x)[:args]
_srepr(x) = sympy_meth(:srepr, x)

## Mapping of SymPy Values into julia values
val_map = Dict(
               "Zero"             => :(0),
               "One"              => :(1),
               "NegativeOne"      => :(-1),
               "Half"             => :(1/2),
               "Pi"               => :pi,
               "Exp1"             => :e,
               "Infinity"         => :Inf,
               "NegativeInfinity" => :(-Inf),
               "ComplexInfinity"  => :Inf, # error?
               "ImaginaryUnit"    => :im,
               "BooleanTrue"      => :true,
               "BooleanFalse"     => :false
               )

## Mapping of Julia function names into julia ones
## most are handled by Symbol(fnname), this catches exceptions
_heaviside(x) = 1//2 * (1 + sign(x))
function _piecewise(args...)
    as = copy([args...])
    val, cond = pop!(as)
    ex = Expr(:call, :ifelse, cond, val, :nothing)
    while length(as) > 0
        val, cond = pop!(as)
        ex = Expr(:call, :ifelse, cond, val, ex)
    end
    ex
end

fn_map = Dict(
              "Add" => :+,
              "Sub" => :-,
              "Mul" => :*,
              "Div" => :/,
              "Pow" => :^,
              "re"  => :real,
              "im"  => :imag,
              "Abs" => :abs,
              "Min" => :min,
              "Max" => :max,
              "Poly" => :identity,
              "Heaviside" => :(_heaviside),
              "Piecewise" => :(_piecewise),
              "And" => :(&),
              "Or" => :(|),
              "Less" => :(<),
              "LessThan" => :(<=),
              "Equal" => :(==),
              "GreaterThan" => :(>=),
              "Greater" => :(>)
              )

map_fn(key, fn_map) = haskey(fn_map, key) ? fn_map[key] : @compat(Symbol(key))

## TODO? In SymPy, one can pass in dictionary keys to replace functions
## http://docs.sympy.org/dev/tutorial/manipulation.html
function walk_expression(ex; values=Dict(), fns=Dict())

    fns_map = merge(fn_map, fns)
    vals_map = merge(val_map, values)

    fn = _funcname(ex)

    if fn == "Symbol"
        return @compat(Symbol(string(ex)))
    elseif fn in ["Integer" , "Float"]
        return N(ex)
    elseif fn == "Rational"
        return convert(Int, numer(ex))//convert(Int, denom(ex))
        ## piecewise requires special treatment
    elseif fn == "Piecewise"
        return _piecewise([walk_expression(cond) for cond in _args(ex)]...)
    elseif fn == "ExprCondPair"
        val, cond = _args(ex)
        return (val, walk_expression(cond))
    elseif haskey(vals_map, fn)
        return vals_map[fn]
    end

    as = _args(ex)

    Expr(:call, map_fn(fn, fns_map), [walk_expression(a) for a in as]...)
end

"""

`lambidfy(ex::Sym,[vars])`: Lambidfy an expression returning a native
Julia function.  SymPy's
[lambdify](http://docs.sympy.org/dev/modules/utilities/lambdify.html)
function translates code into Python, this translates an expression
into a `Julia` function.

Evaluating the function does not call into SymPy, so should be much faster.

The optional `[vars]` specifies the order of the variables when more
than one is in `ex`. The default is to use the ordering of
`free_symbols(ex)`.

The keyword aruguments allow for the passing of expressions that are
not covered by the default ones. These are dictionaries whose keys are
strings with a SymPy name and whose values are symbols representing
`Julia` values. For examples `Dict("sin"=>:sin)` could be used to map
a function, were that not already done.

Not all expressions can be lambdified. If not, an error is thrown.

Some simple examples

```
@vars x y
lambdify(x^2)(2)       # 4
lambdify(x*y^2)(2,3)   # 2*3^2 using default ordering
lambdify(x*y^2, [y, x])(2,3) # 3*2^2, as function is (y,x) -> x*y^2 equivalent in Julia
```

Compare times
```
xs = rand(1000)
@vars x
ex = sin(x)*cos(2x) * exp(x^2/2)
map(u -> N(ex(u)), xs)   # 3.435850 seconds
SymPy.mapsubs(ex, x, xs) # 0.008569 seconds (does calculuations in Python)
map(lambdify(ex), xs)    # 0.007085 seconds
```

This is a *temporary* solution. The proper fix is to do this in SymPy.

"""
function lambdify(ex::Sym, vars=free_symbols(ex); fns=Dict(), values=Dict())
    # if :julia_code printer is there, use it
    if haskey(sympy, :julia_code)
        body = parse(sympy_meth(:julia_code, ex))
    else
        body = walk_expression(ex, fns=fns, values=values)
    end
    try
        eval(Expr(:function,
                  Expr(:call, gensym(), map(Symbol,vars)...),
                  body))
    catch err
        throw(ArgumentError("Expression does not lambdify"))
    end
end

export(lambdify)

function init_lambdify()
    if haskey(sympy, :julia_code)
        global julia_code(ex::Sym; assign_to=nothing, kwargs...) = sympy_meth(:julia_code, ex; assign_to=assign_to, kwargs...)
        eval(Expr(:export, :julia_code))
    end
end


text:
using Printf
import Base: getindex, *, ndims
import LinearAlgebra: ×

abstract type AbstractDomain end
abstract type AbstractBoundary end

export Interval, Circle
export PointDomain, IntervalDomain, ProductDomain, AbstractDomain
export info, ndims, periodicity_function
export ×
export ∂

struct DomainBoundary <: AbstractBoundary
    closure
end

struct PointDomain{S} <: AbstractDomain
    point::S
end
struct IntervalDomain{AT, BT, PT} <: AbstractDomain
    a::AT
    b::BT
    periodic::PT
end

function IntervalDomain(a, b; periodic=false)
    @assert a < b
    return IntervalDomain(a, b, periodic)
end

function Circle(a, b)
    @assert a < b
    return IntervalDomain(a, b, periodic = true)
end

S¹ = Circle

function Interval(a, b)
    @assert a < b
    return IntervalDomain(a, b)
end

function Base.show(io::IO, Ω::IntervalDomain)
    a = Ω.a
    b = Ω.b
    printstyled(io, "[", color = 226)
    printstyled("$a, $b", color = 7)
    Ω.periodic ? printstyled(io, ")", color = 226) : printstyled(io, "]", color = 226)
 end

function Base.show(io::IO, o::PointDomain)
    printstyled("{",o.point,"}", color = 201)
end

# Product Domains
struct ProductDomain{DT} <: AbstractDomain
    domains::DT
end

function Base.show(io::IO, Ω::ProductDomain)
    for (i,domain) in enumerate(Ω.domains)
        print(domain)
        if i != length(Ω.domains)
            printstyled(io, "×", color = 118)
        end
    end
 end

ndims(p::PointDomain) = 0
ndims(Ω::IntervalDomain) = 1
ndims(Ω::ProductDomain) = +(ndims.(Ω.domains)...)

×(arg1::AbstractDomain, arg2::AbstractDomain) = ProductDomain((arg1, arg2))
×(args::ProductDomain, arg2::AbstractDomain) = ProductDomain((args.domains..., arg2))
×(arg1::AbstractDomain, args::ProductDomain) = ProductDomain((arg1, args.domains...))
×(arg1::ProductDomain, args::ProductDomain) = ProductDomain((arg1.domains..., args.domains...))
×(args::AbstractDomain) = ProductDomain(args...)
*(arg1::AbstractDomain, arg2::AbstractDomain) = ProductDomain((arg1, arg2))
*(args::ProductDomain, arg2::AbstractDomain) = ProductDomain((args.domains..., arg2))
*(arg1::AbstractDomain, args::ProductDomain) = ProductDomain((arg1, args.domains...))
*(arg1::ProductDomain, args::ProductDomain) = ProductDomain((arg1.domains..., args.domains...))
*(args::AbstractDomain) = ProductDomain(args...)

function info(Ω::ProductDomain)
    println("This is a ", ndims(Ω),"-dimensional tensor product domain.")
    print("The domain is ")
    println(Ω, ".")
    for (i,domain) in enumerate(Ω.domains)
        domain_string = domain.periodic ? "periodic" : "wall-bounded"
        length = @sprintf("%.2f ", domain.b-domain.a)
        println("The dimension $i domain is ", domain_string, " with length ≈ ", length)
    end
    return nothing
end

function check_full_periodicity(Ω::ProductDomain)
    b = [Ω.domains[i].periodic for i in eachindex(Ω.domains)]
    return prod(b)
end

function periodicity_function(Ω::ProductDomain)
    periodicity = ones(Bool, ndims(Ω))
    for i in 1:ndims(Ω)
        periodicity[i] = Ω[i].periodic
    end
    return Tuple(periodicity)
end

getindex(Ω::ProductDomain, i::Int) = Ω.domains[i]

# Boundaries
struct Boundaries{S}
    boundaries::S
end

getindex(∂Ω::Boundaries, i) = ∂Ω.boundaries[i]

function Base.show(io::IO, ∂Ω::Boundaries)
    for (i,boundary) in enumerate(∂Ω.boundaries)
        printstyled("boundary ", i, ": ", color = 13)
        println(boundary)
    end
end

function ∂(a::IntervalDomain)
    if a.periodic
        return (nothing)
    else
        return Boundaries((PointDomain(a.a), PointDomain(a.b)))
    end
    return nothing
end

function ∂(Ω::ProductDomain)
    ∂Ω  = []
    for domain in Ω.domains
        push!(∂Ω, ∂(domain))
    end
    splitb = []
    for (i, boundary) in enumerate(∂Ω)
        tmp = Any[]
        push!(tmp, Ω.domains...)
        if boundary != nothing
            tmp1 = copy(tmp)
            tmp2 = copy(tmp)
            tmp1[i] = boundary[1]
            push!(splitb, ProductDomain(Tuple(tmp1)))
            tmp2[i] = boundary[2]
            push!(splitb, ProductDomain(Tuple(tmp2)))
        end
    end
    return Boundaries(Tuple(splitb))
end


text:
"""
## hexahedron_xz!

This subroutine generates nodal coordinates and numbering for
8, 14 or 20-node "bricks" counting x-z planes in the y-direction.

### Method
```julia
hexahedron_xz!(iel, x_coords, y_coords, z_coords, coord, num)
```

### Arguments
```julia
* iel::Int                       : Element number
* x_coords::FloatRange{Float64}    : x coordinates
* y_coords::FloatRange{Float64}    : y coordinates
* z_coords::FloatRange{Float64}    : y coordinates
* coord::Matrix{Float64}           : Nodal coordinates (Updated)
* num::Vector{Int}               : Node numbers (Updated)
* dir::Symbol                      : Node numbering direction
```
"""
function hexahedron_xz!(iel::Int, x_coords::Vector{Float64},
  y_coords::Vector{Float64}, z_coords::Vector{Float64},
  coord::LinearAlgebra.Adjoint{Float64,Array{Float64,2}}, num::Vector{Int})
  #
  # This subroutine generates nodal coordinates and numbering for
  # 8, 14 or 20-node "bricks" counting x-z planes in the y-direction.
  #
  # Updates coord, num
  #
  nxe = size(x_coords, 1) - 1
  nze = size(z_coords, 1) - 1
  nod = size(num, 1)
  iq = (iel-1)÷(nxe*nze) + 1
  iplane = iel-(iq-1)*nxe*nze
  is = (iplane-1)÷nxe + 1
  ip = iplane-(is-1)*nxe
  if nod == 8
   fac1=(nxe+1)*(nze+1)*(iq-1)
   num[1] = fac1+is*(nxe+1) + ip
   num[2] = fac1+(is-1)*(nxe+1) + ip
   num[3] = num[2] + 1
   num[4] = num[1] + 1
   num[5] = (nxe+1)*(nze+1)*iq+is*(nxe+1) + ip
   num[6] = (nxe+1)*(nze+1)*iq+(is-1)*(nxe+1) + ip
   num[7] = num[6] + 1
   num[8] = num[5] + 1

   coord[1:2,1] = x_coords[ip]
   coord[5:6,1] = x_coords[ip]
   coord[3:4,1] = x_coords[ip+1]
   coord[7:8,1] = x_coords[ip+1]

   coord[1:4,2] = y_coords[iq]
   coord[5:8,2] = y_coords[iq+1]

   coord[2:3,3] = z_coords[is]
   coord[6:7,3] = z_coords[is]
   coord[1:3:4,3] = z_coords[is+1]
   coord[5:3:8,3] = z_coords[is+1]

  elseif nod == 14
   fac1=(2*nxe+1)*(2*nze+1)*(iq-1)
   fac2=(2*nxe+1)*(2*nze+1)*iq
   num[1] = fac1+is*(2*nxe+1) + ip
   num[2] = num[1] - (2*nxe+1)
   num[3] = num[2] + 1
   num[4] = num[1] + 1
   num[5] = num[2] + nxe + 1
   num[6] = fac1+(nxe+1)*(nze+1) + nxe*nze+(is-1)*(2*nxe+1) + nxe+ip
   num[7] = num[6] - nxe
   num[8] = num[6] + 1
   num[9] = num[8] + nxe
   num[10] = fac2+is*(2*nxe+1) + ip
   num[11] = num[10] - (2*nxe+1)
   num[12] = num[11] + 1
   num[13] = num[10] + 1
   num[14] = num[11] + nxe + 1

   coord[1:2,1] = x_coords[ip]
   coord[6,1] = x_coords[ip]
   coord[10:11,1] = x_coords[ip]
   coord[5:9:2,1] = 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[14,1] = 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[3:4,1] = x_coords[ip+1]
   coord[8,1] = x_coords[ip+1]
   coord[12:13,1] = x_coords[ip+1]

   coord[1:5,2] = y_coords[iq]
   coord[6:9,2] = 0.5*(y_coords[iq] + y_coords[iq+1])
   coord[10:14,2] = y_coords[iq+1]

   coord[2:3,3] = z_coords[is]
   coord[7,3] = z_coords[is]
   coord[11:12,3] = z_coords[is]
   coord[5:6,3] = 0.5*(z_coords[is] + z_coords[is+1])
   coord[8:6:14,3] = 0.5*(z_coords[is] + z_coords[is+1])
   coord[1:4:3,3] = z_coords[is+1]
   coord[9,3] = z_coords[is+1]
   coord[10:3:13,3] = z_coords[is+1]

  elseif nod == 20
   fac1=((2*nxe+1)*(nze+1) + (2*nze+1)*(nxe+1))*(iq-1)
   fac2=((2*nxe+1)*(nze+1) + (2*nze+1)*(nxe+1))*iq
   num[1] = fac1+(3*nxe+2)*is+2*ip-1
   num[2] = fac1+(3*nxe+2)*is-nxe+ip-1
   num[3] = num[1] - 3*nxe-2
   num[4] = num[3] + 1
   num[5] = num[4] + 1
   num[6] = num[2] + 1
   num[7] = num[1] + 2
   num[8] = num[1] + 1
   num[9] = fac2-(nxe+1)*(nze+1) + (nxe+1)*is+ip
   num[10] = num[9]-nxe-1
   num[11] = num[10] + 1
   num[12] = num[9] + 1
   num[13] = fac2+(3*nxe+2)*is+2*ip-1
   num[14] = fac2+(3*nxe+2)*is-nxe+ip-1
   num[15] = num[13] - 3*nxe-2
   num[16] = num[15] + 1
   num[17] = num[16] + 1
   num[18] = num[14] + 1
   num[19] = num[13] + 2
   num[20] = num[13] + 1

   coord[1:3,1] .= x_coords[ip]
   coord[9:10,1] .= x_coords[ip]
   coord[13:15,1] .= x_coords[ip]
   coord[5:7,1] .= x_coords[ip+1]
   coord[11:12,1] .= x_coords[ip+1]
   coord[17:19,1] .= x_coords[ip+1]
   coord[4:4:8,1] .= 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[16:4:20,1] .= 0.5*(x_coords[ip] + x_coords[ip+1])

   coord[1:8,2] .= y_coords[iq]
   coord[13:20,2] .= y_coords[iq+1]
   coord[9:12,2] .= 0.5*(y_coords[iq] + y_coords[iq+1])

   coord[1,3] = z_coords[is+1]
   coord[7:9,3] .= z_coords[is+1]
   coord[12:13,3] .= z_coords[is+1]
   coord[19:20,3] .= z_coords[is+1]
   coord[3:5,3] .= z_coords[is]
   coord[10:11,3] .= z_coords[is]
   coord[15:17,3] .= z_coords[is]
   coord[2:4:6,3] .= 0.5*(z_coords[is] + z_coords[is+1])
   coord[14:4:18,3] .= 0.5*(z_coords[is] + z_coords[is+1])

  else
    println("Wrong number of nodes for hexahedral fin_el.")
    exit(1)
  end
end

function hexahedron_xz!(iel::Int, x_coords::Vector{Float64},
  y_coords::Vector{Float64}, z_coords::Vector{Float64},
  coord::Array{Float64,2}, num::Vector{Int})
  #
  # This subroutine generates nodal coordinates and numbering for
  # 8, 14 or 20-node "bricks" counting x-z planes in the y-direction.
  #
  # Updates coord, num
  #
  nxe = size(x_coords, 1) - 1
  nze = size(z_coords, 1) - 1
  nod = size(num, 1)
  iq = (iel-1)÷(nxe*nze) + 1
  iplane = iel-(iq-1)*nxe*nze
  is = (iplane-1)÷nxe + 1
  ip = iplane-(is-1)*nxe
  if nod == 8
   fac1=(nxe+1)*(nze+1)*(iq-1)
   num[1] = fac1+is*(nxe+1) + ip
   num[2] = fac1+(is-1)*(nxe+1) + ip
   num[3] = num[2] + 1
   num[4] = num[1] + 1
   num[5] = (nxe+1)*(nze+1)*iq+is*(nxe+1) + ip
   num[6] = (nxe+1)*(nze+1)*iq+(is-1)*(nxe+1) + ip
   num[7] = num[6] + 1
   num[8] = num[5] + 1

   coord[1:2,1] = x_coords[ip]
   coord[5:6,1] = x_coords[ip]
   coord[3:4,1] = x_coords[ip+1]
   coord[7:8,1] = x_coords[ip+1]

   coord[1:4,2] = y_coords[iq]
   coord[5:8,2] = y_coords[iq+1]

   coord[2:3,3] = z_coords[is]
   coord[6:7,3] = z_coords[is]
   coord[1:3:4,3] = z_coords[is+1]
   coord[5:3:8,3] = z_coords[is+1]

  elseif nod == 14
   fac1=(2*nxe+1)*(2*nze+1)*(iq-1)
   fac2=(2*nxe+1)*(2*nze+1)*iq
   num[1] = fac1+is*(2*nxe+1) + ip
   num[2] = num[1] - (2*nxe+1)
   num[3] = num[2] + 1
   num[4] = num[1] + 1
   num[5] = num[2] + nxe + 1
   num[6] = fac1+(nxe+1)*(nze+1) + nxe*nze+(is-1)*(2*nxe+1) + nxe+ip
   num[7] = num[6] - nxe
   num[8] = num[6] + 1
   num[9] = num[8] + nxe
   num[10] = fac2+is*(2*nxe+1) + ip
   num[11] = num[10] - (2*nxe+1)
   num[12] = num[11] + 1
   num[13] = num[10] + 1
   num[14] = num[11] + nxe + 1

   coord[1:2,1] = x_coords[ip]
   coord[6,1] = x_coords[ip]
   coord[10:11,1] = x_coords[ip]
   coord[5:9:2,1] = 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[14,1] = 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[3:4,1] = x_coords[ip+1]
   coord[8,1] = x_coords[ip+1]
   coord[12:13,1] = x_coords[ip+1]

   coord[1:5,2] = y_coords[iq]
   coord[6:9,2] = 0.5*(y_coords[iq] + y_coords[iq+1])
   coord[10:14,2] = y_coords[iq+1]

   coord[2:3,3] = z_coords[is]
   coord[7,3] = z_coords[is]
   coord[11:12,3] = z_coords[is]
   coord[5:6,3] = 0.5*(z_coords[is] + z_coords[is+1])
   coord[8:6:14,3] = 0.5*(z_coords[is] + z_coords[is+1])
   coord[1:4:3,3] = z_coords[is+1]
   coord[9,3] = z_coords[is+1]
   coord[10:3:13,3] = z_coords[is+1]

  elseif nod == 20
   fac1=((2*nxe+1)*(nze+1) + (2*nze+1)*(nxe+1))*(iq-1)
   fac2=((2*nxe+1)*(nze+1) + (2*nze+1)*(nxe+1))*iq
   num[1] = fac1+(3*nxe+2)*is+2*ip-1
   num[2] = fac1+(3*nxe+2)*is-nxe+ip-1
   num[3] = num[1] - 3*nxe-2
   num[4] = num[3] + 1
   num[5] = num[4] + 1
   num[6] = num[2] + 1
   num[7] = num[1] + 2
   num[8] = num[1] + 1
   num[9] = fac2-(nxe+1)*(nze+1) + (nxe+1)*is+ip
   num[10] = num[9]-nxe-1
   num[11] = num[10] + 1
   num[12] = num[9] + 1
   num[13] = fac2+(3*nxe+2)*is+2*ip-1
   num[14] = fac2+(3*nxe+2)*is-nxe+ip-1
   num[15] = num[13] - 3*nxe-2
   num[16] = num[15] + 1
   num[17] = num[16] + 1
   num[18] = num[14] + 1
   num[19] = num[13] + 2
   num[20] = num[13] + 1

   coord[1:3,1] .= x_coords[ip]
   coord[9:10,1] .= x_coords[ip]
   coord[13:15,1] .= x_coords[ip]
   coord[5:7,1] .= x_coords[ip+1]
   coord[11:12,1] .= x_coords[ip+1]
   coord[17:19,1] .= x_coords[ip+1]
   coord[4:4:8,1] .= 0.5*(x_coords[ip] + x_coords[ip+1])
   coord[16:4:20,1] .= 0.5*(x_coords[ip] + x_coords[ip+1])

   coord[1:8,2] .= y_coords[iq]
   coord[13:20,2] .= y_coords[iq+1]
   coord[9:12,2] .= 0.5*(y_coords[iq] + y_coords[iq+1])

   coord[1,3] = z_coords[is+1]
   coord[7:9,3] .= z_coords[is+1]
   coord[12:13,3] .= z_coords[is+1]
   coord[19:20,3] .= z_coords[is+1]
   coord[3:5,3] .= z_coords[is]
   coord[10:11,3] .= z_coords[is]
   coord[15:17,3] .= z_coords[is]
   coord[2:4:6,3] .= 0.5*(z_coords[is] + z_coords[is+1])
   coord[14:4:18,3] .= 0.5*(z_coords[is] + z_coords[is+1])

  else
    println("Wrong number of nodes for hexahedral fin_el.")
    exit(1)
  end
end


text:
include(joinpath("implicit_trapezoid", "integrator.jl"))
include(joinpath("implicit_trapezoid", "body.jl"))
# include(joinpath("implicit_trapezoid", "joint.jl"))
include(joinpath("implicit_trapezoid", "bound.jl"))
include(joinpath("implicit_trapezoid", "test.jl"))