text:
const PRECHELT_REF = "[Prechelt, Lutz (1998): \"Early Stopping"*
    "- But When?\", in *Neural Networks: Tricks of the Trade*, "*
    "ed. G. Orr, Springer.](https://link.springer.com/chapter"*
    "/10.1007%2F3-540-49430-8_3)"

const STOPPING_DOC = "An early stopping criterion for loss-reporting "*
    "iterative algorithms. "

const CUSTOM_ALTERNATIVE_DOC = "For a customizable loss-based stopping "*
    "criterion, use [`WithLossDo`](@ref) or [`WithTrainingLossesDo`](@ref) "*
    "with the `stop_if_true=true` option. "

## NEVER

"""
    Never()

$STOPPING_DOC

Indicates early stopping is to be disabled.

See also [`NotANumber`](@ref), for stopping on encountering `NaN`.

"""
struct Never <: StoppingCriterion end

## OUT OF BOUNDS

"""
    InvalidValue()

$STOPPING_DOC

Stop if a loss (or training loss) is `NaN`, `Inf` or `-Inf` (or, more
precisely, if `isnan(loss)` or `isinf(loss)` is `true`).

$CUSTOM_ALTERNATIVE_DOC

"""
struct InvalidValue <: StoppingCriterion end

# state = `true` when `NaN`, `Inf` or `-Inf` has been encountered

update(::InvalidValue, loss, state=false) =
    state || isinf(loss) || isnan(loss)
update_training(::InvalidValue, loss, state=false) =
    state || isinf(loss) || isnan(loss)

done(::InvalidValue, state) = state

message(::InvalidValue, state) = "Stopping early as `NaN`, "*
    "`Inf` or `-Inf` encountered. "

needs_loss(::Type{<:InvalidValue}) = true

## TIME LIMIT

"""
    TimeLimit(; t=0.5)

$STOPPING_DOC

Stopping is triggered after `t` hours have elapsed since the stopping
criterion was initiated.

Any Julia built-in `Real` type can be used for `t`. Subtypes of
`Period` may also be used, as in `TimeLimit(t=Minute(30))`.

Internally, `t` is rounded to nearest millisecond.
``
"""
struct TimeLimit <: StoppingCriterion
    t::Millisecond
    function TimeLimit(t::Millisecond)
        t > Millisecond(0) ||
            throw(ArgumentError("Time limit `t` must be positive. "))
        return new(t)
    end
end
TimeLimit(t::T) where T <: Period = TimeLimit(convert(Millisecond, t))
# for t::T a "numeric" time in hours; assumes `round(Int, ::T)` implemented:
TimeLimit(t) = TimeLimit(round(Int, 3_600_000*t) |> Millisecond)
TimeLimit(; t =Minute(30)) = TimeLimit(t)

# state = time at initialization

update(::TimeLimit, loss) = now()
update_training(::TimeLimit, loss) = now()
update(::TimeLimit, loss, state) = state
done(criterion::TimeLimit, state) = begin
    criterion.t < now() - state
end

## GL

# helper:

generalization_loss(E, E_opt) =  100*(E/abs(E_opt) - one(E_opt))

"""
    GL(; alpha=2.0)

$STOPPING_DOC

A stop is triggered when the (rescaled) generalization loss exceeds
the threshold `alpha`.

**Terminology.** Suppose ``E_1, E_2, ..., E_t`` are a sequence of
losses, for example, out-of-sample estimates of the loss associated
with some iterative machine learning algorithm. Then the
*generalization loss* at time `t`, is given by

`` GL_t = 100 (E_t - E_{opt}) \\over |E_{opt}|``

where ``E_{opt}`` is the minimum value of the sequence.

Reference: $PRECHELT_REF.

"""
struct GL <: StoppingCriterion
    alpha::Float64
    function GL(alpha)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        return new(alpha)
    end
end
GL(; alpha=2.0) = GL(alpha)

update(::GL, loss) = (loss=loss, min_loss=loss)
update(::GL, loss, state) = (loss=loss, min_loss=min(loss, state.min_loss))
# in case first loss consumed was a training loss:
update(criterion::GL, loss, ::Nothing) = update(criterion, loss)
done(criterion::GL, state) =
    generalization_loss(state.loss, state.min_loss) > criterion.alpha

needs_loss(::Type{<:GL}) = true

## PQ

# helpers:

function prepend(v, x, max_length)
    if length(v) < max_length
        vcat([x, ], v)
    else
        vcat([x, ], v[1:end-1])
    end
end

progress(losses::AbstractVector{T}) where T =
    1000*(mean(losses)/minimum(losses) - one(T))

_min(x, ::Nothing) = x
_min(::Nothing, x) = x
_min(x, y) = min(x, y)

"""
    PQ(; alpha=0.75, k=5, tol=eps(Float64))

A stopping criterion for training iterative supervised learners.

A stop is triggered when Prechelt's progress-modified generalization
loss exceeds the threshold `alpha`, or if the training progress drops
below `tol`. Here `k` is the maximum number of training (in-sample)
losses to be used to estimate the training progress.

**Context and explanation of terminology.** The progress-modified loss
is defined in the following scenario: Estimates, ``E_1, E_2, ...,
E_t``, of the out-of-sample loss of an iterative supervised learner
are being computed, but not necessarily at every iteration. However,
training losses for every iteration *are* being made available
(usually as a by-product of training) which can be used to quantify
recent training progress, as follows.

Fix a time ``j``, corresponding to some out-of-sample loss ``E_j``,
and let ``F_1`` be the corresponding training loss, ``F_2`` the
training loss in the previous interation of the model, ``F_3``, the
training loss two iterations previously, and so on. Let ``K`` denote
the number of model iterations since the last out-of-sample loss
``E_{j-1}`` was computed, or `k`, whichever is the smaller.  Then
the *training progress* at time ``j`` is defined by

`` P_j = 1000 |M - m|/|m| ``

where ``M`` is the mean of the training losses ``F_1, F_2, \\ldots ,
F_K`` and ``m`` the minimum value of those losses.

The *progress-modified generalization loss* at time ``t`` is given by

`` PQ_t = GL_t / P_t``

where ``GL_t`` is the generalization loss at time ``t``; see
[`GL`](@ref).

Reference: $PRECHELT_REF.

"""
struct PQ <: StoppingCriterion
    alpha::Float64
    k::Union{Int,Float64} # could be Inf
    tol::Float64
    function PQ(alpha, k, tol)
        alpha > 0 ||
            throw(ArgumentError("Threshold `alpha` must be positive. "))
        if k < 2 || !(k isa Int || isinf(k))
            throw(ArgumentError("`k` must `Int` or `Inf` and `k > 1`. "))
        end
        tol > 0 ||
            throw(ArgumentError("`tol` must be positive. "))

        return new(alpha, k, tol)
    end
end
PQ(; alpha=0.75, k=5, tol=eps(Float64)) = PQ(alpha, k, tol)

struct PQState{T}
    training_losses::Vector{T}
    waiting_for_out_of_sample::Bool
    loss::Union{Nothing,T}
    min_loss::Union{Nothing,T}
end

function update_training(criterion::PQ, loss)
    training_losses = [loss, ]
    return PQState(training_losses,
                   true,
                   nothing,
                   nothing)
end

update(::PQ, loss) = error("First loss reported to the PQ early stopping "*
                           "algorithm must be a training loss. ")

function update_training(criterion::PQ, loss, state)
    training_losses = if state.waiting_for_out_of_sample
        prepend(state.training_losses, loss, criterion.k)
    else
        [loss, ]
    end
    return PQState(training_losses,
                   true,
                   state.loss,
                   state.min_loss)
end

function update(::PQ, loss, state)
    length(state.training_losses) > 1 ||
        error("The PQ stopping criterion requires at least two training "*
              "losses between out-of-sample loss updates. ")
    return PQState(state.training_losses,
                   false,
                   loss,
                   _min(loss, state.min_loss))
end

function done(criterion::PQ, state)
    state.waiting_for_out_of_sample && return false
    GL = generalization_loss(state.loss, state.min_loss)
    P = progress(state.training_losses)
    P > criterion.tol || return true
    PQ = GL/P
    return  PQ > criterion.alpha
end

needs_loss(::Type{<:PQ}) = true
needs_training_losses(::Type{<:PQ}) = true

## PATIENCE

# This is UP_s in Prechelt 1998

"""
    Patience(; n=5)

$STOPPING_DOC

A stop is triggered by `n` consecutive increases in the loss.

Denoted "_UP_s" in $PRECHELT_REF.

$CUSTOM_ALTERNATIVE_DOC

"""
struct Patience <: StoppingCriterion
    n::Int
    function Patience(n::Int)
        n > 0 ||
            throw(ArgumentError("The patience level `n` must be positive. "))
        return new(n)
    end
end
Patience(; n=5) = Patience(n)

# Prechelt alias:
const UP = Patience

update(criterion::Patience, loss) = (loss=loss, n_increases=0)
@inline function update(criterion::Patience, loss, state)
    old_loss, n = state
    if loss > old_loss
        n += 1
    else
        n = 0
    end
    return (loss=loss, n_increases=n)
end
# in case first loss consumed was a training loss:
update(criterion::Patience, loss, ::Nothing) = update(criterion, loss)

done(criterion::Patience, state) = state.n_increases == criterion.n

needs_loss(::Type{<:Patience}) = true

## NUMBER SINCE BEST

"""
    NumberSinceBest(; n=6)

$STOPPING_DOC

A stop is triggered when the number of calls to the control, since the
lowest value of the loss so far, is `n`.

$CUSTOM_ALTERNATIVE_DOC

"""
struct NumberSinceBest <: StoppingCriterion
    n::Int
    function NumberSinceBest(n::Int)
        n > 0 ||
            throw(ArgumentError("`n` must be positive. "))
        return new(n)
    end
end
NumberSinceBest(; n=6) = NumberSinceBest(n)

update(criterion::NumberSinceBest, loss) = (best=loss, number_since_best=0)
@inline function update(criterion::NumberSinceBest, loss, state)
    best, number_since_best = state
    if loss < best
        best = loss
        number_since_best = 0
    else
        number_since_best += 1
    end
    return (best=best, number_since_best=number_since_best)
end

# in case first loss consumed was a training loss:
update(criterion::NumberSinceBest, loss, ::Nothing) = update(criterion, loss)

done(criterion::NumberSinceBest, state) = state.number_since_best == criterion.n

needs_loss(::Type{<:NumberSinceBest}) = true

# # NUMBER LIMIT

"""
    NumberLimit(; n=100)

$STOPPING_DOC

A stop is triggered by `n` consecutive loss updates, excluding
"training" loss updates.

If wrapped in a `stopper::EarlyStopper`, this is the number of calls
to `done!(stopper)`.

"""
struct NumberLimit <: StoppingCriterion
    n::Int
    function NumberLimit(n::Int)
        n > 0 ||
            throw(ArgumentError("`n` must be positive. "))
        return new(n)
    end
end
NumberLimit(; n=100) = NumberLimit(n)

update(criterion::NumberLimit, loss) = 1
@inline function update(criterion::NumberLimit, loss, state)
    return state+1
end
# in case first loss consumed was a training loss:
update(criterion::NumberLimit, loss, ::Nothing) = update(criterion, loss)

done(criterion::NumberLimit, state) = state == criterion.n

# ## THRESHOLD

"""
    Threshold(; value=0.0)

$STOPPING_DOC

A stop is triggered as soon as the loss drops below `value`.

$CUSTOM_ALTERNATIVE_DOC

"""
struct Threshold <: StoppingCriterion
    value::Float64
end
Threshold(; value=0.0) = Threshold(value)

update(criterion::Threshold, loss) = loss
update(criterion::Threshold, loss, state) = loss
# in case first loss consumed was a training loss:
update(criterion::Threshold, loss, ::Nothing) = loss

done(criterion::Threshold, state) = state < criterion.value

needs_loss(::Type{<:Threshold}) = true

"""
    Warmup(c::StoppingCriterion, n)

Wait for `n` updates before checking stopping criterion `c`
"""
struct Warmup{C} <: StoppingCriterion where {C <: StoppingCriterion}
    criterion::C
    n::Int
    function Warmup(criterion::C, n::N) where {C <: StoppingCriterion, N <: Integer}
        n > 0 || throw(ArgumentError("`n` must be positive. "))
        new{C}(criterion, Int(n))
    end
