text:
"""
    classify(x, pairs; lower=(>=), upper=(<), others=nothing)
    classify(x, pairs...; lower, upper, others)

Create a new array with values in `x` classified by the values in `pairs`.

`pairs` can hold tuples fo values `(2, 3)`, a `Fix2` function e.g. `<=(1)`, a `Tuple`
of `Fix2` e.g. `(>=(4), <(7))`, or an IntervalSets.jl interval, e.g. `3..9` or `OpenInterval(10, 12)`.
`pairs` can also be a `n * 3` matrix where each row is lower bounds, upper bounds, replacement.

If if tuples or a `Matrix` are used, the `lower` and `upper` keywords define
how the lower and upper boundaries are chosen.

If `others` is set other values not covered in `pairs` will be set to that values.

# Arguments

- `x`: a `Raster` or `RasterStack`
- `pairs`: each pair contains a value and a replacement, a tuple of lower and upper
    range and a replacement, or a Tuple of `Fix2` like `(>(x), <(y)`.

# Keywords

- `lower`: Which comparison (`<` or `<=`) to use for lower values, if `Fix2` are not used.
- `upper`: Which comparison (`>` or `>=`) to use for upper values, if `Fix2` are not used.
- `others`: A value to assign to all values not included in `pairs`.
    Passing `nothing` (the default) will leave them unchanged.

# Example

```jldoctest
using Rasters, Plots
A = Raster(WorldClim{Climate}, :tavg; month=1)
classes = <=(15) => 10,
          15..25 => 20,
          25..35 => 30,
          >(35) => 40
classified = classify(A, classes; others=0)
plot(classified; c=:magma)

savefig("build/classify_example.png")
# output
```

![classify](classify_example.png)

$EXPERIMENTAL
"""
function classify end
classify(A::AbstractRaster, p1::Pair, pairs::Pair...; kw...) = classify(A, (p1, pairs...); kw...)
function classify(A::AbstractRaster, pairs::Union{Tuple,AbstractArray};
    filename=nothing, suffix=nothing, lower=(>=), upper=(<),
    others=nothing, missingval=missingval(A)
)
    # Make sure we get a concrete type. Broadcast doesn't always work.
    T = promote_type(_pairs_type(pairs), _others_type(others, A), typeof(missingval))
    # We use `Val{T}` to force type stability through the closure
    valT = Val{T}()
    f(x) = _convert_val(valT, _classify(x, pairs, lower, upper, others, Rasters.missingval(A), missingval))
    A1 = create(filename, T, A; suffix, missingval)
    open(A1; write=true) do O
        broadcast!(f, O, A)
    end
    return A1
end
function classify(xs::AbstractRasterStack, pairs...; suffix=keys(xs), kw...)
    mapargs(xs, suffix) do x, s
        classify(x, pairs; suffix=s, kw...)
    end
end
function classify(xs::AbstractRasterSeries, pairs...; kw...)
    map(x -> classify(x, pairs...; kw...), xs)
end

_pairs_type(pairs) = promote_type(map(eltype ∘ last, pairs)...)
_pairs_type(pairs::AbstractArray{T}) where T = T

_others_type(others, A) = typeof(others)
_others_type(others::Nothing, A) = eltype(A)

_convert_val(::Val{T}, x) where T = convert(T, x)

"""
    classify!(x, pairs...; lower, upper, others)
    classify!(x, pairs; lower, upper, others)

Classify the values of `x` in-place, by the values in `pairs`.

If `Fix2` is not used, the `lower` and `upper` keywords

If `others` is set other values not covered in `pairs` will be set to that values.

# Arguments

- `x`: a `Raster` or `RasterStack`
- `pairs`: each pair contains a value and a replacement, a tuple of lower and upper
    range and a replacement, or a Tuple of `Fix2` like `(>(x), <(y)`.

# Keywords

- `lower`: Which comparison (`<` or `<=`) to use for lower values, if `Fix2` are not used.
- `upper`: Which comparison (`>` or `>=`) to use for upper values, if `Fix2` are not used.
- `others`: A value to assign to all values not included in `pairs`.
    Passing `nothing` (the default) will leave them unchanged.

# Example

`classify!` to disk, with key steps:
- copying a tempory file so we don't write over the RasterDataSources.jl version.
- use `open` with `write=true` to open the file with disk-write permissions.
- use `Float32` like `10.0f0` for all our replacement values and `other`, because
    the file is stored as `Float32`. Attempting to write some other type will fail.

```jldoctest
using Rasters, Plots, RasterDataSources
# Download and copy the file
filename = getraster(WorldClim{Climate}, :tavg; month=6)
tempfile = tempname() * ".tif"
cp(filename, tempfile)
# Define classes
classes = (5, 15) => 10,
          (15, 25) => 20,
          (25, 35) => 30,
          >=(35) => 40
# Open the file with write permission
open(Raster(tempfile); write=true) do A
    classify!(A, classes; others=0)
end
# Open it again to plot the changes
plot(Raster(tempfile); c=:magma)

savefig("build/classify_bang_example.png")
# output
```

![classify!](classify_bang_example.png)

$EXPERIMENTAL
"""
classify!(A::AbstractRaster, p1::Pair, pairs::Pair...; kw...) =
    classify!(A, (p1, pairs...); kw...)
function classify!(A::AbstractRaster, pairs;
    lower=(>=), upper=(<), others=nothing, missingval=missingval(A)
)
    T = promote_type(_pairs_type(pairs), _others_type(others, A), typeof(missingval))
    # We use `Val{T}` to force type stability through the closure
    valT = Val{T}()
    out = broadcast!(A, A) do x
        _convert_val(valT, _classify(x, pairs, lower, upper, others, Rasters.missingval(A), missingval))
    end
    return rebuild(out; missingval=missingval)
end
function classify!(xs::RasterSeriesOrStack, pairs...; kw...)
    map(x -> classify!(x, pairs...; kw...),  xs)
    return xs
end

# _classify
# Classify single values
function _classify(x, pairs, lower, upper, others, oldmissingval, newmissingval)
    isequal(x, oldmissingval) && return newmissingval
    # Use a fold instead of a loop, for type stability
    init = (false, nothing)
    found, foundval = foldl(pairs; init) do (found, foundval), (find, replace)
        if !found && _compare(find, x, lower, upper)
            (true, replace)
        else
            (found, foundval)
        end
    end
    if found
        return foundval
    else
        return isnothing(others) ? x : others
    end
end
function _classify(x, pairs::AbstractArray, lower, upper, others, oldmissingval, newmissingval)
    isequal(x, oldmissingval) && return newmissingval
    found = false
    local foundval
    if ndims(pairs) != 2 || !(size(pairs, 2) in (2, 3))
        throw(ArgumentError("pairs must be a N*2 or N*3 matrix, or Pair"))
    elseif size(pairs, 2) == 2
        for i in 1:size(pairs, 1)
            find = pairs[i, 1]
            if _compare(find, x, lower, upper)
                foundval = pairs[i, 2]
                found = true
                break
            end
        end
    elseif size(pairs, 2) == 3
        for i in 1:size(pairs, 1)
            find = pairs[i, 1], pairs[i, 2]
            if _compare(find, x, lower, upper)
                foundval = pairs[i, 3]
                found = true
                break
            end
        end
    end
    if found
        return foundval
    else
        return isnothing(others) ? x : others
    end
end

_compare(find, x, lower, upper) = find === x
_compare(find::Base.Fix2, x, lower, upper) = find(x)
_compare((l, u)::Tuple, x, lower, upper) = lower(x, l) && upper(x, u)
_compare((l, u)::Tuple{<:Base.Fix2,<:Base.Fix2}, x, lower, upper) = l(x) && u(x)
_compare(interval::LA.IntervalSets.Interval, x, lower, upper) = x in interval


text:
"""
    TrainingData

Abstract type for representing training data.
"""
abstract type TrainingData end

"""
    ExactTrainingData(y, x, xm, K, Km, xKinv) <: TrainingData

Create an object that contains the training data when using the full Gram matrix
K.

# Properties
- `y::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Features for
  training data [Q,T] or \\[T\\] (if Q = 1)
- `x::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Latent
  parameters for training data [L,T] or \\[T\\] (if L = 1)
- `xm::Union{<:Real,<:AbstractVector{<:Real}}`: Mean of latent parameters
  \\[L\\] or scalar (if L = 1)
- `K::AbstractMatrix{<:Real}`: De-meaned (both rows and columns) Gram matrix of
  the kernel evaluated on the training data features [T,T]
- `Km::AbstractVector{<:Real}`: Row means of `K` (before de-meaning) \\[T\\]
- `xKinv::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: `x` times
  the regularized inverse of `K` [L,T] or \\[T\\] (if L = 1)
- `Q::Integer`: Number of training features
- `L::Integer`: Number of latent parameters
- `T::Integer`: Number of training points
"""
struct ExactTrainingData{
    T1<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    T2<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    T3<:Union{<:Real,<:AbstractVector{<:Real}},
    T4<:AbstractMatrix{<:Real},
    T5<:AbstractVector{<:Real},
    T6<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}
} <: TrainingData
    y::T1
    x::T2
    xm::T3
    K::T4
    Km::T5
    xKinv::T6
end

Base.getproperty(data::ExactTrainingData, s::Symbol) = begin
    if s == :Q
        y = getfield(data, :y)
        return ndims(y) == 1 ? 1 : size(y, 1)
    elseif s == :L
        x = getfield(data, :x)
        return ndims(x) == 1 ? 1 : size(x, 1)
    elseif s == :T
        return length(getfield(data, :Km))
    else
        return getfield(data, s)
    end
end

"""
    RFFTrainingData(freq, phase, zm, xm, Czz, Cxz, CxzCzzinv) <: TrainingData

Create an object that contains the training data when using an approximation of
the Gram matrix K using random Fourier features.

# Properties
- `freq::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Random
  frequency values for random Fourier features [H,Q] or \\[H\\] (if Q = 1)
- `phase::AbstractVector{<:Real}`: Random phase values for random Fourier
  features \\[H\\]
- `zm::AbstractVector{<:Real}`: Mean of feature maps \\[H\\]
- `xm::Union{<:Real,<:AbstractVector{<:Real}}`: Mean of latent parameters
  \\[L\\] or scalar (if L = 1)
- `Czz::AbstractMatrix{<:Real}`: Auto-covariance matrix of feature maps [H,H]
- `Cxz::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`:
  Cross-covariance matrix between latent parameters and feature maps [L,H] or
  \\[H\\] (if L = 1)
- `CxzCzzinv::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: `Cxz`
  times the regularized inverse of `Czz` [L,H] or \\[H\\] (if L = 1)
- `Q::Integer`: Number of training features
- `L::Integer`: Number of latent parameters
- `H::Integer`: Kernel approximation order
"""
struct RFFTrainingData{
    T1<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    T2<:AbstractVector{<:Real},
    T3<:AbstractVector{<:Real},
    T4<:Union{<:Real,<:AbstractVector{<:Real}},
    T5<:AbstractMatrix{<:Real},
    T6<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}},
    T7<:Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}
} <: TrainingData
    freq::T1
    phase::T2
    zm::T3
    xm::T4
    Czz::T5
    Cxz::T6
    CxzCzzinv::T7
end

