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

# methods related to array printing

# Printing a value requires to take into account the :typeinfo property
# from the IO context; this property encodes (as a type) the type information
# that is supposed to have already been displayed concerning this value,
# so that redundancy can be avoided. For example, when printing an array of
# `Float16` values, the header "Float16" will be printed, and the values
# can simply be printed with the decimal representations:
# show(Float16(1)) -> "Float16(1.0)"
# show([Float16(1)]) -> "Float16[1.0]" (instead of "Float16[Float16(1.0)]")
# Similarly:
# show([[Float16(1)]]) -> "Array{Float16}[[1.0]]" (instead of "Array{Float16}[Float16[1.0]]")
#
# The array printing methods here can be grouped into two categories (and are annotated as such):
# 1) "typeinfo aware" : these are "API boundaries" functions, which will read the typeinfo
#    property from the context, and pass down to their value an updated property
#    according to its eltype; at each layer of nesting, only one "typeinfo aware"
#    function must be called;
# 2) "typeinfo agnostic": these are helper functions used by the first category; hence
#    they don't manipulate the typeinfo property, and let the printing routines
#    for their elements read directly the property set by their callers
#
# Non-annotated functions are even lower level (e.g. print_matrix_row), so they fall
# by default into category 2.
#
# The basic organization of this file is
# 1) printing with `display`
# 2) printing with `show`
# 3) Logic for displaying type information

## printing with `display`

"""
Unexported convenience function used in body of `replace_in_print_matrix`
methods. By default returns a string of the same width as original with a
centered cdot, used in printing of structural zeros of structured matrices.
Accept keyword args `c` for alternate single character marker.
"""
function replace_with_centered_mark(s::AbstractString;c::AbstractChar = '⋅')
    N = length(s)
    return join(setindex!([" " for i=1:N],string(c),ceil(Int,N/2)))
end

const undef_ref_alignment = (3,3)