end

# Constructors for Warmup
Warmup() = Warmup(InvalidValue())   # Default for testing
Warmup(c; n = 1) = Warmup(c, n)     # Provide kwargs interface

# Initialize inner state for type-stability, and record first observation
update(c::Warmup, loss, ::Nothing) = update(c, loss)
update(criterion::Warmup, loss) = (1, update(criterion.criterion, loss))

# Catch uninitialized state
update_training(c::Warmup, loss, ::Nothing) = update_training(c, loss)
update_training(c::Warmup, loss) = (1, update_training(c.criterion, loss))

# Handle update vs update_training
update(c::Warmup, loss, state) = _update(update, c, loss, state)
update_training(c::Warmup, loss, state) = _update(update_training, c, loss, state)
needs_loss(::Type{<:Warmup{C}}) where C = needs_loss(C)
needs_training_losses(::Type{<:Warmup{C}}) where C = needs_training_losses(C)

# Dispatch update and update_training here
function _update(f::Function, criterion::Warmup, loss, state)
    n, inner = state
    n += 1
    if n <= criterion.n
        # Skip inner criterion
        return n, inner
    elseif n == criterion.n+1
        # First step of inner criterion
        return n, f(criterion.criterion, loss)
    else
        # Step inner criterion
        return n, f(criterion.criterion, loss, inner)
    end
end

function done(criterion::Warmup, state)
    # Only check if inner criterion is done after n updates
    return state[1] <= criterion.n ? false : done(criterion.criterion, state[2])
end

message(c::Warmup, state) = message(c.criterion, state[2])

## NOT A NUMBER (deprecated)

"""
    NotANumber()

$STOPPING_DOC

Stop if a loss of `NaN` is encountered.

**Now deprecated** in favour of [`InvalidValue`](@ref).

"""
struct NotANumber <: StoppingCriterion
    function NotANumber()
        Base.depwarn("`NotANumber()` is deprecated. Use `InvalidValue()` "*
                     "to trap `NaN`, `Inf` or `-Inf`. ", :NotANumber)
        return new()
    end
end

# state = `true` when NaN has been encountered

update(::NotANumber, loss) = isnan(loss)
update_training(::NotANumber, loss) = isnan(loss)

update(::NotANumber, loss, state) = state || isnan(loss)
update_training(::NotANumber, loss, state) = state || isnan(loss)

done(::NotANumber, state) = state

message(::NotANumber, state) = "Stopping early as NaN encountered. "

needs_loss(::Type{<:NotANumber}) = true


text:
function loadmnist()
    # (c) David Barber, University College London 2015
    alldata=matread("mnist-original.mat")
    images=(float(alldata["data"]))./255
    return images, alldata["label"]
end

ext: jl

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
### function plot_margin(gain, phase, w, Gm, Pm, wcg, wcp)

Plot the bode diagram with the margin information.

##### Args

* gain: Gain [dB].
* phase: Phase [deg].
* w: Frequencies [rad/s].
* Gm: Gain margin.
* Pm: Phase margin [deg/s].
* wcg: Phase cross frequency for the gain margin computation [rad/s].
* wcp: Gain cross frequency for the phase margin computation [rad/s].

##### Parameters

* tol: Tolerance to check if a value is grater than 0.

"""

function plot_margin(gain, phase, w, Gm, Pm, wcg, wcp; tol = sqrt(eps()))
    # Compute the gain margin in dB.
    Gm_dB = 20*log10(Gm)

    if (plot_system == :pyplot)
        fig, ax = plt.subplots(2, 1, sharex=true)

        # Set window title.
        fig[:canvas][:set_window_title]("Bode Diagram")

        # Plot data.
        ax[1][:semilogx](w, gain,  linewidth=plw, zorder=6)
        ax[2][:semilogx](w, phase, linewidth=plw, zorder=6)
        ax[1][:grid]("on", which="major", color="#AAAAAA", linewidth=glw)
        ax[1][:grid](zorder = 0)
        ax[2][:grid]("on", which="major", color="#AAAAAA", linewidth=glw)
        ax[2][:grid](zorder = 0)

        # Change figure border width.
        ax[1][:spines]["top"][:set_linewidth](blw)
        ax[1][:spines]["bottom"][:set_linewidth](blw)
        ax[1][:spines]["left"][:set_linewidth](blw)
        ax[1][:spines]["right"][:set_linewidth](blw)

        ax[2][:spines]["top"][:set_linewidth](blw)
        ax[2][:spines]["bottom"][:set_linewidth](blw)
        ax[2][:spines]["left"][:set_linewidth](blw)
        ax[2][:spines]["right"][:set_linewidth](blw)

        # Title of the figure.
        title = s"$\bf{GM} = $"

        if !isinf(Gm)
            title *= "$(round(Gm_dB,2)) dB at $(round(wcg,2)) rad/s, "
        else
            title *= "Inf, "
        end

        title *= s"$\bf{PM} = $"

        if !isinf(Pm)
            title *= "$(round(Pm,2)) deg at $(round(wcp,2)) rad/s"
        else
            title *= "Inf"
        end

        fig[:suptitle](title, fontsize=tifsz)

        # Axes labels.
        ax[2][:set_xlabel]("Freq. (rad/s)")

        ax[1][:set_ylabel]("Gain (dB)")
        ax[2][:set_ylabel]("Phase (deg)")

        # Plot a line indicating the margins if applicable.
        if (!isinf(Gm)) && (wcg > tol)
            ax[1][:arrow](wcg, 0, 0, -Gm_dB,
                          linewidth=3.0, length_includes_head=true,
                          head_width=0, head_length=0,
                          facecolor="#F0B300", edgecolor="#F0B300", zorder=4)

            ax[1][:axhline](y = 0.0, color="k", zorder=5)
            ax[2][:axvline](x = wcg, color="k", zorder=3, linestyle="--")
            ax[1][:axvline](x = wcg, color="k", zorder=3, linestyle="--")
        end

        if !isinf(Pm) && (wcp > tol)
            # Find the approximate phase value of the transfer function at the
            # crossing frequency.
            ind_wcp   = indmin(abs.(w-wcp))
            phase_wcp = phase[ind_wcp][1]

            # Find the reference that was used to compute the phase margin.
            pc_ref = round(phase_wcp/180)*180
            (pc_ref == 0.0) && (pc_ref = sign(phase_wcp)*180.0)

            # Avoid problems if phase margin is 0.
            if (Pm != 0.0)
                ax[2][:arrow](wcp, pc_ref, 0, +Pm,
                              linewidth=3.0, length_includes_head=true,
                              head_width=0, head_length=0,
                              facecolor="#F0B300", edgecolor="#F0B300", zorder=3)
            end

            ax[2][:axhline](y = pc_ref, color="k", zorder=5)
            ax[2][:axvline](x = wcp,    color="k", zorder=3, linestyle="--")
            ax[1][:axvline](x = wcp,    color="k", zorder=3, linestyle="--")
        end
    end

    nothing
end


text:
using AstroTime.AstroDates
import Dates

const REFERENCES = (
    (-4713, 12, 31, -2451546),
    (-4712, 01, 01, -2451545),
    ( 0000, 12, 31,  -730122),
    ( 0001, 01, 01,  -730121),
    ( 1500, 02, 28,  -182554),
    ( 1500, 02, 29,  -182553),
    ( 1500, 03, 01,  -182552),
    ( 1582, 10, 04,  -152385),
    ( 1582, 10, 15,  -152384),
    ( 1600, 02, 28,  -146039),
    ( 1600, 02, 29,  -146038),
    ( 1600, 03, 01,  -146037),
    ( 1700, 02, 28,  -109514),
    ( 1700, 03, 01,  -109513),
    ( 1800, 02, 28,   -72990),
    ( 1800, 03, 01,   -72989),
    ( 1858, 11, 15,   -51546),
    ( 1858, 11, 16,   -51545),
    ( 1999, 12, 31,       -1),
    ( 2000, 01, 01,        0),
    ( 2000, 02, 28,       58),
    ( 2000, 02, 29,       59),
    ( 2000, 03, 01,       60),
)

@testset "DateTime" begin
    @testset for ref in REFERENCES
        s = AstroTime.Date(ref[end])
        @test year(s) == ref[1]
        @test month(s) == ref[2]
        @test day(s) == ref[3]
        @test j2000(s) == ref[end]
        @test j2000(AstroTime.Date(ref[1:3]...)) == ref[end]
    end

    @test_throws ArgumentError AstroTime.Date(2018, 2, 29)
    @test_throws ArgumentError AstroTime.Date(2018, 0, 1)
    @test_throws ArgumentError AstroTime.Date(2018, 13, 1)
    @test_throws ArgumentError AstroTime.Time(24, 59, 59.0)
    @test_throws ArgumentError AstroTime.Time(23, 60, 59.0)
    @test_throws ArgumentError AstroTime.Time(23, 59, 61.0)
    @test_throws ArgumentError AstroTime.Time(23, 59, 59, 2.0)
    @test_throws ArgumentError AstroTime.Time(23, 59, 59, -0.2)
    @test_throws ArgumentError AstroTime.Time(86401, 0)

    @test AstroTime.Date(2000, 1) == AstroTime.Date(2000, 1, 1)
    @test AstroTime.Date(-2000, 1) == AstroTime.Date(-2000, 1, 1)
    @test AstroTime.Date(1000, 1) == AstroTime.Date(1000, 1, 1)

    dt = AstroTime.DateTime(2020, 3, 21, 10, 15, 37.245)
    dt1 = AstroTime.DateTime(2020, 3, 21, 10, 15, 37, 0.245)
    d = AstroTime.Date(dt)
    t = AstroTime.Time(dt)

    @test dt ≈ dt1

    @test AstroTime.DateTime(Dates.DateTime(dt)) == dt
    @test AstroTime.Date(Dates.Date(d)) == d
    @test AstroTime.Time(Dates.Time(t)) ≈ t

    @test year(dt) == 2020
    @test year(d) == 2020
    @test month(dt) == 3
    @test month(d) == 3
    @test day(dt) == 21
    @test day(d) == 21
    @test hour(dt) == 10
    @test hour(t) == 10
    @test minute(dt) == 15
    @test minute(t) == 15
    @test second(dt) == 37
    @test second(t) == 37
    @test second(Float64, dt) == 37.245
    @test second(Float64, t) == 37.245
    @test millisecond(dt) == 245
    @test millisecond(t) == 245

    @test secondinday(t) == 36937.245

    jd = Dates.datetime2julian(Dates.DateTime(dt))
    @test julian(dt) ≈ jd
    @test sum(julian_twopart(dt)) ≈ jd
    @test j2000(dt) ≈ jd - AstroTime.AstroDates.J2000

    pt = AstroTime.Time(12, 12, 12, 0.123456789)
    @test millisecond(pt) == 123
    @test microsecond(pt) == 456
    @test nanosecond(pt) == 789
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
using Weave
using Test

# Test rendering of doc chunks
content = """
# Test chunk

Test rendering \$\alpha\$
"""

dchunk = Weave.DocChunk(content, 1, 1)

pformat = Weave.formats["github"]
f = Weave.format_chunk(dchunk, pformat.formatdict, pformat)
@test f == content

docformat = Weave.formats["md2html"]
f_check = "
Test chunk
\n
Test rendering \$\alpha\$
\n"
f = Weave.format_chunk(dchunk, docformat.formatdict, docformat)
@test f_check == f

# Test with actual doc

parsed = Weave.read_doc("documents/chunk_options.noweb")
doc = Weave.run(parsed, doctype = "md2html")

c_check = "
\nx = [12, 10]\nprintln(y)\n
\n"
doc.format.formatdict[:theme] = doc.highlight_theme
c = Weave.format_code(doc.chunks[3].content, doc.format)
@test c_check == c

o_check = "\nprintln(x)\n"
o = Weave.format_output(doc.chunks[4].content, doc.format)
@test o_check == o

doc.template = "templates/mini.tpl"
rendered = Weave.render_doc("Hello", doc, doc.format)
@test rendered == "\nHello\n"

# Tex format
parsed = Weave.read_doc("documents/chunk_options.noweb")
doc = Weave.run(parsed, doctype = "md2tex")

c_check = "\\begin{lstlisting}\n(*@\\HLJLnf{println}@*)(*@\\HLJLp{(}@*)(*@\\HLJLn{x}@*)(*@\\HLJLp{)}@*)\n\\end{lstlisting}\n"
doc.format.formatdict[:theme] = doc.highlight_theme
c = Weave.format_code(doc.chunks[4].content, doc.format)
@test c_check == c

o_check = "\nx = [12, 10]\nprintln(y)\n"
o = Weave.format_output(doc.chunks[3].content, doc.format)
@test o_check == o

doc.template = "templates/mini.tpl"
rendered = Weave.render_doc("Hello", doc, doc.format)
@test rendered == "\nHello\n"

# Test header parsing and stripping
header = """
---
title : Test block
author : Matti Pastell
---

