text:
using Crayons
using StaticArrays

import Base.+,
    Base.-, Base.first, Base.in, Base.intersect, Base.isempty, Base.issubset, Base.union

export SquareSet

export bishopattacks,
    bishopattacksempty,
    filesquares,
    isempty,
    issingleton,
    kingattacks,
    knightattacks,
    onlyfirst,
    pawnattacks,
    pawnshift_n,
    pawnshift_ne,
    pawnshift_nw,
    pawnshift_s,
    pawnshift_sw,
    pawnshift_se,
    pprint,
    ranksquares,
    removefirst,
    rookattacks,
    rookattacksempty,
    queenattacks,
    queenattacksempty,
    shift_e,
    shift_n,
    shift_s,
    shift_w,
    squarecount,
    squares,
    squaresbetween,
    toarray

export SS_EMPTY,
    SS_FILE_A,
    SS_FILE_B,
    SS_FILE_C,
    SS_FILE_D,
    SS_FILE_E,
    SS_FILE_F,
    SS_FILE_G,
    SS_FILE_H,
    SS_RANK_1,
    SS_RANK_2,
    SS_RANK_3,
    SS_RANK_4,
    SS_RANK_5,
    SS_RANK_6,
    SS_RANK_7,
    SS_RANK_8

"""
    SquareSet

A type representing a set of squares on the chess board.

The most common ways of obtaining a square set are:

  * Initializing it with one or more squares, for instance `SquareSet(SQ_D4,
    SQ_E4, SQ_D5, SQ_E5)`.

  * From one of the predefined square set constants, like `SS_FILE_C` (the
    squares on the C file) or `SS_RANK_7` (the squares on the 7th rank).

  * By extracting it from a chess board. See the `Board` type for details
    about this.

  * By performing operations transforming or combining one or more square sets
    to a new square set.

The union or intersection of two sets can be computed by the functions `union`
and `intersect`, or by the corresponding binary operators `∪` and `∩`. The
complement of a square set is denoted by the unary `-` operator. The difference
between two set is obtained by the `setdiff` function or by the binary `-`
operator. Subset relationships can be tested by the `issubset` function or the
binary operator `⊆`.

To add or remove a square to a square set, use the `+` or `-` operators with
the square set as the left operand and the square as the right operand. To test
whether a square set contains a square, use `s in ss` or `s ∈ ss`.
"""
struct SquareSet
    val::UInt64
end

function Base.show(io::IO, ss::SquareSet)
    println(io, "SquareSet:")
    for ri ∈ 1:8
        r = SquareRank(ri)
        for fi ∈ 1:8
            f = SquareFile(fi)
            if Square(f, r) ∈ ss
                print(io, " # ")
            else
                print(io, " - ")
            end
        end
        if ri < 8
            println(io, "")
        end
    end
end

function Base.show(io::IO, ::MIME"text/html", ss::SquareSet)
    print(io, Chess.MIME.html(ss))
end

"""
    SquareSet(ss::Vararg{Square})

Construct a square set with the provided squares.

# Examples

```julia-repl
julia> SquareSet(SQ_A1, SQ_A2, SQ_A3, SQ_A4)
SquareSet:
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 #  -  -  -  -  -  -  -
 #  -  -  -  -  -  -  -
 #  -  -  -  -  -  -  -
 #  -  -  -  -  -  -  -
```
"""
function SquareSet(ss::Vararg{Square})
    result = UInt64(0)
    for s ∈ ss
        result |= UInt64(1) << (s.val - 1)
    end
    SquareSet(result)
end

"""
    SS_EMPTY

An empty square set, containing no squares.
"""
const SS_EMPTY = SquareSet(0x0)

"""
    SS_FILE_A

The square set containing all the squares along the A file.
"""
const SS_FILE_A = SquareSet(0xff)

"""
    SS_FILE_B

The square set containing all the squares along the B file.
"""
const SS_FILE_B = SquareSet(0xff00)

"""
    SS_FILE_C

The square set containing all the squares along the C file.
"""
const SS_FILE_C = SquareSet(0xff0000)

"""
    SS_FILE_D

The square set containing all the squares along the D file.
"""
const SS_FILE_D = SquareSet(0xff000000)

"""
    SS_FILE_E

The square set containing all the squares along the E file.
"""
const SS_FILE_E = SquareSet(0xff00000000)

"""
    SS_FILE_F

The square set containing all the squares along the F file.
"""
const SS_FILE_F = SquareSet(0xff0000000000)

"""
    SS_FILE_G

The square set containing all the squares along the G file.
"""
const SS_FILE_G = SquareSet(0xff000000000000)

"""
    SS_FILE_H

The square set containing all the squares along the H file.
"""
const SS_FILE_H = SquareSet(0xff00000000000000)

"""
    SS_RANK_1

The square set containing all the squares along the 1st rank.
"""
const SS_RANK_1 = SquareSet(0x8080808080808080)

"""
    SS_RANK_2

The square set containing all the squares along the 2nd rank.
"""
const SS_RANK_2 = SquareSet(0x4040404040404040)

"""
    SS_RANK_3

The square set containing all the squares along the 3rd rank.
"""
const SS_RANK_3 = SquareSet(0x2020202020202020)

"""
    SS_RANK_4

The square set containing all the squares along the 4th rank.
"""
const SS_RANK_4 = SquareSet(0x1010101010101010)

"""
    SS_RANK_5

The square set containing all the squares along the 5th rank.
"""
const SS_RANK_5 = SquareSet(0x0808080808080808)

"""
    SS_RANK_6

The square set containing all the squares along the 6th rank.
"""
const SS_RANK_6 = SquareSet(0x0404040404040404)

"""
    SS_RANK_7

The square set containing all the squares along the 7th rank.
"""
const SS_RANK_7 = SquareSet(0x0202020202020202)

"""
    SS_RANK_8

The square set containing all the squares along the 8th rank.
"""
const SS_RANK_8 = SquareSet(0x0101010101010101)

const FILE_SQUARES = SVector(
    SS_FILE_A,
    SS_FILE_B,
    SS_FILE_C,
    SS_FILE_D,
    SS_FILE_E,
    SS_FILE_F,
    SS_FILE_G,
    SS_FILE_H,
)

const RANK_SQUARES = SVector(
    SS_RANK_8,
    SS_RANK_7,
    SS_RANK_6,
    SS_RANK_5,
    SS_RANK_4,
    SS_RANK_3,
    SS_RANK_2,
    SS_RANK_1,
)

"""
    filesquares(f::SquareFile)

The set of all squares on the provided file.

# Examples

```julia-repl
julia> filesquares(FILE_G) == SS_FILE_G
true
```
"""
function filesquares(f::SquareFile)::SquareSet
    FILE_SQUARES[f.val]
end

"""
    ranksquares(r::SquareRank)

The set of all squares on the provided rank.

# Examples

```julia-repl
julia> ranksquares(RANK_2) == SS_RANK_2
true
```
"""
function ranksquares(r::SquareRank)::SquareSet
    RANK_SQUARES[r.val]
end

"""
    isempty(ss::SquareSet)

Determine whether a square set is the empty set.

# Examples

```julia-repl
julia> isempty(SS_RANK_1)
false

julia> isempty(SS_EMPTY)
true

julia> isempty(SS_RANK_1 ∩ SS_RANK_2)
true
```
"""
function isempty(ss::SquareSet)::Bool
    ss == SS_EMPTY
end

"""
    union(ss1::SquareSet, ss2::SquareSet)
    ∪(ss1::SquareSet, ss2::SquareSet)

Compute the union of two square sets.

The binary operator `∪` can be used instead of the named function.

# Examples

```julia-repl
julia> SS_FILE_C ∪ SS_RANK_3
SquareSet:
 -  -  #  -  -  -  -  -
 -  -  #  -  -  -  -  -
 -  -  #  -  -  -  -  -
 -  -  #  -  -  -  -  -
 -  -  #  -  -  -  -  -
 #  #  #  #  #  #  #  #
 -  -  #  -  -  -  -  -
 -  -  #  -  -  -  -  -
```
"""
function union(ss1::SquareSet, ss2::SquareSet)
    SquareSet(ss1.val | ss2.val)
end

"""
    intersect(ss1::SquareSet, ss2::SquareSet)
    ∩(ss1::SquareSet, ss2::SquareSet)

Compute the intersection of two square sets.

The binary operator `∩` can be used instead of the named function.

# Examples

```julia-repl
julia> SS_FILE_D ∩ SS_RANK_7
SquareSet:
 -  -  -  -  -  -  -  -
 -  -  -  #  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
```
"""
function intersect(ss1::SquareSet, ss2::SquareSet)
    SquareSet(ss1.val & ss2.val)
end

"""
    -(ss::SquareSet)

The complement of a square set.

# Examples

```julia-repl
julia> ss = SquareSet(SQ_C4);

julia> SQ_C4 ∈ ss
true

julia> SQ_D4 ∈ ss
false

julia> SQ_C4 ∈ -ss
false

julia> SQ_D4 ∈ -ss
true
```
"""
function -(ss::SquareSet)
    SquareSet(~ss.val)
end

"""
    setdiff(ss1::SquareSet, ss2::SquareSet)
    -(ss1::SquareSet, ss2::SquareSet)

The set of all squares that are in `ss1`, but not in `ss2`.

# Examples

```julia-repl
julia> SquareSet(SQ_A1, SQ_A2, SQ_A3, SQ_B1, SQ_B2, SQ_B3) - SS_RANK_2
SquareSet:
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 #  #  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 #  #  -  -  -  -  -  -
```
"""
function setdiff(ss1::SquareSet, ss2::SquareSet)
    ss1 ∩ -ss2
end, function -(ss1::SquareSet, ss2::SquareSet)
    setdiff(ss1, ss2)
end

"""
    +(ss::SquareSet, s::Square)

Add a square to a square set.

If a square is added that is already a member is added to the set, the set is
returned unchanged.

# Examples

```julia-repl
julia> SquareSet(SQ_A1) + SQ_H8 == SquareSet(SQ_A1, SQ_H8)
true

julia> SS_FILE_A + SQ_A1 == SS_FILE_A
true
```
"""
function +(ss::SquareSet, s::Square)
    SquareSet(ss.val | (UInt64(1) << (s.val - 1)))
end

"""
    -(ss::SquareSet, s::Square)

Remove a square from a square set.

If a non-member square is removed, the set is returned unchanged.

# Examples

```julia-repl
julia> SquareSet(SQ_A1, SQ_B1) - SQ_B1 == SquareSet(SQ_A1)
true

julia> SS_FILE_A - SQ_H8 == SS_FILE_A
true
```
"""
function -(ss::SquareSet, s::Square)
    SquareSet(ss.val & ~(UInt64(1) << (s.val - 1)))
end

"""
    issubset(ss1::SquareSet, ss2::SquareSet)
    ⊆(ss1::SquareSet, ss2::SquareSet)

Determine whether `ss1` is a subset of `ss2`.

# Examples

```julia-repl
julia> SquareSet(SQ_A1, SQ_A2) ⊆ SS_FILE_A
true

julia> SquareSet(SQ_A1, SQ_B1) ⊆ SS_FILE_A
false
```
"""
issubset(ss1::SquareSet, ss2::SquareSet) = ss1 ∩ ss2 == ss1

"""
    in(s::Square, ss::SquareSet)
    ∈(s::Square, ss::SquareSet)

Determine whether a square is a member of a square set.

# Examples

```julia-repl
julia> SQ_D7 ∈ SS_RANK_8
false

julia> SQ_D8 ∈ SS_RANK_8
true
```
"""
function in(s::Square, ss::SquareSet)::Bool
    (ss.val & (UInt64(1) << (s.val - 1))) ≠ 0
end

"""
    toarray(ss::SquareSet, T::Type{<:Number}=Float32)

Convert a square set to a two-dimensional array of the specified element type.

The returned array's columns corresponds to the files, and the rows to the
ranks. The array entries are 1 for the members of the square set, and 0
for non-members.

# Examples

```julia-repl
julia> toarray(SS_FILE_C ∪ SS_RANK_5, Int)
8×8 Array{Int64,2}:
 0  0  1  0  0  0  0  0
 0  0  1  0  0  0  0  0
 0  0  1  0  0  0  0  0
 1  1  1  1  1  1  1  1
 0  0  1  0  0  0  0  0
 0  0  1  0  0  0  0  0
 0  0  1  0  0  0  0  0
 0  0  1  0  0  0  0  0
```
"""
function toarray(ss::SquareSet, T::Type{<:Number} = Float32)::Array{T,2}
    reshape([Square(i) ∈ ss ? one(T) : zero(T) for i ∈ 1:64], 8, 8)
end

"""
    squarecount(ss::SquareSet)

The number of members of a square set.
"""
function squarecount(ss::SquareSet)::Int
    count_ones(ss.val)
end

"""
    first(ss::SquareSet)

The first square in a square set.

Returns SQ_NONE for an empty square set.
"""
function first(ss::SquareSet)::Square
    Square(trailing_zeros(ss.val) + 1)
end

"""
    removefirst(ss::SquareSet)

Remove the first member of a square set.

# Examples

```julia-repl
julia> removefirst(SquareSet(SQ_A4, SQ_D5, SQ_F6)) == SquareSet(SQ_D5, SQ_F6)
true
```
"""
function removefirst(ss::SquareSet)::SquareSet
    SquareSet(ss.val & (ss.val - 1))
end

"""
    issingleton(ss::SquareSet)

Determine whether `ss` contains exactly one square.

# Examples

```julia-repl
julia> issingleton(SquareSet(SQ_D5))
true

julia> issingleton(SquareSet(SQ_D5, SQ_C5))
false

julia> issingleton(SS_EMPTY)
false
```
"""
function issingleton(ss::SquareSet)::Bool
    ss ≠ SS_EMPTY && removefirst(ss) == SS_EMPTY
end

"""
    onlyfirst(ss::SquareSet)

Return a square set with all squares except the first removed.

# Examples

```julia-repl
julia> onlyfirst(SquareSet(SQ_A4, SQ_D5, SQ_F6)) == SquareSet(SQ_A4)
true
```
"""
function onlyfirst(ss::SquareSet)::SquareSet
    SquareSet(ss.val & -ss.val)
