text:
export # Convex Hull Functions
       ch_akl_toussaint,
       ch_bykat,
       ch_eddy,
       ch_graham_andrew,
       ch_jarvis,
       ch_melkman,
       convex_hull_2,
       # Convexity Checking
       is_ccw_strongly_convex_2,
       is_cw_strongly_convex_2,
       # Hull Subsequence Functions
       ch_graham_andrew_scan,
       ch_jarvis_march,
       lower_hull_points_2,
       upper_hull_points_2,
       # Extreme Point Functions
       ch_e_point,
       ch_n_point,
       ch_ns_point,
       ch_nswe_point,
       ch_s_point,
       ch_w_point,
       ch_we_point

@doc raw"""
    ch_akl_toussaint(ps::AbstractVector{Point2})
    ch_akl_toussaint(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

This function uses the algorithm of Akl and Toussaint [^sga78] that requires
``O(n \log n)`` time for ``n`` input points.

[^sga78]: S. G. Akl and G. T. Toussaint. A fast convex hull algorithm. *Inform.
          Process. Lett.*, 7(5):219--222, 1978.

See also: [`ch_bykat()`](@ref), [`ch_eddy()`](@ref),
[`ch_graham_andrew()`](@ref), [`ch_jarvis()`](@ref), [`ch_melkman()`](@ref),
[`convex_hull_2()`](@ref)
"""
ch_akl_toussaint(ps::AbstractVector{Point2})

@doc raw"""
    ch_bykat(ps::AbstractVector{Point2})
    ch_bykat(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

This function implements the non-recursive variation of Eddy's algorithm (see
[`ch_eddy()`](@ref)) described in [^kra78]. This algorithm requires ``O(nh)``
time in the worst case for ``n`` input points with ``h`` extreme points.

[^kra78]: K. R. Anderson. A reevaluation of an efficient algorithm for determining
          the convex hull of a finite planar set. *Inform. Process. Lett.*,
          7(1):53--55, 1978.

See also: [`ch_akl_toussaint()`](@ref), [`ch_eddy()`](@ref),
[`ch_graham_andrew()`](@ref), [`ch_jarvis()`](@ref), [`ch_melkman()`](@ref),
[`convex_hull_2()`](@ref)
"""
ch_bykat(ps::AbstractVector{Point2})

@doc raw"""
    ch_eddy(ps::AbstractVector{Point2})
    ch_eddy(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

The function implements Eddy's algorithm [^wfe77], which is the two-dimensional
version of the quickhull algorithm [^cbb77].

This algorithm requires ``O(nh)`` time in the worst case for ``n`` input points
with ``h`` extreme points.

[^wfe77]: W. F. Eddy. A new convex hull algorithm for planar sets. *ACM Trans.
          Math. Softw.*, 3:398--403 and 411--412, 1977.

[^cbb96]: C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa. The
          Quickhull algorithm for convex hulls. *ACM Trans. Math. Softw.*,
          22(4):469--483, December 1996.

See also: [`ch_akl_toussaint()`](@ref), [`ch_bykat()`](@ref),
[`ch_graham_andrew()`](@ref), [`ch_jarvis()`](@ref), [`ch_melkman()`](@ref),
[`convex_hull_2()`](@ref)
"""
ch_eddy(ps::AbstractVector{Point2})

@doc raw"""
    ch_graham_andrew(ps::AbstractVector{Point2})
    ch_graham_andrew(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

This function implements Andrew's variant of the Graham scan algorithm
[^ama79] and follows the presentation of Mehlhorn [^km84]. This algorithm
requires ``O(n \log n)`` time in the worst case for ``n`` input points.

[^ama79]: A. M. Andrew. Another efficient algorithm for convex hulls in two
    dimensions. *Inform. Process. Lett.*, 9(5):216--219, 1979.

[^km84]: Kurt Mehlhorn. *Data Structures and Algorithms 3: Multi-dimensional
         Searching and Computational Geometry*, volume 3 of *EATCS Monographs
         on Theoretical Computer Science*. Springer-Verlag, Heidelberg,
         Germany, 1984.

See also: [`ch_akl_toussaint()`](@ref), [`ch_bykat()`](@ref),
[`ch_eddy()`](@ref), [`ch_graham_andrew_scan()`](@ref), [`ch_jarvis()`](@ref),
[`ch_melkman()`](@ref), [`convex_hull_2()`](@ref),
[`lower_hull_points_2`](@ref), [`upper_hull_points_2`](@ref)
"""
ch_graham_andrew(ps::AbstractVector{Point2})

@doc raw"""
    ch_jarvis(ps::AbstractVector{Point2})
    ch_jarvis(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

This functions uses the Jarvis march (gift-wrapping) algorithm [^raj73]. This
algorithm requires ``O(nh)`` time in the worst case for ``n`` input points with
``h`` extreme points.

[^raj73]: R. A. Jarvis. On the identification of the convex hull of a finite
          set of points in the plane. *Inform. Process. Lett.*, 2:18--21, 1973.

See also: [`ch_akl_toussaint()`](@ref), [`ch_bykat()`](@ref),
[`ch_eddy()`](@ref), [`ch_graham_andrew()`](@ref), [`ch_jarvis_march()`](@ref),
[`ch_melkman()`](@ref), [`convex_hull_2()`](@ref)
"""
ch_jarvis(ps::AbstractVector{Point2})

@doc raw"""
    ch_melkman(ps::AbstractVector{Point2})
    ch_melkman(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

It uses an implementation of Melkman's algorithm [^am87]. Running time of
this is linear.

[^am87] A. Melkman. On-line construction of the convex hull of a simple
        polyline. *Inform. Process. Lett.*, 25:11--12, 1987.

See also: [`ch_akl_toussaint()`](@ref), [`ch_bykat()`](@ref),
[`ch_eddy()`](@ref), [`ch_graham_andrew()`](@ref), [`ch_jarvis()`](@ref),
[`convex_hull_2()`](@ref)
"""
ch_melkman(ps::AbstractVector{Point2})

@doc raw"""
    convex_hull_2(ps::AbstractVector{Point2})
    convex_hull_2(ps::Point2...)

Returns the counterclockwise sequence of extreme points of points in the vector
`ps`.

It is not specified at which point the cyclic sequence of extreme points is cut
into a linear sequence.

!!! info "Precondition"

    The vector `ps` does not contain the result.

### Implementation

The ``O(n \log n)`` Akl and Toussaint [^sga78] algorithm is used.

See also: [`ch_akl_toussaint()`](@ref), [`ch_bykat()`](@ref),
[`ch_eddy()`](@ref), [`ch_graham_andrew()`](@ref), [`ch_jarvis()`](@ref),
[`ch_melkman()`](@ref)
"""
convex_hull_2(ps::AbstractVector{Point2})

"""
    is_ccw_strongly_convex_2(ps::AbstractVector{Point2})
    is_ccw_strongly_convex_2(ps::Point2...)

The function [`is_ccw_strongly_convex_2()`](@ref) determines if a given
sequence of points defines a counterclockwise-oriented, strongly convex polygon.

It returns `true`, iff the point elements in `ps` form a
counterclockwise-oriented strongly convex polygon.

A set of points is said to be strongly convex if it consists of only extreme
points (*i.e.*, vertices of the convex hull).

### Implementation

The algorithm requires ``O(n)`` time for a set of ``n`` input points.

See also: [`is_cw_strongly_convex_2()`](@ref)
"""
is_ccw_strongly_convex_2(ps::AbstractVector{Point2})

"""
    is_cw_strongly_convex_2(ps::AbstractVector{Point2})
    is_cw_strongly_convex_2(ps::Point2...)

The function [`is_cw_strongly_convex_2()`](@ref) determines if a given sequence
of points defines a clockwise-oriented, strongly convex polygon.

It returns `true`, iff the point elements in `ps` form a clockwise-oriented
strongly convex polygon.

A set of points is said to be strongly convex if it consists of only extreme
points (*i.e.*, vertices of the convex hull).

### Implementation

The algorithm requires ``O(n)`` time for a set of ``n`` input points.

See also: [`is_ccw_strongly_convex_2()`](@ref)
"""
is_cw_strongly_convex_2(ps::AbstractVector{Point2})

for F ∈ (:is_ccw_strongly_convex_2, :is_cw_strongly_convex_2)
    @eval begin
        $F(ps::AbstractVector) = $F(collect(CxxRef{Point2}, CxxRef.(ps)))
        $F(ps::reference_type_union(Point2)...) = $F(collect(CxxRef.(ps)))
    end
end

"""
    ch_graham_andrew_scan(ps::AbstractVector{Point2})
    ch_graham_andrew_scan(p::Point2, q::Point2, ps::Point2...)

Generates the counterclockwise sequence of extreme points from a given
sequence of input points that are not left of the line defined by the first and
last points in this sequence.

More precisely, it generates the counterclockwise sequence of extreme points
from a given sequence of input points that are not left of the line ``pq``
defined by the first (``p``) and last (``q``) points in the sequence (``p`` is
the value of `first(ps)` and ``q`` is the value of `last(ps)`).  The resulting
sequence is returned starting with ``p``; point ``q`` is omitted.

!!! info "Precondition"

    `ps` contains at least two different points.  The points in `ps` are sorted
    with respect to ``pq``, *i.e.*, the sequence of points in `ps` defines a
    counterclockwise polygon, for which the Graham-Sklansky-procedure [^js72]
    works.

### Implementation

This algorithm requires ``O(n)`` time in the worst case for ``n`` input points.

[^js72]: J. Sklansky. Measuring concavity on rectangular mosaic. *IEEE Trans.
         Comput.*, C-21:1355--1364, 1972.

See also: [`ch_graham_andrew()`](@ref), [`lower_hull_points_2()`](@ref),
[`upper_hull_points_2()`](@ref)
"""
ch_graham_andrew_scan(ps::AbstractVector{Point2})

@cxxdereference ch_graham_andrew_scan(ps::AbstractVector) =
    ch_graham_andrew_scan(collect(CxxRef{Point2}, CxxRef.(ps)))
@cxxdereference ch_graham_andrew_scan(p::Point2, q::Point2,
                                      ps::reference_type_union(Point2)...) =
    ch_graham_andrew_scan(CxxRef.([p, q, ps...]))

"""
    ch_jarvis_march(ps::AbstractVector{Point2}, p::Point2, q::Point2)
    ch_jarvis_march(p::Point2, q::Point2, ps::Point2...)

Generates the counterclockwise sequence of extreme points from a given set of
input points that line between two input points.

More precisely, it generates the counterclockwise subsequence of extreme points
between `p` and `q` of the points in `ps`.  Said sequence is returned starting
with point `p`.  The last point generated is the point preceding `q` in the
counterclockwise order of extreme points.

!!! info "Precondition"

    `p` and `q` are extreme points with respect to the points in `ps` and `q`
    belongs to `ps`.

### Implementation

The functions uses the Jarvis march (gift-wrapping) algorithm [^raj73]. This
algorithm requires ``O(nh)`` time in the worst case for ``n`` input points with
``h`` extreme points.

See also: [`ch_jarvis()`](@ref), [`lower_hull_points_2()`](@ref),
[`upper_hull_points_2()`](@ref)
"""
ch_jarvis_march(ps::AbstractVector{Point2}, p::Point2, q::Point2)

@cxxdereference ch_jarvis_march(ps::AbstractVector, p::Point2, q::Point2) =
    ch_jarvis_march(collect(CxxRef{Point2}, CxxRef.(ps)), p, q)
@cxxdereference ch_jarvis_march(p::Point2, q::Point2, ps::Point2...) =
    ch_jarvis_march(collect(ps), p, q)

@doc raw"""
    lower_hull_points_2(ps::AbstractVector{Point2})
    lower_hull_points_2(ps::Point2...)

Generates the counterclockwise sequence of extreme points on the lower hull of a
given set of input points.

More precisely, it generates the counterclockwise sequence of extreme points in
the lower hull of the points in `ps`.  The resulting sequence is returned,
starting with the leftmost point; the rightmost point is not included.  If there
is only one extreme point (*i.e.*, the leftmost and rightmost points are equal),
the extreme point is reported.

!!! info "Precondition"

    `ps` does not contain the result.

The different treatment by [`upper_hull_points_2()`](@ref) of the case that all
points are equal ensures that concatenation of lower and upper hull points gives
the sequence of extreme points.

### Implementation

This function uses Andrew's variant of Graham's scan algorithm [^ama79],
[^km84]. The algorithm has worst-case running time of ``O(n \log n)`` for
``n`` input points.

See also: [`ch_graham_andrew()`](@ref), [`ch_graham_andrew_scan()`](@ref),
[`upper_hull_points_2()`](@ref)
"""
lower_hull_points_2(ps::AbstractVector{Point2})

