text:
module Day08

using AdventOfCode2019

function day08(input::String = readInput(joinpath(@__DIR__, "..", "data", "day08.txt")))
    digits = parse.(UInt8, split(strip(input), ""))
    height = 6
    width = 25
    nPixels = height * width
    nLayers = length(digits) ÷ nPixels
    layerRanges = [1+(i-1)*nPixels:i*nPixels for i=1:nLayers]

    resultLayerIndex = -1
    minSum = nPixels
    for (i, r) in enumerate(layerRanges)
        s = (digits[r] .== 0) |> sum
        if s < minSum
            minSum = s
            resultLayerIndex = i
        end
    end

    nOneDigits = (digits[layerRanges[resultLayerIndex]] .== 1) |> sum
    nTwoDigits = (digits[layerRanges[resultLayerIndex]] .== 2) |> sum
    resultPart1 = nOneDigits * nTwoDigits

    resultImage = UInt8(2) * ones(UInt8, (width, height))
    for r in layerRanges
        for i = 1:nPixels
            if resultImage[i] in (0, 1) || digits[r][i] == 2
                continue
            end
            resultImage[i] = digits[r][i]
        end
    end

    return [resultPart1, resultImage]
end

function generate_image(image)
    block = '\u2588'
    empty = ' '
    output = ""
    for i = 1:size(image, 2)
        row = join(image[:, i])
        row = replace(row, "1" => block)
        row = replace(row, "0" => empty)
        output *= row * "\n"
    end
    return output
end

end # module


text:
export FixedParameterSystem, fix_parameters

"""
    FixedParameterSystem(F:AbstractSystem, parameters)

Construct a system from the given [`AbstractSystem`](@ref) `F` with the given `parameters`
fixed.
"""
struct FixedParameterSystem{S<:AbstractSystem,T} <: AbstractSystem
    system::S
    parameters::Vector{T}
end
FixedParameterSystem(
    F::AbstractSystem,
    p;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = FixedParameterSystem(F, p)
FixedParameterSystem(F::System, p; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[]) =
    FixedParameterSystem(fixed(F; compile = compile), p)
Base.size(F::FixedParameterSystem) = size(F.system)

ModelKit.variables(F::FixedParameterSystem) = variables(F.system)
ModelKit.parameters(F::FixedParameterSystem) = Variable[]
ModelKit.variable_groups(F::FixedParameterSystem) = variable_groups(F.system)

(F::FixedParameterSystem)(x, p = nothing) = F.system(x, F.parameters)

ModelKit.evaluate!(u, F::FixedParameterSystem, x, p = nothing) =
    evaluate!(u, F.system, x, F.parameters)
ModelKit.evaluate_and_jacobian!(u, U, F::FixedParameterSystem, x, p = nothing) =
    evaluate_and_jacobian!(u, U, F.system, x, F.parameters)
ModelKit.taylor!(u, v::Val, F::FixedParameterSystem, tx, p = nothing) =
    taylor!(u, v, F.system, tx, F.parameters)

# Only for certification
ModelKit.jacobian!(U, F::FixedParameterSystem{<:InterpretedSystem}, x, p = nothing) =
    ModelKit.jacobian!(U, F.system, x, F.parameters)
ModelKit.jacobian!(U, F::FixedParameterSystem{<:InterpretedSystem}, x, p, cache) =
    ModelKit.jacobian!(U, F.system, x, F.parameters, cache)

"""
    fix_parameters(F::Union{System,AbstractSystem}, p; compile::Union{Bool,Symbol} = $(COMPILE_DEFAULT[]))

Fix the parameters of the given system `F`. Returns a [`FixedParameterSystem`](@ref).
"""
fix_parameters(
    F::Union{System,AbstractSystem},
    p;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = FixedParameterSystem(F, p; compile = compile)

ext: jl


text:
module MyExample

using ForwardDiff
export my_f, derivative_of_my_f

greet() = print("hello world!")

include("extra_file.jl")# Write your package code here.

end


text:
@rule BIFM(:out, Marginalisation) (m_in::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, m_znext::MultivariateNormalDistributionsFamily, meta::BIFMMeta) = begin

    # fetch information from meta data
    A       = getA(meta)
    B       = getB(meta)
    C       = getC(meta)
    H       = getH(meta)
    ξztilde = getξztilde(meta)
    BHBt    = getBHBt(meta)
    Λz      = getΛz(meta)

    # fetch statistics
    μ_in, Σ_in       = mean_cov(m_in)
    μ_zprev, Σ_zprev = mean_cov(m_zprev)

    # calculate intermediate quantities
    F = I - Λz * BHBt

    μ_ztilde = A * μ_zprev
    Σ_ztilde = A * Σ_zprev * A'

    μ_znext = F' * μ_ztilde + (B * (Σ_in * (B' * ξztilde))) + B * μ_in
    Σ_znext = F' * Σ_ztilde * F + BHBt

    # save required intermediate variables
    setΣu!(meta, Σ_in)

    # calculate outgoing message to output
    μ_out = C * μ_znext
    Σ_out = C * Σ_znext * C'

    # Actual return type depends on meta object as well, so we explicitly cast the result here
    # Should be noop if type matches
    T = promote_type(eltype(m_in), eltype(m_zprev), eltype(m_znext))

    # return outgoing marginal
    return ProdFinal(convert(MvNormalMeanCovariance{T}, MvNormalMeanCovariance(μ_out, Σ_out)))

end


text:
import Base: ==, isapprox
using LinearAlgebra: norm

# Comparison
(==)(a::KnetArray{T},b::KnetArray{T}) where {T}=(size(a)==size(b) && sum(abs2,a-b)==0)
(==)(a::AbstractArray,b::KnetArray)=(size(a)==size(b) && a==Array(b))
(==)(a::KnetArray,b::AbstractArray)=(size(a)==size(b) && Array(a)==b)
# Adapted from base/linalg/generic.jl:589
isapprox(a::AbstractArray,b::KnetArray;o...)=(size(a)==size(b) && isapprox(a,Array(b);o...))
isapprox(a::KnetArray,b::AbstractArray;o...)=(size(a)==size(b) && isapprox(Array(a),b;o...))

function isapprox(a::KnetArray{T}, b::KnetArray{T}; rtol=sqrt(eps(T)), atol=T(0)) where {T}
    (size(a)==size(b) && norm(a-b) <= atol + rtol * max(norm(a), norm(b)))
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# ––––––––––
# Paragraphs
# ––––––––––

type Paragraph
    content
end

Paragraph() = Paragraph([])

function paragraph(stream::IO, md::MD)
    buffer = IOBuffer()
    p = Paragraph()
    push!(md, p)
    skipwhitespace(stream)
    while !eof(stream)
        char = read(stream, Char)
        if char == '\n' || char == '\r'
            if blankline(stream) || parse(stream, md, breaking = true)
                break
            else
                write(buffer, ' ')
            end
        else
            write(buffer, char)
        end
    end
    p.content = parseinline(seek(buffer, 0), md)
    return true
end

# –––––––
# Headers
# –––––––

type Header{level}
    text
end

Header(s, level::Int) = Header{level}(s)
Header(s) = Header(s, 1)

@breaking true ->
function hashheader(stream::IO, md::MD)
    withstream(stream) do
        eatindent(stream) || return false
        level = 0
        while startswith(stream, '#') level += 1 end
        level < 1 || level > 6 && return false

        c = ' '
        # Allow empty headers, but require a space
        !eof(stream) && (c = read(stream, Char); !(c in " \n")) &&
            return false

        if c != '\n' # Empty header
            h = readline(stream) |> strip
            h = match(r"(.*?)( +#+)?$", h).captures[1]
            buffer = IOBuffer()
            print(buffer, h)
            push!(md.content, Header(parseinline(seek(buffer, 0), md), level))
        else
            push!(md.content, Header("", level))
        end
        return true
    end
end

function setextheader(stream::IO, md::MD)
    withstream(stream) do
        eatindent(stream) || return false
        header = readline(stream) |> strip
        header == "" && return false

        eatindent(stream) || return false
        underline = readline(stream) |> strip
        length(underline) < 3 && return false
        u = underline[1]
        u in "-=" || return false
        all(c -> c == u, underline) || return false
        level = (u == '=') ? 1 : 2

        push!(md.content, Header(parseinline(header, md), level))
        return true
    end
end

# ––––
# Code
# ––––

type Code
    language::UTF8String
    code::UTF8String
end

Code(code) = Code("", code)

function indentcode(stream::IO, block::MD)
    withstream(stream) do
        buffer = IOBuffer()
        while !eof(stream)
            if startswith(stream, "    ") || startswith(stream, "\t")
                write(buffer, readline(stream))
            elseif blankline(stream)
                write(buffer, '\n')
            else
                break
            end
        end
        code = takebuf_string(buffer)
        !isempty(code) && (push!(block, Code(rstrip(code))); return true)
        return false
    end
end

# ––––––
# Quotes
# ––––––

type BlockQuote
    content
end

BlockQuote() = BlockQuote([])

# TODO: Laziness
@breaking true ->
function blockquote(stream::IO, block::MD)
    withstream(stream) do
        buffer = IOBuffer()
        empty = true
        while eatindent(stream) && startswith(stream, '>')
            startswith(stream, " ")
            write(buffer, readline(stream))
            empty = false
        end
        empty && return false

        md = takebuf_string(buffer)
        push!(block, BlockQuote(parse(md, flavor = config(block)).content))
        return true
    end
end

# –––––
# Lists
# –––––

type List
    items::Vector{Any}
    ordered::Bool

    List(x::AbstractVector, b::Bool) = new(x, b)
    List(x::AbstractVector) = new(x, false)
    List(b::Bool) = new(Any[], b)
end

List(xs...) = List(vcat(xs...))

const bullets = "*•+-"
const num_or_bullets = r"^(\*|•|\+|-|\d+(\.|\))) "

# Todo: ordered lists, inline formatting
function list(stream::IO, block::MD)
    withstream(stream) do
        eatindent(stream) || return false
        b = startswith(stream, num_or_bullets)
        (b == nothing || b == "") && return false
        ordered = !(b[1] in bullets)
        if ordered
            b = b[end - 1] == '.' ? r"^\d+\. " : r"^\d+\) "
            # TODO start value
        end
        the_list = List(ordered)

        buffer = IOBuffer()
        fresh_line = false
        while !eof(stream)
            if fresh_line
                sp = startswith(stream, r"^ {0,3}")
                if !(startswith(stream, b) in [false, ""])
                    push!(the_list.items, parseinline(takebuf_string(buffer), block))
                    buffer = IOBuffer()
                else
                    # TODO write a newline here, and deal with nested
                    write(buffer, ' ', sp)
                end
                fresh_line = false
            else
                c = read(stream, Char)
                if c == '\n'
                    eof(stream) && break
                    next = peek(stream)
                    if next == '\n'
                        break
                    else
                        fresh_line = true
                    end
                else
                    write(buffer, c)
                end
            end
        end
        push!(the_list.items, parseinline(takebuf_string(buffer), block))
        push!(block, the_list)
        return true
    end
end

# ––––––––––––––
# HorizontalRule
# ––––––––––––––

type HorizontalRule
end