Base.getproperty(data::RFFTrainingData, s::Symbol) = begin
    if s == :Q
        freq = getfield(data, :freq)
        return ndims(freq) == 1 ? 1 : size(freq, 2)
    elseif s == :L
        return length(getfield(data, :xm))
    elseif s == :H
        return length(getfield(data, :zm))
    else
        return getfield(data, s)
    end
end

"""
    train(T, xDists, [νDists,] noiseDist, signalModels, kernel, ρ)

Train PERK using simulated training data.

# Arguments
- `T::Integer`: Number of training points
- `xDists`: Distributions of latent parameters \\[L\\] or scalar (if L = 1);
  `xDists` can be any object such that `rand(xDists, ::Integer)` is defined (or
  a collection of such objects)
- `νDists`: Distributions of known parameters \\[K\\] or scalar (if K = 1);
  `νDists` can be any object such that `rand(νDists, ::Integer)` is defined (or
  a collection of such objects); omit this parameter if K = 0
- `noiseDist`: Distribution of noise (assumes same noise distribution for both
  real and imaginary channels in complex case); `noiseDist` can be any object
  such that `rand(noiseDist, ::Integer)` is defined
- `signalModels::Union{<:Function,<:AbstractVector{<:Function}}`: Signal models
  used to generate noiseless data \\[numSignalModels\\]; each signal model
  accepts as inputs L latent parameters (scalars) first, then K known parameters
  (scalars); user-defined parameters (e.g., scan parameters in MRI) should be
  built into the signal model
- `kernel::Kernel`: Kernel to use
- `ρ::Real`: Tikhonov regularization parameter

# Return
- `trainData::TrainingData`: TrainingData object to be passed to `perk`
"""
function train(
    T::Integer,
    xDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}},
    kernel::Kernel,
    ρ::Real
)

    (y, x) = generatenoisydata(T, xDists, noiseDist, signalModels)

    eltype(y) <: Complex && (y = complex2real(y))

    return krr_train(x, y, kernel, ρ)

end

function train(
    T::Integer,
    xDists,
    νDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractArray{<:Function,1}},
    kernel::Kernel,
    ρ::Real
)

    (y, x, ν) = generatenoisydata(T, xDists, νDists, noiseDist, signalModels)

    eltype(y) <: Complex && (y = complex2real(y))

    q = combine(y, ν) # [D+K,T]

    return krr_train(x, q, kernel, ρ)

end

"""
    complex2real(y)

Split complex data into real and imaginary parts.
"""
function complex2real(
    y::Complex
)

    return reshape([real(y), imag(y)], :, 1) # [D=2,T=1]

end

function complex2real(
    y::AbstractVector{<:Complex} # [T]
)

    return [real.(y) imag.(y)]' # [D=2,T]

end

function complex2real(
    y::AbstractMatrix{<:Complex} # [D,T]
)

    return reduce(vcat, (complex2real(@view(y[i,:])) for i = 1:size(y, 1))) # [2D,T]

end

"""
    combine(y, ν)

Combine the output of the signal models with the known parameters.
"""
function combine(
    y::AbstractVector{<:Real}, # [T]
    ν::AbstractVector{<:Real} # [T]
)

    return transpose([y ν])

end

function combine(
    y::AbstractMatrix{<:Real}, # [D,T]
    ν::AbstractVector{<:Real} # [T]
)

    return [y; transpose(ν)]

end

function combine(
    y::AbstractVector{<:Real}, # [T]
    ν::AbstractMatrix{<:Real} # [K,T]
)

    return [transpose(y); ν]

end

function combine(
    y::AbstractMatrix{<:Real}, # [D,T]
    ν::AbstractMatrix{<:Real} # [K,T]
)

    return [y; ν]

end

function combine(
    y::Real,
    ν::Real
)

    return transpose([y ν])

end

function combine(
    y::AbstractMatrix{<:Real}, # [D,1]
    ν::Real
)

    return [y; ν]

end

function combine(
    y::Real,
    ν::AbstractMatrix{<:Real} # [K,1]
)

    return [y; ν]

end

"""
    generatenoisydata(N, xDists, [νDists,] noiseDist, signalModels)

Generate noisy data from unknown (and possibly known) parameter distributions.

# Arguments
- `N::Integer`: Number of data points
- `xDists`: Distributions of latent parameters \\[L\\] or scalar (if L = 1);
  `xDists` can be any object such that `rand(xDists, ::Integer)` is defined (or
  a collection of such objects)
- `νDists`: Distributions of known parameters \\[K\\] or scalar (if K = 1);
  `νDists` can be any object such that `rand(νDists, ::Integer)` is defined (or
  a collection of such objects); omit this parameter if K = 0
- `noiseDist`: Distribution of noise (assumes same noise distribution for both
  real and imaginary channels in complex case); `noiseDist` can be any object
  such that `rand(noiseDist, ::Integer)` is defined
- `signalModels::Union{<:Function,<:AbstractVector{<:Function}}`: Signal models
  used to generate noiseless data \\[numSignalModels\\]; each signal model
  accepts as inputs L latent parameters (scalars) first, then K known parameters
  (scalars); user-defined parameters (e.g., scan parameters in MRI) should be
  built into the signal model

# Return
- `y::Union{<:AbstractVector{<:Number},<:AbstractMatrix{<:Number}}`: Output
  data of all the (simulated) signals [D,N] or \\[N\\] (if D = 1)
- `x::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Randomly
  generated latent parameters [L,N] or \\[N\\] (if L = 1)
- `ν::Union{<:AbstractVector{<:Real},<:AbstractMatrix{<:Real}}`: Randomly
  generated known parameters [K,N] or \\[N\\] (if K = 1); not returned if
  `νDists` is omitted
"""
function generatenoisydata(
    N::Integer,
    xDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}}
)

    # Sample the distributions of latent parameters
    if xDists isa AbstractVector
        x = rand.(xDists, N) # [L][N]
    else
        x = rand(xDists, N) # [N]
    end

    # Generate the data
    y = _evaluate_signalModels(x, signalModels)
    size(y, 1) == 1 && (y = vec(y)) # [N] (if D = 1)
    addnoise!(y, noiseDist)

    # Reshape x
    xDists isa AbstractVector && (x = transpose(reduce(hcat, x))) # [L,N]

    # Return simulated data and random latent parameters
    return (y, x)

end

function generatenoisydata(
    N::Integer,
    xDists,
    νDists,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}}
)

    # Sample the distributions of latent and known parameters
    if xDists isa AbstractVector
        x = rand.(xDists, N) # [L][N]
    else
        x = rand(xDists, N) # [N]
    end
    if νDists isa AbstractVector
        ν = rand.(νDists, N) # [K][N]
    else
        ν = rand(νDists, N) # [N]
    end

    # Generate the data
    y = _evaluate_signalModels(x, ν, signalModels)
    size(y, 1) == 1 && (y = vec(y)) # [N] (if D = 1)
    addnoise!(y, noiseDist)

    # Reshape x and ν
    xDists isa AbstractVector && (x = transpose(reduce(hcat, x))) # [L,N]
    νDists isa AbstractVector && (ν = transpose(reduce(hcat, ν))) # [K,N]

    # Return simulated data and random latent and known parameters
    return (y, x, ν)

end

function _evaluate_signalModels(
    x,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}}
)

    if signalModels isa AbstractVector
        if x isa AbstractVector{<:AbstractVector}
            y = reduce(vcat, (_hcat(signalModels[i].(x...))
                for i = 1:length(signalModels)))
        else
            y = reduce(vcat, (_hcat(signalModels[i].(x))
                for i = 1:length(signalModels)))
        end
    else
        if x isa AbstractVector{<:AbstractVector}
            y = _hcat(signalModels.(x...))
        else
            y = _hcat(signalModels.(x))
        end
    end

    return y

end

function _evaluate_signalModels(
    x,
    ν,
    signalModels::Union{<:Function,<:AbstractVector{<:Function}}
)

    if signalModels isa AbstractVector
        if x isa AbstractVector{<:AbstractVector}
            if ν isa AbstractVector{<:AbstractVector}
                y = reduce(vcat, (_hcat(signalModels[i].(x..., ν...))
                    for i = 1:length(signalModels)))
            else
                y = reduce(vcat, (_hcat(signalModels[i].(x..., ν))
                    for i = 1:length(signalModels)))
            end
        else
            if ν isa AbstractVector{<:AbstractVector}
                y = reduce(vcat, (_hcat(signalModels[i].(x, ν...))
                    for i = 1:length(signalModels)))
            else
                y = reduce(vcat, (_hcat(signalModels[i].(x, ν))
                    for i = 1:length(signalModels)))
            end
        end
    else
        if x isa AbstractVector{<:AbstractVector}
            if ν isa AbstractVector{<:AbstractVector}
                y = _hcat(signalModels.(x..., ν...))
            else
                y = _hcat(signalModels.(x..., ν))
            end
        else
            if ν isa AbstractVector{<:AbstractVector}
                y = _hcat(signalModels.(x, ν...))
            else
                y = _hcat(signalModels.(x, ν))
            end
        end
    end

    return y

end

# Use _hcat to have type stability, because reduce(hcat, x) is type unstable
# if x is a AbstractVector (returns scalar if length(x) == 1, else 2D array)
_hcat(x::AbstractVector{<:AbstractVector}) = reduce(hcat, x)
_hcat(x::AbstractVector) = x

"""
    addnoise!(y, noiseDist)

Add noise to `y`. If elements of `y` are complex-valued, then add independent
noise to both the real and imaginary parts of `y`.
"""
function addnoise!(y, noiseDist)

    if eltype(y) <: Complex
        n = complex.(rand(noiseDist, size(y)...), rand(noiseDist, size(y)...))
    else
        n = rand(noiseDist, size(y)...) # [D,N]
    end
    y .+= n # [D,N]

end

ext: jl


text:
using Test

using NewsLookout

@test NewsLookout.greet() == "-- The NewsLookout Web Scraping App ---"

@test [1, 2] + [2, 1] == [3, 3]


text:
@testset "env-basic" begin
    mds = raw"""
        \newenvironment{aaa}[5]{pre}{post}
        \begin{aaa}
        bbb
        \end{aaa}
        """
    tokens = F.find_tokens(mds, F.MD_TOKENS, F.MD_1C_TOKENS)

    @test has(tokens, :LX_NEWENVIRONMENT)
    @test has(tokens, :CAND_LX_BEGIN)
    @test has(tokens, :CAND_LX_END)

    blocks, tokens = F.find_all_ocblocks(tokens, F.MD_OCB2)

    num_braces_orig = length(filter(b -> b.name == :LXB, blocks))

    @test num_braces_orig == 5

    F.find_lxenv_delims!(tokens, blocks)

    @test has(tokens, :LX_BEGIN)
    @test has(tokens, :LX_END)
    @test !has(tokens, :CAND_LX_BEGIN)
    @test !has(tokens, :CAND_LX_END)

    num_braces_post = length(filter(b -> b.name == :LXB, blocks))

    @test num_braces_post == 3

    @test tokens[3].name == :LX_BEGIN
    @test F.envname(tokens[3]) == "aaa"
    @test tokens[6].name == :LX_END
    @test F.envname(tokens[6]) == "aaa"

    lxdefs, tokens, braces, blocks = F.find_lxdefs(tokens, blocks)

    @test lxdefs[1].name == "aaa"
    @test lxdefs[1].narg == 5
    @test lxdefs[1].def.first == "pre"
    @test lxdefs[1].def.second == "post"
