text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Implementation of transformations (generics, translations, rotations, ...)

"""
    Transformation

A wrapper around two 4x4 matrices representing a transformation for [`Vec`](@ref), [`Normal`](@ref), and [`Point`](@ref) instances.

A 4x4 matrix is needed to use the properties of homogeneous coordinates in 3D space. Storing the inverse of the transformation
significantly increases performance at the cost of memory space.

# Fields

- `m::SMatrix{4, 4, Float32}`: the homogeneous matrix representation of the transformation.
- `invm::SMatrix{4, 4, Float32}`: the homogeneous matrix representation of the inverse transformation.
"""
struct Transformation
    m::SMatrix{4, 4, Float32}
    invm::SMatrix{4, 4, Float32}

    function Transformation(m::SMatrix{4, 4, Float32} = SMatrix{4, 4, Float32}(I(4)),
                            invm::SMatrix{4, 4, Float32} = inv(m))
        new(m, invm)
    end
end

@doc """
    Transformation(m::SMatrix{4, 4, Float32} = SMatrix{4, 4, Float32}(I(4)),
                   invm::SMatrix{4, 4, Float32} = inv(m))

Constructor for a [`Transformation`](@ref) instance.

If no parameter is specified, then an identity transformation is returned.
If only the direct matrix is specified, then the inverse matrix is automatically computed.

# Examples

```jldoctest
julia> Transformation()
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  1.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  1.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  1.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  1.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
```
""" Transformation(::SMatrix{4, 4, Float32}, invm::SMatrix{4, 4, Float32})

"""
    Transformation(m::AbstractMatrix)
    Transformation(m::AbstractMatrix, invm::AbstractMatrix)

Construct a [`Transformation`](@ref) instance from `m` and `invm`. The elements of the matrix will be casted to `Float32`.

If any argument is an `AbstractMatrix`, it will be implicitly casted to a `StaticArrays.SMatrix` to increase performance.

# Examples

```jldoctest; setup = :(import StaticArrays)
julia> Transformation(StaticArrays.SMatrix{4,4}([1 0 0 0; 0 2 0 0; 0 0 4 0; 0 0 0 1]))
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  4.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0   0.0f0
 0.0f0  0.5f0  0.0f0   0.0f0
 0.0f0  0.0f0  0.25f0  0.0f0
 0.0f0  0.0f0  0.0f0   1.0f0

julia> Transformation([1 0 0 0; 0 2 0 0; 0 0 4 0; 0 0 0 1])
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  4.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0   0.0f0
 0.0f0  0.5f0  0.0f0   0.0f0
 0.0f0  0.0f0  0.25f0  0.0f0
 0.0f0  0.0f0  0.0f0   1.0f0
```

```jldoctest; setup = :(import LinearAlgebra)
julia> Transformation(LinearAlgebra.Diagonal([1,2,4,1]))
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  4.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0   0.0f0
 0.0f0  0.5f0  0.0f0   0.0f0
 0.0f0  0.0f0  0.25f0  0.0f0
 0.0f0  0.0f0  0.0f0   1.0f0
```
"""
Transformation(m::AbstractMatrix) = Transformation(SMatrix{4, 4, Float32}(m))
Transformation(m::AbstractMatrix, invm::AbstractMatrix) = Transformation(SMatrix{4, 4, Float32}(m), SMatrix{4, 4, Float32}(invm))

################
# Miscellaneous

function show(io::IO, ::MIME"text/plain", t::Transformation)
    println(io, "4x4 $(typeof(t)):")
    println(io, "Matrix of type ", typeof(t.m), ":")
    print_matrix(io, t.m)
    println(io, "\nInverse matrix of type ", typeof(t.invm), ":")
    print_matrix(io, t.invm)
end

eltype(::Transformation) = Float32
eltype(::Type{Transformation}) = Float32

"""
    isconsistent(t::Transformation)

Return `true` if `t.m * t.invm` is similar to the identity matrix and so the [`Transformation`](@ref) is consistent.

Mainly used for testing and to verify matrices haven't been mutated.
"""
isconsistent(t::Transformation) = (t.m * t.invm) ≈ I(4)

(≈)(t1::Transformation, t2::Transformation) = t1.m ≈ t2.m && t1.invm ≈ t2.invm

(*)(t1::Transformation, t2::Transformation) = Transformation(t1.m * t2.m, t2.invm * t1.invm)
(*)(t ::Transformation, v ::Vec) = Vec(@view(t.m[1:3,1:3]) * v)
(*)(t ::Transformation, n ::Normal) = Normal(transpose(@view(t.invm[1:3,1:3])) * n)
function (*)(t ::Transformation, p ::Point)
    res = t.m * SVector(p.v..., one(eltype(p)))
    res = res[end] == 1 ? res : res/res[end]
    Point(@view(res[1:3]))
end

##########
# Inverse

"""
    inv(t::Transformation)

Return the inverse [`Transformation`](@ref).

Returns a `Transformation` which has the `m` and `invm` fields swapped.

# Examples

```jldoctest; setup = :(using LinearAlgebra: Diagonal)
julia> t = Transformation(Diagonal([1, 2, 3, 1]))
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  3.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0         0.0f0
 0.0f0  0.5f0  0.0f0         0.0f0
 0.0f0  0.0f0  0.33333334f0  0.0f0
 0.0f0  0.0f0  0.0f0         1.0f0

julia> inv(t)
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0         0.0f0
 0.0f0  0.5f0  0.0f0         0.0f0
 0.0f0  0.0f0  0.33333334f0  0.0f0
 0.0f0  0.0f0  0.0f0         1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  3.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
 ```
"""
inv(t::Transformation) = Transformation(t.invm, t.m)

############
# Rotations

let rotation_matrices = Dict(
        :X => :(@SMatrix(Float32[   1      0      0     0;
                                    0    cos(θ) -sin(θ) 0;
                                    0    sin(θ)  cos(θ) 0;
                                    0      0      0     1])),
        :Y => :(@SMatrix(Float32[ cos(θ)   0     sin(θ) 0;
                                    0      1      0     0;
                                 -sin(θ)   0     cos(θ) 0;
                                    0      0      0     1])),
        :Z => :(@SMatrix(Float32[ cos(θ) -sin(θ)  0     0;
                                  sin(θ)  cos(θ)  0     0;
                                    0      0      1     0;
                                    0      0      0     1]))
    )

    for (ax, mat) ∈ pairs(rotation_matrices)
        quote
            function $(Symbol(:rotation, ax))(θ::Real)
                m = $mat
                Transformation(m, transpose(m))
            end
        end |> eval
    end

    # docstrings
    let docmsg = (ax, mat) -> """
            rotation$ax(θ::Real)

        Return a [`Transformation`](@ref) that rotates a 3D vector field of the given angle in radians around the $ax-axis.

        If an `AbstractVector` is provided as argument it must have a `size = (3,)`.

        # Examples
        ```jldoctest
        julia> rotation$ax(π/4)
        $(replace(repr(MIME("text/plain"), mat), "Raytracer." => "" ))
        ```
        """
        @doc docmsg(:X, rotationX(π/4)) rotationX
        @doc docmsg(:Y, rotationY(π/4)) rotationY
        @doc docmsg(:Z, rotationZ(π/4)) rotationZ
    end
end

###############
# Translations

"""
    translation(v::AbstractVector)
    translation(x::Real, y::Real, z::Real)

Return a [`Transformation`](@ref) that translates a 3D vector field of the given coordinates.

If an `AbstractVector` is provided as argument it must have a `size = (3,)`.

# Examples

```jldoctest
julia> translation(1, 2, 3)
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  1.0f0
 0.0f0  1.0f0  0.0f0  2.0f0
 0.0f0  0.0f0  1.0f0  3.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  -1.0f0
 0.0f0  1.0f0  0.0f0  -2.0f0
 0.0f0  0.0f0  1.0f0  -3.0f0
 0.0f0  0.0f0  0.0f0   1.0f0
```

```jldoctest
julia> translation([1, 2, 3])
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  1.0f0
 0.0f0  1.0f0  0.0f0  2.0f0
 0.0f0  0.0f0  1.0f0  3.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  -1.0f0
 0.0f0  1.0f0  0.0f0  -2.0f0
 0.0f0  0.0f0  1.0f0  -3.0f0
 0.0f0  0.0f0  0.0f0   1.0f0
```
"""
function translation(v::AbstractVector)
    size(v) == (3,) || raise(ArgumentError("argument 'v' has size = $(size(v)) but 'translate' requires an argument of size = (3,)"))

    mat = Diagonal(ones(eltype(v), 4)) |> MMatrix{4, 4}
    mat⁻¹ = copy(mat)
    mat[1:3, end]   =  v
    mat⁻¹[1:3, end] = -v
    Transformation(mat, mat⁻¹)
end

translation(x::Real, y::Real, z::Real) = translation([x,y,z])

##########
# Scaling

"""
    scaling(x::Real, y::Real, z::Real)
    scaling(s::Real)
    scaling(v::AbstractVector)

Return a [`Transformation`](@ref) that scales a 3D vector field of a given factor for each axis.

If a single `Real` is provided as argument then the scaling is considered uniform.
If an `AbstractVector` is provided as argument it must have a `size = (3,)`.

# Examples

```jldoctest
julia> scaling(1, 2, 3)
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  3.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0         0.0f0
 0.0f0  0.5f0  0.0f0         0.0f0
 0.0f0  0.0f0  0.33333334f0  0.0f0
 0.0f0  0.0f0  0.0f0         1.0f0
```

```jldoctest
julia> scaling(2)
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 2.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  2.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 0.5f0  0.0f0  0.0f0  0.0f0
 0.0f0  0.5f0  0.0f0  0.0f0
 0.0f0  0.0f0  0.5f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
```

```jldoctest
julia> scaling([1, 2, 3])
4x4 Transformation:
Matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0  0.0f0
 0.0f0  2.0f0  0.0f0  0.0f0
 0.0f0  0.0f0  3.0f0  0.0f0
 0.0f0  0.0f0  0.0f0  1.0f0
Inverse matrix of type StaticArrays.SMatrix{4, 4, Float32, 16}:
 1.0f0  0.0f0  0.0f0         0.0f0
 0.0f0  0.5f0  0.0f0         0.0f0
 0.0f0  0.0f0  0.33333334f0  0.0f0
 0.0f0  0.0f0  0.0f0         1.0f0
```
"""
function scaling(x::Real, y::Real, z::Real)
    Transformation(Diagonal(        [x, y, z, true]),
                   Diagonal(true ./ [x, y, z, true])) # NOTE: the use of true is to avoid unwanted promotions
end

scaling(s::Real) = scaling(s, s, s)

function scaling(v::AbstractVector)
    size(v) == (3,) || raise(ArgumentError("argument 'v' has size = $(size(v)) but 'scaling' requires an argument of size = (3,)"))
    scaling(v...)
end


text:
function update(c::Client, o::Overwrite, ch::DiscordChannel; kwargs...)
    return edit_channel_permissions(c, ch.id, o.id; kwargs..., type=o.type)
end

function delete(c::Client, o::Overwrite, ch::DiscordChannel)
    return delete_channel_permission(c, ch.id, o.id)
end

ext: jl


text:
function score(player::Player, board::Board)
    cards::Array{Card} = _assemble_cards(player.hole, board)
    _sort_cards!(cards)
    best_hand = _royal_flush_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.ROYAL, best_hand...)
    end
    best_hand = _straight_flush_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.SFLUSH, best_hand...)
    end
    best_hand = _quad_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.QUAD, best_hand...)
    end
    best_hand = _full_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.FULL, best_hand...)
    end
    best_hand = _flush_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.FLUSH, best_hand...)
    end
    best_hand = _straight_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.STRAIGHT, best_hand...)
    end
    best_hand = _trip_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.TRIP, best_hand...)
    end
    best_hand = _dubs_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.DUBS, best_hand...)
    end
    best_hand = _pair_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.PAIR, best_hand...)
    end
    best_hand = _high_draw(cards)
    if best_hand != nothing
        return Hand(player.playerID, Rankings.HIGH, best_hand...)
    end
    return nothing