end

function Base.iterate(ss::SquareSet, state = ss)
    if isempty(state)
        nothing
    else
        (first(state), removefirst(state))
    end
end

"""
    squares(ss::SquareSet)

Convert a square set to a vector of squares.

# Examples

```julia-repl
julia> tostring.(squares(SS_RANK_1))
8-element Array{String,1}:
 "a1"
 "b1"
 "c1"
 "d1"
 "e1"
 "f1"
 "g1"
 "h1"
```
"""
function squares(ss::SquareSet)::Array{Square,1}
    result = Array{Square,1}()
    for s ∈ ss
        push!(result, s)
    end
    result
end

"""
    shift_n(ss::SquareSet)

Shift the square set one step in the 'north' direction.

Squares that are shifted off the edge of the board disappear.

```julia-repl
julia> shift_n(SS_RANK_2) == SS_RANK_3
true

julia> shift_n(SquareSet(SQ_D3, SQ_E4, SQ_F8)) == SquareSet(SQ_D4, SQ_E5)
true
```
"""
function shift_n(ss::SquareSet)::SquareSet
    SquareSet(ss.val >> 1) - SS_RANK_1
end

"""
    shift_s(ss::SquareSet)

Shift the square set one step in the 'south' direction.

Squares that are shifted off the edge of the board disappear.

# Examples

```julia-repl
julia> shift_s(SS_RANK_3) == SS_RANK_2
true

julia> shift_s(SquareSet(SQ_C3, SQ_D2, SQ_E1)) == SquareSet(SQ_C2, SQ_D1)
true
```
"""
function shift_s(ss::SquareSet)::SquareSet
    (SquareSet)(ss.val << 1) - SS_RANK_8
end

"""
    shift_e(ss::SquareSet)

Shift the square set one step in the 'east' direction.

Squares that are shifted off the edge of the board disappear.

# Examples

```julia-repl
julia> shift_e(SS_FILE_F) == SS_FILE_G
true

julia> shift_e(SquareSet(SQ_F5, SQ_G6, SQ_H7)) == SquareSet(SQ_G5, SQ_H6)
true
```
"""
function shift_e(ss::SquareSet)::SquareSet
    SquareSet(ss.val << 8)
end

"""
    shift_w(ss::SquareSet)

Shift the square set one step in the 'west' direction.

Squares that are shifted off the edge of the board disappear.

```julia-repl
julia> shift_w(SS_FILE_C) == SS_FILE_B
true

julia> shift_w(SquareSet(SQ_C5, SQ_B6, SQ_A7)) == SquareSet(SQ_B5, SQ_A6)
true
```
"""
function shift_w(ss::SquareSet)::SquareSet
    SquareSet(ss.val >> 8)
end

"""
    pawnshift_n(ss::SquareSet)

Shift a square set of pawns one step in the 'north' direction.

This is identical to the `shift_n` function except that `pawnshift_n` is a
little faster, but will not work for square sets containing squares on the 1st
or 8th rank.
"""
function pawnshift_n(ss::SquareSet)::SquareSet
    SquareSet(ss.val >> 1)
end

"""
    pawnshift_s(ss::SquareSet)

Shift a square set of pawns one step in the 'south' direction.

This is identical to the `shift_s` function except that `pawnshift_s` is a
little faster, but will not work for square sets containing squares on the 1st
or 8th rank.
"""
function pawnshift_s(ss::SquareSet)::SquareSet
    SquareSet(ss.val << 1)
end

"""
    pawnshift_nw(ss::SquareSet)

Shift a square set of pawns one step in the 'north west' direction.

This is identical to calling `shift_n` followed by `shift_w`, except that
`pawnshift_nw` is a little faster, but will not work for square sets
containing squares on the 1st or 8th rank.
"""
function pawnshift_nw(ss::SquareSet)::SquareSet
    SquareSet(ss.val >> 9)
end

"""
    pawnshift_ne(ss::SquareSet)

Shift a square set of pawns one step in the 'north east' direction.

This is identical to calling `shift_n` followed by `shift_e`, except that
`pawnshift_ne` is a little faster, but will not work for square sets
containing squares on the 1st or 8th rank.
"""
function pawnshift_ne(ss::SquareSet)::SquareSet
    SquareSet(ss.val << 7)
end

"""
    pawnshift_sw(ss::SquareSet)

Shift a square set of pawns one step in the 'south west' direction.

This is identical to calling `shift_s` followed by `shift_w`, except that
`pawnshift_sw` is a little faster, but will not work for square sets
containing squares on the 1st or 8th rank.
"""
function pawnshift_sw(ss::SquareSet)::SquareSet
    SquareSet(ss.val >> 7)
end

"""
    pawnshift_se(ss::SquareSet)

Shift a square set of pawns one step in the 'south east' direction.

This is identical to calling `shift_s` followed by `shift_e`, except that
`pawnshift_se` is a little faster, but will not work for square sets
containing squares on the 1st or 8th rank.
"""
function pawnshift_se(ss::SquareSet)::SquareSet
    SquareSet(ss.val << 9)
end

"""
    bishopattacks(blockers::SquareSet, s::Square)

The squares attacked by a bishop on `s`, with `blockers` being the occupied
squares.
"""
function bishopattacks(blockers::SquareSet, s::Square)::SquareSet
    SquareSet(bishopattacks(blockers.val, s.val))
end

"""
    rookattacks(blockers::SquareSet, s::Square)

The squares attacked by a rook on `s`, with `blockers` being the occupied
squares.
"""
function rookattacks(blockers::SquareSet, s::Square)::SquareSet
    SquareSet(rookattacks(blockers.val, s.val))
end

"""
    queenattacks(blockers::SquareSet, s::Square)

The squares attacked by a queen on `s`, with `blockers` being the occupied
squares.
"""
function queenattacks(blockers::SquareSet, s::Square)::SquareSet
    bishopattacks(blockers, s) ∪ rookattacks(blockers, s)
end

function computestepattacks(s, deltas)
    foldl((acc, d) -> distance(s, s + d) ≤ 2 ? acc + (s + d) : acc, deltas, init = SS_EMPTY)
end

function computeknightattacks(s::Square)::SquareSet
    computestepattacks(
        s,
        [
            2 * DELTA_N + DELTA_E,
            2 * DELTA_N + DELTA_W,
            DELTA_N + 2 * DELTA_E,
            DELTA_N + 2 * DELTA_W,
            DELTA_S + 2 * DELTA_E,
            DELTA_S + 2 * DELTA_W,
            2 * DELTA_S + DELTA_E,
            2 * DELTA_S + DELTA_W,
        ],
    )
end

function computekingattacks(s::Square)::SquareSet
    computestepattacks(
        s,
        [DELTA_NE, DELTA_N, DELTA_NW, DELTA_E, DELTA_W, DELTA_SE, DELTA_S, DELTA_SW],
    )
end

function computewpattacks(s::Square)::SquareSet
    computestepattacks(s, [DELTA_NE, DELTA_NW])
end

function computebpattacks(s::Square)::SquareSet
    computestepattacks(s, [DELTA_SE, DELTA_SW])
end

const N_ATTACKS = @SVector [computeknightattacks(Square(i)) for i ∈ 1:64]
const K_ATTACKS = @SVector [computekingattacks(Square(i)) for i ∈ 1:64]
const WP_ATTACKS = @SVector [computewpattacks(Square(i)) for i ∈ 1:64]
const BP_ATTACKS = @SVector [computebpattacks(Square(i)) for i ∈ 1:64]

"""
    knightattacks(s::Square)

The set of squares attacked by a knight on the square `s`.
"""
function knightattacks(s::Square)::SquareSet
    @inbounds N_ATTACKS[s.val]
end

"""
    kingattacks(s::square)

the set of squares attacked by a king on the square `s`.
"""
function kingattacks(s::Square)::SquareSet
    @inbounds K_ATTACKS[s.val]
end

"""
    pawnttacks(c::PieceColor, s::square)

the set of squares attacked by a pawn of color `c` on the square `s`.
"""
function pawnattacks(c::PieceColor, s::Square)::SquareSet
    if c == WHITE
        @inbounds WP_ATTACKS[s.val]
    else
        @inbounds BP_ATTACKS[s.val]
    end
end

const B_ATTACKS_EMPTY = @SVector [bishopattacks(SS_EMPTY, Square(i)) for i ∈ 1:64]

const R_ATTACKS_EMPTY = @SVector [rookattacks(SS_EMPTY, Square(i)) for i ∈ 1:64]

const Q_ATTACKS_EMPTY = @SVector [queenattacks(SS_EMPTY, Square(i)) for i ∈ 1:64]

"""
    bishopattacksempty(s::Square)

The set of squares a bishop on `s` would attack on an otherwise empty board.
"""
function bishopattacksempty(s::Square)::SquareSet
    @inbounds B_ATTACKS_EMPTY[s.val]
end

"""
    rookattacksempty(s::Square)

The set of squares a rook on `s` would attack on an otherwise empty board.
"""
function rookattacksempty(s::Square)::SquareSet
    @inbounds R_ATTACKS_EMPTY[s.val]
end

"""
    queenattacksempty(s::Square)

The set of squares a queen on `s` would attack on an otherwise empty board.
"""
function queenattacksempty(s::Square)::SquareSet
    @inbounds Q_ATTACKS_EMPTY[s.val]
end

function computesquaresbetween(s1::Square, s2::Square)::SquareSet
    if s2 ∉ queenattacksempty(s1)
        SS_EMPTY
    else
        result = SS_EMPTY
        si1 = min(s1.val, s2.val)
        si2 = max(s1.val, s2.val)
        for si3 ∈ (si1+1):(si2-1)
            s3 = Square(si3)
            ss = SquareSet(s3)
            if s2 ∈ queenattacksempty(s1) && s2 ∉ queenattacks(ss, s1)
                result += s3
            end
        end
        result
    end
end

const SQUARES_BETWEEN = [computesquaresbetween(Square(i), Square(j)) for i ∈ 1:64, j = 1:64]

"""
    squaresbetween(s1::Square, s2::Square)

The set of squares on the line, file or diagonal between `s1` and `s2`.

When a queen on `s1` would attack `s2` on an otherwise empty board, this
function returns the set of squares where a piece would block the queen
on `s1` from attacking `s2`.

# Examples

```julia-repl
julia> squaresbetween(SQ_A4, SQ_D4) == SquareSet(SQ_B4, SQ_C4)
true

julia> squaresbetween(SQ_F7, SQ_A2)
SquareSet:
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  #  -  -  -
 -  -  -  #  -  -  -  -
 -  -  #  -  -  -  -  -
 -  #  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -
```
"""
function squaresbetween(s1::Square, s2::Square)::SquareSet
    @inbounds SQUARES_BETWEEN[s1.val, s2.val]
end

function colorpprint(ss::SquareSet)
    for ri ∈ 1:8
        r = SquareRank(ri)
        for fi ∈ 1:8
            f = SquareFile(fi)
            bg = (ri + fi) % 2 == 0 ? 0x87CEFF : 0x6CA6CD
            bg = (ri + fi) % 2 == 0 ? 0x8accc0 : 0x66b0a3
            if Square(f, r) ∈ ss
                print(Crayon(background = bg, foreground = 0xff1654), " * ")
            else
                print(Crayon(background = bg), "   ")
            end
        end
        println(Crayon(reset = true), "")
    end
end

"""
    pprint(ss::SquareSet)

Pretty-print a square set to the standard output.
"""
function pprint(ss::SquareSet; color = false)
    if color
        colorpprint(ss)
        return
    end
    for ri ∈ 1:8
        r = SquareRank(ri)
        println("+---+---+---+---+---+---+---+---+")
        for fi ∈ 1:8
            f = SquareFile(fi)
            if Square(f, r) ∈ ss
                print("| # ")
            else
                print("|   ")
            end
        end
        println("|")
    end
    println("+---+---+---+---+---+---+---+---+")
end


text:
module eFEM

# Custom types
include("MeshTypes.jl")
include("ParameterTypes.jl")
include("ProblemTypes.jl")
include("SolutionTypes.jl")

# Mesh input/generation
include("MeshGeneration.jl")
include("GMSHreader.jl")
include("MeshTransform.jl")

# Generate Matrices
include("Basis.jl")
include("BilinearForms.jl")
include("MatrixGeneration.jl")
include("BoundaryConditions.jl")

# Solve Problem
include("Assembly.jl")
include("TimeStepping.jl")
include("Solver.jl")

# PostProcessing
include("vtkExport.jl")
include("ErrorAnalysis.jl")
include("DomainQuadrature.jl")
include("PostProcessing.jl")
include("TracerGenerate.jl")

export getBoundaries

# ParameterTypes.jl
export PoissonParam,
       HeatParam,
       AdvDiffParam,
       DarcyParam,
       FluidParam,
       BrinkmanParam,
       BrinkmanMPParam,
       FluidVarParam

# ProblemTypes.jl
export Problem,
       Dirichlet,
       Neumann,
       Robin,
       Forcing

# MeshGeneration.jl
export squareMesh,
       squareMeshFluid,
       Mesh,
       FluidMesh,
       axisymMesh,
       axisymFluid

# MeshTransform.jl
export pointTransform,
       onSegment,
       doIntersect,
       isInsideDomain,
       meshTransform

# Basis.jl
export derivShape2D,
       shapeEval,
       GaussEdge,
       FEMstats

# MatrixGeneration.jl
export Mass2DMatrix

# TimeStepping.jl
export FEMForwardEuler,
       progressBar,
       sToHMS

# Solver.jl
export solve,
       solve!,
       GenerateSystem,
       ApplyBC!

# vtkExport.jl
export vtksave,
       Path,
       ScalarData,
       ScalarNames,
       VectorData,
       VectorNames,
       SURF_UNSTRUCTURED_TO_VTK

# ErrorAnalysis.jl
export DomainNorm,
       hCalc

# DomainQuadrature.jl
export DomainQuad,
       SurfaceQuad,
       SurfaceFlux,
       mySurface,
       VelocityFlux,
       SurfaceInterp

# PostProcessing.jl
export computeStress,
       computeDerivative,
       computeExtension,
       computeCompressibility