"""
`alignment(io, X, rows, cols, cols_if_complete, cols_otherwise, sep)` returns the
alignment for specified parts of array `X`, returning the (left,right) info.
It will look in X's `rows`, `cols` (both lists of indices)
and figure out what's needed to be fully aligned, for example looking all
the way down a column and finding out the maximum size of each element.
Parameter `sep::Integer` is number of spaces to put between elements.
`cols_if_complete` and `cols_otherwise` indicate screen width to use.
Alignment is reported as a vector of (left,right) tuples, one for each
column going across the screen.
"""
function alignment(io::IO, @nospecialize(X::AbstractVecOrMat),
        rows::AbstractVector{T}, cols::AbstractVector{V},
        cols_if_complete::Integer, cols_otherwise::Integer, sep::Integer,
        #= `size(X) may not infer, set this in caller =# ncols::Integer=size(X, 2)) where {T,V}
    a = Tuple{T, V}[]
    for j in cols # need to go down each column one at a time
        l = r = 0
        for i in rows # plumb down and see what largest element sizes are
            if isassigned(X,i,j)
                aij = alignment(io, X[i,j])::Tuple{Int,Int}
            else
                aij = undef_ref_alignment
            end
            l = max(l, aij[1]) # left characters
            r = max(r, aij[2]) # right characters
        end
        push!(a, (l, r)) # one tuple per column of X, pruned to screen width
        if length(a) > 1 && sum(map(sum,a)) + sep*length(a) >= cols_if_complete
            pop!(a) # remove this latest tuple if we're already beyond screen width
            break
        end
    end
    if 1 < length(a) < ncols
        while sum(map(sum,a)) + sep*length(a) >= cols_otherwise
            pop!(a)
        end
    end
    return a
end

"""
`print_matrix_row(io, X, A, i, cols, sep)` produces the aligned output for
a single matrix row X[i, cols] where the desired list of columns is given.
The corresponding alignment A is used, and the separation between elements
is specified as string sep.
`print_matrix_row` will also respect compact output for elements.
"""
function print_matrix_row(io::IO,
        @nospecialize(X::AbstractVecOrMat), A::Vector,
        i::Integer, cols::AbstractVector, sep::AbstractString,
        #= `axes(X)` may not infer, set this in caller =# idxlast::Integer=last(axes(X, 2)))
    for (k, j) = enumerate(cols)
        k > length(A) && break
        if isassigned(X,Int(i),Int(j)) # isassigned accepts only `Int` indices
            x = X[i,j]
            a = alignment(io, x)::Tuple{Int,Int}

            # First try 3-arg show
            sx = sprint(show, "text/plain", x, context=io, sizehint=0)

            # If the output contains line breaks, try 2-arg show instead.
            if occursin('\n', sx)
                sx = sprint(show, x, context=io, sizehint=0)
            end
        else
            a = undef_ref_alignment
            sx = undef_ref_str
        end
        l = repeat(" ", A[k][1]-a[1]) # pad on left and right as needed
        r = j == idxlast ? "" : repeat(" ", A[k][2]-a[2])
        prettysx = replace_in_print_matrix(X,i,j,sx)
        print(io, l, prettysx, r)
        if k < length(A); print(io, sep); end
    end
end

"""
`print_matrix_vdots` is used to show a series of vertical ellipsis instead
of a bunch of rows for long matrices. Not only is the string vdots shown
but it also repeated every M elements if desired.
"""
function print_matrix_vdots(io::IO, vdots::AbstractString,
                            A::Vector, sep::AbstractString, M::Integer, m::Integer,
                            pad_right::Bool = true)
    for k = 1:length(A)
        w = A[k][1] + A[k][2]
        if k % M == m
            l = repeat(" ", max(0, A[k][1]-length(vdots)))
            r = k == length(A) && !pad_right ?
                "" :
                repeat(" ", max(0, w-length(vdots)-length(l)))
            print(io, l, vdots, r)
        else
            (k != length(A) || pad_right) && print(io, repeat(" ", w))
        end
        if k < length(A); print(io, sep); end
    end
end

# typeinfo agnostic
"""
    print_matrix(io::IO, mat, pre, sep, post, hdots, vdots, ddots, hmod, vmod)

Prints a matrix with limited output size. If `io` sets `:limit` to true,
then only the corners of the matrix are printed, separated with vertical,
horizontal, and diagonal ellipses as appropriate.
Optional arguments are string pre (printed before the matrix, e.g. an opening bracket)
which will cause a corresponding same-size indent on following rows, and
string post (printed at the end of the last row of the matrix).
Also options to use different ellipsis characters hdots, vdots, ddots.
These are repeated every hmod or vmod elements.
"""
function print_matrix(io::IO, X::AbstractVecOrMat,
                      pre::AbstractString = " ",  # pre-matrix string
                      sep::AbstractString = "  ", # separator between elements
                      post::AbstractString = "",  # post-matrix string
                      hdots::AbstractString = "  \u2026  ",
                      vdots::AbstractString = "\u22ee",
                      ddots::AbstractString = "  \u22f1  ",
                      hmod::Integer = 5, vmod::Integer = 5)
    _print_matrix(io, inferencebarrier(X), pre, sep, post, hdots, vdots, ddots, hmod, vmod, unitrange(axes(X,1)), unitrange(axes(X,2)))
end

function _print_matrix(io, @nospecialize(X::AbstractVecOrMat), pre, sep, post, hdots, vdots, ddots, hmod, vmod, rowsA, colsA)
    hmod, vmod = Int(hmod)::Int, Int(vmod)::Int
    ncols, idxlast = length(colsA), last(colsA)
    if !(get(io, :limit, false)::Bool)
        screenheight = screenwidth = typemax(Int)
    else
        sz = displaysize(io)::Tuple{Int,Int}
        screenheight, screenwidth = sz[1] - 4, sz[2]
    end
    screenwidth -= length(pre)::Int + length(post)::Int
    presp = repeat(" ", length(pre)::Int)  # indent each row to match pre string
    postsp = ""
    @assert textwidth(hdots) == textwidth(ddots)
    sepsize = length(sep)::Int
    m, n = length(rowsA), length(colsA)
    # To figure out alignments, only need to look at as many rows as could
    # fit down screen. If screen has at least as many rows as A, look at A.
    # If not, then we only need to look at the first and last chunks of A,
    # each half a screen height in size.
    halfheight = div(screenheight,2)
    if m > screenheight
        rowsA = [rowsA[(0:halfheight-1) .+ firstindex(rowsA)]; rowsA[(end-div(screenheight-1,2)+1):end]]
    else
        rowsA = [rowsA;]
    end
    # Similarly for columns, only necessary to get alignments for as many
    # columns as could conceivably fit across the screen
    maxpossiblecols = div(screenwidth, 1+sepsize)
    if n > maxpossiblecols
        colsA = [colsA[(0:maxpossiblecols-1) .+ firstindex(colsA)]; colsA[(end-maxpossiblecols+1):end]]
    else
	    colsA = [colsA;]
    end
    A = alignment(io, X, rowsA, colsA, screenwidth, screenwidth, sepsize, ncols)
    # Nine-slicing is accomplished using print_matrix_row repeatedly
    if m <= screenheight # rows fit vertically on screen
        if n <= length(A) # rows and cols fit so just print whole matrix in one piece
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X,A,i,colsA,sep,idxlast)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        else # rows fit down screen but cols don't, so need horizontal ellipsis
            c = div(screenwidth-length(hdots)::Int+1,2)+1  # what goes to right of ellipsis
            Ralign = reverse(alignment(io, X, rowsA, reverse(colsA), c, c, sepsize, ncols)) # alignments for right
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize, ncols) # alignments for left of ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep,idxlast)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                print_matrix_row(io, X, Ralign, i, (n - length(Ralign)) .+ colsA, sep, idxlast)
                print(io, i == last(rowsA) ? post : postsp)
                if i != last(rowsA); println(io); end
            end
        end
    else # rows don't fit so will need vertical ellipsis
        if n <= length(A) # rows don't fit, cols do, so only vertical ellipsis
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X,A,i,colsA,sep,idxlast)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    print_matrix_vdots(io, vdots, A, sep, vmod, 1, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        else # neither rows nor cols fit, so use all 3 kinds of dots
            c = div(screenwidth-length(hdots)::Int+1,2)+1
            Ralign = reverse(alignment(io, X, rowsA, reverse(colsA), c, c, sepsize, ncols))
            c = screenwidth - sum(map(sum,Ralign)) - (length(Ralign)-1)*sepsize - length(hdots)::Int
            Lalign = alignment(io, X, rowsA, colsA, c, c, sepsize, ncols)
            r = mod((length(Ralign)-n+1),vmod) # where to put dots on right half
            for i in rowsA
                print(io, i == first(rowsA) ? pre : presp)
                print_matrix_row(io, X,Lalign,i,colsA[1:length(Lalign)],sep,idxlast)
                print(io, (i - first(rowsA)) % hmod == 0 ? hdots : repeat(" ", length(hdots)::Int))
                print_matrix_row(io, X,Ralign,i,(n-length(Ralign)).+colsA,sep,idxlast)
                print(io, i == last(rowsA) ? post : postsp)
                if i != rowsA[end] || i == rowsA[halfheight]; println(io); end
                if i == rowsA[halfheight]
                    print(io, i == first(rowsA) ? pre : presp)
                    print_matrix_vdots(io, vdots, Lalign, sep, vmod, 1, true)
                    print(io, ddots)
                    print_matrix_vdots(io, vdots, Ralign, sep, vmod, r, false)
                    print(io, i == last(rowsA) ? post : postsp * '\n')
                end
            end
        end
        if isempty(rowsA)
            print(io, pre)
            print(io, vdots)
            length(colsA) > 1 && print(io, "    ", ddots)
            print(io, post)
        end
    end
end

# typeinfo agnostic
# n-dimensional arrays
show_nd(io::IO, a::AbstractArray, print_matrix::Function, show_full::Bool) =
    _show_nd(io, inferencebarrier(a), print_matrix, show_full, map(unitrange, axes(a)))

function _show_nd(io::IO, @nospecialize(a::AbstractArray), print_matrix::Function, show_full::Bool, axs::Tuple{Vararg{AbstractUnitRange}})
    limit::Bool = get(io, :limit, false)
    if isempty(a)
        return
    end
    tailinds = tail(tail(axs))
    nd = ndims(a)-2
    show_full || print(io, "[")
    Is = CartesianIndices(tailinds)
    lastidxs = first(Is).I
    reached_last_d = false
    for I in Is
        idxs = I.I
        if limit
            for i = 1:nd
                ii = idxs[i]
                ind = tailinds[i]
                if length(ind) > 10
                    if ii == ind[firstindex(ind)+3] && all(d->idxs[d]==first(tailinds[d]),1:i-1)
                        for j=i+1:nd
                            szj = length(axs[j+2])
                            indj = tailinds[j]
                            if szj>10 && first(indj)+2 < idxs[j] <= last(indj)-3
                                @goto skip
                            end
                        end
                        print(io, ";"^(i+2))
                        print(io, " \u2026 ")
                        show_full && print(io, "\n\n")
                        @goto skip
                    end
                    if ind[firstindex(ind)+2] < ii <= ind[end-3]
                        @goto skip
                    end
                end
            end
        end
        if show_full
            _show_nd_label(io, a, idxs)
        end
        slice = view(a, axs[1], axs[2], idxs...)
        if show_full
            print_matrix(io, slice)
            print(io, idxs == map(last,tailinds) ? "" : "\n\n")
        else
            idxdiff = lastidxs .- idxs .< 0
            if any(idxdiff)
                lastchangeindex = 2 + findlast(idxdiff)
                print(io, ";"^lastchangeindex)
                lastchangeindex == ndims(a) && (reached_last_d = true)
                print(io, " ")
            end
            print_matrix(io, slice)
        end
        @label skip
        lastidxs = idxs
    end
    if !show_full
        reached_last_d || print(io, ";"^(nd+2))
        print(io, "]")
    end
end

function _show_nd_label(io::IO, a::AbstractArray, idxs)
    print(io, "[:, :, ")
    for i = 1:length(idxs)-1
        print(io, idxs[i], ", ")
    end
    println(io, idxs[end], "] =")
end

# print_array: main helper functions for show(io, text/plain, array)
# typeinfo agnostic
# Note that this is for showing the content inside the array, and for `MIME"text/plain".
# There are `show(::IO, ::A) where A<:AbstractArray` methods that don't use this
# e.g. show_vector, show_zero_dim
print_array(io::IO, X::AbstractArray{<:Any, 0}) =
    isassigned(X) ? show(io, X[]) : print(io, undef_ref_str)
print_array(io::IO, X::AbstractVecOrMat) = print_matrix(io, X)
print_array(io::IO, X::AbstractArray) = show_nd(io, X, print_matrix, true)

# typeinfo aware
# implements: show(io::IO, ::MIME"text/plain", X::AbstractArray)
function show(io::IO, ::MIME"text/plain", X::AbstractArray)
    if isempty(X) && (get(io, :compact, false) || X isa Vector)
        return show(io, X)
    end
    # 0) show summary before setting :compact
    summary(io, X)
    isempty(X) && return
    print(io, ":")
    show_circular(io, X) && return

    # 1) compute new IOContext
    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        # override usual show method for Vector{Method}: don't abbreviate long lists
        io = IOContext(io, :limit => false)
    end

    if get(io, :limit, false) && displaysize(io)[1]-4 <= 0
        return print(io, " …")
    else
        println(io)
    end

    # 2) update typeinfo
    #
    # it must come after printing the summary, which can exploit :typeinfo itself
    # (e.g. views)
    # we assume this function is always called from top-level, i.e. that it's not nested
    # within another "show" method; hence we always print the summary, without
    # checking for current :typeinfo (this could be changed in the future)
    io = IOContext(io, :typeinfo => eltype(X))

    # 2) show actual content
    recur_io = IOContext(io, :SHOWN_SET => X)
    print_array(recur_io, X)
end

## printing with `show`

### non-Vector arrays

# _show_nonempty & _show_empty: main helper functions for show(io, X)
# typeinfo agnostic

"""
`_show_nonempty(io, X::AbstractMatrix, prefix)` prints matrix X with opening and closing square brackets,
preceded by `prefix`, supposed to encode the type of the elements.
"""
_show_nonempty(io::IO, X::AbstractMatrix, prefix::String) =
    _show_nonempty(io, inferencebarrier(X), prefix, false, axes(X))

function _show_nonempty(io::IO, @nospecialize(X::AbstractMatrix), prefix::String, drop_brackets::Bool, axs::Tuple{AbstractUnitRange,AbstractUnitRange})
    @assert !isempty(X)
    limit = get(io, :limit, false)::Bool
    indr, indc = axs
    nr, nc = length(indr), length(indc)
    rdots, cdots = false, false
    rr1, rr2 = unitrange(indr), 1:0
    cr1 = unitrange(indc)
    cr2 = first(cr1) .+ (0:-1)
    if limit
        if nr > 4
            rr1, rr2 = rr1[1:2], rr1[nr-1:nr]
            rdots = true
        end
        if nc > 4
            cr1, cr2 = cr1[1:2], cr1[nc-1:nc]
            cdots = true
        end
    end
    drop_brackets || print(io, prefix, "[")
    for rr in (rr1, rr2)
        for i in rr
            for cr in (cr1, cr2)
                for j in cr
                    j > first(cr) && print(io, " ")
                    if !isassigned(X,i,j)
                        print(io, undef_ref_str)
                    else
                        el = X[i,j]
                        show(io, el)
                    end
                end
                if last(cr) == last(indc)
                    i < last(indr) && print(io, "; ")
                elseif cdots
                    print(io, " \u2026 ")
                end
            end
        end
        last(rr) != last(indr) && rdots && print(io, "\u2026 ; ")
    end
    if !drop_brackets
        nc > 1 || print(io, ";;")
        print(io, "]")
    end
    return nothing
end

_show_nonempty(io::IO, X::AbstractArray, prefix::String) =
    show_nd(io, X, (io, slice) -> _show_nonempty(io, inferencebarrier(slice), prefix, true, axes(slice)), false)

# a specific call path is used to show vectors (show_vector)
_show_nonempty(::IO, ::AbstractVector, ::String) =
    error("_show_nonempty(::IO, ::AbstractVector, ::String) is not implemented")

_show_nonempty(io::IO, X::AbstractArray{T,0} where T, prefix::String) = print_array(io, X)

# NOTE: it's not clear how this method could use the :typeinfo attribute
function _show_empty(io::IO, X::Array)
    show(io, typeof(X))
    print(io, "(undef, ", join(size(X),", "), ')')
end
_show_empty(io, X::AbstractArray) = summary(io, X)

# typeinfo aware (necessarily)
function show(io::IO, X::AbstractArray)
    ndims(X) == 0 && return show_zero_dim(io, X)
    ndims(X) == 1 && return show_vector(io, X)
    prefix, implicit = typeinfo_prefix(io, X)
    if !implicit
        io = IOContext(io, :typeinfo => eltype(X))
    end
    isempty(X) ?
        _show_empty(io, X) :
        _show_nonempty(io, X, prefix)
end

### 0-dimensional arrays (#31481)
show_zero_dim(io::IO, X::BitArray{0}) = print(io, "BitArray(", Int(X[]), ")")
function show_zero_dim(io::IO, X::AbstractArray{T, 0}) where T
    if isassigned(X)
        print(io, "fill(")
        show(io, X[])
    else
        print(io, "Array{", T, ", 0}(")
        show(io, undef)
    end
    print(io, ")")
end

### Vector arrays

# typeinfo aware
# NOTE: v is not constrained to be a vector, as this function can work with iterables
# in general (it's used e.g. by show(::IO, ::Set))
function show_vector(io::IO, v, opn='[', cls=']')
    prefix, implicit = typeinfo_prefix(io, v)
    print(io, prefix)
    # directly or indirectly, the context now knows about eltype(v)
    if !implicit
        io = IOContext(io, :typeinfo => eltype(v))
    end
    limited = get(io, :limit, false)

    if limited && length(v) > 20
        axs1 = axes1(v)
        f, l = first(axs1), last(axs1)
        show_delim_array(io, v, opn, ",", "", false, f, f+9)
        print(io, "  …  ")
        show_delim_array(io, v, "", ",", cls, false, l-9, l)
    else
        show_delim_array(io, v, opn, ",", cls, false)
    end
end

## Logic for displaying type information

# given type `typeinfo` extracted from context, assuming a collection
# is being displayed, deduce the elements type; in spirit this is
# similar to `eltype` (except that we don't want a default fall-back
# returning Any, as this would cause incorrect printing in e.g. `Vector[Any[1]]`,
# because eltype(Vector) == Any so `Any` wouldn't be printed in `Any[1]`)
typeinfo_eltype(typeinfo) = nothing # element type not precisely known
typeinfo_eltype(typeinfo::Type{<:AbstractArray{T}}) where {T} = eltype(typeinfo)
typeinfo_eltype(typeinfo::Type{<:AbstractDict{K,V}}) where {K,V} = eltype(typeinfo)
typeinfo_eltype(typeinfo::Type{<:AbstractSet{T}}) where {T} = eltype(typeinfo)

# types that can be parsed back accurately from their un-decorated representations
function typeinfo_implicit(@nospecialize(T))
    if T === Float64 || T === Int || T === Char || T === String || T === Symbol ||
        issingletontype(T)
        return true
    end
    return isconcretetype(T) &&
        ((T <: Array && typeinfo_implicit(eltype(T))) ||
         ((T <: Tuple || T <: Pair) && all(typeinfo_implicit, fieldtypes(T))) ||
         (T <: AbstractDict && typeinfo_implicit(keytype(T)) && typeinfo_implicit(valtype(T))))
end

# X not constrained, can be any iterable (cf. show_vector)
function typeinfo_prefix(io::IO, X)
    typeinfo = get(io, :typeinfo, Any)::Type

    if !(X isa typeinfo)
        typeinfo = Any
    end

    # what the context already knows about the eltype of X:
    eltype_ctx = typeinfo_eltype(typeinfo)
    eltype_X = eltype(X)

    if X isa AbstractDict
        if eltype_X == eltype_ctx
            sprint(show_type_name, typeof(X).name), false
        elseif !isempty(X) && typeinfo_implicit(keytype(X)) && typeinfo_implicit(valtype(X))
            sprint(show_type_name, typeof(X).name), true
        else
            string(typeof(X)), false
        end
    else
        # Types hard-coded here are those which are created by default for a given syntax
        if eltype_X == eltype_ctx
            "", false
        elseif !isempty(X) && typeinfo_implicit(eltype_X)
            "", true
        elseif print_without_params(eltype_X)
            sprint(show_type_name, unwrap_unionall(eltype_X).name), false # Print "Array" rather than "Array{T,N}"
        else
            string(eltype_X), false
        end
    end
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
import CSV

# Is this the idiomatic way of specifying types?
# https://csv.juliadata.org/stable/#Providing-Types
file = CSV.File("input_file.txt", header=false; types=[Int])

# Specifying types makes this crash somehow.  Ask Lorenzo later.
complement_dict = Dict()
original_values = Set()

total = 2020

for row1 in file
    for row2 in file
        val1 = row1[1]
        val2 = row2[1]

        two_sum = val1 + val2
        if two_sum < 2020
        # all values are positive so
            complement = total - two_sum
            complement_dict[Set([val1, val2])] = complement
            # println("Two sum $val1 $val2 and $complement")
        end

    push!(original_values, val1)
    end
end

key1 = 0
key2 = 0
key3 = 0

for (k, v) in complement_dict
    if v in original_values
        tmp_set = collect(k)
        global key1 = tmp_set[1]
        global key2 = tmp_set[2]
        global key3 = v
    end
end

if key1 == 0 && key2 == 0 && key3 == 0
    println("Error!!!")
else
    product = key1 * key2 * key3
    println("The two values are $key1 and $key2 and $key3.  The product is $product")
end

ext: jl


text:
using Pkg
Pkg.add.([
          # tabular
          "DataFrames", "CSV", "Query",

          # JuliaPro
          "DataStructures", "LightGraphs", "Calculus", "DataFrames",
          "StatsBase", "Distributions", "HypothesisTests", "GLM",
          "JuMP", "Optim", "Roots", "ODBC", "JDBC",
          "Knet", "Clustering", "DecisionTree", "PyCall",
          "JSON", "HDF5", "JLD", "QuantEcon",
          "SQLite",

          # JuliaPlots
          "Plots", "PyPlot", "PlotlyJS",

          # DSP
          "DSP", "Wavelets", "Deconvolution",

          # "MXNet"  # not supported yet

          # MISC
          "ProgressMeter",
         ])

max_stars_repo_name: deeptimhe/LeinaoPAI

max_issues_repo_name: deeptimhe/LeinaoPAI

max_forks_repo_name: deeptimhe/LeinaoPAI


text:
using DynamicPPL
using Distributions
using Bijectors

@testset "distribution_wrappers.jl" begin
    d = Normal()
    nd = DynamicPPL.NoDist(d)

    # Smoke test
    rand(nd)

    # Actual tests
    @test minimum(nd) == -Inf
    @test maximum(nd) == Inf
    @test logpdf(nd, 15.0) == 0
    @test Bijectors.logpdf_with_trans(nd, 0) == 0
end


text:
context("Hyper Cubes") do
    context("Core") do
        x = centered(HyperCube)
        @fact origin(x) --> Vec3f0(-0.5)
        @fact widths(x) --> Vec3f0(1.0)
        @fact maximum(x) --> Vec3f0(0.5)
        @fact minimum(x) --> Vec3f0(-0.5)
    end
end
context("Hyper Rectangles") do

context("constructors and containment") do
    # HyperRectangle(vals...)
    a = HyperRectangle(0,0,1,1)
    @fact a --> HyperRectangle{2,Int}(Vec(0,0),Vec(1,1))

    b = HyperRectangle(0,0,1,1,1,2)
    @fact b --> HyperRectangle{3,Int}(Vec(0,0,1),Vec(1,1,2))

    a = AABB(0,0,0,1,1,1)
    b = AABB(Vec(0,0,0),Vec(1,1,1))
    @fact a --> HyperRectangle(0,0,0,1,1,1)
    @fact a --> b

    a = HyperRectangle{4, Float64}()
    @fact a --> HyperRectangle{4, Float64}(Vec(Inf,Inf,Inf,Inf), Vec(-Inf,-Inf,-Inf,-Inf))

    a = update(a, Vec(1,2,3,4))

    @fact a --> HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(0.0))
    @fact a == HyperRectangle{4, Float64}(Vec(3.0,2.0,3.0,4.0),Vec(1.0,2.0,3.0,4.0)) --> false

    a = update(a, Vec(5,6,7,8))

    b = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(4.0))
    @fact widths(b) --> Vec(4., 4., 4., 4.)
    @fact a --> b
    @fact isequal(a,b) --> true

    @fact maximum(a) --> Vec(5.0,6.0,7.0,8.0)
    @fact minimum(a) --> Vec(1.0,2.0,3.0,4.0)
    @fact origin(a) --> Vec(1.0,2.0,3.0,4.0)

    @fact_throws MethodError HyperRectangle(Vec(1.0,2.0,3.0), Vec(1.0,2.0,3.0,4.0))

    @fact (in(a,b) && in(b,a) && contains(a,b) && contains(b,a)) --> true

    c = HyperRectangle(Vec(1.1,2.1,3.1,4.1),Vec4f0(3.8))

    @fact (!in(a,c) && in(c,a) && contains(a,c) && !contains(c,a)) --> true

    #conversion
    h = HyperRectangle(0,0,1,1)
    @fact HyperRectangle{2,Float64}(h) --> HyperRectangle(0.,0.,1.,1.)

    # AABB
    a = AABB(0,0,1,1)
    @fact a --> HyperRectangle{2,Int}(Vec(0,0),Vec(1,1))

    centered_rect = centered(HyperRectangle)
    @fact centered_rect --> HyperRectangle{3,Float32}(Vec3f0(-0.5),Vec3f0(1))
    centered_rect = centered(HyperRectangle{2})
    @fact centered_rect --> HyperRectangle{2,Float32}(Vec2f0(-0.5),Vec2f0(1))

    centered_rect = centered(HyperRectangle{2, Float64})
    @fact centered_rect --> HyperRectangle{2,Float64}(Vec(-0.5, -0.5),Vec(1.,1.))

    centered_rect = centered(HyperRectangle{3, Float32})
    @fact centered_rect --> HyperRectangle{3,Float64}(Vec(-0.5,-0.5,-0.5),Vec(1.,1.,1.))