@doc raw"""
    upper_hull_points_2(ps::AbstractVector{Point2})
    upper_hull_points_2(ps::Point2...)

Generates the counterclockwise sequence of extreme points on the upper hull of a
given set of input points.

More precisely, it generates the counterclockwise sequence of extreme points in
the lower hull of the points in `ps`.  The resulting sequence is returned,
starting with the rightmost point; the leftmost point is not included.  If there
is only one extreme point (*i.e.*, the leftmost and rightmost points are equal),
the extreme point is reported.

!!! info "Precondition"

    `ps` does not contain the result.

The different treatment by [`lower_hull_points_2()`](@ref) of the case that all
points are equal ensures that concatenation of lower and upper hull points gives
the sequence of extreme points.

### Implementation

This function uses Andrew's variant of Graham's scan algorithm [^ama79],
[^km84]. The algorithm has worst-case running time of ``O(n \log n)`` for
``n`` input points.

See also: [`ch_graham_andrew()`](@ref), [`ch_graham_andrew_scan()`](@ref),
[`lower_hull_points_2()`](@ref)
"""
upper_hull_points_2(ps::AbstractVector{Point2})

for F ∈ (:ch_akl_toussaint
       , :ch_bykat
       , :ch_eddy
       , :ch_graham_andrew
       , :ch_jarvis
       , :ch_melkman
       , :convex_hull_2
       , :lower_hull_points_2
       , :upper_hull_points_2)
    @eval begin
        $F(ps::AbstractVector) = $F(collect(CxxRef{Point2}, CxxRef.(ps)))
        $F(ps::reference_type_union(Point2)...) = $F(collect(CxxRef.(ps)))
    end
end

@doc raw"""
    ch_e_point(ps::AbstractVector{Point2})
    ch_e_point(ps::Point2...)

The function [`ch_e_point()`](@ref) finds a point of a given set of input points
with maximal ``x`` coordinate.

It traverses `ps`. After execution, the returned value is the first point
``e`` such that ``∀p\! ∈\! ps\; e ≥_{xy} p``.

See also: [`ch_n_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_ns_point()`](@ref), [`ch_s_point()`](@ref), [`ch_w_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_e_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_n_point(ps::AbstractVector{Point2})
    ch_n_point(ps::Point2...)

The function [`ch_n_point()`](@ref) finds a point of a given set of input points
with maximal ``y`` coordinate.

It traverses `ps`. After execution, the returned value is the first point
``n`` such that ``∀p\! ∈\! ps\; n ≥_{yx} p``.

See also: [`ch_e_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_ns_point()`](@ref), [`ch_s_point()`](@ref), [`ch_w_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_n_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_ns_point(ps::AbstractVector{Point2})
    ch_ns_point(ps::Point2...)

The function [`ch_ns_point()`](@ref) finds two points of a given set of input
points with minimal and maximal ``y`` coordinate.

It traverses `ps`. After execution, the returned value is a tuple where the
first point ``n`` is such that ``∀p\! ∈\! ps\; n ≥_{yx} p``, and the second
point ``s`` is such that ``∀p\! ∈\! ps\; s ≤_{yx} p``.

See also: [`ch_e_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_n_point()`](@ref), [`ch_s_point()`](@ref), [`ch_w_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_ns_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_nswe_point(ps::AbstractVector{Point2})
    ch_nswe_point(ps::Point2...)

The function [`ch_nswe_point()`](@ref) finds the four extreme points of a given
set of input points using a linear scan of the input points.

That is, it determines the points with maximal ``y``, minimal ``y``, minimal
``x``, and maximal ``x`` coordinates.

It traverses `ps`. After execution, the returned value is a tuple where the
first point ``n`` is such that ``∀p\! ∈\! ps\; n ≥_{yx} p``, the second
point ``s`` is such that ``∀p\! ∈\! ps\; s ≤_{yx} p``, the third point ``w`` is
such that ``∀p\! ∈\! ps\; w ≤_{xy} p``, and the fourth point ``e`` is such that
``∀p\! ∈\! ps\; e ≥_{xy} p``.

See also: [`ch_e_point()`](@ref), [`ch_n_point()`](@ref),
[`ch_ns_point()`](@ref), [`ch_s_point()`](@ref), [`ch_w_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_nswe_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_s_point(ps::AbstractVector{Point2})
    ch_s_point(ps::Point2...)

The function [`ch_s_point()`](@ref) finds a point of a given set of input points
with minimal ``y`` coordinate.

It traverses `ps`. After execution, the returned value is the first point
``s`` such that ``∀p\! ∈\! ps\; s ≤_{yx} p``.

See also: [`ch_e_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_n_point()`](@ref), [`ch_ns_point()`](@ref), [`ch_w_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_s_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_w_point(ps::AbstractVector{Point2})
    ch_w_point(ps::Point2...)

The function [`ch_w_point()`](@ref) finds a point of a given set of input points
with minimal ``x`` coordinate.

It traverses `ps`. After execution, the returned value is the first point
``w`` such that ``∀p\! ∈\! ps\; w ≤_{xy} p``.

See also: [`ch_e_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_n_point()`](@ref), [`ch_ns_point()`](@ref), [`ch_s_point()`](@ref),
[`ch_we_point()`](@ref)
"""
ch_w_point(ps::AbstractVector{Point2})

@doc raw"""
    ch_we_point(ps::AbstractVector{Point2})
    ch_we_point(ps::Point2...)

The function [`ch_we_point()`](@ref) finds two points of a given set of input
points with minimal and maximal ``x`` coordinate.

It traverses `ps`. After execution, the returned value is a tuple where the
first point ``w`` is such that ``∀p\! ∈\! ps\; w ≤_{xy} p``, and the second
point ``e`` is such that ``∀p\! ∈\! ps\; e ≥_{xy} p``.

See also: [`ch_e_point()`](@ref), [`ch_nswe_point()`](@ref),
[`ch_n_point()`](@ref), [`ch_ns_point()`](@ref), [`ch_s_point()`](@ref),
[`ch_w_point()`](@ref)
"""
ch_we_point(ps::AbstractVector{Point2})

for S ∈ (:e, :n, :ns, :nswe, :s, :w, :we)
    F = Symbol(:ch_, S, :_point)
    @eval begin
        $F(ps::AbstractVector) =
            isempty(ps) ?
                $F([CxxRef(Point2())]) :
                $F(CxxRef.(ps))
        $F(ps::reference_type_union(Point2)...) = $F(collect(CxxRef.(ps)))
    end
end


text:
# Index manipulations
#---------------------
"""
    permute(tsrc::AbstractTensorMap{S}, p1::NTuple{N₁, Int}, p2::NTuple{N₂, Int} = ())
        -> tdst::TensorMap{S, N₁, N₂}

Permute the indices of `tsrc::AbstractTensorMap{S}` such that a new tensor
`tdst::TensorMap{S, N₁, N₂}` is obtained, with indices in `p1` playing the role of the
codomain or range of the map, and indices in `p2` indicating the domain.

To permute into an existing `tdst`, see [`add!`](@ref)
"""
function permute(t::TensorMap{S},
                    p1::IndexTuple{N₁},  p2::IndexTuple{N₂}=();
                    copy::Bool = false) where {S, N₁, N₂}
    cod = ProductSpace{S, N₁}(map(n->space(t, n), p1))
    dom = ProductSpace{S, N₂}(map(n->dual(space(t, n)), p2))
    # share data if possible
    if !copy
        if p1 === codomainind(t) && p2 === domainind(t)
            return t
        elseif has_shared_permute(t, p1, p2)
            return TensorMap(reshape(t.data, dim(cod), dim(dom)), cod, dom)
        end
    end
    # general case
    @inbounds begin
        return add!(true, t, false, similar(t, cod←dom), p1, p2)
    end
end

function permute(t::AdjointTensorMap{S}, p1::IndexTuple, p2::IndexTuple=();
                    copy::Bool = false) where {S}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(permute(adjoint(t), p1′, p2′; copy = copy))
end

function has_shared_permute(t::TensorMap, p1, p2)
    if p1 === codomainind(t) && p2 === domainind(t)
        return true
    elseif sectortype(t) === Trivial
        stridet = i->stride(t[], i)
        sizet = i->size(t[], i)
        canfuse1, d1, s1 = TO._canfuse(sizet.(p1), stridet.(p1))
        canfuse2, d2, s2 = TO._canfuse(sizet.(p2), stridet.(p2))
        return canfuse1 && canfuse2 && s1 == 1 && (d2 == 1 || s2 == d1)
    else
        return false
    end
end

function has_shared_permute(t::AdjointTensorMap, p1, p2)
    p1′ = adjointtensorindices(t, p2)
    p2′ = adjointtensorindices(t, p1)
    return has_shared_permute(t', p1′, p2′)
end

"""
    permute!(tdst::AbstractTensorMap{S, N₁, N₂}, tsrc::AbstractTensorMap{S},
        p1::IndexTuple{N₁}, p2::IndexTuple{N₂}=()
    ) -> tdst

Permute the indices of `tsrc` and write the result into `tdst`, with indices in `p1` playing
the role of the codomain or range of the map `tdst`, and indices in `p2` indicating the
domain of `tdst`.

For more details see [`add!`](@ref), of which this is a special case.
"""
@propagate_inbounds Base.permute!(tdst::AbstractTensorMap{S, N₁, N₂},
                                    tsrc::AbstractTensorMap{S},
                                    p1::IndexTuple{N₁},
                                    p2::IndexTuple{N₂}=()) where {S, N₁, N₂} =
    add_permute!(true, tsrc, false, tdst, p1, p2)

# Braid
function braid(t::TensorMap{S}, levels::IndexTuple,
                    p1::IndexTuple, p2::IndexTuple=();
                    copy::Bool = false) where {S}
    @assert length(levels) == numind(t)
    if BraidingStyle(sectortype(S)) isa SymmetricBraiding
        return permute(t, p1, p2; copy = copy)
    end
    if !copy && p1 == codomainind(t) && p2 == domainind(t)
        return t
    end
    # general case
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    @inbounds begin
        return add_braid!(true, t, false, similar(t, cod←dom), p1, p2, levels)
    end
end
@propagate_inbounds braid!(tdst::AbstractTensorMap{S, N₁, N₂},
                                tsrc::AbstractTensorMap{S},
                                levels::IndexTuple,
                                p1::IndexTuple{N₁},
                                p2::IndexTuple{N₂}=()) where {S, N₁, N₂} =
    add!(true, tsrc, false, tdst, p1, p2, levels)

# Transpose
LinearAlgebra.transpose!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap,
                            p1::IndexTuple = reverse(domainind(tsrc)),
                            p2::IndexTuple = reverse(codomainind(tsrc))) =
    add_transpose!(true, tsrc, false, tdst, p1, p2)

function LinearAlgebra.transpose(t::TensorMap{S},
                                    p1::IndexTuple = reverse(domainind(t)),
                                    p2::IndexTuple = reverse(codomainind(t));
                                    copy::Bool = false) where {S}
    if sectortype(S) === Trivial
        return permute(t, p1, p2; copy = copy)
    end
    if !copy && p1 == codomainind(t) && p2 == domainind(t)
        return t
    end
    # general case
    cod = ProductSpace{S}(map(n->space(t, n), p1))
    dom = ProductSpace{S}(map(n->dual(space(t, n)), p2))
    @inbounds begin
        return add_transpose!(true, t, false, similar(t, cod←dom), p1, p2)
    end
end

function LinearAlgebra.transpose(t::AdjointTensorMap{S},
                                    p1::IndexTuple = reverse(domainind(t)),
                                    p2::IndexTuple = reverse(codomainind(t));
                                    copy::Bool = false) where {S}
    p1′ = map(n->adjointtensorindex(t, n), p2)
    p2′ = map(n->adjointtensorindex(t, n), p1)
    adjoint(transpose(adjoint(t), p1′, p2′; copy = copy))
end

# Twist
twist(t::AbstractTensorMap, i::Int; inv::Bool = false) = twist!(copy(t), i; inv = inv)