# TracerGenerate.jl
export GenerateTracers,
       Tracer,
       TracerInfo,
       TracerLineSource

end # module

ext: jl


text:
using Test
using Profile

@testset "PProf.jl" begin
    include("PProf.jl")
end

@testset "flamegraphs.jl" begin
    include("flamegraphs.jl")
end

@static if isdefined(Profile, :Allocs)  # PR https://github.com/JuliaLang/julia/pull/42768
    @testset "Allocs.jl" begin
        include("Allocs.jl")
    end
end


text:
function ec_map_init()
    opt = ec_init()
    opt["tooltip"]["trigger"] = "item"
    opt["toolbox"]["feature"]["mark"] = false
    opt["dataRange"] = Dict(
        "calculable" => true,
        "text" => ["High", "Low"],
        "color" => ["orangered","yellow","lightskyblue"]
    )

    return opt
end

function ec_map{N<:Real,S<:String,T<:String}(legendlist::Vector{S}, data::Array{N}, maptype::String, labels::Vector{T})
    opt = ec_map_init()
    opt["legend"]["data"] = legendlist

    data_cumsum = cumsum(data, 2)[:, end]
    opt["dataRange"]["min"] = int(0.95 * minimum(data_cumsum))
    opt["dataRange"]["max"] = int(1.05 * maximum(data_cumsum))

    m, n = size(data)
    resize!(opt["series"], n)
    for i = 1:n
        coldata = Dict(
            "name" => legendlist[i],
            "type" => "map",
            "mapType" => maptype,
            "data" => Array(Any, m)
            )

        for j = 1:m
            coldata["data"][j] = Dict("name" => labels[j], "value" => data[j, i])
        end

        opt["series"][i] = coldata
    end

    return opt
end


text:
function verifyConstraints!(mechanism::Mechanism)
    for eqc in mechanism.eqconstraints
        if norm(g(mechanism, eqc)) > 1e-3
            @info string("Bad constraint satisfaction at constraint: ", eqc.id, ", |g| = ", norm(g(mechanism, eqc)))
        end
    end
    return
end

function activateConstraints!(mechanism::Mechanism)
    graph = mechanism.graph

    for (id,eqc) in pairs(mechanism.eqconstraints)
        activate!(eqc)
        activate!(mechanism.graph,id)
    end
    for (id,ineqc) in pairs(mechanism.ineqconstraints)
        activate!(ineqc)
        activate!(mechanism.graph,id)
    end

    return
end

function deactivateConstraints!(mechanism::Mechanism)
    graph = mechanism.graph

    for (id,eqc) in pairs(mechanism.eqconstraints)
        deactivate!(eqc)
        deactivate!(mechanism.graph,id)
    end
    for (id,ineqc) in pairs(mechanism.ineqconstraints)
        deactivate!(ineqc)
        deactivate!(mechanism.graph,id)
    end

    return
end

function nameiddict(mechanism::Mechanism)
    dict = Dict{String,Int64}()
    for (id,body) in pairs(mechanism.bodies)
        if body.name != ""
            dict[body.name] = id
        end
    end
    for (id,eqc) in pairs(mechanism.eqconstraints)
        if eqc.name != ""
            dict[eqc.name] = id
        end
    end
    for (id,ineqc) in pairs(mechanism.ineqconstraints)
        if ineqc.name != ""
            dict[ineqc.name] = id
        end
    end

    return dict
end

function minimalCoordinates(mechanism::Mechanism)
    keys = mechanism.eqconstraints.keys
    values = Vector{SVector}()

    for eqc in mechanism.eqconstraints
        push!(values, minimalCoordinates(mechanism, eqc))
    end

    return UnitDict(keys, values)

end

function setPosition!(mechanism::Mechanism, dict::UnitDict)
    for (id,eqc) in pairs(mechanism.eqconstraints)
        setPosition!(mechanism, eqc, dict[id])
    end

    return
end

function setVelocity!(mechanism::Mechanism, dict::UnitDict)
    for (id,eqc) in pairs(mechanism.eqconstraints)
        setVelocity!(mechanism, eqc, dict[id])
    end

    return
end

function setForce!(mechanism::Mechanism, dict::UnitDict)
    for (id,eqc) in pairs(mechanism.eqconstraints)
        setForce!(mechanism, eqc, dict[id])
    end

    return
end

@inline function currentasknot!(mechanism::Mechanism)
    foreach(currentasknot!, mechanism.bodies)
    return
end


text:
using Bonsai, Test, StructTypes, HTTP, JSON3, HTTP.Messages
using StructTypes: @Struct
using HTTP: Request
using URIs: URI

@Struct struct Payload
    x
end

@Struct struct PayloadTyped
    x::Int
end

@Struct struct PayloadMissing
    x
    y::Union{Float64, Nothing}
end

@testset "Body" begin
    # constructors
    b = Body("ok")
    @test b.val == "ok"
    @test b.t == String

    b = Body(error="test", message="kwconstructor")
    @test b.t == NamedTuple{(:error, :message), Tuple{String, String}}
    @test b.val  == (error ="test", message="kwconstructor")

    b = Body(x=String, y=Float64)
    @test b.t == NamedTuple{(:x, :y), Tuple{String, Float64}}
    @test isnothing(b.val)

    # reading
    io = IOBuffer(JSON3.write(Payload(10)))
    req = Request()
    req.body = take!(io)
    @test Bonsai.read(req, Body(Payload)) == Payload(10)

    io = IOBuffer(JSON3.write(Payload(10)))
    req = Request()
    req.body = take!(io)
    payload = Bonsai.read(req, Body(PayloadMissing))
    @test payload.x  == 10
    @test isnothing(payload.y)
end

@testset "Query" begin
    # constructors
    @test Query(x=String).val |> isnothing
    @test_throws Exception Query(x="x")

    # reading - requires HTTP master
    req = Request()
    req.url = URI("http://localhost?x=10")
    @test Bonsai.read(req, Query(Payload)) isa Payload
    @test Bonsai.convert_numbers!(Dict{Symbol,Any}(:x => "10"), PayloadTyped)[:x] == 10
    @test Bonsai.read(req, Query(PayloadTyped)) isa PayloadTyped
end

@testset "Header" begin

    req = HTTP.Messages.Request()
    req.headers = [
        "X-Test" => "wagwan"
    ]

    @Struct struct A
        x_test::String
    end

    @Struct struct B
        x_missing::Union{String,Missing}
    end

    @test Bonsai.read(req, Headers(A)).x_test == "wagwan"
    @test Bonsai.read(req, Headers(B)).x_missing |> ismissing

    bad_req = HTTP.Messages.Request()
    bad_req.headers = [
        "X-Something-Else" => "wagwan"
    ]

    @test_throws Exception Bonsai.read(bad_req, Headers(A))
end

@testset "Params" begin
    req = HTTP.Messages.Request()
    # Base.match(app, req) should perform this for us
    req.context[:params] = Dict{Any, Any}(:id => "10")
    @test Bonsai.read(req, Params(id=Int)).id == 10
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
mutable struct Simulation{T,W,FC}
	model::Model{T}
	env::Environment{W,FC}
	con::ContactMethods
	res::ResidualMethods
	rz::Any
	rθ::Any
	ϕ::Any
end

function Simulation(model::Model{T}, env::Environment) where T
	Simulation(model, env, ContactMethods(), ResidualMethods(), zeros(0, 0), zeros(0, 0), q -> nothing)
end

function get_simulation(model::String, env::String, sim_name::String;
		model_variable_name::String = model,
		dynamics_name::String = "dynamics",
		gen_base = true,
		gen_dyn = true,
		approx = false)

	#TODO: assert model exists

	dir_model = joinpath(@__DIR__, "..", "dynamics", model)
	dir_sim   = joinpath(@__DIR__, "..", "simulation", model, sim_name)

	dir_base = joinpath(dir_model, dynamics_name, "base.jld2")
	dir_dyn = joinpath(dir_model, dynamics_name, "dynamics.jld2")
	dir_res = joinpath(dir_sim, "residual.jld2")
	dir_jac = joinpath(dir_sim, "jacobians.jld2")

	model = deepcopy(eval(Symbol(model_variable_name)))
	env = deepcopy(eval(Symbol(env)))
	sim = Simulation(model, env)

	instantiate_base!(sim.model, dir_base)

	instantiate_dynamics!(sim.model, dir_dyn,
		derivs = approx)

	instantiate_residual!(sim,
		dir_res, dir_jac,
		jacobians = (approx ? :approx : :full))

	# signed-distance function
	@variables q[1:sim.model.nq]
	ϕ = ϕ_func(sim.model, sim.env, q)
	sim.ϕ = eval(Symbolics.build_function(ϕ, q, checkbounds=true)[2])

	return sim
end

function z_initialize!(z, model::Model, env::Environment{<:World,LinearizedCone}, q)
	z .= 1.0
	iq2 = index_q2(model, env, quat = false)
	z[iq2] = q
end

function z_initialize!(z::Vector{T}, idx::SVector{nq2,Int}, q::Vector{T}) where {T,nq2}
	z .= 1.0
	z[idx] .= q
	return nothing
end

function z_initialize!(z, model::Model, env::Environment{<:World,NonlinearCone}, q1) #TODO: make allocation free
	iq2 = index_q2(model, env, quat = false)
	ib1 = index_b1(model, env, quat = false)
	iψ1 = index_ψ1(model, env, quat = false)
	iη1 = index_η1(model, env, quat = false)
	is2 = index_s2(model, env, quat = false)

    z .= 0.1
    z[iq2] = q1

	# second-order cone initializations
	z[ib1] .= 0.1 # primal cones: vector part # TODO redundant
	z[iψ1] .= 1.0 # primal cones: scalar part
	z[iη1] .= 0.1 # dual cones: vector part # TODO redundant
	z[is2] .= 1.0 # dual cones: scalar part
	return z
end

function z_initialize!(z, idx::SVector{nq2,Int}, model::Model, env::Environment{<:World,NonlinearCone}, q1) where {T,nq2}#TODO: make allocation free
	@warn "soc z initialization not implemented"
end

function z_warmstart!(z, model::Model, env::Environment{<:World,LinearizedCone}, q, a)
	iq2 = index_q2(model, env, quat = false)
	iort = index_ort(model, env, quat = false)
	isoc = index_soc(model, env, quat = false)
	z[iq2] = q
	for i in eachindex(iort)
		z[iort[i]] .+= a * rand(length(iort[i]))
	end
	for i in eachindex(isoc)
		for j in eachindex(isoc[i])
			z[isoc[i][j]] .+= a * rand(length(isoc[i][j]))
		end
	end
	nothing
end

function z_warmstart!(z, model::Model, env::Environment{<:World,NonlinearCone}, q, a)
	# @warn "warm start for second-order cone not implemented"
	z_initialize!(z, model, env, q)
	nothing
end

function θ_initialize!(θ, model::Model, q0, q1, u, w, μ, h)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	off = 0
    θ[1:nq] = q0
	off += nq
    θ[off .+ (1:nq)] = q1
	off += nq
    θ[off .+ (1:nu)] = u
	off += nu
    θ[off .+ (1:nw)] = w
	off += nw
	θ[off .+ (1:1)] .= μ
	off += 1
    θ[off .+ (1:1)] .= h
end

function E_func(model::Model, env::Environment{<:World,LinearizedCone})
	nc = model.nc
	nb = nc * friction_dim(env)
	SMatrix{nc, nb}(kron(Diagonal(ones(nc)), ones(1, Int(nb / nc))))
end

function residual(model::Model, env::Environment{<:World,LinearizedCone}, z, θ, κ)
	nc = model.nc
	nb = nc * friction_dim(env)
	nf = Int(nb / nc)
	np = dim(env)

	q0, q1, u1, w1, μ, h = unpack_θ(model, θ)
	q2, γ1, b1, ψ1, s1, η1, s2 = unpack_z(model, env, z)

	ϕ = ϕ_func(model, env, q2)

	k = kinematics(model, q2)
	λ1 = contact_forces(model, env, γ1, b1, q2, k)
	Λ1 = transpose(J_func(model, env, q2)) * λ1 #@@@@ maybe need to use J_fast
	vT_stack = velocity_stack(model, env, q1, q2, k, h)
	ψ_stack = transpose(E_func(model, env)) * ψ1

	# @warn "define residual order"
	[model.dyn.d(h, q0, q1, u1, w1, Λ1, q2);
	 s1 - ϕ;
	 η1 - vT_stack - ψ_stack;
	 s2 .- (μ[1] * γ1 .- E_func(model, env) * b1);
	 γ1 .* s1 .- κ[1];
	 b1 .* η1 .- κ[1];
	 ψ1 .* s2 .- κ[1]]
end

function residual(model::Model, env::Environment{<:World,NonlinearCone}, z, θ, κ)
	q0, q1, u1, w1, μ, h = unpack_θ(model, θ)
	q2, γ1, b1, ψ1, s1, η1, s2 = unpack_z(model, env, z)
	ϕ = ϕ_func(model, env, q2)
	k = kinematics(model, q2)
	λ1 = contact_forces(model, env, γ1, b1, q2, k)
	vT = velocity_stack(model, env, q1, q2, k, h)
	nc = model.nc
	nf = friction_dim(env)
	ne = dim(env)

	# @warn "define residual order"
	[dynamics(model, env, h, q0, q1, u1, w1, λ1, q2);
	 s1 - ϕ;
	 vcat([η1[(i - 1) * nf .+ (1:nf)] - vT[(i - 1) * nf .+ (1:nf)] for i = 1:nc]...);
	 s2 - μ[1] * γ1;
	 γ1 .* s1 .- κ[1];
	 vcat([
	 	second_order_cone_product( # [ψ1; η1] ∘ [s2; b1] - [κ; 0...0]
			# [ψ1[i]; b1[(i-1) * nf .+ (1:nf)]],
			[ψ1[i]; η1[(i-1) * nf .+ (1:nf)]],
			# [s2[i]; η1[(i-1) * nf .+ (1:nf)]]
			[s2[i]; b1[(i-1) * nf .+ (1:nf)]]
			) - [κ[1]; zeros(nf)]
		for i = 1:nc]...)
	 ]