function horizontalrule(stream::IO, block::MD)
   withstream(stream) do
       n, rule = 0, ' '
       while !eof(stream)
           char = read(stream, Char)
           char == '\n' && break
           isspace(char) && continue
           if n==0 || char==rule
               rule = char
               n += 1
           else
               return false
           end
       end
       is_hr = (n ≥ 3 && rule in "*-")
       is_hr && push!(block, HorizontalRule())
       return is_hr
   end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -6.103907844236108 <= q <= 8.378675476959124)

                     add_NL_constraint(m, :(gelu(-0.8297607126794757 + -0.5064858067335005*gelu(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8619658571743893*gelu(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + 0.8415580445654567*gelu(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + gelu(-0.0809505831099564 + 0.6185835236357913*gelu(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.5288114683363174*gelu(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.42070009030061817*gelu(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) + gelu(0.16011867292674786 + -0.14898607305663658*gelu(-0.6188655598844357 + -0.38627203797486676*$(x[1]) + 0.4878933200630442*$(x[2])) + -0.8220978033005912*gelu(-0.7914291135758744 + 0.8682104207735972*$(x[1]) + 0.3684602934879275*$(x[2])) + -0.8518875556726564*gelu(-0.26405078449033637 + -0.6452849815457378*$(x[1]) + 0.9586530591287103*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
## -- Binning functions, in-place

    """
    ```julia
    nanbinmean!(MU, [N], x, y, xedges::AbstractRange)
    ```
    Ignoring NaNs, fill the array `MU` with the means (and optionally `N` with
    the counts) of non-NAN `y` values that fall into each of `length(xedges)-1`
    equally spaced bins along the `x` axis with bin edges specified by `xedges`.

    The array of `x` data should given as a one-dimensional array (any subtype
    of AbstractVector) and `y` as either a 1-d or 2-d array (any subtype of
    AbstractVecOrMat).

    The output arrays `MU` and `N` must be the same size, and must have the same
    number of columns as `y`; if `y` is a 2-d array (matrix), then each column of
    `y` will be treated as a separate variable.
    """
    function nanbinmean!(MU::AbstractVecOrMat, x::AbstractVector, y::AbstractVecOrMat, xedges::AbstractRange)
        N = Array{Int}(undef, size(MU))
        return nanbinmean!(MU, N, x, y, xedges)
    end
    # As above, but also return an array of counts, N; y, N, and MU as 1D vectors
    function nanbinmean!(MU::AbstractVector, N::AbstractVector, x::AbstractVector, y::AbstractVector, xedges::AbstractRange)
        # What is the size of each bin?
        nbins = length(xedges) - 1
        xmin, xmax = extrema(xedges)
        δ𝑖δx = nbins / (xmax - xmin)

        # Make sure we don't have a segfault by filling beyond the length of N
        # in the @inbounds loop below
        if length(MU) < nbins
            nbins = length(MU)
            @warn "length(MU) < nbins; any bins beyond length(MU) will not be filled"
        end

        # Calculate the means for each bin, ignoring NaNs
        fill!(N, 0)
        fill!(MU, 0) # Fill the output array with zeros to start
        @inbounds for n ∈ eachindex(x)
            𝑖 = (x[n] - xmin) * δ𝑖δx
            yₙ = y[n]
            if (0 < 𝑖 < nbins) && (yₙ==yₙ)
                i = ceil(Int, 𝑖)
                N[i] += 1
                MU[i] += yₙ
            end
        end
        MU ./= N # Divide by N to calculate means. Empty bin = 0/0 = NaN

        return MU
    end
    # In-place binning; as above but with y, N, MU as 2D matrices instead of 1D vectors
    function nanbinmean!(MU::AbstractMatrix, N::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, xedges::AbstractRange)
        # What is the size of each bin?
        nbins = length(xedges) - 1
        xmin, xmax = extrema(xedges)
        δ𝑖δx = nbins / (xmax - xmin)
        ncols = size(y,2)

        # Make sure we don't have a segfault by filling beyond the length of N
        # in the @inbounds loop below
        if size(MU,1) < nbins
            nbins = size(MU,1)
            @warn "size(MU,1) < nbins; any bins beyond size(MU,1) will not be filled"
        end
        if size(MU,2) < ncols
            ncols = size(MU,2)
            @warn "size(MU,2) < size(y,2); any columns beyond size(MU,2) will not be filled"
        end

        # Calculate the means for each bin, ignoring NaNs
        fill!(N, 0)
        fill!(MU, 0) # Fill the output array with zeros to start
        @inbounds for n ∈ eachindex(x)
            𝑖 = (x[n] - xmin) * δ𝑖δx
            if (0 < 𝑖 < nbins)
                i = ceil(Int, 𝑖)
                for j ∈ 1:ncols
                    if y[n,j]==y[n,j]
                        N[i,j] += 1
                        MU[i,j] += y[n,j]
                    end
                end
            end
        end
        MU ./= N # Divide by N to calculate means. Empty bin = 0/0 = NaN

        return MU
    end
    nanbinmean!(MU, x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmean!(MU, x, y, range(xmin,xmax,length=nbins+1))
    nanbinmean!(MU, N, x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmean!(MU, N, x, y, range(xmin,xmax,length=nbins+1))
    export nanbinmean!

    # 2-d binning
    """
    ```julia
    nanbinmean!(MU, N, x, y, z, xedges::AbstractRange, yedges::AbstractRange)
    ```
    Ignoring NaNs, fill the matrix `MU` with the means and `N` with
    the counts of non-NAN `z` values that fall into a 2D grid of x and y bins
    defined by `xedges` and `yedges`. The independent variables `x` and `y`,
    as well as the dependent variable `z`, are all expected as 1D vectors (any
    subtype of AbstractVector).

    The output matrices `MU` and `N` must be the same size, and must each have
    `length(yedges)-1` rows and `length(xedges)-1` columns.
    """
    function nanbinmean!(MU::AbstractMatrix, N::AbstractMatrix, x::AbstractVector, y::AbstractVector, z::AbstractVector, xedges::AbstractRange, yedges::AbstractRange)
        # Calculate bin index from x value
        nxbins = length(xedges)-1
        xmin, xmax = extrema(xedges)
        δjδx = nxbins / (xmax - xmin)

        nybins = length(yedges)-1
        ymin, ymax = extrema(yedges)
        δiδy = nybins / (ymax - ymin)

        # Make sure we don't have a segfault by filling beyond the length of N
        # in the @inbounds loop below
        if size(MU,1) < nybins
            nybins = size(MU,1)
            @warn "size(MU,1) < nybins; any y bins beyond size(MU,1) will not be filled"
        end
        if size(MU,2) < nxbins
            nxbins = size(MU,2)
            @warn "size(MU,2) < nxbins; any x bins beyond size(MU,2) will not be filled"
        end

        # Calculate the means for each bin, ignoring NaNs
        fill!(N, 0)
        fill!(MU, 0) # Fill the output array with zeros to start
        @inbounds for n ∈ eachindex(x)
            𝑖 = (y[n] - ymin) * δiδy
            𝑗 = (x[n] - xmin) * δjδx
            zₙ = z[n]
            if (zₙ==zₙ) && (0 < 𝑖 < nybins) && (0 < 𝑗 < nxbins)
                i = ceil(Int, 𝑖)
                j = ceil(Int, 𝑗)
                N[i,j] += 1
                MU[i,j] += zₙ
            end
        end
        MU ./= N # Divide by N to calculate means. Empty bin = 0/0 = NaN

        return MU
    end

    """
    ```julia
    nanbinmean!(MU, W, x, y, w, xedges::AbstractRange)
    ```
    Ignoring NaNs, fill the array `MU` with the weighted means (and `W` with
    the sum of weights) of non-NAN `y` values that fall into each of
    `length(xedges)-1` equally spaced bins along the `x` axis with bin edges
    specified by `xedges`.

    The array of `x` data should given as a one-dimensional array (any subtype
    of AbstractVector) and `y` as either a 1-d or 2-d array (any subtype of
    AbstractVecOrMat).

    The output arrays `MU` and `N` must be the same size, and must have the same
    number of columns as `y`; if `y` is a 2-d array (matrix), then each column of
    `y` will be treated as a separate variable.
    """
    # In-place binning with weights; y, W, MU as 1D vectors
    function nanbinmean!(MU::AbstractVector, W::AbstractVector, x::AbstractVector, y::AbstractVector, w::AbstractVector, xedges::AbstractRange)
        # What is the size of each bin?
        nbins = length(xedges) - 1
        xmin, xmax = extrema(xedges)
        δ𝑖δx = nbins / (xmax - xmin)

        # Make sure we don't have a segfault by filling beyond the length of N
        # in the @inbounds loop below
        if length(MU) < nbins
            nbins = length(MU)
            @warn "length(MU) < nbins; any bins beyond length(MU) will not be filled"
        end

        # Calculate the means for each bin, ignoring NaNs
        fill!(W, 0)
        fill!(MU, 0) # Fill the output array with zeros to start
        @inbounds for n ∈ eachindex(x)
            𝑖 = (x[n] - xmin) * δ𝑖δx
            yₙ = y[n]
            if (0 < 𝑖 < nbins) && (yₙ==yₙ)
                i = ceil(Int, 𝑖)
                W[i] += w[n]
                MU[i] += yₙ * w[n]
            end
        end
        MU ./= W # Divide by sum of weights to calculate means. Empty bin = 0/0 = NaN

        return MU
    end
    # In-place binning with weights; y, W, MU as 2D matrices
    function nanbinmean!(MU::AbstractMatrix, W::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, w::AbstractVector, xedges::AbstractRange)
        nbins = length(xedges) - 1
        xmin, xmax = extrema(xedges)
        δ𝑖δx = nbins / (xmax - xmin)
        ncols = size(y,2)

        # Make sure we don't have a segfault by filling beyond the length of N
        # in the @inbounds loop below
        if size(MU,1) < nbins
            nbins = size(MU,1)
            @warn "size(MU,1) < nbins; any bins beyond size(MU,1) will not be filled"
        end
        if size(MU,2) < ncols
            ncols = size(MU,2)
            @warn "size(MU,2) < size(y,2); any columns beyond size(MU,2) will not be filled"
        end

        # Calculate the means for each bin, ignoring NaNs
        fill!(W, 0)
        fill!(MU, 0) # Fill the output array with zeros to start
        @inbounds for n ∈ eachindex(x)
            𝑖 = (x[n] - xmin) * δ𝑖δx
            if (0 < 𝑖 < nbins)
                i = ceil(Int, 𝑖)
                for j ∈ 1:ncols
                    if y[n,j]==y[n,j]
                        W[i,j] += w[n]
                        MU[i,j] += y[n,j]*w[n]
                    end
                end
            end
        end
        MU ./= W # Divide by sum of weights to calculate means. Empty bin = 0/0 = NaN

        return MU
    end
    nanbinmean!(MU, W, x, y, w, xmin::Number, xmax::Number, nbins::Integer) = nanbinmean!(MU, W, x, y, w, range(xmin,xmax,length=nbins+1))

    """
    ```julia
    nanbinmedian!(M, [N], x, y, xedges::AbstractRange)
    ```
    Fill the array `M` with the medians (and optionally `N` with the counts) of
    non-NaN `y` values that fall into each of `length(xedges)-1` equally spaced
    bins along the `x` axis with bin edges specified by `xedges`.

    If `y` is a 2-d array (matrix), each column will be treated as a separate variable
    """
    function nanbinmedian!(M::AbstractVector, x::AbstractVector, y::AbstractVector, xedges::AbstractRange)
        nbins = length(xedges) - 1
        t = Array{Bool}(undef, length(x))
        @inbounds for i = 1:nbins
            t .= (xedges[i] .<= x .< xedges[i+1])
            M[i] = _nanmedian!(y[t],:)
        end
        return M
    end
    function nanbinmedian!(M::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, xedges::AbstractRange)
        nbins = length(xedges) - 1
        t = Array{Bool}(undef, length(x))
        @inbounds for i = 1:nbins
            t .= xedges[i] .<= x .< xedges[i+1]
            if any(t)
                for j = 1:size(y,2)
                    M[i,j] = _nanmedian!(y[t,j],:)
                end
            else
                M[i,:] .= float(eltype(y))(NaN)
            end
        end
        return M
    end
    function nanbinmedian!(M::AbstractVector, N::AbstractVector, x::AbstractVector, y::AbstractVector, xedges::AbstractRange)
        nbins = length(xedges) - 1
        t = Array{Bool}(undef, length(x))
        @inbounds for i = 1:nbins
            t .= (xedges[i] .<= x .< xedges[i+1]) .& (y.==y)
            yₜ = y[t]
            M[i] = _nanmedian!(yₜ,:)
            N[i] = countnotnans(yₜ)
        end
        return M
    end
    function nanbinmedian!(M::AbstractMatrix, N::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, xedges::AbstractRange)
        nbins = length(xedges) - 1
        t = Array{Bool}(undef, length(x))
        tj = Array{Bool}(undef, length(x))
        @inbounds for i = 1:nbins
            t .= xedges[i] .<= x .< xedges[i+1]
            if any(t)
                for j = 1:size(y,2)
                    yₜ = y[t,j]
                    M[i,j] = _nanmedian!(yₜ,:)
                    N[i,j] = countnotnans(yₜ)
                end
            else
                M[i,j] = float(eltype(y))(NaN)
                N[i,j] = 0
            end
        end
        return M
    end
    nanbinmedian!(MU, x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmedian!(MU, x, y, range(xmin,xmax,length=nbins+1))
    nanbinmedian!(MU, N, x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmedian!(MU, N, x, y, range(xmin,xmax,length=nbins+1))
    export nanbinmedian!

## -- Binning functions, convenience versions (not in place)

    """
    ```julia
    nanbinmean(x, y, xedges::AbstractRange)
    ```
    Ignoring NaNs, calculate the mean of `y` values that fall into each of
    `length(xedges)-1` equally spaced bins along the `x` axis with bin edges
    specified by `xedges`.

    The array of `x` data should be given as a one-dimensional array (any subtype
    of AbstractVector) and `y` as either a 1-d or 2-d array (any subtype of
    AbstractVecOrMat). If `y` is a 2-d array, then each column of `y` will be
    treated as a separate variable.

    ## Examples
    ```julia
    julia> nanbinmean([1:100..., 1], [1:100..., NaN], 0:25:100)
    4-element Vector{Float64}:
     13.0
     38.0
     63.0
     87.5

    julia> nanbinmean(1:100, reshape(1:300,100,3), 0:25:100)
    4×3 Matrix{Float64}:
     13.0  113.0  213.0
     38.0  138.0  238.0
     63.0  163.0  263.0
     87.5  187.5  287.5
    ```
    """
    function nanbinmean(x::AbstractVector, y::AbstractVecOrMat, xedges::AbstractRange)
        N = Array{Int}(undef, length(xedges)-1, size(y)[2:end]...)
        MU = Array{float(eltype(y))}(undef, length(xedges)-1, size(y)[2:end]...)
        return nanbinmean!(MU, N, x, y, xedges)
    end
    nanbinmean(x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmean(x, y, range(xmin,xmax,length=nbins+1))
    export nanbinmean

    """
    ```julia
    nanbinmean(x, y, z, xedges, yedges)
    ```
    Ignoring NaNs, calculate the mean of `z` values that fall into a 2D grid of
    x and y bins with bin edges defined by `xedges` and `yedges`. The independent
    variables `x` and `y`, as well as the dependent variable `z`, are all expected
    as 1D vectors (any subtype of AbstractVector).

    ## Examples
    ```julia
    julia> x = y = z = 0.5:9.5;

    julia> xedges = yedges = 0:10;

    julia> nanbinmean(x,y,z,xedges,yedges)
    10×10 Matrix{Float64}:
       0.5  NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN
     NaN      1.5  NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN
     NaN    NaN      2.5  NaN    NaN    NaN    NaN    NaN    NaN    NaN
     NaN    NaN    NaN      3.5  NaN    NaN    NaN    NaN    NaN    NaN
     NaN    NaN    NaN    NaN      4.5  NaN    NaN    NaN    NaN    NaN
     NaN    NaN    NaN    NaN    NaN      5.5  NaN    NaN    NaN    NaN
     NaN    NaN    NaN    NaN    NaN    NaN      6.5  NaN    NaN    NaN
     NaN    NaN    NaN    NaN    NaN    NaN    NaN      7.5  NaN    NaN
     NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN      8.5  NaN
     NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN    NaN      9.5
    ```
    """
    function nanbinmean(x::AbstractVector, y::AbstractVector, z::AbstractVector, xedges::AbstractRange, yedges::AbstractRange)
        N = Array{Int}(undef, length(yedges)-1, length(xedges)-1)
        MU = Array{float(eltype(y))}(undef, length(yedges)-1, length(xedges)-1)
        return nanbinmean!(MU, N, x, y, z, xedges, yedges)
    end

    """
    ```julia
    nanbinmean(x, y, xedges::AbstractRange)
    ```
    Ignoring NaNs, calculate the weighted mean of `y` values that
    fall into each of `length(xedges)-1` equally spaced bins along the `x`
    axis with bin edges specified by `xedges`.

    The array of `x` data should given as a one-dimensional array (any subtype
    of AbstractVector) and `y` as either a 1-d or 2-d array (any subtype of
    AbstractVecOrMat). If `y` is a 2-d array, then each column of `y` will be
    treated as a separate variable.
    """
    function nanbinmean(x::AbstractVector, y::AbstractVecOrMat, w::AbstractVector, xedges::AbstractRange)
        W = Array{float(eltype(y))}(undef, length(xedges)-1, size(y)[2:end]...)
        MU = Array{float(eltype(y))}(undef, length(xedges)-1, size(y)[2:end]...)
        return nanbinmean!(MU, W, x, y, w, xedges)
    end
    nanbinmean(x, y, w, xmin::Number, xmax::Number, nbins::Integer) = nanbinmean(x, y, w, range(xmin,xmax,length=nbins+1))

    """
    ```julia
    nanbinmedian(x, y, xedges::AbstractRange)
    ```
    Calculate the median, ignoring NaNs, of y values that fall into each of
    `length(xedges)-1` equally spaced bins along the `x` axis with bin edges
    specified by `xedges`.

    If `y` is a 2-d array (matrix), each column will be treated as a separate variable

    ## Examples
    ```julia
    julia> nanbinmedian([1:100..., 1], [1:100..., NaN], 0:25:100)
    4-element Vector{Float64}:
     12.5
     37.0
     62.0
     87.0

    julia> nanbinmedian(1:100, reshape(1:300,100,3), 0:25:100)
    4×3 Matrix{Float64}:
     12.5  112.5  212.5
     37.0  137.0  237.0
     62.0  162.0  262.0
     87.0  187.0  287.0
    ```
    """
    function nanbinmedian(x::AbstractVector, y::AbstractVecOrMat, xedges::AbstractRange)
        M = Array{float(eltype(y))}(undef, length(xedges)-1, size(y)[2:end]...)
        return nanbinmedian!(M, x, y, xedges)
    end
    nanbinmedian(x, y, xmin::Number, xmax::Number, nbins::Integer) = nanbinmedian(x, y, range(xmin,xmax,length=nbins+1))
    export nanbinmedian

## -- End of File


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
using JLD2;
using Printf;
using StatsPlots;
include("runit.jl"); # for types
include("../utilities/experiment_helpers.jl");
include("../utilities/thresholds.jl");

name = ARGS[1];

@load "$name.dat" dist μ pep srs data δs βs N seed

dump_stats(pep, μ, δs, βs, srs, data, N);

for i in 1:length(δs)
    plot(_boxes(pep, μ, δs[i], βs[i], srs, getindex.(data, i), N));
    savefig("$(name)_$(δs[i])_no_oulier.pdf");
end


text:
using VUMPS
using Test

@testset "VUMPS.jl" begin
    @testset "cuda_patch" begin
        println("cuda_patch tests running...")
        include("cuda_patch.jl")
    end

    @testset "environment" begin
        println("environment tests running...")
        include("environment.jl")
    end

    @testset "fixedpoint" begin
        println("fixedpoint tests running...")
        include("fixedpoint.jl")
    end

    @testset "vumpsruntime.jl" begin
        println("vumpsruntime tests running...")
        include("vumpsruntime.jl")
    end

    @testset "autodiff.jl" begin
        println("autodiff tests running...")
        include("autodiff.jl")
    end
end


text:
# using MatpowerCases, Base.Test
using Test, MatpowerCases

println("Running tests...")

include("misc.jl")
include("brch.jl")
include("bus.jl")
include("gen.jl")
include("gencost.jl")
include("case.jl")


text:
using Markdown: Markdown

# function comment_out(code)
#     "#" * replace(code, "\n" => "\n#")
# end

function convert_to_rcall(code)
    "\n" * "R" * '"'^3  * "\n" * code * "\n" * '"' ^ 3 * "\n"
end

function isrcode(code)
    try
        if code.language == 'r'
            return true
        elseif code.language[2] == 'r'
            return true
        else
            return false
        end
    catch
        return false
    end
end

function weave(ctx, code::Markdown.Code)
   if isrcode(code)
       return Markdown.Code(
           "julia",
           convert_to_rcall(code.code)
       )
   else
       code
   end
end

function weave(ctx, doc::Markdown.MD)
    return Markdown.MD(
        weave.(ctx, doc.content),
        doc.meta
    )
end

weave(ctx, x) = x

function main(ctx, src=joinpath(@__DIR__, "Pest_Control_Example.Rmd"), dst=joinpath(@__DIR__, "notebook.md"))
    parsed = Markdown.parse(read(src, String))
    weaved = string(weave(ctx, parsed))
    __fixed = replace(weaved, r":\$\n\n" => "\$\$\n")
    _fixed = replace(__fixed, r":\$" => "\n\$\$\n")
    fixed = replace(_fixed, r" \$\$" => "\n\$\$\n")
    open(dst, "w") do f
        print(f, fixed)
    end
end


text:
using Turing, Distributions, Test
using MacroTools

model_info = Dict(:main_body_names => Dict(:vi => :vi, :sampler => :sampler))
# unit test model macro
expr = Turing.generate_observe(:x, :y, model_info)
@test expr.head == :block
@test :(vi.logp += Turing.observe(sampler, y, x, vi)) in expr.args

@model testmodel_comp(x, y) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x ~ Normal(m, sqrt(s))
    y ~ Normal(m, sqrt(s))

    return x, y
end
testmodel_comp(1.0, 1.2)

# check if drawing from the prior works
@model testmodel0(x) = begin
    x ~ Normal()
    return x
end
f0_mm = testmodel0()
@test mean(f0_mm() for _ in 1:1000) ≈ 0. atol=0.1

# Test #544
@model testmodel0(x = Vector{Float64}(undef, 2)) = begin
    x[1] ~ Normal()
    x[2] ~ Normal()
    return x
end
f0_mm = testmodel0()
@test all(isapprox.(mean(f0_mm() for _ in 1:1000), 0., atol=0.1))

@model testmodel01(x) = begin
    x ~ Bernoulli(0.5)
    return x
end
f01_mm = testmodel01()
@test mean(f01_mm() for _ in 1:1000) ≈ 0.5 atol=0.1

# test if we get the correct return values
@model testmodel1(x1, x2) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x1 ~ Normal(m, sqrt(s))
    x2 ~ Normal(m, sqrt(s))

    return x1, x2
end
f1_mm = testmodel1(1., 10.)
@test f1_mm() == (1, 10)

# Test for assertions in observe statements.
@model brokentestmodel_observe1(x1, x2) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x1 ~ Normal(m, sqrt(s))
    x2 ~ x1 + 2

    return x1, x2
end

btest = brokentestmodel_observe1(1., 2.)
@test_throws ArgumentError btest()

@model brokentestmodel_observe2(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x = Vector{Float64}(undef, 2)
    x ~ [Normal(m, sqrt(s)), 2.0]

    return x
end

btest = brokentestmodel_observe2([1., 2.])
@test_throws ArgumentError btest()

# Test for assertions in assume statements.
@model brokentestmodel_assume1() = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x1 ~ Normal(m, sqrt(s))
    x2 ~ x1 + 2

    return x1, x2
end

btest = brokentestmodel_assume1()
@test_throws ArgumentError btest()

@model brokentestmodel_assume2() = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))

    x = Vector{Float64}(undef, 2)
    x ~ [Normal(m, sqrt(s)), 2.0]

    return x
end

btest = brokentestmodel_assume2()
@test_throws ArgumentError btest()


text:
module DiffPong

include("data_structures.jl")
export Arena, Ball, Env, Player, Vector_

include("utils.jl")
export render, draw_rectangle!, get_obs,
       get_score, update_paddle_position

include("physics.jl")
export wall_collision!, paddle_collision!, rebound

include("game_logic.jl")
export reset_ball!, human_action, ai_action, step!, reset!

include("gradients.jl")
end # module


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

"""
    prime_range(n)

Return an array consisting of all primes up to and including `n`.
"""
function prime_range(n::Int)

   primes = fill(true, n)
   if n == 0
      return primes
   end

   primes[1] = false
   for p in 2:n
      primes[p] || continue
      for j in 2:div(n, p)
         primes[p*j] = false
       end
   end

   return findall(primes)
end


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

"""
    Leapfrog

Rotation convention of the Leapfrog software.
"""
struct Leapfrog <: RotationConvention end
axesseq(::Type{Leapfrog}) = :ZXZ
orientation(::Type{Leapfrog}) = (:CW,:CW,:CW)
angleunits(::Type{Leapfrog}) = :DEG
mainaxis(::Type{Leapfrog}) = :X
isextrinsic(::Type{Leapfrog}) = false


text:
# Particularly useful or beautiful grammar of graphics invocations.

# A convenience plot function for quickly plotting functions or expressions.
#
# Args:
#   fs: An array in which each object is either a single argument function or an
#       expression computing some value on x.
#   a: Lower bound on x.
#   b: Upper bound on x.
#   elements: One ore more grammar elements.
#
# Returns:
#   A plot objects.
#
function plot(fs::Array, a, b, elements::ElementOrFunction...; mapping...)
    # Catch a common misuse of this function
    if isa(b, ElementOrFunction)
        error(
        """
        Invalid plot usage:
            plot(xs, ys, ...) should be plot(x=xs, y=ys, ...)
        """)
    end

    if isempty(elements)
        elements = ElementOrFunction[]
    elseif isa(elements, Tuple)
        elements = ElementOrFunction[elements...]
    end

    element_types = Set(map(typeof, elements))

    if !in(Guide.xlabel, element_types)
        push!(elements, Guide.xlabel("x"))
    end

    if !in(Guide.ylabel, element_types)
        push!(elements, Guide.ylabel("f(x)"))
    end

    if b < a
        push!(elements, Coord.cartesian(xflip=true))
    end

    mappingdict = @compat Dict{Symbol, Any}(:y => fs, :xmin => [a], :xmax => [b])
    for (k, v) in mapping
        mappingdict[k] = v
    end

    plot(Stat.func, Geom.line, elements...; mappingdict...)
end

# Plot a single function.
function plot(f::Function, a, b, elements::ElementOrFunction...; mapping...)
    plot(Function[f], a, b, elements...; mapping...)
end

# Plot a single function using a contour plot
function plot(f::Function, xmin, xmax, ymin, ymax,
              elements::ElementOrFunction...; mapping...)
    if isempty(elements)
        elements = ElementOrFunction[]
    elseif isa(elements, Tuple)
        elements = ElementOrFunction[elements...]
    end

    element_types = Set(map(typeof, elements))

    if !in(Guide.xlabel, element_types)
        push!(elements, Guide.xlabel("x"))
    end

    if !in(Guide.ylabel, element_types)
        push!(elements, Guide.ylabel("y"))
    end

    if !in(Guide.ColorKey, element_types) && !in(Guide.ManualColorKey, element_types)
        push!(elements, Guide.colorkey("f(x,y)"))
    end

    push!(elements, Coord.cartesian(xflip=xmin > xmax, yflip=ymin > ymax))

    mappingdict = @compat Dict{Symbol, Any}(:z    => f, :xmin => [xmin], :xmax => [xmax],
                                            :ymin => [ymin], :ymax => [ymax])
    for (k, v) in mapping
        mappingdict[k] = v
    end

    plot(Geom.contour, elements...; mappingdict...)
end

# Create a layer from a list of functions or expressions.
function layer(fs::Array, a, b, elements::ElementOrFunction...)
    layer(y=fs, xmin=[a], xmax=[b], Stat.func, Geom.line, elements...)
end

# Create a layer from a single function.
function layer(f::Function, a, b, elements::ElementOrFunction...)
    layer([f], a, b, elements...)
end

# Simple heatmap plots of matrices.
#
# Args:
#   M: A matrix.
#
# Returns:
#   A plot object.
#
function spy(M::AbstractMatrix, elements::ElementOrFunction...; mapping...)
    is, js, values = findnz(M)
    df = DataFrame(i=is, j=js, value=values)
    plot(df, x="j", y="i", color="value",
         Coord.cartesian(yflip=true),
         Scale.color_continuous,
         Scale.x_continuous,
         Scale.y_continuous,
         Geom.rectbin,
         Stat.identity,
         elements...;
         mapping...)
end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test TreeNode" begin

    # Initial a TreeNode with data '1'
    T_e = Int64
    tn = DS.TreeNode{T_e}(1)

    @test DS.TreeNode(1) isa DS.TreeNode{Int64}
    @test DS.value(tn) == 1
    @test DS.parent(tn) isa DS.NullNode
    @test DS.left_child(tn) isa DS.NullNode
    @test DS.right_child(tn) isa DS.NullNode
    @test DS.level(tn) == 0
    @test eltype(tn) == T_e

    @test repr(tn) == "TreeNode{$T_e}(1)"
end

@testset "Test BinaryTree" begin

    # Initial a Binary Tree with data '1'
    bt = DS.BinaryTree(1)
    @test DS.value(DS.root(bt)) == 1
    @test length(bt) == 1
    @test repr(bt) == "BinaryTree(\n\nTreeNode(1)\n)"

    bt = DS.BinaryTree()
    @test DS.value(DS.root(bt)) == nothing
    @test length(bt) == 0
    @test repr(bt) == "BinaryTree(\n)"

    @test eltype(bt) == Int64

    # Set index Test: Push 1:27 value into tree
    for i=1:27
        bt[i] = i
        @test DS.value(bt[i]) == i
    end

    @test length(bt) == 27

    @test DS.parent(bt[100]) isa DS.NullNode
    @test DS.right_child(bt[100]) isa DS.NullNode
    @test DS.left_child(bt[100]) isa DS.NullNode
    @test DS.level(bt[25]) == 4

    @test repr(bt) == "BinaryTree(\n" *
        "\n\t\t\tTreeNode(15)\n" *
        "\n\t\tTreeNode(7)\n" *
        "\n\t\t\tTreeNode(14)\n" *
        "\n\tTreeNode(3)\n" *
        "\n\t\t\t\tTreeNode(27)\n" *
        "\n\t\t\tTreeNode(13)\n" *
        "\n\t\t\t\tTreeNode(26)\n" *
        "\n\t\tTreeNode(6)\n" *
        "\n\t\t\t\tTreeNode(25)\n" *
        "\n\t\t\tTreeNode(12)\n" *
        "\n\t\t\t\tTreeNode(24)\n" *
        "\nTreeNode(1)\n" *
        "\n\t\t\t\tTreeNode(23)\n" *
        "\n\t\t\tTreeNode(11)\n" *
        "\n\t\t\t\tTreeNode(22)\n" *
        "\n\t\tTreeNode(5)\n" *
        "\n\t\t\t\tTreeNode(21)\n" *
        "\n\t\t\tTreeNode(10)\n" *
        "\n\t\t\t\tTreeNode(20)\n" *
        "\n\tTreeNode(2)\n" *
        "\n\t\t\t\tTreeNode(19)\n" *
        "\n\t\t\tTreeNode(9)\n" *
        "\n\t\t\t\tTreeNode(18)\n" *
        "\n\t\tTreeNode(4)\n" *
        "\n\t\t\t\tTreeNode(17)\n" *
        "\n\t\t\tTreeNode(8)\n" *
        "\n\t\t\t\tTreeNode(16)\n" *
    ")"
end


text:
module sem

using Distributions, LinearAlgebra, Optim, Random,
    NLSolversBase, Statistics, SparseArrays, ModelingToolkit,
    NLopt, FiniteDiff

include("types.jl")
include("observed.jl")
include("helper.jl")
include("diff/diff.jl")
include("diff/ML.jl")
include("diff/ML2.jl")
include("diff/fiml.jl")
include("imply/imply.jl")
include("imply/symbolic.jl")
include("loss/loss.jl")
include("loss/ML.jl")
include("loss/fiml.jl")
include("loss/WLS.jl")
include("diff/WLS.jl")
include("hessian/WLS.jl")
include("loss/definition_variables.jl")
include("loss/regularized.jl")
include("loss/constant.jl")
include("model.jl")
include("multigroup.jl")
#include("fiml.jl")
include("collection.jl")
include("optimizer/optim.jl")
include("optimizer/nlopt.jl")

export  AbstractSem,
            Sem, MGSem, SemFIML, computeloss,
        Imply,
            ImplyCommon, ImplySparse, ImplySymbolic,
            ImplyDense, ImplySymbolicDefinition,
        Loss,
            LossFunction, SemML, SemFIML, SemDefinition, SemLasso, SemRidge,
            SemConstant,
        SemDiff,
            SemFiniteDiff, SemForwardDiff, SemReverseDiff,
            SemAnalyticDiff,
                ∇SemML,
        SemObs,
            SemObsCommon, SemObsMissing,
        sem_fit, get_observed, CollectionSem, remove_all_missing

end # module


text:
module EchoStateNetworks
using Random
using LinearAlgebra

abstract type ESNArch end
struct NoTeacherForcing <: ESNArch end
struct TeacherForcing <: ESNArch end

export EchoStateNetwork
export train!, predict!

mutable struct EchoStateNetwork{T<:AbstractFloat}
    Nr::Integer     # the number of reservoir neurons
    Ni::Integer     # the number of input dimension
    No::Integer     # the number of output dimension
    sparsity::T
    spectral_radius::T
    noise_level::T
    leaking_rate::T
    teacher_forcing::Bool
    Wr::Matrix{T}   # reservoir weights
    Wi::Matrix{T}   # input weights
    Wo::Matrix{T}   # readout weights
    Wf::Matrix{T}   # feedback weights
    state::Vector{T}
    input::Vector{T}
    output::Vector{T}
    rng::AbstractRNG
    activation::Function

    function EchoStateNetwork{T}(;Ni::Integer=1,
        No::Integer=1,
        Nr::Integer=100,
        sparsity::AbstractFloat=0.95,
        spectral_radius::AbstractFloat=0.95,
        noise_level::AbstractFloat=0.001,
        leaking_rate::AbstractFloat=1.0,
        teacher_forcing::Bool=true,
        rng::AbstractRNG=MersenneTwister(rand(UInt32)),
        activation::Function=x->tanh(x),
       ) where T<:AbstractFloat
        esn = new()

        esn.Ni = Ni
        esn.No = No
        esn.Nr = Nr
        esn.sparsity = sparsity
        esn.spectral_radius = spectral_radius
        esn.teacher_forcing = teacher_forcing
        esn.noise_level = noise_level
        esn.leaking_rate = leaking_rate
        esn.rng = rng
        esn.activation = activation

        init_weights!(esn)

        return esn
    end
end

function init_weights!(esn::EchoStateNetwork{T}) where T<:AbstractFloat
    # init reservoir weight matrix
    esn.Wr = rand(esn.rng, T, esn.Nr, esn.Nr) .- T(0.5)
    ## reduce connections based on `sparsity`
    Ns = round(Int,length(esn.Wr)*esn.sparsity)
    for i in [rand(esn.rng, 1:length(esn.Wr)) for i in 1:Ns]
        esn.Wr[i] = zero(T)
    end
    ## rescale the matrix to fit the `spectral radius`
    esn.Wr *= esn.spectral_radius/maximum(abs.(eigvals(esn.Wr)))

    # init input weight matrix
    esn.Wi = rand(esn.rng, T, esn.Nr, esn.Ni+1) .- T(0.5)

    # init feedback weight matrix
    esn.Wf = rand(esn.rng, T, esn.Nr, esn.No) .- T(0.5)

    return
end

function activate(f::Function, x::Array)
    return map(f, x)
end

function update(esn::EchoStateNetwork{T}, state::Vector{T},
                input::Vector{T}, output::Vector{T}, arch::NoTeacherForcing) where T<:AbstractFloat
    return activate(esn.activation, esn.Wi*[one(T);input] + esn.Wr*state)
    + esn.noise_level*(rand(esn.rng, T, esn.Nr)-T(0.5))
end

function update(esn::EchoStateNetwork{T}, state::Vector{T},
                                  input::Vector{T}, output::Vector{T}, arch::TeacherForcing) where T<:AbstractFloat
    return activate(esn.activation, esn.Wi*[one(T);input] + esn.Wr*state + esn.Wf*output)
    + esn.noise_level*(rand(esn.rng, T, esn.Nr)-T(0.5))
end

function reservoir_states(esn::EchoStateNetwork{T}, inputs::Matrix{T}, outputs::Matrix{T}) where T<:AbstractFloat
    Nd = size(inputs, 2)
    states = zeros(T, (esn.Nr, Nd))
    arch = esn.teacher_forcing ? TeacherForcing() : NoTeacherForcing()
    for t = 2:size(inputs,2)
        states[:,t] = (one(T)-esn.leaking_rate)*states[:,t-1] + esn.leaking_rate*update(esn, states[:,t-1], inputs[:,t], outputs[:,t-1], arch)
    end
    return states
end

function train!(esn::EchoStateNetwork{T}, inputs::Matrix{T}, outputs::Matrix{T};
                                  discard::Integer=min(div(size(inputs,2),10), 100), reg::AbstractFloat=1e-8) where T<:AbstractFloat

    @assert(size(inputs, 1) == esn.Ni)
    @assert(size(outputs, 1) == esn.No)
    @assert(size(inputs, 2) == size(outputs, 2))

    states = reservoir_states(esn, inputs, outputs)

    # extended system states
    X = zeros(1+esn.Ni+esn.Nr, size(inputs,2))
    X[1,:] .= one(T)
    X[2:1+esn.Ni,:] = inputs
    X[2+esn.Ni:end,:] = states

    # discard initial transient
    Xe = X[:,discard+1:end]
    tXe = Xe'
    O = outputs[:,discard+1:end]

    # calc output weight matrix
    esn.Wo = O*tXe*pinv(Xe*tXe + reg*Matrix{T}(I, size(Xe,1), size(Xe,1)))

    # store last states
    esn.state = states[:,end]
    esn.input = inputs[:,end]
    esn.output = outputs[:,end]

    return esn.Wo * X
end

function predict!(esn::EchoStateNetwork{T}, inputs::Matrix{T}; cont::Bool=true) where T<:AbstractFloat
    Nd = size(inputs, 2)
    outputs = zeros(T, esn.No, Nd)
    state = zeros(T, esn.Nr)

    if cont
        inputs = hcat(esn.input, inputs)
        outputs = hcat(esn.output, outputs)
        state[:] = esn.state
    else
        inputs = hcat(zeros(T,esn.Ni), inputs)
        outputs = hcat(zeros(T,esn.No), outputs)
    end

    arch = esn.teacher_forcing ? TeacherForcing() : NoTeacherForcing()
    for t = 1:Nd
        state[:] = (one(T)-esn.leaking_rate)*state + esn.leaking_rate*update(esn, state, inputs[:,t+1], outputs[:,t], arch)
        outputs[:,t+1] = esn.Wo*[one(T);inputs[:,t];state]
    end

    return outputs[:,2:end]
end

end # module


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

mutable struct CutCollection{T <: AbstractFloat}
    cuts::Vector{SparseOptimalityCut{T}}
    q::T
    considered::Int
    id::Int

    function CutCollection(::Type{T}, id::Integer = 1) where T <: AbstractFloat
        new{T}(Vector{SparseOptimalityCut{T}}(), zero(T), 0, id)
    end
end

collection_size(collection::CutCollection) = length(collection.cuts)
considered(collection::CutCollection) = collection.considered

function aggregate(collection::CutCollection)
    return aggregate(collection.cuts, collection.id)
end

function renew!(collection::CutCollection{T}, id::Integer) where T <: AbstractFloat
    empty!(collection.cuts)
    collection.q = zero(T)
    collection.considered = 0
    collection.id = id
    return nothing
end

function add_to_collection!(collection::CutCollection, cut::HyperPlane, x::AbstractVector)
    collection.considered += 1
    return nothing
end

function add_to_collection!(collection::CutCollection, cut::HyperPlane{OptimalityCut}, x::AbstractVector)
    push!(collection.cuts, cut)
    collection.q += cut(x)
    collection.considered += 1
    return nothing
end


text:
mutable struct Library <: Polyhedra.Library
    precision::Symbol
    solver

    function Library(precision::Symbol=:float, solver=nothing)
        if !(precision in [:float, :exact])
            error("Invalid precision, it should be :float or :exact")
        end
        new(precision, solver)
    end
end

Polyhedra.similar_library(l::Library, ::Polyhedra.FullDim, ::Type{<:AbstractFloat}) = Library(:float)
Polyhedra.similar_library(l::Library, ::Polyhedra.FullDim, ::Type) = Library(:exact)

const HRepT{T} = LiftedHRepresentation{T, Matrix{T}}
const VRepT{T} = LiftedVRepresentation{T, Matrix{T}}

mutable struct Polyhedron{T} <: Polyhedra.Polyhedron{T}
    hrep::Union{HRepT{T}, Nothing}
    vrep::Union{VRepT{T}, Nothing}
    noredundantinequality::Bool
    noredundantgenerator::Bool
    solver
end

function Polyhedron{T}(hrep::HRepresentation, solver) where {T}
    Polyhedron{T}(convert(HRepT{T}, hrep), nothing, false, false, solver)
end
function Polyhedron{T}(vrep::VRepresentation, solver) where {T}
    Polyhedron{T}(nothing, convert(VRepT{T}, vrep), false, false, solver)
end
#function Polyhedron{T}(rep::Representation, solver) where {T}
#    return Polyhedron{T}(Polyhedra.change_coefficient_type(rep, T), solver)
#end

Polyhedra.FullDim(p::Polyhedron) = Polyhedra.FullDim_rep(p.hrep, p.vrep)
Polyhedra.library(p::Polyhedron{<:AbstractFloat}) = Library(:float, p.solver)
Polyhedra.library(p::Polyhedron) = Library(:exact, p.solver)
Polyhedra.default_solver(p::Polyhedron; T = nothing) = p.solver
Polyhedra.supportssolver(::Type{<:Polyhedron}) = true

function Polyhedra.hvectortype(::Type{Polyhedron{T}}) where {T}
    return Polyhedra.hvectortype(HRepT{T})
end
function Polyhedra.vvectortype(::Type{Polyhedron{T}}) where {T}
    return Polyhedra.vvectortype(VRepT{T})
end

Polyhedra.similar_type(::Type{<:Polyhedron}, ::Polyhedra.FullDim, ::Type{T}) where {T} = Polyhedron{T}

function Polyhedron{T}(d::Polyhedra.FullDim, it::Polyhedra.HIt...; solver=nothing) where {T}
    return Polyhedron{T}(HRepT{T}(d, it...), solver)
end
function Polyhedron{T}(d::Polyhedra.FullDim, it::Polyhedra.VIt...; solver=nothing) where {T}
    return Polyhedron{T}(VRepT{T}(d, it...), solver)
end

function clearfield!(p::Polyhedron)
    p.hrep = nothing
    p.vrep = nothing
    p.noredundantinequality = false
    p.noredundantgenerator = false
end

# Implementation of Polyhedron's mandatory interface
function polytypeforprecision(precision::Symbol)
  if !(precision in (:float, :exact))
    error("precision should be :float or :exact, you gave $precision")
  end
  return precision == :float ? Float64 : Rational{BigInt}
end

function Polyhedra.polyhedron(rep::Representation, lib::Library)
    T = polytypeforprecision(lib.precision)
    return Polyhedron{T}(rep, lib.solver)
end

function Base.copy(p::Polyhedron{T}) where {T}
    hrep = nothing
    if !isnothing(p.hrep)
        hrep = copy(p.hrep)
    end
    vrep = nothing
    if !isnothing(p.vrep)
        vrep = copy(p.vrep)
    end
    return Polyhedron{T}(hrep, vrep, p.noredundantinequality,
                         p.noredundantgenerator, p.solver)
end

function Polyhedra.hrepiscomputed(p::Polyhedron)
    !isnothing(p.hrep)
end
function Polyhedra.computehrep!(p::Polyhedron)
    # vrep(p) could trigger an infinite loop if both vrep and hrep are null
    p.hrep = double_description(p.vrep)
    p.noredundantinequality = true
    return
end
function Polyhedra.hrep(p::Polyhedron)
    if !Polyhedra.hrepiscomputed(p)
        Polyhedra.computehrep!(p)
    end
    return p.hrep
end
function Polyhedra.vrepiscomputed(p::Polyhedron)
    !isnothing(p.vrep)
end
function Polyhedra.computevrep!(p::Polyhedron)
    # hrep(p) could trigger an infinite loop if both vrep and hrep are null
    p.vrep = double_description(p.hrep)
    p.noredundantgenerator = true
    return
end
function Polyhedra.vrep(p::Polyhedron)
    if !Polyhedra.vrepiscomputed(p)
        Polyhedra.computevrep!(p)
    end
    return p.vrep
end

function Polyhedra.sethrep!(p::Polyhedron{T}, h::HRepresentation) where T
    p.hrep = convert(HRepT{T}, h)
end
function Polyhedra.setvrep!(p::Polyhedron{T}, v::VRepresentation) where T
    p.vrep = convert(VRepT{T}, v)
end
function Polyhedra.resethrep!(p::Polyhedron, h::HRepresentation)
    clearfield!(p)
    Polyhedra.sethrep!(p, h)
end
function Polyhedra.resetvrep!(p::Polyhedron, v::VRepresentation)
    clearfield!(p)
    Polyhedra.setvrep!(p, v)
end


text:
# this computes the GMM estimator by SA minimization, for
# one of the 1000 data sets. Can be parallelized using
# threads.
using Econometrics, SolveDSGE, DelimitedFiles, Statistics, LinearAlgebra, ForwardDiff
cd(@__DIR__)
include("DSGEmoments.jl")  # computes errors
include("CKlib.jl")

function main()
cd(@__DIR__)
# using first line for the "true" data set, the second to generate new data
data = readdlm("dsgedata.txt")
#data = dgp(TrueParameters())

# define CUE GMM criterion
include("parameters.jl") # load true parameter values
lb, ub = PriorSupport()
n = size(data,1)
moments = theta -> DSGEmoments(theta, data)
m = theta -> vec(mean(moments(theta),dims=1)) # 1Xg
weight = theta -> inv(cov(sqrt(n)*moments(theta)))
obj = theta -> m(theta)'*weight(theta)*m(theta)
thetastart = (ub+lb)/2.0 # prior mean as start
# estimate by simulated annealing
thetahat, objvalue, converged, details = samin(obj, thetastart, lb, ub; ns = 20, verbosity = 1, rt = 0.5)
# compute the estimated standard errors and CIs
D = ForwardDiff.jacobian(m, vec(thetahat))
W = weight(thetahat)
V = inv(D'*W*D)/n
se = sqrt.(diag(V))
println("estimates, st. error, and limits of 95% CI")
prettyprint([thetahat se thetahat-1.96*se thetahat+1.96*se],["estimate", "std. err.", "CI lower", "CI upper"])
# attempt gradient based (often doesn't work)
println("trying gradient-based")
thetahat, objvalue, converged = fmincon(obj, thetastart, [], [], lb, ub; iterlim=10000)
println("fmincon results. objective fn. value: ", objvalue)
println("parameter values:")
prettyprint(thetahat)
return nothing
end
main()

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
module timeevolution_classical

export timeevolution, timeevolution_stochastic

using QuantumOptics
using ..classical

function timeevolution(T, S::MultimodeSystem, state0::ClassicalState; fout=nothing)
    T = float(T)
    Nparticles = state0.particlenumber
    Ncavitymodes = state0.cavitymodenumber
    mass = S.particles.mass
    function f(t, y::Vector{Float64}, dy::Vector{Float64})
        x, v, α = splitstate(Nparticles, Ncavitymodes, y)
        dx, dv, dα = splitstate(Nparticles, Ncavitymodes, dy)
        for j=1:Nparticles
            dx[j] = v[j]
            dv[j] = 0
        end
        for n=1:Ncavitymodes
            kn = S.modes[n].k
            αn = α[2*n-1] + 1im*α[2*n]
            ηn = S.modes[n].eta
            U0n = S.modes[n].U0
            Δn = S.modes[n].delta
            κn = S.modes[n].kappa
            ηn_eff = 2*ηn*real(αn)
            U0n_eff = 2*U0n*abs2(αn)
            θn = 0.
            Ωn = 0.
            for j=1:Nparticles
                dv[j] -= 1./mass*kn*cos(kn*x[j])*(U0n_eff*sin(kn*x[j]) + ηn_eff)
                θn += sin(kn*x[j])
                Ωn += sin(kn*x[j])^2
            end
            dαn = (1im*(Δn - U0n*Ωn) - κn)*αn - 1im*ηn*θn
            dα[2*n-1] = real(dαn)
            dα[2*n] = imag(dαn)
        end
    end
    if fout==nothing
        t_out = Float64[]
        state_out = ClassicalState[]
        function fout_(t, y::Vector{Float64})
            push!(t_out, t)
            push!(state_out, ClassicalState(Nparticles, Ncavitymodes, deepcopy(y)))
        end
        QuantumOptics.ode_dopri.ode(f, T, state0.data; fout=fout_)
        return t_out, state_out
    else
        return QuantumOptics.ode_dopri.ode(f, T, state0.data; fout=(t,y)->fout(t, ClassicalState(Nparticlemodes, Ncavitymodes, y)))
    end
end

function step(f, κ::Vector{Float64}, t0::Float64, dt::Float64, state0::Vector{Float64})
    T, states = QuantumOptics.ode_dopri.ode(f, [t0, t0+dt], state0, display_initialvalue=false)
    state1 = states[1]
    return state1
end

function ode_stochastic(f, κ::Vector{Float64}, dt::Float64, T::Vector{Float64}, state0::Vector{Float64}, fout)
    t = T[1]
    tfinal = T[end]
    for tfixed = T[2:end]
        while t<tfixed
            tnext = min(t+dt, tfixed)
            state0 = step(f, κ, t, tnext-t, state0)
            t = tnext
        end
        fout(t, state0)
    end
end

function timeevolution_stochastic(T, S::MultimodeSystem, state0::ClassicalState; fout=nothing, dt=1e-3)
    T = float(T)
    Nparticles = state0.particlenumber
    Ncavitymodes = state0.cavitymodenumber
    mass = S.particles.mass
    function f(t, y::Vector{Float64}, dy::Vector{Float64})
        x, v, α = splitstate(Nparticles, Ncavitymodes, y)
        dx, dv, dα = splitstate(Nparticles, Ncavitymodes, dy)
        for j=1:Nparticles
            dx[j] = v[j]
            dv[j] = 0
        end
        for n=1:Ncavitymodes
            kn = S.modes[n].k
            αn = α[2*n-1] + 1im*α[2*n]
            ηn = S.modes[n].eta
            U0n = S.modes[n].U0
            Δn = S.modes[n].delta
            κn = S.modes[n].kappa
            ηn_eff = 2*ηn*real(αn)
            U0n_eff = 2*U0n*abs2(αn)
            θn = 0.
            Ωn = 0.
            for j=1:Nparticles
                dv[j] -= 1./mass*kn*cos(kn*x[j])*(U0n_eff*sin(kn*x[j]) + ηn_eff)
                θn += sin(kn*x[j])
                Ωn += sin(kn*x[j])^2
            end
            dαn = (1im*(Δn - U0n*Ωn) - κn)*αn - 1im*ηn*θn
            dα[2*n-1] = real(dαn)
            dα[2*n] = imag(dαn)
        end
    end
    κ = Float64[S.modes[n].kappa for n=1:Ncavitymodes]
    if fout==nothing
        t_out = Float64[]
        state_out = ClassicalState[]
        function fout_(t, y::Vector{Float64})
            push!(t_out, t)
            push!(state_out, ClassicalState(Nparticles, Ncavitymodes, deepcopy(y)))
        end
        ode_stochastic(f, κ, dt, T, state0.data, fout_)
        return t_out, state_out
    else
        return ode_stochastic(f, κ, dt, T, state0.data, fout=(t,y)->fout(t, MFState(Nparticlemodes, Ncavitymodes, y)))
    end
end

end # module


text:
raw = readdlm("unixdict.txt",String)[:]
inter = intersect(raw,map(reverse,raw)) #find the matching strings/revstrings
res = String[b == 1 && a != reverse(a) && a < reverse(a) ? a : reverse(a) for a in inter, b in 1:2] #create pairs
res = res[res[:,1] .!= res[:,2],:] #get rid of duplicates, palindromes

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using OperatorFlux
using Test
using Flux
using LinearAlgebra

@testset "SpectralKernelOperator - FourierTransform" begin
    for M in [32, 33]
        even = Bool(mod(M, 2)) ? false : true
        # test spectral operator forward evaluation
        trafo = FourierTransform(modes = (8, 7, 5), even = even)
        model = Chain(
            Dense(3, 32), # lifting
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            Dense(32, 4), # unlifting
        )
        x = rand(Float32, 3, M, 16, 27, 17)
        @test size(model(x)) == (4, M, 16, 27, 17)
    end

    # integration test using Flux.jl
    trafo = FourierTransform(modes = (8,))
    model = Chain(
        Dense(1, 32),
        SpectralKernelOperator(trafo, 32 => 32),
        Dense(32, 128),
    )
    params_before = deepcopy(Flux.params(model))
    loss(x, y) = Flux.mse(model(x), y)
    x = rand(Float32, 1, 1024, 5)
    data = [(x, rand(Float32, 128, 1024, 5))]
    Flux.train!(loss, Flux.params(model), data, Flux.ADAM())
    params_after = deepcopy(Flux.params(model))
    # test that something actually happened during
    # optimization with Flux.
    for i in 1:length(params_before)
        @test !all(params_after[i] .== params_before[i])
    end
end

@testset "SpectralKernelOperator - ChebyshevTransform" begin
    for M in [32, 33]
        # test spectral operator forward evaluation
        trafo = ChebyshevTransform(modes = (8,))
        model = Chain(
            Dense(3, 32), # lifting
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            Dense(32, 4), # unlifting
        )
        x = rand(Float32, 3, M, 17)
        @test size(model(x)) == (4, M, 17)
    end

    # integration test using Flux.jl
    trafo = ChebyshevTransform(modes = (8,))
    model = Chain(
        Dense(1, 32),
        SpectralKernelOperator(trafo, 32 => 32),
        Dense(32, 128),
    )
    params_before = deepcopy(Flux.params(model))
    loss(x, y) = Flux.mse(model(x), y)
    x = rand(Float32, 1, 1024, 5)
    data = [(x, rand(Float32, 128, 1024, 5))]
    Flux.train!(loss, Flux.params(model), data, Flux.ADAM())
    params_after = deepcopy(Flux.params(model))
    # test that something actually happened during
    # optimization with Flux.
    for i in 1:length(params_before)
        @test !all(params_after[i] .== params_before[i])
    end
end

@testset "SpectralKernelOperator - LegendreTransform" begin
    for M in [32, 33]
        # test spectral operator forward evaluation
        trafo = LegendreTransform(modes = (8, 8))
        model = Chain(
            Dense(3, 32), # lifting
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            SpectralKernelOperator(trafo, 32 => 32),
            Dense(32, 4), # unlifting
        )
        x = rand(Float32, 3, M, 10, 17)
        @test size(model(x)) == (4, M, 10, 17)
    end

    # TODO!: Broken due to mutating array
    # # integration test using Flux.jl
    # trafo = LegendreTransform(modes = (8, 8))
    # model = Chain(
    #     Dense(1, 32),
    #     SpectralKernelOperator(trafo, 32 => 32),
    #     Dense(32, 128),
    # )
    # params_before = deepcopy(Flux.params(model))
    # loss(x, y) = Flux.mse(model(x), y)
    # x = rand(Float32, 1, 8, 10, 5)
    # data = [(x, rand(Float32, 128, 8, 10, 5))]
    # Flux.train!(loss, Flux.params(model), data, Flux.ADAM())
    # params_after = deepcopy(Flux.params(model))
    # # test that something actually happened during
    # # optimization with Flux.
    # for i in 1:length(params_before)
    #     @test !all(params_after[i] .== params_before[i]) broken = true
    # end
end


text:
__precompile__()

module UnitAliases

using Unitful
import Unitful: m, s, Hz, V, N, F, g, A, C

export HasVoltageUnits,
        HasCapacitanceUnits,
        HasCurrentUnits,
        HasChargeUnits,

        HasTimeUnits,
        HasInverseTimeUnits,

        HasLengthUnits,
        HasForceUnits,
        HasMassUnits,

        HasStiffnessUnits,
        HasVelocityUnits,
        HasAccelerationUnits

const HasVoltageUnits{T,U} = Quantity{T, typeof(0.0V).parameters[2], U}
const HasCapacitanceUnits{T,U} = Quantity{T, typeof(0.0F).parameters[2], U}
const HasCurrentUnits{T,U} = Quantity{T, typeof(0.0A).parameters[2], U}
const HasChargeUnits{T,U} = Quantity{T, typeof(0.0C).parameters[2], U}

const HasTimeUnits{T,U} = Quantity{T, typeof(0.0s).parameters[2], U}
const HasInverseTimeUnits{T,U} = Quantity{T, typeof(0.0Hz).parameters[2], U}

const HasLengthUnits{T,U} = Quantity{T, typeof(0.0m).parameters[2], U}
const HasForceUnits{T,U} = Quantity{T, typeof(0.0N).parameters[2], U}
const HasMassUnits{T,U} = Quantity{T, typeof(0.0g).parameters[2], U}
const HasStiffnessUnits{T,U} = Quantity{T, typeof(0.0N/m).parameters[2], U}
const HasVelocityUnits{T,U} = Quantity{T, typeof(0.0m/s).parameters[2], U}
const HasAccelerationUnits{T,U} = Quantity{T, typeof(0.0m/s^2).parameters[2], U}

end # module


text:
## -- Functions used for MCMC estimation of mineral eruption / deposition ages from observed crystallization age spectra

    """
    ```julia
    check_dist_ll(dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray, tmin::Number, tmax::Number)
    ```
    Return the log-likelihood of a set of mineral ages with means `mu` and
    uncertianty `sigma` being drawn from a given source (i.e., crystallization / closure)
    distribution `dist`, with terms to prevent runaway at low N.

    ### Examples
    ```julia
    mu, sigma = collect(100:0.1:101), 0.01*ones(11)
    ll = check_dist_ll(MeltsVolcanicZirconDistribution, mu, sigma, 100, 101)
    ```
    """
    function check_dist_ll(dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray, tmin::Number, tmax::Number)
        # Define some frequently used variables
        loglikelihood = zero(float(eltype(dist)))
        datarows = length(mu)
        distrows = length(dist)
        dist_yave = nanmean(dist)
        nbins = distrows - 1
        dt = abs(tmax-tmin)
        # Cycle through each datum in dataset
        @inbounds for j=1:datarows
            # Find equivalent index position of mu[j] in the `dist` array
            ix = (mu[j] - tmin) / dt * nbins + 1
            # If possible, prevent aliasing problems by interpolation
            if (sigma[j] < dt/nbins) && ix > 1 && ix < distrows
                # Interpolate corresponding distribution value
                f = floor(Int,ix)
                likelihood = (dist[f+1]*(ix-f) + dist[f]*(1-(ix-f))) / (dt * dist_yave)
                # Otherwise, sum contributions from Gaussians at each point in distribution
            else
                likelihood = zero(float(eltype(dist)))
                @inbounds @simd for i=1:distrows
                    distx = tmin + dt*(i-1)/nbins # time-position of distribution point
                    # Likelihood curve follows a Gaussian PDF. Note: dt cancels
                    likelihood += dist[i] / (dist_yave * distrows * sigma[j] * sqrt(2*pi)) *
                            exp( - (distx-mu[j])*(distx-mu[j]) / (2*sigma[j]*sigma[j]) )
                end
            end
            loglikelihood += log(likelihood)
        end
        # Calculate a weighted mean and examine our MSWD
        (wm, wsigma, mswd) = awmean(mu, sigma)
        if datarows == 1 || mswd < 1
            Zf = 1.0
        elseif mswd*sqrt(datarows) > 1000
            Zf = 0.0
        else
            f = datarows - 1
            # Height of MSWD distribution relative to height at MSWD = 1 (see Wendt and Carl, 1991, Chemical geology)
            Zf = exp((f/2-1)*log(mswd) - f/2*(mswd-1))
        end
        # To prevent instability / runaway of the MCMC for small datasets (low N),
        # favor the weighted mean interpretation at high Zf (MSWD close to 1) and
        # the youngest-zircon interpretation at low Zf (MSWD far from one). The
        # penalty factors used here are determined by training against synthetic datasets.
        # In other words, these are just context-dependent prior distributions on tmax and tmin
        return loglikelihood - (2/log(1+datarows)) * (              # Scaling factor that decreases with log number of data points (i.e., no penalty at high N)
        log((abs(tmin - wm)+wsigma)/wsigma)*Zf +                    # Penalty for proposing tmin too far from the weighted mean at low MSWD (High Zf)
        log((abs(tmax - wm)+wsigma)/wsigma)*Zf +                    # Penalty for proposing tmax too far from the weighted mean at low MSWD (High Zf)
        log((abs(tmin - mu[1])+sigma[1])/sigma[1])*(1-Zf) +         # Penalty for proposing tmin too far from youngest zircon at high MSWD (low Zf)
        log((abs(tmax - mu[end])+sigma[end])/sigma[end])*(1-Zf) )   # Penalty for proposing tmax too far from oldest zircon at high MSWD (low Zf)
    end

    """
    ```julia
    metropolis_minmax!(tminDist, tmaxDist, llDist, acceptanceDist, nsteps::Int, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)
    ```
    In-place (non-allocating) version of `metropolis_minmax`, filling existing arrays

    Run a Metropolis sampler to estimate the extrema of a finite-range source
    distribution `dist` using samples drawn from that distribution -- e.g.,
    estimate zircon saturation and eruption ages from a distribution of zircon
    crystallization ages.

    ### Examples
    ```julia
    metropolis_minmax!(tmindist, tmaxdist, lldist, acceptancedist, 2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
    ```
    """
    function metropolis_minmax!(tminDist::AbstractArray, tmaxDist::AbstractArray, llDist::AbstractArray, acceptanceDist::AbstractArray, nsteps::Int, dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
        # standard deviation of the proposal function is stepfactor * last step; this is tuned to optimize accetance probability at 50%
        stepfactor = 2.9
        # Sort the dataset from youngest to oldest
        sI = sortperm(mu)
        mu_sorted = mu[sI] # Sort means
        sigma_sorted = sigma[sI] # Sort uncertainty
        # These quantities will be used more than once
        datarows = length(mu_sorted)
        tmin_obs = nanminimum(mu_sorted)
        tmax_obs = nanmaximum(mu_sorted)
        # Step sigma for Gaussian proposal distributions
        dt = tmax_obs - tmin_obs + sigma_sorted[1] + sigma_sorted[end]
        tmin_step = dt / datarows
        tmax_step = dt / datarows
        # Use oldest and youngest zircons for initial proposal
        tmin = tmin_obs - sigma_sorted[1]
        tmax = tmax_obs + sigma_sorted[end]
        tminₚ = tmin
        tmaxₚ = tmax
        # Log likelihood of initial proposal
        ll = check_dist_ll(dist, mu_sorted, sigma_sorted, tmin, tmax)
        llₚ = ll
        # Burnin
        for i=1:nsteps
            tminₚ = tmin
            tmaxₚ = tmax
            # Adjust either upper or lower bound
            if rand()<0.5
                tminₚ += tmin_step*randn()
            else
                tmaxₚ += tmax_step*randn()
            end
            # Flip bounds if reversed
            if (tminₚ>tmaxₚ)
                tminₚ, tmaxₚ = tmaxₚ, tminₚ
            end
            # Calculate log likelihood for new proposal
            llₚ = check_dist_ll(dist, mu_sorted, sigma_sorted, tminₚ, tmaxₚ)
            # Decide to accept or reject the proposal
            if log(rand()) < (llₚ-ll)
                if tminₚ != tmin
                    tmin_step = abs(tminₚ-tmin)*stepfactor
                end
                if tmaxₚ != tmax
                    tmax_step = abs(tmaxₚ-tmax)*stepfactor
                end

                ll = llₚ
                tmin = tminₚ
                tmax = tmaxₚ
            end
        end
        # Step through each of the N steps in the Markov chain
        @inbounds for i=1:nsteps
            tminₚ = tmin
            tmaxₚ = tmax
            # Adjust either upper or lower bound
            if rand()<0.5
                tminₚ += tmin_step*randn()
            else
                tmaxₚ += tmax_step*randn()
            end
            # Flip bounds if reversed
            if (tminₚ>tmaxₚ)
                tminₚ, tmaxₚ = tmaxₚ, tminₚ
            end
            # Calculate log likelihood for new proposal
            llₚ = check_dist_ll(dist, mu_sorted, sigma_sorted, tminₚ, tmaxₚ)
            # Decide to accept or reject the proposal
            if log(rand()) < (llₚ-ll)
                if tminₚ != tmin
                    tmin_step = abs(tminₚ-tmin)*stepfactor
                end
                if tmaxₚ != tmax
                    tmax_step = abs(tmaxₚ-tmax)*stepfactor
                end

                ll = llₚ
                tmin = tminₚ
                tmax = tmaxₚ
                acceptanceDist[i]=true
            end
            tminDist[i] = tmin
            tmaxDist[i] = tmax
            llDist[i] = ll
        end
        return tminDist, tmaxDist, llDist, acceptanceDist
    end

    """
    ```julia
    metropolis_minmax(nsteps::Int, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)
    ```
    Run a Metropolis sampler to estimate the extrema of a finite-range source
    distribution `dist` using samples drawn from that distribution -- e.g.,
    estimate zircon saturation and eruption ages from a distribution of zircon
    crystallization ages.

    ### Examples
    ```julia
    tmindist, tmaxdist, lldist, acceptancedist = metropolis_minmax(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
    ```
    """
    function metropolis_minmax(nsteps::Int, dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
        # Allocate ouput arrays
        acceptanceDist = falses(nsteps)
        llDist = Array{float(eltype(dist))}(undef,nsteps)
        tmaxDist = Array{float(eltype(mu))}(undef,nsteps)
        tminDist = Array{float(eltype(mu))}(undef,nsteps)
        # Run metropolis sampler
        return metropolis_minmax!(tminDist, tmaxDist, llDist, acceptanceDist, nsteps, dist, mu, sigma; burnin=burnin)
    end

    """
    ```julia
    metropolis_min!(tminDist::Array, nsteps::Int, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)
    ```
    In-place (non-allocating) version of `metropolis_min`, fills existing array `tminDist`.

    Run a Metropolis sampler to estimate the minimum of a finite-range source
    distribution `dist` using samples drawn from that distribution -- e.g., estimate
    zircon eruption ages from a distribution of zircon crystallization ages.

    ### Examples
    ```julia
    metropolis_min!(tminDist, 2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
    ```
    """
    function metropolis_min!(tminDist::AbstractArray, nsteps::Int, dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
        # standard deviation of the proposal function is stepfactor * last step; this is tuned to optimize accetance probability at 50%
        stepfactor = 2.9
        # Sort the dataset from youngest to oldest
        sI = sortperm(mu)
        mu_sorted = mu[sI] # Sort means
        sigma_sorted = sigma[sI] # Sort uncertainty
        # These quantities will be used more than once
        datarows = length(mu_sorted)
        tmin_obs = nanminimum(mu_sorted)
        tmax_obs = nanmaximum(mu_sorted)
        # Step sigma for Gaussian proposal distributions
        dt = tmax_obs - tmin_obs + sigma_sorted[1] + sigma_sorted[end]
        tmin_step = dt / datarows
        tmax_step = dt / datarows
        # Use oldest and youngest zircons for initial proposal
        tmin = tmin_obs - sigma_sorted[1]
        tmax = tmax_obs + sigma_sorted[end]
        tminₚ = tmin
        tmaxₚ = tmax
        # Log likelihood of initial proposal
        ll = check_dist_ll(dist, mu_sorted, sigma_sorted, tmin, tmax)
        llₚ = ll
        # Burnin
        for i=1:burnin
            tminₚ = tmin
            tmaxₚ = tmax
            # Adjust either upper or lower bound
            if rand()<0.5
                tminₚ += tmin_step*randn()
            else
                tmaxₚ += tmax_step*randn()
            end
            # Flip bounds if reversed
            if (tminₚ>tmaxₚ)
                tminₚ, tmaxₚ = tmaxₚ, tminₚ
            end
            # Calculate log likelihood for new proposal
            llₚ = check_dist_ll(dist, mu_sorted, sigma_sorted, tminₚ, tmaxₚ)
            # Decide to accept or reject the proposal
            if log(rand()) < (llₚ-ll)
                if tminₚ != tmin
                    tmin_step = abs(tminₚ-tmin)*stepfactor
                end
                if tmaxₚ != tmax
                    tmax_step = abs(tmaxₚ-tmax)*stepfactor
                end

                ll = llₚ
                tmin = tminₚ
                tmax = tmaxₚ
            end
        end
        # Step through each of the N steps in the Markov chain
        @inbounds for i=1:nsteps
            tminₚ = tmin
            tmaxₚ = tmax
            # Adjust either upper or lower bound
            if rand()<0.5
                tminₚ += tmin_step*randn()
            else
                tmaxₚ += tmax_step*randn()
            end
            # Flip bounds if reversed
            if (tminₚ>tmaxₚ)
                tminₚ, tmaxₚ = tmaxₚ, tminₚ
            end
            # Calculate log likelihood for new proposal
            llₚ = check_dist_ll(dist, mu_sorted, sigma_sorted, tminₚ, tmaxₚ)
            # Decide to accept or reject the proposal
            if log(rand()) < (llₚ-ll)
                if tminₚ != tmin
                    tmin_step = abs(tminₚ-tmin)*stepfactor
                end
                if tmaxₚ != tmax
                    tmax_step = abs(tmaxₚ-tmax)*stepfactor
                end

                ll = llₚ
                tmin = tminₚ
                tmax = tmaxₚ
            end
            tminDist[i] = tmin
        end
        return tminDist
    end

    """
    ```julia
    tminDist = metropolis_min(nsteps::Int, dist::AbstractArray, data::AbstractArray, uncert::AbstractArray; burnin::Integer=0)
    ```
    Run a Metropolis sampler to estimate the minimum of a finite-range source
    distribution `dist` using samples drawn from that distribution -- e.g., estimate
    zircon eruption ages from a distribution of zircon crystallization ages.

    ### Examples
    ```julia
    tmindist = metropolis_min(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
    ```
    """
    function metropolis_min(nsteps::Int, dist::AbstractArray, mu::AbstractArray, sigma::AbstractArray; burnin::Integer=0)
        # Allocate ouput array
        tminDist = Array{float(eltype(mu))}(undef,nsteps)
        # Run Metropolis sampler
        return metropolis_min!(tminDist, nsteps, dist, mu, sigma; burnin=burnin)
    end

## --- Some useful distributions

    #[1.0, 1.0] -- maximum entropy:
    const UniformDistribution = ones(100)

    #[2.0,1.0,0.0]
    const TriangularDistribution = collect(range(2,0,length=100))

    # Best for survivorship processes, apparently including Ar-Ar partial degassing / entrainment / inheritance
    const ExponentialDistribution = [7.94734, 7.38851, 6.86897, 6.38597, 5.93693, 5.51946, 5.13135, 4.77053, 4.43508, 4.12322, 3.83329, 3.56375, 3.31315, 3.08018, 2.8636, 2.66224, 2.47504, 2.301, 2.1392, 1.98878, 1.84893, 1.71892, 1.59805, 1.48568, 1.38122, 1.28409, 1.1938, 1.10986, 1.03181, 0.959261, 0.891808, 0.829099, 0.7708, 0.7166, 0.666211, 0.619365, 0.575813, 0.535324, 0.497681, 0.462686, 0.430152, 0.399905, 0.371785, 0.345642, 0.321338, 0.298742, 0.277736, 0.258206, 0.24005, 0.22317, 0.207478, 0.192889, 0.179325, 0.166716, 0.154993, 0.144094, 0.133962, 0.124542, 0.115785, 0.107643, 0.100074, 0.0930371, 0.0864951, 0.080413, 0.0747586, 0.0695018, 0.0646147, 0.0600712, 0.0558472, 0.0519202, 0.0482693, 0.0448752, 0.0417197, 0.0387861, 0.0360588, 0.0335233, 0.031166, 0.0289745, 0.0269371, 0.025043, 0.023282, 0.0216449, 0.0201229, 0.0187079, 0.0173925, 0.0161695, 0.0150325, 0.0139755, 0.0129928, 0.0120791, 0.0112298, 0.0104401, 0.00970602, 0.00902352, 0.00838902, 0.00779913, 0.00725072, 0.00674087, 0.00626688, 0.00582621, 0.00541653, 0.00503566, 0.00468156, 0.00435237, 0.00404633, 0.0037618, 0.00349728, 0.00325137, 0.00302274, 0.00281019]

    const HalfNormalDistribution =
    [2.65307, 2.6516, 2.64718, 2.63984, 2.6296, 2.61648, 2.60054, 2.58182, 2.5604, 2.53633, 2.5097, 2.48059, 2.4491, 2.41532, 2.37936, 2.34133, 2.30135, 2.25954, 2.21603, 2.17095, 2.12442, 2.07657, 2.02755, 1.97749, 1.92653, 1.87479, 1.82242, 1.76954, 1.71629, 1.66279, 1.60917, 1.55555, 1.50205, 1.44878, 1.39584, 1.34335, 1.29139, 1.24006, 1.18946, 1.13965, 1.09071, 1.04272, 0.995728, 0.9498, 0.904985, 0.861327, 0.818864, 0.777631, 0.737653, 0.698954, 0.66155, 0.625452, 0.590667, 0.557197, 0.52504, 0.494189, 0.464635, 0.436363, 0.409356, 0.383594, 0.359054, 0.335711, 0.313537, 0.292503,0.272576, 0.253725, 0.235916, 0.219112, 0.20328, 0.188382, 0.174383, 0.161244, 0.14893, 0.137403, 0.126628, 0.116568,0.107188, 0.0984536, 0.0903304, 0.0827854, 0.0757864, 0.069302, 0.0633021, 0.0577574, 0.0526399, 0.0479225, 0.0435794, 0.0395859, 0.0359185, 0.0325546, 0.029473, 0.0266534, 0.0240769, 0.0217252, 0.0195815, 0.0176297, 0.0158548, 0.0142428, 0.0127805, 0.0114555, 0.0102566]

    const TruncatedNormalDistribution =
    [1.15708, 1.2038, 1.25037, 1.29662, 1.34239, 1.38751, 1.43181, 1.47511, 1.51724, 1.55802, 1.5973, 1.63489, 1.67064, 1.70439, 1.73598, 1.76527, 1.79213, 1.81643, 1.83805, 1.8569, 1.87289, 1.88593, 1.89596, 1.90294, 1.90682, 1.9076, 1.90527, 1.89983, 1.89133, 1.87978, 1.86526, 1.84783, 1.82758, 1.80461, 1.77901, 1.75092, 1.72046, 1.68777, 1.65301, 1.61632, 1.57786, 1.53781, 1.49633, 1.45359, 1.40977, 1.36504, 1.31958, 1.27355, 1.22712, 1.18045, 1.13371, 1.08704, 1.04059, 0.9945, 0.948905, 0.903921, 0.859665, 0.816241, 0.773749, 0.732271, 0.691887, 0.652662, 0.614655, 0.577921, 0.542496, 0.508411, 0.475689, 0.444348, 0.414395, 0.38583, 0.358649, 0.332839, 0.308382, 0.285256, 0.263435, 0.242885, 0.223573, 0.205461, 0.188509, 0.172673, 0.157909, 0.144172, 0.131416, 0.119592, 0.108655, 0.0985572, 0.0892521, 0.0806935, 0.0728368, 0.0656377, 0.0590535, 0.0530436, 0.0475671, 0.0425867, 0.0380655, 0.0339688, 0.0302636, 0.0269185, 0.0239041, 0.0211926]

    const MeltsZirconDistribution =
    [0.282361, 0.28919, 0.296019, 0.302849, 0.30968, 0.316567, 0.323614, 0.33064, 0.337727, 0.344848, 0.352146, 0.359642,0.367482, 0.375622, 0.384052, 0.392828, 0.401968, 0.411518, 0.421442, 0.43171, 0.44213, 0.45295, 0.464036, 0.47539, 0.486938, 0.498644, 0.51075, 0.523026, 0.535688, 0.548764, 0.562124, 0.575927, 0.590363, 0.604879, 0.620415, 0.636022, 0.652333, 0.669112, 0.686441, 0.704341, 0.72283, 0.742036, 0.761964, 0.782541, 0.803718, 0.825707, 0.848386, 0.871895,0.896139, 0.920462, 0.946071, 0.972964, 0.999905, 1.02776, 1.05664, 1.08637, 1.11731, 1.14919, 1.18202, 1.21582, 1.24956, 1.28342, 1.31828, 1.35427, 1.39153, 1.43006, 1.46879, 1.50812, 1.5477, 1.58888, 1.63149, 1.6748, 1.71724, 1.76126, 1.80668, 1.85101, 1.89546, 1.94144, 1.98379, 2.02785, 2.06738, 2.10669, 2.1377, 2.16306, 2.17843, 2.17924, 2.16073, 2.11744, 2.04444, 1.93323, 1.7923, 1.62527, 1.44425, 1.25401, 1.0528, 0.843628, 0.632687, 0.421876, 0.211064, 0.000252985]

    const MeltsVolcanicZirconDistribution =
    [0.54933, 0.556409, 0.563488, 0.570567, 0.577653, 0.584759, 0.591912, 0.599251, 0.606793, 0.614519, 0.622425, 0.630421, 0.63852, 0.646681, 0.654972, 0.663533, 0.672274, 0.681233, 0.690399, 0.699787, 0.709334, 0.719174, 0.729157, 0.739423, 0.749935, 0.760644, 0.771726, 0.782974, 0.794507, 0.806296, 0.818297, 0.830517, 0.842957, 0.855411, 0.866744, 0.878127, 0.889792, 0.901792, 0.914121, 0.926689, 0.939557, 0.952834, 0.966425, 0.980333, 0.994521, 1.00914, 1.02403, 1.03928, 1.05487, 1.0705, 1.08587, 1.10097, 1.11608, 1.13153, 1.1474, 1.16353, 1.18025, 1.19743, 1.21504, 1.23312, 1.25034, 1.26711, 1.28441, 1.30212, 1.32024, 1.33892, 1.35769, 1.37491, 1.3923, 1.41046, 1.42924, 1.44775, 1.46432, 1.48171, 1.49969, 1.51516, 1.53001, 1.54571, 1.5566, 1.56814, 1.57522, 1.58168, 1.58206, 1.57869, 1.56907, 1.55064, 1.51982, 1.4737, 1.40944, 1.32047, 1.21218, 1.09157, 0.965488, 0.834108, 0.697552, 0.558304, 0.418827, 0.279262, 0.139695, 0.000127237]

    const EllisDistribution =
    [6.80942, 5.36792, 4.45867, 3.83457, 3.28267, 2.77244, 2.33403, 1.98717, 1.72219, 1.52427, 1.37856, 1.27023, 1.18445,1.10697, 1.03176, 0.958823, 0.888329, 0.820435, 0.755302, 0.693089, 0.633956, 0.578063, 0.525569, 0.476635, 0.431421,0.390065, 0.35249, 0.31849, 0.287855, 0.260378, 0.235851, 0.214064, 0.194811, 0.177882, 0.163069, 0.150163, 0.138957,0.129243, 0.120811, 0.113453, 0.106962, 0.101129, 0.095775, 0.0908504, 0.0863393, 0.0822255, 0.0784932, 0.0751262, 0.0721086, 0.0694244, 0.0670576]

    ## --- End of File


text:
@test degree(0) == 1
@test degree(1) == 2
@test detailindex(0,1) == 2

level = 5
@test degree(level) + 1 == first(detailrange(level))
@test detailindex(level, 3) == detailrange(level)[3]

level = 7
signal = zeros(2^level)

@test scales(2^level) == level
@test degree(dyadiclevelmax(2^level)) == 2^(level - 1)
@test scales(signal) == level
@test dyadiclevelmax(signal) == level - 1

dyadiclevel = 8
signalsize = 512
signal = zeros(signalsize)

@test todyadiclevel(signalsize, 1) == dyadiclevel
@test iszero(todyadiclevel(signalsize, dyadiclevel + 1))
@test todyadiclevel(signal, 1) == dyadiclevel
@test iszero(todyadiclevel(signal, dyadiclevel + 1))

signalsize = 64
level = 2
offset = 3

@test levelmax(signalsize) == scales(signalsize)
@test levelmax(signalsize) == overlaplevelmax(signalsize)
@test detailindex(signalsize, level, offset) == detailindex(todyadiclevel(signalsize, level), offset)
@test detailrange(signalsize, level) == detailrange(todyadiclevel(signalsize, level))
@test degree(signalsize, level) == degree(todyadiclevel(signalsize, level))

@test isdyadic(rand(4))
@test isdyadic(4)
@test isdyadic(rand(32,32))
@test !isdyadic(rand(6))
@test !isdyadic(6)


text:
using Ogg
using Compat.Test

@testset "Ogg synthesis/analysis" begin
    # Let's start with building our own Ogg structure, writing it out to an IOBuf,
    # then loading it back in again and checking everything about it we can think of

    # We are going to build three streams, each with 10 packets
    num_packets = 10
    stream_ids = Cint[1, 2, 3]
    packets = Dict{Clong,Vector{Vector{UInt8}}}()
    granulepos = Dict{Int64,Vector{Int64}}()
    for serial in stream_ids
        # The packets are all of different size
        packets[serial] = Vector{UInt8}[UInt8.(mod.(collect(1:100*x), 256)) for x in 1:num_packets]

        # Each packet will have a monotonically increasing granulepos, except for
        # the first two packets which are our "header" packets with granulepos == 0
        granulepos[serial] = Int64[0, 0, [20*x for x in 1:(num_packets - 2)]...]
    end

    # Now we write these packets out to an IOBuffer
    ogg_file = IOBuffer()
    save(ogg_file, packets, granulepos)

    # Rewind to the beginning of this IOBuffer and load the packets back in
    seekstart(ogg_file)

    # This is the verbose way to load an Ogg file that leaves the decoder object
    # around for us to inspect; the load() function returns only `dec.packets`
    dec = OggDecoder()
    Ogg.decode_all_pages(dec, ogg_file)
    Ogg.decode_all_packets(dec, ogg_file)

    # # Are the serial numbers the same?
    @test sort(collect(keys(dec.packets))) == sort(stream_ids)

    # Are the number of packets the same?
    for serial in stream_ids
        @test length(dec.packets[serial]) == length(packets[serial])
    end

    # Are the contents of the packets the same?
    for serial in stream_ids
        for packet_idx in 1:length(dec.packets[serial])
            @test dec.packets[serial][packet_idx] == packets[serial][packet_idx]
        end
    end

    # Let's dig deeper; let's ensure that the first two pages had length equal to
    # our first two packets, proving that our header packets had their own pages:
    for serial in stream_ids
        @test dec.pages[serial][1].body_len == length(packets[serial][1])
        @test dec.pages[serial][2].body_len == length(packets[serial][2])
    end
end

# Next, let's load a known ogg stream and ensure that it's exactly as we expect
@testset "Known .ogg decoding" begin
    testdir = dirname(@__FILE__)
    ogg_packets = load(joinpath(testdir, "zero.ogg"))

    # There is only one stream, and we know its serial number
    @test collect(keys(ogg_packets)) == [1238561138]
    serial = first(keys(ogg_packets))

    # There are four packets, the first starts with \x7fFLAC
    @test length(ogg_packets[serial]) == 4
    @test String(copy(ogg_packets[serial][1][2:5])) == "FLAC"

    # The lengths of the packets are:
    @test [length(x) for x in ogg_packets[serial]] == [51, 55, 13, 0]
end


text:
using Revise
using VCFTools
using MendelImpute
using GeneticVariation
using Random
using StatsBase
using CodecZlib
using ProgressMeter
using JLD2, FileIO, JLSO
using BenchmarkTools
using GroupSlices
using LinearAlgebra

BLAS.set_num_threads(1)
Threads.nthreads() != 10 && error("not 10 threads!")

Random.seed!(2020)
function run()
    chr = 20
    maf = 0.01
    max_d = 1000
    X_complete = convert_gt(UInt8, "target.chr$chr.full.vcf.gz")
    n, p = size(X_complete)
    for overlap in [0.25, 0.5, 1.0]
	tgtfile = "target.chr$chr.typedOnly.masked.vcf.gz"
	reffile = "ref.chr$chr.maxd$max_d.overlap$overlap.excludeTarget.jlso"
	outfile = "mendel.imputed.chr$chr.vcf.gz"

	ph = phase(tgtfile, reffile, outfile=outfile, impute=true, max_d=max_d, phase=true)
        X_mendel = convert_gt(UInt8, outfile)
        println("overlap $overlap error = $(sum(X_mendel .!= X_complete) / n / p) \n")
        flush(stdout)
    end
end
run()

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
@latexrecipe function f(eqs::Vector{ModelingToolkit.Equation})
    # Set default option values.
    env --> :align
    cdot --> false

    # Convert both the left and right hand side to expressions of basic types
    # that latexify can deal with.

    rhs = getfield.(eqs, :rhs)
    rhs = convert.(Expr, rhs)
    rhs = [postwalk(x -> x isa ModelingToolkit.Constant ? x.value : x, eq) for eq in rhs]
    rhs = [postwalk(x -> x isa Expr && length(x.args) == 1 ? x.args[1] : x, eq) for eq in rhs]
    rhs = [postwalk(x -> x isa Expr && x.args[1] == :derivative && length(x.args[2].args) == 2 ? :($(Symbol(:d, x.args[2]))/($(Symbol(:d, x.args[2].args[2])))) : x, eq) for eq in rhs]
    rhs = [postwalk(x -> x isa Expr && x.args[1] == :derivative ? "\\frac{d\\left($(Latexify.latexraw(x.args[2]))\\right)}{d$(Latexify.latexraw(x.args[3]))}" : x, eq) for eq in rhs]

    lhs = getfield.(eqs, :lhs)
    lhs = convert.(Expr, lhs)
    lhs = [postwalk(x -> x isa ModelingToolkit.Constant ? x.value : x, eq) for eq in lhs]
    lhs = [postwalk(x -> x isa Expr && length(x.args) == 1 ? x.args[1] : x, eq) for eq in lhs]
    lhs = [postwalk(x -> x isa Expr && x.args[1] == :derivative && length(x.args[2].args) == 2 ? :($(Symbol(:d, x.args[2]))/($(Symbol(:d, x.args[2].args[2])))) : x, eq) for eq in lhs]
    lhs = [postwalk(x -> x isa Expr && x.args[1] == :derivative ? "\\frac{d\\left($(Latexify.latexraw(x.args[2]))\\right)}{d$(Latexify.latexraw(x.args[3]))}" : x, eq) for eq in lhs]

    return lhs, rhs
end

@latexrecipe function f(sys::ModelingToolkit.AbstractSystem)
    return latexify(equations(sys))
end


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

@doc raw"""CustomResourceSubresourceStatus defines how to serve the status subresource for CustomResources. Status is represented by the `.status` JSON path inside of a CustomResource. When set, * exposes a /status subresource for the custom resource * PUT requests to the /status subresource take a custom resource object, and ignore changes to anything except the status stanza * PUT/POST/PATCH requests to the custom resource ignore changes to the status stanza

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus(;
    )

"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus <: SwaggerModel

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus(;)
        o = new()
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus = Dict{Symbol,Symbol}()
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresourceStatus)
    true
end

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