end

function compare_hands(hands::Hand...)
    if length(hands) < 1
        return nothing
    end
    if length(hands) == 1
        return hands
    end
    best_hands::Vector{Hand} = [hands[1]]
    i = 2
    while i <= length(hands)
        k = compare(hands[i], best_hands[1])
        if k > 0
            best_hands = [hands[i]]
        elseif k == 0
            push!(best_hands, hands[i])
        end
        i += 1
    end
    return best_hands
end

function _assemble_cards(hole::Hole, board::Board)
    cards = [
        hole.card1,
        hole.card2,
        board.flop1,
        board.flop2,
        board.flop3,
        board.turn,
        board.river
    ]
    return cards
end

function _sort_cards!(cards::Array{Card})
    i = 2
    while i <= 7
        j = i
        while j > 1
            if compare(cards[j], cards[j-1], tie_break=true) > 0
                tmp = cards[j]
                cards[j] = cards[j-1]
                cards[j-1] = tmp
                j -= 1
            else
                break
            end
        end
        i += 1
    end
end

function _royal_flush_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    straight_flush = _straight_flush_draw(cards)
    if straight_flush == nothing
        return nothing
    end
    if straight_flush[1].rank == Ranks.ACE
        return straight_flush
    end
    return nothing
end

function _straight_flush_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    clubs::Array{Card} = []
    diamonds::Array{Card} = []
    hearts::Array{Card} = []
    spades::Array{Card} = []
    for card in cards
        if card.suit == Suits.CLUBS
            push!(clubs, card)
        elseif card.suit == Suits.DIAMONDS
            push!(diamonds, card)
        elseif card.suit == Suits.HEARTS
            push!(hearts, card)
        elseif card.suit == Suits.SPADES
            push!(spades, card)
        end
    end
    if length(spades) >= 5
        straight = _straight_draw(spades)
        if straight != nothing
            return straight
        end
    end
    if length(hearts) >= 5
        straight = _straight_draw(hearts)
        if straight != nothing
            return straight
        end
    end
    if length(diamonds) >= 5
        straight = _straight_draw(diamonds)
        if straight != nothing
            return straight
        end
    end
    if length(clubs) >= 5
        straight = _straight_draw(clubs)
        if straight != nothing
            return straight
        end
    end
    return nothing

end

function _quad_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+3 <= length(cards)
        if cards[i].rank == cards[i+1].rank && cards[i+1].rank == cards[i+2].rank && cards[i+2].rank == cards[i+3].rank
            if i == 1
                return cards[1:5]
            else
                return vcat(cards[i:(i+3)], cards[1:1])
            end
        end
        i += 1
    end
    return nothing
end

function _full_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+2 <= length(cards)
        if cards[i].rank == cards[i+1].rank && cards[i+1].rank == cards[i+2].rank
            j = 1
            while j+1 < i
                if cards[j].rank == cards[j+1].rank
                    return vcat(cards[i:(i+2)], cards[j:(j+1)])
                end
                j += 1
            end
            j = i+3
            while j+1 <= length(cards)
                if cards[j].rank == cards[j+1].rank
                    return vcat(cards[i:(i+2)], cards[j:(j+1)])
                end
                j += 1
            end
        end
        i += 1
    end
    return nothing
end

function _flush_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    clubs::Array{Card} = []
    diamonds::Array{Card} = []
    hearts::Array{Card} = []
    spades::Array{Card} = []
    for card in cards
        if card.suit == Suits.CLUBS
            push!(clubs, card)
        elseif card.suit == Suits.DIAMONDS
            push!(diamonds, card)
        elseif card.suit == Suits.HEARTS
            push!(hearts, card)
        elseif card.suit == Suits.SPADES
            push!(spades, card)
        end
    end
    if length(spades) >= 5
        return spades[1:5]
    end
    if length(hearts) >= 5
        return hearts[1:5]
    end
    if length(diamonds) >= 5
        return diamonds[1:5]
    end
    if length(clubs) >= 5
        return clubs[1:5]
    end
    return nothing
end

function _straight_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+4 <= length(cards)
        straight::Array{Card} = [cards[i]]
        j = i+1
        while length(straight) < 5 && j <= length(cards)
            if Int(straight[length(straight)].rank)-1 == Int(cards[j].rank)
                push!(straight, cards[j])
            end
            j += 1
        end
        if length(straight) >= 5
            return straight[1:5]
        end
        i += 1
    end
    return nothing
end

function _trip_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+2 <= length(cards)
        if cards[i].rank == cards[i+1].rank && cards[i+1].rank == cards[i+2].rank
            if i == 1
                return cards[1:5]
            elseif i == 2
                return vcat(cards[2:4], cards[1:1], cards[5:5])
            else
                return vcat(cards[i:(i+2)], cards[1:2])
            end
        end
        i += 1
    end
    return nothing
end

function _dubs_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+1 <= length(cards)
        if cards[i].rank == cards[i+1].rank
            j = i+2
            while j+1 <= length(cards)
                if cards[j].rank == cards[j+1].rank
                    if i == 1
                        if j == 3
                            return cards[1:5]
                        else
                            return vcat(cards[1:2], cards[j:(j+1)], cards[3:3])
                        end
                    else
                        return vcat(cards[i:(i+1)], cards[j:(j+1)], cards[1:1])
                    end
                end
                j += 1
            end
        end
        i += 1
    end
    return nothing
end

function _pair_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    i = 1
    while i+1 <= length(cards)
        if cards[i].rank == cards[i+1].rank
            if i == 1
                return cards[1:5]
            elseif i == 2
                return vcat(cards[2:3], cards[1:1], cards[4:5])
            elseif i == 3
                return vcat(cards[3:4], cards[1:2], cards[5:5])
            else
                return vcat(cards[i:(i+1)], cards[1:3])
            end
        end
        i += 1
    end
    return nothing
end

function _high_draw(cards::Array{Card})
    if length(cards) < 5
        return nothing
    end
    return cards[1:5]
end

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
using GeoInterface
using Documenter

DocMeta.setdocmeta!(GeoInterface, :DocTestSetup, :(using GeoInterface); recursive=true)
cp(joinpath(@__DIR__, "../INTEGRATIONS.md"), joinpath(@__DIR__, "src/reference/integrations.md"); force=true)

makedocs(;
    modules=[GeoInterface],
    authors="JuliaGeo and contributors",
    repo="https://github.com/JuliaGeo/GeoInterface.jl/blob/{commit}{path}#{line}",
    sitename="GeoInterface.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliageo.github.io/GeoInterface.jl",
        assets=String[]
    ),
    pages=[
        "Home" => "index.md",
        "Background" => Any[
            "Simple Features"=>"background/sf.md",
            "History"=>"background/history.md",
        ],
        "Tutorials" => Any[
            "Installation"=>"tutorials/installation.md",
            "Usage"=>"tutorials/usage.md",
        ],
        "Guides" => Any[
            "For developers"=>"guides/developer.md",
            "Defaults"=>"guides/defaults.md",
        ],
        "Reference" => Any[
            "API" => "reference/api.md"
            "Implementations" => "reference/integrations.md"
        ],
    ],
    doctest=true
)

deploydocs(;
    repo="github.com/JuliaGeo/GeoInterface.jl",
    devbranch="main"
)


text:
function lufac(a::Matrix{Float64})
  upper = zeros(size(a))
  upper[1,:] = a[1,:]
  lower = Matrix{Float64}(I, size(a,1), size(a,1))
  n = size(a, 1)
  for k in 1:n-1
    if abs(upper[k,k] > 1e-10)
      for i in k+1:n
        # L component
        for j in 1:i-1
          total = 0.0
          for l in 1:j-1
            total = total - lower[i, l] * upper[l, j]
          end
          lower[i, j] = (a[i, j] + total)/upper[j, j]
        end
        # U component
        for j in 1:n
          total = 0.0
          for l in 1:i-1
            total = total - lower[i, l] * upper[l, j]
          end
          upper[i, j] = a[i, j] +total
        end
      end
    else
      println("Zero pivot found in row $k")
    end
  end
  (lower, upper)
end


text:
using Documenter, SSDPClient