end

function ResidualMethods()
	function f()
		error("Not Implemented: use instantiate_residual!")
		return nothing
	end
	return ResidualMethods(fill(f, 3)...)
end


text:
using ExcelWriter
using Test
using DataFrames

@warn("Tests are currently rudimentary and currently only test whether the functions can be called (without running into any errors)")

@testset "Smoke Tests" begin
    adir=mktempdir()
    fi=joinpath(adir,"testfile1.xlsx")
    fi2=joinpath(adir,"testfile2.xlsx")
    fi3=joinpath(adir,"testfile3.xlsx")

    adf=DataFrame(rand(4,4),:auto)
    xlsheet1=ExcelSheet("mySheet",adf)
    xlsheet2=ExcelSheet("myOtherSheet",DataFrame(rand(6,4),:auto))
    xlData=ExcelData([xlsheet1,xlsheet2])

    writeExcel(xlData,fi,true,false)
    if isfile(fi)
        @test true
    end

    #write Matrix or some DF
    writeExcel(adf,fi2,true,false)
    if isfile(fi2)
        @test true
    end

    #write a matrix
    writeExcel(rand(4,4),fi3,true,false)
    if isfile(fi3)
        @test true
    end

    for sz in [20_000,40_000]
        df2=DataFrame(l=["l$(l)" for l=1:sz],n=["n$(n)" for n=1:sz],n2= ["n$(m)" for m=1:sz])

        xlDataLarge=ExcelData([ExcelSheet("Test",deepcopy(df2)),ExcelSheet("Test2",deepcopy(df2)),ExcelSheet("Test3",deepcopy(df2)),ExcelSheet("Test4",deepcopy(df2))])

        fi2=joinpath(adir,"testfileLarge.xlsx")
        @show sz
        @time writeExcel(xlDataLarge,fi2,true,false)
        @time writeExcel(xlDataLarge,fi2,true,false,printfilename = true)
        if isfile(fi2)
            @test true
        end

        isfile(fi2) && rm(fi2)
        @info("Running tests without printing filename:")
        @time writeExcel(xlDataLarge,fi2,true,false,printfilename = false)
        if isfile(fi2)
            @test true
        end

    end
end


text:
"""
    mutable struct Context
A contextual data object.
"""
Base.@kwdef mutable struct Context
    id::Int=new_id()
    Window::GLFW.Window = GLFW.Window(C_NULL)
    Time::Cfloat = Cfloat(0.0f0)
    MouseJustPressed::Vector{Bool} = fill(false, Int(ImGuiMouseButton_COUNT))
    MouseCursors::Vector{GLFW.Cursor} = fill(GLFW.Cursor(C_NULL), Int(ImGuiMouseCursor_COUNT))
    KeyOwnerWindows::Vector{GLFW.Window} = fill(GLFW.Window(C_NULL), 512)
    InstalledCallbacks::Bool = true
    WantUpdateMonitors::Bool = true
    PrevUserCallbackMousebutton::Ptr{Cvoid} = C_NULL
    PrevUserCallbackScroll::Ptr{Cvoid} = C_NULL
    PrevUserCallbackKey::Ptr{Cvoid} = C_NULL
    PrevUserCallbackChar::Ptr{Cvoid} = C_NULL
    PrevUserCallbackMonitor::Ptr{Cvoid} = C_NULL
    GlslVersion = 150
    FontTexture = GLuint(0)
    ShaderHandle = GLuint(0)
    VertHandle = GLuint(0)
    FragHandle = GLuint(0)
    AttribLocationTex = GLint(0)
    AttribLocationProjMtx = GLint(0)
    AttribLocationVtxPos = GLint(0)
    AttribLocationVtxUV = GLint(0)
    AttribLocationVtxColor = GLint(0)
    VboHandle = GLuint(0)
    ElementsHandle = GLuint(0)
    ui_handle = igCreateContext(C_NULL)
    ui_func = ()->nothing
end

Base.show(io::IO, x::Context) = print(io, "Context(id=$(x.id))")

const __IMGUI_CONTEXTS = Dict{Int,Context}()
const __IMGUI_CONTEXT_COUNTER = Threads.Atomic{Int}(0)

reset_counter() = Threads.atomic_sub!(__IMGUI_CONTEXT_COUNTER,  __IMGUI_CONTEXT_COUNTER[])
add_counter() = Threads.atomic_add!(__IMGUI_CONTEXT_COUNTER, 1)
get_counter() = __IMGUI_CONTEXT_COUNTER[]
new_id() = (add_counter(); get_counter();)

"""
    create_context()
Return a GLFW backend contextual data object.
"""
function create_context()
    screen = GLMakie.global_gl_screen()
    window = GLMakie.to_native(screen)

    ctx = Context(; Window=window)

    config_imgui()

    # store the ctx in a global variable to prevent it from being GC-ed
    __IMGUI_CONTEXTS[ctx.id] = ctx

    init(ctx)

    return ctx
end

function config_imgui()
    # enable docking and multi-viewport
    io = igGetIO()
    io.ConfigFlags = unsafe_load(io.ConfigFlags) | ImGuiConfigFlags_DockingEnable
    io.ConfigFlags = unsafe_load(io.ConfigFlags) | ImGuiConfigFlags_ViewportsEnable

    # set style
    igStyleColorsDark(C_NULL)

    # When viewports are enabled we tweak WindowRounding/WindowBg so platform windows can look identical to regular ones.
    style = Ptr{ImGuiStyle}(igGetStyle())
    if unsafe_load(io.ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
        style.WindowRounding = 5.0f0
        col = c_get(style.Colors, ImGuiCol_WindowBg)
        c_set!(style.Colors, ImGuiCol_WindowBg, ImVec4(col.x, col.y, col.z, 1.0f0))
    end
end

"""
    release_context(ctx::Context)
Relese the ctx so it can be GC-ed.
"""
release_context(ctx::Context) = delete!(__IMGUI_CONTEXTS, ctx.id)

get_window(ctx::Context) = ctx.Window


text:
using Test
using Planets
using Statistics

include("test_core_radius.jl")
include("test_habitable_zone.jl")
include("test_stellar_evolution.jl")


text:
module Utils

using LinearAlgebra: norm
using JuLIP: JVecF
using SHIPs: TransformedJacobi,  inv_transform

import Base: rand

function rand_sphere()
   R = randn(JVecF)
   return R / norm(R)
end

function rand_radial(J::TransformedJacobi)
   # uniform sample from [tl, tu]
   x = J.tl + rand() * (J.tu - J.tl)
   # transform back
   return inv_transform(J.trans, x)
end

rand_radial(J::TransformedJacobi, N::Integer) = [ rand_radial(J) for _=1:N ]

rand(J::TransformedJacobi) = rand_radial(J) *  rand_sphere()

rand(J::TransformedJacobi, N::Integer) =  [ rand(J) for _ = 1:N ]

_get_ll(KL, νz) = getfield.(KL[νz.ν], :l)

end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE
g = TypeDBClient

# Scenario: create one database
@when("connection create database: alice") do context
    client = context[:client]
    db = g.create_database(client, "alice")
    @expect db === true
end

@then("connection has database: alice") do context
    result = g.contains_database(context[:client], "alice")
    @expect result === true
    #until beforescenario works for me :-) the following
    delete_all_databases(context[:client])
end

# Scenario: create many databases
@when("connection create databases:") do context
    db_names = context.datatable
    all_db = Bool[]
    for db in db_names
        push!(all_db, g.create_database(context[:client], db[1]))
    end
    @expect all(all_db) === true
end

@then("connection has databases:") do context
    db_names = [db[1] for db in context.datatable]
    server_dbs = [item.name for item in g.get_all_databases(context[:client])]

    @expect (Set(db_names) == Set(server_dbs))

    delete_all_databases(context[:client])
end

# Scenario: create many databases in parallel
@when("connection create databases in parallel:") do context
    db_names = context.datatable
    @sync @async for db in db_names
        g.create_database(context[:client], db[1])
    end
end

# Scenario: delete one database
@when("connection delete database: alice") do context
    g.delete_database(context[:client],"alice")
end

@then("connection does not have database: alice") do context
    @expect g.contains_database(context[:client], "alice") == false
end

# Scenario: connection can delete many databases
@when("connection delete databases:") do context
    db_names = [db[1] for db in context.datatable]
    for db in db_names
        g.delete_database(context[:client], db)
    end
end

@then("connection does not have databases:") do context
    db_names = [db[1] for db in context.datatable]
    db_there = Bool[]
    for db in db_names
        push!(db_there, g.contains_database(context[:client], db))
    end
    @expect all(db_there) !== true
end

# Scenario: delete many databases in parallel
@when("connection delete databases in parallel:") do context
    db_names = [db[1] for db in context.datatable]
    @sync @async for db in db_names
        g.delete_database(context[:client], db)
    end
    @expect isempty(g.get_all_databases(context[:client])) === true
end

# Scenario: delete a database causes open sessions to fail
@when("connection create database: typedb") do context
    g.create_database(context[:client], "typedb")
end

@when("connection open session for database: typedb") do context
    context[:session] = g.CoreSession(context[:client], "typedb" , g.Proto.Session_Type.DATA, request_timout=Inf)
end

@when("connection delete database: typedb") do context
    g.delete_database(context[:client], "typedb")
end

@then("connection does not have database: typedb") do context
    @expect g.contains_database(context[:client], "typedb") === false
end

@then("session open transaction of type; throws exception: write") do context
    try
        transaction(context[:session], g.Proto.Transaction_Type.WRITE)
    catch ex
        @info "if a database isn't there a transaction can't be open"
    end
end

@then("graql define; throws exception containing \"transaction has been closed\"") do context
    define_string = "define person sub entity;"
    try
        g.define(context[:transaction], define_string)
    catch ex
       # @expect occursin("transaction has been closed", ex.message) === true
       @info ex
    end
end

@then("connection delete database; throws exception: typedb") do context
    try
        g.delete_database(context[:client], "typedb")
    catch ex
        @expect (typeof(ex) == g.TypeDBClientException)
        @expect occursin("The database typedb does not exist", string(ex))
    end
end

##############  utility functions ########################
function delete_all_databases(client::g.CoreClient)
    for (_, session) in client.sessions
        close(session)
    end
    all_db = g.get_all_databases(client)
    for db in all_db
        g.delete_database(client, db.name)
    end
end


text:
include("ALPA.jl")

function main(dynetfile, commfile)
    m = Dict{Int,Int}() # initialize nodes' membership
    g = Main.ALPA.Graph() # start with an empty graph
    fin = open(dynetfile, "r")
    fout = open(commfile, "w")
    for line in eachline(fin)
        items = split(line)
        if length(items) == 3
            if items[1] == "+"
                Main.ALPA.lpa_addedge!(g, m, parse.(Int,items[2:3])...)
            end
            if items[1] == "-"
                Main.ALPA.lpa_deleteedge!(g, m, parse.(Int,items[2:3])...)
            end
        end
        if length(items) == 2
            if items[1] == "+"
                Main.ALPA.lpa_addnode!(g, m, parse(Int,items[2]))
            end
            if items[1] == "-"
                Main.ALPA.lpa_deletenode!(g, m, parse(Int,items[2]))
            end
        end
        for node in keys(m)
            print(fout, node, "=>", m[node], " ")
        end
        print(fout, "\n")
    end
    close(fin)
    close(fout)
end

if length(ARGS) == 1
    dynetfile = ARGS[1]
    commfile = splitext(dynetfile)[1]*".comm"
elseif length(ARGS) == 2
    dynetfile = ARGS[1]
    commfile = ARGS[2]
else
    error("ARGS number error")
end
main(dynetfile, commfile)


text:
struct MOLFiniteDifference{G} <: DiffEqBase.AbstractDiscretization
    dxs
    time
    approx_order::Int
    upwind_order::Int
    grid_align::G
end

# Constructors. If no order is specified, both upwind and centered differences will be 2nd order
function MOLFiniteDifference(dxs, time=nothing; approx_order = 2, upwind_order = 1, grid_align=CenterAlignedGrid())

    if approx_order % 2 != 0
        @warn "Discretization approx_order must be even, rounding up to $(approx_order+1)"
    end
    @assert approx_order >= 1 "approx_order must be at least 1"
    @assert upwind_order >= 1 "upwind_order must be at least 1"

    @assert (time isa Num) | (time isa Nothing) "time must be a Num, or Nothing - got $(typeof(time)). See docs for MOLFiniteDifference."

    return MOLFiniteDifference{typeof(grid_align)}(dxs, time, approx_order, upwind_order, grid_align)
end


text:
######### Stan program example  ###########

using Mamba, Stan

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

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

bernoullidata = [
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
  Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
]

monitor = ["theta", "lp__", "accept_stat__"]

stanmodel = Stanmodel(update=1200, thin=2, name="bernoulli", model=bernoullimodel);

println("\nStanmodel that will be used:")
stanmodel |> display
println()
println("Input observed data dictionary:")
bernoullidata |> display
println()

sim1 = stan(stanmodel, bernoullidata, ProjDir, diagnostics=false, CmdStanDir=CMDSTAN_HOME);

## Subset Sampler Output to variables suitable for describe().
sim = sim1[1:size(sim1, 1), monitor, 1:size(sim1, 3)]
describe(sim)
println()

## Brooks, Gelman and Rubin Convergence Diagnostic
try
  gelmandiag(sim, mpsrf=true, transform=true) |> display
catch e
  #println(e)
  gelmandiag(sim, mpsrf=false, transform=true) |> display
end
println()

## Geweke Convergence Diagnostic
gewekediag(sim) |> display

## Highest Posterior Density Intervals
hpd(sim) |> display

## Cross-Correlations
cor(sim) |> display

## Lag-Autocorrelations
autocor(sim) |> display

## Deviance Information Criterion
#dic(sim) |> display

## Plotting
p = plot(sim, [:trace, :mean, :density, :autocor], legend=true);
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:svg)
draw(p, ncol=4, filename="$(stanmodel.name)-summaryplot", fmt=:pdf)

# Below will only work on OSX, please adjust for your environment.
# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
        for i in 1:4
          isfile("$(stanmodel.name)-summaryplot-$(i).svg") &&
            run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-summaryplot-$(i).svg"`)
        end
      end : println()