end

@testset "env-1" begin
    s = raw"""
        \newenvironment{aaa}{pre}{post}
        \begin{aaa}
        bbb
        \end{aaa}""" |> fd2html
    @test s // "pre bbb post"
    s = raw"""
        \newenvironment{aaa}[1]{pre:#1}{post:#1}
        \begin{aaa}{00}
        bbb
        \end{aaa}""" |> fd2html
    @test s // "pre: 00 bbb post: 00"
end

@testset "env-redef" begin
    s = raw"""
        \newcommand{\abc}{123}
        \abc
        \newcommand{\abc}{321}
        \abc
        """ |> fd2html
    @test s // "123\n321"

    s = raw"""
        \newcommand{\abc}{123}
        \abc
        \newcommand{\abc}[1]{321:#1}
        \abc{aa}
        """ |> fd2html
    @test s // "123\n321: aa"

    s = raw"""
        \newenvironment{aaa}{pre}{post}
        \begin{aaa}
        bbb
        \end{aaa}

        ---

        \newenvironment{aaa}{PRE}{POST}
        \begin{aaa}
        ccc
        \end{aaa}
        """ |> fd2html
    @test s // "pre bbb post\n
\nPRE ccc POST"

    s = raw"""
        \newenvironment{aaa}{pre}{post}
        \begin{aaa}
        bbb
        \end{aaa}

        ---

        \newenvironment{aaa}[1]{pre:#1}{post:#1}
        \begin{aaa}{00}
        bbb
        \end{aaa}
        """ |> fd2html
    @test s // "pre bbb post\n
\npre: 00 bbb post: 00"
    s = raw"""
        \newenvironment{aaa}{pre}{post}
        \newenvironment{bbb}[2]{abc:#1}{def:#2}
        ABC
        \begin{aaa}
        A
        \begin{bbb}{00}{11}
        B
        \end{bbb}
        C
        \end{aaa}
        """ |> fd2html
    @test s // "
ABC
\npre A abc: 00 B def: 11 C post"
end

@testset "env-maths" begin
    s = raw"""
        AA
        \begin{align}
        A &= B \\
        C &= D+E
        \end{align}
        BB
        """ |> fd2html
    @test s // raw"""

AA

        \[\begin{aligned}
        A &= B \\
        C &= D+E
        \end{aligned}\]