makedocs(sitename = "SSDPClient Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
using DNC
using Test

in_Sn(vec) = sum(vec) == 1.0
in_Δn(vec) = sum(vec) >= 0.0 && sum(vec) <= 1.0

@testset "DNC.jl" begin
    @testset "utils.jl" begin include("utils_tests.jl") end
    @testset "addressing.jl" begin include("addressing_tests.jl") end
    @testset "access.jl" begin include("access_tests.jl") end
    @testset "computer.jl" begin include("computer_tests.jl") end
    @testset "Gradients" begin include("gradient_tests.jl") end
end


text:
using BasisChange
using LinearAlgebra
using Test

δ(x,y) = (x == y ? 1. : 0.)

@testset "BasisChange.jl" begin
    @testset "Proper Transformation" begin
        for N in 3:2:11 # even dimensions require minuses to keep handedness
            B = I
            for i in 1:N # couldn't be bothered to figure the formula out
                B′ = [δ(j-1, (N + k - i)% N) for (k, j) in Iterators.product(1:N,1:N)]
                R = rotation_matrix(B,B′)
                Δ = det(R)
                @test abs(Δ) ≈ 1. atol = 1e-15
                @test sign(Δ) == 1
            end
        end
    end
end


text:
using TimeSeries
using Dates

dates  = collect(DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime("1/1/2024  23:00:00", "d/m/y  H:M:S"))

nodes14= [
                Bus(1 , "Bus 1"  , "REF" ,0.0, 1.06  , (min=0.94, max=1.06), 69),
                Bus(2 , "Bus 2"  , "PV" , -0.08691739674931762,  1.045 , (min=0.94, max=1.06), 69),
                Bus(3 , "Bus 3"  , "PV" , -0.22200588085367873,  1.01  , (min=0.94, max=1.06), 69),
                Bus(4 , "Bus 4"  , "PQ" , -0.18029251173101424,  1.019 , (min=0.94, max=1.06), 69),
                Bus(5 , "Bus 5"  , "PQ" , -0.15323990832510212,  1.02  , (min=0.94, max=1.06), 69),
                Bus(6 , "Bus 6"  , "PV" , -0.24818581963359368,  1.07  , (min=0.94, max=1.06), 13.8),
                Bus(7 , "Bus 7"  , "PQ" , -0.23335052099164186,  1.062 , (min=0.94, max=1.06), 13.8),
                Bus(8 , "Bus 8"  , "PV" , -0.2331759880664424,  1.09  , (min=0.94, max=1.06), 18),
                Bus(9 , "Bus 9"  , "PQ" , -0.2607521902479528,  1.056 , (min=0.94, max=1.06), 13.8),
                Bus(10, "Bus 10" , "PQ" , -0.26354471705114374,  1.051 , (min=0.94, max=1.06), 13.8),
                Bus(11, "Bus 11" , "PQ" , -0.2581341963699613,  1.057 , (min=0.94, max=1.06), 13.8),
                Bus(12, "Bus 12" , "PQ" , -0.2630211182755455,  1.055 , (min=0.94, max=1.06), 13.8),
                Bus(13, "Bus 13" , "PQ" , -0.2645919146023404,  1.05  , (min=0.94, max=1.06), 13.8),
                Bus(14, "Bus 14" , "PQ" , -0.27995081201989047,  1.036 , (min=0.94, max=1.06), 13.8)
            ]
branches14_dc = [
                Line("Line1",  true, 0.0, 0.0, Arc(from=nodes14[1],to=nodes14[2]), 0.01938, 0.05917, (from=0.0264, to=0.0264), 18.046, 1.04),
                Line("Line2",  true, 0.0, 0.0, Arc(from=nodes14[1],to=nodes14[5]),   0.05403, 0.22304, (from=0.0246, to=0.0246), 4.896, 1.04),
                HVDCLine("DCLine3", true, 0.0, Arc(from=nodes14[2],to=nodes14[3]), (min = -600.0, max = 600), (min = -600.0, max = 600), (min = -600.0, max = 600), (min = -600.0, max = 600), (l0 = 0.01, l1 = 0.001)),
                HVDCLine("DCLine4", true, 0.0, Arc(from=nodes14[2],to=nodes14[4]), (min = -600.0, max = 600), (min = -600.0, max = 600), (min = -600.0, max = 600), (min = -600.0, max = 600), (l0 = 0.01, l1 = 0.001)),
                #Line("Line3",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[3]),   0.04699, 0.19797, (from=0.0219, to=0.0219), 5.522, 1.04),
                #Line("Line4",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[4]),   0.05811, 0.17632, (from=0.017,  to=0.017), 6.052, 1.04),
                Line("Line5",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[5]),   0.05695, 0.17388, (from=0.0173, to=0.0173), 6.140, 1.04),
                Line("Line6",  true, 0.0, 0.0, Arc(from=nodes14[3],to=nodes14[4]),   0.06701, 0.17103, (from=0.0064, to=0.0064), 6.116, 1.04),
                Line("Line7",  true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[5]),   0.01335, 0.04211, (from=0.0, to=0.0), 25.434, 1.04),
                TapTransformer("Trans3", true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[7]),  0.0    , 0.20912,  0.0, 0.978,  20.0),
                TapTransformer("Trans1", true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[9]),  0.0    , 0.55618,  0.0, 0.969,  20.0),
                TapTransformer("Trans2", true, 0.0, 0.0, Arc(from=nodes14[5],to=nodes14[6]),  0.0    , 0.25202,  0.0, 0.932,  20.0),
                Line("Line8",  true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[11]),  0.09498, 0.19890, (from=0.0, to=0.0), 5.373, 1.04),
                Line("Line9",  true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[12]),  0.12291, 0.25581, (from=0.0, to=0.0), 2.020, 1.04),
                Line("Line10", true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[13]),  0.06615, 0.13027, (from=0.0, to=0.0), 4.458, 1.04),
                Transformer2W("Trans4", true, 0.0, 0.0, Arc(from=nodes14[7],to=nodes14[8]),  0.0      , 0.17615,  0.0,    20.0),
                Line("Line16", true, 0.0, 0.0, Arc(from=nodes14[7],to=nodes14[9]),   0.0,     0.11001, (from=0.0, to=0.0), 12.444, 1.04),
                Line("Line11", true, 0.0, 0.0, Arc(from=nodes14[9],to=nodes14[10]),  0.03181, 0.08450, (from=0.0, to=0.0), 5.097, 1.04),
                Line("Line12", true, 0.0, 0.0, Arc(from=nodes14[9],to=nodes14[14]),  0.12711, 0.27038, (from=0.0, to=0.0), 3.959, 1.04),
                Line("Line13", true, 0.0, 0.0, Arc(from=nodes14[10],to=nodes14[11]), 0.08205, 0.19207, (from=0.0, to=0.0), 7.690, 1.04),
                Line("Line14", true, 0.0, 0.0, Arc(from=nodes14[12],to=nodes14[13]), 0.22092, 0.19988, (from=0.0, to=0.0), 6.378, 1.04),
                Line("Line15", true, 0.0, 0.0, Arc(from=nodes14[13],to=nodes14[14]), 0.17093, 0.34802, (from=0.0, to=0.0), 10.213, 1.04)
            ]

branches14 = [
                Line("Line1",  true, 0.0, 0.0, Arc(from=nodes14[1],to=nodes14[2]),   0.01938, 0.05917, (from=0.0264, to=0.0264), 18.046, 1.04),
                Line("Line2",  true, 0.0, 0.0, Arc(from=nodes14[1],to=nodes14[5]),   0.05403, 0.22304, (from=0.0246, to=0.0246), 4.896, 1.04),
                Line("Line3",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[3]),   0.04699, 0.19797, (from=0.0219, to=0.0219), 5.522, 1.04),
                Line("Line4",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[4]),   0.05811, 0.17632, (from=0.017,  to=0.017), 6.052, 1.04),
                Line("Line5",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[5]),   0.05695, 0.17388, (from=0.0173, to=0.0173), 6.140, 1.04),
                Line("Line6",  true, 0.0, 0.0, Arc(from=nodes14[3],to=nodes14[4]),   0.06701, 0.17103, (from=0.0064, to=0.0064), 6.116, 1.04),
                Line("Line7",  true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[5]),   0.01335, 0.04211, (from=0.0, to=0.0), 25.434, 1.04),
                TapTransformer("Trans3", true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[7]),  0.0    , 0.20912,  0.0, 0.978,  20.0),
                TapTransformer("Trans1", true, 0.0, 0.0, Arc(from=nodes14[4],to=nodes14[9]),  0.0    , 0.55618,  0.0, 0.969,  20.0),
                TapTransformer("Trans2", true, 0.0, 0.0, Arc(from=nodes14[5],to=nodes14[6]),  0.0    , 0.25202,  0.0, 0.932,  20.0),
                Line("Line8",  true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[11]),  0.09498, 0.19890, (from=0.0, to=0.0), 5.373, 1.04),
                Line("Line9",  true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[12]),  0.12291, 0.25581, (from=0.0, to=0.0), 2.020, 1.04),
                Line("Line10", true, 0.0, 0.0, Arc(from=nodes14[6],to=nodes14[13]),  0.06615, 0.13027, (from=0.0, to=0.0), 4.458, 1.04),
                Transformer2W("Trans4", true, 0.0, 0.0, Arc(from=nodes14[7],to=nodes14[8]),  0.0      , 0.17615,  0.0,    20.0),
                Line("Line16", true, 0.0, 0.0, Arc(from=nodes14[7],to=nodes14[9]),   0.0,     0.11001, (from=0.0, to=0.0), 12.444, 1.04),
                Line("Line11", true, 0.0, 0.0, Arc(from=nodes14[9],to=nodes14[10]),  0.03181, 0.08450, (from=0.0, to=0.0), 5.097, 1.04),
                Line("Line12", true, 0.0, 0.0, Arc(from=nodes14[9],to=nodes14[14]),  0.12711, 0.27038, (from=0.0, to=0.0), 3.959, 1.04),
                Line("Line13", true, 0.0, 0.0, Arc(from=nodes14[10],to=nodes14[11]), 0.08205, 0.19207, (from=0.0, to=0.0), 7.690, 1.04),
                Line("Line14", true, 0.0, 0.0, Arc(from=nodes14[12],to=nodes14[13]), 0.22092, 0.19988, (from=0.0, to=0.0), 6.378, 1.04),
                Line("Line15", true, 0.0, 0.0, Arc(from=nodes14[13],to=nodes14[14]), 0.17093, 0.34802, (from=0.0, to=0.0), 10.213, 1.04)
            ]

thermal_generators14 = [ThermalStandard("Bus1", true, nodes14[1], 2.0, -0.169,
                    TechThermal(2.324, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=3.332), (min=0.0, max=0.1), nothing, nothing),
                    ThreePartCost((430.292599, 2000.0), 0.0, 0.0, 0.0)
                    ),
                    ThermalStandard("Bus2", true, nodes14[2], 0.40, 0.42,
                    TechThermal(1.4, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=1.40), (min=-0.4, max=0.5), nothing, nothing),
                    ThreePartCost((2500.0, 2000.0), 0.0, 0.0, 0.0)
                    ),
                    ThermalStandard("Bus3", true, nodes14[3], 0.0, 0.23,
                    TechThermal(1.0, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=1.0), (min=0.0, max=0.4), nothing, nothing),
                    ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0)
                    ),
                    ThermalStandard("Bus6", true, nodes14[6], 0.0, 0.12,
                    TechThermal(1.0, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=1.0), (min=-0.06, max=0.24), nothing, nothing),
                    ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0)
                    ),
                    ThermalStandard("Bus8", true, nodes14[8], 0.0, 0.174,
                    TechThermal(1.0, PowerSystems.ST, PowerSystems.COAL, (min=0.0, max=1.0), (min=-0.06, max=0.24), nothing, nothing),
                    ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0)
                    )
            ];

loadz1_ts = [ 0.792729978
            0.723201574
            0.710952098
            0.677672816
            0.668249175
            0.67166919
            0.687608809
            0.711821241
            0.756320618
            0.7984057
            0.827836527
            0.840362459
            0.84511032
            0.834592803
            0.822949221
            0.816941743
            0.824079963
            0.905735139
            0.989967048
            1
            0.991227765
            0.960842114
            0.921465115
            0.837001437 ]

loadz2_ts = [ 0.831093782
              0.689863228
              0.666058513
              0.627033103
              0.624901388
              0.62858924
              0.650734211
              0.683424321
              0.750876413
              0.828347191
              0.884248576
              0.888523615
              0.87752169
              0.847534405
              0.8227661
              0.803809323
              0.813282799
              0.907575962
              0.98679848
              1
              0.990489904
              0.952520972
              0.906611479
              0.824307054]

loadz3_ts = [ 0.871297342
              0.670489749
              0.642812243
              0.630092987
              0.652991383
              0.671971681
              0.716278493
              0.770885833
              0.810075243
              0.85562361
              0.892440566
              0.910660449
              0.922135467
              0.898416969
              0.879816542
              0.896390855
              0.978598576
              0.96523761
              1
              0.969626503
              0.901212601
              0.81894251
              0.771004923
              0.717847996]

loads14 = [PowerLoad("Bus2", true, nodes14[2], PowerSystems.ConstantPower, 0.217, 0.127, 0.217, 0.127),
          PowerLoad("Bus3", true, nodes14[3], PowerSystems.ConstantPower, 0.942, 0.19, 0.942, 0.19),
          PowerLoad("Bus4", true, nodes14[4], PowerSystems.ConstantPower, 0.478, -0.039, 0.478, -0.039),
          PowerLoad("Bus5", true, nodes14[5], PowerSystems.ConstantPower, 0.076, 0.016, 0.076, 0.016),
          PowerLoad("Bus6", true, nodes14[6], PowerSystems.ConstantPower, 0.112, 0.075, 0.112, 0.075),
          PowerLoad("Bus9", true, nodes14[9], PowerSystems.ConstantPower, 0.295, 0.166, 0.295, 0.166),
          PowerLoad("Bus10", true, nodes14[10], PowerSystems.ConstantPower, 0.09, 0.058, 0.09, 0.058),
          PowerLoad("Bus11", true, nodes14[11], PowerSystems.ConstantPower, 0.035, 0.018, 0.035, 0.018),
          PowerLoad("Bus12", true, nodes14[12], PowerSystems.ConstantPower, 0.061, 0.016, 0.061, 0.016),
          PowerLoad("Bus13", true, nodes14[13], PowerSystems.ConstantPower, 0.135, 0.058, 0.135, 0.058),
          PowerLoad("Bus14", true, nodes14[14], PowerSystems.ConstantPower, 0.149, 0.050, 0.149, 0.050)]