if Pkg.installed("Mamba") > v"0.7.1"
	# Pairwise contour plots
	simmod = ModelChains(sim,Mamba.Model())
	p = plot(simmod, :contour)
	draw(p, nrow=2, ncol=2, filename="$(stanmodel.name)-contourplot.svg")
	draw(p, nrow=2, ncol=2, filename="$(stanmodel.name)-contourplot.pdf", fmt=:pdf)

	# Below will only work on OSX, please adjust for your environment.
	# JULIA_SVG_BROWSER is set from environment variable JULIA_SVG_BROWSER
	@static is_apple() ? if isdefined(Main, :JULIA_SVG_BROWSER) && length(JULIA_SVG_BROWSER) > 0
	        isfile("$(stanmodel.name)-contourplot.svg") &&
	          run(`open -a $(JULIA_SVG_BROWSER) "$(stanmodel.name)-contourplot.svg"`)
	      end : println()
end

end # cd


text:
struct BoundingBoxDetection{T} <: AbstractObjectDetection{T}
    rect::Rect2{T}
    core::ObjectDetectionCore

    function BoundingBoxDetection(top_left::Point2{T}, width::T, height::T, core_args...) where T
        core = ObjectDetectionCore(core_args...)
        new{T}(Rect2{T}(top_left.data[1], top_left.data[2], width, height), core)
    end

    function BoundingBoxDetection(top_left::Point2{T}, bottom_right::Point2{T}, core_args...) where T
        width = bottom_right.data[1] - top_left.data[1]
        height = bottom_right.data[2] - top_left.data[2]
        BoundingBoxDetection(top_left, width, height, core_args...)
    end

    function BoundingBoxDetection(vertices::Vector{Point2{T}}, core::ObjectDetectionCore) where T
        min_x = minimum([x for (x,_) in vertices])
        max_x = maximum([x for (x,_) in vertices])
        min_y = minimum([y for (_,y) in vertices])
        max_y = maximum([y for (_,y) in vertices])

        width = max_x - min_x
        height = max_y - min_y

        new{T}(Rect2{T}(min_x, min_y, width, height), core)
    end

    function BoundingBoxDetection(vertices::Vector{Point2{T}}, core_args...) where T
        core = ObjectDetectionCore(core_args...)
        BoundingBoxDetection(vertices, core)
    end

    function BoundingBoxDetection(box::Rect2{T}, core_args...) where T
        core = ObjectDetectionCore(core_args...)
        new{T}(box, core)
    end
end

function create_with_center(center::Point2{T}, width::T, height::T, core_args...) where T
    top_left = center - Point2{T}(width, height) / 2
    return BoundingBoxDetection(top_left, width, height, core_args...)
end

function centroid(detection::BoundingBoxDetection{T}) :: Point2{T} where T
    return bounding_box(detection).origin + bounding_box(detection).widths / 2
end

function bounding_box_detection(detection::BoundingBoxDetection{T}) :: BoundingBoxDetection{T} where T
    return detection
end


text:
# In this exercise we will try to predict the quality class of wines given some chemical characteristics

# In detail, the attributes of this dataset are:
#   1) Alcohol
#   2) Malic acid
#   3) Ash
#   4) Alcalinity of ash
#   5) Magnesium
#   6) Total phenols
#   7) Flavanoids
#   8) Nonflavanoid phenols
#   9) Proanthocyanins
#   10) Color intensity
#   11) Hue
#   12) OD280/OD315 of diluted wines
#   13) Proline

# Further information concerning this dataset can be found online on the [UCI Machine Learning Repository dedicated page](https://archive.ics.uci.edu/ml/datasets/wine) or in particular on [this file](https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.names)

# Our prediction concerns the quality class of the wine (1, 2 or 3) that is given as first column of the data.

# 1) Start by setting the working directory to the directory of this file and activate it. If you have the provided `Manifest.toml` file in the directory, just run `Pkg.instantiate()`, otherwise manually add the packages Pipe, HTTP, Plots and BetaML.
# Also, seed the random seed with the integer `123`.

# 2) Load the packages/modules DelimitedFiles, Pipe, HTTP, Plots, BetaML

# 3) Load from internet or from local file the input data as a Matrix.
# You can use `readdlm`` using the comma as field separator.
dataURL="https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"

# 4) Now create the X matrix of features using the second to final columns of the data you loaded above and the Y vector by taking the 1st column. Transform the Y vector to a vector of integers using the `Int()` function (broadcasted). Make shure you have a 178×13 matrix and a 178 elements vector

# 5) Partition the data in (xtrain,xval) and (ytrain,yval) keeping 80% of the data for training and reserving 35% for testing. Keep the default option to shuffle the data, as the input data isn't.

# 6) As the output is multinomial we need to encode ytrain. We use the `oneHotEncoder()` function to make `ytrain_oh`

# 7) Define a Neural Network model with the following characteristics:
#   - 3 dense layers with respectively 13, 20 and 3 nodes and activation function relu
#   - a `VectorFunctionLayer` with 3 nodes and `softmax` as activation function
#   - `crossEntropy` as the neural network cost function

# 8) Train your model using `ytrain` and a scaled version of `xtrain`` (where all columns have zero mean and 1 standard deviaiton) for 100 epochs and use a batch size of 6 records.
# Save the output of your training function to `trainingLogs`

# 9) Predict the training labels ŷtrain and the validation labels ŷval. Recall you did the training on the scaled features!

# 10) Compute the train and test accuracies using the function `accuracy`

# 11) Compute and print a ConfutionMatrix of the validation data true vs. predicted

# 12) Run the following commands to plots the average loss per epoch
plot(trainingLogs.ϵ_epochs)

# 13) (Optional) Run the same workflow without scaling the data or using `squaredCost` as cost function. How this affect the quality of your predictions ?

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
function gather_stdlib_uuids()
    return Set{UUIDs.UUID}(x for x in keys(Pkg.Types.stdlibs()))
end


text:
#=
[SSDPClient.jl]
version =  0.1
author = "William Herrera"
copyright = "Copyright 2018 William Herrera"
created = "27 January 2018"
purpose = "Detect devices on the local network using the SSDP protocol"
=#

module SSDPClient

export ssdpquery

# Julia > v0.7 needs this
if !isdefined(Base, :UDPSocket)
    using Sockets
end

"""
    ssdpquery(matchtext="")

Query local network for devices responding to query with the
Simple Service Discovery Protocol (SSDP).

The optional argument is a text to match, such as "Roku" or "Arduino".

If you specify a capture in the match text that will be returned,
otherwise the entry that matched will be returned in entirety.

NOTE: If you do not specify an argument all replies will be printed and
the function will loop until interrupted. This is useful when searching a network.
"""
function ssdpquery(matchtxt="")
    MULTICAST = ip"239.255.255.250"
    SERVERPORT = 1900
    SEARCHHEADER = """M-SEARCH * HTTP/1.1\r\nHOST: 239.255.255.250:1900\r\nMAN: "ssdp:discover"\r\nST: ssdp:all\r\nMX: 10\r\n\r\n"""
    try
        while true
            udpsock = UDPSocket()
            send(udpsock, MULTICAST,SERVERPORT, SEARCHHEADER)
            reply = join([Char(ch) for ch in recv(udpsock)], "")
            if matchtxt == ""
                println("Received reply:\r\n$reply")
            else
                regex = Regex(matchtxt, "s")
                if (mat = match(regex, reply)) != nothing
                    if length(mat.captures) > 0
                        return mat.captures
                    end
                    return reply
                end
            end
            close(udpsock)
            sleep(3)
        end
    catch y
        println("Caught error: $y")
        return ""
    end
end

end # module


text:
export Model2Fit_Transformed_Parameters
export get_model_θ, get_model
export create_model_transform_src_insert_dest

@doc raw"""

Create a new model, where the  parameter vector ``\theta`` is computed using an
[`Abstract_Map`](@ref):
```math
 \hat{\theta}\mapsto \theta = f_{\hat{\theta}_f}(\hat{\theta}_m)
```

The resulting model is thus:
```math
\hat{m}(X,\hat{\theta}=[\hat{\theta}_m,\hat{\theta}_f]) = m(X,\theta = f_{\hat{\theta}_f}(\hat{\theta}_m))
```

# Constructor

## General case

```julia
Model2Fit_Transformed_Parameters(model::Abstract_Model2Fit,
                                 map_domain_size::Int,
                                 map::Abstract_Map)
```

- `map_domain_size` is the expected ``\hat{\theta}_m`` size. This
  quantity is used to implement model `parameter_size()` method. The
  reason is that the current [`Abstract_Map`](@ref) interface only
  provides ``\hat{\theta}_f`` size (its `parameter_size()` method) but
  tells nothing about ``\hat{\theta}_m`` size.
- Also, please note that ``\theta`` size is not necessary equal to
  ``\hat{\theta}_m`` size

## Using `Transform_Src_Insert_Dest_Map`

Very often this structure is used with a map of type
[`Transform_Src_Insert_Dest_Map`](@ref). In that case it is
simpler to use a dedicated function:
[`create_model_transform_src_insert_dest`](@ref) (follow this link to
see an example).

# Example

TODO: add an example using the general constructor

# Extra methods

Beside the regular [`Abstract_Model2Fit`](@ref) interface, you can use these extra methods:

- [`get_model(model::Model2Fit_Transformed_Parameters)`](@ref)
- [`get_model_θ(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)`](@ref)

## Not exported

These methods are not exported but may be useful:
- [`get_model_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)`](@ref)
- [`get_map_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)`](@ref)
"""
struct Model2Fit_Transformed_Parameters{MODEL <: Abstract_Model2Fit,
                                        MAP <: Abstract_Map} <: Abstract_Model2Fit
    _model::MODEL
    _map_domain_size::Int # length(X) where map(X,θ̂f), mandatory to define parameter_size
    # _map_codomain_size must be equal to the wrapped _model parameter_size
    _map::MAP

    function Model2Fit_Transformed_Parameters(model::MODEL,
                                              map_domain_size::Int,
                                              map::MAP) where {MODEL <: Abstract_Model2Fit,
                                                               MAP <: Abstract_Map}
        new{MODEL,MAP}(model, map_domain_size, map)
    end

end

@doc raw"""

This function creates an [`Model2Fit_Transformed_Parameters`](@ref)
instance from a map of type
[`Transform_Src_Insert_Dest_Map`](@ref).

```julia
function create_model_transform_src_insert_dest(model, map, src=>dest) -> model
```

# Example

```jldoctest
using NLS_Fit

# Two Gaussian peaks with σ(μ) an affine function
#
model = Gaussian_Peak() + Gaussian_Peak()

# model θ is:
# idx:  1   2  3   4   5   6
#  θ : h1, μ1, σ1, h2, μ2, σ2
#
# hence:
#
src  = [2,5]
dest = [3,6]

# The σ(μ) relation: σ(μ=1)=1 and σ(μ=100)=10 for [θ̂A, θ̂B]=[1, 1]
#
f_σ_μ = Map_Affine(1.0=>1.0,100.0=>10.0)

# Now the model_σ_μ parameters are:
#
# [ h1, μ1, h2, μ2, θ̂A, θ̂B ]
#
# the initial model is called with:
#
# [ h1, μ1, σ1 = affine(μ1, θ̂A, θ̂B), h2, μ2, , σ2 = affine(μ2, θ̂A, θ̂B) ]
#
model_σ_μ = create_model_transform_src_insert_dest(model,f_σ_μ,src=>dest)

# Examine that the compute θ takes into account the σ(μ) dependance:
#
θ̂ = Float64[1,5, 2, 95, 1, 1]
θ = get_model_θ(model_σ_μ, θ̂)

# output
6-element Vector{Float64}:
  1.0
  5.0
  1.3636363636363638
  2.0
 95.0
  9.545454545454547

```

"""
function create_model_transform_src_insert_dest(model::Abstract_Model2Fit,
                                                f_map::Abstract_Map,
                                                src_dest::Pair{SOURCE,DEST}) where {SOURCE<:AbstractVector{Int},DEST<:AbstractVector{Int}}
    @assert length(first(src_dest)) == length(last(src_dest))

    # Initialize g_map (see Transform_Src_Insert_Dest_Map
    #
    g_map = NLS_Fit.Transform_Src_Insert_Dest_Map(f_map, src_dest)

    # compute map_domain_size:
    # -> this is model θ length minus the number of inserted element
    #
    map_domain_size = parameter_size(model) - length(first(src_dest))
    @assert map_domain_size ≥ 0

    Model2Fit_Transformed_Parameters(model, map_domain_size, g_map)
end

# Specific methods  ================
#
@doc raw"""
```julia
get_model_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)
```

Return a view on ``\hat{\theta}_m`` knowing ``\hat{\theta}=[\hat{\theta}_m,\hat{\theta}_f]``
"""
function get_model_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)
    @assert parameter_size(hat_model) == length(hat_θ)

    @view hat_θ[1:hat_model._map_domain_size]
end

@doc raw"""
```julia
get_model_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)
```

Return a view on ``\hat{\theta}_f`` knowing ``\hat{\theta}=[\hat{\theta}_m,\hat{\theta}_f]``
"""
function get_map_hat_θ_view(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)
    @assert parameter_size(hat_model) == length(hat_θ)

    @view hat_θ[(hat_model._map_domain_size+1):end]
end

@doc raw"""
```julia
get_model(mp::Model2Fit_Transformed_Parameters)::Abstract_Model2Fit
```

Get back the wrapped model

See [`Model2Fit_Transformed_Parameters`](@ref)
"""
get_model(model::Model2Fit_Transformed_Parameters) = model._model

@doc raw"""
```julia
get_model_θ(mp::Model2Fit_Transformed_Parameters, hat_θ::AbstractVector) -> θ::AbstractVector
```

Retrieve the parameter vector θ associated to the wrapped model [`get_model`](@ref).
"""
function get_model_θ(hat_model::Model2Fit_Transformed_Parameters,hat_θ::AbstractVector)
    model_hat_θ = get_model_hat_θ_view(hat_model,hat_θ)
    map_hat_θ = get_map_hat_θ_view(hat_model,hat_θ)

    θ = eval_map(hat_model._map,model_hat_θ,map_hat_θ)

    # check that map codomain = model parameter size
    @assert parameter_size(get_model(hat_model)) == length(θ)

    θ