end

# Testing split function
context("Testing split function") do
    d = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(1.0))
    d1, d2 = split(d, 3, 3.5)

    @fact maximum(d1)[3] --> 3.5
    @fact minimum(d1)[3] --> 3.0
    @fact maximum(d2)[3] --> 4.0
    @fact minimum(d2)[3] --> 3.5

    d = HyperRectangle(0,0,2,2)
    d1, d2 = split(d, 1, 1)
    @fact d1 --> HyperRectangle(0,0,1,2)
    @fact d2 --> HyperRectangle(1,0,1,2)
end

context("Test distance functions") do
    a = HyperRectangle(Vec(0.0,0.0), Vec(1.0, 1.0))
    b = HyperRectangle(Vec(2.0,3.0), Vec(1.0, 1.0))
    p = Vec(2.5, 1.5)

    # Rect - Rect
    @fact min_dist_dim(a, b, 1) --> 1.0
    @fact min_dist_dim(a, b, 2) --> 2.0
    @fact max_dist_dim(a, b, 1) --> 3.0
    @fact max_dist_dim(a, b, 2) --> 4.0

    @fact min_euclideansq(a, b) --> 5.0
    @fact max_euclideansq(a, b) --> 25.0
    @fact minmax_euclideansq(a, b) --> (5.0, 25.0)

    @fact min_euclidean(a, b) --> sqrt(5.0)
    @fact max_euclidean(a, b) --> sqrt(25.0)
    @fact minmax_euclidean(a, b) --> (sqrt(5.0), sqrt(25.0))

    # Rect - Point
    @fact min_dist_dim(a, p, 1) --> 1.5
    @fact max_dist_dim(a, p, 1) --> 2.5
    @fact minmax_dist_dim(a, p, 1) --> (1.5, 2.5)

    @fact min_dist_dim(a, p, 2) --> 0.5
    @fact max_dist_dim(a, p, 2) --> 1.5
    @fact minmax_dist_dim(a, p, 2) --> (0.5, 1.5)

    @fact min_euclideansq(a, p) --> 2.5
    @fact max_euclideansq(a, p) --> 8.5
    @fact minmax_euclideansq(a, p) --> (2.5, 8.5)

    @fact min_euclidean(a, p) --> sqrt(2.5)
    @fact max_euclidean(a, p) --> sqrt(8.5)
    @fact minmax_euclidean(a, p) --> (sqrt(2.5), sqrt(8.5))

    p2 = Vec(0.75, 0.75)
    @fact min_dist_dim(a, p2, 1) --> 0.0
    @fact min_dist_dim(a, p2, 2) --> 0.0

    b2  = HyperRectangle(Vec(0.25, 0.25), Vec(0.75,0.75))
    @fact min_dist_dim(a, b2, 1) --> 0
    @fact min_dist_dim(a, b2, 2) --> 0
end

context("set operations") do
    a = HyperRectangle(Vec(0,0),Vec(1,1))
    b = HyperRectangle(Vec(1,1),Vec(1,1))

    @fact union(a,b) --> union(b,a)
    @fact union(a,b) --> HyperRectangle(Vec(0,0), Vec(2,2))
    @fact intersect(a,b) --> intersect(b,a)
    @fact intersect(a,b) --> HyperRectangle(Vec(1,1), Vec(0,0))
    @fact diff(a,b) --> a
    @fact diff(b,a) --> b

    c = HyperRectangle(Vec(0,0),Vec(2,2))
    d = HyperRectangle(Vec(1,1),Vec(2,2))
    @fact union(c,d) --> union(d,c)
    @fact union(c,d) --> HyperRectangle(Vec(0,0), Vec(3,3))
    @fact intersect(c,d) --> intersect(d,c)
    @fact intersect(c,d) --> HyperRectangle(Vec(1,1), Vec(1,1))
    @fact diff(c,d) --> c
    @fact diff(d,c) --> d

end

# fact relations
context("relations") do
    a = HyperRectangle(Vec(0,0),Vec(1,1))
    b = HyperRectangle(Vec(1,1),Vec(1,1))
    c = HyperRectangle(Vec(0,0),Vec(2,2))
    d = HyperRectangle(Vec(1,1),Vec(2,2))
    e = HyperRectangle(Vec(3,3),Vec(1,1))
    d = HyperRectangle(Vec(0.25,0.25),Vec(0.5,0.5))
    f = HyperRectangle(Vec(0.9,0.9), Vec(1.1,1.1))

    @fact (finishes(b,c) && !finishes(c,b)) --> true
    @fact (!finishes(a,b)) --> true
    @fact (meets(a,b) && !meets(b,a)) --> true
    @fact (before(a,e) && !before(e,a)) --> true
    @fact (during(d,a) && !during(a,d)) --> true
    @fact (starts(a,c) && !starts(c,a)) --> true
    @fact (!starts(a,b)) --> true
    @fact (!overlaps(a,b) && !overlaps(b,a)) --> true
    @fact (overlaps(a,f) && !overlaps(f,a)) --> true
end

context("point membership") do
    a = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(4.0))
    @fact (in(Vec(4,5,6,7),a) && contains(a,Vec(4,5,6,7))) --> true
    @fact (in(Vec(5,6,7,8),a) && contains(a,Vec(5,6,7,8))) --> true
    @fact (!in(Vec(6,7,8,9),a) && !contains(a,Vec(6,7,8,9))) --> true
end

context("from Points") do
    a = HyperRectangle([Point(1,1), Point(2,3), Point(4,5), Point(0,-1)])
    @fact a --> HyperRectangle(0,-1,4,6)
    a = HyperRectangle{3,Int}([Point(1,1), Point(2,3), Point(4,5), Point(0,-1)])
    @fact a --> HyperRectangle(0,-1,0,4,6,0)
end

context("transforms") do
    t = Mat((1,0,0),(0,1,0),(1,2,0))
    h = t*HyperRectangle(0,0,1,1)
    @fact h --> HyperRectangle(1,2,1,1)
    t = Mat((0,1),(1,0))
    h = t*HyperRectangle(0,0,1,2)
    @fact h --> HyperRectangle(0,0,2,1)
end

end


text:
using Base.Test
using SpikeSorter

include("io.jl")


text:
@doc raw"""
    LinearAffineMetric <: Metric

The linear affine metric is the metric for symmetric positive definite matrices, that employs
matrix logarithms and exponentials, which yields a linear and affine metric.
"""
struct LinearAffineMetric <: RiemannianMetric end

default_metric_dispatch(::SymmetricPositiveDefinite, ::LinearAffineMetric) = Val(true)

@doc raw"""
    distance(M::SymmetricPositiveDefinite, p, q)
    distance(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, q)

Compute the distance on the [`SymmetricPositiveDefinite`](@ref) manifold between `p` and `q`,
as a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

```math
d_{\mathcal P(n)}(p,q)
= \lVert \operatorname{Log}(p^{-\frac{1}{2}}qp^{-\frac{1}{2}})\rVert_{\mathrm{F}}.,
```
where $\operatorname{Log}$ denotes the matrix logarithm and
$\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
"""
function distance(M::SymmetricPositiveDefinite{N}, p, q) where {N}
    cq = cholesky(q)
    s = eigvals(Symmetric(cq.L \ p / cq.U))
    return any(s .<= eps()) ? zero(eltype(p)) : sqrt(sum(abs.(log.(s)) .^ 2))
end

@doc raw"""
    exp(M::SymmetricPositiveDefinite, p, X)
    exp(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, p, X)

Compute the exponential map from `p` with tangent vector `X` on the
[`SymmetricPositiveDefinite`](@ref) `M` with its default [`MetricManifold`](@ref) having the
[`LinearAffineMetric`](@ref). The formula reads

```math
\exp_p X = p^{\frac{1}{2}}\operatorname{Exp}(p^{-\frac{1}{2}} X p^{-\frac{1}{2}})p^{\frac{1}{2}},
```

where $\operatorname{Exp}$ denotes to the matrix exponential.
"""
exp(::SymmetricPositiveDefinite, ::Any...)

function exp!(M::SymmetricPositiveDefinite{N}, q, p, X) where {N}
    e = eigen(Symmetric(p))
    U = e.vectors
    S = e.values
    Ssqrt = Diagonal(sqrt.(S))
    SsqrtInv = Diagonal(1 ./ sqrt.(S))
    pSqrt = Symmetric(U * Ssqrt * transpose(U))
    pSqrtInv = Symmetric(U * SsqrtInv * transpose(U))
    T = Symmetric(pSqrtInv * X * pSqrtInv)
    eig1 = eigen(T) # numerical stabilization
    Se = Diagonal(exp.(eig1.values))
    Ue = eig1.vectors
    pUe = pSqrt * Ue
    return copyto!(q, pUe * Se * transpose(pUe))
end

@doc raw"""
    [Ξ,κ] = get_basis(M::SymmetricPositiveDefinite, p, B::DiagonalizingOrthonormalBasis)
    [Ξ,κ] = get_basis(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, p, B::DiagonalizingOrthonormalBasis)

Return a orthonormal basis `Ξ` as a vector of tangent vectors (of length
[`manifold_dimension`](@ref) of `M`) in the tangent space of `p` on the
[`MetricManifold`](@ref) of [`SymmetricPositiveDefinite`](@ref) manifold `M` with
[`LinearAffineMetric`](@ref) that diagonalizes the curvature tensor $R(u,v)w$
with eigenvalues `κ` and where the direction `B.frame_direction` has curvature `0`.
"""
function get_basis(
    M::SymmetricPositiveDefinite{N},
    p,
    B::DiagonalizingOrthonormalBasis,
) where {N}
    pSqrt = sqrt(p)
    eigv = eigen(B.frame_direction)
    V = eigv.vectors
    Ξ = [
        (i == j ? 1 / 2 : 1 / sqrt(2)) *
        (V[:, i] * transpose(V[:, j]) + V[:, j] * transpose(V[:, i])) for i in 1:N
        for j in i:N
    ]
    λ = eigv.values
    κ = [-1 / 4 * (λ[i] - λ[j])^2 for i in 1:N for j in i:N]
    return CachedBasis(B, κ, Ξ)
end