# Actual header

and some text

"""

dchunk = Weave.DocChunk(header, 1, 1)
h = Weave.parse_header(dchunk)
h_ref = Dict("author" => "Matti Pastell", "title" => "Test block")
@test h_ref == h

htext = Weave.strip_header(dchunk)
h_ref = """
# Actual header

and some text

"""
@test htext.content[1].content == h_ref

# Test wrapping

cows = repeat("🐄", 100)
testcows = """
🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄
🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄🐄"""

wcows = Weave.wrapline(cows)

@test wcows == testcows
@test length(split(wcows, "\n")[1]) == 75
@test length(split(wcows, "\n")[2]) == 25

tfied = "\\ensuremath{\\bm{\\mathrm{L}}} \\ensuremath{\\bm{\\mathfrak{F}}} \\ensuremath{\\bm{\\iota}} \\ensuremath{\\mathfrak{A}} \\ensuremath{\\bm{\\varTheta}}"

@test Weave.uc2tex("𝐋 𝕱 𝛊 𝔄 𝚹") == tfied

# Test markdown output from chunks
parsed = Weave.read_doc("documents/markdown_output.jmd")
doc = Weave.run(parsed, doctype = "md2html")
@test doc.chunks[1].rich_output == "\n

Small markdown sample
\n
Hello from code block.
\n
"
@test doc.chunks[2].rich_output == "\n

\n

one
\n
\n

two
\n
\n

three
\n
\n
\n
"

ldoc = Weave.run(parsed, doctype = "md2tex")
@test ldoc.chunks[1].rich_output == "\n\\subsubsection{Small markdown sample}\n\\textbf{Hello} from \\texttt{code} block.\n\n"
@test ldoc.chunks[2].rich_output == "\n\\begin{itemize}\n\\item one\n\n\n\\item two\n\n\n\\item three\n\n\\end{itemize}\n"

mdoc = Weave.run(parsed, doctype = "github")
@test mdoc.chunks[1].rich_output == "\n\n### Small markdown sample\n\n**Hello** from `code` block.\n\n"
@test mdoc.chunks[2].rich_output == "\n\n* one\n* two\n* three\n\n"


text:
function equality_test(dist, d, A, B)
    for i in CartesianIndices(A)
        for j in CartesianIndices(B)
            dist[i, j] != d(A[i], B[j]) && return false
        end
    end
    return true
end
function equality_test(dist, index_map, d, A, B)
    for i in CartesianIndices(A)
        for j in CartesianIndices(B)
            p = A[index_map(i)]
            q = B[index_map(j)]
            dist[i, j] != d(p, q) && return false
        end
    end
    return true
end
function equality_test(dist, index_map, getindex_op, d, A, B)
    for i in CartesianIndices(A)
        for j in CartesianIndices(B)
            p = getindex_op(A, index_map(i))
            q = getindex_op(B, index_map(j))
            dist[i, j] != d(p, q) && return false
        end
    end
    return true
end

@testset "pairwise" begin
    @testset "NullCache" begin
        A, B = rand(1:5, 6), rand(1:5, 4)
        d = Euclidean()

        # (_, _, f)
        dist = PairwiseDistance(d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, d, A, B)

        # {T}(_, _, f)
        dist = PairwiseDistance{Float32}(d, (A, B))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, d, A, B)

        # (index_map, _, f)
        index_map = i->CartesianIndex(1)
        dist = PairwiseDistance(index_map, d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, d, A, B)

        # (index_map_list, _, f)
        index_map = (
            i->CartesianIndex(1),
            j->CartesianIndex(2)
        )
        dist = PairwiseDistance(index_map, d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test all(isequal(d(A[1], B[2])), dist)

        # {T}(index_map, _, f)
        index_map = i->CartesianIndex(1)
        dist = PairwiseDistance{Float32}(index_map, d, (A, B))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, d, A, B)

        # {T}(index_map_list, _, f)
        index_map = (
            i->CartesianIndex(1),
            j->CartesianIndex(2)
        )
        dist = PairwiseDistance{Float32}(index_map, d, (A, B))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test all(isequal(d(A[1], B[2])), dist)

        # (index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map_list, getindex_op_list, f)
        index_map = (
            i->i[1],
            j->j[1]
        )
        getindex_op = (
            (A, i)->i,
            (A, j)->j
        )
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map[1], getindex_op[1], d, A, B)

        # {T}(index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map_list, getindex_op_list, f)
        index_map = (
            i->i[1],
            j->j[1]
        )
        getindex_op = (
            (A, i)->i,
            (A, j)->j
        )
        dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map[1], getindex_op[1], d, A, B)

        @testset "Invalid cases" begin
            # all inputs should be tuple

            index_map = [
                i->i[1],
                j->j[1]
            ]
            getindex_op = [
                (A, i)->i,
                (A, j)->j
            ]
            dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B));
            @test_throws MethodError dist[1, 1]

            @test_throws MethodError PairwiseDistance(d, A, B);
            @test_throws MethodError PairwiseDistance{Float32}(d, A, B);
            @test_throws MethodError PairwiseDistance(index_map, getindex_op, d, A, B);
            @test_throws MethodError PairwiseDistance{Float32}(index_map, getindex_op, d, A, B);

            # Currently we only support two arrays
            @test_throws ArgumentError PairwiseDistance{Float32}(d, (A, B, A));
        end
    end

    @testset "LocalWindowCache" begin
        A, B = rand(1:5, 6), rand(1:5, 4)
        d = Euclidean()

        # (_, _, f)
        dist = PairwiseDistance(d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, d, A, B)

        # {T}(_, _, f)
        dist = PairwiseDistance{Float32}(d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, d, A, B)

        # (index_map, _, f)
        index_map = i->CartesianIndex(1)
        dist = PairwiseDistance(index_map, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, d, A, B)

        # (index_map_list, _, f)
        index_map = (
            i->CartesianIndex(1),
            j->CartesianIndex(2)
        )
        dist = PairwiseDistance(index_map, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test all(isequal(d(A[1], B[2])), dist)

        # {T}(index_map, _, f)
        index_map = i->CartesianIndex(1)
        dist = PairwiseDistance{Float32}(index_map, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, d, A, B)

        # {T}(index_map_list, _, f)
        index_map = (
            i->CartesianIndex(1),
            j->CartesianIndex(2)
        )
        dist = PairwiseDistance{Float32}(index_map, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test all(isequal(d(A[1], B[2])), dist)

        # (index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map_list, getindex_op_list, f)
        index_map = (
            i->i[1],
            j->j[1]
        )
        getindex_op = (
            (A, i)->i,
            (A, j)->j
        )
        dist = PairwiseDistance(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == result_type(d, A, B)
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map[1], getindex_op[1], d, A, B)

        # {T}(index_map, getindex_op, f)
        index_map = i->i[1]
        getindex_op = (A, i)->i
        dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map, getindex_op, d, A, B)

        # (index_map_list, getindex_op_list, f)
        index_map = (
            i->i[1],
            j->j[1]
        )
        getindex_op = (
            (A, i)->i,
            (A, j)->j
        )
        dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )))
        @test eltype(dist) == Float32
        @test axes(dist) == (axes(A, 1), axes(B, 1))
        @test equality_test(dist, index_map[1], getindex_op[1], d, A, B)

        @testset "Invalid cases" begin
            # all inputs should be tuple

            index_map = [
                i->i[1],
                j->j[1]
            ]
            getindex_op = [
                (A, i)->i,
                (A, j)->j
            ]
            dist = PairwiseDistance{Float32}(index_map, getindex_op, d, (A, B), LocalWindowCache((4, )));
            @test_throws MethodError dist[1, 1]

            @test_throws MethodError PairwiseDistance(d, A, B, LocalWindowCache((4, )));
            @test_throws MethodError PairwiseDistance{Float32}(d, A, B, LocalWindowCache((4, )));
            @test_throws MethodError PairwiseDistance(index_map, getindex_op, d, A, B, LocalWindowCache((4, )));
            @test_throws MethodError PairwiseDistance{Float32}(index_map, getindex_op, d, A, B, LocalWindowCache((4, )));

            # Currently we only support two arrays
            @test_throws ArgumentError PairwiseDistance{Float32}(d, (A, B, A), LocalWindowCache((4, )));
        end
    end
end


text:
using Documenter, IOUtils

makedocs(;
    modules=[IOUtils],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jagot/IOUtils.jl/blob/{commit}{path}#L{line}",
    sitename="IOUtils.jl",
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>",
    doctest=false
)

deploydocs(;
    repo="github.com/jagot/IOUtils.jl",
)


text:
function test_mat2mat_base(f::Function, M, N)
    srand(98765)
    srcmat = repmat(rand(M), 1, N)
    dstlen = length(f(srcmat[:,1]))
    dstmat = Array{Float64}(dstlen, N)
    for i = 1:N
        dstmat[:,i] = f(srcmat[:,i])
    end
    dstmat2 = f(srcmat)
    @test dstmat2 ≈ dstmat
end

function test_mat2mat()
    src = zeros(20)
    srcmat = zeros(20, 2)

    # all extended functions should accept function calls like:
    # 1. f(x) where x is a input vector
    # 2. f(x) where x is a input matrix
    for f in SPTK.vec2vec
        @show f
        @test applicable(eval(f), src)
        @test applicable(eval(f), srcmat)
    end

    # test_mat2mat_base(mcep, 1024, 2)
    # test_mat2mat_base(mfcc, 512, 2)
    test_mat2mat_base(lpc2c, 20, 2)
    test_mat2mat_base(mc2b, 20, 2)
    test_mat2mat_base(gc2gc, 20, 2)
    test_mat2mat_base(mgc2mgc, 20, 2)
end

println("-- test mat2mat conversion applicable (col-wise)")
test_mat2mat()


text:
"""
```
Module      : FinanceLib.Derivatives.Forwards
Description : Implement Derivatives modules for the FinanceLib library
Copyright   : (c) 2022 Kishaloy Neogi
License     : MIT
Maintainer  : Kishaloy Neogi
Email       : nkishaloy@yahoo.com
```
The module describes the base modules of Forward Derivatives.