end

# Visit  ================
#
visit_submodel_size(hat_model::Model2Fit_Transformed_Parameters) = 1
visit_get_submodel(hat_model::Model2Fit_Transformed_Parameters,submodel_idx::Int) = get_model(hat_model)
visit_get_X(hat_model::Model2Fit_Transformed_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,hat_θ::AbstractVector) = X
visit_get_Y(hat_model::Model2Fit_Transformed_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,hat_θ::AbstractVector) = Y
visit_get_θ(hat_model::Model2Fit_Transformed_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,hat_θ::AbstractVector) = get_model_θ(hat_model,hat_θ)

# Interface ================
#
parameter_size(hat_model::Model2Fit_Transformed_Parameters) = hat_model._map_domain_size + parameter_size(hat_model._map)

function accumulate_y!(hat_model::Model2Fit_Transformed_Parameters,Y::AbstractVector,X::AbstractVector,hat_θ::AbstractVector)
    model = get_model(hat_model)
    θ = get_model_θ(hat_model,hat_θ)

    accumulate_y!(model,Y,X,θ)
end


text: include("solvers/cem.jl")


text:
import Pkg
Pkg.activate(".")
Pkg.add("Plots")
Pkg.add("AbstractTrees")
Pkg.add("Cascadia")
Pkg.add("Clustering")
Pkg.add("CodecZlib")
Pkg.add("ColorTypes")
Pkg.add("Convex")
Pkg.add("CoordinateTransformations")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("DataFramesMeta")
Pkg.add("Dates")
Pkg.add("DecisionTree")
Pkg.add("Distributed")
Pkg.add("Distributions")
Pkg.add("ECOS")
Pkg.add("ExcelFiles")
Pkg.add("ExcelReaders")
Pkg.add("FFTW")
Pkg.add("FredData")
Pkg.add("Gadfly")
Pkg.add("GeometryTypes")
Pkg.add("GLM")
Pkg.add("Glob")
Pkg.add("Gumbo")
Pkg.add("HDF5")
Pkg.add("HTTP")
Pkg.add("JLD2")
Pkg.add("JSON")
Pkg.add("JSONTables")
Pkg.add("JuliaDB")
Pkg.add("Lazy")
Pkg.add("LIBSVM")
Pkg.add("LightXML")
Pkg.add("LinearAlgebra")
Pkg.add("MacroTools")
Pkg.add("MeshCat")
Pkg.add("Plots")
Pkg.add("Random")
Pkg.add("RDatasets")
Pkg.add("ScikitLearn")
Pkg.add("SCS")
Pkg.add("SharedArrays")
Pkg.add("Statistics")
Pkg.add("StatsBase")
Pkg.add("TimeZones")
Pkg.add("PackageCompiler")

max_stars_repo_name: AkillesAILimited/Thetis

max_issues_repo_name: AkillesAILimited/Thetis

max_forks_repo_name: AkillesAILimited/Thetis


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

module GeoStatsBase

using CSV
using Optim
using Random: randperm, shuffle
using Combinatorics: multiexponents
using Distributed: CachingPool, pmap, myid
using LinearAlgebra: Diagonal, normalize, norm, ⋅
using StatsBase: Histogram, Weights, AbstractWeights
using Distances: Metric, Euclidean, Mahalanobis, pairwise
using Distributions: ContinuousUnivariateDistribution, median, mode
using CategoricalArrays: CategoricalValue, CategoricalArray
using CategoricalArrays: levels, isordered, pool, levelcode
using NearestNeighbors: KDTree, BallTree, knn, inrange
using ReferenceFrameRotations: angle_to_dcm
using DataFrames: DataFrame, DataFrame!
using StaticArrays: SVector, MVector, SOneTo
using AverageShiftedHistograms: ash
using Transducers: Map, foldxt
using SpecialFunctions: gamma
using DensityRatioEstimation
using ScientificTypes
using LossFunctions
using RecipesBase
using Parameters

import Tables
import MLJModelInterface
import Base: values, ==
import Base: in, filter, map, split, error
import StatsBase: fit, sample, varcorrection
import Statistics: mean, var, quantile
import Distributions: quantile, cdf
import ScientificTypes: Scitype, scitype
import Distances: evaluate
import DataFrames: groupby
import NearestNeighbors: MinkowskiMetric

# aliases
const MI = MLJModelInterface
const Vec{N,T} = Union{SVector{N,T},MVector{N,T}}

# convention of scientific types
include("convention.jl")

function __init__()
  ScientificTypes.set_convention(GeoStats())
end

include("variables.jl")
include("geotraits.jl")
include("domains.jl")
include("data.jl")
include("ensembles.jl")
include("georef.jl")
include("macros.jl")
include("paths.jl")
include("trends.jl")
include("geometries.jl")
include("distances.jl")
include("neighborhoods.jl")
include("neighborsearch.jl")
include("distributions.jl")
include("estimators.jl")
include("partitioning.jl")
include("weighting.jl")
include("discretizing.jl")
include("sampling.jl")
include("geoops.jl")
include("learning.jl")
include("mappings.jl")
include("problems.jl")
include("solvers.jl")
include("folding.jl")
include("errors.jl")
include("statistics.jl")
include("plotrecipes.jl")
include("utils.jl")

export
  # spatial variable
  Variable,
  name, mactype,

  # geotraits
  nelms,
  ncoords,
  coordtype,
  coordinates,
  coordinates!,
  domain,
  values,

  # spatial domains
  AbstractDomain,
  PointSet,
  RegularGrid,
  StructuredGrid,
  origin, spacing,

  # spatial data
  AbstractData,
  SpatialData,
  variables,
  georef,

  # spatial ensembles
  Ensemble,

  # mapping
  MappingMethod,
  NearestMapping,
  CopyMapping,

  # learning tasks
  AbstractLearningTask,
  SupervisedLearningTask,
  UnsupervisedLearningTask,
  RegressionTask,
  ClassificationTask,
  ClusteringTask,
  issupervised,
  inputvars,
  outputvars,
  features,
  label,

  # learning models
  issupervised,
  isprobabilistic,
  learn, perform,

  # learning losses
  defaultloss,

  # problems
  AbstractProblem,
  EstimationProblem,
  SimulationProblem,
  LearningProblem,
  data,
  domain,
  sourcedata,
  targetdata,
  task,
  variables,
  coordinates,
  hasdata,
  nreals,

  # solvers
  AbstractSolver,
  AbstractEstimationSolver,
  AbstractSimulationSolver,
  AbstractLearningSolver,
  SeqSim,
  SeqSimParam,
  CookieCutter,
  CookieCutterParam,
  PointwiseLearn,
  variables,
  covariables,
  preprocess,
  solve, solvesingle,

  # folding
  FoldingMethod,
  RandomFolding,
  PointFolding,
  BlockFolding,
  BallFolding,
  folds,

  # errors
  ErrorEstimationMethod,
  LeaveOneOut,
  LeaveBallOut,
  CrossValidation,
  BlockCrossValidation,
  WeightedCrossValidation,
  DensityRatioValidation,

  # helper macros
  @estimsolver,
  @simsolver,

  # paths
  AbstractPath,
  LinearPath,
  RandomPath,
  SourcePath,
  ShiftedPath,
  traverse,

  # geometries
  AbstractGeometry,
  Rectangle,
  sides,
  center,
  diagonal,
  volume,

  # distances
  aniso2distance,
  evaluate,

  # neighborhoods
  AbstractNeighborhood,
  BallNeighborhood,
  isneighbor,
  radius,
  metric,

  # neighborhood search
  NeighborSearchMethod,
  BoundedNeighborSearchMethod,
  KNearestSearch,
  NeighborhoodSearch,
  BoundedSearch,
  KBallSearch,
  search!, search,
  maxneighbors,
  object,

  # distributions
  EmpiricalDistribution,
  transform!, quantile, cdf,

  # estimators
  fit, predict, status,

  # partitioning
  SpatialPartition,
  PartitionMethod,
  PredicatePartitionMethod,
  SPredicatePartitionMethod,
  RandomPartition,
  FractionPartition,
  SLICPartition,
  BlockPartition,
  BisectPointPartition,
  BisectFractionPartition,
  BallPartition,
  PlanePartition,
  DirectionPartition,
  PredicatePartition,
  SPredicatePartition,
  VariablePartition,
  ProductPartition,
  HierarchicalPartition,
  partition,
  subsets,
  metadata,
  →,

  # weighting
  SpatialWeights,
  WeightingMethod,
  UniformWeighting,
  BlockWeighting,
  DensityRatioWeighting,
  weight,

  # discretizing
  DiscretizationMethod,
  BlockDiscretization,
  discretize,

  # sampling
  SamplingMethod,
  UniformSampling,
  WeightedSampling,
  BallSampling,
  sample,

  # operations
  uniquecoords,
  inside,

  # trends
  polymat,
  trend,

  # statistics
  EmpiricalHistogram,
  mean, var, quantile,

  # plot recipes
  cornerplot,

  # utilities
  readgeotable,
  groupby,
  boundbox,
  slice,
  spheredir

end


text:
"""
    get_setup(SetupType)

Get a preconfigured setup instance of type `SetupType`.
"""
function get_setup end

get_setup(S::Type{PlateSetup{T}}) where {T} = S(;
    name = "plates/someplates",
    width = 20.0,
    depth = 20.0,
    heights = fill(5.0, 4),
    bend = 0.0,
    twist = 0.0,
    # refinement = 0.5,
)

get_setup(S::Type{CylinderSetup{T}}) where {T} = S(;
    name = "cylinders/somecylinders",
    ncell = 1,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 1.0,
    bend = 0.0,
    twist = 0.0,
    include_in = true,
    in_ratio = 0.6,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
    # refinement = 0.5,
)

get_setup(S::Type{SphereSetup{T}}) where {T} = S(;
    name = "spheres/somespheres",
    ncell = 4,
    rmin = 3.0,
    rmax = 5.0,
    dmin = 0.2,
    dmax = 0.3,
    include_in = false,
    in_ratio = 0.7,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.3,
    # refinement = 0.5,
)

get_setup(S::Type{NeuronSetup{T}}) where {T} = S(;
    name = "neuron",
    ecs_shape = :no_ecs,
    ecs_ratio = 0.3,
    # refinement = 0.5,
)

"""
    get_coeffs(setup)

Get preconfigured compartment coefficients for `setup`.
"""
function get_coeffs end

get_coeffs(setup::PlateSetup) = coefficients(
    setup;
    D = [0.002 * I(3) for _ = 1:length(setup.heights)],
    T₂ = [Inf for _ = 1:length(setup.heights)],
    ρ = [1.0 for _ = 1:length(setup.heights)],
    κ = (;
        interfaces = [1e-4 for _ = 1:length(setup.heights)-1],
        boundaries = [0.0 for _ = 1:length(setup.heights)],
    ),
    γ = 2.67513e-4,
)

get_coeffs(setup::CylinderSetup) = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

get_coeffs(setup::SphereSetup) = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