function get_coordinates!(
    M::SymmetricPositiveDefinite{N},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis,
) where {N}
    dim = manifold_dimension(M)
    @assert size(Y) == (dim,)
    @assert size(X) == (N, N)
    @assert dim == div(N * (N + 1), 2)
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds Y[k] = X[i, j] * scale
        k += 1
    end
    return Y
end

function get_vector!(
    M::SymmetricPositiveDefinite{N},
    Y,
    p,
    X,
    B::DefaultOrthonormalBasis,
) where {N}
    dim = manifold_dimension(M)
    @assert size(X) == (div(N * (N + 1), 2),)
    @assert size(Y) == (N, N)
    k = 1
    for i in 1:N, j in i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds Y[i, j] = X[k] * scale
        @inbounds Y[j, i] = X[k] * scale
        k += 1
    end
    return Y
end

@doc raw"""
    inner(M::SymmetricPositiveDefinite, p, X, Y)
    inner(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, X, Y)

Compute the inner product of `X`, `Y` in the tangent space of `p` on
the [`SymmetricPositiveDefinite`](@ref) manifold `M`, as
a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

````math
g_p(X,Y) = \operatorname{tr}(p^{-1} X p^{-1} Y),
````
"""
function inner(M::SymmetricPositiveDefinite, p, X, Y)
    F = cholesky(Symmetric(p))
    return tr((F \ Symmetric(X)) * (F \ Symmetric(Y)))
end

@doc raw"""
    log(M::SymmetricPositiveDefinite, p, q)
    log(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, q)

Compute the logarithmic map from `p` to `q` on the [`SymmetricPositiveDefinite`](@ref)
as a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

```math
\log_p q =
p^{\frac{1}{2}}\operatorname{Log}(p^{-\frac{1}{2}}qp^{-\frac{1}{2}})p^{\frac{1}{2}},
```
where $\operatorname{Log}$ denotes to the matrix logarithm.
"""
log(::SymmetricPositiveDefinite, ::Any...)

function log!(M::SymmetricPositiveDefinite{N}, X, p, q) where {N}
    e = eigen(Symmetric(p))
    U = e.vectors
    S = e.values
    Ssqrt = Diagonal(sqrt.(S))
    SsqrtInv = Diagonal(1 ./ sqrt.(S))
    pSqrt = Symmetric(U * Ssqrt * transpose(U))
    pSqrtInv = Symmetric(U * SsqrtInv * transpose(U))
    T = Symmetric(pSqrtInv * q * pSqrtInv)
    e2 = eigen(T)
    Se = Diagonal(log.(max.(e2.values, eps())))
    pUe = pSqrt * e2.vectors
    return mul!(X, pUe, Se * transpose(pUe))
end