You may see the github repository at <https://github.com/n-kishaloy/FinanceLib.jl>
"""
module Forwards

using Dates: Date
using FinanceLib: yearFrac

"""
A struct defining Forward contracts for Assets, with periods given as Double
```
struct Forward
  r     = risk-free rate of return per period
  T     = Forward expiry date in Double
  FT    = Forward rate @ T = S0*(1+r)^T
  benf  = Dividends and other benefits - Cost of holding the asset
end
```
"""
struct Forward
  r     :: Float64 # risk-free rate of return per period
  T     :: Float64 # Forward expiry date in Double
  FT    :: Float64 # Forward rate @ T = S0*(1+r)^T
  benf  :: Float64 # Dividends and other benefits - Cost of holding the asset
end

"""
A struct defining Forward contracts for Assets, with periods given as Date
  ```
struct XForward
  r     # risk-free rate of return per year
  T0    # Contract start date as Date
  T     # Forward expiry date as Date
  FT    # Forward rate @ T = S0*(1+r)^T
  benf  # Dividends and other benefits - Cost of holding the asset
end
```
"""
struct XForward
  r     ::  Float64     # risk-free rate of return per year
  T0    ::  Date        # Contract start date as Date
  T     ::  Date        # Forward expiry date as Date
  FT    ::  Float64     # Forward rate @ T = S0*(1+r)^T
  benf  ::  Float64     # Dividends and other benefits - Cost of holding the asset
end

Forward_S0(;r, T, S0, benf=0.0) = Forward(r, T, (S0 - benf)*(1+r)^T, benf)

XForward_S0(;r, T0, T, S0, benf=0.0) = XForward(r, T0, T, (S0 - benf)*(1+r)^(yearFrac(T0,T)), benf)

getS0(x::Forward) = x.FT/(1 + x.r)^x.T + x.benf

getS0(x::XForward) = x.FT/(1 + x.r)^(yearFrac(x.T0,x.T)) + x.benf

valuation(x::Forward, St, t) = St - x.benf*(1+x.r)^t - x.FT*(1 + x.r)^(t - x.T)

valuation(x::XForward, St, t) = St - x.benf*(1+x.r)^yearFrac(x.T0,t) - x.FT*(1 + x.r)^yearFrac(x.T, t)

end


text:
# info("Running Test: testmonitor.jl")
# include("testmonitor.jl")
info(" XXX modified XXX")
info("Running Test: testTAP.jl")
include("testTAP.jl")


text:
#=
P06 (*) Find out whether a list is a palindrome.
A palindrome can be read forward or backward; e.g. (x a m a x).
=#

function is_palindrome(coll)
    # With built-in `reverse` function
    coll == reverse(coll)
end

function is_palindrome_2(coll)
    # With tail recursion
    if isempty(coll)
        true
    else
        coll[1] == coll[end] && is_palindrome_2(coll[2:end - 1])
    end
end

@assert is_palindrome((:x, :a, :m, :a, :x))
@assert !is_palindrome((:a, :b, :c, :d, :b, :a))

@assert is_palindrome_2((:x, :a, :m, :a, :x))
@assert !is_palindrome_2((:a, :b, :c, :d, :b, :a))

println("Tests passed: JL-06.jl")


text:
function real_gdp(interval::String = "annual"; client=GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(interval, ["quarterly", "annual"])
    @argcheck in(datatype, ["json", "csv", nothing])
    params = Dict(
        "function" => "REAL_GDP",
        "interval" => interval,
        "datatype" => isnothing(datatype) ? "csv" : datatype,
        "apikey" => key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

function treasury_yield(interval::String = "monthly", maturity::String = "10year"; client=GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(interval, ["daily", "weekly", "monthly", nothing])
    @argcheck in(maturity, ["3month", "5year", "10year", "30year", nothing])
    @argcheck in(datatype, ["json", "csv", nothing])
    params = Dict(
        "function" => "TREASURY_YIELD",
        "interval" => interval,
        "maturity" => maturity,
        "datatype" => isnothing(datatype) ? "csv" : datatype,
        "apikey" => key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

function federal_fund_rate(interval::String = "monthly"; client=GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(interval, ["daily", "weekly", "monthly", nothing])
    @argcheck in(datatype, ["json", "csv", nothing])
    params = Dict(
        "function" => "FEDERAL_FUNDS_RATE",
        "interval" => interval,
        "datatype" => isnothing(datatype) ? "csv" : datatype,
        "apikey" => key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

function cpi(interval::String = "monthly"; client=GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(interval, ["monthly", "semiannual", nothing])
    @argcheck in(datatype, ["json", "csv", nothing])

    params = Dict(
        "function" => "CPI",
        "datatype" => isnothing(datatype) ? "csv" : datatype,
        "apikey" => key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

for func in (:real_gdp_per_capita, :inflation,
             :inflation_expectation, :consumer_sentiment,
             :retail_sales, :durables,
             :unemployment, :nonfarm_payroll)

    @eval begin
        function $(func)(; client = GLOBAL[], datatype::Union{String, Nothing}=nothing, parser = "default")
            @argcheck in(datatype, ["json", "csv", nothing])
            params = Dict(
                "function"=>uppercase($(String(func))),
                "datatype"=>isnothing(datatype) ? "csv" : datatype,
                "apikey"=>key(client)
            )
            uri = _build_uri(client.scheme, client.host, "query", params)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
        export $func
    end

end


text:
using Boltz, Lux

models_available = Dict(alexnet => [(:alexnet, true), (:alexnet, false)],
                        convmixer => [(:small, false), (:base, false), (:large, false)],
                        densenet => [
                            (:densenet121, false),
                            (:densenet161, false),
                            (:densenet169, false),
                            (:densenet201, false),
                        ],
                        googlenet => [(:googlenet, false)],
                        mobilenet => [
                            (:mobilenet_v1, false),
                            (:mobilenet_v2, false),
                            (:mobilenet_v3_small, false),
                            (:mobilenet_v3_large, false),
                        ],
                        resnet => [
                            (:resnet18, true),
                            (:resnet18, false),
                            (:resnet34, true),
                            (:resnet34, false),
                            (:resnet50, true),
                            (:resnet50, false),
                            (:resnet101, true),
                            (:resnet101, false),
                            (:resnet152, true),
                            (:resnet152, false),
                        ],
                        resnext => [
                            (:resnext50, false),
                            (:resnext101, false),
                            (:resnext152, false),
                        ],
                        vgg => [
                            (:vgg11, false),
                            (:vgg11, true),
                            (:vgg11_bn, false),
                            (:vgg11_bn, true),
                            (:vgg13, false),
                            (:vgg13, true),
                            (:vgg13_bn, false),
                            (:vgg13_bn, true),
                            (:vgg16, false),
                            (:vgg16, true),
                            (:vgg16_bn, false),
                            (:vgg16_bn, true),
                            (:vgg19, false),
                            (:vgg19, true),
                            (:vgg19_bn, false),
                            (:vgg19_bn, true),
                        ],
                        vision_transformer => [
                            (:tiny, false),
                            (:small, false),
                            (:base, false),
                            # CI cant handle these
                            # (:large, false),
                            # (:huge, false),
                            # (:giant, false),
                            # (:gigantic, false),
                        ])

@testset "$model_creator" for (model_creator, config) in pairs(models_available)
    @time begin @testset "name = $name & pretrained = $pretrained" for (name, pretrained) in config
        if VERSION <= v"1.7" && pretrained
            @warn "Skipping pretrained models in Julia < 1.7"
            continue
        end
        model, ps, st = model_creator(name; pretrained)
        st = Lux.testmode(st)

        imsize = string(model_creator) == "vision_transformer" ? (256, 256) : (224, 224)
        x = randn(Float32, imsize..., 3, 1)

        @test size(first(model(x, ps, st))) == (1000, 1)

        GC.gc(true)
    end end
end


text:
@testset "public interface" begin
    @test SGFNode === SmartGameFormat.SGFNode
    @test SGFGameTree === SmartGameFormat.SGFGameTree
end

@testset "SGFNode" begin
    @testset "constructor" begin
        n = @inferred SGFNode(:a => 2)
        @test n isa SGFNode
        @test n.properties == Dict(:a => [2])
        @test n.properties isa Dict{Symbol,Vector{Any}}
        @test n == @inferred SGFNode(:a => [2])
        n = @inferred SGFNode()
        @test n isa SGFNode
        @test n.properties == Dict()
        @test n.properties isa Dict{Symbol,Vector{Any}}
        n = @inferred SGFNode(:a => 2, :b => "test")
        @test n == @inferred SGFNode(:a => 2, :b => ["test"])
        @test n == @inferred SGFNode(:a => [2], :b => ["test"])
        @test n isa SGFNode
        @test n.properties == Dict(:a => [2], :b => ["test"])
        @test n.properties isa Dict{Symbol,Vector{Any}}
        n = @inferred SGFNode(Dict(:a => [2], :b => [3]))
        @test n isa SGFNode
        @test n.properties == Dict(:a => [2], :b => [3])
        @test n.properties isa Dict{Symbol,Vector{Any}}
    end

    @testset "equality" begin
        @test SGFNode(:a => 3) == SGFNode(:a => 3)
        @test hash(SGFNode(:a => 3)) == hash(SGFNode(:a => 3))
        @test SGFNode(:a => "3") == SGFNode(:a => "3")
        @test hash(SGFNode(:a => "3")) == hash(SGFNode(:a => "3"))
        @test SGFNode(:a => 3) != SGFNode(:a => "3")
        @test hash(SGFNode(:a => 3)) != hash(SGFNode(:a => "3"))
        @test SGFNode(:a => [3], :b => ["c"]) == SGFNode(:a => [3], :b => ["c"])
        @test hash(SGFNode(:a => [3], :b => ["c"])) == hash(SGFNode(:a => [3], :b => ["c"]))
    end

    @testset "interface" begin
        n = SGFNode()
        @test @inferred(keytype(n)) == Symbol
        @test @inferred(valtype(n)) == Vector{Any}
        @test @inferred(length(n)) == 0
        @test @inferred(get(n, :a, Any[3])) == Any[3]
        @test @inferred(getkey(n, :a, :b)) == :b
        @test @inferred(haskey(n, :a)) == false
        n[:a] = 2
        @test @inferred(length(n)) == 1
        @test @inferred(get(n, :a, Any[3])) == Any[2]
        @test @inferred(getkey(n, :a, :b)) == :a
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(n[:a]) == Any[2]
        n[:a] = 4
        @test @inferred(length(n)) == 1
        @test @inferred(get(n, :a, Any[3])) == Any[4]
        @test @inferred(getkey(n, :a, :b)) == :a
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(n[:a]) == Any[4]
        n[:a] = [5]
        @test @inferred(length(n)) == 1
        @test @inferred(get(n, :a, Any[3])) == Any[5]
        @test @inferred(getkey(n, :a, :b)) == :a
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(n[:a]) == Any[5]
        push!(n, :b=>2)
        @test @inferred(length(n)) == 2
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(haskey(n, :b)) == true
        @test @inferred(n[:b]) == Any[2]
        @test collect(keys(n)) == [:a, :b] || collect(keys(n)) == [:b, :a]
        @test collect(values(n)) == [[5], [2]] || collect(keys(n)) == [[2], [5]]
        push!(n, :b=>3)
        @test @inferred(length(n)) == 2
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(haskey(n, :b)) == true
        @test @inferred(n[:b]) == Any[2, 3]
        push!(n, :b=>[4])
        @test @inferred(length(n)) == 2
        @test @inferred(haskey(n, :a)) == true
        @test @inferred(haskey(n, :b)) == true
        @test @inferred(n[:b]) == Any[4]
    end

    @testset "show" begin
        if Int == Int64
            @test_reference "ref/show_node0.txt" @io2str show(::IO, MIME"text/plain"(), SGFNode())
            @test_reference "ref/show_node1.txt" @io2str show(::IO, MIME"text/plain"(), SGFNode(:a => [2]))
            @test_reference "ref/show_node2.txt" @io2str show(::IO, MIME"text/plain"(), SGFNode(:a => [2], :b => [1.2, 3.2], :c => ["test"]))
            @test_reference "ref/show_node_cpct0.txt" @io2str show(IOContext(::IO, :compact=>true), MIME"text/plain"(), SGFNode())
            @test_reference "ref/show_node_cpct1.txt" @io2str show(IOContext(::IO, :compact=>true), MIME"text/plain"(), SGFNode(:a => [2]))
            @test_reference "ref/show_node_cpct2.txt" @io2str show(IOContext(::IO, :compact=>true), MIME"text/plain"(), SGFNode(:a => [2], :b => [1.2, 3.2], :c => ["test"]))
            @test_reference "ref/show_node_cpct3.txt" @io2str show(IOContext(::IO, :compact=>true), MIME"text/plain"(), SGFNode((Symbol("a$a") => [a] for a in 1:10)...))
        else
            warn("skipping show tests on 32 bit systems")
        end
    end
end

@testset "SGFGameTree" begin
    @test SGFGameTree <: AbstractVector

    @testset "constructor" begin
        t = @inferred SGFGameTree()
        @test t isa SGFGameTree
        @test t.sequence isa Vector{SGFNode}
        @test t.variations isa Vector{SGFGameTree}
        @test length(t.sequence) == 0
        @test length(t.variations) == 0
        t = @inferred SGFGameTree(SGFNode(:a=>2))
        @test t isa SGFGameTree
        @test t.sequence isa Vector{SGFNode}
        @test t.variations isa Vector{SGFGameTree}
        @test length(t.sequence) == 1
        @test length(t.variations) == 0
        @test t.sequence == [SGFNode(:a=>2)]
        t = @inferred SGFGameTree([SGFNode(:b=>3),SGFNode(:c=>"4")])
        @test t isa SGFGameTree
        @test t.sequence isa Vector{SGFNode}
        @test t.variations isa Vector{SGFGameTree}
        @test length(t.sequence) == 2
        @test length(t.variations) == 0
        @test t.sequence == [SGFNode(:b=>3),SGFNode(:c=>"4")]
        t2 = @inferred SGFGameTree([SGFNode(:a=>2)], [t])
        @test t2 isa SGFGameTree
        @test t2.sequence isa Vector{SGFNode}
        @test t2.variations isa Vector{SGFGameTree}
        @test length(t2.sequence) == 1
        @test length(t2.variations) == 1
        @test t2.sequence == [SGFNode(:a=>2)]
        @test t2.variations == [t]
    end

    @testset "equality" begin
        @test SGFGameTree() == SGFGameTree()
        @test hash(SGFGameTree()) == hash(SGFGameTree())
        @test SGFGameTree(SGFNode(:a=>2)) == SGFGameTree(SGFNode(:a=>2))
        @test hash(SGFGameTree(SGFNode(:a=>2))) == hash(SGFGameTree(SGFNode(:a=>2)))
        @test SGFGameTree(SGFNode(:a=>2)) != SGFGameTree(SGFNode(:a=>3))
        @test hash(SGFGameTree(SGFNode(:a=>2))) != hash(SGFGameTree(SGFNode(:a=>3)))
        @test SGFGameTree([SGFNode(:b=>"a")], [SGFGameTree(SGFNode(:a=>2))]) == SGFGameTree([SGFNode(:b=>"a")], [SGFGameTree(SGFNode(:a=>2))])
        @test hash(SGFGameTree([SGFNode(:b=>"a")], [SGFGameTree(SGFNode(:a=>2))])) == hash(SGFGameTree([SGFNode(:b=>"a")], [SGFGameTree(SGFNode(:a=>2))]))
    end

    @testset "interface" begin
        t = SGFGameTree()
        @test @inferred(length(t)) == 0
        @test @inferred(size(t)) == (0,)
        @test_throws BoundsError t[0]
        @test_throws BoundsError t[1]
        push!(t.sequence, SGFNode(:a=>2))
        @test @inferred(length(t)) == 1
        @test @inferred(size(t)) == (1,)
        @test @inferred(t[1]) == SGFNode(:a=>2)
        @test_throws BoundsError t[2]
        t.variations = [SGFGameTree(SGFNode(:b=>"c")), SGFGameTree(SGFNode(:c=>"d"))]
        @test @inferred(length(t)) == 2
        @test @inferred(size(t)) == (2,)
        @test t.sequence == [SGFNode(:a=>2)]
        @test @inferred(t[1]) == SGFNode(:a=>2)
        @test @inferred(t[2]) == SGFNode(:b=>"c")
        t[1] = SGFNode(:AB=>3)
        @test @inferred(length(t)) == 2
        @test @inferred(t[1]) == SGFNode(:AB=>3)
        @test @inferred(t[2]) == SGFNode(:b=>"c")
        @test t.sequence == [SGFNode(:AB=>3)]
        @test t.variations == [SGFGameTree(SGFNode(:b=>"c")), SGFGameTree(SGFNode(:c=>"d"))]
        t[2] = SGFNode(:BC=>3.4)
        @test @inferred(t[1]) == SGFNode(:AB=>3)
        @test @inferred(t[2]) == SGFNode(:BC=>3.4)
        @test t.sequence == [SGFNode(:AB=>3)]
        @test t.variations == [SGFGameTree(SGFNode(:BC=>3.4)), SGFGameTree(SGFNode(:c=>"d"))]
        @test_throws BoundsError t[3]
        @test_throws BoundsError t[3] = SGFNode()
        # include simple show test right here
        if Int == Int64
            @test_reference "ref/show_gametree.txt" @io2str show(::IO, MIME"text/plain"(), t)
        else
            warn("skipping show tests on 32 bit systems")
        end
    end
end


text:
N = 10
NG = N#5

maxW = 2000
minW = 10
tmpStrIO = round.([ Vector(linspace(minW,maxW,N));Vector(linspace(minW,maxW ,N))])
N, StrI,StrO = splitVec(tmpStrIO)
sum(StrI) == sum(StrO)?():error()
Amat = trues(N,N);for n=1:N Amat[n,n] =false end
Amat[StrI.==0,:] = false
Amat[:,StrO.==0] = false
dgpStrIO  =   DegSeq2graphDegSeq(fooErgmDirW1Afixed,tmpStrIO;Amat=Amat)
N, StrI,StrO = splitVec(dgpStrIO)
sum(StrI) == sum(StrO)?():error()
Amat = trues(N,N);for n=1:N Amat[n,n] =false end
Amat[StrI.==0,:] = false
Amat[:,StrO.==0] = false

Mod =  ErgmDirW1Afixed(dgpStrIO,Amat)
S = 0.5*sum(dgpStrIO,1)[1]
Par = ones(2N) * 0.5 /(S/(N^2-N))
expMatTest=expMatrix(Mod,Par,Amat)
firstOrderCond(Mod; strIO =dgpStrIO, parNodesIO = Par,Amat=Amat)

(dgpStrIO .== 0 )
##
estPar,estIt,estMod = estimate(Mod)
StrO .==0
degsO = sumSq(Amat,1)
(degsO.==0)
## primi test
uBndPar = bndPar2uBndPar(Mod,estPar)
sum(abs.(estPar - uBndPar2bndPar(Mod,uBndPar)))

##
using Plots
plotly()
plot(estPar[abs.(estPar).<100])

##
N,parI,parO = splitVec(estPar)
minVal = minimum(parI)
parI = parI - minVal
parO = parO + minVal

sum(parI) - sum(parO)

(sum(parI.<0) + sum(parO.<0)) == 0 ? ():error()

##
expMatEst= expMatrix(Mod,estPar)

sam = sampl(estMod,10000)

tmp = sumSq(sam[N,:,:] .- expMatEst[N,:],1)./sqrt(sumSq( expMatEst[N,:],1).^2)
histogram( tmp )
meanMatSam = squeeze(mean(sam,3),3)
strSam = expValStatsFromMat(estMod,meanMatSam)
strGroupsSam = zeros(2NG) ; for i =1:N strGroupsSam[groupsInds[i]] += strSam[i]; strGroupsSam[NG + groupsInds[i]] += strSam[ N + i] end
relErrSam = maximum(abs.(( strGroupsSam - strGroups )./strGroups )[strGroups.!=0])
absErrZeros = maximum(abs.(( strGroupsSam - strGroups ) )[strGroups.==0])

## Test sui dati
using JLD2
 load_fold =   "/home/Domenico/Dropbox/Dynamic_Networks/data/emid_data/juliaFiles/"
# file_name = "Weekly_eMid_Estimates.jld"
# load_path = load_fold*file_name#
# @load(load_path,estSS,estGas,estSS_1GW,estGas_1GW,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)
#
 N,~,T = size(YeMidWeekly_T)
# ##
# estPar_T = zeros(T,2N)
#
# for t = 1:T
#     Y = YeMidWeekly_T[:,:,t]
#     strIOeMid_t = expValStatsFromMat(fooErgmDirW1,Y)
#     eMidMod = ErgmDirW1(strIOeMid_t)
#
#     estPar,estIt,estMod = estimate(eMidMod)
#     estPar_T[t,:] = estPar
#     N,parI,parO = splitVec(estPar)
#     minVal = minimum(parI)
#     parI = parI - minVal
#     parO = parO + minVal
#
#     println(t)
#      (sum(parI.<0) + sum(parO.<0)) == 0 ? ():error()
# end
#
# estContW_SS_T = estPar_T
# file_name = "Weekly_eMid_Estimates_ContinuousWeights.jld"
# load_path = load_fold*file_name#
# @save(load_path,estContW_SS_T,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)

file_name = "Weekly_eMid_Estimates_ContinuousWeights.jld"
load_path = load_fold*file_name#
@load(load_path,estContW_SS_T,AeMidWeekly_T, YeMidWeekly_T ,weekInd,degsIO_T)
 N,~,T = size(YeMidWeekly_T)

##
estParPos_T=zeros(estContW_SS_T)
for t=1:T
    estPar = estContW_SS_T[t,:]
    N,parI,parO = splitVec(estPar)
    minVal = minimum(parI)
    parI = parI - minVal
    parO = parO + minVal
    estParPos_T[t,:] = [parI;parO]
end

estParPos_T

##
t=1
par = estParPos_T[t,:]
N,parI,parO = splitVec(par)

minimum(parI)
minimum(parO)
sum(parI) - sum(parO)

##
using Plots
plotly()
pdegsI = plot(degsIO_T[:,1:N],title = "  $(round.(squeeze(mean(degsIO_T[:,1:N],1),1),1))")
pparO = plot(estParPos_T[:,N+1:2N])#,title = round.(squeeze(mean(Fitness_T[:,N+1:2N],1),1),1))  $(round.(UMdgp[N+1:2N],1))")
pdegsO = plot(degsIO_T[:,N+1:2N],title =  "  $(round.(squeeze(mean(degsIO_T[:,N+1:2N],1),1),1))")

plot(pdegsI,pdegsO,layout = (2,1),legend=:none,size = (1200,600))

## test IPF function
YallT = meanSq(YeMidWeekly_T,3)
strIOallT = expValStatsFromMat(fooErgmDirW1,YallT)
#zeros on the diagonal
A =  ones(N,N); for i=1:N A[i,i] = 0 end  #Int.(YallT .> 0) #
nnzInds = strIO.!=0
~,nnzIndsI,nnzIndsO = splitVec(nnzInds)

@time  ipfMat,~ = estimateIPFMat(fooErgmDirW1,strIOallT,A)
eMidMod = ErgmDirW1(strIOallT)
@time  estPar,estIt,estMod = estimate(eMidMod)
gammaMat = expMatrix(estMod,estPar)

gammaMat = sortMat(gammaMat )
ipfMat = sortMat(ipfMat)
gammaMatA = (gammaMat.>0)
ipfMatA = ipfMat.>0
sum(gammaMatA) - sum(ipfMatA)
sum(.!ipfMatA.& .!gammaMatA)
plotMat = log.(gammaMat./ipfMat)
 contour(plotMat,fill=true,yflip = true)
##

contour(sortMat(log.(gammaMat)),fill=true,yflip = true)
function sortMat(M::Array{<:Real,2})
    strI = sum(M,2)
    N = length(strI)
    M1 = [strI Vector(1:N) M ]

    tmp = sortrows(M1,rev=true)
    inds = Int.(tmp[:,2])
    Msorted = tmp[:,3:end]
    #M2 = [strI Msorted' ]
    #Msorted = sortrows(M2,rev=true)[:,2:end]'
    Msorted = Msorted[:,inds]

    return Msorted
end

sortMat(gammaMat[1:4,1:4,])
##

#

#


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module Encounter

export
    AbstractEncounterModel,
    AddObserver,

    update,

    generateEncounter,
    getInitialState,
    getNextCommand,
    getTrajectory,

    CorrAEM,
    CorrAEMInitialState,
    CorrAEMCommand,
    generateEncountersToFile,
    setInitialDistributions,
    validate,
    getInitialSample,

    LLAEM,
    LLAEMInitialState,
    LLAEMCommand

import CommonInterfaces: initialize, update
using AbstractEncounterModelImpl

using CorrAEMImpl
using LLAEMImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
"""
Function that returns Utility for given prices and coordinates. Requires precompilation of Utility Cube.
CPU version. Requires precompilation of utils cube in RAM. Accesses values directly.
# Arguments
- `ind_vec::NTuple{M,<:Integer}`: coordinates of Utility.
- `mdl::Model`: mutable struct Model.
"""

function cube_CPU( ind_vec::Vector{Int64}, mdl::Model, additional_args=nothing)

    grid = mdl.GridDef
    status, index = (ind_vec[length(ind_vec)], CartesianIndex(Base.front(Tuple(ind_vec))))
    return(mdl.Utils.qb[ status ][ index ])

end
# function cube_CPU( ind_vec::CartesianIndex, mdl::Model, additional_args=nothing)

#     grid = mdl.GridDef
#     status, index = (ind_vec[length(ind_vec)], CartesianIndex(Base.front(Tuple(ind_vec))))
#     return(mdl.Utils.qb[ status ][ index ])

# end
function cube_CPU( ind_vec::CartesianIndex, status::Int, mdl::Model, additional_args=nothing)

    return(mdl.Utils.qb[ status ][ ind_vec ])

end

"""
Function building cube. CPU version. Requires precompilation of utils cube in RAM. Accesses values directly.
# Arguments
- `mdl::ModelInit`: struct ModelInit.
"""

function cube_CPU_qb(mdl::Model)

    n_dim = tuple(length.(mdl.GridDef.Grids)[1:(end-1)]...)
    cube = Vector( [repeat([-Inf::Float64], n_dim... ), repeat([-Inf::Float64], n_dim... )] )

    for status = mdl.GridDef.Grids[end]

        Threads.@threads for idx = CartesianIndices(n_dim)

            ind_vec = CartesianIndex( idx, status )
            cube[status][idx] =  mdl.HAutility( mdl.HAconstraint( Tuple(ind_vec) , mdl )... )

        end

    end

    return(cube)

end


text:
using Documenter
using BinaryTraits

makedocs(
    modules = [BinaryTraits],
    sitename="BinaryTraits.jl",
    pages = [
        "Motivation" => "index.md",
        "User Guide" => "guide.md",
        "Concepts" => "concepts.md",
        "Reference" => "reference.md",
        "Under the hood" => "design.md",
    ]
)

deploydocs(
    repo = "github.com/tk3369/BinaryTraits.jl.git",
)


text:
function path(parentmx::AbstractMatrix{U}, orig::U, dest::U) where {U<:Integer}
    path = Vector{Tuple{U,U}}()
    currnode = dest
    while currnode != orig
        prevnode = parentmx[orig,currnode]
        push!(path, (prevnode, currnode))
        currnode = prevnode
    end
    reverse(path)
end

function pathto(parents::AbstractVector, dest::U) where {U<:Integer}
    path = Vector{Tuple{U,U}}()
    currnode = dest
    while true #currnode != orig
        prevnode = parents[currnode]
        if currnode == prevnode
            break
        end
        push!(path, (prevnode, currnode))
        currnode = prevnode
    end
    reverse(path)
end


text:
# CPU cross-correlation using a single core
# this was tested on a machine with 32 GB of memory
# This will run out of memory with 16 GB of memory

using SeisIO, SeisNoise, Dates, Glob
fs = 1.
cc_len, cc_step = 2^15, 26768
freqmin = 0.003
freqmax = 0.005
maxlag = 12000.
YEARSEIS = expanduser("~/SeisNoise_paper/DATA/YEARSEIS")
SINGLECORR = expanduser("~/SeisNoise_paper/SINGLECORR")
XML = expanduser("~/SeisNoise_paper/DATA/XML")
files = glob("*seisio",YEARSEIS)
N = length(files)

if !isdir(SINGLECORR)
    mkpath(SINGLECORR)
end

function seis2fft(file::String,XML::String,cc_len::Int,cc_step::Int)
    S = SeisData()
    println("Reading file ",basename(file[1:end-7]))
    read_data!(S,"sac",file)
    # detrend!(S)
    # taper!(S,t_max=400.)
    filtfilt!(S,fl=0.01,fh=0.4,np=3)

    # get xmlfile
    net,sta,loc,chan = split(S[1].id,'.')
    xmlfile = joinpath(XML,net*'.'*sta*".xml")
    Resp = read_meta("sxml",xmlfile)
    translate_resp!(S,Resp)
    remove_resp!(S)
    filtfilt!(S,fh=0.05,rt="Lowpass")
    R = RawData(S[1],cc_len,cc_step)
    demean!(R)
    taper!(R,max_length=1000.)
    return rfft(R)
end

function corrmapper(A::Array{FFTData,1},maxlag::Float64,OUTDIR::String)
        N = size(A,1)
        # copy the current FFT and correlate against all remaining
        for ii = 1:N-1
            for jj = ii+1:N
                mapcc(A[ii],A[jj],maxlag,OUTDIR)
            end
            A[ii] = FFTData()
        end
end

function mapcc(FFT1::FFTData,FFT2::FFTData,maxlag::Float64,OUTDIR::String)
    println("Correlation $(FFT1.name), $(FFT2.name)")
    C = compute_cc(FFT1,FFT2,maxlag)
    if !isnothing(C)
        abs_max!(C)
        stack!(C,allstack=true)
        save_corr(C,OUTDIR)
    else
        println("FAILED $(FFT1.name), $(FFT2.name)")
    end
    return nothing
end

# precompile
Fs = map(seis2fft,files[1:24],fill(XML,24),fill(cc_len,24),fill(cc_step,24))
corrmapper(Fs,maxlag,SINGLECORR)
toremove = glob("*",SINGLECORR)
rm.(toremove)

println("Beginning XCORR $(today())")
println("$(length(files)) total stations to process...")
println("$(length(files) * (length(files)-1) / 2) total correlation to process...")
tALL = @elapsed begin
Fs = map(seis2fft,files,fill(XML,N),fill(cc_len,N),fill(cc_step,N))
corrmapper(Fs,maxlag,SINGLECORR)
end

println("XCORR took $tALL seconds")


text:
module SimpleHypergraphs

using LightGraphs

export Hypergraph, getvertices, gethyperedges, hg_load, hg_save
export add_vertex!, add_hyperedge!, remove_vertex!
export set_vertex_meta!, get_vertex_meta
export set_hyperedge_meta!, get_hyperedge_meta
export BipartiteView, shortest_path
export TwoSectionView

export nhv, nhe
export modularity, randompartition
export AbstractCommunityFinder, CFModularityRandom, CFModularityCNMLike
export findcommunities

include("hypergraph.jl")
include("bipartite.jl")
include("io.jl")
include("twosection.jl")
include("modularity.jl")

end # module


text:
using Base.Dates
using BusinessDays
using Base.Test

bd = BusinessDays

# lets find out minimum and maximum month for easter
y = 1582
easter_1582 = bd.easter_date(Year(1582))
month_min = month(easter_1582)
date_min = easter_1582
month_max = month(easter_1582)
date_max = easter_1582

while y <= 2100
    e_date = bd.easter_date(Year(y))
    m_e_date = month(e_date)
    month_min = min(month_min, m_e_date)
    month_max = max(month_max, m_e_date)

    if month_min == m_e_date
        date_min = e_date
    end

    if month_max == m_e_date
        date_max = e_date
    end

    y += 1
end

print("easter minimum month is $(month_min) on date $(date_min) \neaster maximum month is $(month_max) on date $(date_max)\n")

@test month_min >= 3
@test month_max <= 4


text:
# https://www.codewars.com/kata/56bc28ad5bdaeb48760009b0

removechar(str) = str[2:end-1]

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
using PixelArt
using Test
using TestImages

@testset "PixelArt.jl" begin
    img = testimage("l")
    @test typeof(pixel(img)) <: AbstractArray
    # A wrong path often causes passing nothing.
    @test_throws MethodError pixel(nothing)
end


text:
"""
ARD Transform
```
    v = rand(3)
    tr = ARDTransform(v)
```
Multiply every vector of observation by `v` element-wise
"""
struct ARDTransform{T,N} <: Transform
    v::Vector{T}
end

function ARDTransform(s::T,dims::Integer) where {T<:Real}
    @check_args(ARDTransform, s, s > zero(T), "s > 0")
    ARDTransform{T,dims}(fill(s,dims))
end

function ARDTransform(v::AbstractVector{T}) where {T<:Real}
    @check_args(ARDTransform, v, all(v.>zero(T)), "v > 0")
    ARDTransform{T,length(v)}(v)
end

function set!(t::ARDTransform{T},ρ::AbstractVector{T}) where {T<:Real}
    @assert length(ρ) == dim(t) "Trying to set a vector of size $(length(ρ)) to ARDTransform of dimension $(dim(t))"
    t.v .= ρ
end

dim(t::ARDTransform) = length(t.v)

function apply(t::ARDTransform,X::AbstractMatrix{<:Real};obsdim::Int = defaultobs)
    @boundscheck if dim(t) != size(X,feature_dim(obsdim))
        throw(DimensionMismatch("Array has size $(size(X,!Bool(obsdim-1)+1)) on dimension $(!Bool(obsdim-1)+1)) which does not match the length of the scale transform length , $(dim(t)).")) #TODO Add test
    end
    _transform(t,X,obsdim)
end
apply(t::ARDTransform,x::AbstractVector{<:Real};obsdim::Int=defaultobs) = t.v .* x
_transform(t::ARDTransform,X::AbstractMatrix{<:Real},obsdim::Int=defaultobs) = obsdim == 1 ? t.v'.*X : t.v .* X

Base.isequal(t::ARDTransform,t2::ARDTransform) = isequal(t.v,t2.v)


text:
A = [true false false
     false true false
     false false true]
B = copy(A); B[1,2] = true;
sz_img = (3,3)
sz_img_3 = (3, 3, 3)
m, n = 3, 5

type_list = generate_test_types([Bool], [Gray])

@testset "GenericHausdorff" begin
    @info "test: GenericHausdorff"

    dist = Hausdorff()
    @testset "Hausdorff" begin
        for T in type_list
            @testset "$T" begin
                test_Metric(dist, sz_img, T)
                test_ndarray(dist, sz_img_3, T)
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)

                a = A .|> T
                b = B .|> T
                test_numeric(dist, a, b, T; filename="references/Hausdorff_2d_Gray")
            end
        end
    end
    test_cross_type(dist, A, B, type_list)

    dist = ModifiedHausdorff()
    @testset "ModifiedHausdorff" begin
        for T in type_list
            @testset "$T" begin
                test_SemiMetric(dist, sz_img, T)
                test_ndarray(dist, sz_img_3, T)
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)

                a = A .|> T
                b = B .|> T
                test_numeric(dist, a, b, T; filename="references/ModifiedHausdorff_2d_Gray")
            end
        end
    end
    test_cross_type(dist, A, B, type_list)
end


text:
module TapFiltersTests

using TapFilters, Base.Test

@test size(dx) == (5,5)
@test size(dy) == (5,5)
@test size(dxx) == (5,5)
@test size(dyy) == (5,5)
@test size(dxy) == (5,5)

end


text:
#
# io.jl -
#
# Read and write OI-FITS files.
#
#------------------------------------------------------------------------------

# Union of possible FITS keyword types.
const KeywordTypes = Union{Bool,Int,Cdouble,String}

"""
    OIFITS.column_type(sym)