forecast_DA14 = [Deterministic(loads14[1], "scalingfactor", TimeArray(dates, loadz1_ts)),
               Deterministic(loads14[2], "scalingfactor", TimeArray(dates, loadz1_ts)),
               Deterministic(loads14[3], "scalingfactor", TimeArray(dates, loadz3_ts)),
               Deterministic(loads14[4], "scalingfactor", TimeArray(dates, loadz1_ts)),
               Deterministic(loads14[5], "scalingfactor", TimeArray(dates, loadz2_ts)),
               Deterministic(loads14[6], "scalingfactor", TimeArray(dates, loadz3_ts)),
               Deterministic(loads14[7], "scalingfactor", TimeArray(dates, loadz2_ts)),
               Deterministic(loads14[8], "scalingfactor", TimeArray(dates, loadz2_ts)),
               Deterministic(loads14[9], "scalingfactor", TimeArray(dates, loadz2_ts)),
               Deterministic(loads14[10], "scalingfactor", TimeArray(dates, loadz2_ts)),
               Deterministic(loads14[11], "scalingfactor", TimeArray(dates, loadz2_ts))
];

forecasts14 = Dict{Symbol,Vector{<:Forecast}}(:DA=>forecast_DA14);

battery14 = [GenericBattery(name = "Bat",
primemover = PowerSystems.BA,
available = true,
bus = nodes14[1],
energy = 5.0,
capacity = (min = 5.0, max = 100.0),
rating = 70,
activepower = 10.0,
inputactivepowerlimits = (min = 0.0, max = 50.0),
outputactivepowerlimits = (min = 0.0, max = 50.0),
reactivepower = 0.0,
reactivepowerlimits = (min = -50.0, max = 50.0),
efficiency = (in = 0.80, out = 0.90)
)]


text:
# Set random seed:
srand(1234)

#uniform.jl

"""
Create a NxM uniformly distributed matrix using the method rand()*(b-a)-a where a,b are the Unif Parameters.

Arguments
---------
a
    Lower Bound
b
    Upper Bound
N
    Row Dimension
M
    Column Dimension

"""
function uniform(a,b,N,M)
    #Returns a NxM uniformly distributed matrix
    #for random number in range [a,b] we note that it can be generated by: rand() * (b-a) - a
    #where rand() prints real numbers:
    rand(N,M)*(b-a)+a
end

#sigmoid.jl

"""
Return the sigmoid of x:

Arguments
---------
x
    input
h
    hyperparemeter

"""

function sigmoid(x,h)
    return 1/(1+exp(-h*x))
end

#sigmoid_deriv.jl

"""
Return the sigmoid of x:

Arguments
---------
x
    input
h
    hyperparemeter

"""

function sigmoid_deriv(x,h)
    return (1/(1+exp(-h*x)))*(1-1/(1+exp(-h*x)))
end

#Train_Neural_Net.jl

"""
Train the XOR neural net for a given number of epochs:

Arguments
---------
epochs
    number of times the network will foward and back propagate
Layer_1
    weighting for hidden layer
Layer_2
    weighting for output layer
learning_rate
    Learning Rate for Gradient Descent
"""

function Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)

    #Initialise XOR truth values:
    X=[ 0 0 ; 0 1 ; 1 0 ; 1 1] #When selected X[i,:] we have column vectors so need to transpose for matmult
    Y=[ 0;1;1;0]

    MSE=0

    for i=1:epochs #For the weight matrix wij is the weight for input i going to Neuron j:

        s_1=map(node_function,X*Layer_1) #Applies Sigmoid to X*Layer_1 the weighting matrix
        s_2=map(node_function, s_1*Layer_2) #Applies Sigmoid to a2 * second weighting matrix
        direct_error=s_2-Y
        MSE=0.5*sum(direct_error.*direct_error) #element wise squares direct_error then summs up and divides by length

    #Begin updating:
    delta_outer=-1.0*direct_error.*map(node_deriv,s_1*Layer_2) #this is the delta (which is the raw loss element-wise
                                                        #multiplied by z(3) see explanation in readme

    delta_inner=delta_outer*transpose(Layer_2).*map(node_deriv,X*Layer_1)
    # map(sigmoid_deriv,X*Layer_1) is f'(z2)

    #Update the weights:

    Layer_2 +=learning_rate*transpose(s_1)*delta_outer #Outer layer has been updated.
    Layer_1 += learning_rate*transpose(X)*delta_inner

    print(string(MSE,"\r"))

    end
    return Layer_1, Layer_2, MSE
    end

#Train_Neural_Net.jl

"""
Train the XOR neural net for a given mean square error it needs to reach:

Arguments
---------
MSE_Min
    Error the network trains till
Layer_1
    weighting for hidden layer
Layer_2
    weighting for output layer
learning_rate
    Learning Rate for Gradient Descent
"""

function Train_Neural_Net_MSE(MSE_Min,Layer_1,Layer_2,learning_rate,node_function,node_deriv)

    #Initialise XOR truth values:
    X=[ 0 0 ; 0 1 ; 1 0 ; 1 1] #When selected X[i,:] we have column vectors so need to transpose for matmult
    Y=[ 0;1;1;0]

    MSE=MSE_Min+1

    while MSE>=MSE_Min #For the weight matrix wij is the weight for input i going to Neuron j:

        s_1=map(node_function,X*Layer_1) #Applies Sigmoid to X*Layer_1 the weighting matrix
        s_2=map(node_function, s_1*Layer_2) #Applies Sigmoid to a2 * second weighting matrix
        direct_error=s_2-Y
        MSE=0.5*sum(direct_error.*direct_error) #element wise squares direct_error then summs up and divides by length

    #Begin updating:
    delta_outer=-1.0*direct_error.*map(node_deriv,s_1*Layer_2) #this is the delta (which is the raw loss element-wise
                                                        #multiplied by z(3) see explanation in readme

    delta_inner=delta_outer*transpose(Layer_2).*map(node_deriv,X*Layer_1)
    # map(sigmoid_deriv,X*Layer_1) is f'(z2)

    #Update the weights:

    Layer_2 +=learning_rate*transpose(s_1)*delta_outer #Outer layer has been updated.
    Layer_1 += learning_rate*transpose(X)*delta_inner

    print(string(MSE,"\r"))

    end
    return Layer_1, Layer_2, MSE
    end

#XOR_Net.jl

"""
Output the values of a given XOR_Net

Arguments
---------
M
    Input for the XOR
w1
    Layer 1
w2
    Layer 2

"""

function XOR_Net(M,w1,w2,node_function)
    #Where M is a matrix that
    s_1=map(node_function,M*w1) #Applies Sigmoid to X*Layer_1 the weighting matrix
    s_2=map(node_function,s_1*w2)
    return s_2

end

#-----------------------------Training Example--------------------------------#

# #Initialise Layers of Neuron Weights:

# Layer_1=uniform(0,1,2,2) #(In video above W1)
# Layer_2=uniform(0,1,2,1) #Column vector for the outer layer (in video above W2)

# function hyper_curry(h)
#     return (x->sigmoid(x,h))
# end

# function hyper_curry_deriv(h)
#     return (x->sigmoid_deriv(x,h))
# end

# node_function=hyper_curry(1)
# node_deriv=hyper_curry_deriv(1)

# w1, w2, MSE1 =Train_Neural_Net_Loop(1000000,Layer_1,Layer_2,0.01,node_function,node_deriv)

# print(XOR_Net([0 1],w1,w2,node_function))


text:
const FIELDS = Dict(
    :name       => (String,   0, 100),
    :mode       => (Int,    100,   8),
    :uid        => (Int,    108,   8),
    :gid        => (Int,    116,   8),
    :size       => (Int,    124,  12),
    :mtime      => (Int,    136,  12),
    :chksum     => (Int,    148,   8),
    :typeflag   => (Char,   156,   1),
    :linkname   => (String, 157, 100),
    :magic      => (String, 257,   6),
    :version    => (String, 263,   2),
    :uname      => (String, 265,  32),
    :gname      => (String, 297,  32),
    :devmajor   => (Int,    329,   8),
    :devminor   => (Int,    337,   8),
    :prefix     => (String, 345, 155),
)

function field_data(buf::Vector{UInt8}, field::Symbol)
    type, offset, len = FIELDS[field]
    return String(buf[offset .+ (1:len)])
end

function parse_header(buf::Vector{UInt8})
    Dict(f => field_data(buf, f) for f in keys(FIELDS))
end


text:
module DynamicEnergyBudgets
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end DynamicEnergyBudgets

using ConstructionBase,
      DimensionalData,
      DocStringExtensions,
      FieldMetadata,
      FieldDefaults,
      FieldDocTables,
      Flatten,
      Mixers,
      Setfield,
      SimpleRoots,
      Requires,
      Unitful

using Unitful: °C, K, Pa, kPa, MPa, J, kJ, W, L, g, kg, cm, m, s, hr, d, mol, mmol, μmol, σ, R
using Base: tail

import FieldMetadata: @default, @description, @units, @bounds, @logscaled, @flattenable, @plottable, @selectable,
                      default, description, units, bounds, logscaled, flattenable, plottable, selectable

export AbstractAssim,
       AbstractNAssim, NH4_NO3Assim, NAssim, ConstantNAssim,
       AbstractCAssim, C3Photosynthesis,
       KooijmanPhotosynthesis, KooijmanSLAPhotosynthesis, KooijmanWaterPotentialPhotosynthesis,
       KooijmanNH4_NO3Assim, ConstantCAssim

export AbstractSynthesizingUnit, ParallelComplementarySU, MinimumRuleSU, KfamilySU

export AbstractRate, SimpleRate, FZeroRate

export AbstractResorption, LosslessResorption, StructuralLossResorption, DissipativeResorption

export AbstractTemperatureCorrection, TempCorr, TempCorrLower, TempCorrLowerUpper, ParentTardieu

export ManualTemperature

export AbstractScaling, KooijmanArea

export AbstractMaturity, Maturity

export AbstractMaintenance, Maintenance

export AbstractDEBCore, DEBCore

export AbstractGermination, ThresholdGermination

export AbstractProduction, Production

export AbstractScaling, Isomorph, V0morph, V1morph, V1V0morph, Plantmorph

export PassiveTranslocation, LosslessPassiveTranslocation, DissipativePassiveTranslocation

export ActiveTranslocation, LosslessActiveTranslocation, DissipativeActiveTranslocation

export AbstractCatabolism, Catabolism, CatabolismCN, CatabolismCNshared, CatabolismCNE

export AbstractAllometry, Allometry, SqrtAllometry, FixedAllometry

export AbstractParams, Params, AbstractSharedParams, SharedParams

export AbstractVars, Vars, PlottableVars, CarbonVars, NitrogenVars

export AbstractOrgan, Organ, AbstractOrganism, Plant

const FIELDDOCTABLE = FieldDocTable((:Description, :Default, :Bounds),
                                    (description, default, bounds);
                                    truncation=(100, 40, 100))

const DEAD, ALIVE = false, true

# Auto docstrings
@template TYPES =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

# Field metadata columns
@chain columns @udefault_kw @units @bounds @logscaled @description

const BI_XTOL = 1e-10
const BI_MAXITER = 100

include("traits.jl")
include("components/synthesizing_units.jl")
include("components/temperature_correction.jl")
include("components/scaling.jl")
include("components/allometry.jl")
include("components/resorption.jl")
include("components/assimilation.jl")
include("components/germination.jl")
include("components/maturity.jl")
include("components/production.jl")
include("components/translocation.jl")
include("components/core.jl")
include("components/rate.jl")
include("components/catabolism.jl")
include("organism.jl")
include("components/environment.jl")
include("functions.jl")
include("setup.jl")
include("model.jl")

function __init__()
    @require Photosynthesis="5e10c064-2706-53a3-a67d-d473e313a663" include("components/fvcb.jl")
    @require Microclimate="e0eb800d-4a9f-54ae-b0f8-217228d9d7c3" include("components/microclimate.jl")