get_coeffs(setup::NeuronSetup) = coefficients(
    setup;
    D = (; neuron = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; neuron = Inf, ecs = Inf),
    ρ = (; n = 1.0, neuron = 1.0, ecs = 1.0),
    κ = (; neuron_ecs = 1e-4, neuron = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)


text:
using Test
using NCDatasets
using CalibrateEDMF.LESUtils

@testset "LESUtils" begin
    scm_names = ["total_flux_qt", "total_flux_h", "u_mean", "foo"]
    tmpdir = mktempdir()
    filename = joinpath(tmpdir, "stats", "test.nc")
    mkdir(dirname(filename))
    NCDataset(filename, "c") do ds

        defDim(ds, "t", 10)
        defDim(ds, "z", 20)

        defGroup(ds, "profiles")
        defVar(ds.group["profiles"], "qt_flux_z", Float32, ("t", "z"))
        defVar(ds.group["profiles"], "resolved_z_flux_thetali", Float32, ("t", "z"))
        defVar(ds.group["profiles"], "u_translational_mean", Float32, ("t", "z"))
        defGroup(ds, "timeseries")
        defVar(ds.group["timeseries"], "foo", Float32, ("t",))
    end

    @test find_alias(("u_mean", "u_translational_mean"), filename) == "u_translational_mean"
    @test find_alias(("foo",), filename) == "foo"
    @test_throws ErrorException find_alias(("vorticity",), filename)
    @test get_les_names(scm_names, filename) == ["qt_flux_z", "resolved_z_flux_thetali", "u_translational_mean", "foo"]

    @testset "cfSite_getter" begin
        @test isa(get_cfsite_les_dir(2), String)
        @test_throws AssertionError get_cfsite_les_dir(75)
        @test_throws AssertionError get_cfsite_les_dir(2, forcing_model = "foo")
    end

    @testset "artifact_getter" begin
        art_dir = get_path_to_artifact()
        @test isdir(art_dir)
        @test_throws KeyError get_path_to_artifact("foo")
    end
end


text:
# Basic polynomials

function poly1(x)
    # 0 (2 times)
    return x^2
end

function poly2(x)
    # 0 (5 times)
    return x^5
end

function poly3(x)
    # -1.222 + i (2 times)
    return (x+1.222-im)^2
end

function poly4(x)
    # 2, -1
    return (x - 2) * (x + 1)
end

function poly5(x)
    # i, -i
    return x^2 + 1
end

function poly6(x)
    # 2i, -3i
    return (x - 2im) * (x + 3im)
end

function poly7(x)
    # i, -i, -10i
    return (x - im) * (x + im) * (x + 10im)
end

function poly8(x)
    # i, -i, -1.5i
    return (x - im) * (x + im) * (x + 1.5im)
end

function poly9(x)
    # 2, -2, 3
    return (x - 2) * (x + 2) * (x - 3)
end

function poly10(x)
    # 2i, -2i, 3
    return (x - 2 * im) * (x + 2 * im) * (x - 3)
end

function poly11(x)
    # 3 roots
    return x^3 - 1
end

# Irrational roots
function poly12(x)
    return x + 3 * (pi / 8 + pi / 8 * im)
end

function poly13(x)
    return x - (pi / 4 + pi / 8 * im)
end

function poly14(x)
    return x - (exp(1) + exp(1) * im)
end

function poly15(x)
    return x - exp(1)
end

function poly16(x)
    return x + pi * im
end

function poly17(x)
    return (x - exp(1))^3
end

function poly18(x)
    return (x - exp(1) * im)^3
end

# Others

function poly19(x)
    return (x - exp(1) * im)^3 * (x - exp(1))^2
end

function poly20(x)
    return (x - 2pi * im) * (x - pi)^3 * (x + 3pi * im -2.1)^2
end

function poly21(x)
    return (x - 0.001 + 0.5*im)^2
end

function poly22(x)
    return (x - 3.4656) * (x + 3.99)^2
end

function poly100(x)
    return x^8
end

function poly101(x)
    return (x^8) * (x - 101) * 9
end

function poly102(x)
    return (x - pi)^21(x + pi * im)^18
end

function poly103(x)
    return x^25
end

function poly104(x)
    return x^17
end

function poly105(x)
    return (x^9) * (x - pi * im)^5
end

function poly106(x)
    return x^10
end

function poly107(x)
    return (x - 2.185935 * im)^10
end

function poly109(x)
    return x^59
end

function poly110(x)
    return x^60
end


text:
"""
    _merge!(x::AbstractArray, p::Int, q::Int, r::Int)

Merges tow sorted arrays inplace of size `q-p+1` and `r-q` into one
sorted array. For private use only.
"""
function _merge!(x::AbstractArray, p::Int, q::Int, r::Int)
    @inbounds x1 = x[p:q]
    @inbounds x2 = x[(q+1):r]

    i = p
    j = q+1
    k = p

    while i ≤ q && j ≤ r
        @inbounds if x1[i-p+1] < x2[j-q]
            @inbounds x[k] = x1[i-p+1]
            k += 1
            i += 1
        else
            @inbounds x[k] = x2[j-q]
            k += 1
            j += 1
        end
    end

    while i ≤ q
        @inbounds x[k] = x1[i-p+1]
        k += 1
        i += 1
    end

    while j ≤ r
        @inbounds x[k] = x2[j-q]
        k += 1
        j += 1
    end

    return nothing
end

function mergesort!(x::AbstractArray, p::Int, r::Int)
    if p ≥ r
        return nothing
    end
    q::Int = floor(( p+r ) / 2)
    mergesort!(x, p, q)
    mergesort!(x, q+1, r)
    _merge!(x, p, q, r)
    return nothing
end;


text:
type DoubleDescription{T<:Real}
    A::Matrix{T}
    R::Vector
    K::Set{Int}
    adj::Dict{(Int,Int),Bool}
    num_rays::Int
end

type CountedVector{T<:Real}
    v::Vector{T}
    Av::Vector{T}
    dd::DoubleDescription{T}
    id::Int

    function CountedVector(v::Vector{T},dd::DoubleDescription{T})
        dd.num_rays += 1
        canonicalize!(v)
        new(v,dd.A*v,dd,dd.num_rays)
    end
end
CountedVector{T<:Real}(v::Vector{T},dd::DoubleDescription{T}) = CountedVector{T}(v,dd)

Base.vec(c::CountedVector) = c.v

function initial_description{T<:Real}(A::Matrix{T})
    m,n = size(A)
    B = rref(A')
    # find pivots
    r = 1
    K = Set{Int}()
    for i in 1:m
        comp = zeros(n)
        comp[r] = 1
        if B[:,i] == comp
            r += 1
            push!(K, i)
        end
        r > n && break
    end
    cK = sort(collect(K))
    Ak = A[cK,:]
    R = Ak \ eye(n,n)
    dd = DoubleDescription(A,CountedVector{T}[],K,Dict{(Int,Int),Bool}(),n)
    Rk = [CountedVector{T}(R[:,i],dd) for i in 1:n]
    dd.R = kR
    for i in 1:n
        # Ar = Rk[i].Av[cK]
        for j in (i+1):n
            # As = Rk[j].Av[cK]
            id = extrema([Rk[i].id,Rk[j].id])
            cache_adjacency!(dd, n, Rk[i].Av[cK], Rk[j].Av[cK], id)
        end
    end
    return dd
end

function double_description{T<:Real}(A::Matrix{T})
    A = [zeros(T,1,size(A,2)); A]
    A[1,1] = one(T)
    m, n = size(A)
    dd = initial_description(A)
    Kᶜ = setdiff(1:m, dd.K)
    while !isempty(Kᶜ)
        i = pop!(Kᶜ)
        println("iteration $(length(dd.K)+1) of $m")
        # length(dd.K) == 31 && break
        update!(dd, i)
    end
    return dd
end

function is_approx_included(haystack, needle)
    n = length(needle)
    diff = zeros(n)
    for h in haystack
        diff = norm(vec(h)-needle)
        diff < n*ε && return true
    end
    return false
end

function canonicalize!{T<:Real}(v::Vector{T})
    n = length(v)
    val = abs(v[1])
    if val < ε
        val = abs(v[findfirst(abs(v) .> ε)])
    end
    for i in 1:n
        v[i] = v[i] / val
    end
    v
end

# use Lemma 8 from Fukuda (1996) to update the double description
function update!{T<:Real}(dd::DoubleDescription{T}, i)
    m, n = size(dd.A)
    Aᵢ = reshape(dd.A[i,:], (n,))
    Rⁿᵉʷ = CountedVector{T}[]
    R⁺, R⁰, R⁻ = partition_rays(dd.R, Aᵢ)
    for r in R⁺, s in R⁻
        if isadjacent(dd,r,s)
            w = dot(Aᵢ,vec(r))*vec(s) - dot(Aᵢ,vec(s))*vec(r)
            v = CountedVector(w,dd)
            if sum(abs(w)) > n*ε &&
               !is_approx_included(R⁰,   vec(w)) &&
               !is_approx_included(Rⁿᵉʷ, vec(w))
                dd.num_rays += 1
                push!(Rⁿᵉʷ, v)
            end
        end
    end
    dd.R = vcat(R⁺, R⁰, Rⁿᵉʷ)
    push!(dd.K, i)
    cK = sort(collect(dd.K))
    Ak = dd.A[cK,:]
    # should really add a test right about here to ensure
    # that old rays do not become adjacent...I think this
    # can only happen if both v,w ∈ R⁰
    d = rank(Ak)
    for s in Rⁿᵉʷ
        # As = s.Av[cK]#Ak*vec(s)
        for r in dd.R
            r.id == s.id && continue
            # Ar = r.Av[cK]#Ak*vec(r)
            id = extrema([r.id, s.id])
            cache_adjacency!(dd, d, r.Av[cK], s.Av[cK], id)
        end
    end
    nothing
end

function partition_rays{T<:Real}(R::Vector{CountedVector{T}}, a::Vector{T})
    R⁺, R⁰, R⁻ = CountedVector{T}[], CountedVector{T}[], CountedVector{T}[]
    n = length(a)
    for r in R
        if sum(abs(vec(r))) < n*ε
            println("have a zero vector!")
            continue
        end
        val = dot(a,vec(r))
        if val > ε
            push!(R⁺, r)
        elseif val < -ε
            push!(R⁻, r)
        else
            push!(R⁰, r)
        end
    end
    return R⁺, R⁰, R⁻
end

isadjacent(dd, r, s) = dd.adj[extrema([r.id,s.id])]

function cache_adjacency!(dd, d, Ar, As, id::(Int,Int))
    Z = active_sets(dd, Ar, As)
    if length(Z) < d - 2
        val = false
    elseif length(intersect(Z,dd.K)) ≥ d - 2
        val = true
    else
        val = (rank(dd.A[sort(collect(dd.K)),:][Z,:]) == d - 2)
    end
    dd.adj[id] = val
    return val
end

function active_sets(dd, Ar, As)
    Z = Int[]
    m = length(Ar)
    sizehint!(Z,m)
    for i in 1:m
        if abs(Ar[i]) ≤ ε && abs(As[i]) ≤ ε
            push!(Z, i)
        end
    end
    return Z
end


text:
using DelimitedFiles

"""
Insert transcribtion generated from speech to text software to a TextGrid file.

Input is the text file generated by speech to text tool which contains the words and corresponding start and end time of word occurance in seconds.
Output is the tartget TextGrid file to include the transcription on as an "IntervalTier" interval.
"""
function insert(input,output)
    words = readdlm(input)

    # set item info
    class = "IntervalTier"
    name = "Words"
    size1 = size(words)[1]
    xmin = words[1,:][2]
    xmax = words[end,:][3]

    open(output, "a") do file
        write(file, " "^4*"item [10]:\n")
        write(file, " "^8*"class = \""*class*"\" \n")
        write(file, " "^8*"name = \""*name*"\" \n")
        write(file, " "^8*"xmin = "*string(xmin)*" \n")
        write(file, " "^8*"xmax = "*string(xmax)*" \n")
        write(file, " "^8*"intervals: size = "*string(size1)*" \n")
        for i in 1:size1
            write(file, " "^8*"intervals ["*string(i)*"]:\n")
            write(file, " "^12*"xmin = "*string(words[i,:][2])*" \n")
            write(file, " "^12*"xmax = "*string(words[i,:][3])*" \n")
            write(file, " "^12*"text = \""*string(words[i,:][1])*"\" \n")
        end
    end
end

input = raw"C:\Users\hemad\Desktop\Master\GoogleSpeech\Example1\file2\transcription.txt"
output = raw"C:\Users\hemad\Desktop\Master\GoogleSpeech\Example1\file2\A032_Sarah_edited.TextGrid"
insert(input,output)


text:
# export delaunayWall;
# #===============================================================================
# #
# #	src/deWall.jl
# #
# #	This File Contains:
# #
# #	 - delaunayWall(
# #			P::Matrix,
# #			ax = 1,
# #			Pblack = Float64[],
# #			AFL = Array{Int64,1}[],
# #			tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
# #			DEBUG = false
# #		)::Array{Array{Int64,1},1}
# #
# #	 - firstDeWallSimplex(
# #			P::Matrix,
# #			ax::Int64,
# #			off::Float64;
# #			DEBUG = false
# #		)::Array{Int64,1}
# #
# #	 - findWallSimplex(
# #			P::Matrix,
# #			blackidx::Int64,
# #			face::Array{Int64,1},
# #			oppoint::Array{Float64,1};
# #			DEBUG = false
# #		)::Union{Array{Int64,1}, Nothing}
# #
# #	 - recursiveDelaunayWall(
# # 			P::Matrix,
# #			Pblack::Matrix,
# # 			tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# # 			AFL::Array{Array{Int64,1},1},
# # 			ax::Int64,
# # 			off::Float64,
# # 			positive::Bool;
# # 			DEBUG = false
# # 		)::Array{Array{Int64,1},1}
# #
# #	 - updateAFL!(
# #			P::Matrix
# #			new::Array{Int64,1}[],
# #			AFLα = Array{Int64,1}[],
# #			AFLplus = Array{Int64,1}[],
# #			AFLminus = Array{Int64,1}[],
# #			ax::Int64,
# #			off::Float64;
# # 			DEBUG = false
# #		)::Bool
# #
# #	 - updatelist!(
# #			list,
# #			element
# #		)::Bool
# #
# #	 - updateTetraDict!(
# #			P::Matrix,
# #			tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# #			AFL::Array{Array{Int64,1},1},
# #			σ::Array{Int64,1}
# #		)::Nothing
# #
# ===============================================================================#
#
# """
# 	delaunayWall(
# 		P::Matrix, ax = 1, Pblack::Float64[], AFL = Array{Int64,1}[],
# 		tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
# 		DEBUG = false
# 	)::Array{Array{Int64,1},1}
#
# Return the Delaunay Triangulation of sites `P` via Delaunay Wall algorithm.
# The optional argument `ax` specify on wich axis it will build the Wall.
# The optional argument `AFL` is used in recursive call.
# If the keyword argument `DEBUG` is set to true than all the procedure is shown.
#
# # Examples
# ```jldoctest
#
# julia> P = [
#  0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
#  0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
#  0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
# ];
#
# julia> DT = AlphaStructures.delaunayWall(P)
# 6-element Array{Array{Int64,1},1}:
#  [1, 2, 4, 5]
#  [1, 3, 4, 5]
#  [2, 4, 5, 6]
#  [3, 4, 5, 7]
#  [4, 5, 6, 7]
#  [4, 6, 7, 8]
#
# ```
# """
# function delaunayWall(
# 		P::Matrix,
# 		ax = 1,
# 		Pblack = Float64[],
# 		AFL = Array{Int64,1}[],
# 		tetraDict = DataStructures.Dict{Array{Int64,1},Array{Float64,1}}();
# 		DEBUG = false
# 	)::Array{Array{Int64,1},1}
#
# 	if DEBUG @show "Delaunay Wall with parameters" P ax AFL tetraDict end
#
# 	# 0 - Data Reading and Container definition
# 	DT = Array{Int64,1}[]		# Delaunay Triangulation
# 	AFLα = Array{Int64,1}[]		# (d-1)faces intersecting the Wall
# 	AFLplus = Array{Int64,1}[]  # (d-1)faces in positive Wall half-space
# 	AFLminus = Array{Int64,1}[] # (d-1)faces in positive Wall half-space
# 	off = AlphaStructures.findMedian(P, ax)
# 	if !isempty(Pblack) Pext = [P Pblack] else Pext = copy(P) end
#
# 	# 1 - Determine first simplex (if necessary)
# 	if isempty(AFL)
# 		@assert isempty(Pblack) "delaunayWall: If AFL is empty => Pblack must be"
# 		@assert isempty(tetraDict) "delaunayWall: If AFL is empty => tetraDict must be"
# 		σ = sort(AlphaStructures.firstDeWallSimplex(P, ax, off, DEBUG = DEBUG))
# 		push!(DT, σ)
# 		AFL = AlphaStructures.simplexFaces(σ)
# 		AlphaStructures.updateTetraDict!(P, tetraDict, AFL, σ)
# 	else
# 		@assert !isempty(Pblack) "delaunayWall: Data missing - Pblack"
# 		@assert !isempty(AFL) "delaunayWall: Data missing - AFL"
# 		@assert !isempty(tetraDict) "delaunayWall: Data missing - tetraDict"
# 	end
#
# 	# 2 - Build `AFL*` according to the axis `ax` with constant term `off`
# 	AlphaStructures.updateAFL!(
# 		P, AFL, AFLα, AFLplus, AFLminus, ax, off, DEBUG = DEBUG
# 	)
#
# 	# 4 - Build simplex Wall
# 	while !isempty(AFLα)
# 		# if face ∈ keys(tetraDict) oppoint = tetraDict[face]
# 		# else Pselection = setdiff([i for i = 1 : n], face) end
# 		σ = AlphaStructures.findWallSimplex(
# 				Pext, AFLα[1], tetraDict[AFLα[1]], size(P, 2), DEBUG = DEBUG
# 			)
# 		if σ != nothing && σ ∉ DT
# 			push!(DT, σ)
# 			AFL = AlphaStructures.simplexFaces(σ)
# 			AlphaStructures.updateTetraDict!(P, tetraDict, AFL, σ)
# 			# Split σ's Faces according in semi-spaces
# 			AlphaStructures.updateAFL!(
# 				P, AFL, AFLα, AFLplus, AFLminus, ax, off, DEBUG=DEBUG
# 			)
# 		else
# 			@assert AlphaStructures.updatelist!(AFLα, AFLα[1]) == false "delaunayWall:
# 				Something unespected happends while removing a face."
# 		end
# 	end
#
# 	# 5 - Change the axis `ax` and repeat until there are no faces but exposed.
# 	#      A.K.A. Divide & Conquer phase.
# 	if !isempty(AFLminus)
# 		union!(DT, recursiveDelaunayWall(
# 			P, Pblack, tetraDict, AFLminus, ax, off, false; DEBUG = DEBUG
# 		))
# 	end
# 	if !isempty(AFLplus)
# 		union!(DT, recursiveDelaunayWall(
# 			P, Pblack, tetraDict, AFLplus, ax, off, true; DEBUG = DEBUG
# 		))
# 	end
#
# 	return DT
# end
#
# #-------------------------------------------------------------------------------
# """
# 	findWallSimplex(
# 		P::Matrix,
# 		face::Array{Int64,1}, oppoint::Array{Float64,1},
# 		blackidx = size(P, 2);
# 		DEBUG = false
# 	)::Union{Array{Int64,1}, Nothing}
#
# Returns the simplex 'σ' build with `face` and a point from `P[:, 1:blackidx]`
# such that it is in the opposite half plane of `oppoint`.
# If such a simplex do not exists it returns `nothing`.
# If `blackidx` is not specified all the points `P` are treaten as valid.
# If the keyword argument `DEBUG` is set to true than all the procedure is shown.
#
# # Examples
# ```jldoctest
#
# julia> P = [
#  0. 1. 0. 0. 2.;
#  0. 0. 1. 0. 2.;
#  0. 0. 0. 1. 2.
# ];
#
# julia> newtetra = AlphaStructures.findWallSimplex(P,[2,3,4],[0., 0., 0.])
# 4-element Array{Int64,1}:
#  2
#  3
#  4
#  5
#
# ```
# """
# function findWallSimplex(
# 		P::Matrix,
# 		face::Array{Int64,1},
# 		oppoint::Array{Float64,1},
# 		blackidx = size(P, 2);
# 		DEBUG = false
# 	)::Union{Array{Int64,1}, Nothing}
#
# 	if DEBUG @show "find Wall Simplex of" face oppoint end
# 	# Find the points in the halfspace defined by `face` that do not
# 	#  containsother the other point of the simplex.
# 	Pselection =
# 		AlphaStructures.oppositeHalfSpacePoints(P, P[:, face], oppoint)
#
# 	if DEBUG @show Pselection end
#
# 	# If there are no such points than the face is part of the convex hull.
# 	if isempty(Pselection)
# 		return nothing
# 	end
#
# 	# Find the Closest Point in the other halfspace with respect to σ
# 	#  according to dd-distance.
# 	idxbase = Pselection[ AlphaStructures.findClosestPoint(
# 		P[:, face], P[:, Pselection], metric = "dd"
# 	) ]
#
# 	# @assert !isnothing(idxbase)
# 	# if isnothing(idxbase)
# 	# 	return nothing
# 	# end
#
# 	# It prevent from adding the same simplex again (cause it has been
# 	#  determined in a previous recursive call in the stacktrace).
# 	if idxbase > blackidx
# 		if DEBUG println("Excluding $face cause simplex already inside.") end
# 		return nothing
# 	end
#
# 	σ = sort([face; idxbase])
# 	if DEBUG @show "Found face" σ end
#
# 	# Check the simplex correctness
# 	radius, center = AlphaStructures.findRadius(P[:, σ], true)
# 	for i = 1 : size(P, 2)
# 		if LinearAlgebra.norm(center - P[:, i]) < radius - 1.e-14
# 			# @assert i ∉ Pselection "ERROR: Numerical error
# 			# 	evaluating minimum radius for $σ"
# 			if DEBUG println("$σ discarded due to a closer point.") end
# 			return nothing
# 		end
# 	end
#
# 	return σ
# end
#
#
# #-------------------------------------------------------------------------------
#
# """
#  	firstDeWallSimplex(
# 		P::Matrix, ax::Int64, off::Float64;
# 		DEBUG = false
# 	)::Array{Int64,1}
#
# Returns the indices array of the points in `P` that form the first thetrahedron
# built over the Wall if the `ax` axes with contant term `off`.
# If the keyword argument `DEBUG` is set to true than all the procedure is shown.
#
# # Examples
# ```jldoctest
#
# julia> V = [
# 	0.0 1.0 0.0 0.0 2.0
# 	0.0 0.0 1.0 0.0 2.0
# 	0.0 0.0 0.0 1.0 2.0
# ];
#
# julia> firstDeWallSimplex(V, 1, AlphaStructures.findMedian(V,1))
# 4-element Array{Int64,1}:
#  1
#  2
#  3
#  4
#
# ```
# """
# function firstDeWallSimplex(
# 		P::Matrix,
# 		ax::Int64,
# 		off::Float64;
# 		DEBUG = false
# 	)::Array{Int64,1}
#
# 	dim = size(P, 1)
#     n = size(P, 2)
#
# 	if DEBUG println("Determine first Simplex with ax = $ax") end
#     # the first point of the simplex is the one with coordinate `ax` maximal
#     #  such that it is less than `off` (closer to α from minus)
#
# 	Pselection = findall(x -> x < off, P[ax, :])
#
# 	# it gives an error if no point are less than `off`
# 	#  in fact it means that all the points are located on the median,
# 	#  with respect to `ax`.
# 	@assert !isempty(Pselection) "firstDeWallSimplex: not able to build the first Delaunay
# 		dimplex; all the points have the same `ax` coordinate."
#     newidx = Pselection[findmax(P[ax, Pselection])[2]]
#     # indices will store the indices of the simplex ...
#     indices = [newidx]                      #Array{Int64,1}
#     # ... and `Psimplex` will store the corresponding points
#     Psimplex = P[:, newidx][:,:]    #Array{Float64,2}
#
#     # the second point must be seeken across those with coordinate `ax`
#     #  grater than `off`
#     Pselection = findall(x -> x > off, P[ax, :])
#
#     for d = 1 : dim
# 		idxbase = AlphaStructures.findClosestPoint(Psimplex, P[:, Pselection])
# 		@assert !isnothing(idxbase) "firstDeWallSimplex:
# 			not able to determine first Delaunay Simplex"
#         newidx = Pselection[idxbase]
#         indices = [indices; newidx]
#         Psimplex = [Psimplex P[:, newidx]]
#         Pselection = [i for i = 1 : n if i ∉ indices]
#     end
#
#     # Correctness check
# 	radius, center = AlphaStructures.findRadius(Psimplex, true)
#     for i = 1 : n
# 		@assert LinearAlgebra.norm(center - P[:, i]) >= radius "firstDeWallSimplex:
# 			Unable to find first Simplex."
# 	end
#
# 	if DEBUG println("First Simplex = $indices") end
#
# 	return indices
# end
#
# #-------------------------------------------------------------------------------
#
# """
# 	recursiveDelaunayWall(
# 		P::Matrix,
# 		Pblack::Matrix,
# 		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# 		AFL::Array{Array{Int64,1},1},
# 		ax::Int64,
# 		off::Float64,
# 		positive::Bool;
# 		DEBUG = false
# 	)::Array{Array{Int64,1},1}
#
# Utility function that prepeares the Divide phase for Delaunay Wall.
# Returns the Delaunay Triangulation for the positve or negative subspace of `P`
# (according to `positive`) determined by the hyperplane with normal `ax` and
# constant term `off`.
# If the keyword argument `DEBUG` is set to true than all the procedure is shown.
# """
# function recursiveDelaunayWall(
# 		P::Matrix,
# 		Pblack::Array{Float64},
# 		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# 		AFL::Array{Array{Int64,1},1},
# 		ax::Int64,
# 		off::Float64,
# 		positive::Bool;
# 		DEBUG = false
# 	)::Array{Array{Int64,1},1}
#
# 	#DEBUG = true
#
# 	dim, n = size(P)
# 	newaxis = mod(ax, dim) + 1
#
# 	if DEBUG println("Divide Plus/Minus $positive") end
#
# 	Psubset = findall(x -> (x > off) == positive, P[ax, :])
# 	blacklist = setdiff(unique([(keys(tetraDict)...)...]), Psubset)
# 	if !isempty(Pblack)
# 		Pblack = [Pblack P[:, blacklist]]
# 	else
# 	 	Pblack = P[:, blacklist]
# 	end
#
# 	if DEBUG println("Step In") end
#
# 	DT = AlphaStructures.delaunayWall(
# 			P[:, Psubset],
# 			newaxis,
# 			Pblack,
# 			[[findall(Psubset.==p)[1] for p in σ] for σ in AFL],
# 			Dict([
# 				[findall(Psubset.==p)[1] for p in k] => v
# 				for (k,v) in tetraDict
# 					if k ⊆ Psubset
# 			]),
# 			DEBUG = DEBUG
# 		)
#
# 	if DEBUG @show "Step Out with " DT end
#
# 	return [[Psubset[i] for i in σ] for σ in DT]
# end
#
# #-------------------------------------------------------------------------------
#
# """
# 	updateAFL!(
# 		P::Matrix
# 		new::Array{Int64,1}[],
# 		AFLα = Array{Int64,1}[],
# 		AFLplus = Array{Int64,1}[],
# 		AFLminus = Array{Int64,1}[],
# 		ax::Int64, off::Float64;
# 		DEBUG = false
# 	)::Bool
#
# Modify the `AFL*` lists of faces by adding to them the faces inside `new`
# (that refers to the points `P`) according to their position with respect
# to the axis defined by the normal direction `ax` and the contant term `off`.
# The function returns a Bool value that states if the operation was succesfully.
# If the keyword argument `DEBUG` is set to true than all the procedure is shown.
# """
# function updateAFL!(
# 		P::Matrix,
# 		newσ::Array{Array{Int64,1},1},
# 		AFLα::Array{Array{Int64,1},1},
# 		AFLplus::Array{Array{Int64,1},1},
# 		AFLminus::Array{Array{Int64,1},1},
# 		ax::Int64, off::Float64;
# 		DEBUG = false
# 	)::Bool
#
# 	for face in newσ
# 		inters = AlphaStructures.planarIntersection(P, face, ax, off)
#     	if inters == 0 # intersected by plane α
# 			AlphaStructures.updatelist!(AFLα, face)
# 		elseif inters == -1 # in NegHalfspace(α)
#         	AlphaStructures.updatelist!(AFLminus, face)
#     	elseif inters == 1 # in PosHalfspace(α)
#         	AlphaStructures.updatelist!(AFLplus, face)
#     	else
# 			return false
# 		end
# 	end
#
# 	if DEBUG @show AFLα AFLminus AFLplus end
#
# 	return true
#
# end
#
# #-------------------------------------------------------------------------------
#
# """
# 	updatelist!(list, element)::Bool
#
# If `element` is in `list` then it is removed (returns `false`);
# If `element` is not in `list` then it is added (returns `true`).
#
# # Examples
# ```jldoctest
# julia> list = [[1, 2, 3, 4], [2, 3, 4, 5]]
# 2-element Array{Array{Int64,1},1}:
#  [1, 2, 3, 4]
#  [2, 3, 4, 5]
#
# julia> updatelist!(list, [1, 2, 4, 5])
#  true
#
# julia> list
# 3-element Array{Array{Int64,1},1}:
#  [1, 2, 3, 4]
#  [2, 3, 4, 5]
#  [1, 2, 4, 5]
#
# julia> updatelist!(list, [1, 2, 4, 5])
#  false
#
# julia> list
# 2-element Array{Array{Int64,1},1}:
#  [1, 2, 3, 4]
#  [2, 3, 4, 5]
#
# ```
# """
#
# function updatelist!(list, element)::Bool
# 	if element ∈ list
# 		setdiff!(list, [element])
# 		return false
# 	else
# 		push!(list, element)
# 		return true
# 	end
# end
#
# #-------------------------------------------------------------------------------
#
# """
# 	updateTetraDict!(
# 		P::Matrix,
# 		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# 		AFL::Array{Array{Int64,1},1},
# 		σ::Array{Int64,1}
# 	)::Nothing
#
# Update the content of `tetraDict` by adding the outer points of the faces of `σ`
# in the dictionary.
# """
# function updateTetraDict!(
# 		P::Matrix,
# 		tetraDict::DataStructures.Dict{Array{Int64,1},Array{Float64,1}},
# 		AFL::Array{Array{Int64,1},1},
# 		σ::Array{Int64,1}
# 	)::Nothing
# 	for cell in AFL
# 		point = setdiff(σ, cell)
# 		@assert length(point) == 1 "updateTetraDict!: Error during update of TetraDict $σ, $cell"
# 		tetraDict[ cell ] = P[:, point[1]]
# 	end
# end