yields the exact element type in a column of a FITS table whose type is
specified by the letter `sym`.

"""
column_type(sym::Symbol) = (
    sym === :A ? String :
    sym === :I ? Int16 :
    sym === :J ? Int32 :
    sym === :K ? Int64 :
    sym === :L ? Bool :
    sym === :E ? Float32 :
    sym === :D ? Float64 :
    sym === :C ? Complex{Float32} :
    sym === :M ? Complex{Float64} :
    error("unknown FITS column type `", sym, "`"))

# Exceptions.

MissingColumn(col::AbstractString, src::Union{TableHDU,OIDataBlock}) =
    MissingColumn(col, extname(src))

MissingKeyword(key::AbstractString, src::Union{TableHDU,OIDataBlock}) =
    MissingKeyword(key, extname(src))

show(io::IO, ::MIME"text/plain", e::MissingKeyword) =
    print(io, "keyword \"", e.key, "\" not found in FITS extension ", e.ext)

show(io::IO, ::MIME"text/plain", e::MissingColumn) =
    print(io, "column \"", e.col, "\" not found in FITS extension ", e.ext)

#------------------------------------------------------------------------------
# WRITING OF OI-FITS FILES

function write(filename::AbstractString, ds::OIDataSet;
               overwrite::Bool=false, kwds...)
    if !overwrite && isfile(filename)
        error("file \"", filename, "\" already exists, ",
              "use keyword `overwrite=true` to overwrite")
    end
    FITS(filename, "w") do f
        write(f, ds; kwds...)
    end
    return nothing
end

# Writing an OI-FITS file is easy: just write all datablocks.  To make reading
# easier, dependencies are written first.
function write(f::FITS, ds::OIDataSet; quiet::Bool=false)
    write(f, ds.target)
    for db in ds.array
        write(f, db)
    end
    for db in ds.instr
        write(f, db)
    end
    for db in ds.correl
        write(f, db)
    end
    for db in ds.vis
        write(f, db)
    end
    for db in ds.vis2
        write(f, db)
    end
    for db in ds.t3
        write(f, db)
    end
    for db in ds.flux
        write(f, db)
    end
    for db in ds.inspol
        write(f, db)
    end
end

function write(f::FITS, db::OIDataBlock)
    hdr_names = String[]
    hdr_values = Any[]
    hdr_comments = String[]
    col_names = String[]
    col_data = Any[]
    col_units = Dict{String,String}()
    check(db) # recheck!
    for spec in get_format(db; throw_errors=true)
        # Skip undefined fields.
        isdefined(db, spec.symb) || continue
        if ndims(spec) == 0
            # Header keyword.
            push!(hdr_names, spec.name)
            push!(hdr_values, getfield(db, spec.symb))
            if spec.units == ""
                push!(hdr_comments, spec.descr)
            else
                push!(hdr_comments, "["*spec.units*"] "*spec.descr)
            end
        else
            # Column keyword.
            push!(col_names, spec.name)
            push!(col_data, with_eltype(column_type(spec.type),
                                        getfield(db, spec.symb)))
            col_units[spec.name] = spec.units
        end
    end
    write(f, col_names, col_data;
          hdutype = TableHDU,
          name = db.extname,
          header = FITSHeader(hdr_names, hdr_values, hdr_comments),
          units = col_units)
end

function write(f::FITS, db::OI_TARGET)
    hdr_names = String[]
    hdr_values = Any[]
    hdr_comments = String[]
    col_names = String[]
    col_data = Any[]
    col_units = Dict{String,String}()
    check(db) # recheck!
    for spec in get_format(db)
        if ndims(spec) == 0
            # Header keyword.
            push!(hdr_names, spec.name)
            push!(hdr_values, getfield(db, spec.symb))
            if spec.units == empty_string
                push!(hdr_comments, spec.descr)
            else
                push!(hdr_comments, "["*spec.units*"] "*spec.descr)
            end
        else
            # Column keyword.
            push!(col_names, spec.name)
            push!(col_data, get_column(column_type(spec.type),
                                       db, spec.symb))
            col_units[spec.name] = spec.units
        end
    end
    write(f, col_names, col_data;
          hdutype = TableHDU,
          name = db.extname,
          header = FITSHeader(hdr_names, hdr_values, hdr_comments),
          units = col_units)
end

with_eltype(::Type{T}, A::AbstractArray{T,N}) where {T,N} = A
with_eltype(::Type{T}, A::AbstractArray{<:Any,N}) where {T,N} =
    convert(Array{T,N}, A)

#------------------------------------------------------------------------------
# READING OF OI-FITS FILES

# Yields whether an exception was due to a missing FITS keyword.
missing_keyword(ex::CFITSIO.CFITSIOError) = ex.errcode == 202
missing_keyword(ex::Exception) = false

# Yields whether an exception was due to a missing FITS column.
missing_column(ex::CFITSIO.CFITSIOError) = ex.errcode == 219
missing_column(ex::Exception) = false

"""
    OIFITS.read_keyword(T, hdu, key, def=OIFITS.unspecified) -> val