end

end # module


text:
"""
    Variable(name::String, value :: Vector{<:Real})

Construct a Variable type
"""
struct Variable <: DataItem
    name :: String
    value :: Vector{<:Real}
end

"""
    standardize(v::Variable)

Standardize the values of the Variable.

# Implementation

The Variable values are standardized by substracting the empirical mean
and dividing by the empirical standard deviation. A [`VariableStd`](@ref) type
is returned.

See also [`Variable`](@ref), [`VariableStd`](@ref) and [`reconstruct`](@ref).
"""
function standardize(v::Variable)::VariableStd

    x = v.value

    x̄ = mean(x)
    s = std(x)

    offset = isapprox(x̄,0) ? zero(typeof(x[1])) : x̄
    scale = (isapprox(s,1) || isapprox(s,0.0)) ? one(typeof(x[1])) : s

    z = (x .- offset) ./ scale

    vstd = VariableStd(v.name, z, offset, scale)

    return vstd

end


text:
r5(x) = round(x,5)
Σ½(x) = sqrt(sum(x))

function meshgrid(a,b)
    grid_a = [i for i in a, j in b]
    grid_b = [j for i in a, j in b]
    grid_a, grid_b
end

function meshgrid2(a,b,c)
    grid_a = [i for i in a, j in b, k in c]
    grid_b = [j for i in a, j in b, k in c]
    grid_c = [k for i in a, j in b, k in c]
    grid_a, grid_b, grid_c
end