BB

        """
end

@testset "env-nest" begin
    s = raw"""
        \newenvironment{aaa}{AA}{AB}
        \newenvironment{bbb}{BA}{BB}
        00
        \begin{aaa}
        11
        \begin{bbb}
        22
        \end{bbb}
        33
        \end{aaa}
        44
        """ |> fd2html
    @test s // "
00
\nAA 11 BA 22 BB 33 AB\n
44
"
end

@testset "env-errors" begin
    s = raw"""
        \newenvironment{aaa}
        hello
        """
    @test_throws F.LxDefError (s |> fd2html)
    s = raw"""
        \newenvironment{aaa}[d]{a#1}{b#1}
        hello
        """
    @test_throws F.LxDefError (s |> fd2html)
    s = raw"""
        \begin{error}
        A
        \end{error}
        """
    @test_throws F.LxObjError (s |> fd2html)
    s = raw"""
        \begin A \end
        """
    @test_throws F.LxObjError (s |> fd2html)
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module Export

import ACE, JuLIP
using ACE.RPI: BasicPSH1pBasis, PSH1pBasisFcn
using ACE: PIBasis, PIBasisFcn, PIPotential
using ACE.OrthPolys: TransformedPolys
using ACE: rand_radial, cutoff, numz
using JuLIP: energy, bulk, i2z, z2i, chemical_symbol
using ACE.PairPotentials: RepulsiveCore, PolyPairPot

function export_ace(fname::AbstractString, args...;  kwargs...)
   fptr = open(fname; write=true)
   export_ace(fptr, args...; kwargs...)
   close(fptr)
end

function export_ace(fptr::IOStream, Pr::TransformedPolys; ntests=0, kwargs...)
   p = Pr.trans.p
   r0 = Pr.trans.r0
   xr = Pr.J.tr
   xl = Pr.J.tl
   pr = Pr.J.pr
   pl = Pr.J.pl
   rcut = cutoff(Pr)
   maxn = length(Pr)

   println(fptr, "transform parameters: p=$(p) r0=$(r0)")
   println(fptr, "cutoff parameters: rcut=$rcut xl=$xl xr=$xr pl=$pl pr=$pr")
   println(fptr, "recursion coefficients: maxn=$(maxn)")
   for n = 1:maxn
      println(fptr, " $(Pr.J.A[n]) $(Pr.J.B[n]) $(Pr.J.C[n])")
   end

   # save some tests
   println(fptr, "tests: ntests=$(ntests)")
   for itest = 1:ntests
      r = ACE.rand_radial(Pr)
      P = ACE.evaluate(Pr, r)
      dP = ACE.evaluate_d(Pr, r)
      println(fptr, " r=$(r)")
      for n = 1:length(P)
         println(fptr, " $(P[n]) $(dP[n])")
      end
   end
end

function export_ace(fptr::IOStream, Vpair::PolyPairPot; kwargs...)
   println(fptr, "begin polypairpot")
   # this exports everything except the coefficients
   export_ace(fptr, Vpair.basis.J; kwargs...)
   println(fptr, "coefficients")
   c = Vpair.coeffs
   for n = 1:length(c)
      println(fptr, "$(c[n])")
   end
   println(fptr, "end polypairpot")
end

function export_ace(fptr::IOStream, Vrep::RepulsiveCore; kwargs...)
   println(fptr, "begin repulsive potential")
   # export the pair potential part
   export_ace(fptr, Vrep.Vout; kwargs...)
   println(fptr, "spline parameters")
   println(fptr, "   e_0 + B  exp(-A*(r/ri-1)) * (ri/r)")
   Vin = Vrep.Vin[1]
   println(fptr, "ri=$(Vin.ri)")
   println(fptr, "e0=$(Vin.e0)")
   println(fptr, "A=$(Vin.A)")
   println(fptr, "B=$(Vin.B)")
   println(fptr, "end repulsive potential")
end

function export_ace(fptr::IOStream, V::PIPotential, Vpair=nothing, V0=nothing; kwargs...)
   @assert numz(V) == 1
   inner = V.pibasis.inner[1]
   coeffs = V.coeffs[1]
   # sort the basis functions into groups the way the ace evaluator wants it
   groups = _basis_groups(inner, coeffs)
   lmax = maximum(maximum(g["l"]) for g in groups)

   sym = chemical_symbol(i2z(V, 1))
   haspair = (Vpair == nothing) ? "f" : "t"

   # header
   println(fptr, "nelements=1")
   println(fptr, "elements: $sym")
   println(fptr, "")
   println(fptr, "lmax=$lmax")
   println(fptr, "")
   println(fptr, "2 FS parameters:  1.000000 1.000000")
   println(fptr, "core energy-cutoff parameters: 100000.000000000000000000 1.000000000000000000")

   # write E0 values
   if V0 == nothing
      println(fptr, "E0: 0.0")
   else
      println(fptr, "E0: $(V0(sym))")
   end

   println(fptr, "")
   println(fptr, "radbasename=ACE.jl.Basic")
   println(fptr, "haspair: $(haspair)")
   println(fptr, "")

   # export_ace the radial basis
   export_ace(fptr, V.pibasis.basis1p.J; kwargs...)
   println(fptr, "")

   if haspair == "t"
      export_ace(fptr, Vpair)
      println(fptr, "")
   end

   # header
   rankmax = maximum(length(g["l"]) for g in groups)
   println(fptr, "rankmax=$rankmax")
   println(fptr, "ndensitymax=1")
   println(fptr, "")

   # header for basis list pair contributions
   println(fptr, "num_c_tilde_max=$(length(groups))")
   num_ms_combinations_max = maximum( length(g["M"]) for g in groups )
   println(fptr, "num_ms_combinations_max=$(num_ms_combinations_max)")

   # write the pair basis groups
   total_basis_size_rank1 = sum( (length(g["l"]) ==  1) for g in groups )
   println(fptr, "total_basis_size_rank1: $(total_basis_size_rank1)")
   for i = 1:total_basis_size_rank1
      g = groups[i]
      _write_group(fptr, g)
   end

   # write the rest
   total_basis_size = length(groups) - total_basis_size_rank1
   println(fptr, "total_basis_size: $(total_basis_size)")
   for i = (total_basis_size_rank1+1):length(groups)
      g = groups[i]
      _write_group(fptr, g)
   end

end

function _write_group(fptr, g)
   order = length(g["l"])
   println(fptr, "ctilde_basis_func: rank=$(order) ndens=1 mu0=0 mu=(" * " 0 "^order * ")")
   println(fptr, "n=(" * prod(" $(ni) " for ni in g["n"]) * ")")
   println(fptr, "l=(" * prod(" $(li) " for li in g["l"]) * ")")
   println(fptr, "num_ms=$(length(g["M"]))")
   for (m, c) in zip(g["M"], g["C"])
      c_ace = c / (4*π)^(order/2)
      println(fptr, "<" * prod(" $(mi) " for mi in m) * ">:  $(c_ace)")
   end
end

function export_ace_tests(fname::AbstractString, V, ntests = 1;
                          nrepeat = 3, pert=0.05,
                          s = JuLIP.chemical_symbol(i2z(V, 1)))
   at = bulk(s, cubic=true) * nrepeat
   JuLIP.set_pbc!(at, false)
   r0 = JuLIP.rnn(s)
   for n = 1:ntests
      JuLIP.rattle!(at, pert * r0)
      E = energy(V, at)
      _write_test(fname * "_$n.dat", JuLIP.positions(at), E)
   end
   return at
end

function _write_test(fname, X, E)
   fptr = open(fname; write=true)
   println(fptr, "E = $E")
   println(fptr, "natoms = $(length(X))")
   println(fptr, "# type x y z")
   for n = 1:length(X)
      println(fptr, "0 $(X[n][1]) $(X[n][2]) $(X[n][3])")
   end
   close(fptr)
end

function _basis_groups(inner, coeffs)
   NL = []
   M = []
   C = []
   for b in keys(inner.b2iAA)
      if coeffs[ inner.b2iAA[b] ] != 0
         push!(NL, ( [b1.n for b1 in b.oneps], [b1.l for b1 in b.oneps] ))
         push!(M, [b1.m for b1 in b.oneps])
         push!(C, coeffs[ inner.b2iAA[b] ])
      end
   end
   ords = length.(M)
   perm = sortperm(ords)
   NL = NL[perm]
   M = M[perm]
   C = C[perm]
   @assert issorted(length.(M))
   bgrps = []
   alldone = fill(false, length(NL))
   for i = 1:length(NL)
      if alldone[i]; continue; end
      nl = NL[i]
      Inl = findall(NL .== Ref(nl))
      alldone[Inl] .= true
      Mnl = M[Inl]
      Cnl = C[Inl]
      pnl = sortperm(Mnl)
      Mnl = Mnl[pnl]
      Cnl = Cnl[pnl]
      push!(bgrps, Dict("n" => nl[1], "l" => nl[2],
                        "M" => Mnl, "C" => Cnl))
   end
   return bgrps
end

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/reverse_pass.jl
# Functions used to compute reverse pass of nonlinear functions.
#############################################################################

function r_init!(t::Relax, g::DAT, b::RelaxCache{V,N,T}) where {V,N,T<:RelaxTag}
    if !is_num(b, 1)
        b[1] = set(b, 1) ∩ g.sink_bnd
    end
    return !isempty(z)
end

function rprop!(t::Relax, v::Variable, g::DAT, c::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(c, i)
    z = var_set(MC{N,T}, rev_sparsity(g, i, k), x, x, lbd(c, i), ubd(c, i))
    if first_eval(c)
        z = z ∩ interval(c, k)
    end
    c[k] = z
    return !isempty(z)
end

function rprop!(t::Relax, v::Subexpression, g::DAT, c::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    store_subexpression!(c, set(c, k), first_index(g, k))
    return true
end

const MAX_ASSOCIATIVE_REVERSE = 6

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = x + y` which updates x & y.
"""
function rprop_2!(t::Relax, v::Val{PLUS}, g::DAT, c::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    is_num(c, k) && return true
    x = child(g, 1, k)
    y = child(g, 2, k)
    if xset_ynum(b, x, y)
        b, a, q = IntervalContractors.plus_rev(interval(c, k), interval(c, x), num(c, y))
    elseif xnum_yset(b, x, y)
        b, a, q = IntervalContractors.plus_rev(interval(c, k), num(c, x), interval(c, y))
    else
        b, a, q = IntervalContractors.plus_rev(interval(c, k), interval(c, x), interval(c, y))
    end
        if !is_num(c, x)
            isempty(a) && return false
            b[x] = MC{N,T}(a)
        end
        if !is_num(c, y)
            isempty(q) && return false
            b[y] = MC{N,T}(q)
        end
    return true
end

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = +(x,y,z...)` which updates x, y, z and so on.
"""
function rprop_n!(t::Relax, v::Val{PLUS}, g::DAT, c::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    # out loops makes a temporary sum (minus one argument)
    # a reverse is then compute with respect to this argument
    count = 0
    children_idx = children(g, k)
    for i in children_idx
        is_num(c, i) && continue                     # don't contract a number valued argument
        (count >= MAX_ASSOCIATIVE_REVERSE) && break
        tsum = zero(MC{N,T})
        count += 1
        for j in children_idx
            if j != i
                if is_num(c, j)
                    tsum += num(c, j)
                else
                    tsum += set(c, j)
                end
            end
        end
        _, w, _ = IntervalContractors.plus_rev(interval(c, k), interval(c, i), interval(tsum))
        isempty(w) && (return false)
        c[i] = MC{N,T}(w)
    end
    return true
end

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = x * y` which updates x & y.
"""
function rprop_2!(t::Relax, v::Val{MULT}, g::DAT, c::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}

    is_num(c,k) && (return true)
    x = child(g, 1, k)
    y = child(g, 2, k)

    if xset_ynum(b, x, y)
        c, a, q = IntervalContractors.mul_rev(interval(c, k), interval(c, x), num(c, y))
    elseif xnum_yset(b, x, y)
        c, a, q = IntervalContractors.mul_rev(interval(c, k), num(c, x), interval(c, y))
    else
        c, a, q = IntervalContractors.mul_rev(interval(c, k), interval(c, x), interval(c, y))
    end

    if !is_num(c, x)
        isempty(a) && (return false)
        c[x] = MC{N,T}(a)
    end
    if !is_num(c, x)
        isempty(q) && (return false)
        c[y] = MC{N,T}(q)
    end
    return true
end

"""
$(FUNCTIONNAME)

Updates storage tapes with reverse evalution of node representing `n = *(x,y,z...)` which updates x, y, z and so on.
"""
function rprop_n!(t::Relax, v::Val{MULT}, g::DAT, c::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    # a reverse is then compute with respect to this argument
    count = 0
    children_idx = children(g, k)
    for i in children_idx
        is_num(b, i) && continue                     # don't contract a number valued argument
        (count >= MAX_ASSOCIATIVE_REVERSE) && break
        tmul = one(MC{N,T})
        count += 1
        for j in children_idx
            if i != j
                if is_num(b, j)
                    tmul *= num(b, j)
                else
                    tmul *= set(b, j)
                end
            end
        end
        _, w, _ = IntervalContractors.mul_rev(interval(b, k), interval(b, c), interval(tmul))
        isempty(w) && (return false)
        c[i] = MC{N,T}(w)
    end
    return true
end

for (f, fc, F) in ((-, MINUS, IntervalContractors.minus_rev),
                   (^, POW, IntervalContractors.power_rev),
                   (/, DIV, IntervalContractors.div_rev))
    @eval function rprop!(t::Relax, v::Val{$fc}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        is_num(b,k) && (return true)
        x = child(g, 1, k)
        y = child(g, 2, k)

        if xset_ynum(b, x, y)
            z, u, v = ($F)(interval(b, k), interval(b, x), num(b, y))
        elseif xnum_yset(b, x, y)
            z, u, v = ($F)(interval(b, k), num(b, x), interval(b, y))
        else
            z, u, v = ($F)(interval(b, k), interval(b, x), interval(b, y))
        end
        if !is_num(c, x)
            isempty(a) && (return false)
            b[x] = MC{N,T}(u)
        end
        if !is_num(c, y)
            isempty(b) && (return false)
            b[y] = MC{N,T}(v)
        end
        return true
    end
end

rprop!(t::Relax, v::Val{USER}, g::DAT, b::RelaxCache, k::Int) = true
rprop!(t::Relax, v::Val{USERN}, g::DAT, b::RelaxCache, k::Int) = true

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval rprop!(t::Relax, v::Val{$ft}, g::DAT, b::RelaxCache, k::Int) = true
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using CUDAnative: device!, CuDevice
device!(CuDevice(3))
using CuArrays
CuArrays.allowscalar(false)

# define the target model
using Yao, Yao.ConstGate, BitBasis
using CuYao
using YaoArrayRegister: u1rows!, mulrow!
using QuantumPEPS
using BenchmarkTools

function run_benchmark(nx=6, ny=4; usecuda, nv=2, depth=1)
    model = J1J2(nx, ny; J2=0.5, periodic=false)
    config = QPEPSConfig(ny, nv, nx-nv, depth)

    reg0 = zero_state(nqubits(config); nbatch=1024)
    usecuda && (reg0 = reg0 |> cu)
    qpeps = QPEPSMachine(config, reg0)
    dispatch!(qpeps.runtime.circuit, :random)
    display(@benchmark energy($qpeps, $model))
end

#run_benchmark(;usecuda=false)
#run_benchmark(;usecuda=true)
run_benchmark(6, 6; usecuda=true, nv=1, depth=2)

# REPORT
# Titan-V: 640ms
# CPU: 59s


text:
export SGD, solve,getSGDsolver

"""
Stochastic Gradient Descent
"""
mutable struct SGD{T}
    maxEpochs::Int
    miniBatch::Int
    out::Bool
    learningRate::T
    momentum::T
    nesterov::Bool
	ADAM::Bool
end

function getSGDsolver(TYPE::Type ;maxEpochs=10,miniBatch=16,out=true,learningRate=0.1,momentum=0.9,nesterov=false,ADAM=false)
	if ADAM && nesterov
		warn("sgd(): ADAM and nestrov together - choosing ADAM");
		nesterov  = false;
	end
	return SGD{TYPE}(maxEpochs,miniBatch,out,convert(TYPE,learningRate),convert(TYPE,momentum), nesterov, ADAM)
end

Base.display(this::SGD)=println("SGD(maxEpochs=$(this.maxEpochs),miniBatch=$(this.miniBatch),learningRate=$(this.learningRate),momentum=$(this.momentum),nesterov=$(this.nesterov),ADAM=$(this.ADAM))")

function solve(this::SGD{T},objFun::dnnObjFctn,xc::Array{T},Y::Array{T},C::Array{T},Yv::Array{T},Cv::Array{T}) where {T<:Number}

    # evaluate training and validation
    epoch = 1;
    xOld = copy(xc);
    dJ = zeros(T,size(xc));
    mJ = zeros(T,size(xc));
    vJ = zeros(T,size(xc));
    if this.ADAM
        mJ = zeros(T,size(xc));
        vJ = zeros(T,size(xc));
    end
    beta2 = convert(T,0.999);
    beta1 = this.momentum;

    lr    = this.learningRate

    if this.out; display(this); end

    # Declare tmp - We know nothing about its shape or datatypes
    tmp = Array{Any}(0,0)

    while epoch <= this.maxEpochs
        nex = size(Y,2)
        ids = randperm(nex)

        for k=1:ceil(Int64,nex/this.miniBatch)
            idk = ids[(k-1)*this.miniBatch+1: min(k*this.miniBatch,nex)]
            if this.nesterov && !this.ADAM
                Jk,dummy,dJk,tmp = evalObjFctn(objFun,xc-this.momentum*dJ,Y[:,idk],C[:,idk],tmp);
            else
                Jk,dummy,dJk,tmp = evalObjFctn(objFun,xc,Y[:,idk],C[:,idk],tmp);
            end

            if this.ADAM
               mJ = beta1*mJ + (one(T)-beta1)*(dJk)
               vJ = beta2*vJ + (one(T)-beta2)*(dJk.^2)
#              dJ = lr*((mJ./(1-beta1^(epoch)))./sqrt.((vJ./(1-beta2^(epoch)))+1e-8))
            else
               dJ = lr*dJk + this.momentum*dJ
            end
            xc = xc - dJ
        end
        # we sample 2^12 images from the training set for displaying the objective.
        idt     = ids[1:min(nex,2^12)]
        Jtrain,ptrain   = getMisfit(objFun,xc,Y[:,idt],C[:,idt],tmp,false);
        Jval,pVal = getMisfit(objFun,xc,Yv,Cv,tmp,false);

        if this.out;
            @printf "%d\t%1.2e\t%1.2f\t%1.2e\t%1.2e\t%1.2f\n" epoch Jtrain 100*(1-ptrain[3]/ptrain[2]) norm(xOld-xc) Jval 100*(1-pVal[3]/pVal[2])
        end

        xOld       = copy(xc);
        epoch = epoch + 1;
    end
    return xc
end


text:
"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 6/28/2017, Last Modified: 6/28/2017 \n
-------------------------------------------------------------------------------------\n
"""
macro controlVariable(n,name,num)
  n=esc(n)
  return quote
      $name=$n.r.ocp.u[:,num]
    end
end

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 5/31/2017, Last Modified: 6/14/2017 \n
-------------------------------------------------------------------------------------\n
"""
macro DiffEq(n,expr)
  n=esc(n)
  return quote
      $n.stateEquations=$(Expr(:quote,expr))
    end
end

function dx_EQ{T<:Any}(n::NLOpt,x::Array{T,2},u::Array{T,2})
  dx_expr=parse_DiffEq(n,x,u)
  for st in 1:n.ocp.state.num
    for j in 1:n.ocp.state.pts
      dx[j,st]=create_DiffEq(n,dx_expr[j,st])
    end
  end
  return dx
end

"""
#TODO consider appending a zero so it is not just a jump variable
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 6/11/2017, Last Modified: 12/06/2019 \n
-------------------------------------------------------------------------------------\n
"""

function parse_DiffEq{T<:Any}(n::NLOpt,x::Array{T,2},u::Array{T,2})
  if n.s.ocp.integrationMethod==:tm; L=size(x)[1]; else; L=size(x)[1]-1; end

  dx=Array{Expr}(L,n.ocp.state.num);
  for qq in 1:L
    dx[qq,:]=deepcopy(n.ocp.DXexpr)
  end

  #TODO use the statenames here! -> currently running too often!
  state=Array{Symbol}(undef, n.ocp.state.num,1)
  for st in 1:n.ocp.state.num
    state[st]=Symbol(string("x",:($(st))))
  end
  control=Array{Symbol}(undef, n.ocp.control.num,1)
  for ctr in 1:n.ocp.control.num
    control[ctr]=Symbol(string("u",:($(ctr))))
  end

  for st_high in 1:n.ocp.state.num
    for qq in 1:L #TODO this will not start at 1 for the :ps method...unless you go bac to x and u
      # states
      for st in 1:n.ocp.state.num
        dx[qq,st_high]=rePlace(dx[qq,st_high],state[st],:($(n.r.ocp.x[qq,st])))
      end

      # controls
      for ctr in 1:n.ocp.control.num
        dx[qq,st_high]=rePlace(dx[qq,st_high],control[ctr],:($(n.r.ocp.u[qq,ctr])))
      end
    end
  end
  return dx
end
#dx=parse_DiffEq(n,n.r.ocp.x,n.r.ocp.u)

"""
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 6/12/2017, Last Modified: 6/12/2017 \n
Citations: \n
----------\n
Based off roughly of:
Source: [located here](https://github.com/JuliaLang/julia/blob/master/base/cartesian.jl#L278-L353)
-------------------------------------------------------------------------------------\n
"""
type Replace
    sym::Symbol
    val::JuMP.Variable
end

function rePlaceEX(ex::Expr,r::Replace)
  for i in 1:length(ex.args)
    if ex.args[i]==r.sym
      ex.args[i]=r.val
    end
  end
  return ex
end

function dig(ex)
  branch=falses(length(ex.args));
  for i in 1:length(ex.args)
    if typeof(ex.args[i])==Expr;
      branch[i]=true;
    end
  end
  return branch
end

function rePlace(ex::Expr,sym::Symbol,val::JuMP.Variable)
  r=Replace(sym,val)
  ex_tmp=[deepcopy(ex)]
  unwrap=true
  level=1
  while unwrap
    expr_arr=[];startCT=true;
    for j in 1:length(ex_tmp)
      ex_tmp[j]=rePlaceEX(ex_tmp[j],r);   # replace any on this branch
      branch=dig(ex_tmp[j]);
      if any(branch) && startCT
        expr_arr=ex_tmp[j].args[branch];
        tree=[any(branch)];
        startCT=false;
      elseif any(branch)
        expr_arr=[expr_arr,ex_tmp[j].args[branch]];
        tree=[tree,any(branch)];
      end
    end

    if level==1; ex=ex_tmp; end

    if isempty(expr_arr)
      unwrap=false
    else
      ex_tmp=expr_arr
      level+=1;
      if level>100
        error("\n Are the variables in the expression nested > 100 levels deep? If not report a bug, if so, consider increasing this limit \n ")
      end
    end
  end
  return ex[1]
end

#= TODO make test function
using JuMP
m=Model()
@variable(m,x[1:2,1:4])
expr=:($:x3*3*sin($:x3)*cos(cos($:x3)))
expr=rePlace(expr,:x3,:($(x[2,3])))
=#


text:
using NeRCA
using Test

# make_cherenkov_calculator()
# ccalc = NeRCA.make_cherenkov_calculator(NeRCA.Track([0,0,-1], [0,0,1], 0))
# @test 3.33564095 ≈ ccalc([0, 0, 0])


text:
function listFunctions(_module::Module)
    _names = names(_module; all=true)

    _functions = Function[]
    for _name in _names
        _function = getfield(_module, _name)
        if (
            typeof(_function) <: Function
         && !occursin("#",string(_name))
         && !(string(_function) in ["eval","include"])
         )
            push!(_functions, _function)
        end
    end

    _functions
end

function listMethods(_module::Module)
    _functions = listFunctions(_module)
    _methods = Method[]
    for _f in _functions
        _methodList = methods(_f)
        push!(_methods,_methodList.ms...)
    end
    _methods
end

function createDefaultArguments(method)
    @info "Create default arguments for method[$(method.name)]"
    sig = try
            collect(method.sig.types)
        catch e
            @warn e
            return nothing
        end
    popfirst!(sig)

    argsAsVector = []
    for _type in sig
        arg = try
           createDefaultObject(getNonMissingTypeOfUnionType(_type))
        catch e
            @warn e
            return nothing
        end
        push!(argsAsVector,
              arg)
    end
    argsAsVector
end

function createDefaultObject(_type::DataType)

    if _type <: Number
        return _type(0)
    end
    if _type <: String
        return "lorem ipsum"
    end
    if _type <: Date
        return today()
    end
    if _type <: Enum
        return first(instances(_type))
    end

    return _type()
end

function createDefaultObject(::Type{Dict})
    return Dict()
end

function invokeMethod(_function::Function, args::Vector, procID::Int64)

    argsTypesForPrinting = join(string.(typeof.(args)),", ")
    @info "# Invoke $(_function)($argsTypesForPrinting) on procID[$procID]"

    try
        future = @spawnat procID _function(args...)
    catch e
        # try-catch on the @spawnat is not enough because Exceptions on remote
        #   computations are captured and rethrown locally. Therefore the calling
        #   method needs to try-catch the call to this function.
        # see https://docs.julialang.org/en/v1/stdlib/Distributed/index.html#Distributed.RemoteException
        error(e)
    end
end

function invokeMethod(_method::Method, args::Vector, procID::Int64)
    fct = getfield(_method.module, _method.name)
    invokeMethod(fct, args, procID)
end

function invokeMethodOnAllProcs(_method::Method)
    args = createDefaultArguments(_method)
    if isnothing(args)
        @warn "Unable to create default arugments => Skip method[$(_method)]"
        return
    end
    for procID in 1:nprocs()
        try
            invokeMethod(_method, args, procID)
        # try-catch on the @spawnat is not enough because Exceptions on remote
        #   computations are captured and rethrown locally
        # see https://docs.julialang.org/en/v1/stdlib/Distributed/index.html#Distributed.RemoteException
        catch e
            @debug "RemoteException was caught so that we can carry on."
        end
    end
end

function invokeMethodsOfModule(_module::Module
                              ;excludeMethods::Vector{Symbol} = Symbol[])

    methods = listMethods(_module)
    filter!(x-> !(x.name in excludeMethods),methods)

    for m in methods
        invokeMethodOnAllProcs(m)
    end
end

# This method exists so that we can call it even without having to test if the
#   argument is a Union
function getNonMissingTypeOfUnionType(arg::Any)
    return arg
end

# NOTE: We cannot say that this method returns a DataType (::DataType) because
#         arrays are not datatypes
function getNonMissingTypeOfUnionType(arg::Union)
    if arg.a != Missing
        return arg.a
    else
        return arg.b
    end
end


text:
using Test

const repository_root = dirname(dirname(@__DIR__))

for filename in ARGS
    path = joinpath(repository_root, filename)
    @info "Starting $(filename)"
    try
        @testset "$(filename)" begin
            include(path)
        end
    catch ex
        @error "" exception=(ex, catch_backtrace())
    end
    @info "Finished $(filename)"
end

max_stars_repo_name: chengchingwen/julia

max_issues_repo_name: chengchingwen/julia

max_forks_repo_name: chengchingwen/julia


text:
struct _Edge
    vertex::Int
    weight::Int
end
_Edge() = _Edge(-1, -1)
_getweight(e::_Edge) = e.weight
_getvertex(e::_Edge) = e.vertex

mutable struct AdjacencyList <: AbstractGraph
    cntedge::Int
    cntvertex::Int
    vertex::Vector # Should eventually be of type Vertex
    mark::Vector # 0...unvisited; 1...visited
    function AdjacencyList(n::Int)
        cntedge = 0
        cntvertex = n
        vertex = Vector{LinkedList{_Edge}}(undef, cntvertex)
        for i = 1:cntvertex
            vertex[i] = LinkedList{_Edge}()
        end
        mark = zeros(Int, cntvertex)
        new(cntedge, cntvertex, vertex, mark)
    end
end

countedges(g::AdjacencyList) = g.cntedge
countvertices(g::AdjacencyList) = g.cntvertex

function first(g::AdjacencyList, v::Vertex)
    isempty(g.vertex[v]) && return g.cntvertex + 1
    movetostart!(g.vertex[v])
    edge = getvalue(g.vertex[v])
    return _getvertex(edge)
end

function next(g::AdjacencyList, v::Vertex, w::Vertex)
    if isedge(g, v, w)
        if currentposition(g.vertex[v]) < length(g.vertex[v])
            next!(g.vertex[v])
            edge = getvalue(g.vertex[v])
            return _getvertex(edge)
        end
    end
    return g.cntvertex + 1
end

function setedge!(g::AdjacencyList, v::Vertex, w::Vertex, weight)
    (weight > 0) || throw(ArgumentError("Illegal weight value."))
    edge = _Edge(w, weight)
    if isedge(g, v, w)
        remove!(g.vertex[v])
        insert!(g.vertex[v], edge)
    else
        g.cntedge += 1
        movetostart!(g.vertex[v])
        while currentposition(g.vertex[v]) <= length(g.vertex[v])
            currentedge = getvalue(g.vertex[v])
            (_getvertex(currentedge) > w) && break
            next!(g.vertex[v])
        end
        insert!(g.vertex[v], edge)
    end
end

function remove!(g::AdjacencyList, v::Vertex, w::Vertex)
    if isedge(g, v, w)
        remove!(g.vertex[v])
        g.cntedge -= 1
    end
end

function isedge(g::AdjacencyList, v::Vertex, w::Vertex)
    movetostart!(g.vertex[v])
    while currentposition(g.vertex[v]) <= length(g.vertex[v])
        edge = getvalue(g.vertex[v])
        if _getvertex(edge) == w
            return true
        end
        next!(g.vertex[v])
    end
    return false
end

function getweight(g::AdjacencyList, v::Vertex, w::Vertex)
    if isedge(g, v, w)
        edge = getvalue(g.vertex[v])
        return _getweight(edge)
    else
        return 0
    end
end

getmark(g::AdjacencyList, v::Vertex) = g.mark[v]
setmark!(g::AdjacencyList, v::Vertex, value) = (g.mark[v] = value)


text:
function FlowDirection(wn_dict::Dict{Any,Any}, Qmin::Dict{String, Float64}, Qmax::Dict{String,Float64})
    PositiveFlowLinks = Array{String,1}(0)
    NegativeFlowLinks = Array{String,1}(0)
    ReversibleFlowLinks = Array{String,1}(0)

    for link in wn_dict["links"]
        name = link["name"]
        if Qmin[name] > 1e-9
            push!(PositiveFlowLinks, name)
        elseif Qmax[name] <-1e-9
            push!(NegativeFlowLinks, name)
        else
            push!(ReversibleFlowLinks, name)
        end
    end
    return PositiveFlowLinks, NegativeFlowLinks, ReversibleFlowLinks
end


text:
function sendtosimple(p::Int, nm, val)
  ref = @spawnat(p, eval(Main, Expr(:(=), nm, val)))
end

macro sendto(p, nm, val)
  return :( sendtosimple($p, $nm, $val) )
end

macro broadcast(nm, val)
  quote
  @sync for p in workers()
    @async sendtosimple(p, $nm, $val)
  end
  end
end

# same as np.allclose
@everywhere function allclose(a::Any, b::Any; rtol::Float64 = 1e-5, atol::Float64 = 1e-8)
  all(abs(a - b) .<= (atol + rtol * abs(b)))
end

function separateAxis2D(A::Array{Float64,2})
  x = Array(Float64, size(A,1))
	y = Array(Float64, size(A,1))
	for j=1:size(A,1)
	  x[j] = A[j,1]
	  y[j] = A[j,2]
	end
	x,y
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function get_all_png_files_alpha(folder)
    imagelist = filter(f -> !startswith(f, ".") && endswith(f, "png"), readdir(folder))
    return map(f -> string(folder, f), imagelist)
end

function alphaimages(fname)
    imagelist = get_all_png_files_alpha(dirname(@__FILE__) * "/../docs/src/assets/figures/")
    width, height = 3000, 3000
    Drawing(width, height, fname)
    origin()
    background("grey50")
    for imgfile in imagelist
        img = readpng(imgfile)
        rand(Bool) ? placeimage(img, rand(-width/2:width/2), rand(-height/2:height/2), rand(0.25:0.1:0.75)) : placeimage(img, Point(rand(-width/2:width/2), rand(-height/2:height/2)), rand(0.25:0.1:0.75))
    end
    @test finish() == true
end

fname = "images-with-alpha.pdf"
alphaimages(fname)
println("...finished test: output in $(fname)")


text:
import Base: *, +, -, /

function *(d1::Dict, d2::Dict)
    """
    entrywise multiplication. Will error if any keys are present in d1 and not in d2
    """
    d3 = Dict(p => 0 for p in keys(d1))
    for p in keys(d1)
        d3[p] = d1[p]*d2[p]
    end
    return d3
end

function +(d1::Dict, d2::Dict)
    """
    entrywise addition. Will error if any keys are present in d1 and not in d2
    """
    d3 = Dict(p => 0 for p in keys(d1))
    for p in keys(d1)
        d3[p] = d1[p]+d2[p]
    end
    return d3
end

function -(d1::Dict, d2::Dict)
    """
    entrywise subtraction. Will error if any keys are present in d1 and not in d2
    """
    d3 = Dict(p => 0 for p in keys(d1))
    for p in keys(d1)
        d3[p] = d1[p]-d2[p]
    end
    return d3
end

function /(d1::Dict, d2::Dict)
    """
    entrywise division. Will error if any keys are present in d1 and not in d2
    """
    d3 = Dict(p => 0 for p in keys(d1))
    for p in keys(d1)
        d3[p] = d1[p]/d2[p]
    end
    return d3
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
function constraint_hash(constraint::BasicConstraint)
    return hash([nameof(constraint.fct), constraint.indices])
end

function constraint_hash(constraint::LinearConstraint)
    return hash([nameof(constraint.fct), constraint.indices, constraint.coeffs, constraint.operator, constraint.rhs])
end

function constraint_hash(constraint::SingleVariableConstraint)
    return hash([nameof(constraint.fct), constraint.rhs, constraint.lhs])
end


text:
using MultivariateStats
using LinearAlgebra, StatsBase, SparseArrays
using Test
import Statistics: mean, cov
import Random

@testset "Whitening" begin

    Random.seed!(34568)

    ## data

    d = 3
    n = 5

    X = rand(d, n)
    mval = vec(mean(X, dims=2))
    C = cov(X, dims=2)
    C0 = copy(C)

    emax = maximum(eigvals(C))
    rc = 0.01
    Cr = C + (emax * rc) * Matrix(I, d, d)

    ## cov_whitening

    cf = cholesky(Hermitian(C, :U))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    cf = cholesky(Hermitian(C, :L))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, 0.0)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, rc)
    @test istriu(W)
    @test W'Cr * W ≈ Matrix(I, d, d)

    # Whitening from data

    f = fit(Whitening, X)
    W = f.W
    @test isa(f, Whitening{Float64})
    @test mean(f) === f.mean
    @test length(f) == d
    @test size(f) == (d,d)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)
    @test transform(f, X) ≈ W' * (X .- f.mean)

    f = fit(Whitening, X; regcoef=rc)
    W = f.W
    @test W'Cr * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=mval)
    W = f.W
    @test W'C * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=0)
    Cx = (X * X') / (n - 1)
    W = f.W
    @test W'Cx * W ≈ Matrix(I, d, d)

    # invsqrtm

    R = invsqrtm(C)
    @test C == C0
    @test R ≈ inv(sqrt(C))

    # mixing types
    X = rand(Float64, 5, 10)
    XX = convert.(Float32, X)

    M = fit(Whitening, X)
    MM = fit(Whitening, XX)

    # mixing types should not error
    transform(M, XX)
    transform(MM, X)

    # type consistency
    @test eltype(mean(M)) == Float64
    @test eltype(mean(MM)) == Float32

    # sparse arrays
    SX = sprand(Float32, d, n, 0.75)
    SM = fit(Whitening, SX; mean=sprand(Float32, 3, 0.75))
    Y = transform(SM, SX)
    @test eltype(Y) == Float32

    # different dimensions
    @test_throws DomainError fit(Whitening, X'; dims=3)
    M1 = fit(Whitening, X'; dims=1)
    M2 = fit(Whitening, X; dims=2)
    @test M1.W == M2.W
    @test_throws DimensionMismatch transform(M1, rand(6,4))
    @test_throws DimensionMismatch transform(M2, rand(4,6))
    Y1 = transform(M1,X')
    Y2 = transform(M2,X)
    @test Y1' == Y2
    @test_throws DimensionMismatch transform(M1, rand(7))
    V1 = transform(M1,X[:,1])
    V2 = transform(M2,X[:,1])
    @test V1 == V2
end


text:
"""
    autocovariance(v::Vector,h::Int; corrected::Bool=true)
``\\hat{\\gamma}(h) =\\frac{1}{n}\\sum_{t=1}^{n-h}(v_{t+h}-\\bar{v})(v_t-\\bar{v})^*``
Calculate the autocovariance of dataset `v` with a delay `h`. If `corrected` is `true`
(the default) then the sum is scaled with `n-h`, whereas the sum
is scaled with `n` if corrected is `false` where `n = length(v)`.
"""
function autocovariance(v::AbstractVector,h::Int; corrected::Bool=true, mean = mean(v))
    n = length(v)
    covsum = zero(mean)
    for i in 1:n-h
        @inbounds covsum += (v[i]-mean)*conj(v[i+h]-mean)
    end
    gamma = covsum/(n - Int(corrected))
    return gamma
end

@doc raw"""
    pseudo_cov(x, y; xmean = mean(x), ymean = mean(y), corrected = true)
Compute the pseudo covariance between collections `x` and `y` returning a scalar:
```math
\frac{1}{n}\sum_{i=1}^{n} (x_i - \bar{x})(y_{i} - \bar{y})
```
Optionally,
precomputed means can be passed as keyword arguments. `pseudo_cov(x,y)` is functionally
equivalent to `Statistics.cov(x, conj(y); corrected = false)` but it is found to be
significantly faster and avoids allocations.
"""
@inline function pseudo_cov(x, y; xmean = mean(x), ymean = mean(y), corrected = true)
    n = length(x)
    @assert length(y) == n
    res = zero(promote_type(eltype(x), eltype(y))) / 1
    for i = 1:n
        @inbounds res += (x[i] - xmean) * (y[i] - ymean)
    end
    return res / (n - Int(corrected))
end


text:
@doc raw"""
    UnbiasedSKCE(k)

Unbiased estimator of the squared kernel calibration error (SKCE) with kernel `k`.

Kernel `k` on the product space of predictions and targets has to be a `Kernel` from the
Julia package
[KernelFunctions.jl](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl)
that can be evaluated for inputs that are tuples of predictions and targets.

# Details

The estimator is unbiased and not guaranteed to be non-negative. Its sample complexity
is ``O(n^2)``, where ``n`` is the total number of samples.

Let ``(P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and corresponding
targets. The estimator is defined as
```math
\frac{2}{n(n-1)} \sum_{1 \leq i < j \leq n} h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big),
```
where
```math
\begin{aligned}
h_k\big((μ, y), (μ', y')\big) ={}&   k\big((μ, y), (μ', y')\big)
                                   - 𝔼_{Z ∼ μ} k\big((μ, Z), (μ', y')\big) \\
                                 & - 𝔼_{Z' ∼ μ'} k\big((μ, y), (μ', Z')\big)
                                   + 𝔼_{Z ∼ μ, Z' ∼ μ'} k\big((μ, Z), (μ', Z')\big).
\end{aligned}
```

# References

Widmann, D., Lindsten, F., & Zachariah, D. (2019). [Calibration tests in multi-class
classification: A unifying framework](https://proceedings.neurips.cc/paper/2019/hash/1c336b8080f82bcc2cd2499b4c57261d-Abstract.html).
In: Advances in Neural Information Processing Systems (NeurIPS 2019) (pp. 12257–12267).

Widmann, D., Lindsten, F., & Zachariah, D. (2021). [Calibration tests beyond
classification](https://openreview.net/forum?id=-bxf89v3Nx).

See also: [`BiasedSKCE`](@ref), [`BlockUnbiasedSKCE`](@ref)
"""
struct UnbiasedSKCE{K<:Kernel} <: SKCE
    """Kernel of estimator."""
    kernel::K
end

function (skce::UnbiasedSKCE)(predictions::AbstractVector, targets::AbstractVector)
    @unpack kernel = skce

    # obtain number of samples
    nsamples = check_nsamples(predictions, targets, 2)

    @inbounds begin
        # evaluate the kernel function for the first pair of samples
        hij = unsafe_skce_eval(
            kernel, predictions[1], targets[1], predictions[2], targets[2]
        )

        # initialize the estimate
        estimate = hij / 1

        # for all other pairs of samples
        n = 1
        for j in 3:nsamples
            predictionj = predictions[j]
            targetj = targets[j]

            for i in 1:(j - 1)
                predictioni = predictions[i]
                targeti = targets[i]

                # evaluate the kernel function
                hij = unsafe_skce_eval(kernel, predictioni, targeti, predictionj, targetj)

                # update the estimate
                n += 1
                estimate += (hij - estimate) / n
            end
        end
    end

    return estimate
end

@doc raw"""
    BlockUnbiasedSKCE(k[, blocksize = 2])

Unbiased estimator of the squared kernel calibration error (SKCE) with kernel `k` that
uses blocks with `blocksize` samples.

Kernel `k` on the product space of predictions and targets has to be a `Kernel` from the
Julia package
[KernelFunctions.jl](https://github.com/JuliaGaussianProcesses/KernelFunctions.jl)
that can be evaluated for inputs that are tuples of predictions and targets.

The number of samples per block must be at least 2 and at most the total number of samples.
Note that samples in the last block are discarded if it is incomplete (see details below).

# Details

The estimator is unbiased and not guaranteed to be non-negative. Its sample complexity
is ``O(Bn)``, where ``B`` is the block size and ``n`` is the total number of samples.

Let ``(P_{X_i}, Y_i)_{i=1,\ldots,n}`` be a data set of predictions and corresponding
targets. The estimator with block size ``B`` is defined as
```math
{\bigg\lfloor \frac{n}{B} \bigg\rfloor}^{-1} \sum_{b=1}^{\lfloor n/B \rfloor}
\frac{2}{B(B-1)} \sum_{(b - 1) B < i < j \leq bB} h_k\big((P_{X_i}, Y_i), (P_{X_j}, Y_j)\big),
```
where
```math
\begin{aligned}
h_k\big((μ, y), (μ', y')\big) ={}&   k\big((μ, y), (μ', y')\big)
                                   - 𝔼_{Z ∼ μ} k\big((μ, Z), (μ', y')\big) \\
                                 & - 𝔼_{Z' ∼ μ'} k\big((μ, y), (μ', Z')\big)
                                   + 𝔼_{Z ∼ μ, Z' ∼ μ'} k\big((μ, Z), (μ', Z')\big).
\end{aligned}
```

# References

Widmann, D., Lindsten, F., & Zachariah, D. (2019). [Calibration tests in multi-class
classification: A unifying framework](https://proceedings.neurips.cc/paper/2019/hash/1c336b8080f82bcc2cd2499b4c57261d-Abstract.html).
In: Advances in Neural Information Processing Systems (NeurIPS 2019) (pp. 12257–12267).

Widmann, D., Lindsten, F., & Zachariah, D. (2021). [Calibration tests beyond
classification](https://openreview.net/forum?id=-bxf89v3Nx).

See also: [`BiasedSKCE`](@ref), [`UnbiasedSKCE`](@ref)
"""
struct BlockUnbiasedSKCE{K<:Kernel} <: SKCE
    """Kernel of estimator."""
    kernel::K
    """Number of samples per block."""
    blocksize::Int

    function BlockUnbiasedSKCE{K}(kernel::K, blocksize::Int) where {K}
        blocksize ≥ 2 || error("there must be at least two samples per block")
        return new{K}(kernel, blocksize)
    end
end

function BlockUnbiasedSKCE(kernel::Kernel, blocksize::Int=2)
    return BlockUnbiasedSKCE{typeof(kernel)}(kernel, blocksize)
end

function (skce::BlockUnbiasedSKCE)(predictions::AbstractVector, targets::AbstractVector)
    @unpack kernel, blocksize = skce

    # obtain number of samples
    nsamples = check_nsamples(predictions, targets, blocksize)

    # compute number of blocks
    nblocks = nsamples ÷ blocksize

    # create iterator of partitions
    blocks = Iterators.take(
        zip(
            Iterators.partition(predictions, blocksize),
            Iterators.partition(targets, blocksize),
        ),
        nblocks,
    )

    # compute average u-statistic
    unbiased = UnbiasedSKCE(kernel)
    estimate = mean(unbiased(_predictions, _targets) for (_predictions, _targets) in blocks)

    return estimate
end


text:
# Numerical operations
1 + 1
2 * 2
2 / 4
4 - 1

dump(3)

dump(9223372036854775807+1)

dump(3.4)

dump(44 // 23)

dump("Helllo")

#d
"Hallo, " * " Julia" ^ 3

false

true

true

!true

!false

a = 3
dump(a)

a = 3.4
dump(a)

# Type assertion

a = 3.0
b = 4

c = a + b
c = a % b
dump(c)

staslijads = "Hi this is a string \n Test"

staslijads

multi = """
Hi
thi
is a
multi line
string
"""

print(staslijads)

println("Kickarz")
println("Nextrow")

grades = [3.0, 1.3, 4.0, 1.0, 1.0]
grades2 = [3.0, 2.0, 5.0, 1.3, 1.0]

grades[4]

using(Statistics)
mean(grades)
var(grades)

std(grades)

cor(grades, grades2)

friend_count = 10:5:200
length(friend_count)
friend_count[39]

large_range = 1:1000000000000

large_range * 2
large_range[1000000]

### do a mini simulation
using Statistics
function simulate(agent_count::Integer, steps::Integer)

    agents = [rand() for i in 1:agent_count]

    for step in 1:steps
        #do simulation step
        agents[1] = agents[1] + 1
    end
    mean(agents)
end

println("---------------------")
for agent_count in 1:5:100
    println(simulate(agent_count, 8))
end


text:
#__precompile__(false)

module Gen

const generated_functions = []

"""
    load_generated_functions(__module__=Main)

Permit use of generative functions written in the static modeling language up
to this point. Functions are loaded into Main by default.
"""
function load_generated_functions(__module__::Module=Main)
    for function_defn in generated_functions
        Core.eval(__module__, function_defn)
    end
end

"""
    @load_generated_functions

Permit use of generative functions written in the static modeling language up
to this point. Functions are loaded into the calling module.
"""
macro load_generated_functions()
    for function_defn in generated_functions
        Core.eval(__module__, function_defn)
    end
end

export load_generated_functions, @load_generated_functions

# built-in extensions to the reverse mode AD
include("backprop.jl")

# addresses and address selections
include("address.jl")

# abstract and built-in concrete choice map data types
include("choice_map.jl")

# a homogeneous trie data type (not for use as choice map)
include("trie.jl")

# generative function interface
include("gen_fn_interface.jl")

# built-in data types for arg-diff and ret-diff values
include("diff.jl")

# built-in probability disributions
include("modeling_library/modeling_library.jl")

# optimization of trainable parameters
include("optimization.jl")

# dynamic embedded generative function
include("dynamic/dynamic.jl")

# static IR generative function
include("static_ir/static_ir.jl")

# optimization for built-in generative functions (dynamic and static IR)
include("builtin_optimization.jl")

# DSLs for defining dynamic embedded and static IR generative functions
# 'Dynamic DSL' and 'Static DSL'
include("dsl/dsl.jl")

# inference and learning library
include("inference/inference.jl")

end # module Gen


text:
@testset "experiment.jl" begin

end


text:
using Test
using Example
using BigProject

@test BigProject.f() == 1

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using Primes
using Test
using DataStructures: SortedDict

import Primes: isprime, primes, primesmask, factor, ismersenneprime, isrieselprime, Factorization

@test primes(10000) == primes(2, 10000) == [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
    157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
    239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
    331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
    421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503,
    509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
    613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701,
    709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
    821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911,
    919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
    1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091,
    1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
    1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
    1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361,
    1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
    1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531,
    1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609,
    1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699,
    1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789,
    1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
    1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
    1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083,
    2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161,
    2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273,
    2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
    2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441,
    2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551,
    2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663,
    2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729,
    2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
    2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917,
    2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023,
    3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137,
    3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251,
    3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
    3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,
    3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533,
    3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617,
    3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
    3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
    3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917,
    3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013,
    4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111,
    4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219,
    4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
    4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423,
    4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519,
    4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639,
    4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729,
    4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
    4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951,
    4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023,
    5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147,
    5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261,
    5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
    5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,
    5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563,
    5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659,
    5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779,
    5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
    5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981,
    5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089,
    6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199,
    6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287,
    6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
    6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491,
    6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607,
    6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709,
    6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827,
    6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
    6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013,
    7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129,
    7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243,
    7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369,
    7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
    7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577,
    7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681,
    7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
    7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901,
    7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017,
    8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123,
    8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237,
    8243, 8263, 8269, 8273, 8287, 8291, 8293, 8297, 8311, 8317, 8329, 8353,
    8363, 8369, 8377, 8387, 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461,
    8467, 8501, 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597,
    8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, 8681, 8689,
    8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, 8747, 8753, 8761, 8779,
    8783, 8803, 8807, 8819, 8821, 8831, 8837, 8839, 8849, 8861, 8863, 8867,
    8887, 8893, 8923, 8929, 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001,
    9007, 9011, 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109,
    9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, 9203, 9209,
    9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, 9293, 9311, 9319, 9323,
    9337, 9341, 9343, 9349, 9371, 9377, 9391, 9397, 9403, 9413, 9419, 9421,
    9431, 9433, 9437, 9439, 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511,
    9521, 9533, 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631,
    9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, 9739, 9743,
    9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, 9817, 9829, 9833, 9839,
    9851, 9857, 9859, 9871, 9883, 9887, 9901, 9907, 9923, 9929, 9931, 9941,
    9949, 9967, 9973 ]

for n = 100:100:1000
    @test primes(n, 10n) == primes(10n)[(length(primes(n)) + 1):end]
    @test primesmask(n, 10n) == primesmask(10n)[n:end]
end

for T in [Int, BigInt], n = [1:1000;1000000]
    n = convert(T, n)
    f = factor(n)
    @test n == prod(T[p^k for (p, k) = f])
    prime = n != 1 && length(f) == 1 && get(f, n, 0) == 1
    @test isprime(n) == prime

    s = primesmask(n)
    for k = 1:n
        @test s[k] == isprime(k)
        @test s[k] == primesmask(k, k)[1]
    end
end

# Issue 25
@test primes(-2, 7) == primes(0, 7) == primes(2, 7)

@test !isprime(1000000003)
@test !isprime(1000000005)
@test  isprime(1000000007)
@test  isprime(1000000009)
@test !isprime(1000000011)
@test !isprime(1000000013)

@test !isprime(10000000015)
@test !isprime(10000000017)
@test  isprime(10000000019)
@test !isprime(10000000021)
@test !isprime(10000000023)

@test !isprime(9223372036854775779)
@test !isprime(9223372036854775781)
@test  isprime(9223372036854775783)
@test !isprime(9223372036854775785)
@test !isprime(9223372036854775787)

@test !isprime(0xffffffffffffffc1)
@test !isprime(0xffffffffffffffc3)
@test  isprime(0xffffffffffffffc5)
@test !isprime(0xffffffffffffffc7)
@test !isprime(0xffffffffffffffc9)

for T in [Int8, UInt8, Int16, UInt16, Int128, UInt128]
    @test isprime(T(2))
    @test !isprime(T(4))
end

# Base issue #5210
@test prod([k^v for (k, v) in factor(typemax(UInt32))]) == typemax(UInt32)
@test prod([k^v for (k, v) in factor(typemax(Int8))]) == typemax(Int8)

# factorization of factors > 2^16
@test factor((big(2)^31 - 1)^2) == Dict(big(2^31 - 1) => 2)
@test factor((big(2)^31 - 1) * (big(2)^17 - 1)) == Dict(big(2^31 - 1) => 1, big(2^17 - 1) => 1)

# fast factorization of Int128 (Base issue #11477)
@test factor((Int128(2)^39 - 7)^2) == Dict(Int128(2)^39 - 7 => 2)

# Base issue #9611
@test factor(Int128(2)^101 + 1) == Dict(3 => 1, 845100400152152934331135470251 => 1)

# test second branch, after all small primes in list have been searched
@test factor(10009 * Int128(1000000000000037)) == Dict(10009 => 1, 1000000000000037 => 1)

for n = 1:100
    m = 1
    for (p,k) in factor(n)
        m *= p^k
    end
    @test n == m
end

let i = rand(1:2^(3 * min(Sys.WORD_SIZE,64) ÷ 4))
    @test primes(i, i + 300) == filter(isprime, i:(i + 300)) == filter(isprime, big(i):big(i + 300))
end

@test isprime(BigInt(1000000007))
@test isprime(BigInt(1000000007), 1)
@test isprime(BigInt(10000000019))
@test isprime(parse(BigInt,"359334085968622831041960188598043661065388726959079837"))
@test !isprime(BigInt(1))
@test !isprime(BigInt(10000000020))

# project euler 3: 6857
@test maximum(keys(factor(600851475143))) == 6857

# project euler 7: 104743
euler7(n) = primes(floor(Int, n * log(n * log(n))))[n]
@test euler7(10001) == 104743

# project euler 10: 142913828922
@test sum(map(Int64,primes(2000000))) == 142913828922

# factor(Vector, n)
for V in (Vector, Vector{Int}, Vector{Int128})
    @test factor(V, 1) == Int[]
    @test factor(V, 3) == [3]
    @test factor(V, 4) == [2, 2]
end

# factor with non-default associative containers
@test factor(SortedDict, 100) == factor(Dict, 100) == factor(100)

# factor sets
@test factor(Set, 100) == Set([2, 5])
@test factor(BitSet, 100) == BitSet([2, 5])

# factor other things and fail
@test_throws MethodError factor(Int, 10)
@test_throws MethodError factor(Any, 10)
@test_throws MethodError factor(Tuple, 10)

# factor non-positive numbers:
@test factor(0)  == Dict( 0 => 1)
@test factor(-1) == Dict(-1 => 1)
@test factor(-9) == Dict(-1 => 1, 3 => 2)

@test factor(typemin(Int32))  == Dict(-1 => 1, 2 => 31)
@test factor(typemin(Int64))  == Dict(-1 => 1, 2 => 63)
@test factor(typemin(Int128)) == Dict(-1 => 1, 2 => 127)

@test factor(1) == Dict{Int,Int}()

# factor returns a sorted dict
@test all([issorted(collect(factor(rand(Int)))) for x in 1:100])

# Lucas-Lehmer
@test !ismersenneprime(2047)
@test ismersenneprime(8191)
@test_throws ArgumentError ismersenneprime(9)
@test_throws ArgumentError ismersenneprime(9, check=true)
# test the following does not throw
ismersenneprime(9, check=false)

#  Lucas-Lehmer-Riesel
@test_throws ArgumentError isrieselprime(1000, 511)
@test_throws ArgumentError isrieselprime(0, 1)
@test isrieselprime(1, 8191) == ismersenneprime(8191)  # Case 1
@test isrieselprime(3, BigInt(2)^607 - 1)              # Case 2
@test_throws ErrorException isrieselprime(20, 31)      # Case `else`

# @testset "Factorization{$T} as an AbstractDict" for T = (Int, UInt, BigInt)
for T = (Int, UInt, BigInt)
    d = Dict(map(Pair, rand(T(1):T(100), 30), 1:30))
    f = Factorization{T}(d)
    @test f == d == Dict(f) == Factorization(d) == convert(Factorization, d)
    @test collect(f) == sort!(collect(d)) # test start/next/done
    @test length(f) == length(d)
    @test get(f, T(101), nothing) == nothing
    @test f[101] == 0
    @test f[0] == 0
    f[0] = 1
    @test get(f, T(0), 0) == 1
    @test f[0] == 1
    @test f[101] == 0
    f = Factorization{T}()
    for i in 100:-1:1
        f[T(i)] = i
    end
    @test length(f) == 100
    @test issorted(f.pe)
end

# dumb implementation of Euler totient for correctness tests
ϕ(n) = count(m -> gcd(n, m) == 1, 1:n)

@testset "totient(::$T)" for T in [Int16, Int32, Int64, BigInt]
    for n in 1:1000
        @test ϕ(T(n)) == totient(T(n)) == totient(-T(n))
    end
    @test_throws ArgumentError @inferred(totient(T(0)))
end

# check some big values
@testset "totient() correctness" begin
    @test totient(2^4 * 3^4 * 5^4) == 216000
    @test totient(big"2"^1000) == big"2"^999

    some_coprime_numbers = BigInt[
        450000000, 1099427429702334733, 200252151854851, 1416976291499, 7504637909,
        1368701327204614490999, 662333585807659, 340557446329, 1009091
    ]

    for i in some_coprime_numbers
        for j in some_coprime_numbers
            if i ≠ j
                @test totient(i*j) == totient(i) * totient(j)
            end
        end
        # can use directly with Factorization
        @test totient(i) == totient(factor(i))
    end
end

# check copy property for big primes relied upon in nextprime/prevprime
for n = rand(big(-10):big(10), 10)
    @test n+0 !== n
end

@testset "nextprime(::$T)" for T in (Int64, Int32, BigInt)
    for N in zip(T[rand(-1000:1), 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
                   2^20, 2^30],
                 T[2, 2, 3, 5, 5, 7, 7, 11, 11, 11, 11, 13,
                   1048583, 1073741827],
                 T[3, 3, 5, 7, 7, 11, 11, 13, 13, 13, 13, 17,
                   1048589, 1073741831])
        @test nextprime(N[1]) == N[2]
        @test nextprime(N[1], 1) == N[2]
        @test nextprime(N[1], 2) == N[3]
    end
    @test nextprime(Int64(2)^60) == 1152921504606847009
    @test nextprime(Int64(2)^60, 2) == 1152921504606847067
    @test_throws DomainError nextprime(rand(-100:100), 0)
    for i = rand(-100:-1, 100),
        n = rand(600:2^20)
        @test nextprime(n, i) == prevprime(n, -i)
    end

    # interval
    @test nextprime(0, interval=2) == 2
    @test nextprime(0, interval=3) == 3
    @test_throws DomainError nextprime(0, 2, interval=2)
    @test_throws DomainError nextprime(0, interval=4)
    @test nextprime(-20, interval=2) == 2
    @test nextprime(4, interval=5) == 19
    @test nextprime(4, 2, interval=5) == 29
    @test nextprime(4, 3, interval=5) == 59
    @test gcd(nextprime(2^17-1024+1; interval=1024) - 1, 1024) == 1024
    @test gcd(nextprime(1024*rand(1:2^6)+1; interval=1024) - 1, 1024) == 1024
end

@testset "prevprime(::$T)" for T in (Int64, Int32, BigInt)
    for N in zip(T[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
                   2^20, 2^30],
                 T[3, 5, 5, 7, 7, 7, 7, 11, 11, 13, 13,
                   1048573, 1073741789],
                 T[2, 3, 3, 5, 5, 5, 5, 7, 7, 11, 11,
                   1048571, 1073741783])
        @test prevprime(N[1]) == N[2]
        @test prevprime(N[1], 1) == N[2]
        @test prevprime(N[1], 2) == N[3]
    end
    @test prevprime(Int64(2)^60) == 1152921504606846883
    @test prevprime(Int64(2)^60, 2) == 1152921504606846869
    @test prevprime(3) == 3
    @test prevprime(3, 2) == 2
    @test prevprime(2) == 2
    @test prevprime(typemax(UInt8)) == prevprime(big(typemax(UInt8)))
    @test typeof(prevprime(typemax(UInt8))) === UInt8
    @test_throws ArgumentError prevprime(2, 2)

    @test_throws DomainError prevprime(rand(-100:100), 0)
    for n = rand(-100:1, 10)
        @test_throws ArgumentError prevprime(n)
        @test_throws ArgumentError prevprime(n, 1)
        @test_throws ArgumentError prevprime(n, 2)
    end

    for i = rand(-100:-1, 100),
        n = rand(600:2^20)
        @test prevprime(n, i) == nextprime(n, -i)
    end

    # interval
    @test gcd(prevprime(2^17-1024+1; interval=1024) - 1, 1024) == 1024
    @test gcd(prevprime(1024*rand(12:2^6)+1; interval=1024) - 1, 1024) == 1024
    @test prevprime(10, interval=4) == 2
    @test [prevprime(11, i, interval=2) for i=1:4] == [11, 7, 5, 3]
    @test [prevprime(11, i, interval=3) for i=1:3] == [11, 5, 2]
end

@testset "prime(::$T)" for T = (Int64, Int32, BigInt)
    for (n, p) = zip([1, 2, 3, 4, 5, 6, 7, 100, 1000, 10000],
                     T[2, 3, 5, 7, 11, 13, 17, 541, 7919, 104729])
        @test prime(n) == p
        @test prime(T, n) == p
    end
end

for T in (Int, UInt, BigInt)
    for n in [T(1); rand(T(2):T(100000), 10)]
        # for n=T(1), must not error out (#51)
        for C = (Factorization, Vector, Dict)
            @test prodfactors(factor(C, n)) == n
        end
        if Primes.radical(n) == n
            for C = (Set, BitSet)
                @test prodfactors(factor(C, n)) == n
            end
        end
    end
    @test prodfactors(factor(Set, T(123456))) == 3858
    @test prod(factor(T(123456))) == 123456
end

@testset "nextprimes(::$T)" for T = (Int32, Int64, BigInt)
    for (i, p) in enumerate(nextprimes(T))
        @test nextprime(0, i) == p
        i > 20 && break
    end
    @test nextprimes() == nextprimes(Int)
    for (i, p) in enumerate(nextprimes(T(5)))
        @test nextprime(T(5), i) == p
        i > 20 && break
    end
    @test nextprimes(T(5), 10) == [nextprime(T(5), i) for i=1:10]
    @test nextprimes(1, 1)[1] == nextprimes(2, 1)[1] == 2
    @test nextprimes(3, 1)[1] == 3
    @test nextprimes(4, 1)[1] == nextprimes(5, 1)[1] == 5
    @test eltype(nextprimes(10)) == Int
    @test eltype(nextprimes(big(10))) == BigInt
    @test Base.IteratorEltype(nextprimes(10)) == Base.HasEltype()
    @test Base.IteratorSize(nextprimes(10)) == Base.IsInfinite()

end

@testset "prevprimes(::$T)" for T = (Int32, Int64, BigInt)
    for (i, p) in enumerate(prevprimes(T(500)))
        @test prevprime(T(500), i) == p
        i > 20 && break
    end
    @test prevprimes(T(500), 10) == [prevprime(T(500), i) for i=1:10]
    @test prevprimes(6, 1)[1] == prevprimes(5, 1)[1] == 5
    @test prevprimes(4, 1)[1] == prevprimes(3, 1)[1] == 3
    @test prevprimes(2, 1)[1] == 2
    @test isempty(prevprimes(1, 1))
    let p8 = collect(prevprimes(typemax(Int8)))
        @test length(p8) == 31
        @test p8[end] == 2
        @test p8[1] == 127
        @test eltype(p8) == Int8
    end
    @test eltype(prevprimes(10)) == Int
    @test eltype(prevprimes(big(10))) == BigInt
    @test Base.IteratorEltype(prevprimes(10)) == Base.HasEltype()
    @test Base.IteratorSize(prevprimes(10)) == Base.SizeUnknown()
end