yields the value of keyword `key` in FITS header of `hdu` converted to type
`T`.  If the keyword is not part of the header, then the default value `def` is
returned if specified, otherwise a `MissingKeyword` exception is thrown.  This
method provides some type-stability.

"""
function read_keyword(T::Type{<:KeywordTypes}, hdu::HDU, key::String,
                      def = unspecified)
    try
        val, com = FITSIO.read_key(hdu, key)
        return _convert_keyword(T, key, val)
    catch ex
        if missing_keyword(ex)
            def === unspecified && throw(MissingKeyword(key, hdu))
            return def
        end
        rethrow()
    end
end

for (T, S) in ((Bool, Bool),
               (Int, Integer),
               (Cdouble, AbstractFloat),
               (String, AbstractString))
    @eval _convert_keyword(::Type{T}, key::String, val::$S) where {T<:$T} =
        convert(T, val)
end

@noinline _convert_keyword(::Type{T}, key::String, val::S) where {T,S} =
    error("FITS keyword \"", key, "\" of type ", S,
          " cannot be converted to type ", T)

"""
    OIFITS.read_column(T=Array, hdu, col, def=OIFITS.unspecified) -> val

yields the contents of column `col` in FITS table `hdu` and converted to array
type `T`.  If the column is not part of the table, then the default value `def`
is returned if specified, otherwise a `MissingColumn` exception is thrown.
This method provides some type-stability and add missing leading dimensions as
needed.