"""
`C, Σ = get_centers_multi(centers,σvec)`\n
Returns multidimensional centers and covariance matrices by gridding.
´centers´ is a matrix with n_basis×n_signals center coordinates\n
`σvec` is a matrix with the corresponding n_basis×n_signals widths
`C, Σ` are the corresponding multidimensional center locations and shapes, both are n_signals × nbasis^n_signals
"""
function get_centers_multi(centers,σvec)
    (nbasis,Nsignals) = size(centers)
    Nbasis::Int64 = nbasis^Nsignals
    Centers = zeros(Nsignals, Nbasis)
    Σ = zeros(Nsignals, Nbasis)
    v = Nbasis
    h = 1
    for i = 1:Nsignals
        v = convert(Int64,v / nbasis)
        Centers[i,:] = vec(repmat(centers[:,i]',v,h))'
        Σ[i,:] = vec(repmat(σvec[:,i]',v,h))'
        h *= nbasis
    end
    Centers, Σ
end

function RLS!(Θ, y, ϕ, P, λ)
    Pϕ = P*ϕ
    P[:,:] = 1/λ*(P - (Pϕ*Pϕ')./(λ + ϕ'*Pϕ))
    yp = (ϕ'Θ)[1]
    e = y-yp
    Θ += Pϕ*e
    return nothing
end

function RLS(Θ, y, ϕ, P, λ)
    Pϕ = P*ϕ
    P = 1/λ*(P - (Pϕ*Pϕ')./(λ + ϕ'*Pϕ))
    yp = (ϕ'Θ)[1]
    e = y-yp
    Θ = Θ + Pϕ*e
    return Θ, P
end

function kalman(R1,R2,R12,Θ, y, ϕ, P)
    ϕTP = ϕ'P
    K = (P*ϕ+R12)/(R2+ϕTP*ϕ)
    P = P - (P*ϕ+R12)/(R2 + ϕTP*ϕ)*(ϕTP+R12') + R1
    yp = (ϕ'Θ)[1]
    e = y-yp
    Θ = Θ + K*e
    return Θ, P
end

function quadform(a,Q)
    vecdot(a,(Q*a))
end


text:
using DataUtils
using Test

@testset "countmap" begin
    d = Dict(rand(1:1000) => rand(1:5) for i in 1:10);
    cm = DataUtils.countmap(values(d))
    @test isa(cm, Dict)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type AbstractDensitySample end
export AbstractDensitySample

# ToDo: Make DensitySample immutable
mutable struct DensitySample{
    P<:Real,
    T<:Real,
    W<:Real,
    PA<:AbstractVector{P}
} <: AbstractDensitySample
    params::PA
    log_value::T
    weight::W
end

export DensitySample

DensitySample(params::PA, log_value::T, weight::W) where {
    P<:Real,
    T<:Real,
    W<:Real,
    PA<:AbstractVector{P}
} = DensitySample{P,T,W,PA}(params, log_value, weight)

import Base.==
==(A::DensitySample, B::DensitySample) =
    A.params == B.params && A.log_value == B.log_value && A.weight == B.weight

# ToDo: remove?
Base.length(s::DensitySample) = length(s.params)

function Base.similar(s::DensitySample{P,T,W}) where {P,T,W}
    params = oob(s.params)
    log_value = convert(T, NaN)
    weight = zero(W)
    PA = typeof(params)
    DensitySample{P,T,W,PA}(params, log_value, weight)
end

function Base.copy!(dest::DensitySample, src::DensitySample)
    copy!(dest.params, src.params)
    dest.log_value = src.log_value
    dest.weight = src.weight
    dest
end

nparams(s::DensitySample) = length(s)

# ToDo: Make immutable
struct DensitySampleVector{
    P<:Real,T<:AbstractFloat,W<:Real,
    PA<:AbstractArray{P,2},TA<:AbstractArray{T,1},WA<:AbstractArray{W,1}
} <: BATDataVector{DensitySample{P,T,W,Vector{P}}}
    params::PA
    log_value::TA
    weight::WA
end

export DensitySampleVector

DensitySampleVector(params::PA, log_value::TA, weight::WA) where {
    P<:Real,T<:AbstractFloat,W<:Real,
    PA<:AbstractArray{P,2},TA<:AbstractArray{T,1},WA<:AbstractArray{W,1}
} = DensitySampleVector{P,T,W,PA,TA,WA}(params, log_value, weight)

DensitySampleVector{P,T,W}(nparams::Integer) where {P<:Real,T<:AbstractFloat,W<:Real} =
    DensitySampleVector(ElasticArray{P}(nparams, 0), Vector{T}(0), Vector{W}(0))

DensitySampleVector(::Type{S}, nparams::Integer) where {P<:Real,T<:AbstractFloat,W<:Real,S<:DensitySample{P,T,W}} =
    DensitySampleVector{P,T,W}(nparams)

Base.size(xs::DensitySampleVector) = size(xs.log_value)

Base.getindex(xs::DensitySampleVector, i::Integer) =
    DensitySample(view(xs.params, :, i), xs.log_value[i], xs.weight[i])

Base.@propagate_inbounds Base._getindex(l::IndexStyle, xs::DensitySampleVector, idxs::AbstractVector{<:Integer}) =
    DensitySampleVector(xs.params[:, idxs], xs.log_value[idxs], xs.weight[idxs])

Base.IndexStyle(xs::DensitySampleVector) = IndexStyle(xs.log_value)

function Base.push!(xs::DensitySampleVector, x::DensitySample)
    append!(xs.params, x.params)
    push!(xs.log_value, x.log_value)
    push!(xs.weight, x.weight)
    xs
end

function Base.append!(A::DensitySampleVector, B::DensitySampleVector)
    append!(A.params, B.params)
    append!(A.log_value, B.log_value)
    append!(A.weight, B.weight)
    A
end

Base.@propagate_inbounds function Base.view(A::DensitySampleVector, idxs)
    DensitySampleVector(
        view(A.params, :, idxs),
        view(A.log_value, idxs),
        view(A.weight, idxs)
    )
end

function _swap!(A::DensitySampleVector, i_A::SingleArrayIndex, B::DensitySampleVector, i_B::SingleArrayIndex)
    _swap!(view(A.params, :, i_A), view(A.params, :, i_B))  # Memory allocation!
    _swap!(A.log_value, i_A, B.log_value, i_B)
    _swap!(A.weight, i_A, B.weight, i_B)
    A
end

function read_fom_hdf5(input, ::Type{DensitySampleVector})
    DensitySampleVector(
        input["params"][:,:],
        input["log_value"][:],
        input["weight"][:]
    )
end

function write_to_hdf5(output, samples::DensitySampleVector)
    output["params"] = Array(samples.params)
    output["log_value"] = samples.log_value
    output["weight"] = samples.weight
    nothing
end


text:
using LinearAlgebra
using Logging
using Test

LinearAlgebra.BLAS.set_num_threads(parse(Int, ENV["OMP_NUM_THREADS"]))

include("benchmarker/Benchmarker.jl")
include("kernel_invocations_syrk.jl")
include("kernel_invocations_gemm.jl")
include("kernel_invocations_syr2k.jl")
include("transposition.jl")
include("common_subexpressions.jl")
include("solve_linear_systems.jl")
include("matrix_chain.jl")
include("diagonal_elements.jl")
include("index_problems.jl")
include("partitioned_matrices.jl")
include("properties_solve.jl")
include("partial_operand.jl")
include("loop_translation.jl")
include("add_scal.jl")

csv = Benchmarker.CSV(string(ENV["LAMP_OUTPUT_DIR"], "julia_", string(ENV["OMP_NUM_THREADS"])))

n = parse(Int32, ENV["LAMP_N"])
println(Threads.nthreads())

add_scal(n)
properties_solve(n)
kernel_invocations_syrk(n)
kernel_invocations_gemm(n)
kernel_invocations_syr2k(n)
transposition(n)
common_subexpressions(n)
solve_linear_systems(n)
matrix_chain(n)
diagonal_elements(n)
partial_operand(n)
index_problems(n)
loop_translation(n)
partitioned_matrices(n)


text:
@testset "I481" begin
    @test isempty(lintstr("""
    s = Symbol("abc")
    """))

    @test Set([:I481]) ⊆ messageset(lintstr("""
    s = symbol("end")
    """))
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"""
    𝐶

`𝐶` (𝐶 = `\\itC+tab`) is a default clock, ready to use.

# Examples

```jldoctest
julia> using DiscreteEvents

julia> resetClock!(𝐶)
"clock reset to t₀=0.0, sampling rate Δt=0.01."

julia> 𝐶  # default clock
Clock 1: state=:idle, t=0.0, Δt=0.01, prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
const 𝐶 = Clock()

# ----------------------------------------------------
# setting and getting clock parameters
# ----------------------------------------------------

"""
    setUnit!(clk::Clock, new::FreeUnits)

set a clock to a new time unit in `Unitful`. If necessary convert
current clock times to the new unit.

# Arguments
- `clk::Clock`
- `new::FreeUnits`: new is one of `ms`, `s`, `minute` or `hr` or another Unitful
    `Time` unit.

# Examples
```jldoctest
julia> using DiscreteEvents, Unitful

julia> import Unitful: Time, s, minute, hr

julia> c = Clock(t0=60)     # setup a new clock with t0=60
Clock 1: state=:idle, t=60.0, Δt=0.01, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> tau(c)               # current time is 60.0 NoUnits
60.0

julia> setUnit!(c, s)       # set clock unit to Unitful.s
60.0 s

julia> tau(c)               # current time is now 60.0 s
60.0 s

julia> setUnit!(c, minute)  # set clock unit to Unitful.minute
1.0 minute

julia> tau(c)               # current time is now 1.0 minute
1.0 minute

julia> isa(tau(c), Time)
true

julia> uconvert(s, tau(c))  # ... which can be converted to other time units
60.0 s

julia> tau(c).val           # it has a value of 1.0
1.0

julia> c.time               # internal clock time is set to 1.0 (a Float64)
1.0

julia> c.unit               # internal clock unit is set to Unitful.minute
minute
```
"""
function setUnit!(clk::Clock, new::FreeUnits)
    if isa(1new, Time)
        if clk.unit == new
            println("clock is already set to $new")
        elseif clk.unit == NoUnits
            clk.unit = new
        else
            old = clk.unit
            clk.unit = new
            fac = uconvert(new, 1*old).val
            clk.time *= fac
            clk.end_time *= fac
            clk.tev *= fac
            clk.Δt *= fac
            clk.tn *= fac
        end
    else
        clk.unit = NoUnits
    end
    tau(clk)
end

"""
    tau(clk::Clock=𝐶)

Return the current simulation time.

# Examples
```jldoctest
julia> using DiscreteEvents

julia> resetClock!(𝐶)   # reset the default clock
"clock reset to t₀=0.0, sampling rate Δt=0.01."
julia> tau()            # gives the default clock's time
0.0
```
"""
tau(clk::Clock=𝐶) = clk.unit == NoUnits ? clk.time : clk.time*clk.unit

"""
    sync!(clk::Clock, to::Clock=𝐶)

Force a synchronization of two clocks. Change all registered times of
`clk` accordingly. Convert or force clk.unit to to.unit.
"""
function sync!(clk::Clock, to::Clock=𝐶)
    if (clk.unit == NoUnits) | (clk.unit == to.unit)
        fac = 1
    elseif to.unit == NoUnits
        println(stderr, "Warning: deleted time unit without conversion")
        fac = 1
    else
        fac = uconvert(to.unit, 1clk.unit).val
    end
    Δt = to.time - clk.time*fac
    clk.time = clk.time*fac + Δt
    clk.unit = to.unit
    clk.tn  = clk.tn*fac + Δt
    clk.tev  = clk.tev*fac + Δt
    clk.end_time = clk.end_time*fac + Δt
    clk.Δt = to.Δt
    evq = PriorityQueue{DiscreteEvent,Float64}()
    for (ev, t) ∈ pairs(clk.sc.events)
        evq[ev] = t*fac + Δt
    end
    clk.sc.events = evq
    clk
end

"""
    resetClock!(clk::Clock, Δt::T=0.01; t0::U=0; <keyword arguments>) where {T<:Number, U<:Number}

Reset a clock.

# Arguments
- `clk::Clock`
- `Δt::T=0.01`: sample rate
- `t0::Float64=0` or `t0::Time`: start time
- `hard::Bool=true`: time is reset, all scheduled events and sampling are
    deleted. If hard=false, then only time is reset, event and
    sampling times are adjusted accordingly.
- `unit=NoUnits`: the Time unit for the clock after reset. If a `Δt::Time` is
    given, its Time unit goes into the clock Time unit. If only t0::Time is given,
    its Time unit goes into the clock time unit.

# Examples

```jldoctest
julia> using DiscreteEvents, Unitful

julia> import Unitful: s

julia> c = Clock(1s, t0=60s)
Clock 1: state=:idle, t=60.0s, Δt=1.0s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> resetClock!(c)
"clock reset to t₀=0.0, sampling rate Δt=0.01."

julia> c
Clock 1: state=:idle, t=0.0, Δt=0.01, prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
function resetClock!(clk::Clock, Δt::T=0.01;
                t0::U=0, hard::Bool=true, unit=NoUnits) where {T<:Number,U<:Number}
    if (Δt == 0) && !isempty(clk.ac)
        Δt = clk.Δt
    end
    if  isa(1unit, Time)
        Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(Δt, Time)
        unit = Unitful.unit(1Δt)
        Δt = Δt.val
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif isa(t0, Time)
        unit = Unitful.unit(t0)
        t0 = t0.val
    else
        nothing
    end
    if hard
        clk.state = Idle()
        clk.time = t0
        clk.unit = unit
        clk.tn = t0
        clk.tev = t0
        clk.end_time = t0
        clk.evcount = 0
        clk.scount = 0
        clk.Δt = Δt
        clk.sc.events = PriorityQueue{DiscreteEvent,Float64}()
        clk.sc.cevents = DiscreteCond[]
        clk.processes = Dict{Any, Prc}()
        clk.channels = Channel[]
        clk.sc.samples = Sample[]
    else
        sync!(clk, Clock(Δt, t0=t0, unit=unit))
    end
    if threadid() == 1
        foreach(ac->put!(ac.forth, Reset(hard)), clk.ac)
        foreach(ac->take!(ac.back), clk.ac)
        "clock reset to t₀=$(float(t0*unit)), sampling rate Δt=$(float(Δt*unit))."
    end
end

# initialize a clock.
init!(clk::Clock) = step!(clk, clk.state, Init(""))

# Adjust/convert `t` given according to clock settings and return a Float64
function _tadjust(clk::Clock, t::Unitful.Time) :: Float64
    if clk.unit == NoUnits
        println(stderr, "Warning: clock has no time unit, ignoring units")
        return t.val
    else
        return uconvert(clk.unit, t).val
    end
end
_tadjust(clk::Clock, t::T) where {T<:Number} = float(t)
_tadjust(ac::ActiveClock, t::T) where {T<:Number} = _tadjust(ac.clock, t)
_tadjust(rtc::RTClock, t::T) where {T<:Number} = _tadjust(rtc.clock, t)

"""
    sample_time!([clk::Clock], Δt::N) where {N<:Number}

Set the clock's sample rate starting from now (`tau(clk)`).

# Arguments
- `clk::Clock`: if not supplied, set the sample rate on 𝐶,
- `Δt::N`: sample rate, time interval for sampling
"""
function sample_time!(clk::Clock, Δt::T)  where {T<:Number}
    clk.Δt = Δt isa Unitful.Time ? _tadjust(clk, Δt) : Δt
    clk.tn = clk.time + clk.Δt
end
sample_time!(Δt::T) where {T<:Number} = sample_time!(𝐶, Δt)

# Is a Clock busy?
_busy(clk::Clock) = clk.state == Busy()

# set internal clock times for next event or sampling action.
# If sampling rate Δt==0, c.tn is set to 0
# If no events are present, c.tev is set to c.end_time
function _setTimes(clk::Clock)
    clk.tev = isempty(clk.sc.events) ? clk.end_time : _nextevtime(clk)
    clk.tn  = isempty(clk.sc.samples) ? 0.0 : clk.time + clk.Δt
end

# ----------------------------------------------------
# step! transition functions for Clocks
# ----------------------------------------------------

step!(clk::Clock, ::Undefined, ::Init) = ( clk.state = Idle() )

# if uninitialized, initialize and then Step or Run.
function step!(clk::Clock, ::Undefined, σ::Σ) where {Σ<:Union{Step,Run}}
    step!(clk, clk.state, Init(0))
    step!(clk, clk.state, σ)
end

step!(c::Clock, ::Q, ::Step) where {Q<:Union{Idle,Halted}} = _step!(c)

# Run a clock for a time Δt.
function _run!(c::Clock, Δt::Float64)
    c.end_time = c.time + Δt
    _setTimes(c)
    while any(i->(c.time ≤ i < c.end_time), (c.tn, c.tev))
        _step!(c) == 0 || break
        if c.state == Halted()
            return c.end_time
        end
    end
    !isempty(c.sc.events) && (c.tev == c.end_time) && _event!(c)
    !isempty(c.sc.samples) && (c.tn == c.end_time) && _tick!(c)
    c.time = c.end_time
end

# Finish a run at time tend, catch remaining events scheduled for tend.
function _finish!(c::Clock, tend::Float64)
    while (length(c.sc.events) ≥ 1) && (_nextevtime(c) ≤ tend + Base.eps(tend)*10)
        _event!(c)
        # step!(c, c.state, Step())
        tend = nextfloat(tend)
    end
    c.time = c.end_time
end

# ----------------------------------------
# Run a simulation for a given duration.
#
# The duration is given with `Run(duration)`. Call scheduled events and
# evaluate sampling expressions at each tick in that timeframe.
# ----------------------------------------
function step!(clk::Clock, ::Idle, σ::Run)

    function handle_response(ix::Int)
        while true
            token = take!(clk.ac[ix].back)
            if token isa Done
                clk.ac[ix].load += token.t
                break
            elseif token isa Forward
                assign(clk, token.ev, token.id)
            elseif token isa Error
                return nothing
            else
                error("invalid response ac $ix: $token")
            end
        end
    end

    clk.evcount = 0
    clk.scount = 0
    sync = false
    foreach(ac->put!(ac.forth, Start()), clk.ac)
    if length(clk.ac) == 0
        tend = _run!(clk, σ.duration)
    else
        tend = clk.end_time = clk.time + σ.duration
        while clk.time < tend
            Δt = min(clk.Δt, abs(tend-clk.time))       # sync yet missing
            foreach(ac->put!(ac.forth, Run(Δt, sync)), clk.ac) # run pclocks
            _cycle!(clk, Δt, sync)                             # run yourself
            foreach(x->handle_response(x), eachindex(clk.ac))
        end
    end
    if clk.state == Halted()
        return "run! halted with $(clk.evcount) clock events, $(clk.scount) sample steps, simulation time: $(round(clk.time, digits=2))"
    end
    clk.time = clk.end_time
    _finish!(clk, tend)
    foreach(ac->put!(ac.forth, Finish(tend)), clk.ac)
    for ac in clk.ac
        c1, c2 = take!(ac.back).x
        clk.evcount += c1
        # clk.scount  += c2
    end

    "run! finished with $(clk.evcount) clock events, $(clk.scount) sample steps, simulation time: $(clk.time)"
end

# Stop the clock.
function step!(clk::Clock, ::Busy, ::Stop)
    clk.state = Halted()
    "Halted after $(clk.evcount) events, simulation time: $(clk.time)"
end

# Resume a halted clock.
function step!(clk::Clock, ::Halted, ::Resume)
    clk.state = Idle()
    step!(clk, clk.state, Run(clk.end_time - clk.time, false))
end

# catch all step!-function.
function step!(clk::Clock, q::ClockState, σ::ClockEvent)
    println(stderr, "Warning: undefined transition ",
            "$(typeof(clk)), ::$(typeof(q)), ::$(typeof(σ)))\n",
            "maybe, you should resetClock!")
end

# --------------------------------------------------------
# Clock user interface functions
# --------------------------------------------------------

"""
    run!(clk::Clock, duration::N) where {N<:Number}

Run a simulation for a given duration.
"""
function run!(clk::Clock, duration::T) where {T<:Number}
    yield()      # allow eventually processes to start
    duration = duration isa Unitful.Time ? _tadjust(clk, duration) : duration
    step!(clk, clk.state, Run(duration, false))
end

"""
    incr!(clk::Clock)

Take one simulation step, execute the next tick or event.
"""
incr!(clk::Clock) = step!(clk, clk.state, Step())

"""
    stop!(clk::Clock)

Stop a running simulation.
"""
stop!(clk::Clock) = step!(clk, clk.state, Stop())

"""
    resume!(clk::Clock)

Resume a halted simulation.
"""
resume!(clk::Clock) = step!(clk, clk.state, Resume())


text:
abstract type Curve end

struct Patch{C<:Curve}
    curve::C
    n::Int
end

function ncurves(patch::Patch, points)
    nₚ = length(points)
    1 + Int((nₚ - patch.n) / 2)
end

function (patch::Patch)(t, points)
    nc = ncurves(patch, points)
    curve_offset = clamp(Int(t ÷ (1/nc)), 0, nc - 1)
    remap = Scaling(nc) ∘ Translation(-curve_offset/nc)
    t = remap(Point(t))[]
    patch.curve(t, curve_points(patch, points, curve_offset))
end

startindex(patch::Patch, curve_offset) = 1 + (patch.n - 1) * curve_offset

function curve_points(patch::Patch, points, curve_offset=0)
    start = startindex(patch, curve_offset)
    @view points[start:start + patch.n - 1]
end

function Base.split(points, patch::Patch)
    nc = ncurves(patch, points)
    map(offset -> curve_points(patch, points, offset), 0:(nc - 1))
end


text:
# SPD-License-Identifier: MIT

using ArgParse
using JLD2
using MarkovModels

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "lexicon"
            required = true
            help = "input lexicon file"
        "lexiconfsm"
            required = true
            help = "output lexicon in JLD2 format"
    end
    s.description = """
    Build the lexicon FSM. The lexicon input file should be formatted
    as:

    word1 p r o n u n

    word1 p R o N u n

    word2 w o r d
    """
    parse_args(s)
end

function LinearFSM(seq)
    SF = LogSemifield{Float32}
    fsm = VectorFSM{SF}()

    prev = nothing
    for (i, label) in enumerate(seq)
        initweight = i == 1 ? one(SF) : zero(SF)
        finalweight = i == length(seq) ? one(SF) : zero(SF)
        s = addstate!(fsm, label; initweight, finalweight)
        i > 1 && addarc!(fsm, prev, s)
        prev = s
    end

    fsm
end

function main(args)
    pronuns = Dict()
    open(args["lexicon"], "r") do f
        for line in eachline(f)
            tokens = split(line)
            word = tokens[1]
            pronun = tokens[2:end]
            list = get(pronuns, word, [])
            push!(list, pronun)
            pronuns[word] = list
        end
    end

    jldopen(args["lexiconfsm"], "w") do f
        for word in keys(pronuns)
            fsms = [LinearFSM(pronun) for pronun in pronuns[word]]
            fsm = union(fsms...)
            f[word] = fsm |> minimize |> renormalize

        end
    end
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
using Base.Test
# fourierlayer.jl
import WaveletScattering: InputLayer

inputlayer = InputLayer(
        tops = [:signal],
        symbols = [:time, :chunk],
        data = rand(Float32, 256, 2))

@test !Mocha.can_do_bp(inputlayer)
@test !Mocha.has_neuron(inputlayer)
@test !Mocha.has_param(inputlayer)
@test Mocha.is_source(inputlayer)


text:
# To ensure that tests can run on travis we have to do a little
# hackadoodle here. The tests require a license file. We include
# a license file that is only valid for one day (the day when
# change is submitted).
# If there is no valid license file, we default to that file.

using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const FALLBACK_URL = "mosek://solve.mosek.com:30080"

using MosekTools
const optimizer = Mosek.Optimizer(fallback = FALLBACK_URL)
MOI.set(optimizer, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Mosek"
end

@testset "Parameters" begin
    optimizer = Mosek.Optimizer(fallback = FALLBACK_URL)
    @testset "Double Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("INTPNT_CO_TOL_DFEAS"), 1e-7)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-7
        MOI.set(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1e-8)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1e-8
        @testset "with integer value" begin
            MOI.set(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS"), 1)
            @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_INTPNT_CO_TOL_DFEAS")) == 1
        end
    end
    @testset "Integer Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS"), 100)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 100
        MOI.set(optimizer, MOI.RawParameter("INTPNT_MAX_ITERATIONS"), 200)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_INTPNT_MAX_ITERATIONS")) == 200
        @testset "with enum value" begin
            MOI.set(optimizer, MOI.RawParameter("MSK_IPAR_OPTIMIZER"), MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
            @test MOI.get(optimizer, MOI.RawParameter("MSK_IPAR_OPTIMIZER")) == convert(Int32, MosekTools.Mosek.MSK_OPTIMIZER_DUAL_SIMPLEX)
        end
    end
    @testset "String Parameter" begin
        MOI.set(optimizer, MOI.RawParameter("PARAM_WRITE_FILE_NAME"), "foo.txt")
        # Needs https://github.com/JuliaOpt/Mosek.jl/pull/174
        #@test MOI.get(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "foo.txt"
        MOI.set(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME"), "bar.txt")
        #@test MOI.get(optimizer, MOI.RawParameter("MSK_SPAR_PARAM_WRITE_FILE_NAME")) == "bar.txt"
    end
    @testset "TimeLimitSec" begin
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
        MOI.set(optimizer, MOI.TimeLimitSec(), 1.0)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == 1.0
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === 1.0
        MOI.set(optimizer, MOI.TimeLimitSec(), nothing)
        @test MOI.get(optimizer, MOI.RawParameter("MSK_DPAR_OPTIMIZER_MAX_TIME")) == -1
        @test MOI.get(optimizer, MOI.TimeLimitSec()) === nothing
    end
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(optimizer, false)
    @test MOIU.supports_default_copy_to(optimizer, true)
end

const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3)

@testset "Basic" begin
    @testset "Linear" begin
        MOIT.basic_constraint_tests(
            optimizer, config,
            include=[
                (MOI.SingleVariable, MOI.EqualTo{Float64}),
                (MOI.SingleVariable, MOI.LessThan{Float64}),
                (MOI.SingleVariable, MOI.GreaterThan{Float64}),
                (MOI.SingleVariable, MOI.Interval{Float64}),
                (MOI.SingleVariable, MOI.ZeroOne),
                (MOI.SingleVariable, MOI.Integer),
                (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}),
                (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64})
        ])
    end
    @testset "Conic" begin
        MOIT.basic_constraint_tests(
            optimizer, config,
            include=[
                (MOI.VectorOfVariables, MOI.ExponentialCone),
                (MOI.VectorOfVariables, MOI.DualExponentialCone),
                (MOI.VectorOfVariables, MOI.PowerCone{Float64}),
                (MOI.VectorOfVariables, MOI.DualPowerCone{Float64}),
                (MOI.VectorOfVariables, MOI.SecondOrderCone),
                (MOI.VectorOfVariables, MOI.RotatedSecondOrderCone)
        ])
    end
end

const bridged = MOIB.full_bridge_optimizer(optimizer, Float64)

# Mosek errors during `MOI.set` instead of `MOI.get` when there are duplicates.
#@testset "Name" begin
#    MOIT.nametest(bridged)
#end

@testset "Copy" begin
    model = MOIB.full_bridge_optimizer(Mosek.Optimizer(), Float64)
    MOIT.copytest(bridged, model)
end

@testset "Unit" begin
    # Mosek does not support names
    MOIT.unittest(bridged, config, [
        # TODO
        "number_threads",
        # Find objective bound of 0.0 which is lower than 4.0
        "solve_objbound_edge_cases",
        # Cannot put multiple bound sets of the same type on a variable
        "solve_integer_edge_cases",
        # Cannot mix `ZeroOne` with `GreaterThan`/`LessThan`
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end

@testset "Continuous Linear" begin
    MOIT.contlineartest(bridged, config)
end

@testset "Continuous Quadratic" begin
    MOIT.contquadratictest(bridged, config, [
        # Non-convex
        "ncqcp",
        # QuadtoSOC does not work as the matrix is not SDP
        "socp"
    ])
end

@testset "Continuous Conic" begin
    MOIT.contconictest(bridged, config, ["rootdets", "logdet"])
end

@testset "Integer Linear" begin
    MOIT.intlineartest(optimizer, config, ["int2", "indicator1", "indicator2", "indicator3"])
end

# Test that objective and constraint data are copied over correctly when
# a scalar variable is transformed to a matrix one
@testset "SDP $add_before" for add_before in [true, false]
    atol = 1e-4
    rtol = 1e-3
    # Problem SDP1 - sdo1 from MOSEK docs
    # See sdp1 of MOI contconic tests for how to get the analytical solution
    δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    y2 = 1 - ε*δ
    y1 = 1 - √2*y2
    obj = y1 + y2/2
    k = -2*δ/ε
    x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    α = √(3-2obj-4x2)/2
    β = k*α

    MOI.empty!(bridged)

    X = MOI.add_variables(bridged, 6)
    x = MOI.add_variables(bridged, 3)

    cx = MOI.add_constraint(bridged, MOI.VectorOfVariables(x), MOI.SecondOrderCone(3))
    if add_before
        cX = MOI.add_constraint(bridged, MOI.VectorOfVariables(X), MOI.PositiveSemidefiniteConeTriangle(3))
    end

    c1 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 1, 1, 1], [X[1], X[3], X[end], x[1]]), 0.), MOI.EqualTo(1.))
    c2 = MOI.add_constraint(bridged, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]]), 0.), MOI.EqualTo(1/2))

    objXidx = [1:3; 5:6]
    objXcoefs = 2*ones(5)
    MOI.set(bridged, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
    MOI.set(bridged, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    if !add_before
        cX = MOI.add_constraint(bridged, MOI.VectorOfVariables(X), MOI.PositiveSemidefiniteConeTriangle(3))
    end

    @test MOI.get(bridged, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    MOI.optimize!(bridged)
    @test MOI.get(bridged, MOI.TerminationStatus()) == MOI.OPTIMAL

    @test MOI.get(bridged, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(bridged, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(bridged, MOI.ObjectiveValue()) ≈ obj atol=atol rtol=rtol

    Xv = [α^2, α*β, β^2, α^2, α*β, α^2]
    xv = [√2*x2, x2, x2]
    @test MOI.get(bridged, MOI.VariablePrimal(), X) ≈ Xv atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.VariablePrimal(), x) ≈ xv atol=atol rtol=rtol

    @test MOI.get(bridged, MOI.ConstraintPrimal(), cX) ≈ Xv atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintPrimal(), cx) ≈ xv atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintPrimal(), c1) ≈ 1. atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintPrimal(), c2) ≈ .5 atol=atol rtol=rtol

    cX0 = 1+(√2-1)*y2
    cX1 = 1-y2
    cX2 = -y2
    cXv = [cX0, cX1, cX0, cX2, cX1, cX0]
    @test MOI.get(bridged, MOI.ConstraintDual(), cX) ≈ cXv atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintDual(), cx) ≈ [1-y1, -y2, -y2] atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintDual(), c1) ≈ y1 atol=atol rtol=rtol
    @test MOI.get(bridged, MOI.ConstraintDual(), c2) ≈ y2 atol=atol rtol=rtol

    @test MosekTools.Mosek.getnumvar(optimizer.task) == (add_before ? 3 : 9)
end

#include("test_jump.jl")


text:
"""
# SeisRequests

Request seismic data via the FDSN or IRIS web services interfaces.

Currently implemented schemes are the FDSN Web Services specification, and IRIS's
web services specification for data (timeseries).  Further specifications will
be added in the future.

Use one of the following high-level functions to get information on seismic stations,
seismic events and data:
- [`get_events`](@ref)
- [`get_stations`](@ref)
- [`get_data`](@ref)

See the docstrings for more information.

### Servers
Inbuilt servers can be listed using [`server_list`](@ref).  Not all servers
will return all kinds of information; for instance `"ISC"` only holds earthquake
catalogues.

### Specifications

FDSN Web Services request specification:
    http://www.fdsn.org/webservices/FDSN-WS-Specifications-1.1.pdf

Details on IRISWS timeseries can be found at:
    https://service.iris.edu/irisws/timeseries/1/
"""
module SeisRequests

export
    FDSNDataSelect,
    FDSNEvent,
    FDSNStation,
    IRISTimeSeries,
    add_server!,
    get_data,
    get_events,
    get_request,
    get_stations,
    server_list

using Dates

using HTTP
using DataStructures: OrderedDict
using Parameters
import QuakeML
import Seis
using Seis: GeogEvent, GeogStation
import StationXML
using StationXML: FDSNStationXML

"URIs of servers to which requests can be sent by key"
const SERVERS = Dict{String,String}(
    "Geofon" => "http://geofon.gfz-potsdam.de",
    "INGV" => "http://webservices.ingv.it",
    "IPGP" => "http://ws.ipgp.fr",
    "IRIS" => "http://service.iris.edu",
    "ISC" => "http://www.isc.ac.uk",
    "NCEDC" => "http://service.ncedc.org",
    "NEIC" => "http://earthquake.usgs.gov",
    "Orfeus" => "http://www.orfeus-eu.org",
    "SCEDC" => "http://service.scedc.caltech.edu",
    )

const DEFAULT_SERVER = "IRIS"

"Dict describing the common HTTP status codes returned by FDSN services"
const STATUS_CODES = Dict{Int,String}(
    200 => "Successful request, results follow",
    204 => "Request was properly formatted and submitted but no data matches the selection",
    400 => "Bad request due to improper specification, unrecognized parameter, parameter value out of range, etc.",
    401 => "Unauthorized, authentication required",
    403 => "Authentication failed or access blocked to restricted data",
    404 => "Alternate to 204 (set via the ‘nodata’ parameter), normally used for results returned to a web browser.",
    413 => "Request would result in too much data being returned or the request itself is too large, returned error message should include the service limitations in the detailed description. Service limits should also be documented in the service WADL.",
    414 => "Request URI too large",
    500 => "Internal server error",
    503 => "Service temporarily unavailable, used in maintenance and error conditions")
const CODE_SUCCESS = 200
const CODES_FAILURE = (400, 401, 403, 404, 413, 414, 500, 503)

# Types which accept missing values
const MBool = Union{Missing,Bool}
const MDateTime = Union{Missing,DateTime}
const MFloat = Union{Missing,Float64}
const MInt = Union{Missing,Int}
const MString = Union{Missing,String}

"""
    SeisRequest

An abstract type representing any kind of web request for seismic data.

Current subtypes of `SeisRequest`:
- `FDSNRequest`
- `IRISRequest`
"""
abstract type SeisRequest end

include("compat.jl")
include("util.jl")
include("fdsnws.jl")
include("fdsnws-formats.jl")
include("irisws.jl")
include("high-level.jl")
include("io.jl")

"""Default version string for all requests."""
version_string(::SeisRequest) = "1"

"""
    request_uri(r::SeisRequest; server=$(DEFAULT_SERVER)) -> uri

Return a URI for the request `r`, which can then be obtained via HTTP GET.

`server` can either be a URI or one of the available servers.  (See [`server_list`](@ref).)
"""
function request_uri(r::SeisRequest; server=DEFAULT_SERVER)
    uri = base_uri(r, server=server) * "?"
    firstfield = true
    for f in fieldnames(typeof(r))
        v = getfield(r, f)
        if !ismissing(v)
            if firstfield
                uri *= "$(f)=$(v)"
                firstfield = false
            else
                uri = join((uri, "$(f)=$(v)"), "&")
            end
        end
    end
    uri
end

"""
    post_uri(rs::AbstractArray{SeisRequest}; server=$(DEFAULT_SERVER)) -> uri

Return a URI for the requests `rs`, which can then be obtained via HTTP POST.

`server` can either be a URI or one of the available servers.  (See [`server_list`](@ref).)
"""
post_uri(r::SeisRequest; server=DEFAULT_SERVER) = base_uri(r, server=server)

"""
    base_uri(request::SeisRequest, server=$(DEFAULT_SERVER)) -> uri

Return the base URI for the `request`, which will look something like
`http://service.iris.edu/fdsnws/station/1/query`.

Note that the `query` part is included.
"""
function base_uri(request::SeisRequest; server=DEFAULT_SERVER)
    server = server in keys(SERVERS) ? SERVERS[server] : server
    protocol = protocol_string(request)
    service = service_string(request)
    version = version_string(request)
    uri = join((server, protocol, service, version, "query"), "/")
    uri
end

"""
    add_server!(label, uri)

Add a server with key `label` to the global list of servers.

See also: [`server_list`](@ref)
"""
function add_server!(label, uri)
    label ∈ keys(SERVERS) && throw(ArgumentError("A server with label `$label` already exists"))
    SERVERS[label] = uri
end

"""
    server_list() -> servers

Return a list of available servers for use with a `SeisRequest`.

See also: [`add_server!`](@ref)
"""
server_list() = collect(keys(SERVERS))

"""
    get_request(r::SeisRequest; server=$(DEFAULT_SERVER), verbose=true) -> response::HTTP.Message.Response

Return `response`, the result of requesting `r` via the HTTP GET command.  Optionally specify
the server either by URI or one of the available servers.  (See [`server_list`](@ref).)
"""
function get_request(r::SeisRequest; server=DEFAULT_SERVER, verbose=true)
    uri = request_uri(r; server=server)
    response = HTTP.request("GET", uri)
    status_text = STATUS_CODES[response.status]
    if verbose
        @info("Request status: " * status_text)
    else
        response.status in CODES_FAILURE && @warn("Request status: " * status_text)
    end
    response
end

"""
    post_request(requests::AbstractArray{<:SeisRequest}; server=$(DEFAULT_SERVER), verbose=true) -> response::HTTP.Message.Response

Send a set of `requests` to `server` using the HTTP POST method, returning
the `response`.

This requires that all the options apart from channel, start time and
end time are the same for all the `requests`.
"""
function post_request(rs::AbstractArray{T};
        server=DEFAULT_SERVER, verbose=true) where {T<:FDSNRequest}
    if !requests_can_be_posted(rs)
        throw(ArgumentError("all requests must be identical apart from their " *
            "channel, start time and end time"))
    end
    body = post_string(rs)
    uri = post_uri(first(rs); server=server)
    # FIXME: Currently IRIS doesn't support a `nodata` line in POST requests
    #        so do not include that for compatibility and simply warn
    #        if nodata != 204.
    #        Other datacentres *do* allow this.
    if server == "IRIS" || server == SERVERS["IRIS"]
        first(rs).nodata != 204 &&
            @warn("IRIS does not support setting the nodata field in POST " *
                  "requests, but it has a non-default value for this request.")
        body = replace(body, r"^\s*nodata=.*\n"=>"")
    end
    response = HTTP.request("POST", uri, [], body)
    status_text = STATUS_CODES[response.status]
    if verbose
        @info("Request status: " * status_text)
    else
        response.status in CODES_FAILURE && @warn("Request status: " * status_text)
    end
    response
end

post_request(::AbstractArray{T}; kwargs...) where {T<:IRISTimeSeries} =
    throw(ArgumentError("$T requests cannot be posted"))

"List of fields which are written in the body of a POST request"
const POST_FIELDS = (:network, :station, :location, :channel, :starttime, :endtime)

"""
    requests_can_be_posted(requests::AbstractArray{<:SeisRequest}) -> ::Bool

Return `true` if the `requests` have options which are consistent (i.e., they
are all the same apart from the channel code, start time and end time),
and hence can be sent via a POST request.
"""
function requests_can_be_posted(rs::AbstractArray{T}) where {T<:SeisRequest}
    for field in fieldnames(T)
        field in POST_FIELDS && continue
        all(getfield(x, field) === getfield(first(rs), field) for x in rs) ||
            return false
    end
    true
end

"Return the full string to be sent as part of a POST request."
post_string(rs) = join((post_string_header(first(rs)), post_string_body(rs)), '\n')

"Return the header part of a POST request."
function post_string_header(r::T) where {T<:SeisRequest}
    join([String(f)*"="*string(getfield(r, f)) for f in fieldnames(T)
          if f ∉ POST_FIELDS && getfield(r, f) !== missing], '\n')
end

"Return the body part of a POST request for a set of `SeisRequest`s."
function post_string_body(rs::AbstractArray{<:SeisRequest})
    s = ""
    for r in rs
        for field in POST_FIELDS
            if field !== :network
                s *= " "
            end
            val = getfield(r, field)
            if val === missing
                throw(ArgumentError("$field field cannot be missing"))
            end
            if field === :location && val == ""
                val = "--"
            end
            s *= string(val)
        end
        s *= "\n"
    end
    s
end

# FDSNEvents can be POSTed in theory, though I don't know of any
# servers which allow this.
post_string_body(rs::AbstractArray{FDSNEvent}) = ""

end # module


text:
"""
    filter(Y, P)

"Filter" data `Y` with iterator `P`
calling (handling of `nothing`s omitted)
```
(t, y), state = iterate(Y, state)
(s => x) = dyniterate(P, s => x, (observation = t => y,))
```

"""
struct Filtered{T,S} <: DynamicIterator
    Y::S
    P::T
end
#filter(Y, P, ::Nothing) = Filtered(Y, P)
filter(Y, P) = Filtered(Y, P)

function dyniterate(M::Filtered, start::Start)
    v, q = @returnnothing dyniterate(M.Y, nothing)
    u, p = @returnnothing dyniterate(M.P, Condition(start, 0.0), v)
    u, (q, p)
end

function dyniterate(M::Filtered, (q, p)::Tuple)
    v, q = @returnnothing iterate(M.Y, q)
    u, p = @returnnothing dyniterate(M.P, p, v)
    u, (q, p)
end

#=
struct ControlledFilter{R,T,S} <: DynamicIterator
    C::R
    Y::S
    P::T
end
filter(Y, P; control = nothing) = filter(Y, P, control)
filter(Y, P, C) = ControlledFilter(C, Y, P)

function dyniterate(M::ControlledFilter, ::Nothing, (value,)::Value)
    w, r = @returnnothing dyniterate(M.C, nothing, (control = value,))
    v, q = @returnnothing iterate(M.Y)
    u, p = @returnnothing dyniterate(M.P, nothing, (value = value, control = w, observation = v))
    u, (u, r, q, p)
end

function iterate(M::ControlledFilter, (u, r, q, p)::Tuple)
    w, r = @returnnothing iterate(M.C, r, (control = u,))
    v, q = @returnnothing iterate(M.Y, q)
    u, p = @returnnothing dyniterate(M.P, p, (control = w, observation = v))
    u, (u, r, q, p)
end
=#


text:
include("./lib/read_settings.jl")
include("./sdm16_sim2.jl")
S  = readJSON()
# include("./sdm16_sim3.jl")
Sim2.run()
# function run_sum3()
#
#     s=Dict("nchps"=>15,"mu"=>10.0,"sigma"=>0.5, "errRate"=>1,
#            "tvFrom"=>29,"tvTo"=>30,"evFrom"=>28,"evTo"=>29)
#     Sim3.run(s)
# end

# run_sim2()

# module ExpJrnVer
#
# export roc_adwin1, roc_bayeschp1
#
# function roc_adwin1() end
#
# function roc_bayeschp1()
#     #
# end
#
# function roc_cusum1_pccf()
#     #
# end
#
# end # module ExpJrnVer

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
abstract type AbstractModel{T} end

abstract type AbstractParams end

function k(m::M) where {M <: AbstractModel}
    return sum(m.ϕ .!= 0) + sum(m.θ .!= 0) + m.μ !=0 + 1
end

function fit(params::AbstractVector{AbstractParams}, z::AbstractVector; criterium=aicc)
    criteria = []
    models = []
    for param in params
        model = fit(param, z)
        perf = criterium(model, residuals(model, z))
        push!(criteria, perf)
        push!(models, model)
        println("param : $param")
        println("criterium : $perf")
    end
    return models[findmin(criteria)[2]]
end

forecast(model::M) where {T,M <: AbstractModel{T}} = forecast(model, T[])

function residuals!(α::AbstractVector{T}, model::M, z::AbstractVector{T}) where {T,M <: AbstractModel{T}}
    length(α) == length(z) || throw(DimensionMismatch("α should have the same length as z"))
    N = length(z)
    if !isempty(α) α[1] = z[1] - forecast(model) end
    for i in 2:N
        α[i] = z[i] - forecast(model, view(z, 1:i - 1), view(α, 1:i-1))
    end
    return α
end

residuals(model::M, z::AbstractVector{T}) where {T,M <: AbstractModel{T}} = residuals!(Vector{T}(undef, length(z)), model, z)