function twist!(t::AbstractTensorMap, i::Int; inv::Bool = false)
    if i > numind(t)
        msg = "Can't twist index $i of a tensor with only $(numind(t)) indices."
        throw(ArgumentError(msg))
    end
    BraidingStyle(sectortype(t)) == Bosonic() && return t
    N₁ = numout(t)
    for (f1, f2) in fusiontrees(t)
        θ = i <= N₁ ? twist(f1.uncoupled[i]) : twist(f2.uncoupled[i-N₁])
        inv && (θ = θ')
        rmul!(t[f1, f2], θ)
    end
    return t
end

# Fusing and splitting

ext: jl


text:
using BayesOpt
using Test
include("branin.jl")


text:
module JSONWebTokens

#=
Algorithms

HS256   HMAC using SHA-256 hash algorithm
HS384   HMAC using SHA-384 hash algorithm
HS512   HMAC using SHA-512 hash algorithm
RS256   RSASSA using SHA-256 hash algorithm
RS384   RSASSA using SHA-384 hash algorithm
RS512   RSASSA using SHA-512 hash algorithm
ES256   ECDSA using P-256 curve and SHA-256 hash algorithm
ES384   ECDSA using P-384 curve and SHA-384 hash algorithm
ES512   ECDSA using P-521 curve and SHA-512 hash algorithm
none    No digital signature or MAC value included
=#

import JSON, SHA, MbedTLS
import Base64
using Random

abstract type Encoding end
function alg end
Base.show(io::IO, encoding::Encoding) = print(io, alg(encoding))

include("errors.jl")

include("base64url/Base64URL.jl")

include("jws.jl")
include("none.jl")
include("hs.jl")
include("rs.jl")

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 48d87103-4c23-4144-a121-1e33d2bb87f3
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using MLCourse, MLJ, DataFrames, MLJMultivariateStatsInterface, OpenML,
          LinearAlgebra, Statistics, Random, MLJClusteringInterface, StatsPlots,
          Distributions, Distances
end

# ╔═╡ 7b013132-0ee2-11ec-1dd2-25a9f16f0568
begin
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ da8f8ce7-c3ad-4e1a-bb9e-c7be15646d72
md"# K-means Clustering"

# ╔═╡ 49a1493b-ba68-4e7a-8c3d-47423f7a8204
begin
    iris = DataFrame(OpenML.load(61))
    iris = rename!(iris, replace.(names(iris), "'" => ""))
    coerce!(iris, :class => Multiclass)
end

# ╔═╡ 67b8f27f-b5b5-4dc0-ab51-83ac3dcaa2f7
iris.class

# ╔═╡ 1ebd74bc-4479-41d2-aba4-934cdd2b778a
@df iris cornerplot([:sepallength :sepalwidth :petallength :petalwidth],
                    compact = true)

# ╔═╡ 472f0d57-cf50-4c3d-b251-dcf2f7c121c2
md"Number of clusters $(@bind nclust Slider(2:8, show_value = true))"

# ╔═╡ d1c88a44-be52-4b0e-bc23-cca00d10ffb6
mach1 = fit!(machine(KMeans(k = nclust), select(iris, Not(:class)))); # nclust is defined by the slider value above

# ╔═╡ 4df24b22-47a7-446a-8e27-9f2031a75764
prediction = MLJ.predict(mach1, select(iris, Not(:class)));

# ╔═╡ 0e7f34b9-e24f-447b-840b-e2750d2e778b
let
    p1 = scatter(iris.petallength, iris.petalwidth,
                 legend = false, c = Int.(int(prediction)),
                 title = "prediction", xlabel = "petal length", ylabel = "petal width")
    p2 = scatter(iris.petallength, iris.petalwidth,
                 legend = false, c = Int.(int(iris.class)),
                 title = "ground truth", xlabel = "petal length", ylabel = "petal width")
    plot(p1, p2, layout = (1, 2), size = (700, 500))
end

# ╔═╡ 70b3f1bb-7c47-4bb0-aa17-cda6fdbe0469
md"## KMeans Clustering Applied to MNIST"

# ╔═╡ 8b458fe0-8de4-46a1-b14b-1ef430a3c0f3
Markdown.parse("
Let us try to apply K-Means Clustering to the MNIST data set.

```julia
mnist_x, mnist_y = let df = OpenML.load(554) |> DataFrame
    coerce!(df, :class => Multiclass)
    coerce!(df, Count => MLJ.Continuous)
    select(df, Not(:class)),
    df.class
end

mnist_mach = fit!(machine(KMeans(k = 10), mnist_x));

mnist_pred = predict(mnist_mach)

let
    idxs = [findall(==(k), mnist_pred)[1:15] for k in 1:10]
    img = vcat([[hcat([[reshape(Array(mnist_x[idx, :]), 28, :)' fill(255, 28)]
                       for idx in row]...); fill(255, 15*29)']
                for row in idxs]...)./255
    plot(Gray.(img), showaxis = false, ticks = false,
         xlabel = \"examples\", ylabel = \"cluster 1 to 10\")
end
```
$(MLCourse.embed_figure("mnist_kmeans.png"))

Some clusters (rows in the figure above) seem indeed to have specialised on some
digit classes. This can further be seen by looking at the confusion_matrix.

```julia
confusion_matrix(mnist_pred, mnist_y) |> x -> DataFrame(x.mat, x.labels)
```

$(MLCourse.embed_figure("mnist_kmeans_table.png"))

We see that the ninth cluster (9th row in this table and in the figure above) has indeed specialized on images with class label \"0\" and the eigth cluster has specialized on images with class label \"2\", but e.g. cluster 5 contains similar amounts of 0s, 5s and 6s.
")

# ╔═╡ 12317841-c1f8-4022-970e-8ef613c79b78
md"# Hierarchical Clustering"

# ╔═╡ 20297924-5e09-40d2-9f7e-27e0e1b8a968
hc = MLJ.transform(HierarchicalClustering(linkage = :complete,
                                          metric = Euclidean()),
                     select(iris, Not(:class)));

# ╔═╡ 6a603cb6-7a1b-4a1d-9358-c4c811efa2bb
md"h = $(@bind hclusth Slider(range(minimum(hc.dendrogram.heights), stop = maximum(hc.dendrogram.heights), length = 100), default = maximum(hc.dendrogram.heights)/2))"

# ╔═╡ e266ada3-ba4d-4177-8129-f1220f293c72
let
    pred = hc(h = hclusth, k = 1)
    p1 = scatter(iris.petallength, iris.petalwidth,
                 legend = false, c = pred,
                 title = "prediction", xlabel = "petal length",
                 ylabel = "petal width")
    p2 = plot(hc.dendrogram)
    hline!([hclusth], c = :red, w = 3)
    plot(p2, p1, layout = (1, 2), size = (700, 500))
end

# ╔═╡ 85d574c2-b823-4dcf-b711-efc755e724b7
md"# Exercises
## Conceptual

**Exercise 1.** In this exercise, you will perform clustering manually, with ``K = 2``, on a small example with ``n = 6`` observations and ``p = 2`` features. The observations are as follows.

 ````\begin{array}{c|cc}
 \hline
 \mbox{Obs.} &X_1 &X_2 \cr
 \hline
 1 &1 &4 \cr
 2 &1 &3 \cr
 3 &0 &4 \cr
 4 &5 &1 \cr
 5 &6 &2 \cr
 6 &4 &0 \cr
 \hline
 \end{array}````

(a) Plot the observations.

(b) Randomly assign a cluster label to each observation. You can use `rand(1:2, 6)` to do this, if you do not want to do this manually. Report the cluster labels for each observation.

(c) Compute the centroid for each cluster.

(d) Assign each observation to the centroid to which it is closest, in terms of Euclidean distance. Report the cluster labels for each observation.

(e) Repeat (c) and (d) until the answers obtained stop changing.

(f) In your plot from (a), color the observations according to the clusters labels obtained.

(g) Apply hierarchical clustering with complete linkage and Euclidean metric to the same data set.  Draw the dendrogram with the links at the correct heights.

## Applied

**Exercise 1.** In this exercise, you will generate simulated data and perform K-means clustering on the data.

(a) Generate a simulated data set with 20 observations in each of three classes (i.e. 60 observations total), and 50 variables with the following generator
```julia
function data_generator(d; n = 20)
    cluster1 = rand(MvNormal([0; -d; fill(0, 48)], [3; fill(.5, 49)]), n)
    cluster2 = rand(MvNormal([0; 0; fill(0, 48)], [3; fill(.5, 49)]), n)
    cluster3 = rand(MvNormal([0; d; fill(0, 48)], [3; fill(.5, 49)]), n)
    (data = DataFrame(vcat(cluster1', cluster2', cluster3'), :auto),
     true_labels = [fill(1, n); fill(2, n); fill(3, n)])
end
```

(b) PCA is an unsupervised machine learning method (you will learn more about it next week) that allows to visualize high-dimensional data in lower dimensions. You can use the command
```julia
MLJ.transform(fit!(machine(PCA(maxoutdim = 2), data)), data)
```
to fit an unsupervised machine called `PCA` (with at most 2 output dimensions) to the data and transform the data according to this machine to 2 dimensions. Plot the result. Use a different color to indicate the observations in each of the three classes.

(c) Perform K-means clustering of the observations with ``K = 3``. Repeat K-means clustering for multiple random initializations and keep the best result. *Hint:* for reproducibility you can use `Random.seed!(SEED)`, where SEED is positive integer. How well do the clusters that you obtained in K-means clustering compare to the true class labels? *Hint:* a `confusion_matrix` may be helpful to comare the result to the true labels.

(d) Perform K-means clustering with ``K = 2``. Describe your results.

(e) Now perform K-means clustering with ``K = 4``, and describe your results.

(f) Now perform K-means clustering with ``K = 3`` on the 2-dimensional data you obtained with PCA, rather than on the raw data. Comment on the results.

(g) Perform K-means clustering with ``K = 3`` on the data after scaling each variable to have standard deviation one. How do these results compare to those obtained in (c)? Explain.

(h) *Optional:* Create other artificial data sets where scaling does not affect the result of clustering or where it improves the result of clustering.

**Exercise 2**

On the website of the text book, there is a gene expression data set that consists of 40 tissue samples with measurements on 1000 genes. Some samples are from healthy patients, while others are from a diseased group.

(a) Download the data set with
```julia
CSV.read(download(\"https://www.statlearning.com/s/Ch12Ex13.csv\"),
         DataFrame, header = false, transpose = true)
```

(b) Apply hierarchical clustering and plot the result.

(c) Do the results change qualitatively with different types of linkage?

(d) How many patients do you think are in each group?

"

# ╔═╡ 15830699-57c5-4bc2-bc92-54105597ab26
MLCourse.list_notebooks(@__FILE__)

# ╔═╡ 35ac2056-ab72-44b0-9972-723a0887a622
MLCourse.footer()

# ╔═╡ Cell order:
# ╠═48d87103-4c23-4144-a121-1e33d2bb87f3
# ╟─da8f8ce7-c3ad-4e1a-bb9e-c7be15646d72
# ╠═49a1493b-ba68-4e7a-8c3d-47423f7a8204
# ╠═67b8f27f-b5b5-4dc0-ab51-83ac3dcaa2f7
# ╟─1ebd74bc-4479-41d2-aba4-934cdd2b778a
# ╟─472f0d57-cf50-4c3d-b251-dcf2f7c121c2
# ╠═d1c88a44-be52-4b0e-bc23-cca00d10ffb6
# ╠═4df24b22-47a7-446a-8e27-9f2031a75764
# ╟─0e7f34b9-e24f-447b-840b-e2750d2e778b
# ╟─70b3f1bb-7c47-4bb0-aa17-cda6fdbe0469
# ╟─8b458fe0-8de4-46a1-b14b-1ef430a3c0f3
# ╟─12317841-c1f8-4022-970e-8ef613c79b78
# ╠═20297924-5e09-40d2-9f7e-27e0e1b8a968
# ╟─6a603cb6-7a1b-4a1d-9358-c4c811efa2bb
# ╟─e266ada3-ba4d-4177-8129-f1220f293c72
# ╟─85d574c2-b823-4dcf-b711-efc755e724b7
# ╟─15830699-57c5-4bc2-bc92-54105597ab26
# ╟─7b013132-0ee2-11ec-1dd2-25a9f16f0568
# ╟─35ac2056-ab72-44b0-9972-723a0887a622

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
"""
    smooth(x_curr::Gaussian, x_next_smoothed::Gaussian, Ah::AbstractMatrix, Qh::AbstractMatrix)

SMOOTH step of the (extended) Kalman smoother, or (extended) Rauch-Tung-Striebel smoother.
It is implemented in Joseph Form:
```math
m_{n+1}^P = A(h)*m_n
P_{n+1}^P = A(h)*P_n*A(h) + Q(h)

G = P_n * A(h)^T * (P_{n+1}^P)^{-1}
m_n^S = m_n + G * (m_{n+1}^S - m_{n+1}^P)
P_n^S = (I - G*A(h)) P_n (I - G*A(h))^T + G * Q(h) * G + G * P_{n+1}^S * G
```
"""
function smooth(
    x_curr::Gaussian,
    x_next_smoothed::Gaussian,
    Ah::AbstractMatrix,
    Qh::AbstractMatrix,
)
    x_pred = predict(x_curr, Ah, Qh)

    P_p = x_pred.Σ
    P_p_inv = inv(P_p)

    G = x_curr.Σ * Ah' * P_p_inv

    smoothed_mean = x_curr.μ + G * (x_next_smoothed.μ - x_pred.μ)
    smoothed_cov =
        (X_A_Xt(x_curr.Σ, (I - G * Ah)) + X_A_Xt(Qh, G) + X_A_Xt(x_next_smoothed.Σ, G))
    x_curr_smoothed = Gaussian(smoothed_mean, smoothed_cov)
    return x_curr_smoothed, G
end
function smooth(
    x_curr::SRGaussian,
    x_next_smoothed::SRGaussian,
    Ah::AbstractMatrix,
    Qh::SRMatrix,
)
    x_pred = predict(x_curr, Ah, Qh)

    P_p = x_pred.Σ
    P_p_inv = inv(P_p)

    G = x_curr.Σ * Ah' * P_p_inv

    smoothed_mean = x_curr.μ + G * (x_next_smoothed.μ - x_pred.μ)

    _R = [
        x_curr.Σ.squareroot' * (I - G * Ah)'
        Qh.squareroot' * G'
        x_next_smoothed.Σ.squareroot' * G'
    ]
    _, P_s_R = qr(_R)
    smoothed_cov = SRMatrix(P_s_R')

    x_curr_smoothed = Gaussian(smoothed_mean, smoothed_cov)
    return x_curr_smoothed, G
end

function smooth!(x_curr, x_next, Ah, Qh, integ, diffusion=1)
    # x_curr is the state at time t_n (filter estimate) that we want to smooth
    # x_next is the state at time t_{n+1}, already smoothed, which we use for smoothing
    @unpack d, q = integ.cache
    @unpack x_pred = integ.cache
    @unpack C1, G1, G2, C2 = integ.cache

    # Prediction: t -> t+1
    predict_mean!(x_pred, x_curr, Ah)
    predict_cov!(x_pred, x_curr, Ah, Qh, C1, diffusion)

    # Smoothing
    # G = x_curr.Σ * Ah' * P_p_inv
    P_p_chol = Cholesky(x_pred.Σ.squareroot, :L, 0)
    G = rdiv!(_matmul!(G1, x_curr.Σ.mat, Ah'), P_p_chol)

    x_curr.μ .+= G * (x_next.μ .- x_pred.μ)

    # Joseph-Form:
    M, L = C2.mat, C2.squareroot
    D = length(x_pred.μ)

    _matmul!(G2, G, Ah)
    copy!(view(L, 1:D, 1:D), x_curr.Σ.squareroot)
    _matmul!(view(L, 1:D, 1:D), G2, x_curr.Σ.squareroot, -1.0, 1.0)

    _matmul!(view(L, 1:D, D+1:2D), _matmul!(G2, G, sqrt.(diffusion)), Qh.squareroot)
    _matmul!(view(L, 1:D, 2D+1:3D), G, x_next.Σ.squareroot)

    # _matmul!(M, L, L')
    # chol = cholesky!(Symmetric(M), check=false)

    QL =
        false && issuccess(chol) ? Matrix(chol.U)' :
        eltype(L) <: Union{Float16,Float32,Float64} ? lq!(L).L : qr(L').R'
    copy!(x_curr.Σ.squareroot, QL)
    _matmul!(x_curr.Σ.mat, QL, QL')

    return nothing
end


text:
using BinaryBuilder

# Collection of pre-build Bitwarden CLI binaries
name = "bitwarden_cli"
version = v"1.17.1"

url_prefix = "https://github.com/bitwarden/cli/releases/download/v$(version)"
sources = [
    ArchiveSource("$(url_prefix)/bw-linux-$(version).zip", "4704297B438041D39774AA7B077DB72A184A50223FAAE906D2C238D14E2056E9"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)/bw-macos-$(version).zip", "9d5c5a997c73b84aeb43db4c7be93d3fa6443f83ade35a4953b0f1c6862c00c2"; unpack_target = "x86_64-apple-darwin14"),
    ArchiveSource("$(url_prefix)/bw-windows-$(version).zip", "38FE9F5126BC723FB3C0FC00DC15B013826030D6A9F54539B18DBB56EB6FB5EE"; unpack_target = "x86_64-w64-mingw32"),
    FileSource("https://raw.githubusercontent.com/bitwarden/cli/v$(version)/LICENSE.txt", "b98fbb37db5b23bc5cfdcd16793206a5a7120a7b01f75374e5e0888376e4691c")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p "${bindir}"
cp ${target}/bw${exeext} ${bindir}
chmod +x ${bindir}/*
install_license LICENSE.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; cxxstring_abi="cxx11"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows")
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("bw", :bw)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("definitions.jl")

function readData(dataFilename::AbstractString,dictionaryFilename::AbstractString,nTrain,nTest)
# read a corpus of one line per document, word indices (with one-based indexing) -  space separated
    f = open(dictionaryFilename);
        dictionary = readlines(f); # an array of strings
    close(f);
    docCounter = 0;
    f = open(dataFilename);
    documents = Array{Document}(0);
    test_documents = Array{Document}(0);
    for ln in eachline(f)
        docCounter+=1;
        splitLn = split(ln);
        wordVector = Array{Int64}(length(splitLn));
        for i = 1:length(splitLn)
            wordVector[i] = parse(Int64,splitLn[i]);
        end
        if docCounter>(nTrain+nTest)
            break;
        elseif docCounter>nTrain && docCounter<=(nTrain+nTest)
            push!(test_documents, Document(wordVector));
        elseif docCounter<=nTrain
            push!(documents, Document(wordVector));
        end
    end
    close(f);
    return documents,test_documents, dictionary;
end

#convert the non-sparse data structure for documents to the sparse data structure (token index, token count pairs)
function wordVectorToSparseCounts(documents::Array{Document,1})
    sparse_documents = Array{SparseDocument}(length(documents))

    for i = 1:length(documents)
        tokensFound = Dict{Int64, Int64}()

        numTokensFound = 0
        for j = 1:length(documents[i].wordVector)
            word = documents[i].wordVector[j];
            if !haskey(tokensFound,word)
                numTokensFound += 1
                tokensFound[word] = numTokensFound;
            end
        end

        sparse_documents[i] = SparseDocument(Int64[], zeros(Int64, numTokensFound), zeros(Int64, numTokensFound), length(documents[i].wordVector))

        for j = 1:length(documents[i].wordVector)
            word = documents[i].wordVector[j];
            token = tokensFound[word]
            sparse_documents[i].tokenIndex[token] = word
            sparse_documents[i].tokenCount[token] += 1;
        end
        sparse_documents[i].wordVector = copy(documents[i].wordVector)
    end
    return sparse_documents;
end

function splitTestSet(test)
#Split test documents into halves to learn theta on, and halves to predict
#on for a document completion task
    testTrain = deepcopy(test);
    testTest = deepcopy(test);
    for i = 1:length(test);
        test[i].wordVector = test[i].wordVector[randperm(length(test[i].wordVector))];
        lastTrain = Int64(round(length(test[i].wordVector) ./ 2));
        testTrain[i].wordVector = test[i].wordVector[1:lastTrain];
        testTest[i].wordVector = test[i].wordVector[lastTrain+1:end];
    end
    return testTrain, testTest;
end


text:
## TODO: Remove
const _LSEP = UInt8('\n')
const _EOL = UInt8('}')
const _BOL = UInt8('{')
const _ABOL = UInt8('[')
const _INT_MAX = typemax(Int)

const _RowType = SubArray{UInt8,1,Array{UInt8,1},Tuple{UnitRange{Int}},true}
## Detect space in UInt8
import Base: isspace
@inline Base.isspace(c::UInt8) =
c == 0x20 || 0x09 <= c <= 0x0d || c == 0x85 || c == 0xa0

## Raw row utils
@inline _equaleol(c::UInt8) = (c == 0x0a) # 0x0a == UInt8('\n')
function _findnexteol(buf::Vector{UInt8}, size::Int, index::Int)
    while true
        if _equaleol(@inbounds(buf[index])) || isequal(index, size)
            return index
        else
            index += 1
        end
    end
end

function _findpreveol(buf::Vector{UInt8}, filestart::Int, index::Int)
    while true
        if _equaleol(@inbounds(buf[index])) || isequal(index, filestart)
            return index
        else
            index -= 1
        end
    end
end

function skiprows(buf::Vector{UInt8}, size::Int, n::Int, from::Int)
    isequal(from, size) && (return size)
    for _ in 1:n
        from = _findnexteol(buf, size, from+1)
        isequal(from, size) && (return size)
    end
    return from
end

function partitionrows(nrows::Int, nworkers::Int)
    partitionsize = (nrows + 1) ÷ nworkers
    parts = Vector{UnitRange{Int}}(undef, nworkers)
    parts[1] = 1:partitionsize
    for i in 2:(nworkers-1)
        parts[i] = (((i-1) * partitionsize)+1):(i * partitionsize)
    end
    parts[nworkers] = (last(parts[nworkers-1])+1):nrows
    return parts
end

function rowindex(rows::T, row::Int) where T
    return @inbounds((rows[row]+1):rows[row+1])
end

function indexrows(buf::T, size::Int, nrows::Int, start::Int) where {T}
    rows = Vector{Int}(undef, 1)
    rows[1] = start
    for i in 1:nrows
        start = _findnexteol(buf, size, start+1)
        push!(rows, start)
        isequal(start, size) && (return rows)
    end
    return rows
end

function tindexrows(buf::T, size::Int, nworkers::Int) where {T}
    parts = partitionbuf(buf, size, nworkers)
    nparts = length(parts)
    out = Vector{Task}(undef, nparts)
    for i in 1:nparts
        out[i] = @spawn indexrows(buf, last(parts[i]), typemax(Int), _findnexteol(buf, size, first(parts[i])+1))
    end
    ret = mapfoldl(fetch, vcat, out)
    pushfirst!(ret, 0)
    return ret
end

function partitionbuf(buf::T, size::Int, nworkers::Int) where {T}
    guessize = cld(size, nworkers)
    parts = Vector{UnitRange{Int}}(undef, nworkers)
    parts[1] = 1:_findnexteol(buf, size, guessize)
    for i in 2:(nworkers - 1)
        parts[i] = (last(parts[i-1])):_findnexteol(buf, size, (i * guessize))
    end
    parts[nworkers] = (last(parts[nworkers-1])):size
    return parts[length.(parts) .> 1]
end

## Reading Utils
function _eatwhitespace(row::_RowType)
    for (i, c) in pairs(row)
        !isspace(c) && (return @inbounds(@view(row[i:end])))
    end
    @inbounds(@view(row[end:end]))
end

function parserow(row::_RowType, structtype::Nothing)
    return JSON3.read(_eatwhitespace(row))
end

function parserow(row::_RowType, structtype::DataType)
    return JSON3.read(_eatwhitespace(row), structtype)
end

function parserows(buf::Vector{UInt8}, rows::Vector{Int}, indices, RT, structtype, nworkers::Int)
   if nworkers > 1
       return _tparserrows(buf, rows, indices, RT, structtype, nworkers)
   else
       return _parserows(buf, rows, indices, RT, structtype)
   end
end

function _tparserrows(buf::Vector{UInt8}, rows::Vector{Int}, indices, RT, structtype, nworkers::Int)
    parts = partitionrows(length(indices), nworkers)
    prows = Vector{RT}(undef, length(indices))
    plen = length(parts[1])
    @sync for i in 1:nworkers
        @spawn begin
            crows = indices[parts[i]]
            for (j, row) in pairs(crows)
                rindices = rowindex(rows, row)
                prows[(i-1) * plen + j] = parserow(@inbounds(@view(buf[rindices])), structtype)
            end
        end
    end
    return prows
end

function _parserows(buf::Vector{UInt8}, rows::Vector{Int}, indices, RT,  structtype)
    prows = Vector{RT}(undef, length(indices))
    for (i, row) in pairs(indices)
        rindices = rowindex(rows, row)
        prows[i] = parserow(@inbounds(@view(buf[rindices])), structtype)
    end
    return prows
end

function _tmaterialize(buf::Vector{UInt8}, rows::Vector{Int}, indices, names, structtype, nworkers::Int)
    parts = partitionrows(length(indices), nworkers)
    prows = Vector{NamedTuple{names, T} where T<:Tuple}(undef, length(indices))
    plen = length(parts[1])
    @sync for i in 1:nworkers
        @spawn begin
            crows = indices[parts[i]]
            for (j, row) in pairs(crows)
                rindices = rowindex(rows, row)
                row = parserow(@inbounds(@view(buf[rindices])), structtype)
                rnames = propertynames(row)
                prows[(i-1) * plen + j] = (;zip(rnames, [getproperty(row, n) for n in rnames])...)
            end
        end
    end
    return prows
end

function _materialize(buf::Vector{UInt8}, rows::Vector{Int}, indices, names, structtype)
    prows = Vector{NamedTuple{names, T} where T<:Tuple}(undef, length(indices))
    for (j, row) in pairs(indices)
        rindices = rowindex(rows, row)
        row = parserow(@inbounds(@view(buf[rindices])), structtype)
        rnames = propertynames(row)
        prows[j] = (;zip(rnames, [getproperty(row, n) for n in rnames])...)
    end
    return prows
end

function _filter(buf::Vector{UInt8}, f, rowtype, lineindex::Vector{Int}, structtype)
    prows = Vector{rowtype}(undef, 0)
    for row in 1:length(lineindex)-1
        rindices = rowindex(lineindex, row)
        row = parserow(@inbounds(@view(buf[rindices])), structtype)
        f(row) && push!(prows, row)
    end
    return prows
end

function _findall(lines::LineIndex, f)
    inds = Int[]
    for (i, row) in enumerate(lines)
        f(row) && push!(inds, i)
    end
    return inds
end

function _tfindall(buf::Vector{UInt8}, f, lineindex::Vector{Int}, nworkers::Int, structtype)
    tvec = Vector{Task}(undef, nworkers)
    parts = partitionrows(length(lineindex), nworkers)
    @sync for (i, part) in enumerate(parts)
        tvec[i] = @spawn begin
            out_i = Int[]
            for row in part
                index = rowindex(lineindex, row)
                rp = parserow(@inbounds(@view(buf[index])), structtype)
                f(rp) && push!(out_i, rp)
            end
            out_i
        end
    end
    return mapfoldl(fetch, vcat, tvec)
end

function _findnext(lines, f, i)
    len = length(lines)
    while !f(@inbounds(lines[i])) && i < len
        i += 1
    end
    return f(@inbounds(lines[i])) ? i : nothing
end

function _findprev(lines, f, i)
    while !f(@inbounds(lines[i])) && i > 1
        i -=1
    end
    return f(@inbounds(lines[i])) ? i : nothing
end

function _ftmaterialize(buf::Vector{UInt8}, f::Function, outtype, rows::Vector{Int}, indices, structtype, nworkers::Int)
    parts = partitionrows(length(indices), nworkers)
    prows = Vector{outtype}(undef, length(indices))
    plen = length(parts[1])
    @sync for i in 1:nworkers
        @spawn begin
            crows = indices[parts[i]]
            for (j, row) in pairs(crows)
                rindices = rowindex(rows, row)
                rp = parserow(@inbounds(@view(buf[rindices])), structtype)
                prows[(i-1) * plen + j] = f(rp)
            end
        end
    end
    return prows
end

function _tfilter(buf::Vector{UInt8}, f, rowtype, lineindex::Vector{Int}, structtype, nworkers::Int)
    parts = partitionrows(length(lineindex)-1, nworkers)
    rowbuf = [rowtype[] for _ in 1:nworkers]
    @sync for i in 1:nworkers
        @spawn begin
            for row in parts[i]
                rindices = rowindex(lineindex, row)
                rp = parserow(@inbounds(@view(buf[rindices])), structtype)
                f(rp) && push!(rowbuf[i], rp)
            end
        end
    end
    return foldl(vcat, rowbuf)
end

function _fmaterialize(buf::Vector{UInt8}, f::Function, outtype, rows::Vector{Int}, indices, structtype)
    prows = Vector{outtype}(undef, length(indices))
    for (j, row) in pairs(indices)
        rindices = rowindex(rows, row)
        rp = parserow(@inbounds(@view(buf[rindices])), structtype)
        prows[j] = f(rp)
    end
    return prows
end

function _columnwise(lines, coltypes)
    lookup = Dict{Symbol, Int}()
    columns = Vector{AbstractVector}(undef, 0)
    ncols = 0
    # Allocate known columns
    if !isnothing(coltypes)
        for (k, v) in coltypes
            ncols += 1
            lookup[k] = ncols
            push!(columns, Vector{v}(undef,length(lines)))
        end # for
    end # if
    # Add rows
    for (i, row) in enumerate(lines)
        rnames = propertynames(row)
        # Add values
        for name in rnames
            # Known columns
            if haskey(lookup, name)
                columns[lookup[name]][i] = row[name]
            # Unknonw columns
            else
                ncols += 1
                push!(columns, Vector{Union{Missing, Any}}(undef, length(lines)))
                @inbounds columns[ncols][1:(i-1)] .= missing
                lookup[name] = ncols
                columns[ncols][i] = row[name]
            end # if
        end # for rnames
        # Missing in current row
        mnames = filter(x -> x ∉ rnames, collect(keys(lookup)))
        for name in mnames
           columns[lookup[name]][i] = missing
        end # for mnames
    end # for lines
    fcolnames = sort(collect(keys(lookup)), by = k -> lookup[k])
    if !isnothing(coltypes)
        untyped = filter(x -> x ∉ collect(keys(coltypes)), fcolnames)
    else
        untyped = fcolnames
    end
    for col in untyped
        columns[lookup[col]] = [promote(columns[lookup[col]]...)...]
    end
    return (;zip(fcolnames, columns)...)
end # fun

JSON3.StructType(::Type{CategoricalArrays.CategoricalValue{T, M}} where {T, M}) = JSON3.StringType()


text:
# Example of MLE estimation. The data is really drawn from the
# Logit DGP, so the model is well specified, and the MLE has
# the properties discussed in class. E.g., if you make n very
# large you should see that the estimator is very close to the
# true value of theta used to generate data

using Econometrics
include("LogitDGP.jl")
n = 30 # sample size
theta = [0, 0.5] # true theta for generating data
(y, x) = LogitDGP(n, theta) # generate the data

# now define things for estimation
model = theta -> logit(theta, y, x)
theta = zeros(size(x,2)) # start values for estimation

# Perform the estimation - Make sure that you examine
# the MLE estimation programs so that you see how this works
thetahat, objvalue, V, converged = mleresults(model, theta, "estimate logit model");

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
#
# Copyright 2021, Andrew Lin.
#
# This software is licensed under the MIT License.
# See LICENSE.txt at the root of the project or
# https://opensource.org/licenses/MIT
#
include("generate_catan.jl")
GenerateCatan.print_board(GenerateCatan.generate_board())

max_stars_repo_name: AndrewWasHere/setup-catan

max_issues_repo_name: AndrewWasHere/setup-catan

max_forks_repo_name: AndrewWasHere/setup-catan


text:
"
Return the k'th moment of a truncated normal distribution.
"
function moment(d::Truncated{Normal{T},Continuous}, k::Int) where T
    k == 0 ? one(T) : last(moments(d,k))
end

"
Compute the 1'st to k'th moment of a truncated normal distribution. Uses a recursive formula.
"
function moments(d::Truncated{Normal{T},Continuous}, k::Int) where T
    k == 0 && return 1
    m = Array{T}(undef,k+2) #Array of moments with index i being the i-2's moment
                            #(treating the -1's moment as 0 and 0'ths moment as 1)
    m[1], m[2] = 0, 1
    pars = params(d)
    μ, σ, σ² =pars[1], pars[2], pars[2]^2
    L, U = pars[3], pars[4]
    zL, zU = (L - μ)/σ, (U - μ)/σ
    ϕL, ϕU = pdf.(Normal(),(zL, zU))
    ΦUL = cdf(Normal(),zU) - cdf(Normal(),zL)
    for i in 3:(k+2)
        kk = i-2
        #recursive formula for kk'th moment as a function of the previous two moments (if kk=-1 it uses 0)
        m[i] = (kk-1)*σ²*m[i-2] + μ*m[i-1] - σ*(U^(kk-1)*ϕU - L^(kk-1)*ϕL)/ΦUL
    end
    return m[3:(k+2)]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ScheduledEventsProfile <: SwaggerModel
    terminateNotificationProfile::Any # spec type: Union{ Nothing, TerminateNotificationProfile } # spec name: terminateNotificationProfile

    function ScheduledEventsProfile(;terminateNotificationProfile=nothing)
        o = new()
        validate_property(ScheduledEventsProfile, Symbol("terminateNotificationProfile"), terminateNotificationProfile)
        setfield!(o, Symbol("terminateNotificationProfile"), terminateNotificationProfile)
        o
    end
end # type ScheduledEventsProfile

const _property_map_ScheduledEventsProfile = Dict{Symbol,Symbol}(Symbol("terminateNotificationProfile")=>Symbol("terminateNotificationProfile"))
const _property_types_ScheduledEventsProfile = Dict{Symbol,String}(Symbol("terminateNotificationProfile")=>"TerminateNotificationProfile")
Base.propertynames(::Type{ ScheduledEventsProfile }) = collect(keys(_property_map_ScheduledEventsProfile))
Swagger.property_type(::Type{ ScheduledEventsProfile }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ScheduledEventsProfile[name]))}
Swagger.field_name(::Type{ ScheduledEventsProfile }, property_name::Symbol) =  _property_map_ScheduledEventsProfile[property_name]

function check_required(o::ScheduledEventsProfile)
    true
end

function validate_property(::Type{ ScheduledEventsProfile }, name::Symbol, val)
end


text:
using CombinedUncertainDiffEq
using ComponentArrays
using ConcreteStructs
using DifferentialEquations
using DiffEqUncertainty
using IntervalArithmetic
using LinearAlgebra
using Plots; plotlyjs()
using Statistics
using UnPack

## Functions
IntervalArithmetic.mid(x::Float64) = x

range_width(x) = 0
range_width(x::Interval) = x.hi - x.lo

function soft_step(x; halfwidth=0.1)
    if x < -halfwidth
        zero(x)
    elseif x > halfwidth
        one(x)
    else
        typeof(x)((sinpi(x/halfwidth/2) + 1)/2)
    end
end

function soft_between(x, lo, hi; kwargs...)
    sgn_lo = soft_step(x-lo; kwargs...)
    sgn_hi = soft_step(hi-x; kwargs...)
    return sgn_lo * sgn_hi
end

# ODE function for 2D lander problem
function lander_2D!(D, vars, params, t; T=0.0, θ=0.0)
    @unpack r, v, m_prop = vars
    @unpack g, Isp, m_dry = params

    m = m_dry + m_prop
    Fx = T*sin(θ)
    Fy = T*cos(θ)

    D.r.x = v.x
    D.r.y = v.y
    D.v.x = Fx/m
    D.v.y = -g + Fy/m
    D.m_prop = -T / (Isp * g)

    return nothing
end

# Thrust control
constant_trust(vars, p, t) = (vars.m_prop > 0) * p.control.thrust_mag

# Angle control
function angle_pd(vars, params, t)
    @unpack r, v = vars
    @unpack kp, kd = params.control
    return -kd*v.x - kp*r.x
end

## Build problem
function build_prob()
    # Initial conditions
    ic = ComponentArray{Any}(
        r = (
            x = 0.0,
            y = 120..150,
        ),
        v = (
            x = 0 ± 20,
            y = -45 ± 5,
        ),
        m_prop = 20.0,
    )

    # Parameters
    p = ComponentArray{Any}(;
        g = 9.80665,
        m_dry = 70 ± 3,
        Isp = 60 ± 5,
        control = (
            kp = 0.0, #1,
            kd = 0.05,
            thrust_mag = 1400..1450,
        ),
    )

    # Callback functions
    hits_ground = ContinuousCallback((vars, t, integrator) -> vars.r.y, terminate!)
    burns_out = ContinuousCallback((vars, t, integrator) -> vars.m_prop, identity)

    f! = WithControls(lander_2D!; T=constant_trust, θ=angle_pd)
    cb = CallbackSet(burns_out, hits_ground)
    return ic, p, ODEProblem(f!, mid.(ic), (0.0, 100.0), mid.(p); callback=hits_ground)
end

##
const ic, p, ODE_prob = build_prob();

##
# Requirement functions
const requirements = ComponentArray(
    x = -0.5..0.5,
    y = -10..0,
)
velocity_req(v) = (v[1] ∈ requirements[1]) && (v[2] ∈ requirements[2])
meets_requirements(sol) = velocity_req(sol.u[end].v)
function soft_meets_requirements(sol; halfwidth=0.1)
    v = sol.u[end].v
    return mapreduce(*, v, requirements) do v, req
        Δ = (req.hi - req.lo) * halfwidth
        soft_between(v, req.lo-Δ, req.hi+Δ; halfwidth=Δ)
    end
    return prod(soft_between.(v, first.(requirements).v.-halfwidth, last.(requirements).v.+halfwidth; halfwidth))
end

# widths = [1, 20] #range_width.(requirements)
# g_single = (v,w) -> (v/w)^2
# g = sol -> mapreduce(g_single, +, sol.u[end].v, widths)
# g_vec = (vx,vy) -> g((;u=[(;v=[vx,vy])]))
# plotify_obj(f) = (vx,vy) -> f((;u=[(;v=[vx,vy])]))

## Simulate
@time mc_sol = mc_solve(ODE_prob, ic, p, Tsit5(), EnsembleSerial(); trajectories=10000);
vx = [sol[3, end] for sol in mc_sol.u]
vy = [sol[4, end] for sol in mc_sol.u]

## Solve
# @time mc_exp = mc_expectation(meets_requirements, ODE_prob, ic, p, Tsit5(), EnsembleThreads(); trajectories=100000)
@time koop_exp = koopman_expectation(meets_requirements, ODE_prob, ic, p, Tsit5())

## Plot Simulation
plot(mc_sol, vars=(1,2))
# color = 2 .- Int.(meets_requirements.(mc_sol.u))
# scatter(vx, vy, markercolor=color, legend=false, markersize=1)

##
run_with_trajectories(trajectories) = mc_expectation(exp_fun, ODE_prob, ic, p, Tsit5(), EnsembleThreads(); trajectories)
function running_expectation(trajectories)
    mc_sol = mc_solve(ODE_prob, ic, p, Tsit5(), EnsembleThreads(); trajectories)
    outs = meets_requirements.(mc_sol.u)
    return cumsum(outs) ./ (1:length(outs))
end

## Plot
plot()
# contour!(-1.5:0.1:1.5, -60:0, plotify_obj(exp_fun), colorbar=false, colorbartitle="g(x)")
plot!(reduce(CombinedUncertainDiffEq.:×, requirements), opacity=0, color=3, linecolor=3, linealpha=1, linewidth=2)
scatter!(vx, vy, markersize=1, color=:ivory, markerstrokewidth=0, legend=false, size=(600,400))
xlims!(-1.5, 1.5)
ylims!(-60, 0)
xlabel!("Horizontal Velocity (m/s)")
ylabel!("Vertical Velocity (m/s)")
title!("Velocity Requirements")

## Plot mc vs koopman
i = 200:200:200_000
plot(i, reduce(hcat, running_expectation(i[end])[i] for j in 1:20), label="Monte Carlo", primary=false)

# for j in 1:3
#     plot!(i, running_expectation(i[end])[i], label="mc$j")
# end
hline!([koop_exp.u], lw=2, label="koopman", legend=:bottomright)
# ylims!(0.095, 0.11)
ylims!(0.065, 0.1)
title!("Convergence of Monte Carlo Solution")

##
using DiffEqFlux


text:
module_path = splitdir(@__FILE__)[1]            # Path to this file

# NOTE: In Julia v0.7, JLD throws the error "Cannot `convert` an object of type
# getfield(JLD, Symbol("##JLD.AssociativeWrapper" unless, imported outside of
# WFVisual through the following line. I'm not sure, but I think it's because
# JLD must be imported before GeometricTools to be able to read Grid JLDs.
# Also, notice that `generate_windfarm(...)` takes an abnormal long time to run
# in v0.7.
import JLD

# Load WFVisual module
import WFVisual
wfv=WFVisual

# Load GeometricTools: https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
gt=GeometricTools

# Data path with geometry JLDs
if Int(VERSION.major)==0 && Int(VERSION.minor)==6   # Case of Julia v0.6
    data_path = joinpath(module_path, "../data/")
else                                                # Case of Julia v0.7
    data_path = joinpath(module_path, "../datav07/")
end

# Create save path where to store vtks
save_path = "temps/windfarmgauss00/"         # Save path of this example
gt.create_path(save_path, true)

# --------------------- WIND FARM LAYOUT ---------------------------------------
turbine_x = [  909.98349606,  1005.0174903 ,   900.40238835,  -479.57607866,
                 937.92551703,   461.20472344,   123.06165965, -1073.3529325 ,
                -698.76200523, -1083.53094471,  1365.9338773 ,   204.10557045,
                 306.99725222,  -585.35767689,    10.9706543 ,   543.54927959,
                 -20.50212678, -1313.07762466, -1327.00852816,   137.06316521,
                1378.99410072,  -645.66250771,  -535.7876181 ,  -100.66538419,
               -1331.29278587]

turbine_y = [  500.59889721,  -974.65186327,    76.51586507,  1315.29789208,
                1039.37304649, -1321.85187113,  -300.57817461,  -765.82650713,
               -1213.14956771,   886.54973742,  -306.96866239,   194.51410928,
                 735.69027357,    81.95252777,  -892.14222889,  -930.41798152,
                 970.9117144 ,   398.75277565,    -3.52236011,  1393.27798468,
                 241.61586052,  -791.92492366,   919.26560593, -1396.37638078,
                -433.19849727]
turbine_z = [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
                0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]

hub_height = [ 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
                100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
                100.,  100.,  100.,  100.,  100.,  100.,  100.]

rotor_diameter = [ 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
                100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
                100.,  100.,  100.,  100.,  100.,  100.,  100.]

nBlades = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
               3, 3]

wind_direction = 228.0

yaw = [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.] - wind_direction

# --------------------- PERIMETER AND FLUID DOMAIN -----------------------------
NDIVSx = 70              # Cells in the parametric x-direction
NDIVSy = 70              # Cells in the parametric y-direction
NDIVSz = 70              # Cells in the geometric z-direction

# Circular wind farm perimeter
Rper = norm(maximum.([turbine_x, turbine_y]) - minimum.([turbine_x, turbine_y]))/2*3/4
perimeter_points = Rper*[ [cos(a), sin(a), 0] for a in linspace(0, 2*pi, 179)]

# Guassian wake function
magVinf = 12.0                          # (m/s) wind speed
Vinf_Oaxis = gt.rotation_matrix2(0, 0, wind_direction) # Free-stream coordinate system
invVinf_Oaxis = Vinf_Oaxis'             # Inverse transformation
refH = mean(hub_height)                 # (m) reference height for shear layer

function wake(X)

    k = 0.0325
    CT = 8.0/9.0
    loss = 0
    alpha = 0.12

    # Calculate loss in freestream direction
    for i = 1:length(turbine_x)     # Iterate over turbines

        # Distance from this turbine in global coordinates
        DX = X - [turbine_x[i], turbine_y[i], turbine_z[i]+hub_height[i]]

        # Distance in free-stream coordinates
        dX = invVinf_Oaxis*DX
        dx = dX[1]

        if dx > 0.
            dy = dX[2]
            dz = dX[3]
            sigma = k*dx + rotor_diameter[i]/sqrt(8.0)
            this_loss = (1.0-sqrt(1.0-CT/(8*sigma^2/(rotor_diameter[i]^2))))*
                           exp(-1.0/2.0*(dy/sigma)^2)*exp(-1.0/2.0*(dz/sigma)^2)
            loss += this_loss^2
        end

    end
    loss = sqrt(loss)

    windspeed= magVinf
    # Uncomment this line adds the shear layer
    # windspeed *= (X[3]/refH)^alpha

    # Return velocity vector in global coordinate system
    return Vinf_Oaxis*[ windspeed*(1-loss), 0, 0 ]
end

# NOTE: For some reason Paraview 5.6.0 runs into a glitch rendering the fluid
# domain as a volume. In order to solve that, create a calculator with the
# magnitude of the velocity and apply the volume rendering over the calculator.

# --------------------- GENERATE WIND FARM -------------------------------------

wfv.generate_windfarm(rotor_diameter, hub_height, nBlades,
                                turbine_x, turbine_y, turbine_z,
                                yaw,
                                perimeter_points, wake;
                                NDIVSx=NDIVSx, NDIVSy=NDIVSy, NDIVSz=NDIVSz,
                                save_path=save_path, spl_s=0.01,
                                data_path=data_path);

nothing

max_issues_repo_name: byuflowlab/wind-farm-utilities

max_forks_repo_name: byuflowlab/wind-farm-utilities


text:
import PredictMD
import Test

Test.@test( isdir(PredictMD.package_directory()) )

Test.@test( isdir(PredictMD.package_directory("ci")) )

Test.@test( isdir(PredictMD.package_directory("ci", "travis")) )

Test.@test( isdir(PredictMD.package_directory(TestModuleA)) )

Test.@test( isdir(PredictMD.package_directory(TestModuleB)) )

Test.@test( isdir( PredictMD.package_directory(
            TestModuleB, "directory2",
            ) ) )

Test.@test( isdir( PredictMD.package_directory(
            TestModuleB, "directory2", "directory3",
            ) ) )

Test.@test_throws(
    ErrorException,
    PredictMD.package_directory(TestModuleC),
    )


text:
function getlaydown(h, g::AirDefenseGame, C, A)
    len = size(g.X, 1)
    cities = DataFrame(id = 1:len, x = g.X[:,1], y = g.X[:,2], value = g.v,
        full_coverage = fill(0, len), partial_coverage = fill(0, len),
        iads = fill(false, len), physical_attack = fill(false, len),
        cyber_defense = fill(false, len), cyber_attack = fill(false, len),
        iads_status = fill("none", len))
    for i in g.iads # physical defenses
        cities.iads[i] = true
    end
    for i in A[6][h[6]] # physically attacked cities
        cities.physical_attack[i] = true
    end
    for i in getindex(g.iads, getindex(A[1][h[1]], A[5][h[5]])) # cyber-defended iads
        cities.cyber_defense[i] = true
    end
    for i in getindex(g.iads, getindex(A[2][h[2]], A[3][h[3]])) # cyber-attacked iads
        cities.cyber_attack[i] = true
    end
    for i in 1:len
        if cities.iads[i] == true && cities.cyber_attack[i] == false
            cities.iads_status[i] = "full"
        elseif cities.iads[i] == true && cities.cyber_attack[i] == true && cities.cyber_defense[i] == true
            cities.iads_status[i] = "partial"
        elseif cities.iads[i] == true && cities.cyber_attack[i] == true && cities.cyber_defense[i] == false
            cities.iads_status[i] = "inoperable"
        end
    end
    full_defenders = @where(cities, :iads_status .== "full")
    partial_defenders = @where(cities, :iads_status .== "partial")
    for i in 1:len
        for j in 1:size(full_defenders, 1)
            if C[full_defenders.id[j], cities.id[i]] == 1
                cities.full_coverage[i] += 1
            end
        end
        for pd in 1:size(partial_defenders, 1)
            if C[partial_defenders.id[pd], cities.id[i]] == 1
                cities.partial_coverage[i] += 1
            end
        end
    end
    cities
end

function plot_laydown(df::DataFrame)
    clibrary(:Plots)
    df_dc = @where(df, :cyber_defense .== true)
    df_ac = @where(df, :cyber_attack .== true)
    df_fc = @where(df, :iads_status .== "full")
    df_pc = @where(df, :iads_status .== "partial")
    df_nc = @where(df, :iads_status .== "inoperable")
    df_ap = @where(df, :physical_attack .== true)
    # plot cities
    p = plot(df.x, df.y, seriestype = :scatter, markershape = :circle, markercolor = :black,
        markersize = 1, markeralpha = 0.5, legend = :none,
        xlabel = "Longitude (Regularized)", ylabel = "Latitude (Regularized)",
        label = "Cities")
    # plot cyber defenses
    plot!(p, df_dc.x, df_dc.y, seriestype = :scatter, markershape = :circle,
        markeralpha = 0.2, markercolor = :blue, markersize = 10,
        label = "Cyber Defense")
    # plot cyber attacks
    plot!(p, df_ac.x, df_ac.y, seriestype = :scatter, markershape = :+,
        markeralpha = 0.0, markercolor = :purple, markersize = 7,
        markerstrokealpha = 1.0, markerstrokecolor = :purple, markerstrokewidth = 0.0001,
        markerstrokestyle = :dash, label = "Cyber Attack")
    # label cities with value
    for i in 1:size(df, 1)
        annotate!(p, df.x[i], df.y[i], text(string("  ", round(df.value[i], digits = 2)), 10, :left, :darkgray))
        annotate!(p, df.x[i], df.y[i], text(string(df.id[i], ""), 10, :bottom, :black))
    end
    # plot coverage
    for i in 1:size(df, 1)
        if df.iads[i] == true
            cx, cy = circlecoord(df.x[i], df.y[i], g.r)
            style = df.iads_status[i] == "full" ? :solid : df.iads_status[i] == "partial" ? :dash : :dot
            plot!(p, cx, cy, linestyle = style, linecolor = :black)
        end
    end
    # plot physical attack
    colors = Vector{Symbol}(undef, size(df_ap, 1))
    for i in 1:size(df_ap, 1)
        if df_ap.full_coverage[i] > 0
            colors[i] = :green
        elseif df_ap.partial_coverage[i] > 0
            colors[i] = :orange
        else
            colors[i] = :red
        end
    end
    plot!(p, df_ap.x, df_ap.y, seriestype = :scatter, markershape = :x,
        markeralpha = 0.0, markerstrokecolor = colors, markersize = 7,
        markerstrokealpha = 1.0, markerstrokewidth = 1,
        markerstrokestyle = :dash, label = "Physical Attack")
    p
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
import Printf
import CSV, DataFrames
import CSV

include("exp.jl")

function gen_row(row, d)
    one = [Printf.@sprintf("%.3f / %.3f", item.prec, item.recall) for item in row]
    join(["model-"*string(d), one...], ", ")
end

function transfer_modelX_dY(modelX, dY, Y)
    @info "creating model .."
    tmp_model = Chain(x->reshape(x, Y, Y, 1, :),
                      modelX[2:end-1]...,
                      x->reshape(x, Y, Y, :))
    # FIXME use multile batches
    x, y = next_batch!(dY) |> cpu
    @info "testing .."
    m = sup_graph_metrics(tmp_model(x), y)
    @info "result" Y m.prec, m.recall
    m
end

function transfer_modelX(fname, d5, d10, d15, d20)
    @info "loading model .."
    @load ("trained_model/" * fname) model
    [transfer_modelX_dY(model, d5, 5),
     transfer_modelX_dY(model, d10, 10),
     transfer_modelX_dY(model, d15, 15),
     transfer_modelX_dY(model, d20, 20)]
end

function test_transfer()
    # bson_files = filter(x->occursin("EQ-deep", x), readdir("trained_model"))
    # fd = parse(Int, match(r"d=(\d+)", bson_files[1]).captures[1])
    f5 = "NEW-EQ-deep-d=5-ng=5000-N=20-2020-03-04T23:48:40.762.bson"
    f10 = "NEW-EQ-deep-d=10-ng=5000-N=20-2020-03-05T02:08:38.134.bson"
    f15 = "NEW-EQ-deep-d=15-ng=5000-N=20-2020-03-05T06:13:58.645.bson"
    f20 = "NEW-EQ-deep-d=20-ng=5000-N=20-2020-03-05T13:08:06.271.bson"
    # identify the d=X
    d5, _ = gen_sup_ds_cached(ng=10000, N=20, d=5, batch_size=100)
    d10, _ = gen_sup_ds_cached(ng=10000, N=20, d=10, batch_size=100)
    d15, _ = gen_sup_ds_cached(ng=10000, N=20, d=15, batch_size=100)
    d20, _ = gen_sup_ds_cached(ng=10000, N=20, d=20, batch_size=100)

    row5 = transfer_modelX(f5, d5, d10, d15, d20)
    row10 = transfer_modelX(f10, d5, d10, d15, d20)
    row15 = transfer_modelX(f15, d5, d10, d15, d20)
    row20 = transfer_modelX(f20, d5, d10, d15, d20)

    # generate table
    println(join(["model\\data", "graph-5", "graph-10", "graph-15", "graph-20"], ", "))
    println(gen_row(row5, 5))
    println(gen_row(row10, 10))
    println(gen_row(row15, 15))
    println(gen_row(row20, 20))
end

function test_universal()
    @load "saved_models/EQ-deep-mixed-d=[10,15,20]-ONE/step-100000.bson" model
    @load "saved_models/EQ-deep-mixed-d=[10,15,20]-ONE/step-17671.bson" model
    @load "saved_models/EQ-deep-mixed-d=[10,15,20]-COV/step-50000.bson" model

    x, y = next_batch!(d15)
    sup_graph_metrics(model(x), y)

    # for d in [5,7,10,13,15,17,20,23,25,27,30,35,40]
    # for d in [5,7,10,13,15]
    # for d in [17,20,23,25]
    # for d in [27,30,35,40]
    for d in [7, 10, 15, 20, 25, 30]
        ds, _ = gen_sup_ds_cached(ng=2000, N=2, d=d, batch_size=100)
        prec = MeanMetric{Float64}()
        recall = MeanMetric{Float64}()
        shd = MeanMetric{Float64}()
        @showprogress 0.1 "evaluating .." for i in 1:10
            x, y = next_batch!(ds)
            m = sup_graph_metrics(model(x), y)
            add!(prec, m.prec)
            add!(recall, m.recall)
            add!(shd, m.shd)
        end
        @info "result" d get!(prec) get!(recall) get!(shd)
    end
end

function test_syn()
    g = gen_ER_dag(10)
    W = gen_weights(g, ()->((rand() * 1.5 + 0.5) * rand([1,-1])))
    # (1000, 10)
    Matrix(W)
    X = gen_data(g, W, 10000)
    X2 = gen_data2(W, 10000)
    mean(X, dims=1)
    maximum(cov(X))
    μ
    Σ
    size(X)

    # FIXME these are not equal
    # Σ = (I-W)-1 (I-W)-T     (wrong)
    # Σ = (I-W)-T (I-W)-1
    inv(myeye(10) - W) * inv((myeye(10) - W)')
    SS = inv((myeye(10) - W)') * inv(myeye(10) - W)
    (inv((myeye(10) - W)))' * inv(myeye(10) - W)
    # cor(X)
    cov(X)
    cov(X2)
    var(X, dims=1)

    cor(X .- μ)

    out = reshape(model(reshape(cov(X), size(cov(X))..., 1)), size(cov(X)))
    out = reshape(model(reshape(cor(X), size(cor(X))..., 1)), size(cor(X)))
    out = reshape(model(reshape(SS, size(SS)..., 1)), size(SS))
    Wout = threshold(out, 0.05, true)
    DiGraph(Wout)
    display_graph_with_label(DiGraph(Wout), names(df))
end

function test_save()
    # generate and save as csv file
    # d10, _ = gen_sup_ds_cached(ng=10000, N=20, d=10, batch_size=100)
    # gen a graph
    g = gen_ER_dag(100)
    W = gen_weights(g, ()->((rand() * 1.5 + 0.5) * rand([1,-1])))
    X = gen_data(g, W, 1000)
    # save
    df = DataFrames.DataFrame(X)
    CSV.write("g100.csv", df)
end

function test_notears()
    include("notears.jl")
    Ŵ = notears(X)
end

function test()
    # read d=20 bson
    model_file = "./trained_model/NEW-EQ-deep-d=20-ng=5000-N=20-2020-03-05T13:08:06.271.bson"
    @load model_file model
    # generate d=20 data and monitor results
    ds, test_ds = gen_sup_ds_cached(ng=10000, N=20, d=20, batch_size=100)
    x, y = next_batch!(ds) |> cpu
    model(x)
    # gpu(model)(gpu(x))
    m = sup_graph_metrics(model(x), y)
    m.prec
    m.recall
    # adapt to d=10
    d10_model = Chain(x->reshape(x, 10, 10, 1, :),
                      model[2:end-1]...,
                      x->reshape(x, 10, 10, :))

    # generate d=10 data and monitor results
    ds10, _ = gen_sup_ds_cached(ng=10000, N=20, d=10, batch_size=100)
    x10, y10 = next_batch!(ds10) |> cpu
    m10 = sup_graph_metrics(d10_model(x10), y10)
    m10.prec
    m10.recall

end

using Distributions
using CausalInference
using LightGraphs

using DelimitedFiles, LinearAlgebra, Statistics

function test_PC_demo()
    p = 0.01

    # Download data
    run(`wget http://nugget.unisa.edu.au/ParallelPC/data/real/NCI-60.csv`)

    # Read data and compute correlation matrix
    X = readdlm("NCI-60.csv", ',')
    d, n = size(X)
    C = Symmetric(cor(X, dims=2))

    # Compute skeleton `h` and separating sets `S`
    h, S = skeleton(d, gausscitest, (C, 100), quantile(Normal(), 1-p/2))

    # Compute the CPDAG `g`
    g = pcalg(d, gausscitest, (C, 100), quantile(Normal(), 1-p/2))
end

function test_PC()
    # 1. load data
    ds, test_ds = spec2ds(DataSpec(d=10, k=1, gtype=:ER, noise=:Gaussian))
    # 2. get X
    X, Y = next_batch!(test_ds) |> gpu
    # 3. run PC algorithm
    d, n = size(X)
    C = Symmetric(cor(X, dims=2))

    # Compute skeleton `h` and separating sets `S`
    h, S = skeleton(d, gausscitest, (C, n), quantile(Normal(), 1-p/2))

    # Compute the CPDAG `g`
    g = pcalg(d, gausscitest, (C, n), quantile(Normal(), 1-p/2))

end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
@testset "Test common utility functions" begin
    @test DualOptimalFiltering.normalise(1:4) == (1:4)/10

    @test_nowarn DualOptimalFiltering.lgamma_local(50)
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# The decimal number, 585 = 10010010012 (binary), is palindromic in both bases.
#
# Find the sum of all numbers, less than one million, which are palindromic in
# base 10 and base 2.
#
# (Please note that the palindromic number, in either base, may not include
# leading zeros.)

using ProjectEulerSolutions

# Simple using Julia's builtin bitstring and parse functions (with BigInt).
function p036solution(n::Integer=10)::Integer
    sums = 0
    for i in 1:n
        if ispalindrome_integer(i) && ispalindrome(parse(BigInt, bitstring(i)))
            sums += i
        end
    end
    return sums
end

p036 = Problems.Problem(p036solution)

Problems.benchmark(p036, 1_000_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


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

__precompile__(true)

"""
    ValueShapes

Provides a Julia API to describe the shape of values, like scalars, arrays
and structures.
"""
module ValueShapes

using ArgCheck
using ArraysOfArrays
using Distributions
using ElasticArrays
using FillArrays
using Random
using Statistics
using StatsBase

import ChainRulesCore
import IntervalSets
import Tables
import TypedTables

# Long-term, ChainRulesCore should be sufficient:
import ZygoteRules

using ChainRulesCore: AbstractTangent, Tangent, NoTangent, ZeroTangent, AbstractThunk, unthunk

include("value_shape.jl")
include("value_accessor.jl")
include("scalar_shape.jl")
include("array_shape.jl")
include("const_value_shape.jl")
include("named_tuple_shape.jl")
include("functions.jl")
include("distributions.jl")
include("scalar_dist.jl")
include("const_value_dist.jl")
include("named_tuple_dist.jl")
include("reshaped_dist.jl")

end # module


text:
module WhatIsHappening

using Logging
using MacroTools  # (This is fast).

export @withfeedback, prettify_logging_in_notebook!, get_pretty_notebook_logger

"""
    @withfeedback [message] expression

Print what is happening during `expression`, then execute the expression. If no `message` is
specified, the expression itself is used as message.
"""
macro withfeedback(message::String, expr)
    quote
        @info $message
        _flush_logging_stream()
        $(esc(expr))
    end
end

macro withfeedback(expr)
    :( @withfeedback $(_content_as_string(rmlines(expr))) $(esc(expr)))
end

function _content_as_string(expr)
    if expr.head == :block
        # For a `begin` block, only print contents (not "begin … end").
        return "\n" * join(rmlines(expr).args, " jf\n")
        #    Leading newline, so that first code line comes on its own line after "[ Info:".
    else
        return string(expr)
    end
end

function _flush_logging_stream()
    logger = current_logger()
    # Handle the wrapping loggers of LoggingExtras.jl.
    if hasproperty(logger, :logger)
        logger = logger.logger
    end
    #   [Not done: deeper unwrapping (`while`); and flushing multiple for teed/multiplexed]
    if hasproperty(logger, :stream)
        io = logger.stream
        # Emulate what the default loggers do both in IJulia and in the REPL (one is a
        # SimpleLogger, the other a ConsoleLogger). They both check `isopen` and otherwise
        # use `stderr`.
        if isopen(io)
            flush(io)
        else
            flush(stderr)
        end
    else
        # This is the case in vscode. We could check whether `logger isa VSCodeLogger`, but
        # that adds a dependency; and the defining package is not registered anyway. In any
        # case: no action needed, VSCodeLogger flushes after every message:
        # https://github.com/julia-vscode/julia-vscode/blob/master/scripts/packages/VSCodeServer/src/progress.jl#L7
    end
end

"""Prettier @info display in IJulia than the default (which is noisy, and on a red bg)."""
prettify_logging_in_notebook!() = global_logger(get_pretty_notebook_logger())

get_pretty_notebook_logger() = ConsoleLogger(stdout)
#   A function, as `stdout` is modified, e.g. by IJulia.

end # module


text:
# Basic Corpus type
@auto_hash_equals mutable struct Corpus{S, T<:AbstractDocument{S}}
    documents::Vector{T}
    total_terms::Int
    lexicon::OrderedDict{S, Int}
    inverse_index::OrderedDict{S, Vector{Int}}
    h::TextHashFunction
end

# Constructors
Corpus(docs::Vector{T};
       hash_function::Function = DEFAULT_HASH_FUNCTION,
       cardinality::Int=DEFAULT_CARDINALITY
      ) where T<:AbstractDocument{S} where S<:AbstractString =
    Corpus(
        docs,
        0,
        OrderedDict{S, Int}(),
        OrderedDict{S, Vector{Int}}(),
        TextHashFunction(hash_function, cardinality)
    )

Corpus(docs::Vector{AbstractDocument{S}};
       hash_function::Function = DEFAULT_HASH_FUNCTION,
       cardinality::Int=DEFAULT_CARDINALITY) where S<:AbstractString =
    Corpus(Vector{GenericDocument{S}}(docs),
           hash_function=hash_function,
           cardinality=cardinality)

Corpus(docs::Vector{AbstractDocument};
       hash_function::Function = DEFAULT_HASH_FUNCTION,
       cardinality::Int=DEFAULT_CARDINALITY) where S<:AbstractString = begin
    # Force-convert to String ;)
    n = length(docs)
    T = String
    new_docs = Vector{GenericDocument{T}}(undef, n)
    for i in 1:n
        if docs[i] isa FileDocument
            new_docs[i] = convert(FileDocument{T}, docs[i])
        elseif docs[i] isa StringDocument
            new_docs[i] = convert(StringDocument{T}, docs[i])
        elseif docs[i] isa TokenDocument
            new_docs[i] = convert(TokenDocument{T}, docs[i])
        elseif docs[i] isa NGramDocument
            new_docs[i] = convert(NGramDocument{T}, docs[i])
        else
            # Some other unknown AbstractDocument type
            new_docs[i] = abstract_convert(docs[i], T)
        end
    end
    return Corpus(new_docs; hash_function=hash_function, cardinality=cardinality)
end

# Construct a Corpus from a directory of text files
function DirectoryCorpus(dirname::AbstractString;
                         hash_function::Function = DEFAULT_HASH_FUNCTION,
                         cardinality::Int=DEFAULT_CARDINALITY)
    # Recursive descent of directory
    # Add all non-hidden files to Corpus
    docs = GenericDocument{String}[]

    function add_files(dirname::AbstractString)
        if !isdir(dirname)
            error("DirectoryCorpus() can only be called on directories")
        end
        starting_dir = pwd()
        cd(dirname)
        for filename in readdir(".")
            if isfile(filename) && !occursin(r"^\.", filename)
                push!(docs, FileDocument(abspath(filename)))
            end
            if isdir(filename) && !islink(filename)
                add_files(filename)
            end
        end
        cd(starting_dir)
    end

    add_files(dirname)
    return Corpus(docs, hash_function=hash_function, cardinality=cardinality)
end

# Basic Corpus properties
documents(c::Corpus) = c.documents

# Treat a Corpus as an iterable
function Base.iterate(crps::Corpus, ind=1)
    ind > length(crps.documents) && return nothing
    crps.documents[ind], ind+1
end

Base.eltype(::Type{Corpus{S,T}}) where {S,T} = T

Base.length(crps::Corpus) = length(crps.documents)

Base.size(crps, i) = size(crps.documents, i)

# Treat a Corpus as a container
Base.push!(crps::Corpus, d::AbstractDocument) = push!(crps.documents, d)

Base.pop!(crps::Corpus) = pop!(crps.documents)

Base.pushfirst!(crps::Corpus, d::AbstractDocument) = pushfirst!(crps.documents, d)

Base.popfirst!(crps::Corpus) = popfirst!(crps.documents)

function Base.insert!(crps::Corpus, index::Int, d::AbstractDocument)
    insert!(crps.documents, index, d)
end

Base.deleteat!(crps::Corpus, index::Integer) = deleteat!(crps.documents, index)

# Indexing into a Corpus
#
# (a) Numeric indexing just provides the n-th document
# (b) String indexing is effectively a trivial search engine
Base.getindex(crps::Corpus, ind::Integer) = crps.documents[ind]
Base.getindex(crps::Corpus, inds::Vector{T}) where {T <: Integer} = crps.documents[inds]
Base.getindex(crps::Corpus, r::AbstractRange) = crps.documents[r]
Base.getindex(crps::Corpus, term::AbstractString) = begin
    isempty(crps.inverse_index) && @warn "Inverse index is empty."
    get(crps.inverse_index, term, Int[])
end
# Assignment into a Corpus
function Base.setindex!(crps::Corpus, d::AbstractDocument, ind::Real)
    crps.documents[ind] = d
    return d
end

# Lexicon and inverse index
lexicon(crps::Corpus) = crps.lexicon

function update_lexicon!(crps::Corpus,
                         doc::AbstractDocument,
                         ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY)
    ngs = ngrams(doc, ngram_complexity)
    for (ngram, counts) in ngs
        crps.total_terms += counts
        crps.lexicon[ngram] = get(crps.lexicon, ngram, 0) + counts
    end
end

function update_lexicon!(crps::Corpus,
                         ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY)
    crps.total_terms = 0
    crps.lexicon = OrderedDict{String,Int}()
    for doc in crps
        update_lexicon!(crps, doc, ngram_complexity)
    end
end

function create_lexicon(docs::Union{AbstractVector{S}, Corpus{S,T}},
                        ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY
                       ) where {S,T}
    lexicon = OrderedDict{S, Int}()
    for doc in docs
        ngs = ngrams(doc, ngram_complexity)
        for (ngram, counts) in ngs
            lexicon[ngram] = get(lexicon, ngram, 0) + counts
        end
    end
    return lexicon
end

lexicon_size(crps::Corpus) = length(keys(crps.lexicon))

lexical_frequency(crps::Corpus, term::AbstractString) =
    (get(crps.lexicon, term, 0) / crps.total_terms)

# Work with the Corpus's inverse index
inverse_index(crps::Corpus) = crps.inverse_index

function create_inverse_index(docs::Union{AbstractVector{S}, Corpus{S,T}},
                              ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY
                             ) where {S,T}
    idx = OrderedDict{S, Vector{Int}}()
    @inbounds for i in 1:length(docs)
        ngram_arr = collect(S, keys(ngrams(docs[i], ngram_complexity)))
        for ngram in ngram_arr
            if haskey(idx, ngram)
                push!(idx[ngram], i)
            else
                idx[ngram] = [i]
            end
        end
    end
    for key in keys(idx)
        idx[key] = unique(idx[key])
    end
    return idx
end

function update_inverse_index!(crps::Corpus,
                               ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY)
    crps.inverse_index = create_inverse_index(crps, ngram_complexity)
    return nothing
end

index_size(crps::Corpus) = length(crps.inverse_index)

# Every Corpus prespecifies a hash function for hash trick analysis
hash_function(crps::Corpus) = crps.h

hash_function!(crps::Corpus, f::TextHashFunction) = (crps.h = f; nothing)

# Standardize the documents in a Corpus to a common type
function standardize!(crps::Corpus, ::Type{T}) where T<:AbstractDocument
    for i in 1:length(crps)
        crps.documents[i] = convert(T, crps.documents[i])
    end
end


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

using Documenter
using SpeechFeatures

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

makedocs(
    sitename="SpeechFeatures",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Home" => "index.md",
        "Extracting Features" => "feaextract.md"
    ]
)

deploydocs(
    repo = "github.com/lucasondel/SpeechFeatures.jl.git",
)

max_stars_repo_name: lucasondel/SpeechFeatures

max_issues_repo_name: BUTSpeechFIT/SpeechFeatures

max_forks_repo_name: BUTSpeechFIT/SpeechFeatures


text:
export SSDFun, HuberFun, MVFun, MVFunTotal

"""
	mis,dmis,d2mis = SSDFun(dc,dobs,Wd)

	Input:

		dc::Array   -  simulated data
		dobs::Array -  measured data
		Wd::Array   -  diagonal weighting

	Output:

		mis::Real   -  misfit, 0.5*|dc-dobs|_Wd^2
		dmis        -  gradient
		d2mis       -  diagonal of Hessian

"""
function SSDFun(dc::Union{Array{Float64},Array{Float32}},dobs::Union{Array{Float64},Array{Float32}},Wd::Union{Array{Float64},Array{Float32}})
	res   = vec(dc)-vec(dobs) # predicted - observed data
	Wd    = vec(Wd)
	mis   = .5*real(dot(Wd.*res,Wd.*res))  # data misfit
	dmis  = Wd.*(Wd.*res)
	d2mis = Wd.*Wd
	return mis, dmis, d2mis
end # function SSDFun

"""
	For complex data misfit is computed as 0.5*|real(dc)-(dobs)|_Wd^2 +  0.5*|complex(dc)-complex(dobs)|_W^2
"""
function SSDFun(dc::Array{ComplexF64},dobs::Array{ComplexF64},Wd::Array{ComplexF64})

	wdr   = vec(real(Wd)); wdi = vec(imag(Wd))
	# wdr.*dRe + im*wdi.*dIm
	res   = vec(dc)-vec(dobs)
	resw  = wdr.*real(res) + im*wdi.*imag(res)

	mis   = .5*real(dot(resw,resw))
	dmis  = sdiag(wdr.*wdr)*real(res) + 1im*sdiag(wdi.*wdi)*imag(res)
	d2mis =  wdr.*wdr + im*wdi.*wdi
	return mis,dmis,d2mis
end

"""
	mis,dmis,d2mis = HuberFun(dc,dobs,Wd,C)

	Computes misfit via

		misfit(dc,dobs) = sqrt(abs(Wd*res).^2 + eps)

	Input:
		dc::Array   -  simulated data
		dobs::Array -  measured data
		Wd::Array   -  diagional weighting
		eps         -  conditioning parameter (default=1e-3)

	Output:
		mis::Real   -  misfit
		dmis        -  gradient
		d2mis       -  diagonal of Hessian

"""
function HuberFun(dc::Array{Float64},dobs::Array{Float64},Wd::Array{Float64},eps=1e-3)
	# compute Huber distance
	res   = vec(dc-dobs)
	G     = sqrt.( abs.(Wd.*res).^2 .+ eps)
	mis   = sum(G)
	dmis  = sdiag(Wd./G)*(Wd.*res)
	d2mis = (Wd.*Wd)./G
	return mis,dmis,d2mis
end
# Old Huber fun
# function HuberFun(dc::Array{Float64},dobs::Array{Float64},Wd::Array{Float64},eps=1e-3)
# 	# compute Huber distance
# 	res   = dc-dobs
# 	G     = sqrt.( abs.(Wd.*res).^2 .+ eps)
# 	mis   = sum(G)
# 	dmis  = sdiag(Wd./G)*(Wd.*res)
# 	d2mis = (Wd.^Wd)./G
# 	return mis,dmis,d2mis
# end


text:
using DataDeps
using Test
using ExpectationStubs

# HACK: todo, work out how ExpectationStubs should be changed to make this make sense
Base.open(stub::Stub, t::Any, ::AbstractString) = stub(t)
path_current = @__DIR__

withenv("DATADEPS_ALWAYS_ACCEPT"=>"true") do
    @testset "automatic download" begin
        @stub dummyhash
        @expect(dummyhash(::Any) = [0x12, 0x34])

        @stub dummydown
        @expect(dummydown("http://www.example.com/eg.zip", ::String) = joinpath(path_current, "eg.zip"))
        # HACK: this gives a directory, because we can't trivially mock out the cd

        register(DataDep( "Test1",
         "A dummy message",
         "http://www.example.com/eg.zip",
         (dummyhash, "1234"),#
         fetch_method=dummydown
        ))

        @test endswith(datadep"Test1", "Test1") || endswith(datadep"Test1", "Test1/") ||  endswith(datadep"Test1", "Test1\\")

        @test all_expectations_used(dummyhash)
        @test all_expectations_used(dummydown)

        rm(datadep"Test1"; recursive=true) # delete the directory
    end

    @testset "sanity check the macro's behaviour with variables" begin
        var = "foo/bar"
        macroexpand(@__MODULE__, :(@datadep_str var)) # this line would throw an error if the varibles were being handle wrong
        @test true
    end

    @testset "Ensure when errors occur the datadep will still retrydownloading" begin
        @testset "error in checksum" begin
            dummy_file_path = @__FILE__
            @stub dummydown
            @expect dummydown(::Any, ::Any) = dummy_file_path # give path to an actual file so `open` works

            register(DataDep("TestErrorChecksum", "dummy message", "http://example.void",
                             (error, "1234"); # this will throw an error
                             fetch_method=dummydown))
            @test_throws ErrorException datadep"TestErrorChecksum"
            @test @usecount(dummydown(::Any, ::Any)) == 1

            @test_throws ErrorException datadep"TestErrorChecksum"
            @test @usecount(dummydown(::Any, ::Any)) == 2 # it should have tried to download again
        end

        @testset "error in post fetch" begin
            @stub dummydown
            @expect dummydown(::Any, ::Any) = joinpath(path_current, "eg.zip")

            register(DataDep("TestErrorPostFetch", "dummy message", "http://example.void", Any,
                             fetch_method=dummydown,
                             post_fetch_method = error
                            ))
            @test_throws ErrorException datadep"TestErrorPostFetch"
            @test @usecount(dummydown(::Any, ::Any)) == 1

            @test_throws ErrorException datadep"TestErrorPostFetch"
            @test @usecount(dummydown(::Any, ::Any)) == 2 # it should have tried to download again
        end

        @testset "error in fetch" begin
            use_count = 0
            function error_down(rp,lp)
                use_count += 1
                error("no download for you")
            end

            register(DataDep("TestErrorFetch", "dummy message", "http://example.void", Any,
                             fetch_method = error_down
                            ))
            @test_throws ErrorException datadep"TestErrorFetch"
            @test use_count == 1

            @test_throws ErrorException datadep"TestErrorFetch"
            @test use_count == 2 # it should have tried to download again
        end
    end
end

@testset "Ensure disabled if in CI #70" begin
    withenv("CI"=> "true") do
        @assert !haskey(ENV, "DATADEPS_ALWAYS_ACCEPT")
        @stub dummydown
        @expect dummydown(::Any, ::Any) = @__FILE__ # give path to an actual file so `open` works

        register(DataDep("TestErrorInCI", "dummy message", "http://example.void", Any,
                         fetch_method=dummydown))
        @test_throws DataDeps.DisabledError datadep"TestErrorInCI"
        @test @usecount(dummydown(::Any, ::Any)) == 0

    end
end


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

name = "MPFR"
version = v"4.1.0"

# Collection of sources required to build MPFR
sources = [
    ArchiveSource("https://www.mpfr.org/mpfr-$(version)/mpfr-$(version).tar.xz",
                  "0c98a3f1732ff6ca4ea690552079da9c597872d30e96ec28414ee23c95558a7f"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mpfr-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-shared --disable-static --with-gmp=${prefix} --enable-thread-safe --enable-shared-cache --disable-float128 --disable-decimal-float
make -j${nproc}
make install

# On Windows, make sure non-versioned filename exists...
if [[ ${target} == *mingw* ]]; then
    cp -v ${prefix}/bin/libmpfr-*.dll ${prefix}/bin/libmpfr.dll
fi
"""

# We enable experimental platforms as this is a core Julia dependency
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmpfr", :libmpfr),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.2.0"),
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
version = v"4.1.1"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5", julia_compat="1.6")