"""
read_column(hdu::TableHDU, col::String, def = unspecified) =
    read_column(Array, hdu, col, def)

function read_column(T::Type{<:Array}, hdu::TableHDU, col::String,
                     def = unspecified)
    try
        val = read(hdu, col; case_sensitive=false)
        return _convert_column(T, col, val)
    catch ex
        if missing_column(ex)
            def === unspecified && throw(MissingColumn(col, hdu))
            return def
        end
        rethrow()
    end
end

# Convert column data, provide missing dimensions.
function _convert_column(::Type{T}, col::String, val::T) where {T}
    return val
end

for T in (Bool, Integer, AbstractFloat, Complex{<:AbstractFloat},
          AbstractString)
    @eval begin
        function _convert_column(::Type{Array{T,N}}, col::String,
                                 val::AbstractArray{<:$T,N}) where {T<:$T,N}
            return convert(Array{T,N}, val)
        end
        function _convert_column(::Type{Array{T,2}}, col::String,
                                 val::AbstractArray{<:$T,1}) where {T<:$T}
            # FIXME: add trailing dimension(s) instead?
            return copyto!(Array{T,2}(undef, 1, length(val)), val)
        end
        function _convert_column(::Type{Array{T,3}}, col::String,
                                 val::AbstractArray{<:$T,1}) where {T<:$T}
            # FIXME: add trailing dimension(s) instead?
            return copyto!(Array{T,3}(undef, 1, 1, length(val)), val)
        end
    end

end
@noinline _convert_column(::Type{T}, col::String, val::S) where {T,S} =
    error("FITS table column \"", col, "\" of type ", S,
          " cannot be converted to type ", T)

# To read into an existing data-set, first read the OI-FITS file (to make sure
# it is a consistent data-set), then merge.
function read!(ds::OIDataSet, args::Union{AbstractString,FITS}...; kwds...)
    for arg in args
        merge!(ds, read(OIDataSet, arg; kwds...))
    end
    return ds
end

OIDataSet(args::Union{AbstractString,FITS}...; kwds...) =
    read(OIDataSet, args...; kwds...)

function read(::Type{OIDataSet}, arg::Union{AbstractString,FITS},
              args::Union{AbstractString,FITS}...; kwds...)
    read!(read(OIDataSet, arg; kwds...),  args...; kwds...)
end

read(::Type{OIDataSet}, filename::AbstractString; kwds...) =
    FITS(filename, "r") do f
        read(OIDataSet, f; kwds...)
    end

# Thanks to the implemented methods, reading an OI-FITS file is not too
# difficult.  The dependencies must however be read first and the OI-FITS file
# must be a consistent data-set in itself.
function read(::Type{OIDataSet}, f::FITS; kwds...)
    # Starting with an empty data-set, first read all dependencies, then all
    # data.
    ds = OIDataSet()
    for pass in 1:2
        for i in 2:length(f)
            _read!(ds, f[i], pass; kwds...)
        end
    end
    return ds
end

# Type-stable methods to read a single HDU.
for T in (:OI_TARGET, :OI_ARRAY, :OI_WAVELENGTH, :OI_CORR, :OI_VIS,
          :OI_VIS2, :OI_T3, :OI_FLUX, :OI_INSPOL)
    @eval begin
        $T(hdu::HDU; kwds...) = read($T, hdu; kwds...)
        function read(::Type{$T}, hdu::TableHDU; kwds...)
            extname(hdu) == extname($T) || error(
                "FITS table is not an ", extname($T), " extension")
            return _read($T, hdu; kwds...)
        end
    end
end
read(T::Type{<:OIDataBlock}, hdu::HDU; kwds...) =
    error("FITS HDU is not an ", extname(T), " extension")

# Type-instable methods to read a single HDU. For debugging purposes.
OIDataBlock(hdu::HDU; kwds...) = read(OIDataBlock, hdu; kwds...)
function read(::Type{OIDataBlock}, hdu::TableHDU; kwds...)
    extn = extname(hdu)
    extn == "" && error("FITS HDU has no keyword \"EXTNAME\"")
    return _read(datablock_type(extn), hdu; kwds...)
end
read(::Type{OIDataBlock}, hdu::HDU; kwds...) =
    error("FITS HDU is not a FITS table")

datablock_type(extn::AbstractString) =
    (extn == "OI_TARGET"     ? OI_TARGET     :
     extn == "OI_ARRAY"      ? OI_ARRAY      :
     extn == "OI_WAVELENGTH" ? OI_WAVELENGTH :
     extn == "OI_CORR"       ? OI_CORR       :
     extn == "OI_VIS"        ? OI_VIS        :
     extn == "OI_VIS2"       ? OI_VIS2       :
     extn == "OI_T3"         ? OI_T3         :
     extn == "OI_FLUX"       ? OI_FLUX       :
     extn == "OI_INSPOL"     ? OI_INSPOL     :
     error("\"", extn, "\" is not the name of an OI-FITS extension"))

# skip non-table HDUs
_read!(::OIDataSet, ::HDU, ::Integer; kwds...) = nothing

function _read!(ds::OIDataSet, hdu::TableHDU, pass::Integer; kwds...)
    extn = extname(hdu)
    if extn != ""
        if pass == 1
            # Read dependencies.
            if extn == "OI_TARGET"
                isempty(ds.target.list) || error(
                    "only one OI_TARGET data-block is allowed in ",
                    "a compliant OI-FITS file")
                push!(ds, _read(OI_TARGET, hdu; kwds...))
            elseif extn == "OI_ARRAY"
                push!(ds, _read(OI_ARRAY, hdu; kwds...))
            elseif extn == "OI_WAVELENGTH"
                push!(ds, _read(OI_WAVELENGTH, hdu; kwds...))
            elseif extn == "OI_CORR"
                push!(ds, _read(OI_CORR, hdu; kwds...))
            end
        elseif pass == 2
            # Read data.
            if extn == "OI_VIS"
                push!(ds, _read(OI_VIS, hdu; kwds...))
            elseif extn == "OI_VIS2"
                push!(ds, _read(OI_VIS2, hdu; kwds...))
            elseif extn == "OI_T3"
                push!(ds, _read(OI_T3, hdu; kwds...))
            elseif extn == "OI_FLUX"
                push!(ds, _read(OI_FLUX, hdu; kwds...))
            elseif extn == "OI_INSPOL"
                push!(ds, _read(OI_INSPOL, hdu; kwds...))
            end
        end
    end
end

function _read(T::Type{<:OIDataBlock}, hdu::TableHDU; hack_revn = undef)
    db = T(undef)
    db.revn = _read_revn(T, hdu, hack_revn)
    for spec in get_format(db)
        if spec.rank == 0
            if spec.symb !== :revn
                _read_keyword!(db, Val(spec.symb), hdu, spec)
            end
        else
            _read_column!(db, Val(spec.symb), hdu, spec)
        end
    end
    return db
end

function _read(T::Type{<:OI_TARGET}, hdu::TableHDU; hack_revn = undef)
    # Read keywords.
    revn = _read_revn(T, hdu, hack_revn)
    nrows = read_keyword(Int, hdu, "NAXIS2")

    # Read columns.
    target_id = read_column(Vector{Int16  }, hdu, "TARGET_ID")
    target    = read_column(Vector{String }, hdu, "TARGET")
    raep0     = read_column(Vector{Float64}, hdu, "RAEP0")
    decep0    = read_column(Vector{Float64}, hdu, "DECEP0")
    equinox   = read_column(Vector{Float32}, hdu, "EQUINOX")
    ra_err    = read_column(Vector{Float64}, hdu, "RA_ERR")
    dec_err   = read_column(Vector{Float64}, hdu, "DEC_ERR")
    sysvel    = read_column(Vector{Float64}, hdu, "SYSVEL")
    veltyp    = read_column(Vector{String }, hdu, "VELTYP")
    veldef    = read_column(Vector{String }, hdu, "VELDEF")
    pmra      = read_column(Vector{Float64}, hdu, "PMRA")
    pmdec     = read_column(Vector{Float64}, hdu, "PMDEC")
    pmra_err  = read_column(Vector{Float64}, hdu, "PMRA_ERR")
    pmdec_err = read_column(Vector{Float64}, hdu, "PMDEC_ERR")
    parallax  = read_column(Vector{Float32}, hdu, "PARALLAX")
    para_err  = read_column(Vector{Float32}, hdu, "PARA_ERR")
    spectyp   = read_column(Vector{String }, hdu, "SPECTYP")
    if revn ≥ 2
        category = read_column(Vector{String}, hdu, "CATEGORY")
    end
    list = Vector{OITargetEntry}(undef, nrows)
    for i in 1:nrows
        list[i] = OITargetEntry(
            target_id[i],
            target[i],
            raep0[i],
            decep0[i],
            equinox[i],
            ra_err[i],
            dec_err[i],
            sysvel[i],
            veltyp[i],
            veldef[i],
            pmra[i],
            pmdec[i],
            pmra_err[i],
            pmdec_err[i],
            parallax[i],
            para_err[i],
            spectyp[i],
            (revn ≥ 2 ? category[i] : empty_string),
        )
    end
    return OI_TARGET(list; revn=revn)
end

# Methods to allow for hacking the revision number.
_read_revn(T::Type{<:OIDataBlock}, hdu::TableHDU, hack::Integer) = hack
_read_revn(T::Type{<:OIDataBlock}, hdu::TableHDU, ::typeof(undef)) =
   read_keyword(Int, hdu, "OI_REVN")
_read_revn(T::Type{<:OIDataBlock}, hdu::TableHDU, hack) =
    if applicable(hack, hdu)
        hack(hdu)
    else
        hack(T, read_keyword(Int, hdu, "OI_REVN"))
    end

function _read_keyword!(db::T, ::Val{S}, hdu::TableHDU,
                        spec::FieldDefinition) where {T<:OIDataBlock,S}
    try
        val = read_keyword(fieldtype(T, S), hdu, spec.name)
        setfield!(db, S, val)
    catch ex
        (spec.optional && isa(ex, MissingKeyword)) || rethrow()
    end
    nothing
end

function _read_column!(db::T, ::Val{S}, hdu::TableHDU,
                       spec::FieldDefinition) where {T<:OIDataBlock,S}
    try
        val = read_column(fieldtype(T, S), hdu, spec.name)
        setfield!(db, S, val)
    catch ex
        (spec.optional && isa(ex, MissingColumn)) || rethrow()
    end
    nothing
end


text:
#PhysicalCommunications: Pseudo-Random Bit Sequence Generators/Checkers
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#Integer representation of polynomial x^p1 + x^p2 + x^p3 + x^p4 + 1
_poly(p1::Int) = one(UInt64)<<(p1-1)
_poly(p1::Int, p2::Int) = _poly(p1) + _poly(p2)
_poly(p1::Int, p2::Int, p3::Int, p4::Int) = _poly(p1) + _poly(p2) + _poly(p3) + _poly(p4)

#==Maximum-length Linear-Feedback Shift Register (LFSR) polynomials/taps (XNOR form)
Ref: Alfke, Efficient Shift Registers, LFSR Counters, and Long Pseudo-Random
     Sequence Generators, Xilinx, XAPP 052, v1.1, 1996.==#
const MAXLFSR_POLYNOMIAL = [
	_poly(64,64) #1: not supported
	_poly(64,64) #2: not supported
	_poly(3,2) #3
	_poly(4,3) #4
	_poly(5,3) #5
	_poly(6,5) #6
	_poly(7,6) #7
	_poly(8,6,5,4) #8
	_poly(9,5) #9
	_poly(10,7) #10
	_poly(11,9) #11
	_poly(12,6,4,1) #12
	_poly(13,4,3,1) #13
	_poly(14,5,3,1) #14
	_poly(15,14) #15
	_poly(16,15,13,4) #16
	_poly(17,14) #17
	_poly(18,11) #18
	_poly(19,6,2,1) #19
	_poly(20,17) #20
	_poly(21,19) #21
	_poly(22,21) #22
	_poly(23,18) #23
	_poly(24,23,22,17) #24
	_poly(25,22) #25
	_poly(26,6,2,1) #26
	_poly(27,5,2,1) #27
	_poly(28,25) #28
	_poly(29,27) #29
	_poly(30,6,4,1) #30
	_poly(31,28) #31
	_poly(32,22,2,1) #32
]

#==Types
===============================================================================#
abstract type SequenceGenerator end #Defines algorithm used by sequence() to create a bit sequence
abstract type PRBSGenerator <: SequenceGenerator end #Specifically a pseudo-random bit sequence

#Define supported algorithms:
struct MaxLFSR{LEN} <: PRBSGenerator; end #Identifies a "Maximum-Length LFSR" algorithm

#Define iterator & state objects:
struct MaxLFSR_Iter{LEN,TRESULT} #LFSR "iterator" object
	seed::UInt64 #Initial state (easier to define here than creating state in parallel)
	mask::UInt64 #Store mask value since it cannot easily be statically evaluated.
	len::Int
end
mutable struct MaxLFSR_State{LEN}
	reg::UInt64 #Current state of LFSR register
	bitsrmg::Int #How many bits left to generate
end

#==Constructors
===============================================================================#
"""
    MaxLFSR(reglen::Int)