@doc raw"""
    vector_transport_to(M::SymmetricPositiveDefinite, p, X, q, ::ParallelTransport)
    vector_transport_to(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, X, y, ::ParallelTransport)

Compute the parallel transport of `X` from the tangent space at `p` to the
tangent space at `q` on the [`SymmetricPositiveDefinite`](@ref) as a
[`MetricManifold`](@ref) with the [`LinearAffineMetric`](@ref).
The formula reads

```math
\mathcal P_{q←p}X = p^{\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}p^{-\frac{1}{2}}\log_p(q)p^{-\frac{1}{2}}
\bigr)
p^{-\frac{1}{2}}X p^{-\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}p^{-\frac{1}{2}}\log_p(q)p^{-\frac{1}{2}}
\bigr)
p^{\frac{1}{2}},
```

where $\operatorname{Exp}$ denotes the matrix exponential
and `log` the logarithmic map on [`SymmetricPositiveDefinite`](@ref)
(again with respect to the [`LinearAffineMetric`](@ref)).
"""
vector_transport_to(::SymmetricPositiveDefinite, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(
    M::SymmetricPositiveDefinite{N},
    Y,
    p,
    X,
    q,
    ::ParallelTransport,
) where {N}
    distance(M, p, q) < 2 * eps(eltype(p)) && copyto!(Y, X)
    e = eigen(Symmetric(p))
    U = e.vectors
    S = e.values
    Ssqrt = sqrt.(e.values)
    SsqrtInv = Diagonal(1 ./ Ssqrt)
    Ssqrt = Diagonal(Ssqrt)
    pSqrt = Symmetric(U * Ssqrt * transpose(U)) # p^1/2
    pSqrtInv = Symmetric(U * SsqrtInv * transpose(U)) # p^(-1/2)
    tv = Symmetric(pSqrtInv * X * pSqrtInv) # p^(-1/2)Xp^{-1/2}
    ty = Symmetric(pSqrtInv * q * pSqrtInv) # p^(-1/2)qp^(-1/2)
    e2 = eigen(ty)
    Se = Diagonal(log.(e2.values))
    Ue = e2.vectors
    logty = Symmetric(Ue * Se * transpose(Ue)) # nearly log_pq without the outer p^1/2
    e3 = eigen(logty) # since they cancel with the pInvSqrt in the next line
    Sf = Diagonal(exp.(e3.values / 2)) # Uf * Sf * Uf' is the Exp
    Uf = e3.vectors
    pUe = pSqrt * Uf * Sf * transpose(Uf) # factors left of tv (and transposed right)
    vtp = Symmetric(pUe * tv * transpose(pUe)) # so this is the documented formula
    return copyto!(Y, vtp)
end


text:
using CoordinateTransformations
using StaticArrays
using Base: @propagate_inbounds

inside_annulus(rmin, rmax, center, idx...) = inside_annulus(rmin, rmax, center, SVector(idx[end-1:end]))
inside_annulus(rmin, rmax, center, idx::CartesianIndex{2}) = inside_annulus(rmin, rmax, center, SVector(idx.I))
function inside_annulus(rmin, rmax, center::AbstractVector, position::AbstractVector)
    Δ = center - position
    r = sqrt(sum(abs2, Δ))
    return rmin ≤ r ≤ rmax
end

include("annulus.jl")
include("multiannulus.jl")


text:
module UnregisteredWithoutProject
export f
f(x) = x

end

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
include("filtering.jl");
## Number of pokes during the session

Df = combine(groupby(streaks,:Phase)) do dd
    summarize(dd,:Trial,:Num_pokes)
end
Drug_colors!(Df)
@df Df plot(:Xaxis,:Mean,
    group = :Phase,
    color = :color,
    ribbon = :SEM,
    linewidth = 2,
    fillalpha = 0.1,
    xticks = 0:10:100)
savefig(joinpath(figs_loc,"Fig1/trials_selected.pdf"))


text:
export crammer_svm

using Flux

"""
Implement multi-class kernel-based SVM

# References

- http://jmlr.csail.mit.edu/papers/volume2/crammer01a/crammer01a.pdf
"""
function crammer_svm(embed_support, labels_support, embed_query)
    # construct dual QP problem
    # solve QP
    # compute compatibility
    return nothing
end


text:
#RSDeltaSigmaPort SNR calculations
#-------------------------------------------------------------------------------

"""`snr = calculateSNR(hwfft,f,nsig=1)`

Estimate the signal-to-noise ratio, given the in-band bins of a (Hann-windowed)
fft and the location of the input signal (f>0).

 - For nsig=1, the input tone is contained in `hwfft(f:f+2)`;
   this range is appropriate for a Hann-windowed fft.
 - Each increment in `nsig` adds a bin to either side.
 - The SNR is expressed in dB.
"""
function calculateSNR(hwfft, f::Int, nsig::Int = 1)
	signalBins = f-nsig+1:f+nsig+1
	signalBins = signalBins[signalBins .> 0]
	signalBins = signalBins[signalBins .<= length(hwfft)]
	s = norm(hwfft[signalBins]) #*4/(N*sqrt(3)) for true rms value
	noiseBins = collect(1:length(hwfft))
	deleteat!(noiseBins, signalBins)
	n = norm(hwfft[noiseBins])
	snr = (n!=0) ? dbv(s/n) : Inf
	return snr
end

"""`(peak_snr,peak_amp) = peakSNR(snr,amp)`

Estimate the snr peak by threading a line of slope 1 through the (amp,snr) data.
Both amp and snr are expressed in dB.
"""
function peakSNR(snr, amp)
	snr = deepcopy(snr); amp = deepcopy(amp)
	#Delete garbage data
	i = abs.(snr) .== Inf;   deleteat!(snr, i); deleteat!(amp, i)
	i = isnan.(snr);         deleteat!(snr, i); deleteat!(amp, i)
	i = snr .< 3;            deleteat!(snr, i); deleteat!(amp, i)
	n = length(amp)

	#Sort by amplitude
	i = sortperm(amp)
	amp = amp[i]; snr = snr[i]

	i = [true]
	local m = 0.0
	while any(i) && n > 3
		#Draw a 45-degree median line through the data
		tmp = sort(snr .- amp)
		if mod(n,2) == 0
			m = mean( tmp[div(n,2) .+ [1, 0]] )
		else
			m = tmp[div(n+1,2)]
		end
		#Discard data that is more than 6dB off
		i = abs.(amp .- snr .+ m) .> 6
		deleteat!(snr, i); deleteat!(amp, i)
		n = length(amp)
	end

	peak_amp = maximum(amp)
	peak_snr = peak_amp + m
	return (peak_snr, peak_amp)
end

"""`(snr,amp,k0,k1,sigma_e2) = predictSNR(ntf,R=64,amp=...,f0=0)`

Predict the SNR curve of a binary delta-sigma modulator by using
the describing function method of Ardalan and Paulos.

The modulator is specified by a noise transfer function (ntf).
The band of interest is defined by the oversampling ratio (R)
and the center frequency (f0).
The input signal is characterized by the amp vector.
amp defaults to [-120 -110...-20 -15 -10 -9 -8 ... 0]dB, where 0 dB means
a full-scale (peak value = 1) sine wave.

The algorithm assumes that the amp vector is sorted in increasing order;
once instability is detected, the remaining SNR values are set to -Inf.

# Output:
 - snr: a vector of SNR values (in dB)
 - amp: a vector of amplitudes (in dB)
 - k0: the quantizer signal gain
 - k1: the quantizer noise gain
 - sigma_e2: the power of the quantizer noise (not in dB)

The describing function method of A&P assumes that the quantizer processes
signal and noise components separately. The quantizer is modelled as two
(not necessarily equal) linear gains, k0 and k1, and an additive white
gaussian noise source of power sigma_e2. k0, k1 and sigma_e2 are calculated
as functions of the input.
The modulator's loop filter is assumed to have nearly infinite gain at
the test frequency.

Future versions may accommodate STFs.
"""
function predictSNR(ntf, R=64, amp=vcat(-120:10:-20, -15, -10:0); f0::Float64=0.0)
	Nb = 100
	XTAB = range(-2, stop=0, length=21)
	#The following code was used to create the YTAB matrix
	#YTAB = [];
	#for xi=XTAB
	#	YTAB = [YTAB; hypergeo(0.5,1,xi) hypergeo(0.5,2,xi)];
	#end
	YTAB = [0.46575960516930   0.67366999387741;
		0.47904652357101   0.68426650762558;
		0.49316295981407   0.69527947902679;
		0.50817364454269   0.70673173666000;
		0.52414894104004   0.71864765882492;
		0.54116523265839   0.73105299472809;
		0.55930554866791   0.74397552013397;
		0.57866013050079   0.75744456052780;
		0.59932720661163   0.77149158716202;
		0.62141352891922   0.78615015745163;
		0.64503526687622   0.80145609378815;
		0.67031890153885   0.81744754314423;
		0.69740217924118   0.83416539430618;
		0.72643494606018   0.85165339708328;
		0.75758063793182   0.86995816230774;
		0.79101717472076   0.88912981748581;
		0.82693856954575   0.90922164916992;
		0.86555624008179   0.93029111623764;
		0.90710091590881   0.95239937305450;
		0.95182400941849   0.97561222314835;
		1.00000000000000   1.00000000000000;
	]

	if f0==0
		band_of_interest = range(0, stop=pi/R, length=Nb)
	else
		band_of_interest = range(2*pi*(f0-0.25/R), stop=2*pi*(f0+0.25/R), length=Nb)
	end

	num, den = _zp2tf(ntf.z, ntf.p, 1)
	num1 = num - den

	N = length(amp)
	snr = zeros(1,N) .- Inf
	k0 = zeros(1,N)
	k1 = zeros(1,N)
	sigma_e2 = zeros(1,N)

	u = 10 .^ (amp/20)

	Nimp = 100
	unstable = false
	f_prev = 0.0; fprime = 0.0
	m0 = 0.0; rho = 0.0; erfinvu = 0.0
	for n=1:N
		#Calculate sigma_e2
		if f0==0
			erfinvu = erfinv(u[n])
			sigma_e2[n] = 1 - u[n]^2 - 2/pi*exp(-2*erfinvu^2)
		else #Sinusoidal input.
			#Solve sqrt(pi)*u/2 = rho * hypergeo(0.5,2,-rho^2);
			#Formulate as solve f(rho)=0, f = rho*M(0.5,2,-rho^2)-K
			#and use the secant method.
			K = 0.5 * sqrt(pi) * u[n]
			if n==1
				rho = u[n]^2 #Initial guess; otherwise use previous value.
				fprime = 1.0
			end
			drho = 1.0
			for itn = 1:20
				m0 = interp1_cubic(XTAB,YTAB[:,2],-rho^2)
				f = rho*m0 - K
				if( itn >1 )
					fprime = max((f-f_prev)/drho,0.5) #Secant approx.
				end
				if abs(f) < 1e-8; break; end #!Converged
				drho = -f/fprime
				if abs(drho) > 0.2; drho = sign(drho)*0.2; end
				if abs(drho) < 1e-6; break; end #!Converged
				rho = rho + drho
				f_prev = f
			end
			m1 = interp1_cubic(XTAB,YTAB[:,1], -rho^2)
			sigma_e2[n] = 1 - u[n]^2/2 - 2/pi*m1^2
		end #Sinusoidal input

		#Iterate to solve for k1 and sigma_1.
		#Using one of MATLAB's nonlinear equation solvers would be more efficient,
		#but this function code would then require the optimization toolbox.
		#!Future work: put in 2-D BFGS code.
		k1[n] = (n>1) ? k1[n-1] : 1.2 #Use the previous value of k1 as the initial guess.
		k1_prev = 0
		itn = 0
		k1sigma1 = (f0==0) ? sqrt(2/pi)*exp(-erfinvu^2) : sqrt(2/pi)*m1
		local sigma_1 = 0.0
		while abs(k1[n]-k1_prev) > 1e-6*(1+k1[n]) && itn < 100
			#Create the function: H_hat = L1/(1-k1*L1)=(H-1)/(H*(1-k1)+k1).
			den1 = (1-k1[n])*num + den*k1[n]
			#Calculate pGain, the square of the 2-norm of H_hat.
			pGain, Nimp = powerGain(num1, den1, Nimp)
			if isinf(pGain)
				unstable = true
				break
			end

			sigma_1 = sqrt(pGain*sigma_e2[n])
			k1_prev = k1[n]
			k1[n] = k1sigma1/sigma_1
			itn = itn+1
		end
		if unstable
			break
		end

		if f0==0
			y0 = sqrt(2)*erfinvu*sigma_1
			k0[n] = u[n]/y0
		else
			k0[n] = sqrt(2/pi)*m0/sigma_1
		end

		h = _freqz(num, (1-k1[n])*num + k1[n]*den, band_of_interest)
		#For both DC and sine wave inputs, use u^2/2 as the signal
		#power since true DC measurements are usually impossible.
		snr[n] = dbp( 0.5*u[n]^2 / (sum(h.^2)/(R*Nb)*sigma_e2[n]) )
	end

	return (snr, amp, k0, k1, sigma_e2)
end

"""`(snr,amp) = simulateSNR(ΔΣmodel,osr,amp; f0=0,nlev=2,f=1/(4*osr),k=13,quadrature=false)`

Determine the SNR for a delta-sigma modulator by using simulations.

`ΔΣmodel` describes the modulator by a noise transfer function (ntf) and the
number of quantizer levels (`nlev`). Alternatively, `ΔΣmodel` may be a
function taking the input signal as its sole argument:
 - `ΔΣmodel::ZPKData`: z, p, k transfer function
 - `ΔΣmodel::Array` ABCD matrix
 - `ΔΣmodel::Function` funΔΣ(u) generating ΔΣ output.

The band of interest is defined by the oversampling ratio (`osr`)
and the center frequency (`f0`).
The input signal is characterized by the `amp` vector and the `f` variable.

 - `amp` defaults to [-120 -110...-20 -15 -10 -9 -8 ... 0]dB, where 0 dB means
   a full-scale (peak value = nlev-1) sine wave.
 - `f` is the input frequency, normalized such that 1 -> fs;
 - `f` is rounded to an FFT bin.

Using sine waves located in FFT bins, the SNR is calculated as the ratio
of the sine wave power to the power in all in-band bins other than those
associated with the input tone. Due to spectral smearing, the input tone
is not allowed to lie in bins 0 or 1. The length of the FFT is 2^k.

 - If ntf is complex, simulateQDSM (which is slow) is called.
 - If ABCD is complex, simulateDSM is used with the real equivalent of ABCD
   in order to speed up simulations.

Future versions may accommodate STFs.
"""
function simulateSNR(arg1, osr::Int=64, amp=vcat(-120:10:-20, -15, -10:0);
	f0::Float64=0.0, nlev::Int=2, f::Float64=NaN, k::Int=13, quadrature::Bool=false
)
	is_function = isa(arg1, Function)

	#Look at arg1 and decide if the system is quadrature
	quadrature_ntf = false
	if !is_function
		if isa(arg1, ZPKData)
			pz = [arg1.p, arg1.z]
			for i=1:2
				if any( abs.(imag.(poly(pz[i]))) .> 1e-4 )
					quadrature = true
					quadrature_ntf = true
				end
			end
		else #ABCD matrix
			if !all(all(imag.(arg1) .== 0))
				quadrature = true
			end
		end
	end

	osr_mult = 2
	if f0!=0 && !quadrature
   	osr_mult = 2*osr_mult
	end
	if isnan(f)
		f = f0 + 0.5/(osr*osr_mult) #Halfway across the band
	end
	M = nlev-1
	if quadrature && !quadrature_ntf
		#Modify arg1 (ABCD) and nlev so that simulateDSM can be used
		nlev = [nlev; nlev]
		arg1 = mapQtoR(arg1)
	end

	if abs(f-f0) > 1/(osr*osr_mult)
		@warn("The input tone is out-of-band.")
	end

	N = 2^k
	if N < 8*2*osr #Require at least 8 bins to be "in-band"
		@warn("Increasing k to accommodate a large oversampling ratio.")
		k = ceil(Int, log2(8*2*osr))
		N = 2^k
	end
	F = round(Int, f*N)
	if abs(F)<=1
		@warn("Increasing k to accommodate a low input frequency.")
		#Want f*N >= 1
		k = ceil(Int, log2(1/f))
		N = 2^k
		F = 2
	end

	Ntransient = 100
	soft_start = 0.5*(1 .- cos.(2π/Ntransient * (0:div(Ntransient,2)-1)))
	if !quadrature
		tone = M * sin.(2π*F/N * (0:(N+Ntransient-1)))
		tone[1:div(Ntransient,2)] = tone[1:div(Ntransient,2)] .* soft_start
	else
		tone = M * exp.(2π*j*F/N * (0:(N+Ntransient-1)) )
		tone[1:div(Ntransient,2)] = tone[1:div(Ntransient,2)] .* soft_start
		if !quadrature_ntf
			tone = [real(tone); imag(tone)]
		end
	end
	window = .5*(1 .- cos.(2π*(0:N-1)/N) ) #Hann window
	if f0==0
		#Exclude DC and its adjacent bin
		inBandBins = div(N,2) .+ (3:round(Int, N/(osr_mult*osr)))
		F = F-2
	else
		f1 = round(Int, N*(f0-1/(osr_mult*osr)))
		f2 = round(Int, N*(f0+1/(osr_mult*osr)))
		inBandBins = div(N,2) .+ (f1:f2) #Should exclude DC
		F = F-f1+1
	end

	snr = zeros(size(amp))
	i=1
	for A = 10 .^ (amp/20)
		if is_function
			v = arg1(A*tone)
		elseif quadrature_ntf
			v = simulateQDSM(A*tone, arg1, nlev=nlev).v
		else
			v = simulateDSM(A*tone, arg1, nlev=nlev).v
			if quadrature
				v = v[1,:] + im*v[2,:]
			end
		end
		hwfft = fftshift(fft(applywnd(v[1+Ntransient:N+Ntransient], window)))
		snr[i] = calculateSNR(hwfft[inBandBins], F)
		i=i+1
	end

	return (snr,amp)
end

#Not meant to be called externally:
function _calcSNRInfo(NTF, OSR::Int, M::Int, f0::Float64)
	nlev=M+1
	SNR, amp = simulateSNR(NTF, OSR, f0=f0, nlev=nlev)
	pk_snr, pk_amp = peakSNR(SNR, amp) #Single values

	vs_amp_predicted = nothing
	if nlev == 2 #predictSNR() only supports 2 levels
		SNR_pred, amp_pred = predictSNR(NTF, OSR, f0=f0)
		amp_pred = conv2seriesmatrix2D(amp_pred*1.0) #Ensure amplitudes are Float64
		vs_amp_predicted=(amp_pred, SNR_pred)
	end

	return (vs_amp_sim=(amp, SNR), peak=(pk_amp, pk_snr),
		vs_amp_predicted=vs_amp_predicted, OSR=OSR
	)
end

"""`calcSNRInfo(sig, dsm::RealDSM; NTF=nothing)`
"""
function calcSNRInfo(dsm::RealDSM; NTF=nothing)
	isnothing(NTF) && (NTF = synthesizeNTF(dsm))
	return _calcSNRInfo(NTF, dsm.OSR, dsm.M, dsm.f0)
end

#Last line


text:
using YAML, JLD2, HDF5, FileIO
using DataFrames
using Glob
using Base.Threads
import Agtor: generate_agparams

function load_yaml(data_dir::String, ext::String=".yml")::Dict{String, Dict}

    if !endswith(data_dir, ext)
        if startswith(ext, ".") == false
            throw(ArgumentError("Extension must include period, e.g. '.yml'"))
        end

        if endswith(data_dir, "/") == false || endswith(data_dir, "\\") == false
            data_dir *= '/'
        end

        data_files = glob("$(data_dir)" * "*$ext")
    else
        data_files = [data_dir]
    end

    return load_yaml(data_files)
end

function load_yaml(file_list::Array{String})::Dict{String, Dict}
    loaded_dataset::Dict{String, Dict} = Dict{String, Dict}()
    Threads.@threads for fn in file_list
        data = YAML.load(open(fn))
        loaded_dataset[data["name"]] = data
    end

    return loaded_dataset
end

function load_spec(data_dir::String)::Dict{Symbol, Dict}
    specs::Dict{String, Dict} = load_yaml(data_dir)
    return Dict(Symbol(k) => generate_agparams("", v) for (k,v) in specs)
end

function load_climate(data_fn::String)::Climate
    # Expect only CSV for now...
    if endswith(data_fn, ".csv")
        use_threads = Threads.nthreads() > 1
        climate_seq = DataFrame!(CSV.File(data_fn, threaded=use_threads, dateformat="YYYY-mm-dd"))
    else
        throw(ArgumentError("Currently, climate data can only be provided in CSV format"))
    end

    return Climate(climate_seq)
end

function determine_file_mode(fn::String)::String
    if isfile(fn)
        return "a+"
    end

    return "w"
end

"""Save state to JLD"""
function save_state!(fn::String, obj::Any, group::String)::Nothing
    mode = determine_file_mode(fn)
    jldopen(fn, mode) do file
        file[group] = obj
    end

    return nothing
end

"""Save arbitrary results to JLD."""
function save_results!(fn::String, results::Any, group::Tuple{String, String})::Nothing
    mode = determine_file_mode(fn)
    jldopen(fn, mode) do file
        idx, name = group
        file["$(idx)/$(name)"] = results
    end

    return nothing
end

"""Save arbitrary results to JLD."""
function save_results!(fn::String, results::NamedTuple)::Nothing
    mode = determine_file_mode(fn)
    jldopen(fn, mode) do file
        file["results"] = results
    end

    return nothing
end

"""Save results to JLD2.

# Example

```julia
output_fn = "example.jld2"
results = run_scenarios!(samples, z1, run_timestep!; post=allocation_callback!)
save_results!(output_fn, results)
```
"""
function save_results!(fn::String, results::Dict)::Nothing
    mode = determine_file_mode(fn)
    jldopen(fn, mode) do file
        for (i, res) in results
            file[i] = res
        end
    end

    return nothing
end

# """Save results to JLD2."""
# function save_results!(fn::String, results::Dict)::Nothing
#     mode = determine_file_mode(fn)
#     jldopen(fn, mode) do file
#         for (i, res) in results
#             if res isa NamedTuple
#                 # Handle Zone level results
#                 file["$(i)/zone_results"] = res.zone_results
#                 file["$(i)/field_results"] = res.field_results
#             elseif res isa Dict
#                 # handle Basin level results
#                 for (k, v) in res
#                     file["$(i)/$(k)/zone_results"] = v.zone_results
#                     file["$(i)/$(k)/field_results"] = v.field_results
#                 end
#             end
#         end
#     end

#     return nothing
# end

"""Save results for a single zone to JLD."""
function save_results!(fn::String, idx::Union{String, Int64}, results::NamedTuple)::Nothing
    mode = determine_file_mode(fn)

    try
        jldopen(fn, mode) do file
            file["$(idx)/zone_results"] = results.zone_results  # zone_results
            file["$(idx)/field_results"] = results.field_results  # field_results
        end
    catch e
        if isa(e, TaskFailedException)
            jldopen(fn, "a+") do file
                file["$(idx)/zone_results"] = results.zone_results  # zone_results
                file["$(idx)/field_results"] = results.field_results  # field_results
            end
        else
            println("Could not save results ($(idx)) to $(fn)!")
        end
    end

    return nothing
end

"""Save results for a scenario, basin, and zone combination"""
function save_results!(fn, sid, bid, zid, results; mode="w")::Nothing
    jldopen(fn, mode) do file
        prefix = "Scenario_$(sid)/$(bid)/$(zid)"
        file["$(prefix)/zone_results"] = results.zone_results  # zone_results
        file["$(prefix)/field_results"] = results.field_results  # field_results
    end

    return nothing
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
abstract type AbstractFittable end
abstract type AbstractEstimator <: AbstractFittable end
abstract type AbstractPipeline <: AbstractFittable end
abstract type AbstractTransformer <: AbstractFittable end

abstract type AbstractFeatureContrasts end
abstract type AbstractNonExistentFeatureContrasts <: AbstractFeatureContrasts end

abstract type AbstractNonExistentUnderlyingObject end

abstract type AbstractBackend end

abstract type AbstractPlot end


text:
module AtomixTests

include("test_core.jl")
include("test_sugar.jl")
include("test_doctest.jl")

if isdefined(Base, :replaceproperty!)  # 1.7 or later
    include("test_properties.jl")
end

end  # module AtomixTests


text:
# Utilities for making it easier
# to define solvers (like ODE solvers)
# for TDVP

"""
    to_vec(x)
Transform `x` into a `Vector`, and return the vector, and a closure which inverts the
transformation.

Modeled after FiniteDifferences.to_vec:

https://github.com/JuliaDiff/FiniteDifferences.jl/blob/main/src/to_vec.jl
"""
to_vec(x) = error("Not implemented")

function to_vec(x::ITensor)
  function ITensor_from_vec(x_vec)
    return itensor(x_vec, inds(x))
  end
  return vec(array(x)), ITensor_from_vec
end

# Represents a time-dependent sum of terms:
#
# H(t) = f[1](t) * H0[1] + f[2](t) * H0[2] + …
#
struct TimeDependentSum{S,T}
  f::Vector{S}
  H0::T
end
TimeDependentSum(f::Vector, H0::ProjMPOSum) = TimeDependentSum(f, H0.pm)
Base.length(H::TimeDependentSum) = length(H.f)

function Base.:*(c::Number, H::TimeDependentSum)
  return TimeDependentSum([t -> c * fₙ(t) for fₙ in H.f], H.H0)
end
Base.:*(H::TimeDependentSum, c::Number) = c * H

# Calling a `TimeDependentOpSum` at a certain time like:
#
# H(t)
#
# Returns a `ScaledSum` at that time.
(H::TimeDependentSum)(t::Number) = ScaledSum([fₙ(t) for fₙ in H.f], H.H0)

# Represents the sum of scaled terms:
#
# H = coefficient[1] * H[1] + coefficient * H[2] + …
#
struct ScaledSum{S,T}
  coefficients::Vector{S}
  H::T
end
Base.length(H::ScaledSum) = length(H.coefficients)

# Apply the scaled sum of terms:
#
# H(ψ₀) = coefficient[1] * H[1](ψ₀) + coefficient[2] * H[2](ψ₀) + …
#
# onto ψ₀.
function (H::ScaledSum)(ψ₀)
  ψ = ITensor(inds(ψ₀))
  for n in 1:length(H)
    ψ += H.coefficients[n] * apply(H.H[n], ψ₀)
  end
  return permute(ψ, inds(ψ₀))
end


text:
export color, chromatic_number

"""
`color(G,k)`: Return a `k`-coloring of `G` (or error if none exists).
"""
function color(G::SimpleGraph, k::Int)
    VV = vlist(G)
    EE = elist(G)
    n = NV(G)
    m = NE(G)

    if k<1
        error("Number of colors must be positive")
    end

    err_msg = "This graph is not " * string(k) * " colorable"

    # Special cases that don't require integer programming

    # see if we know the chromatic number already
    if cache_check(G,:chromatic_number)
        chi = cache_recall(G,:chromatic_number)
        if k<chi
            error(err_msg)
        end
    end

    result = Dict{vertex_type(G),Int}()

    if k==1
        if NE(G) > 0
            error(err_msg)
        end
        for v in VV
            result[v] = 1
        end
        return result
    end

    if k==2
        return two_color(G)
    end

    MOD = Model(get_solver())

    @variable(MOD, x[VV,1:k], Bin)

    for v in VV
        @constraint(MOD, sum(x[v,i] for i=1:k) == 1)
    end

    for e in EE
        u = e[1]
        v = e[2]
        for i=1:k
            @constraint(MOD, x[u,i] + x[v,i] <= 1)
        end
    end

    optimize!(MOD)
    status = Int(termination_status(MOD))

    if status != 1
        error(err_msg)
    end

    X = value.(x)

    for v in VV
        for c = 1:k
            if X[v,c] > 0
                result[v] = c
            end
        end
    end

    return result
end

function chromatic_number(G::SimpleGraph, verb::Bool=false)::Int
    if cache_check(G,:chromatic_number)
        return cache_recall(G,:chromatic_number)
    end

    if NV(G) == 0
        return 0
    end
    if NE(G) == 0
        return 1
    end

    # lower bound: larger of clique size or n/alpha
    n = NV(G)
    alf = length(max_indep_set(G))
    lb1 = Int(floor(n/alf))
    lb2 = length(max_clique(G))
    lb = max(lb1,lb2)

    # upper bound: try a random greedy coloring
    f = greedy_color(G)
    ub = maximum(values(f))

    k =  chromatic_number_work(G,lb,ub,verb)
    cache_save(G,:chromatic_number,k)
    return k
end

function chromatic_number_work(G::SimpleGraph, lb::Int, ub::Int, verb::Bool)::Int
    if verb
        print("$lb <= chi(G) <= $ub")
    end
    if lb == ub
        if verb
            println("\tand we're done")
        end
        return lb
    end

    mid = Int(floor((ub+lb)/2))

    if verb
        print("\tlooking for a $mid coloring")
    end

    try
        f = color(G,mid)  # success
        if verb
            println("\tfound")
        end
        return chromatic_number_work(G,lb,mid,verb)
    catch
    end
    if verb
        println("\tno such coloring")
    end
    return chromatic_number_work(G,mid+1,ub,verb)
end


text:
include("model.jl")

# particle marginal metropolis-hastings

@compiled @gen function var_proposal(prev)
    var_x::Float64 = get_assignment(prev)[:var_x]
    var_y::Float64 = get_assignment(prev)[:var_y]
    #@addr(normal(var_x, sqrt(0.15)), :var_x)
    #@addr(normal(var_y, sqrt(0.08)), :var_y)
    @addr(normal(var_x, sqrt(0.5)), :var_x)
    @addr(normal(var_y, sqrt(0.5)), :var_y)
end

@gen function observer(ys)
    for (i, y) in enumerate(ys)
        @addr(dirac(y), :hmm => :y => i)
    end
end

load_generated_functions()

##########
function strip_lineinfo(expr::Expr)
    @assert !(expr.head == :line)
    new_args = []
    for arg in expr.args
        if (isa(arg, Expr) && arg.head == :line) || isa(arg, LineNumberNode)
        elseif isa(arg, Expr) && arg.head == :block
            stripped = strip_lineinfo(arg)
            append!(new_args, stripped.args)
        else
            push!(new_args, strip_lineinfo(arg))
        end
    end
    Expr(expr.head, new_args...)
end

function strip_lineinfo(expr)
    expr
end

println("\n######################################################################\n")
obs = get_assignment(simulate(obs_sub, (1.2,)))
#println(strip_lineinfo(
    #Gen.codegen_generate(typeof(kernel), Tuple{Int,State,Params}, typeof(obs), Nothing)))
import InteractiveUtils
    InteractiveUtils.code_warntype(
        generate,
        (typeof(kernel), Tuple{Int,State,Params}, typeof(obs), Nothing))
println("\n######################################################################\n")

function initial_collapsed_trace(ys)
    T = length(ys)
    constraints = get_assignment(simulate(observer, (ys,)))
    (trace, weight) = generate(model_collapsed, (T,), constraints)
    trace
end

import Random
Random.seed!(1)

# generate synthetic dataset
T = 100 # was 500
(xs_sim, ys_sim) = hmm(10., 1., T)

# do inference
function do_inference(n)
    trace = initial_collapsed_trace(ys_sim)
    for iter=1:n
        score = get_call_record(trace).score
        println("score: $score")
        trace = mh(model_collapsed, var_proposal, (), trace)
        choices = get_assignment(trace)
	    println("var_x: $(choices[:var_x]), var_y: $(choices[:var_y])")
    end
end

import Profile

@time do_inference(2)
@time do_inference(100)

#Profile.@profile do_inference(17)

#Profile.print(format=:flat, sortedby=:count)
#Profile.print(mincount=10)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Statistics
using DataFrames

"""
Collate individually saved `.jld` files from distributed runs into a single JLD2 file.
"""
function collate_results!(fn_pattern::String, main_fn::String)::Nothing
    all_fns = glob(fn_pattern)
    jldopen(main_fn, "w") do file
        for fn in all_fns
            data = load(fn)
            for (grp_id, grp) in data
                file[grp_id] = grp
            end
        end
    end

    return nothing
end

"""
Collate a given list of names from across scenario runs into a individual DataFrames.
"""
function collate_results(fn::String, group::String, target::Array{String})::DataFrame
    @assert endswith(fn, ".jld2")
    x = load(fn)

    return collate_results(x, group, target)
end

"""
Collate a given result name from across scenario runs into a single DataFrame.
"""
function collate_results(fn::String, group::String, target::String)::DataFrame
    @assert endswith(fn, ".jld2")
    x = load(fn)

    return collate_results(x, group, target)
end

"""
Collate a given list of names from across scenario runs into individual DataFrames.
"""
function collate_results(data::Dict, group::String, target::Array{String})::DataFrame
    res = Dict()
    for tgt in target
        res[tgt] = collate_results(data, group, tgt)
    end

    return res
end

"""
Collate a given result name from across scenario runs into a single DataFrame.
"""
function collate_results(data::Dict, group::String, target::String)::DataFrame

    res = Dict()
    tgt_sym = Symbol(target)
    for (k, dataset) in data
        if !occursin(group, k)
            continue
        end

        if dataset isa NamedTuple
            for (fid, v) in dataset
                res[k*"_"*fid] = v[:, tgt_sym]
            end
        elseif dataset isa DataFrame
            res[k] = dataset[:, tgt_sym]
        elseif dataset isa Dict
            # collate basin level results
            compiled_results = for (idx, res) in dataset
                result_set = results_df(res, target; res_type="zone_results")
                mean(result_set[:total])
            end
            res[k] = compiled_results
        end
    end

    return DataFrame(res)
end

"""Filter DataFrame down to a subset of results based on partial column name match."""
function select_results(data::DataFrame, needle::String)::DataFrame
    return data[:, [c for c in names(data) if contains(String(c), needle)]]
end

"""Get aggregate statistics for an entire scenario run."""
function scenario_stats(data::DataFrame, needle::String)::NamedTuple
    df = select_results(data, needle)
    total = sum(sum.(eachrow(df)))
    avg = mean(mean.(eachrow(df)))
    med = median(median.(eachrow(df)))

    return (total=total, mean=avg, median=med)

end

function results_df(results::Union{Dict, String}, metric::String; res_type::String="zone_results")
    res_df = collate_results(results, res_type, metric)
    n_samples = ncol(res_df)
    res_stats = [scenario_stats(res_df, "$(i)/$(res_type)") for i in 1:n_samples]
    result_set = DataFrame(res_stats)

    return result_set
end


text:
# Uses new formulation from limbdark paper
using Limbdark
import Limbdark: Transit_Struct, s2

#include("s2_stable.jl")
#include("../../../src/transit_structure.jl")
include("../../../test/loglinspace.jl")
#include("../../../src/s2.jl")
#include("../../../src/define_constants.jl")

s2diff(b,r) = s2(r,b)-convert(Float64,s2(big(r),big(b)))

# Run some tests:

using PyPlot

# Found solution for using mathfrak here:
# https://github.com/JuliaPlots/Plots.jl/issues/997
#
using PyCall
PyDict(pyimport("matplotlib")["rcParams"])["text.latex.preamble"] = ["\\usepackage{eufrak}",]

nb = 1024; nr = 1024
fig,axes = subplots(3,2, figsize=(13, 8))

# Compute r = b+-\epsilon
epsilon = 1e-8
fepsp= zeros(Float64,nb)
fepsbigp= zeros(Float64,nb)
fepsm= zeros(Float64,nb)
fepsbigm= zeros(Float64,nb)
#b=range(epsilon,stop=2.0,length=nr)
b=linearspace(epsilon,2.0,nr)
for ib=1:nb
  fepsbigp[ib] = float(s2(big(b[ib]+epsilon),big(b[ib])))
  fepsp[ib] = float(s2(b[ib]+epsilon,b[ib]))
  fepsbigm[ib] = float(s2(big(b[ib]-epsilon),big(b[ib])))
  fepsm[ib] = float(s2(b[ib]-epsilon,b[ib]))
end

ax = axes[1]
ax[:plot](b,s2.(b,b),"-", lw=2,color="k",alpha=0.3,label=L"$r=b$")
ax[:plot](b,fepsp,"--", lw=2,color="C0",label=L"$r=b+\epsilon$")
ax[:plot](b,fepsm,":", lw=2,color="C1",label=L"$r=b-\epsilon$")
ax[:set_ylabel](L"$\mathfrak{s}_1$", fontsize=14)
ax[:set_title](L"$\mathfrak{s}_1(r = b \pm \epsilon,b)$", fontsize=16)
ax[:legend](loc="upper right",fontsize=10)

ax = axes[2]
ax[:plot](b,(fepsp-s2.(b,b))/epsilon,lw=2,color="C0",label=L"$\mathrm{double}$")
ax[:plot](b,(fepsbigp-s2.(b,b))/epsilon,linewidth=2,linestyle="--", color="C1",label=L"$\mathrm{double}$")
ax[:plot](b,(fepsm-s2.(b,b))/epsilon,lw=2,color="C0")
ax[:plot](b,(fepsbigm-s2.(b,b))/epsilon,linewidth=2,linestyle="--",color="C1")
ax[:set_ylabel](L"$\Delta \mathfrak{s}_1 / \epsilon$", fontsize=14)

ax[:annotate](L"$r = b - \epsilon$", xy=(1.8, 1.8), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0)

ax[:annotate](L"$r = b + \epsilon$", xy=(1.8, -1.8), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=0)

ax = axes[3]
ax[:plot](b,1e15 * s2diff.(b,b .+ epsilon),label=L"$r=b+\epsilon$")
ax[:plot](b,1e15 * s2diff.(b,b .- epsilon),linewidth=2,linestyle="--",label=L"$r=b-\epsilon$")
ax[:set_ylabel](L"$\mathrm{error}\ \times 10^{15}$", fontsize=14)
ax[:set_xlabel](L"$b$", fontsize=14)
ax[:legend](loc="lower right",fontsize=10)

# Compute b+r = 1+-\epsilon
epsilon = 1e-8
#b=range(0.0,stop=2.0,length=nr)
b=linearspace(0.0,2.0,nr)
for ib=1:nr
  fepsbigp[ib] = float(s2(big(1.0 .- b[ib] .+ epsilon),big(b[ib])))
  fepsp[ib] = float(s2(1.0 .- b[ib] .+ epsilon,b[ib]))
  fepsbigm[ib] = float(s2(big(1.0 .- b[ib] .- epsilon),big(b[ib])))
  fepsm[ib] = float(s2(1.0 .- b[ib] .- epsilon,b[ib]))
end

ax = axes[4]
ax[:plot](b,s2.(1.0 .- b,b),"k-",label=L"$r=1-b$",lw=2,alpha=0.3)
ax[:plot](b,fepsp,"--",label=L"$r=1-b+\epsilon$",lw=2,color="C0")
ax[:plot](b,fepsm,":",label=L"$r=1-b-\epsilon$",lw=2,color="C1")
ax[:legend](loc="lower right",fontsize=10)
ax[:set_title](L"$\mathfrak{s}_1(r = 1 - b \pm \epsilon,b)$", fontsize=16)

ax = axes[5]
ax[:plot](b,(fepsp-s2.(1.0 .- b,b))/epsilon,lw=2,label=L"$\mathrm{double}$", color="C0")
ax[:plot](b,(fepsbigp-s2.(1.0 .- b,b))/epsilon,lw=2,linestyle="--",label=L"$\mathrm{BigFloat}$", color="C1")
ax[:plot](b,(fepsm-s2.(1.0 .- b,b))/epsilon,lw=2, color="C0")
ax[:plot](b,(fepsbigm-s2.(1.0 .- b,b))/epsilon,lw=2,linestyle="--", color="C1")
ax[:legend](loc="lower right",fontsize=10)

ax[:annotate](L"$r = 1-b - \epsilon$", xy=(0.75, 1.5), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation=-35)

ax[:annotate](L"$r = 1-b + \epsilon$", xy=(0.75, -1.5), xycoords="data",
              xytext=(0, 0), textcoords="offset points",
              ha="center", va="center", rotation= 35)

ax = axes[6]
ax[:plot](b,1e15*s2diff.(1.0 .- b .+ epsilon,b),label=L"$r=1-b+\epsilon$")
ax[:plot](b,1e15*s2diff.(1.0 .- b .- epsilon,b),linewidth=2,linestyle="--",label=L"$r=1-b-\epsilon$")
ax[:legend](loc="lower right",fontsize=10)
ax[:set_xlabel](L"$b$", fontsize=14)

savefig("s2_machine.pdf", bbox_inches="tight")

max_issues_repo_name: rodluger/limbdark


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

using Markdown
using InteractiveUtils

# ╔═╡ 7da308ac-f282-11ea-3f24-95dc931f0456
using Plots

# ╔═╡ ea01ac60-f27d-11ea-13ee-f568175af295
md"""
# TP - F
Bon à savoir:
* où suis-je en train de travailler?
* où se trouvent mes fichiers?
* comment fonctionne mon ordi?
"""

# ╔═╡ 0db200ec-f27e-11ea-1950-a9d7e61afa63
pwd()

# ╔═╡ 5973e55e-f27e-11ea-353c-41dab383b09d
readdir(pwd())

# ╔═╡ 900ea8e2-f27e-11ea-3761-b108c219efc0
# path relatif
readdir("./data/")

# ╔═╡ a84c5c4c-f27e-11ea-304d-b54ef4bfdd80
# path absolu
readdir("/Users/bart/Documents/Stack/ES313/ES313.jl/Exercises/data/")

# ╔═╡ c11d3ce6-f27e-11ea-3e1d-bb7e6080bd82
# lire un fichier tct

# ╔═╡ 8d5bd6e0-f280-11ea-05a6-bbaa78cfaecd

# ╔═╡ 0cf36016-f27f-11ea-2165-bda0d94fa9da
function readrules(path::String)
	r = readlines(path)
	d = Dict()
	for line in r
		# split line
		intermed = parse.(Int,split(line,""))
		actuel = intermed[1]
		voisins = intermed[2:5]
		suivant = intermed[6]
		for i in 0:3
			get!(d, (actuel, circshift(voisins,i)), suivant)
		end
	end
	d
end

# ╔═╡ f11e8d0c-f281-11ea-32bd-1d9a3426aba8
d = readrules("./data/Langtonsrules.txt")

# ╔═╡ 6af82be0-f281-11ea-317d-d349315bfe91
d[(3,[0,0,0,7])]

# ╔═╡ 8a4a4332-f281-11ea-23bd-c35f9ef897c9
X = [1,2,3,4]

# ╔═╡ 92a12de8-f281-11ea-1d6c-97f1edfa3156
readlines("./data/Langtonstart.txt")

# ╔═╡ 07913c56-f282-11ea-0e39-b950958fbe77
begin
	r = readlines("./data/Langtonstart.txt")
	depart = permutedims(hcat([parse.(Int, split(line,"")) for line in r]...))
end

# ╔═╡ e7e5143e-f27e-11ea-3335-378eb2a0ffa1
val = r[173]

# ╔═╡ 804f6ffa-f282-11ea-1d17-ed258755405c
heatmap(depart, yflip=true)

# ╔═╡ 0f2dea90-f7ff-11ea-3a2a-3b4ff3176ff4
	mycmap = cgrad([RGBA(0/255,0/255,0/255),
    RGBA(0/255,0/255,255/255),
    RGBA(255/255,0/255,0/255),
    RGBA(0/255,255/255,0/255),
    RGBA(255/255,255/255,0/255),
    RGBA(255/255,0/255,255/255),
    RGBA(255/255,255/255,255/255),
    RGBA(0/255,255/255,255/255)])

# ╔═╡ a7902b40-f7eb-11ea-24e8-49b97db4e0f8
md"""
Choses à faire aujourd'hui:
* rectifier les couleurs
* incorporer la situation intiale dans un monde plus grand
* actualiser l'état
* générer une animation
"""

# ╔═╡ fb6968b0-f7fc-11ea-3b8c-9f8aeb89eb49
# créer le monde
function genworld(EI, dims)
	ssize = size(EI)
	#println(ssize .> dims)
	#println(typeof(ssize .> dims))
	if any(ssize .> dims)
		error("le monde est trop petit...")
	end
	# initialisation
	monde = zeros(Int,dims...)
	# incorporer EI
	i = round(Int,(dims[1] - ssize[1])/2) +1
	j = round(Int,(dims[2] - ssize[2])/2) +1
	monde[i:i+ssize[1]-1, j:j+ssize[2]-1] = EI

	return monde
end

# ╔═╡ faa26f2e-f7fc-11ea-0689-35a9fe9702ca
genworld(depart, (14,25))

# ╔═╡ fde2e240-f7fe-11ea-0ea2-377c6e99f379
heatmap(genworld(depart, (40,60)), yflip=true, color=mycmap, size=(250,200))

# ╔═╡ dd0c8d90-f7ff-11ea-0e50-d5eea9756d0c
# mise à jour
function update(monde, regles)
	neuf = zeros(Int, size(monde)...)
	for i in 2:size(monde,1)-1
		for j in 2:size(monde,2)-1
			neuf[i,j] = regles[(monde[i,j],[monde[i-1,j], monde[i,j+1], monde[i+1,j], monde[i,j-1]])]
		end
	end
	return neuf
end

# ╔═╡ 7645dd3e-f800-11ea-1dbd-55d7af40fb49
monde = genworld(depart, (12,20));

# ╔═╡ 8b03ec40-f800-11ea-1885-8b78f155da83
monde_v2 = update(monde,d)

# ╔═╡ ab35a710-f800-11ea-1e6f-8fdeffc51470
monde .== monde_v2

# ╔═╡ 03be2470-f801-11ea-153d-4fa68edc60f0
# combiner le tout
mutable struct Langton
	etat::Array{Int,2}
	regles::Dict
end

# ╔═╡ 236f210e-f802-11ea-1c0a-a96757d124ee
begin
	import Base.show
	function show(io,IO,L::Langton)
		print(io,"Langton instance de taille $(size(L.etat)) avec $(length(L.regles)) règles")
	end
end

# ╔═╡ 37fc5180-f801-11ea-0280-211d18e7a5c0
function langton(fei::String, fr::String, dims::Tuple=(12,20))
	regles = readrules(fr)
	r = readlines(fei)
	depart = permutedims(hcat([parse.(Int, split(line,"")) for line in r]...))
	monde = genworld(depart, dims)
	return Langton(monde, regles)
end

# ╔═╡ 9e6cd610-f801-11ea-3928-35eadc58071e
L = langton("./data/Langtonstart.txt","./data/Langtonsrules.txt")

# ╔═╡ 1fc5a110-f802-11ea-3546-ffd6cfd1ac18
L

# ╔═╡ 021d2930-f802-11ea-3e0c-673994053432
L

# ╔═╡ 260e5270-f801-11ea-03d7-27e618b8da3c
Langton(monde, d)

# ╔═╡ Cell order:
# ╟─ea01ac60-f27d-11ea-13ee-f568175af295
# ╠═0db200ec-f27e-11ea-1950-a9d7e61afa63
# ╠═5973e55e-f27e-11ea-353c-41dab383b09d
# ╠═900ea8e2-f27e-11ea-3761-b108c219efc0
# ╠═a84c5c4c-f27e-11ea-304d-b54ef4bfdd80
# ╠═c11d3ce6-f27e-11ea-3e1d-bb7e6080bd82
# ╟─e7e5143e-f27e-11ea-3335-378eb2a0ffa1
# ╠═8d5bd6e0-f280-11ea-05a6-bbaa78cfaecd
# ╠═0cf36016-f27f-11ea-2165-bda0d94fa9da
# ╠═f11e8d0c-f281-11ea-32bd-1d9a3426aba8
# ╠═6af82be0-f281-11ea-317d-d349315bfe91
# ╠═8a4a4332-f281-11ea-23bd-c35f9ef897c9
# ╠═92a12de8-f281-11ea-1d6c-97f1edfa3156
# ╠═07913c56-f282-11ea-0e39-b950958fbe77
# ╠═7da308ac-f282-11ea-3f24-95dc931f0456
# ╠═804f6ffa-f282-11ea-1d17-ed258755405c
# ╠═0f2dea90-f7ff-11ea-3a2a-3b4ff3176ff4
# ╟─a7902b40-f7eb-11ea-24e8-49b97db4e0f8
# ╠═fb6968b0-f7fc-11ea-3b8c-9f8aeb89eb49
# ╠═faa26f2e-f7fc-11ea-0689-35a9fe9702ca
# ╠═fde2e240-f7fe-11ea-0ea2-377c6e99f379
# ╠═dd0c8d90-f7ff-11ea-0e50-d5eea9756d0c
# ╠═7645dd3e-f800-11ea-1dbd-55d7af40fb49
# ╠═8b03ec40-f800-11ea-1885-8b78f155da83
# ╠═ab35a710-f800-11ea-1e6f-8fdeffc51470
# ╠═03be2470-f801-11ea-153d-4fa68edc60f0
# ╠═1fc5a110-f802-11ea-3546-ffd6cfd1ac18
# ╠═236f210e-f802-11ea-1c0a-a96757d124ee
# ╠═37fc5180-f801-11ea-0280-211d18e7a5c0
# ╠═021d2930-f802-11ea-3e0c-673994053432
# ╠═9e6cd610-f801-11ea-3928-35eadc58071e
# ╠═260e5270-f801-11ea-03d7-27e618b8da3c

max_stars_repo_path: Exercises/TPF (endstate).jl

max_issues_repo_path: Exercises/TPF (endstate).jl

max_forks_repo_path: Exercises/TPF (endstate).jl


text:
using Pkg

diffractorspec = PackageSpec(url="https://github.com/JuliaDiff/Diffractor.jl", rev="main")
Pkg.add([diffractorspec])
checkpointingspec = PackageSpec(path=joinpath(dirname(@__FILE__), ".."))
Pkg.develop(checkpointingspec)

# when first running instantiate
Pkg.instantiate()

using Documenter
using Checkpointing

makedocs(
    sitename = "Checkpointing.jl",
    format = Documenter.HTML(
        prettyurls = Base.get(ENV, "CI", nothing) == "true",
        mathengine = Documenter.KaTeX()
    ),
    modules = [Checkpointing],
    repo = "https://github.com/Argonne-National-Laboratory/Checkpointing.jl/blob/{commit}{path}#{line}",
    strict = true,
    checkdocs = :exports,
    pages = [
        "Home" => "index.md",
        "Quick Start" => "quickstart.md",
        "Library" => [
        "Checkpointing" => "lib/checkpointing.md",
        ],
    ]
)

deploydocs(
    repo = "github.com/Argonne-National-Laboratory/Checkpointing.jl.git",
    target = "build",
    devbranch = "main",
    devurl = "main",
    push_preview = true,
)


text:
function ObjectType(conn::Connection, name::AbstractString)
    obj_ref = Ref{Ptr{Cvoid}}()
    result = dpiConn_getObjectType(conn.handle, name, obj_ref)
    error_check(context(conn), result)
    return ObjectType(conn, obj_ref[])
end

object_type_schema(obj::ObjectType) = unsafe_string(obj.type_info.schema, obj.type_info.schema_length)
object_type_name(obj::ObjectType) = unsafe_string(obj.type_info.name, obj.type_info.name_length)
object_type_num_attributes(obj::ObjectType) = Int(obj.type_info.num_attributes)
object_type_is_collection(obj::ObjectType) = obj.type_info.is_collection != 0

function object_type_element_type_info(obj::ObjectType) :: Union{Nothing, OraDataTypeInfo}
    if object_type_is_collection(obj)
        return obj.element_type_info
    else
        return nothing
    end
end


text:
const errmax = 1e-15

#Compute tail of noncentral Beta distribution
#Russell Lenth, Algorithm AS 226: Computing Noncentral Beta Probabilities,
#Applied Statistics,Volume 36, Number 2, 1987, pages 241-244

"""
	ncbeta_tail(x,a,b,lambda)

Compute tail of the noncentral beta distribution.
Uses the recursive relation
```math
I_{x}(a,b+1;0) = I_{x}(a,b;0) - \\Gamma(a+b)/\\Gamma(a+1)\\Gamma(b)x^{a}(1-x)^{b}
```
and ``\\Gamma(a+1) = a\\Gamma(a)`` given in https://dlmf.nist.gov/8.17.21.
"""
function ncbeta_tail(a::Float64, b::Float64, lambda::Float64, x::Float64)
    if x <= 0.0
        return 0.0
    elseif x >= 1.0
        return 1.0
    end

    c = 0.5*lambda
    #Init series

    beta = logabsbeta(a,b)[1]
    temp = beta_inc(a,b,x)[1]
    gx = (beta_integrand(a,b,x,1.0-x))/a
    q = exp(-c)
    xj = 0.0
    ax = q*temp
    sumq = 1.0 - q
    ans = ax

    while true
        xj += 1.0
        temp -= gx
        gx *= x*(a+b+xj-1.0)/(a+xj)
        q *= c/xj
        sumq -= q
        ax = temp*q
        ans += ax

        #Check convergence
        errbd = abs((temp-gx)*sumq)
        if xj > 1000 || errbd < 1e-10
            break
        end
    end
    return ans
end

"""
    ncbeta_poisson(a,b,lambda,x)

Compute CDF of noncentral beta if lambda >= 54 using:
First ``\\lambda/2`` is calculated and the Poisson term is calculated using ``P(j-1)=j/\\lambda P(j)`` and ``P(j+1) = \\lambda/(j+1) P(j)``.
Then backward recurrences are used until either the Poisson weights fall below `errmax` or `iterlo` is reached.
```math
I_{x}(a+j-1,b) = I_{x}(a+j,b) + \\Gamma(a+b+j-1)/\\Gamma(a+j)\\Gamma(b)x^{a+j-1}(1-x)^{b}
```
Then forward recurrences are used until error bound falls below `errmax`.
```math
I_{x}(a+j+1,b) = I_{x}(a+j,b) - \\Gamma(a+b+j)/\\Gamma(a+j)\\Gamma(b)x^{a+j}(1-x)^{b}
```
"""
function ncbeta_poisson(a::Float64, b::Float64, lambda::Float64, x::Float64)
    c = 0.5*lambda
    xj = 0.0
    m = round(Int, c)
    mr = float(m)
    iterlo = m - trunc(Int, 5.0*sqrt(mr))
    iterhi = m + trunc(Int, 5.0*sqrt(mr))
    t = -c + mr*log(c) - logabsgamma(mr + 1.0)[1]
    q = exp(t)
    r = q
    psum = q

    beta = logabsbeta(a+mr,b)[1]
    gx = beta_integrand(a+mr,b,x,1.0-x)/(a + mr)
    fx = gx
    temp = beta_inc(a+mr,b,x)[1]
    ftemp = temp
    xj += 1.0

    sm = q*temp
    iter1 = m

    #Iterations start from M and goes downwards

    for iter1 = m:-1:iterlo
        if q < errmax
            break
        end

        q *= iter1/c
        xj += 1.0
        gx *= (a + iter1)/(x*(a+b+iter1-1.0))
        iter1 -= 1
        temp += gx
        psum += q
        sm += q*temp
    end

    t0 = logabsgamma(a+b)[1] - logabsgamma(a+1.0)[1] - logabsgamma(b)[1]
    s0 = a*log(x) + b*log1p(-x)

    s = 0.0
    for j = 0:iter1-1
        s += exp(t0+s0+j*log(x))
        t1 = log(a+b+j) - log(a+j+1.0) + t0
        t0 = t1
    end
    #Compute first part of error bound

    errbd = (gamma_inc(float(iter1),c,0)[2])*(temp+s)
    q = r
    temp = ftemp
    gx = fx
    iter2 = m
    #Iterations for the higher part

    for iter2 = m:iterhi-1
        ebd = errbd + (1.0 - psum)*temp
        if ebd < errmax
            return sm
        end
        iter2 += 1
        xj += 1.0
        q *= c/iter2
        psum += q
        temp -= gx
        gx *= x*(a+b+iter2-1.0)/(a+iter2)
        sm += q*temp
    end
    return sm
end

#R Chattamvelli, R Shanmugam, Algorithm AS 310: Computing the Non-central Beta Distribution Function,
#Applied Statistics, Volume 46, Number 1, 1997, pages 146-156

"""
	ncbeta(a,b,lambda,x)

Compute the CDF of the noncentral beta distribution given by
```math
I_{x}(a,b;\\lambda ) = \\sum_{j=0}^{\\infty}q(\\lambda/2,j)I_{x}(a+j,b;0)
```
For ``\\lambda < 54`` : algorithm suggested by Lenth(1987) in `ncbeta_tail(a,b,lambda,x)`.
Else for ``\\lambda >= 54`` : modification in Chattamvelli(1997) in `ncbeta_poisson(a,b,lambda,x)` by using both forward and backward recurrences.
"""
function ncbeta(a::Float64, b::Float64, lambda::Float64, x::Float64)
    ans = x
    if x <= 0.0
        return 0.0
    elseif x >= 1.0
        return 1.0
    end

    if lambda < 54.0
        return ncbeta_tail(a,b,lambda,x)
    else
        return ncbeta_poisson(a,b,lambda,x)
    end
end

"""
    ncF(x,v1,v2,lambda)

Compute CDF of noncentral F distribution given by:
```math
F(x, v1, v2; lambda) = I_{v1*x/(v1*x + v2)}(v1/2, v2/2; \\lambda)
```
where ``I_{x}(a,b; lambda)`` is the noncentral beta function computed above.

Wikipedia: https://en.wikipedia.org/wiki/Noncentral_F-distribution
"""
function ncF(x::Float64, v1::Float64, v2::Float64, lambda::Float64)
    return ncbeta(v1/2, v2/2, lambda, (v1*x)/(v1*x + v2))
end

function ncbeta(a::T,b::T,lambda::T,x::T) where {T<:Union{Float16,Float32}}
	T.(ncbeta(Float64(a),Float64(b),Float64(lambda),Float64(x)))
end

function ncF(x::T,v1::T,v2::T,lambda::T) where {T<:Union{Float16,Float32}}
	T.(ncF(Float64(x),Float64(v1),Float64(v2),Float64(lambda)))
end

ncbeta(a::Real,b::Real,lambda::Real,x::Real) = ncbeta(promote(float(a),float(b),float(lambda),float(x))...)
ncbeta(a::T,b::T,lambda::T,x::T) where {T<:AbstractFloat} = throw(MethodError(ncbeta,(a,b,lambda,x,"")))
ncF(x::Real,v1::Real,v2::Real,lambda::Real) = ncF(promote(float(x),float(v1),float(v2),float(lambda))...)
ncF(x::T,v1::T,v2::T,lambda::T) where {T<:AbstractFloat} = throw(MethodError(ncF,(x,v1,v2,lambda,"")))


text:
function call_bbox(bb_methods,all_bounds,func_evals,j,func,max_time)

    opt1 = bbsetup(@eval $func; Method=bb_methods[j], SearchRange = all_bounds, TraceMode = :silent, MaxFuncEvals = func_evals)# , MaxTime = max_time
    res1 = bboptimize(opt1)
    opti1=res1.archive_output.best_fitness
    sol1=res1.archive_output.best_candidate
    hist_num = length(opt1.runcontrollers[1].evaluator.archive.fitness_history)
    hist_tmp=zeros(hist_num)
    for i=1:hist_num
            hist_tmp[i]=opt1.runcontrollers[1].evaluator.archive.fitness_history[i].fitness
    end
    matt=repeat(hist_tmp', trunc(Int, func_evals/hist_num)+1, 1)[:]
    matt[end]=copy(opti1)

    return matt,opti1,sol1
end

#


text:
module TestSweep

using TimeZoneLookup
using TimeZoneLookup: V
using Test

@testset "Points comparison" begin
    @test V(1, 2) > V(1, 3)
    @test V(1, 2) < V(2, 2)
end

end


text:
using POMDPs
using LaserTag
using ParticleFilters
using ProfileView
using QMDP
using POMDPSimulators
using Profile

p = gen_lasertag()

solver = QMDPSolver(max_iterations=1000, verbose=true)
@time pol = solve(solver, p)

# @show value(pol, initialstate_distribution(p))

Profile.clear()
@profile solve(solver, p)
ProfileView.view()

# sim = HistoryRecorder(max_steps=5, rng=MersenneTwister(2), show_progress=true)
#
# filter = SIRParticleFilter(p, 10000)
#
# #=
# b = initialize_belief(filter, initial_state_distribution(p))
# o = rand(MersenneTwister(3), observation(p, initial_state(p, MersenneTwister(4))))
# @code_warntype update(filter, b, 1, o)
# =#
#
# simulate(sim, p, pol, filter)
#
# # @time simulate(sim, p, pol, filter)
# # simulate(sim, p, pol, filter)
#
# Profile.clear()
# @profile simulate(sim, p, pol, filter)
#
# ProfileView.view()