Construct an object used to identify the Maximum-length LFSR algorithm of a given shift register length, `reglen`.
"""
MaxLFSR(reglen::Int) = MaxLFSR{reglen}()

#==Helper functions:
===============================================================================#
#Find next bit & update state:
function _nextbit(state::MaxLFSR_State{LEN}, polymask::UInt64) where LEN
	msb = UInt64(1)<<(LEN-1) #Statically compiles if LEN is known

	#Mask out all "non-tap" bits:
	reg = state.reg | polymask
	bit = msb
	for j in 1:LEN
		bit = ~xor(reg, bit)
		reg <<= 1
	end
	bit = UInt64((bit & msb) > 0) #Convert resultant MSB to an integer

	state.reg = (state.reg << 1) | bit #Leaves garbage @ bits above LEN
	state.bitsrmg -= 1
	return bit
end

#Core algorithm for sequence() function (no kwargs):
function _sequence(::MaxLFSR{LEN}, seed::UInt64, len::Int, output::DataType) where LEN
	ensure(in(LEN, 3:32), ArgumentError("Invalid LFSR register length, $LEN: 3 <= length <= 32"))
	ensure(LEN < 64,
		OverflowError("Cannot build sequence for MaxLFSR{LEN} with LEN=$LEN >= 64."))
	availbits = (UInt64(1)<<LEN)-UInt64(1) #Available LFSR bits
	ensure((seed & availbits) == seed,
		OverflowError("seed=$seed does not fit in LFSR with register length of $LEN."))
	ensure(len>=0, ArgumentError("Invalid sequence length.  len must be non-negative"))

	poly = UInt64(MAXLFSR_POLYNOMIAL[LEN])
	mask = ~poly
	#==Since `1 XNOR A => A`, we can ignore all taps that are not part of the
	polynomial, simply by forcing all non-tap bits to 1 (OR-ing with `~poly`)
	Thus, we build a mask from `~poly`.
	==#

	return MaxLFSR_Iter{LEN, output}(seed, mask, len)
end

#==Iterator interface:
===============================================================================#

Base.length(iter::MaxLFSR_Iter) = iter.len
Base.eltype(iter::MaxLFSR_Iter{LEN, TRESULT}) where {LEN, TRESULT} = TRESULT
Iterators.IteratorSize(iter::MaxLFSR_Iter) = Base.HasLength()

function Iterators.iterate(iter::MaxLFSR_Iter{LEN, TRESULT}, state::MaxLFSR_State{LEN}) where {LEN, TRESULT}
	if state.bitsrmg < 1
		return nothing
	end
	bit = _nextbit(state, iter.mask)

	return (TRESULT(bit), state)
end

function Iterators.iterate(iter::MaxLFSR_Iter{LEN}) where LEN
	state = MaxLFSR_State{LEN}(iter.seed, iter.len)
	return iterate(iter, state)
end

#==High-level interface
===============================================================================#
"""
    sequence(t::SequenceGenerator; seed::Integer=11, len::Int=-1, output::DataType=Int)

Create an iterable object that defines a bit sequence of length `len`.

Inputs:
  - t: Instance defining type of algorithm used to generate bit sequence.
  - seed: Initial value of register used to build sequence.
  - len: Number of bits in sequence.
  - output: DataType used for sequence elements (typical values are `Int` or `Bool`).

Example returning the first `1000` bits of a PRBS-`31` pattern constructed with the Maximum-length LFSR algorithm seeded with an initial register value of `11`.:

    pattern = collect(sequence(MaxLFSR(31), seed=11, len=1000, output=Bool)).
"""
sequence(t::MaxLFSR; seed::Integer=11, len::Int=-1, output::DataType=Int) =
	_sequence(t, UInt64(seed), len, output)

"""
    sequence_detecterrors(t::SequenceGenerator, v::Array)

Tests validity of bit sequence using sequence generator algorithm.

NOTE: Seeded from first bits of sequence in v.
"""
function sequence_detecterrors(t::MaxLFSR{LEN}, v::Vector{T}) where {LEN, T<:Number}
	ensure(length(v) > LEN, ArgumentError("Pattern vector too short to test validity (must be at least > $LEN)"))

	if Bool == T
		#Will be valid
	else
		for i in 1:length(v)
			ensure(v[i]>=0 && v[i]<=1, ArgumentError("Sequence value not ∉ [0,1] @ index $i.)"))
		end
	end

	#Build seed register for Max LFSR algorithm:
	_errors = similar(v)
	seed = UInt64(0)
	for i in 1:LEN
		seed = (seed << 1) | UInt64(v[i])
		_errors[i] = 0
	end

	#Test for errors in remaining sequence:
	iter = _sequence(t, seed, length(v)-LEN, T)
	state = MaxLFSR_State{LEN}(iter.seed, iter.len)
	for i in (LEN+1):length(_errors)
		(b, state) = iterate(iter, state)
		_errors[i] = convert(T, b!=v[i])
	end

	return _errors
end

#Last line


text:
using Gnuplot
function ribbon(x, y, yerr; scolor = "black", linecolor = "black", leg = "y")
    @gp x y "w l lw 1 lc '$(linecolor)' t '$(leg)'"
    @gp :- x y y .+ yerr "w filledcu lc '$(scolor)' t ''"
    @gp :- "set style fill transparent solid 0.2 noborder"
    @gp :- x y y .- yerr "w filledcu lc '$(scolor)' t ''"
end
# example
x = LinRange(-π,3π,100)
y = sin.(x)
yerr =0.15
ribbon(x, y, yerr; scolor = "orange", linecolor = "red")
save(term ="pngcairo size 600,400", output ="ribbon.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using TupleTools, Base.Cartesian
export loop_einsum, loop_einsum!

"""
    loop_einsum(::EinCode, xs, size_dict)

evaluates the eincode specified by `EinCode` and the tensors `xs` by looping
over all possible indices and calculating the contributions ot the result.
Scales exponentially in the number of distinct index-labels.
"""
function loop_einsum(code::EinCode{ixs, iy}, xs::NTuple{N, AbstractArray{<:Any,M} where M},
                size_dict) where {N,T, ixs, iy}
    size = getindex.(Ref(size_dict), iy)
    loop_einsum!(code, xs, get_output_array(xs, size), size_dict)
end

"""
    loop_einsum!(::EinCode, xs, y, size_dict)

inplace-version of `loop_einsum`, saving the result in a preallocated tensor
of correct size `y`.
"""
function loop_einsum!(code::EinCode{ixs, iy},
                xs::NTuple{N, AbstractArray{<:Any,M} where M},
                y::AbstractArray{T,L}, size_dict) where {N,L,T,IT <: Union{AbstractChar,Integer}, ixs, iy}
    A = einarray(code, xs, size_dict)
    reduce_einarray!(A, y)
end

function reduce_einarray!(A::EinArray, y)
    @inbounds for ind_y in A.OCIS
        iy = subindex(A.y_indexer,ind_y)
        for ind_x in A.ICIS
            y[iy] += map_prod(A.xs, TupleTools.vcat(ind_x.I,ind_y.I), A.x_indexers)
        end
    end
    y
end

function get_output_array(xs::NTuple{N, AbstractArray{<:Any,M} where M}, size) where N
    zeros(promote_type(map(eltype,xs)...), size)
end