text:
using Tables
using DataAPI
using Thrift
using Snappy
using CodecZstd: ZstdCompressor
using CodecZlib: GzipCompressor
#using CodecLz4: LZ4HCCompressor # wating for CodecLz4.jl devs to fix a bug
using LittleEndianBase128
using Base.Iterators: partition
using CategoricalArrays: CategoricalArray, CategoricalValue

using Base: SkipMissing

if VERSION < v"1.3"
    using Missings: nonmissingtype
end

# a mapping of Julia types to _Type codes in Parquet format
const COL_TYPE_CODE = Dict{DataType, Int32}(
    Bool => PAR2._Type.BOOLEAN,
    Int32 => PAR2._Type.INT32,
    Int64 => PAR2._Type.INT64,
    #INT96 => 3,  // deprecated, only used by legacy implementations. # not supported by Parquet.jl
    Float32 => PAR2._Type.FLOAT,
    Float64 => PAR2._Type.DOUBLE,
    String => PAR2._Type.BYTE_ARRAY, # BYTE_ARRAY
    # FIXED_LEN_BYTE_ARRAY => 7, # current there is no Julia type that we support that maps to this type
    )

function write_thrift(fileio, thrift_obj)
    """write thrift definition to file"""
    pos_before_write = position(fileio)
    p = TCompactProtocol(TFileTransport(fileio))
    Thrift.write(p, thrift_obj)
    pos_after_write = position(fileio)

    size_of_written = pos_after_write - pos_before_write

    size_of_written
end

function compress_using_codec(colvals::AbstractArray, codec::Integer)::Vector{UInt8}
    """Compress `isbits` column types using codec"""
    uncompressed_byte_data = reinterpret(UInt8, colvals) |> collect

    if codec == PAR2.CompressionCodec.UNCOMPRESSED
        return uncompressed_byte_data
    elseif codec == PAR2.CompressionCodec.SNAPPY
        compressed_data = Snappy.compress(uncompressed_byte_data)
    elseif codec == PAR2.CompressionCodec.GZIP
        compressed_data = transcode(GzipCompressor, uncompressed_byte_data)
    elseif codec == PAR2.CompressionCodec.LZ4
        error("lz4 is not supported as data compressed with https://github.com/JuliaIO/CodecLz4.jl can't seem to be read by R or Python. If you know how to fix it please help out.")
        #compressed_data = transcode(LZ4HCCompressor, uncompressed_byte_data)
    elseif codec == PAR2.CompressionCodec.ZSTD
        compressed_data = transcode(ZstdCompressor, uncompressed_byte_data)
    else
        error("not yet implemented")
    end

    return compressed_data
end

function compress_using_codec(colvals::AbstractVector{String}, codec::Int)::Vector{UInt8}
    """Compress `String` column using codec"""
    # the output
    io = IOBuffer()

    # write the values
    for val in colvals
        # for string it needs to be stored as BYTE_ARRAY which needs the length
        # to be the first 4 bytes UInt32
        write(io, val |> sizeof |> UInt32 |> htol)
        # write each of the strings one after another
        write(io, val)
    end

    uncompressed_bytes = take!(io)
    return compress_using_codec(uncompressed_bytes, codec)
end

function write_defn_levels(data_to_compress_io, colvals::AbstractVector{Union{Missing, T}}) where T
    """ A function to write definition levels for `Union{Missing, T}`"""
    # if there is missing
    # use the bit packing algorithm to write the
    # definition_levels
    bytes_needed = ceil(Int, length(colvals) / 8sizeof(UInt8))
    tmp = UInt32((UInt32(bytes_needed) << 1) | 1)
    bitpacking_header = LittleEndianBase128.encode(tmp)

    tmpio = IOBuffer()
    not_missing_bits::BitArray = .!ismissing.(colvals)
    write(tmpio, not_missing_bits)
    seek(tmpio, 0)

    encoded_defn_data = read(tmpio, bytes_needed)

    encoded_defn_data_length = length(bitpacking_header) + bytes_needed
    # write the definition data
    write(data_to_compress_io, UInt32(encoded_defn_data_length) |> htol)
    write(data_to_compress_io, bitpacking_header)
    write(data_to_compress_io, encoded_defn_data)
end

function write_defn_levels(data_to_compress_io, colvals::AbstractVector)
    """ A function to write definition levels for NON-missing data
    """
    # if there is no missing can just use RLE of one
    # using rle
    rle_header = LittleEndianBase128.encode(UInt32(length(colvals)) << 1)
    repeated_value = UInt8(1)
    encoded_defn_data_length = sizeof(rle_header) + sizeof(repeated_value)

    # write the definition data
    write(data_to_compress_io, UInt32(encoded_defn_data_length) |> htol)
    write(data_to_compress_io, rle_header)
    write(data_to_compress_io, repeated_value)
end

# TODO turn this on when writing dictionary is necessary
# function write_col_dict(fileio, colvals::AbstractArray{T}, codec) where T
#     """ write the column dictionary page """
#     # note: `level`s does not return `missing` as a level
#     uvals = DataAPI.levels(colvals)
#
#     # do not support dictionary with more than 127 levels
#     # TODO relax this 127 restriction
#     if length(uvals) > 127
#         @warn "More than 127 levels in dictionary. Parquet.jl does not support this at this stage."
#         return (offset = missing, uncompressed_size = 0, compressed_size = 0)
#     end
#
#     if nonmissingtype(T) == String
#         # the raw bytes of made of on UInt32 to indicate string length
#         # and the content of the string
#         # so the formula for dict size is as below
#         uncompressed_dict_size = sizeof(UInt32)*length(uvals) + sum(sizeof, uvals)
#     else
#         uncompressed_dict_size = length(uvals)*sizeof(eltype(uvals))
#     end
#
#     compressed_uvals::Vector{UInt8} = compress_using_codec(uvals, codec)
#     compressed_dict_size = length(compressed_uvals)
#
#     # TODO do the CRC properly
#     crc = 0
#
#     # construct dictionary metadata
#     dict_page_header = PAR2.PageHeader()
#
#     dict_page_header._type = PAR2.PageType.DICTIONARY_PAGE
#     dict_page_header.uncompressed_page_size = uncompressed_dict_size
#     dict_page_header.compressed_page_size = compressed_dict_size
#     dict_page_header.crc = crc
#
#     dict_page_header.dictionary_page_header = PAR2.DictionaryPageHeader()
#     dict_page_header.dictionary_page_header.num_values = Int32(length(uvals))
#     dict_page_header.dictionary_page_header.encoding = PAR2.Encoding.PLAIN_DICTIONARY
#     dict_page_header.dictionary_page_header.is_sorted = false
#
#     before_write_page_header_pos = position(fileio)
#
#     dict_page_header_size = write_thrift(fileio, dict_page_header)
#
#     # write the dictionary data
#     write(fileio, compressed_uvals)
#
#     return (offset = before_write_page_header_pos, uncompressed_size = uncompressed_dict_size + dict_page_header_size, compressed_size = compressed_dict_size + dict_page_header_size)
# end

write_encoded_data(data_to_compress_io, colvals::AbstractVector{Union{Missing, T}}) where T =
    write_encoded_data(data_to_compress_io, skipmissing(colvals))

function write_encoded_data(data_to_compress_io, colvals::Union{AbstractVector{String}, SkipMissing{S}}) where S <: AbstractVector{Union{Missing, String}}
    """ Write encoded data for String type """
    # write the values
    for val in colvals
        # for string it needs to be stored as BYTE_ARRAY which needs the length
        # to be the first 4 bytes UInt32
        write(data_to_compress_io, val |> sizeof |> UInt32 |> htol)
        # write each of the strings one after another
        write(data_to_compress_io, val)
    end
end

function write_encoded_data(data_to_compress_io, colvals::Union{AbstractVector{Bool}, SkipMissing{S}}) where S <: AbstractVector{Union{Missing, Bool}}
    """ Write encoded data for Bool type """
    # write the bitacpked bits
    # write a bitarray seems to write 8 bytes at a time
    # so write to a tmpio first
    no_missing_bit_vec =  BitArray(colvals)
    bytes_needed = ceil(Int, length(no_missing_bit_vec) / 8sizeof(UInt8))
    tmpio = IOBuffer()
    write(tmpio, no_missing_bit_vec)
    seek(tmpio, 0)
    packed_bits = read(tmpio, bytes_needed)
    write(data_to_compress_io, packed_bits)
end

function write_encoded_data(data_to_compress_io, colvals::AbstractArray)
    """ Efficient write of encoded data for `isbits` types"""
    @assert isbitstype(eltype(colvals))
    write(data_to_compress_io, colvals |> htol)
end

function write_encoded_data(data_to_compress_io, colvals::SkipMissing)
    """ Write of encoded data for skipped missing types"""
    for val in colvals
        write(data_to_compress_io, val |> htol)
    end
end

function write_encoded_data(data_to_compress_io, colvals)
    """ Write of encoded data for the most general type.
    The only requirement is that colvals has to be iterable
    """
    for val in skipmissing(colvals)
        write(data_to_compress_io, val |> htol)
    end
end

# TODO set the encoding code into a dictionary
function write_col_page(fileio, colvals::AbstractArray, codec, ::Val{PAR2.Encoding.PLAIN})
    """
    Write a chunk of data into a data page using PLAIN encoding where the values
    are written back-to-back in memory and then compressed with the codec.
    For `String`s, the values are written with length (UInt32), followed by
    content; it is NOT null terminated.
    """

    # generate the data page header
    data_page_header = PAR2.PageHeader()

    # set up an IO buffer to write to
    data_to_compress_io = IOBuffer()

    # write repetition level data
    ## do nothing
    ## this seems to be related to nested columns
    ## and hence is not needed here as we only supported unnested column write

    # write definition levels
    write_defn_levels(data_to_compress_io, colvals)

    # write the encoded data
    write_encoded_data(data_to_compress_io, colvals)

    data_to_compress::Vector{UInt8} = take!(data_to_compress_io)

    compressed_data::Vector{UInt8} = compress_using_codec(data_to_compress, codec)

    uncompressed_page_size = length(data_to_compress)
    compressed_page_size = length(compressed_data)

    data_page_header._type = PAR2.PageType.DATA_PAGE
    data_page_header.uncompressed_page_size = uncompressed_page_size
    data_page_header.compressed_page_size = compressed_page_size

    # TODO proper CRC
    data_page_header.crc = 0

    data_page_header.data_page_header = PAR2.DataPageHeader()
    data_page_header.data_page_header.num_values = Int32(length(colvals))
    data_page_header.data_page_header.encoding = PAR2.Encoding.PLAIN
    data_page_header.data_page_header.definition_level_encoding = PAR2.Encoding.RLE
    data_page_header.data_page_header.repetition_level_encoding = PAR2.Encoding.RLE

    position_before_page_header_write = position(fileio)

    size_of_page_header_defn_repn = write_thrift(fileio, data_page_header)

    # write data
    write(fileio, compressed_data)

    return (
        offset = position_before_page_header_write,
        uncompressed_size = uncompressed_page_size + size_of_page_header_defn_repn,
        compressed_size = compressed_page_size + size_of_page_header_defn_repn,
    )
end

function write_col_page(fileio, colvals::AbstractArray, codec, ::Val{PAR2.Encoding.PLAIN_DICTIONARY})
    """write Dictionary encoding data page"""
    error("PLAIN_DICTIONARY encoding not implemented yet")

    # TODO finish the implementation
    rle_header = LittleEndianBase128.encode(UInt32(length(colvals)) << 1)
    repeated_value = UInt8(1)

    encoded_defn_data_length = sizeof(rle_header) + sizeof(repeated_value)

    ## write the encoded data length
    write(fileio, encoded_defn_data_length |> UInt32 |> htol)

    write(fileio, rle_header)
    write(fileio, repeated_value)

    position(fileio)

    # write the data

    ## firstly, bit pack it

    # the bitwidth to use
    bitwidth = ceil(UInt8, log(2, length(uvals)))
    # the max bitwidth is 32 according to documentation
    @assert bitwidth <= 32
    # to do that I have to figure out the Dictionary index of it
    # build a JuliaDict
    val_index_dict = Dict(zip(uvals, 1:length(uvals)))

    bitwidth_mask = UInt32(2^bitwidth-1)

    bytes_needed = ceil(Int, bitwidth*length(colvals) / 8)

    bit_packed_encoded_data = zeros(UInt8, bytes_needed)
    upto_byte = 1

    bits_written = 0
    bitsz = 8sizeof(UInt8)

    for val in colvals
        bit_packed_val = UInt32(val_index_dict[val]) & bitwidth_mask
        if bitwidth_mask <= bitsz - bits_written
            bit_packed_encoded_data[upto_byte] = (bit_packed_encoded_data[upto_byte] << bitwidth_mask) | bit_packed_val
        else
            # this must mean
            # bitwidth_mask > bitsz - bits_written
            # if the remaining bits is not enough to write a packed number
            42
        end
    end
end

function write_col_page(fileio, colvals::AbstractArray{T}, codec, encoding) where T
    error("Page encoding $encoding is yet not implemented.")
end

write_col(fileio, colvals::CategoricalArray, args...; kwars...) = begin
    throw("Currently CategoricalArrays are not supported.")
end

function write_col(fileio, colvals::AbstractArray{T}, colname, encoding, codec; nchunks = 1) where T
    """Write a column to a file"""
    # TODO turn writing dictionary on
    # Currently, writing the dictionary page is not turned on for any type.
    # Normally, for Boolean data, dictionary is not supported. However for other
    # data types, dictionary page CAN be supported. However, since Parquet.jl
    # only supports writing PLAIN encoding data, hence there is no need to write
    # a dictionary page until other dictionary-based encodings are supported
    dict_info = (offset = missing, uncompressed_size = 0, compressed_size = 0)

    num_vals_per_chunk = ceil(Int, length(colvals) / nchunks)

    chunk_info = [write_col_page(fileio, val_chunk, codec, Val(encoding)) for val_chunk in partition(colvals, num_vals_per_chunk)]

    sizes = reduce(chunk_info; init = dict_info) do x, y
        (
            uncompressed_size = x.uncompressed_size + y.uncompressed_size,
            compressed_size = x.compressed_size + y.compressed_size
        )
    end

    # write the column metadata
    # can probably write the metadata right after the data chunks
    col_meta = PAR2.ColumnMetaData()

    col_meta._type = COL_TYPE_CODE[eltype(colvals) |> nonmissingtype]
    # these are all the fields
    # TODO collect all the encodings used
    if eltype(colvals) == Bool
        col_meta.encodings = Int32[0, 3]
    else
        col_meta.encodings = Int32[2, 0, 3]
    end
    col_meta.path_in_schema = [colname]
    col_meta.codec = codec
    col_meta.num_values = length(colvals)

    col_meta.total_uncompressed_size = sizes.uncompressed_size
    col_meta.total_compressed_size = sizes.compressed_size

    col_meta.data_page_offset = chunk_info[1].offset
    if !ismissing(dict_info.offset)
        col_meta.dictionary_page_offset = dict_info.offset
    end

    # write the column meta data right after the data
    # keep track of the position so it can put into the column chunk
    # metadata
    col_meta_offset = position(fileio)
    write_thrift(fileio, col_meta)

    # Prep metadata for the filemetadata
    ## column chunk metadata
    col_chunk_meta = PAR2.ColumnChunk()

    col_chunk_meta.file_offset = col_meta_offset
    col_chunk_meta.meta_data = col_meta

    return (
        data_page_offset = chunk_info[1].offset,
        dictionary_page_offset =  dict_info.offset,
        col_chunk_meta = col_chunk_meta,
        col_meta_offset = col_meta_offset
    )
end

function create_schema_parent_node(ncols)
    """Create the parent node in the schema tree"""
    schmea_parent_node = PAR2.SchemaElement()
    schmea_parent_node.name = "schema"
    schmea_parent_node.num_children = ncols
    schmea_parent_node
end

function create_col_schema(type, colname)
    """Create a column node in the schema tree for non-strings"""
    schema_node = PAR2.SchemaElement()
    # look up type code
    schema_node._type = COL_TYPE_CODE[type |> nonmissingtype]
    schema_node.repetition_type = 1
    schema_node.name = colname
    schema_node.num_children = 0

    schema_node
end

function create_col_schema(type::Type{String}, colname)
    """create col schema for string"""
    schema_node = PAR2.SchemaElement()
    # look up type code
    schema_node._type = COL_TYPE_CODE[type]
    schema_node.repetition_type = 1
    schema_node.name = colname
    schema_node.num_children = 0

    # for string set converted type to UTF8
    schema_node.converted_type = PAR2.ConvertedType.UTF8

    logicalType = PAR2.LogicalType()
    logicalType.STRING = PAR2.StringType()

    schema_node.logicalType = logicalType

    schema_node
end

"""
    Write a parquet file from a Tables.jl compatible table e.g DataFrame

io                  -   A writable IO stream
tbl                 -   A Tables.jl columnaccessible table e.g. a DataFrame
compression_code    -   Default "SNAPPY". The compression codec. The supported
                        values are "UNCOMPRESSED", "SNAPPY", "ZSTD", "GZIP"
"""
function write_parquet(io::IO, x; compression_codec = "SNAPPY")
    tbl = Tables.columns(x)

    # check that all types are supported
    sch = Tables.schema(tbl)
    err_msgs = String[]
    for type in sch.types
        if type <: CategoricalValue
            push!(err_msgs, "CategoricalArrays are not supported at this stage. \n")
        elseif !(nonmissingtype(type) <: Union{Int32, Int64, Float32, Float64, Bool, String})
            push!(err_msgs, "Column whose `eltype` is $type is not supported at this stage. \n")
        end
    end

    err_msgs = unique(err_msgs)
    if length(err_msgs) > 0
        throw(reduce(*, err_msgs))
    end

    # set the data page encoding
    # currently only PLAIN is supported
    # TODO add support for other encodings see
    # https://github.com/apache/parquet-format/blob/master/Encodings.md
    encoding = Encoding.PLAIN

    # convert a string or symbol compression codec into the numeric code
    codec = getproperty(PAR2.CompressionCodec, Symbol(uppercase(string(compression_codec))))

    # figure out the right number of chunks
    # TODO test that it works for all supported table
    nrows = Tables.rowcount(tbl)
    sample_size = min(100, nrows)
    rs = collect(Iterators.take(Tables.namedtupleiterator(tbl), sample_size))
    table_size_bytes = Base.summarysize(rs) / sample_size * nrows

    approx_raw_to_parquet_compression_ratio = 6
    approx_post_compression_size = (table_size_bytes / 2^30) / approx_raw_to_parquet_compression_ratio

    # if size is larger than 64mb and has more than 6 rows
    if (approx_post_compression_size > 0.064) & (nrows > 6)
        recommended_chunks = ceil(Int, approx_post_compression_size / 6) * 6
    else
        recommended_chunks = 1
    end

    colnames = String.(Tables.columnnames(tbl))
    _write_parquet(
        io,
        Tables.Columns(tbl),
        Tables.columnnames(tbl),
        recommended_chunks;
        encoding    =   Dict(col => encoding for col in colnames),
        codec       =   Dict(col => codec for col in colnames)
    )
end

"""
    Write a parquet file from a Tables.jl compatible table e.g DataFrame

path                -   The file path
tbl                 -   A Tables.jl columnaccessible table e.g. a DataFrame
compression_code    -   Default "SNAPPY". The compression codec. The supported
                        values are "UNCOMPRESSED", "SNAPPY", "ZSTD", "GZIP"
"""
function write_parquet(path, x; compression_codec = "SNAPPY")
    open(path, "w") do io
        write_parquet(io, x; compression_codec=compression_codec)
    end
end

function _write_parquet(io::IO, itr_vectors, colnames, nchunks; ncols = length(itr_vectors), encoding::Dict{String, Int32}, codec::Dict{String, Int32})

    """Internal method for writing parquet

    itr_vectors -   An iterable of `AbstractVector`s containing the values to be
                    written
    colnames    -   Column names for each of the vectors
    path        -   The output parquet file path
    nchunks     -   The number of chunks/pages to write for each column
    ncols       -   The number of columns. This is provided as an argument for
                    the case where the `length(itr_vectors)` is not defined,
                    e.g. lazy loading of remote resources.
    encoding    -   A dictionary mapping from column names to encoding
    codec       -   A dictionary mapping from column names to compression codec
    """
    write(io, "PAR1")

    # the + 1 comes from the fact that schema is a tree and there is an extra
    # parent node
    schemas = Vector{PAR2.SchemaElement}(undef, ncols + 1)
    schemas[1] = create_schema_parent_node(ncols)
    col_chunk_metas = Vector{PAR2.ColumnChunk}(undef, ncols)
    row_group_file_offset = missing

    # write the columns one by one
    # TODO parallelize this
    nrows = -1 # initialize it
    for (coli, (colname_sym, colvals)) in enumerate(zip(colnames, itr_vectors))
        colname = String(colname_sym)

        col_encoding = encoding[colname]
        col_codec = codec[colname]
        # write the data including metadata
        col_info = write_col(io, colvals, colname, col_encoding, col_codec; nchunks = nchunks)

        # the `row_group_file_offset` keeps track of where the data starts, so
        # keep it at the dictonary of the first data
        if coli == 1
            nrows = length(colvals)
            if ismissing(col_info.dictionary_page_offset)
                row_group_file_offset = col_info.data_page_offset
            else
                row_group_file_offset = col_info.dictionary_page_offset
            end
        end

        col_chunk_metas[coli] = col_info.col_chunk_meta

        # add the schema
        schemas[coli + 1] = create_col_schema(eltype(colvals) |> nonmissingtype, colname)
    end

    # now all the data is written we write the filemetadata
    # finalise it by writing the filemetadata
    filemetadata = PAR2.FileMetaData()
    filemetadata.version = 1
    filemetadata.schema = schemas
    filemetadata.num_rows = nrows
    filemetadata.created_by = "Parquet.jl $(Parquet.PARQUET_JL_VERSION)"

    # create row_groups
    # TODO do multiple row_groups
    row_group = PAR2.RowGroup()

    row_group.columns = col_chunk_metas
    row_group.total_byte_size = Int64(sum(x->x.meta_data.total_compressed_size, col_chunk_metas))
    row_group.num_rows = nrows
    if ismissing(row_group_file_offset)
        error("row_group_file_offset is not set")
    else
        row_group.file_offset = row_group_file_offset
    end
    row_group.total_compressed_size = Int64(sum(x->x.meta_data.total_compressed_size, col_chunk_metas))

    filemetadata.row_groups = [row_group]

    filemetadata_size = write_thrift(io, filemetadata)

    write(io, UInt32(filemetadata_size) |> htol)
    write(io, "PAR1")
end


text:
#include("../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

@testset "ttv_cartesian" begin

# This routine takes derivative of transit times with respect
# to the initial cartesian coordinates of bodies. [x]
#n = 8
n = 3
H = [3,1,1]
t0 = 7257.93115525-7300.0
#h  = 0.12
h  = 0.04
#h  = 0.02
#tmax = 600.0
#tmax = 1000.0
#tmax = 100.0
tmax = 10.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')
elements[:,3] .-= 7300.0

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
#println("ntt: ",ntt)
tt  = zeros(n,maximum(ntt))
tt1 = zeros(n,maximum(ntt))
tt2 = zeros(n,maximum(ntt))
tt3 = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12
#dq = ttv_elements!(H,t0,h,tmax,elements,tt1,count1,0.0,0,0,rstar)
# Now call with half the timestep:
count2 = zeros(Int64,n)
count3 = zeros(Int64,n)
dq = ttv_elements!(H,t0,h/2,tmax,elements,tt2,count2,0.0,0,0,rstar)

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtdq0 = zeros(n,maximum(ntt),7,n)
dtdelements = ttv_elements!(H,t0,h,tmax,elements,tt,count,dtdq0,rstar)
#read(STDIN,Char)

# Check that this is working properly:
for i=1:n
  for j=1:count2[i]
#    println(i," ",j," ",tt[i,j]," ",tt2[i,j]," ",tt[i,j]-tt2[i,j]," ",tt1[i,j]-tt2[i,j])
  end
end
#read(STDIN,Char)

# Compute derivatives numerically:
# Compute the numerical derivative:
dtdq0_num = zeros(BigFloat,n,maximum(ntt),7,n)
itdq0 = zeros(Int64,n,maximum(ntt),7,n)
dlnq = big(1e-15)
hbig = big(h); t0big = big(t0); tmaxbig=big(tmax); tt2big = big.(tt2); tt3big = big.(tt3)
for jq=1:n
  for iq=1:7
    elements2  = big.(elements)
    dq_plus = ttv_elements!(H,t0big,hbig,tmaxbig,elements2,tt2big,count2,dlnq,iq,jq,big(rstar))
    elements3  = big.(elements)
    dq_minus = ttv_elements!(H,t0big,hbig,tmaxbig,elements3,tt3big,count3,-dlnq,iq,jq,big(rstar))
    for i=1:n
      for k=1:count2[i]
        # Compute double-sided derivative for more accuracy:
        dtdq0_num[i,k,iq,jq] = (tt2big[i,k]-tt3big[i,k])/(dq_plus-dq_minus)
#        println(i," ",k," ",iq," ",jq," ",tt2big[i,k]," ",tt3big[i,k]," ")
      end
    end
  end
end

nbad = 0
ntot = 0
diff_dtdq0 = zeros(n,maximum(ntt),7,n)
mask = zeros(Bool, size(dtdq0))
for i=2:n, j=1:count[i], k=1:7, l=1:n
  if abs(dtdq0[i,j,k,l]-dtdq0_num[i,j,k,l]) > 0.1*abs(dtdq0[i,j,k,l]) && ~(abs(dtdq0[i,j,k,l]) == 0.0  && abs(dtdq0_num[i,j,k,l]) < 1e-3)
#    println(i," ",j," ",k," ",l," ",dtdq0[i,j,k,l]," ",dtdq0_num[i,j,k,l]," ",itdq0[i,j,k,l])
    nbad +=1
  end
  diff_dtdq0[i,j,k,l] = abs(dtdq0[i,j,k,l]-dtdq0_num[i,j,k,l])
  if k != 2 && k != 5
    mask[i,j,k,l] = true
  end
  ntot +=1
end

tt_big = big.(tt); elementsbig = big.(elements); rstarbig = big(rstar)
dqbig = ttv_elements!(H,t0big,hbig,tmaxbig,elementsbig,tt_big,count,big(0.0),0,0,rstarbig)
# Now halve the time steps:
tt_big_half = copy(tt_big)
dqbig = ttv_elements!(H,t0big,hbig/2,tmaxbig,elementsbig,tt_big_half,count1,big(0.0),0,0,rstarbig)

# Compute the derivatives in BigFloat precision to see whether finite difference
# derivatives or Float64 derivatives are imprecise at the start:
dtdq0_big = zeros(BigFloat,n,maximum(ntt),7,n)
hbig = big(h); tt_big = big.(tt); elementsbig = big.(elements); rstarbig = big(rstar)
dtdelements_big = ttv_elements!(H,t0big,hbig,tmaxbig,elementsbig,tt_big,count,dtdq0_big,rstarbig)

#=
using PyPlot

clf()
# Plot the difference in the TTVs:
for i=2:3
#  diff1 = abs.(tt1[i,2:count1[i]]./tt_big[i,2:count1[i]]-1.0);
  diff1 = convert(Array{Float64,1},abs.(tt1[i,2:count1[i]].-tt_big[i,2:count1[i]])/elements[i,2]);
  dtt=tt[i,2:count1[i]].-tt[i,1]
  loglog(dtt,diff1);
#  diff2 = abs.(tt2[i,2:count1[i]]./tt_big_half[i,2:count1[i]]-1.0);
#  diff2 = abs.(tt2[i,2:count1[i]].-tt_big_half[i,2:count1[i]])/elements[i,2];
#  loglog(tt[i,2:count[i]]-tt[i,1],diff2);
end
loglog([1.0,1024.0],2e-15*[1,2^15],":")
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdelements_big[i,2:count[i],k,l])) > 0
    diff1 = convert(Array{Float64,1},abs.(dtdelements_big[i,2:count[i],k,l]./dtdelements[i,2:count[i],k,l].-1));
    diff3 = convert(Array{Float64,1},abs.(asinh.(dtdelements_big[i,2:count[i],k,l])-asinh.(dtdelements[i,2:count[i],k,l])));
#    loglog(tt[i,2:count[i]]-tt[i,1],diff1);
    dtt = tt[i,2:count[i]].-tt[i,1]
    loglog(dtt,diff3);
    println(i," ",k," ",l," frac error: ",convert(Float64,maximum(diff1))," asinh error: ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
  if maximum(abs.(dtdq0_big[i,2:count[i],k,l])) > 0
    diff1 = convert(Array{Float64,1},abs.(dtdq0[i,2:count[i],k,l]./dtdq0_big[i,2:count[i],k,l].-1.0));
#    diff2 = abs.(asinh.(dtdq0_big[i,2:count[i],k,l]).-asinh.(dtdq0_num[i,2:count[i],k,l]));
    diff3 = convert(Array{Float64,1},abs.(asinh.(dtdq0_big[i,2:count[i],k,l]).-asinh.(dtdq0[i,2:count[i],k,l])));
#    loglog(tt[i,2:count[i]].-tt[i,1],diff1);
    dtt = tt[i,2:count[i]].-tt[i,1]
    loglog(dtt,diff3,linestyle=":");
#    loglog(tt[i,2:count[i]].-tt[i,1],diff2,".");
    println(i," ",k," ",l," frac error: ",convert(Float64,maximum(diff1))," asinh error: ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
mederror = zeros(size(tt))
for i=2:3
  for j=1:count1[i]
    data_list = Float64[]
    for k=1:7, l=1:3
      if abs(dtdq0_num[i,j,k,l]) > 0
        push!(data_list,abs(dtdq0[i,j,k,l]/dtdq0_num[i,j,k,l]-1.0))
      end
    end
    mederror[i,j] = median(data_list)
  end
end

# Plot a line that scales as time^{3/2}:

loglog([1.0,1024.0],1e-12*[1,2^15],":",linewidth=3)

=#
#println("Max diff log(dtdq0): ",maximum(abs.(dtdq0_num[mask]./dtdq0[mask].-1.0)))
#println("Max diff asinh(dtdq0): ",maximum(abs.(asinh.(dtdq0_num[mask]).-asinh.(dtdq0[mask]))))
#unit = ones(dtdq0[mask])
#@test isapprox(dtdq0[mask]./convert(Array{Float64,4},dtdq0_num)[mask],unit;norm=maxabs)
#@test isapprox(dtdq0[mask],convert(Array{Float64,4},dtdq0_num)[mask];norm=maxabs)
@test isapprox(asinh.(dtdq0[mask]),asinh.(convert(Array{Float64,4},dtdq0_num)[mask]);norm=maxabs)
#end

#
#nderiv = n^2*7*maximum(ntt)
#loglog(abs.(reshape(dtdq0,nderiv)),abs.(reshape(convert(Array{Float64,4},dtdq0_num),nderiv)),".")
#loglog(abs.(reshape(dtdq0,nderiv)),abs.(reshape(convert(Array{Float64,4},diff_dtdq0),nderiv)),".")

## Make a plot of some TTVs:
#
#fig,axes = subplots(4,2)
#
#for i=2:8
#  ax = axes[i-1]
#  fn = zeros(Float64,2,count1[i])
#  sig = ones(count1[i])
#  tti1 = tt1[i,1:count1[i]]
#  tti2 = tt2[i,1:count2[i]]
#  for j=1:count1[i]
#    fn[1,j] = 1.0
#    fn[2,j] = round(Int64,(tti1[j]-elements[i,3])/elements[i,2])
#  end
#  coeff,cov = regress(fn,tti1,sig)
#  tt_ref1 = coeff[1]+coeff[2]*fn[2,:]
#  ttv1 = (tti1-tt_ref1)*24.*60.
#  coeff,cov = regress(fn,tti2,sig)
#  tt_ref2 = coeff[1]+coeff[2]*fn[2,:]
#  ttv2 = (tti2-tt_ref2)*24.*60.
#  ax[:plot](tti1,ttv1)
##  ax[:plot](tti2,ttv2)
#  ax[:plot](tti2,((ttv1-ttv2)-mean(ttv1-ttv2)))
#  println(i," ",coeff," ",elements[i,2:3]," ",coeff[1]-elements[i,3]," ",coeff[2]-elements[i,2])
#  println(i," ",maximum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec ", minimum(ttv1-ttv2-mean(ttv1-ttv2))*60.," sec" )
end

ext: jl


text:
"""
Created on Sun Apr 26 2020

@author: Yoann Pradat

    CentraleSupelec
    MICS laboratory
    9 rue Juliot Curie, Gif-Sur-Yvette, 91190 France

Rank selection method in NMF based on multiple metrics computed from the consensus matrices. These metrics include
    - the cophenetic correlation
        J-P. Brunet et al, Metagenes and molecular pattern discovery using matrix factorization. 2009, PNAS.
    - the dispersion
        H. Kim and H. Park, Sparse non-negative matrix factorizations via alternating non-negativity-constrained least
        squares for microarray data analysis. 2007, Bioinformatics
    - the silhouette width
        Rousseeuw PJ. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Jour- nal of Computational and Applied Mathematics. 1987
"""

abstract type AbstractRCParams <: AbstractNMFParams end

"""
    RCParams([,K_min, K_max, n_iter, pert, seed])

Parameters passed to `rank_by_consensus`.

# Arguments
- `K_min::Integer`     : Minimum rank evaluated.
- `K_max::Integer`     : Maximum rank evaluated.
- `n_iter::Integer`    : Total number of iterations, each iteration corresponding to one run of NMF on one random
    initialization of the matrices W and H.
- `seed::Integer`      : Seed for the random number generator.
"""

mutable struct RCParams <: AbstractRCParams
    K_min::Integer
    K_max::Integer
    n_iter::Integer
    seed::Integer

    function RCParams(;
        K_min  = 2,
        K_max  = 10,
        n_iter = 100,
        seed   = 0)

        K_min > 1 || throw(ArgumentError("K_min must be superior or equal to 2"))
        K_max >= K_min || throw(ArgumentError("K_max must be superior or equal to K_min"))
        n_iter > 0 || throw(ArgumentError("n_iter must be superior or equal to 1"))

        new(K_min, K_max, n_iter, seed)
    end
end

abstract type AbstractRCNMFResults <: AbstractNMFResults end

"""
    RCNMFResults(grid_nmf_results, list_rank, list_coph, list_idxs, list_cons)

Struct for storing results returned by `rank_by_consensus`.
"""
mutable struct RCNMFResults <: AbstractRCNMFResults
    df_metrics::DataFrame

    list_idxs::Vector{Vector{Int64}}
    list_cons::Vector{Array{Float64, 2}}

    function RCNMFResults(
        df_metrics::DataFrame,
        list_idxs::Vector{Vector{Int64}},
        list_cons::Vector{Array{Float64, 2}})

        new(df_metrics, list_idxs, list_cons)
    end
end

"""
    _rank_by_consensus_one((V, nmf))

Helper function for `rank_by_consensus`.
"""
function _rank_by_consensus_one((V, nmf))
    nmf_results = nmf.solver(V, nmf.global_params, nmf.local_params)
end

"""
    rank_by_consensus(V, rc_params, nmf)

To run in parallel the arguments must have been made available to all workers using `@everywhere` macro.
Returns in a tuple the ranks, cophenetics correlations, ordering indices and consensus matrices.
"""

function rank_by_consensus(V::Matrix{T}, rc_params::RCParams, nmf::NMF) where T <: Real
    #### fix random number generator
    rng = MersenneTwister(rc_params.seed)

    #### build grid
    list_nmf    = Vector{NMF}()

    for rank in rc_params.K_min:rc_params.K_max
        for seed in randperm(rng, rc_params.n_iter)
            nmf_rank_seed = copy(nmf)
            nmf_rank_seed.global_params.rank = rank
            nmf_rank_seed.global_params.seed = seed

            push!(list_nmf, nmf_rank_seed)
        end
    end

    grid = vec(collect(Iterators.product([V], list_nmf)))

    #### run in parallel
    grid_res = @showprogress pmap(grid) do g
        _rank_by_consensus_one(g)
    end

    grid_rank     = (x -> x.global_params.rank).(grid_res)
    grid_W        = (x -> x.W).(grid_res)
    grid_H        = (x -> x.H).(grid_res)
    list_rank     = Vector{Int64}()
    list_coph     = Vector{Float64}()
    list_disp     = Vector{Float64}()
    list_silh     = Vector{Float64}()
    list_idxs     = Vector{Vector{Int64}}()
    list_cons     = Vector{Array{Float64, 2}}()
    list_sparse_W = Vector{Float64}()
    list_sparse_H = Vector{Float64}()

    for rank in rc_params.K_min:rc_params.K_max
        #### compute consensus matrix
        cons = mean((H -> _connectivity_matrix(H)).(grid_H[grid_rank .== rank]))

        #### perform hierarchical clustering
        hc = hclust(1 .- cons, linkage=:average, branchorder=:r)

        #### compute cophenetic correlation
        coph = cons_cophenetic(cons, hc)

        #### compute dispersion
        disp = cons_dispersion(cons)

        #### compute avg silhouette width
        silh = cons_avg_silhouette(cons, hc, rank)

        #### reorder consensus matrix
        cons = cons[hc.order, hc.order]

        ### compute average sparseness of rows of W
        sparse_W = mean((W -> matr_avg_row_sparseness(W)).(grid_W[grid_rank .== rank]))

        ### compute average sparseness of columns of H
        sparse_H = mean((H -> matr_avg_col_sparseness(H)).(grid_H[grid_rank .== rank]))

        push!(list_rank, rank)
        push!(list_coph, coph)
        push!(list_disp, disp)
        push!(list_silh, silh)
        push!(list_idxs, hc.order)
        push!(list_cons, cons)
        push!(list_sparse_W, sparse_W)
        push!(list_sparse_H, sparse_H)
    end

    df_metrics = DataFrame(
        rank       = list_rank,
        cophenetic = list_coph,
        dispersion = list_disp,
        silhouette = list_silh,
        sparse_W   = list_sparse_W,
        sparse_H   = list_sparse_H
       )

    RCNMFResults(
        df_metrics,
        list_idxs,
        list_cons
    )
end


text:
#!/usr/bin/julia
using ArgParse

function empty(input::String)
    println(input)
end

# selects a parser based on extension, if there is a match, then that parser is
# handed the file's raw text
function selectParser(file_names)

    # get all the file extensions and select a parser based on them
    for name in file_names
        ext = getExtension(name)
        parse_function = empty
        if ext == "cpp"

        # place holder for later parsers
        elseif false

        else
            println("parser not found for $name, skipping $name .")

        end

        # open the file and hand it off to the parser
        if parse_function != empty

            try
                src_file = open(name,"r")
                contents = read(src_file,String)
                close(src_file)
                parse_function(contents)
            catch e

                println(sprint(showerror, e, catch_backtrace()))

            finally
            end

        end
    end
end

function getExtension(file_name::String)::String

    # get the last . and then return if the extension matches
    dot_index = findlast(isequal('.'),file_name)
    if typeof(dot_index) == Nothing
        return ""
    end

    return SubString(file_name,dot_index+1,length(file_name))
end

# iterate through the array from right to left, returning all text before the . of the file extension
function extensionMatches(filename::String,extenstion::String)::Bool

    file_ext = getExtension(filename)

    return file_ext == extenstion
end

function main(args)

    settings = ArgParseSettings("test text")

    @add_arg_table settings begin
        "-r"
        help = "Step through all directories, scanning all files in them"
        action = :store_true

        "-o"
        help = "Specify the output file name like this -o outputFile.svg"
        arg_type = String
        default = "controlFlow.svg"

        "filenames"
        nargs = '*'
        help = " a list of files to parse"
        required = true

        "--ext"
        nargs = 1
        help = "file extensions to search, such as cpp, c, py, pl, jl, js, etc. Example usage: \"--ext cpp \". Only accepts one argument"
        arg_type = String
    end

    parsed_args = parse_args(args, settings)

    # recursive mode
    if parsed_args["r"]

        if parsed_args["ext"][1] == ""
            println("The --ext argument is required for recursive mode.")
            exit()
        end

        file_list = []
        for file_name in parsed_args["filenames"]
            if extensionMatches(file_name,parsed_args["ext"])
                push!(file_list,file_name)
            end
        end
        if length(file_list) == 0
            println("No files with the extension $(parsed_args["ext"]) were found")
        end

    # not recursive mode
    else
        selectParser(parsed_args["filenames"])
    end
end

# run the program
main(ARGS)

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
module Htop

import Htop_jll, Ncurses_jll
export htop
const htop = addenv(Htop_jll.htop(), "TERMINFO_DIRS" => joinpath(Ncurses_jll.artifact_dir, "share", "terminfo"))
Htop_jll.is_available() ||
    (__init__() = @warn("Htop is not available on this system ($((Base.BinaryPlatforms.HostPlatform())))."))

end # module


text:
a = collect(1:10)
b = Array{Float64}(undef,10)
b .= NaN
Plots.plot(a,b)


text:
using Romberg
using Test
using Trapz

@testset "Romberg.jl" begin
    # Test interfaces
    x = range(0, π, length=2^8+1)
    y = sin.(x)

    @test romberg(x, y) == romberg(x, y, 8)

    max_steps = 7
    R = zeros(max_steps+1, max_steps+1)
    @test romberg!(R, x, y)[end] == romberg(x, y, max_steps)

    # Test ispow2(length(x) - 1) == false
    x = range(0, π, length=2^8)
    y = sin.(x)
    @test romberg(x, y) ≈ 2

    x = range(0, π, length=2^8-1)
    y = sin.(x)
    @test romberg(x, y) ≈ 2

    # Test `max_steps` too large
    @test_throws DomainError romberg(x, y, 9)

    # Test length(x) == 1
    x = 0:0
    y = sin.(x)
    R = zeros(1, 1)
    @test romberg(x, y) == 0
    @test romberg!(R, x, y) == zeros(1,1)

    # Test length(x) == 2
    x = range(0, 1, length=2)
    y = x.^2
    @test romberg(x,y) ≈ trapz(x, y)  # integration is inaccurate

    # Test length(x) == 3
    x = range(0, 1, length=3)
    y = x.^2
    @test romberg(x, y) ≈ 1/3

    # Integrate different functions
    x = range(0, π, length=2^8+1)
    y = sin.(x)
    @test romberg(x, y) ≈ 2

    x = range(0, 1, length=2^8+1)
    y = x.^3
    @test romberg(x, y) ≈ 1/4

    x = range(0, π/2, length=2^5+1)
    y = sin.(x).^2
    @test romberg(x, y) ≈ π/4

    m = 3
    n = 4
    x = range(0, π, length=2^8+1)
    y = sin.(m*x).*cos.(n*x)
    @test romberg(x, y) ≈ 2*m/(m^2 - n^2)

    # this one requires lots of samples...
    a = 15
    x = range(0, a, length=2^16+1)
    y = sqrt.(a^2 .- x.^2)
    @test romberg(x, y) ≈ π*a^2/4

    # tricky to integrate
    x = range(1e-15, 1, length=2^16+1)
    y = log.(x)./(1 .+ x)
    @test romberg(x, y) ≈ -π^2/12 atol=1e-4
end


text:
# ---
# title: 1438. Longest Continuous Subarray With Absolute Diff Less Than or Equal to Limit
# id: problem1438
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sliding Window
# link: <https://leetcode.com/problems/longest-continuous-subarray-with-absolute-diff-less-than-or-equal-to-limit/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `limit`, return the size of
# the longest **non-empty** subarray such that the absolute difference between
# any two elements of this subarray is less than or equal to `limit` _._
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [8,2,4,7], limit = 4
#     Output: 2
#     Explanation: All subarrays are:
#     [8] with maximum absolute diff |8-8| = 0 <= 4.
#     [8,2] with maximum absolute diff |8-2| = 6 > 4.
#     [8,2,4] with maximum absolute diff |8-2| = 6 > 4.
#     [8,2,4,7] with maximum absolute diff |8-2| = 6 > 4.
#     [2] with maximum absolute diff |2-2| = 0 <= 4.
#     [2,4] with maximum absolute diff |2-4| = 2 <= 4.
#     [2,4,7] with maximum absolute diff |2-7| = 5 > 4.
#     [4] with maximum absolute diff |4-4| = 0 <= 4.
#     [4,7] with maximum absolute diff |4-7| = 3 <= 4.
#     [7] with maximum absolute diff |7-7| = 0 <= 4.
#     Therefore, the size of the longest subarray is 2.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [10,1,2,4,7,2], limit = 5
#     Output: 4
#     Explanation: The subarray [2,4,7,2] is the longest since the maximum absolute diff is |2-7| = 5 <= 5.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [4,2,2,2,4,4,2,2], limit = 0
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10^5`
#   * `1 <= nums[i] <= 10^9`
#   * `0 <= limit <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""

This function fills a pre-allocated square matrix `D` with the taxonomic
distance between taxa in a vector `tax`. The keywords areguments are `d` (the
distance as a function to closest matching level), `strict` (a boolean to decide
whether the distance of a taxon with itself should be read from the distance
dictionary, or set to 0), and other keywords passed to the `lineage` function.

Because the distances are symetrical, there are only (n(n-1))/2 measurements to
do.

By default, this function uses the distances from Shimatani (2001):

| equal rank  | distance |
|-------------|----------|
| species     | 0        |
| genus       | 1        |
| family      | 2        |
| subclass    | 3        |
| `:fallback` | 4        |

Shimatani, K. 2001. On the measurement of species diversity incorporating
species differences. Oikos 93:135–147.
"""
function taxonomicdistance!(D::Matrix{Float64}, tax::Vector{NCBITaxon}; d::Dict{Symbol,Float64}=Dict(:species=>0.0, :genus=>1.0, :family=>2.0, :subclass=>3.0, :fallback=>4.0), strict::Bool=true, kwargs...)
    if size(D,1) != size(D,2)
        throw(DimensionMismatch("The matrix `D` ($(size(D))) must be square" ))
    end
    if size(D,1) != length(tax)
        throw(DimensionMismatch("The matrix `D` ($(size(D,1))) must have the same sizes as the length of `tax` ($(length(tax)))."))
    end
    ancestry = lineage.(tax; kwargs...)
    for i in 1:(length(tax)-1)
        D[i,i] = strict ? get(d, rank(tax[i]), d[:fallback]) : 0.0
        ai = ancestry[i]
        for j in (i+1):length(tax)
            aj = ancestry[j]
            shared = rank.(ai ∩ aj)
            D[i,j] = D[j,i] = minimum([get(d, s, d[:fallback]) for s in shared])
        end
    end
    return D
end

"""
    taxonomicdistance(tax::Vector{NCBITaxon}; kwargs...)

See `taxonomicdistance!`.
"""
function taxonomicdistance(tax::Vector{NCBITaxon}; kwargs...)
    D = zeros(Float64, length(tax), length(tax))
    taxonomicdistance!(D, tax; kwargs...)
    return D
end


text:
module FV3Kernels

export State

include("states.jl")

#=
export SillyType

struct SillyType{T}
    a :: T
    sin_a :: T
end

function SillyType(a)
    sin_a = sin(a)
    return SillyType(a, sin_a)
end
=#

end # module


text:
using SafeTestsets
@safetestset "Sniffles test" begin include("testsets.jl") end


text:
using Test

include("robot-simulator.jl")

@testset "constructor" begin
    r = Robot((0, 0), NORTH)
    @test position(r) == Point(0, 0)
    @test heading(r) == NORTH

    r = Robot((-1, -1), SOUTH)
    @test position(r) == Point{Int}(-1, -1)
    @test heading(r) == SOUTH
end

@testset "mutating functions should return robot" begin
    r = Robot((0, 0), NORTH)
    @test r == turn_right!(r)
    @test r == turn_left!(r)
    @test r == advance!(r)
    @test r == move!(r, "A")
end

@testset "rotate +π/2" begin
    r = Robot((0, 0), NORTH)
    turn_right!(r)
    @test position(r) == Point(0, 0)
    @test heading(r) == EAST
    turn_right!(r)
    @test heading(r) == SOUTH
    turn_right!(r)
    @test heading(r) == WEST
    turn_right!(r)
    @test heading(r) == NORTH
end

@testset "rotate -π/2" begin
    r = Robot((0, 0), NORTH)
    turn_left!(r)
    @test position(r) == Point(0, 0)
    @test heading(r) == WEST
    turn_left!(r)
    @test heading(r) == SOUTH
    turn_left!(r)
    @test heading(r) == EAST
    turn_left!(r)
    @test heading(r) == NORTH
end

@testset "advance" begin
    r = Robot((0, 0), NORTH)
    advance!(r)
    @test heading(r) == NORTH
    @test position(r) == Point(0, 1)

    @test position(advance!(Robot((0, 0), SOUTH))) == Point(0, -1)
    @test position(advance!(Robot((0, 0), EAST))) == Point(1, 0)
    @test position(advance!(Robot((0, 0), WEST))) == Point(-1, 0)
end

@testset "instructions" begin
    @testset "moving east and north from README" begin
        r = Robot((7, 3), NORTH)
        move!(r, "RAALAL")
        @test position(r) == Point(9, 4)
        @test heading(r) == WEST
    end

    @testset "move west and north" begin
        r = Robot((0, 0), NORTH)
        move!(r, "LAAARALA")
        @test position(r) == Point(-4, 1)
        @test heading(r) == WEST
    end

    @testset "move west and south" begin
        r = Robot((2, -7), EAST)
        move!(r, "RRAAAAALA")
        @test position(r) == Point(-3, -8)
        @test heading(r) == SOUTH
    end

    @testset "move east and north" begin
        r = Robot((8, 4), SOUTH)
        move!(r, "LAAARRRALLLL")
        @test position(r) == Point(11, 5)
        @test heading(r) == NORTH
    end
end


text:
export RedPitayaCluster, master, readDataPeriods, numChan, readDataSlow, readFrames, readPeriods, readPipelinedSamples, startPipelinedData, collectSamples!, readData, readDataPeriods, convertSamplesToFrames!, convertSamplesToFrames, SampleChunk

import Base: length, iterate, getindex, firstindex, lastindex

struct RedPitayaCluster
  rp::Vector{RedPitaya}
end

#TODO: set first RP to master
function RedPitayaCluster(hosts::Vector{String}, port=5025)
  # the first RP is the master
  rps = RedPitaya[ RedPitaya(host, port, i==1) for (i,host) in enumerate(hosts) ]

  @sync for rp in rps
    @async triggerMode(rp, "EXTERNAL")
  end

  return RedPitayaCluster(rps)
end

length(rpc::RedPitayaCluster) = length(rpc.rp)
numChan(rpc::RedPitayaCluster) = 2*length(rpc)
master(rpc::RedPitayaCluster) = rpc.rp[1]

# Indexing Interface
function getindex(rpc::RedPitayaCluster, index::Integer)
  1 <= index <= length(rpc) || throw(BoundsError(rpc.rp, index))
  return rpc.rp[index]
end
firstindex(rpc::RedPitayaCluster) = start_(rpc)
lastindex(rpc::RedPitayaCluster) = length(rpc)

# Iterable Interface
start_(rpc::RedPitayaCluster) = 1
next_(rpc::RedPitayaCluster,state) = (rpc[state],state+1)
done_(rpc::RedPitayaCluster,state) = state > length(rpc)
iterate(rpc::RedPitayaCluster, s=start_(rpc)) = done_(rpc, s) ? nothing : next_(rpc, s)

function RPInfo(rpc::RedPitayaCluster)
  return RPInfo([RPPerformance([]) for i = 1:length(rpc)])
end

const _currentFrame = Ref(0)

function currentFrame(rpc::RedPitayaCluster)
  _currentFrame[] = currentFrame(rpc[1]) #[ currentFrame(rp) for rp in rpc.rp ]
  @debug "Current frame: $(_currentFrame[])"
  #return minimum(currentFrames)
  return _currentFrame[]
end

function currentPeriod(rpc::RedPitayaCluster)
  currentPeriods = currentPeriod(rpc[1])  #[ currentPeriod(rp) for rp in rpc.rp ]
  @debug "Current period: $currentPeriods"
  #return minimum(currentPeriods)
  return currentPeriods
end

function currentWP(rpc::RedPitayaCluster)
  return currentWP(master(rpc))
end

for op in [:periodsPerFrame, :samplesPerPeriod, :decimation, :keepAliveReset, :sequenceRepetitions,
           :triggerMode, :samplesPerSlowDACStep, :slowDACStepsPerFrame,
           :rampUpTime, :rampUpFraction]
  @eval $op(rpc::RedPitayaCluster) = $op(master(rpc))
  @eval begin
    function $op(rpc::RedPitayaCluster, value)
      @sync for rp in rpc
        @async $op(rp, value)
      end
    end
  end
end

for op in [:connectADC, :stopADC, :disconnect, :connect, :clearSequence]
  @eval begin
    function $op(rpc::RedPitayaCluster)
      @sync for rp in rpc
        @async $op(rp)
      end
    end
  end
end

function startADC(rpc::RedPitayaCluster)
  @sync for rp in rpc
    @async startADC(rp)
  end
end

function masterTrigger(rpc::RedPitayaCluster, val::Bool)
    if val
        masterTrigger(master(rpc), val)
    else
        keepAliveReset(rpc, true)
        masterTrigger(master(rpc), false)
        ramWriterMode(rpc, "TRIGGERED")
        keepAliveReset(rpc, false)
    end
end
masterTrigger(rpc::RedPitayaCluster) = masterTrigger(master(rpc))
bufferSize(rpc::RedPitayaCluster) = bufferSize(master(rpc))

# "TRIGGERED" or "CONTINUOUS"
function ramWriterMode(rpc::RedPitayaCluster, mode::String)
  @sync for rp in rpc
    @async ramWriterMode(rp, mode)
  end
end

for op in [:amplitudeDAC, :amplitudeDACNext, :frequencyDAC, :phaseDAC, :modulusFactorDAC, :modulusDAC]
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, component::Integer)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc[idxRP], chanRP, component)
  end
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, component::Integer, value)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc[idxRP], chanRP, component, value)
  end
end

for op in [:signalTypeDAC,  :DCSignDAC, :jumpSharpnessDAC]
  @eval function $op(rpc::RedPitayaCluster, chan::Integer)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc[idxRP], chanRP)
  end
  @eval function $op(rpc::RedPitayaCluster, chan::Integer, value)
    idxRP = div(chan-1, 2) + 1
    chanRP = mod1(chan, 2)
    return $op(rpc[idxRP], chanRP, value)
  end
end

function setSlowDAC(rpc::RedPitayaCluster, chan, value)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod1(chan, 2)
  setSlowDAC(rpc[idxRP], chanRP, value)
end

function getSlowADC(rpc::RedPitayaCluster, chan::Integer)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod1(chan, 2)
  getSlowADC(rpc[idxRP], chanRP)
end

function offsetDAC(rpc::RedPitayaCluster, chan, value)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod1(chan, 2)
  offsetDAC(rpc[idxRP], chanRP, value)
end

function offsetDAC(rpc::RedPitayaCluster, chan::Integer)
  idxRP = div(chan-1, 2) + 1
  chanRP = mod1(chan, 2)
  offsetDAC(rpc[idxRP], chanRP)
end

function numSlowADCChan(rpc::RedPitayaCluster)
  tmp = [0 for rp in rpc]
  @sync for (d, rp) in enumerate(rpc)
    @async tmp[d] = numSlowADCChan(rp)
  end
  return sum(tmp)
end

function numSlowADCChan(rpc::RedPitayaCluster, num)
  @sync for rp in rpc
    @async numSlowADCChan(rp, num)
  end
  return
end

function passPDMToFastDAC(rpc::RedPitayaCluster, val::Vector{Bool})
  @sync for (d,rp) in enumerate(rpc)
    @async passPDMToFastDAC(rp, val[d])
  end
end

function passPDMToFastDAC(rpc::RedPitayaCluster)
  result = [false for rp in rpc]
  @sync for (d, rp) in enumerate(rpc)
    @async result[d] = passPDMToFastDAC(rp)
  end
  return result
end

function appendSequence(rpc::RedPitayaCluster, seq::AbstractSequence)
  @sync for rp in rpc
    @async appendSequence(rp, seq)
  end
end

function appendSequence(rpc::RedPitayaCluster, index, seq::AbstractSequence)
  appendSequence(rpc[index], seq)
end

function prepareSequence(rpc::RedPitayaCluster)
  success = [false for rp in rpc]
  @sync for (i, rp) in enumerate(rpc)
    @async success[i] = prepareSequence(rp)
  end
  all(success)
end

computeRamping(rpc::RedPitayaCluster, stepsPerSeq, time, fraction) = computeRamping(master(rpc), stepsPerSeq, time, fraction)

modeDAC(rpc::RedPitayaCluster) = modeDAC(master(rpc))

function modeDAC(rpc::RedPitayaCluster, mode::String)
  @sync for rp in rpc
    @async modeDAC(rp, mode)
  end
end

function slowDACInterpolation(rpc::RedPitayaCluster, enable::Bool)
  @sync for rp in rpc
    @async slowDACInterpolation(rp, enable)
  end
end

include("ClusterView.jl")

struct SampleChunk
  samples::Matrix{Int16}
  performance::Vector{PerformanceData}
end

function startPipelinedData(rpc::RedPitayaCluster, reqWP::Int64, numSamples::Int64, chunkSize::Int64)
  @sync for rp in rpc
    @async startPipelinedData(rp, reqWP, numSamples, chunkSize)
  end
end

function readSamples(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, wpStart::Int64, numOfRequestedSamples::Int64; chunkSize::Int64 = 25000, rpInfo=nothing)
  numOfReceivedSamples = 0
  index = 1
  rawData = zeros(Int16, numChan(rpu), numOfRequestedSamples)
  chunkBuffer = zeros(Int16, chunkSize * 2, length(rpu))

  while numOfReceivedSamples < numOfRequestedSamples
    wpRead = wpStart + numOfReceivedSamples
    wpWrite = currentWP(rpu)
    chunk = min(numOfRequestedSamples - numOfReceivedSamples, chunkSize)

    # Wait for data to be written
    while wpRead + chunk >= wpWrite
      wpWrite = currentWP(rpu)
      @debug wpWrite
    end

    # Collect data
    collectSamples!(rpu, wpRead, chunk, rawData, chunkBuffer, index, rpInfo=rpInfo)
    index += chunk
    numOfReceivedSamples += chunk
  end
  return rawData
end

function readPipelinedSamples(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, wpStart::Int64, numOfRequestedSamples::Int64; chunkSize::Int64 = 25000, rpInfo=nothing)
  numOfReceivedSamples = 0
  index = 1
  rawData = zeros(Int16, numChan(rpu), numOfRequestedSamples)
  chunkBuffer = zeros(Int16, chunkSize * 2, length(rpu))

  startPipelinedData(rpu, wpStart, numOfRequestedSamples, chunkSize)
  while numOfReceivedSamples < numOfRequestedSamples
    wpRead = wpStart + numOfReceivedSamples
    chunk = min(numOfRequestedSamples - numOfReceivedSamples, chunkSize)
    collectSamples!(rpu, wpRead, chunk, rawData, chunkBuffer, index, rpInfo=rpInfo)
    index += chunk
    numOfReceivedSamples += chunk
  end
  return rawData

end

function collectSamples!(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, wpRead::Int64, chunk::Int64, rawData::Matrix{Int16}, chunkBuffer, index; rpInfo=nothing)
  done = zeros(Bool, length(rpu))
  iterationDone = Condition()
  timeoutHappened = false

  for (d, rp) in enumerate(rpu)
    @async begin
      buffer = view(chunkBuffer, 1:(2 * chunk), d)
      (u, perf) = readSamplesChunk_(rp, Int64(wpRead), Int64(chunk), buffer)
      samples = reshape(u, 2, chunk)
      rawData[2*d-1, index:(index + chunk - 1)] = samples[1, :]
      rawData[2*d, index:(index + chunk - 1)] = samples[2, :]

      # Status
      if perf.status.overwritten
        @error "RP $d: Requested data from $wpRead until $(wpRead+chunk) was overwritten"
      end
      if perf.status.corrupted
        @error "RP $d: Requested data from $wpRead until $(wpRead+chunk) might have been corrupted"
      end

      # Performance
      if !isnothing(rpInfo)
        push!(rpInfo.performances[d].data, perf)
      end

      done[d] = true
      if (all(done))
        notify(iterationDone)
      end
    end
  end

  # Setup Timeout
  t = Timer(_timeout)
  @async begin
    wait(t)
    notify(iterationDone)
    timeoutHappened = true
  end

  wait(iterationDone)
  close(t)
  if timeoutHappened
    error("Timout reached when reading from sockets")
  end

end

function readPipelinedSamples(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, wpStart::Int64, numOfRequestedSamples::Int64, channel::Channel; chunkSize::Int64 = 25000)
  numOfReceivedSamples = 0
  chunkBuffer = zeros(Int16, chunkSize * 2, length(rpu))

  startPipelinedData(rpu, wpStart, numOfRequestedSamples, chunkSize)
  while numOfReceivedSamples < numOfRequestedSamples
    wpRead = wpStart + numOfReceivedSamples
    chunk = min(numOfRequestedSamples - numOfReceivedSamples, chunkSize)
    collectSamples!(rpu, wpRead, chunk, channel, chunkBuffer)
    numOfReceivedSamples += chunk
  end

end

function collectSamples!(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, wpRead::Int64, chunk::Int64, channel::Channel, chunkBuffer)
  done = zeros(Bool, length(rpu))
  iterationDone = Condition()
  timeoutHappened = false
  result = zeros(Int16, numChan(rpu), chunk)
  performances = Vector{PerformanceData}(undef, length(rpu))

  for (d, rp) in enumerate(rpu)
    @async begin
      buffer = view(chunkBuffer, 1:(2 * chunk), d)
      (u, perf) = readSamplesChunk_(rp, Int64(wpRead), Int64(chunk), buffer)
      samples = reshape(u, 2, chunk)
      result[2*d-1, :] = samples[1, :]
      result[2*d, :] = samples[2, :]
      performances[d] = perf

      done[d] = true
      if (all(done))
        notify(iterationDone)
      end
    end
  end

  # Setup Timeout
  t = Timer(_timeout)
  @async begin
    wait(t)
    notify(iterationDone)
    timeoutHappened = true
  end

  wait(iterationDone)
  close(t)
  if timeoutHappened
    error("Timout reached when reading from sockets")
  end

  put!(channel, SampleChunk(result, performances))

end

# High level read. numFrames can adress a future frame. Data is read in
# chunks
function readFrames(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, startFrame, numFrames, numBlockAverages=1, numPeriodsPerPatch=1; rpInfo=nothing, chunkSize = 50000)
  numSampPerPeriod = samplesPerPeriod(rpu)
  numPeriods = periodsPerFrame(rpu)
  numSampPerFrame = numSampPerPeriod * numPeriods

  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end

  wpStart = startFrame * numSampPerFrame
  numOfRequestedSamples = numFrames * numSampPerFrame

  # rawSamples Int16 numofChan(rpc) x numOfRequestedSamples
  rawSamples = readPipelinedSamples(rpu, Int64(wpStart), Int64(numOfRequestedSamples), chunkSize = chunkSize, rpInfo = rpInfo)

  # Reshape/Avg Data
  data = convertSamplesToFrames(rawSamples, numChan(rpu), numSampPerPeriod, numPeriods, numFrames, numBlockAverages, numPeriodsPerPatch)

  return data
end

function convertSamplesToFrames(samples, numChan, numSampPerPeriod, numPeriods, numFrames, numBlockAverages=1, numPeriodsPerPatch=1)
  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end
  numTrueSampPerPeriod = div(numSampPerPeriod,numBlockAverages*numPeriodsPerPatch)
  frames = zeros(Float32, numTrueSampPerPeriod, numChan, numPeriods*numPeriodsPerPatch, numFrames)
  convertSamplesToFrames!(samples, frames, numChan, numSampPerPeriod, numPeriods, numFrames, numTrueSampPerPeriod, numBlockAverages, numPeriodsPerPatch)
  return frames
end

function convertSamplesToFrames!(samples, frames, numChan, numSampPerPeriod, numPeriods, numFrames, numTrueSampPerPeriod, numBlockAverages=1, numPeriodsPerPatch=1)
  temp = reshape(samples, numChan, numSampPerPeriod, numPeriods, numFrames)
  for d = 1:div(numChan,2)
    u = temp[2*d-1:2*d, :, :, :]
    utmp1 = reshape(u,2,numTrueSampPerPeriod,numBlockAverages, size(u,3)*numPeriodsPerPatch,size(u,4))
    utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=3) : utmp1
    frames[:,2*d-1,:,:] = utmp2[1,:,1,:,:]
    frames[:,2*d,:,:] = utmp2[2,:,1,:,:]
  end
end

function readPeriods(rpu::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, startPeriod, numPeriods, numBlockAverages=1; rpInfo=nothing, chunkSize = 50000)
  numSampPerPeriod = samplesPerPeriod(rpu)

  if rem(numSampPerPeriod,numBlockAverages) != 0
    error("block averages has to be a divider of numSampPerPeriod")
  end

  numAveragedSampPerPeriod = div(numSampPerPeriod,numBlockAverages)

  data = zeros(Float32, numAveragedSampPerPeriod, numChan(rpu), numPeriods)
  wpStart = startPeriod * numSampPerPeriod
  numOfRequestedSamples = numPeriods * numSampPerPeriod

  # rawSamples Int16 numofChan(rpc) x numOfRequestedSamples
  rawSamples = readPipelinedSamples(rpu, Int64(wpStart), Int64(numOfRequestedSamples), chunkSize = chunkSize, rpInfo = rpInfo)

  # Reshape/Avg Data
  convertSamplesToPeriods!(rawSamples, data, numChan(rpu), numSampPerPeriod, numPeriods, numBlockAverages)

  return data
end

function convertSamplesToPeriods!(samples, periods, numChan, numSampPerPeriod, numPeriods, numBlockAverages=1)
  temp = reshape(samples, numChan, numSampPerPeriod, numPeriods)
  for d = 1:div(numChan,2)
    u = temp[2*d-1:2*d, :, :]
    utmp1 = reshape(u,2,div(numSampPerPeriod,numBlockAverages), numBlockAverages, size(u,3))
    utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=3) : utmp1
    periods[:,2*d-1,:] = utmp2[1,:,1,:]
    periods[:,2*d,:] = utmp2[2,:,1,:]
  end
end

function readData(rpc::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, startFrame, numFrames, numBlockAverages=1, numPeriodsPerPatch=1; chunkSize = 50000)
  data = readFrames(rpc, startFrame, numFrames, numBlockAverages, numPeriodsPerPatch, chunkSize = chunkSize)
  return data
end

function readDataPeriods(rpc::Union{RedPitaya,RedPitayaCluster, RedPitayaClusterView}, startPeriod, numPeriods, numBlockAverages=1; chunkSize = 50000)
  data = readPeriods(rpc, startPeriod, numPeriods, numBlockAverages, chunkSize = chunkSize)
  return data
end

# High level read. numFrames can adress a future frame.
function readDataSlow(rpc::RedPitayaCluster, startFrame, numFrames)
  numPeriods = master(rpc).periodsPerFrame
  numChanTotal = numSlowADCChan(rpc)
  numSampPerFrame = numPeriods * numChanTotal

  data = zeros(Float32, numChanTotal, numPeriods, numFrames)
  wpRead = startFrame
  l=1

  # This is a wild guess for a good chunk size
  chunkSize = max(1,  round(Int, 1000000 / numSampPerFrame)  )
  @debug "chunkSize = $chunkSize"
  while l<=numFrames
    wpWrite = currentFrame(rpc)
    while wpRead >= wpWrite # Wait that startFrame is reached
      wpWrite = currentFrame(rpc)
      @debug wpWrite
    end
    chunk = min(wpWrite-wpRead,chunkSize) # Determine how many frames to read
    @debug chunk
    if l+chunk > numFrames
      chunk = numFrames - l + 1
    end

    @debug "Read from $wpRead until $(wpRead+chunk-1), WpWrite $(wpWrite), chunk=$(chunk)"

    p = 1
    for (d,rp) in enumerate(rpc)
      u = readDataSlow_(rp, Int64(wpRead), Int64(chunk))
      numChan = size(u,1)

      #utmp1 = reshape(u, numChan, numBlockAverages,
      #                    div(size(u,2),numBlockAverages),size(u,3))
      #utmp2 = numBlockAverages > 1 ? mean(utmp1,dims=2) : utmp1

      data[p:(p+numChan-1),:,l:(l+chunk-1)] = u #utmp2[:,1,:,:]
      p += numChan
    end

    l += chunk
    wpRead += chunk
  end

  return data
end

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
using LinearAlgebra
using Zygote: @adjoint

""" maximum eigenvalue of a symmetric square matrix """
function maxeig(symmat::Array{T,2})::T where {T}
    eig = eigen(Symmetric(symmat))
    max_idx = findmax(eig.values)[2]
    eigmax = eig.values[max_idx]
    return eigmax
end

@adjoint function maxeig(symmat::Array{T,2}) where {T}
    eig = eigen(Symmetric(symmat))
    max_idx = findmax(eig.values)[2]
    eigmax = eig.values[max_idx]
    vecmax = eig.vectors[:, max_idx]
    jac_adj = vecmax * (vecmax')
    return eigmax, c -> Tuple(jac_adj * c)
end

""" minimum eigenvalue of a symmetric square matrix """
function mineig(symmat::Array{T,2})::T where {T}
    eig = eigen(Symmetric(symmat))
    min_idx = findmin(eig.values)[2]
    eigmin = eig.values[min_idx]
    return eigmin
end

@adjoint function mineig(symmat::Array{T,2}) where {T}
    eig = eigen(Symmetric(symmat))
    min_idx = findmin(eig.values)[2]
    eigmin = eig.values[min_idx]
    vecmin = eig.vectors[:, min_idx]
    jac_adj = vecmin * (vecmin')
    return eigmax, c -> Tuple(jac_adj * c)
end

""" eigenvalue regularization """
function softplus_eigreg(symmat)
    eig = eigen(Symmetric(symmat))
    return mean(softplus.(eig.values))
end

@adjoint function softplus_eigreg(symmat)
    eig = eigen(Symmetric(symmat))
    return sum(softplus.(eig.values)),
    c -> Tuple(
        sum((
            sigmoid(eig.values[i]) * eig.vectors[:, i] * (eig.vectors[:, i]') for
            i = 1:size(symmat, 1)
        )) * c,
    )
end

function relu_eigreg(symmat)
    eig = eigen(Symmetric(symmat))
    return sum(relu.(eig.values))
end

@adjoint function relu_eigreg(symmat)
    eig = eigen(Symmetric(symmat))
    return sum(relu.(eig.values)),
    c -> Tuple(
        sum((
            0.5 * (sign(eig.values[i]) + 1) * eig.vectors[:, i] * (eig.vectors[:, i]')
            for i = 1:size(symmat, 1)
        )) * c,
    )
end

""" compute all eigenvalues """
function eigenvalues(symmat)
    eig = eigen(Symmetric(symmat))
    return eig.values
    # return eigvals(Symmetric(symmat))
end

@adjoint function eigenvalues(symmat)
    eig = eigen(Symmetric(symmat))
    return eig.values,
    c -> Tuple((eig.vectors[:, i] * (eig.vectors[:, i]')) * c for i = 1:size(symmat, 1))
    # vals = eigvals(Symmetric(symmat))
    # vecs = eigvecs(Symmetric(symmat))
    # return vals, c -> Tuple((vecs[:,i]*(vecs[:,i]'))*c for i=1:size(symmat,1))
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
@testset "size_pyramid" begin
    @test SinGAN.size_pyramid(4 / 3, (32, 32), (400, 300)) ==
        [(32, 32), (42, 42), (56, 56), (75, 75), (101, 101),
         (134, 134), (179, 179), (239, 239), (319, 300), (400, 300)]
    @test SinGAN.size_pyramid(4 / 3, (32, 32), (128, 128)) ==
        [(32, 32), (42, 42), (56, 56), (75, 75), (101, 101), (128, 128)]
    @code_warntype SinGAN.size_pyramid(4 / 3, (32, 32), (128, 128))
end

@testset "channel_pyramid" begin
    @test SinGAN.channel_pyramid(5) == [32, 32, 32, 64, 64]
    @code_warntype SinGAN.channel_pyramid(5)
end

@testset "resize_and_padding" begin
    x = randn(Float32, 2, 2, 3, 1)
    @info x
    @info SinGAN.resize_and_padding(x, (3, 3), (5, 5))
    # @code_warntype SinGAN.resize_and_padding(x, (3, 3), (5, 5))
end

@testset "build_zero_pyramid" begin
    xs = Float32[1 2; 3 4] |> gpu
    ts = [(1, 1), (2, 2)]
    @test SinGAN.build_zero_pyramid(xs, ts) == [zeros(Float32, 1, 1, 3, 1), zeros(Float32, 2, 2, 3, 1)] |> gpu
    # @code_warntype SinGAN.build_zero_pyramid(xs, ts)
end

@testset "build_noise_pyramid" begin
    xs = Float32[1 2; 3 4] |> gpu
    ts = [(1, 1), (2, 2)]
    as = [1f0, 0.5f0]
    @info SinGAN.build_noise_pyramid(xs, ts, as)
    # @code_warntype SinGAN.build_noise_pyramid(xs, ts, as)
end

@testset "rec_vector_generator" begin
    xs = Float32[1 2; 3 4] |> gpu
    ts = [(1, 1), (2, 2)]
    a = 1f0
    @info SinGAN.build_rec_pyramid(xs, ts, a)
    # @code_warntype SinGAN.build_rec_pyramid(xs, ts, a)
end

max_stars_repo_name: matsueushi/SinGAN

max_issues_repo_name: matsueushi/SinGAN

max_forks_repo_name: matsueushi/SinGAN


text:
module CreateTableBooks

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:books) do
    [
      primary_key()
      column(:title, :string)
      column(:author, :string)
    ]
  end

  add_index(:books, :title)
  add_index(:books, :author)
end

function down()
  drop_table(:books)
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
"""

    Soliton(K::Integer, M::Integer, δ::Real, atol::Real=0) <: Distribution{Univariate, Discrete}

The Robust Soliton distribution of length `K`, mode `M` (i.e., the
location of the robust component spike), peeling process failure
probability `δ`, and minimum non-zero probability mass `atol`. More
specifically, degrees `i` for which `pdf(Ω, i)<atol` are set to
`0`. Letting `atol=0` yields the regular robust Soliton distribution.

```julia
Soliton(K, M, δ)        # Robust Soliton distribution (with atol=0)
Soliton(K, M, δ, atol)  # Robust Soliton distribution with minimum non-zero probability mass atol

params(Ω)               # Get the parameters ,i.e., (K, M, δ, atol)
degrees(Ω)              # Return a vector composed of the degrees with non-zero probability mass
pdf(Ω, i)               # Evaluate the pdf at i
cdf(Ω, i)               # Evaluate the pdf at i
rand(Ω)                 # Sample from Ω
rand(Ω, n)              # Draw n samples from Ω
```

External links:

* [Soliton distribution on Wikipedia](https://en.wikipedia.org/wiki/Soliton_distribution)

"""
struct Soliton <: DiscreteUnivariateDistribution
    K::Int # Number of input symbols
    M::Int # Location of the robust component spike
    δ::Float64 # Peeling process failure probability
    atol::Float64 # Minimum non-zero probability assigned to a degree
    degrees::Vector{Int} # Degrees with non-zero probability
    CDF::Vector{Float64} # CDF evaluated at each element in degrees

    function Soliton(K::Integer, M::Integer, δ::Real, atol::Real=0)
        0 < K || throw(DomainError(K, "Expected 0 < K."))
        0 < δ < 1 || throw(DomainError(δ, "Expected 0 < δ < 1."))
        0 < M <= K || throw(DomainError(M, "Expected 0 < M <= K."))
        0 <= atol < 1 || throw(DomainError(atol, "Expected 0 <= atol < 1."))
        PDF = [soliton_τ(K, M, δ, i)+soliton_ρ(K, i) for i in 1:K]
        PDF ./= sum(PDF)
        degrees = [i for i in 1:K if PDF[i] > atol]
        CDF = cumsum([PDF[i] for i in degrees])
        CDF ./= CDF[end]
        new(K, M, δ, atol, degrees, CDF)
    end
end

Base.show(io::IO, Ω::Soliton) = print(io, "Soliton(K=$(Ω.K), M=$(Ω.M), δ=$(Ω.δ), atol=$(Ω.atol))")

"""
    degrees(Ω)

Return a vector composed of the degrees with non-zero probability.
"""
degrees(Ω::Soliton) = copy(Ω.degrees)

"""
Robust component of the Soliton distribution.
"""
function soliton_τ(K::Integer, M::Integer, δ::Real, i::Integer)
    i <= K || throw(DomainError(i, "Expected i <= K, but got i=$i, K=$K."))
    T = promote_type(typeof(δ), Float64)
    R = K / M
    if i < M
        return T(inv(i * M))
    elseif i == M
        return T(log(R / δ) / M)
    else # i <= K
        return zero(T)
    end
end

"""
Ideal component of the Soliton distribution.
"""
function soliton_ρ(K::Integer, i::Integer)
    i <= K || throw(DomainError(i, "Expected i <= K, but got i=$i, K=$K."))
    if i == 1
        return 1 / K
    else # i <= K
        return 1 / (i * (i - 1))
    end
end

StatsBase.params(Ω::Soliton) = (Ω.K, Ω.M, Ω.δ, Ω.atol)

function pdf(Ω::Soliton, i::Real)
    j = searchsortedfirst(Ω.degrees, i)
    (j > length(Ω.degrees) || Ω.degrees[j] != i) && return zero(eltype(Ω.CDF))
    rv = Ω.CDF[j]
    if j > 1
        rv -= Ω.CDF[j-1]
    end
    return rv
end
logpdf(Ω::Soliton, i::Real) = log(pdf(Ω, i))

function cdf(Ω::Soliton, i::Integer)
    i < Ω.degrees[1] && return 0.0
    i > Ω.degrees[end] && return 1.0
    j = searchsortedfirst(Ω.degrees, i)
    Ω.degrees[j] == i && return Ω.CDF[j]
    return Ω.CDF[j-1]
end

Statistics.mean(Ω::Soliton) = sum(i -> i*pdf(Ω, i), Ω.degrees)

function Statistics.var(Ω::Soliton)
    μ = mean(Ω)
    rv = 0.0
    for d in Ω.degrees
        rv += pdf(Ω, d)*(d-μ)^2
    end
    return rv
end

function Statistics.quantile(Ω::Soliton, v::Real)
    0 <= v <= 1 || throw(DomainError(v, "Expected 0 <= v <= 1."))
    j = searchsortedfirst(Ω.CDF, v)
    j = min(length(Ω.degrees), j)
    return Ω.degrees[j]
end

Base.minimum(Ω::Soliton) = 1
Base.maximum(Ω::Soliton) = Ω.K


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

test_resource_manager = load_op_and_grad("./build/libTestResourceManager","test_resource_manager")
# TODO: specify your input parameters
k = constant(0, dtype=Int32)
u1 = test_resource_manager(k)
u2 = test_resource_manager(k)
u3 = test_resource_manager(k)

control_dependencies([u1, u2, u3]) do
    global z = test_resource_manager(k)  # use z to read the summation
end

# z = test_resource_manager(k)
sess = Session(); init(sess)
run(sess, z)

# output
# Create a new container
# Current Value=1
# Current Value=2
# Current Value=3
# Current Value=4
# 4


text:
module CellValuesTests

using Test

include("CellValuesMocks.jl")

include("MapsMocks.jl")

@testset "CellValuesMocks" begin include("CellValuesMocksTests.jl") end

@testset "MapsMocks" begin include("MapsMocksTests.jl") end

@testset "CellNumbers" begin include("CellNumbersTests.jl") end

@testset "CellArrays" begin include("CellArraysTests.jl") end

@testset "CellMaps" begin include("CellMapsTests.jl") end

@testset "Kernels" begin include("KernelsTests.jl") end

@testset "CellNumberApply" begin include("CellNumberApplyTests.jl") end

@testset "CellArrayApply" begin include("CellArrayApplyTests.jl") end

@testset "MapApply" begin include("MapApplyTests.jl") end

@testset "CellMapApply" begin include("CellMapApplyTests.jl") end

@testset "CellValuesOperations" begin include("CellValuesOperationsTests.jl") end

@testset "CellValuesGallery" begin include("CellValuesGalleryTests.jl") end

@testset "CellValuesAppend" begin include("CellValuesAppendTests.jl") end

@testset "CellValuesReindex" begin include("CellValuesReindexTests.jl") end

@testset "CellNumberReindex" begin include("CellNumberReindexTests.jl") end

@testset "CellArrayReindex" begin include("CellArrayReindexTests.jl") end

@testset "FindLocalIndex" begin include("FindLocalIndexTests.jl") end

@testset "ConstantCellValues" begin include("ConstantCellValuesTests.jl") end

@testset "CompressedCellValues" begin include("CompressedCellValuesTests.jl") end

@testset "NonIterableCellMaps" begin include("NonIterableCellMapsTests.jl") end

@testset "IdentityCellNumbers" begin include("IdentityCellNumbersTests.jl") end

end # module


text:
abstract type PRAlphaModel <: AlphaModel end

struct PRAlphaParam <: EoSParam
    acentricfactor::SingleParam{Float64}
end

@newmodelsimple PRAlpha PRAlphaModel PRAlphaParam

export PRAlpha
function PRAlpha(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = PRAlphaParam(acentricfactor)
    model = PRAlpha(packagedparams, verbose=verbose)
    return model
end

function α_function(model::CubicModel,V,T,z,alpha_model::PRAlphaModel)
    Tc = model.params.Tc.values
    ω  = alpha_model.params.acentricfactor.values
    α = zeros(typeof(T),length(Tc))
    for i in @comps
        ωi = ω[i]
        Tr = T/Tc[i]
        m = evalpoly(ωi,(0.37464,1.54226,-0.26992))
        α[i] = (1+m*(1-√(Tr)))^2
    end
    return α
end


text:
# Autogenerated wrapper script for GAP_pkg_edim_jll for x86_64-linux-gnu
export ediv

using GAP_jll
using GAP_lib_jll
JLLWrappers.@generate_wrapper_header("GAP_pkg_edim")
JLLWrappers.@declare_file_product(ediv)
function __init__()
    JLLWrappers.@generate_init_header(GAP_jll, GAP_lib_jll)
    JLLWrappers.@init_file_product(
        ediv,
        "lib/gap/ediv.so",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import SemialgebraicSets
import PolyJuMP

@testset "Non-symmetric matrix SOS constraint" begin
    @polyvar x
    m = SOSModel()
    # Throws an Argument Error because it should be symmetric to be an SOS matrix
    @test_throws ArgumentError PolyJuMP.addpolyconstraint!(m, [1 x; -x 0], SOSMatrixCone(), SemialgebraicSets.BasicSemialgebraicSet{Int, polynomialtype(x, Int)}(), PolyJuMP.MonomialBasis)
end


text:
using LightGraphsFlows, LightGraphs

flow_graph = DiGraph(8)
# (node,node,flow)
flow_edges = [(1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
              (2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
              (5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)];
capacity_matrix = zeros(Int ,nv(flow_graph) ,nv(flow_graph));
for e in flow_edges
    u,v,f = e
    LightGraphs.add_edge!(flow_graph,u,v)
    capacity_matrix[u,v] = f
end
xxx = LightGraphs.DiGraph(LightGraphs.Graph(flow_graph))
LightGraphsFlows.boykov_kolmogorov_impl(xxx, 1, 8, capacity_matrix)

# pkg> activate .
using BKMaxflow

weights, neighbors = create_graph(JuliaImpl{Float64,Int}, 8)

BKMaxflow.add_edge!(weights, neighbors, 1, 2, 10., 0.)
BKMaxflow.add_edge!(weights, neighbors, 1, 3, 5., 0.)
BKMaxflow.add_edge!(weights, neighbors, 1, 4, 15., 0.)
BKMaxflow.add_edge!(weights, neighbors, 2, 3, 4., 0.)
BKMaxflow.add_edge!(weights, neighbors, 2, 5, 9., 0.)
BKMaxflow.add_edge!(weights, neighbors, 2, 6, 15., 0.)
BKMaxflow.add_edge!(weights, neighbors, 3, 4, 4., 0.)
BKMaxflow.add_edge!(weights, neighbors, 3, 6, 8., 0.)
BKMaxflow.add_edge!(weights, neighbors, 4, 7, 16., 0.)
BKMaxflow.add_edge!(weights, neighbors, 5, 6, 15., 0.)
BKMaxflow.add_edge!(weights, neighbors, 5, 8, 10., 0.)
BKMaxflow.add_edge!(weights, neighbors, 6, 7, 15., 0.)
BKMaxflow.add_edge!(weights, neighbors, 6, 8, 10., 0.)
BKMaxflow.add_edge!(weights, neighbors, 7, 3, 6., 0.)
BKMaxflow.add_edge!(weights, neighbors, 7, 8, 10., 0.)

flow, label = boykov_kolmogorov(1, 8, neighbors, weights)

# c wrapper
g = create_graph(CImpl{Cdouble}, 8, 15)

add_node(CImpl{Cdouble}, g, 6)

add_tweights(CImpl, g, 0, 10., 0.)
add_tweights(CImpl, g, 1, 5., 0.)
add_tweights(CImpl, g, 2, 15., 0.)
add_tweights(CImpl, g, 3, 0., 10.)
add_tweights(CImpl, g, 4, 0., 10.)
add_tweights(CImpl, g, 5, 0., 10.)

add_edge(CImpl, g, 0, 1, 4., 0.)
add_edge(CImpl, g, 0, 3, 9., 0.)
add_edge(CImpl, g, 0, 4, 15., 0.)
add_edge(CImpl, g, 1, 2, 4., 0.)
add_edge(CImpl, g, 1, 4, 8., 0.)
add_edge(CImpl, g, 2, 5, 16., 0.)
add_edge(CImpl, g, 3, 4, 15., 0.)
add_edge(CImpl, g, 4, 5, 15., 0.)
add_edge(CImpl, g, 5, 1, 6., 0.)

flow = maxflow(CImpl{Cdouble}, g)  #-> 3

what_segment(CImpl{Cdouble}, g, 0)
what_segment(CImpl{Cdouble}, g, 1)
what_segment(CImpl{Cdouble}, g, 2)
what_segment(CImpl{Cdouble}, g, 3)
what_segment(CImpl{Cdouble}, g, 4)
what_segment(CImpl{Cdouble}, g, 5)

delete_graph(CImpl{Cdouble}, g)

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
using Test
using ShowSet

A = BitSet([1,3,9])
@test string(A) == "{1,3,9}"

B = Set{String}()
push!(B,"Hello")
push!(B,"Bye")
@test "{Bye,Hello}" == string(B)

set_empty("{}")
Z = Set([])
@test string(Z) == get_empty()

A = Set([-im, im])
B = Set([im, -im])
@test string(A) == string(B)


text:
using Winston
p = plot(abs(sin(0:.1:10)))
setattr(p, "ylog", true)
setattr(p, "yrange", (5e-6,0.2))


text:
# This file investigates integration error with respect to number of points in 1D integration

using FractalTools
using Makie

# Construct interpolation data
f(x) = sin(2π * x) + 1
xi = 0.
xf = 1.
freevar = 0.001
npts = 10 : 2 : 50
ntpts = 10 * npts[end]

# Compute theoretical value
fval = (cos(2π * xi) - cos(2π * xf)) / (2π) + (xf - xi)

# Compute errors
mse = map(npts) do npt
    # Construct interpolation data
    x = collect(range(xi, xf, length=npt))
    y = f.(x)
    pts = collect.(zip(x, y))

    # Compute integration
    ival = integrate(pts, Interp1D(freevar))

    # Compute error
    abs(fval - ival)
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Number of Points", ylabel="Error", title="1D Integration Error")
stem!(ax, npts, mse, color=:black)
save(joinpath(@__DIR__, "integ1d_error.png"), fig)
display(fig)


text:
using PQEcolaPoint
using Dates

# build a PUT contract -
put_options = Dict{String,Any}()
put_options["ticker"] = "AMD"
put_options["expiration_date"] = Date(2022, 04, 14)
put_options["strike_price"] = 105.0
put_options["premium"] = 3.59
put_options["number_of_contracts"] = 1
put_options["direction"] = -1
put_options["current_price"] = 2.87
put_contract_model = build(PutContractModel, put_options)

# build a CALL contract -
call_options = Dict{String,Any}()
call_options["ticker"] = "AMD"
call_options["expiration_date"] = Date(2022, 04, 14)
call_options["strike_price"] = 130.0
call_options["premium"] = 302.0
call_options["number_of_contracts"] = 1
call_options["direction"] = -1
call_options["current_price"] = 6.28
call_contract_model = build(CallContractModel, call_options)

# build an EQUITY model -
equity_options = Dict{String,Any}()
equity_options["ticker"] = "AMD"
equity_options["purchase_price"] = 106.0
equity_options["current_price"] = 122.99
equity_options["direction"] = 1
equity_options["number_of_shares"] = 1
equity_contract_model = build(EquityModel, equity_options)


text:
using DeepCFR
using Test
using CounterfactualRegret
using CounterfactualRegret.Games
const CFR = CounterfactualRegret
using StaticArrays
using JLD2
using Flux
using Random

Random.seed!(1337)

include(joinpath(@__DIR__, "solver.jl"))

include(joinpath(@__DIR__, "fitting.jl"))

include(joinpath(@__DIR__, "memory.jl"))

include(joinpath(@__DIR__, "check.jl"))

include(joinpath(@__DIR__, "callback.jl"))


text:
@testset "User defined models" begin
    y = [1 2; 1 2.0]
    Z = [1 2.0; 2.0 1.0]
    T = [1 2.0]
    R = [1 2.0]
    @test_throws ErrorException StateSpaceModel(y, Z, T, R)

    Z = Array{Float64, 3}(undef, 2, 2, 2)
    Z[:, :, 1] = [1 2.0; 2.0 1.0]
    Z[:, :, 2] = [1 2.0; 2.0 1.0]
    T = [1 2.0; 2.0 1.0]
    R = [1 2.0; 2.0 1.0]
    ssm = StateSpaceModel(y, Z, T, R)
    @test ssm.mode == "time-variant"

    Z = Array{Float64, 3}(undef, 2, 2, 2)
    Z[:, :, 1] = [1 2.0; 2.0 1.0]
    Z[:, :, 2] = [1 2.0; 2.0 1.0]
    T = [1 2.0]
    R = [1 2.0; 2.0 1.0]
    @test_throws ErrorException StateSpaceModel(y, Z, T, R)

    y = ones(5, 1)
    Z = ones(1, 3)
    T = ones(3, 3)
    R = ones(3, 3)
    d = zeros(5, 1) # correct d dimension
    c = zeros(5, 3) # correct c dimension
    H = SSM.build_H(1, Float64)
    Q = SSM.build_Q(3, 1, Float64)
    ssm = StateSpaceModel(y, Z, T, R, d, c, H, Q)
    @test ssm.dim == StateSpaceDimensions(5, 1, 3, 3)

    d = zeros(4, 1) # incorrect d dimension
    @test_throws AssertionError StateSpaceModel(y, Z, T, R, d, c, H, Q)

    d = zeros(5, 1) # correct d dimension
    c = zeros(5, 2) # incorrect c dimension
    @test_throws AssertionError StateSpaceModel(y, Z, T, R, d, c, H, Q)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE
abstract type AbstractTypeDBStub end

mutable struct Core_TypeDBStub <: AbstractTypeDBStub
    blockingStub::Proto.TypeDBBlockingStub
    asyncStub::Proto.TypeDBStub
end

function Core_TypeDBStub(channel::gRPCClient.gRPCChannel)
    blockingStub = Proto.TypeDBBlockingStub(channel)
    asyncStub = Proto.TypeDBStub(channel)
    return Core_TypeDBStub(blockingStub,asyncStub)
end

mutable struct Cluster_TypeDBStub <: AbstractTypeDBStub
    blockingStub::Proto.TypeDBBlockingStub
end

function Cluster_TypeDBStub(channel::gRPCClient.gRPCChannel)
    blockingStub = Proto.TypeDBBlockingStub(channel)
    return Cluster_TypeDBStub(blockingStub)
end

function ensure_connected(stub::T) where {T<:AbstractTypeDBStub}
    throw(TypeDBClientException(GENERAL_UNKOWN_ERROR,"Isn't implemented yet"))
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export typeofweight
"""
`typeofweight(fst_or_path)`

Returns the weight type of a WFST or a path.
"""
typeofweight(fst::Union{WFST{W},Path{W}}) where {W} = W

export get_ialphabet
"""
`get_ialphabet(label::fst)`

Returns the alphabet of the input labels of `fst`.
The alphabet is defined as the set of symbols composing the input or output labels.
"""
get_ialphabet(fst::WFST) = Set(keys(get_isym(fst)))

export get_oalphabet
"""
`get_oalphabet(label::fst)`

Returns the alphabet of the output labels of `fst`.
The alphabet is defined as the set of symbols composing the input or output labels.
"""
get_oalphabet(fst::WFST) = Set(keys(get_osym(fst)))

export count_eps
"""
`count_eps(label::Function, fst::FST)`

Returns number of epsilon transition of either input or output labels.
To specify input or output plug in `ilabel` and `olabel` respectively in `label`.
"""
function count_eps(label::Function, fst::WFST)
  cnt = 0
  for state in fst
    for arc in state
      if iseps(label(arc))
        cnt += 1
      end
    end
  end
  return cnt
end

export has_eps
"""
`has_eps([label::Function,] fst::FST)`

Check if `fst` has epsilon transition in either input or output labels.
To specify input or output plug in either `get_ilabel` or `get_olabel` respectively in `label`.
"""
function has_eps(label::Function, fst::WFST)
  for state in fst
    for arc in state
      if iseps(label(arc))
        return true
      end
    end
  end
  return false
end

function has_eps(fst::WFST)
  for state in fst
    for arc in state
      if iseps(get_ilabel(arc))
        return true
      end
      if iseps(get_olabel(arc))
        return true
      end
    end
  end
  return false
end

export is_acceptor
"""
`is_acceptor(fst::WFST)`

Returns `true` if `fst` is an acceptor, i.e. if all arcs have equal input and output labels.
"""
function is_acceptor(fst::WFST)
  acceptor = true
  for s in fst
    for a in s
      acceptor = acceptor && (get_ilabel(a) == get_olabel(a))
      if acceptor == false
        break
      end
    end
  end
  return acceptor
end

export is_deterministic
"""
`is_deterministic(fst::WFST; get_label=get_ilabel)`

Returns `true` if `fst` is deterministic in the input labels.
A input label deterministic WFST must have a single initial state and arcs leaving any state do not share the same input label.

Change the keyword `get_label` to `get_olabel` in order to check determinism in the output labels.
"""
function is_deterministic(fst::WFST; get_label=get_ilabel)
  if length(get_initial(fst; single=false)) > 1
    return false
  end
  for s in fst
    for i in 1:length(s)
      for ii in i+1:length(s)
        if get_label(s[i]) == get_label(s[ii])
          return false
        end
      end
    end
  end
  return true
end

export is_acyclic
"""
  `is_acyclic(fst::WFST[,i = get_initial(fst;single=true)]; kwargs...)`

Returns the number of states of the fst. For `kwargs` see [`DFS`](@ref).
"""
function is_acyclic(fst::WFST, i = get_initial(fst;single=true); kwargs...)
  if isempty(fst)
    return true
  else
    dfs = DFS(fst,i; kwargs...)

    for (p,s,n,d,e,a) in dfs
      if !d # already discovered states
        if !e && !dfs.is_completed[n]
          return false
        end
      end
    end
    return true
  end
end

"""
`length(fst::WFST)`

Returns the number of states of the fst.
"""
Base.length(fst::WFST) = length(get_states(fst))

"""
`size(fst::WFST,[i])`

Returns number of states and total number of arcs in a tuple.
If 'i=1' returns the number of states and if `i=2` the number of arcs.
"""
Base.size(fst::WFST) = (length(get_states(fst)), count_arcs(fst))
Base.size(fst::WFST,i::Int) = size(fst)[i]


text:
const BACKEND = :CUDA
using Pkg,KernelAbstractions, Adapt, OffsetArrays,CUDAKernels
Pkg.activate(string(BACKEND, "Env"))
using CUDA, CUDAKernels

@kernel function saxpy!(z, α, x, y)
    I = @index(Global)
    @inbounds z[I] = α * x[I] + y[I]
end

"""
The workgroupsize is a local block of threads that are co-executed. The ndrange specifies the global index space
. This index space will be subdivided by the workgroupsize and each group will be executed in parallel.
"""
ndrange = (128,)
workgroupsize = (16,)
blocks, workgroupsize, dynamic = KernelAbstractions.NDIteration.partition(ndrange, workgroupsize)

"""
Initialization
"""
kernel = saxpy!(KernelAbstractions.Device())
x = adapt(ArrayT, rand(64, 32))
y = adapt(ArrayT, rand(64, 32))
z = similar(x)


text:
using Test
using Spandex

@testset "from_csr 1" begin
    local a = from_csr(
        3,
        4,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    @test 3 == a.row_count
    @test 4 == a.column_count
    @test 11 == a.nnz

    @test 10 == get_rowwise(a, 1, 1)
    @test 11 == get_rowwise(a, 1, 2)
    @test 12 == get_rowwise(a, 1, 3)
    @test 13 == get_rowwise(a, 1, 4)
    @test 23 == get_rowwise(a, 2, 4)
    @test 20 == get_rowwise(a, 2, 1)
    @test 21 == get_rowwise(a, 2, 2)
    @test 33 == get_rowwise(a, 3, 4)
    @test 32 == get_rowwise(a, 3, 3)
    @test 31 == get_rowwise(a, 3, 2)
    @test 30 == get_rowwise(a, 3, 1)
end

@testset "from_graph 1" begin
    local a = from_csr(
        3,
        4,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    local graph = DirectedGraph{Int64}(3)
    graph[1, 1] = 10
    graph[1, 2] = 11
    graph[1, 3] = 12
    graph[1, 4] = 13
    graph[2, 4] = 23
    graph[2, 1] = 20
    graph[2, 2] = 21
    graph[3, 4] = 33
    graph[3, 3] = 32
    graph[3, 2] = 31
    graph[3, 1] = 30

    local b = from_graph(graph, 3, 4)

    @test equals(a, b)
    @test equals(b, a)
end

@testset "transpose 1" begin
    local a = from_csr(
        3,
        4,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    local at = transpose(a)
    local att = transpose(at)
    local attt = transpose(att)

    @test equals(a, att)
    @test equals(at, attt)
    @test !equals(a, at)
end

@testset "get_row 1" begin
    local a = from_csr(
        4,
        4,
        [1, 5, 8, 12, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    local row = get_row(a, 1)
    @test 4 == row.size
    @test 4 == row.nnz
    @test 10 == row[1]
    @test 11 == row[2]
    @test 12 == row[3]
    @test 13 == row[4]

    row = get_row(a, 2)
    @test 4 == row.size
    @test 3 == row.nnz
    @test 20 == row[1]
    @test 21 == row[2]
    @test 23 == row[4]

    row = get_row(a, 3)
    @test 4 == row.size
    @test 4 == row.nnz
    @test 30 == row[1]
    @test 31 == row[2]
    @test 32 == row[3]
    @test 33 == row[4]

    row = get_row(a, 4)
    @test 4 == row.size
    @test 0 == row.nnz
end

@testset "get_column 1" begin
    local a = from_csr(
        3,
        5,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    local column = get_column(a, 1)
    @test 3 == column.size
    @test 3 == column.nnz
    @test 10 == column[1]
    @test 20 == column[2]
    @test 30 == column[3]

    column = get_column(a, 2)
    @test 3 == column.size
    @test 3 == column.nnz
    @test 11 == column[1]
    @test 21 == column[2]
    @test 31 == column[3]

    column = get_column(a, 3)
    @test 3 == column.size
    @test 2 == column.nnz
    @test 12 == column[1]
    @test 32 == column[3]

    column = get_column(a, 4)
    @test 3 == column.size
    @test 3 == column.nnz
    @test 13 == column[1]
    @test 23 == column[2]
    @test 33 == column[3]

    column = get_column(a, 5)
    @test 3 == column.size
    @test 0 == column.nnz
end

@testset "get 1" begin
    local a = from_csr(
        3,
        4,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [10, 11, 12, 13, 20, 21, 23, 30, 31, 32, 33],
    )

    @test 10 == get_rowwise(a, 1, 1)
    @test 11 == get_rowwise(a, 1, 2)
    @test 12 == get_rowwise(a, 1, 3)
    @test 13 == get_rowwise(a, 1, 4)
    @test 23 == get_rowwise(a, 2, 4)
    @test 20 == get_rowwise(a, 2, 1)
    @test 21 == get_rowwise(a, 2, 2)
    @test 33 == get_rowwise(a, 3, 4)
    @test 32 == get_rowwise(a, 3, 3)
    @test 31 == get_rowwise(a, 3, 2)
    @test 30 == get_rowwise(a, 3, 1)

    @test get_columnwise(a, 1, 1) == get_rowwise(a, 1, 1)
    @test get_columnwise(a, 1, 2) == get_rowwise(a, 1, 2)
    @test get_columnwise(a, 1, 3) == get_rowwise(a, 1, 3)
    @test get_columnwise(a, 1, 4) == get_rowwise(a, 1, 4)
    @test get_columnwise(a, 2, 4) == get_rowwise(a, 2, 4)
    @test get_columnwise(a, 2, 1) == get_rowwise(a, 2, 1)
    @test get_columnwise(a, 2, 2) == get_rowwise(a, 2, 2)
    @test get_columnwise(a, 3, 4) == get_rowwise(a, 3, 4)
    @test get_columnwise(a, 3, 3) == get_rowwise(a, 3, 3)
    @test get_columnwise(a, 3, 2) == get_rowwise(a, 3, 2)
    @test get_columnwise(a, 3, 1) == get_rowwise(a, 3, 1)
end

@testset "set! 1" begin
    local a = from_csr(
        3,
        4,
        [1, 5, 8, 12],
        [1, 2, 3, 4, 1, 2, 4, 1, 2, 3, 4],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    )

    set_rowwise!(a, 1, 1, 10)
    set_rowwise!(a, 1, 2, 11)
    set_rowwise!(a, 1, 3, 12)
    set_rowwise!(a, 1, 4, 13)
    set_rowwise!(a, 2, 4, 23)
    set_rowwise!(a, 2, 1, 20)
    set_rowwise!(a, 2, 2, 21)
    set_rowwise!(a, 3, 4, 33)
    set_rowwise!(a, 3, 3, 32)
    set_rowwise!(a, 3, 2, 31)
    set_rowwise!(a, 3, 1, 30)

    @test 10 == get_rowwise(a, 1, 1)
    @test 11 == get_rowwise(a, 1, 2)
    @test 12 == get_rowwise(a, 1, 3)
    @test 13 == get_rowwise(a, 1, 4)
    @test 23 == get_rowwise(a, 2, 4)
    @test 20 == get_rowwise(a, 2, 1)
    @test 21 == get_rowwise(a, 2, 2)
    @test 33 == get_rowwise(a, 3, 4)
    @test 32 == get_rowwise(a, 3, 3)
    @test 31 == get_rowwise(a, 3, 2)
    @test 30 == get_rowwise(a, 3, 1)

    set_columnwise!(a, 1, 1, 100)
    set_columnwise!(a, 1, 2, 110)
    set_columnwise!(a, 1, 3, 120)
    set_columnwise!(a, 1, 4, 130)
    set_columnwise!(a, 2, 4, 230)
    set_columnwise!(a, 2, 1, 200)
    set_columnwise!(a, 2, 2, 210)
    set_columnwise!(a, 3, 4, 330)
    set_columnwise!(a, 3, 3, 320)
    set_columnwise!(a, 3, 2, 310)
    set_columnwise!(a, 3, 1, 300)

    @test 100 == get_rowwise(a, 1, 1)
    @test 110 == get_rowwise(a, 1, 2)
    @test 120 == get_rowwise(a, 1, 3)
    @test 130 == get_rowwise(a, 1, 4)
    @test 230 == get_rowwise(a, 2, 4)
    @test 200 == get_rowwise(a, 2, 1)
    @test 210 == get_rowwise(a, 2, 2)
    @test 330 == get_rowwise(a, 3, 4)
    @test 320 == get_rowwise(a, 3, 3)
    @test 310 == get_rowwise(a, 3, 2)
    @test 300 == get_rowwise(a, 3, 1)
end

@testset "contains 1" begin
    local a = from_csr(4, 3, [1, 2, 3, 5, 5], [2, 1, 2, 3], [1, 1, 1, 1])

    @test contains(a, 1, 2)
    @test contains(a, 2, 1)
    @test contains(a, 3, 2)
    @test contains(a, 3, 3)

    @test !contains(a, 1, 1)
    @test !contains(a, 1, 3)
    @test !contains(a, 2, 2)
    @test !contains(a, 2, 3)
    @test !contains(a, 3, 1)
    @test !contains(a, 4, 1)
    @test !contains(a, 4, 2)
    @test !contains(a, 4, 3)
end


text:
d = @destruct [a, b] = Dict(:a => 1, :b => 2)
@test d == Dict(:a => 1, :b => 2)
@test (a, b) == (1, 2)

@destruct [a] = Dict("a" => "foo")
@test a == "foo"

@destruct [foo = :a || 5, b = :b || 6, c || 7] = Dict(:a => 1)
@test (foo, b, c) == (1, 6, 7)

import Base: ==

struct S
    data::Vector{UInt8}
    S(s) = new(codeunits(s))
end
==(s1::S, s2::S) = s1.data == s2.data

@destruct x.data = S("foo")
@test x == S("foo")
@test data == S("foo").data

@destruct _.(re, im) = Complex(1,2)
@test (re, im) == (1, 2)

@destruct [s.data = :a] = Dict(:a => S("foo"))
@test s == S("foo")
@test data === s.data

@destruct x[a, [c, d] = b] = Dict(:a => 1, :b => Dict(:c => 2, :d => 3))
@test x == Dict(:a => 1, :b => Dict(:c => 2, :d => 3))
@test (a, c, d) == (1, 2, 3)
@test b == Dict(:c => 2, :d => 3)


text:
function addsub(x, y)
  return x + y, x - y
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractDijkstraState <: AbstractPathState end

struct DijkstraHeapEntry{T}
    vertex::Int
    dist::T
end

isless(e1::DijkstraHeapEntry, e2::DijkstraHeapEntry) = e1.dist < e2.dist

mutable struct DijkstraState{T}<: AbstractDijkstraState
    parents::Vector{Int}
    dists::Vector{T}
    predecessors::Vector{Vector{Int}}
    pathcounts::Vector{Int}
end

"""
    dijkstra_shortest_paths(g, s, distmx=weights(g); allpaths=false)
    dijkstra_shortest_paths(g, sources, distmx=weights(g); allpaths=false)

Performs [Dijkstra's algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)
on a graph, computing shortest distances between a source vertex `s` (or a vector
`sources`)  and all other veritces.
Returns a `DijkstraState` that contains various traversal information.

With `allpaths=true`, returns a `DijkstraState` that keeps track of all
predecessors of a given vertex.
"""
function dijkstra_shortest_paths(
        g::AGraphOrDiGraph,
        srcs::Vector{Int},
        distmx::AEdgeMap=weights(g);
        allpaths=false
    )
    T = valtype(distmx)
    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(Int, nvg)
    preds = fill(Vector{Int}(),nvg)
    visited = zeros(Bool, nvg)
    pathcounts = zeros(Int, nvg)
    H = Vector{DijkstraHeapEntry{T}}()  # this should be Vector{T}() in 0.4, I think.
    dists[srcs] .= 0
    pathcounts[srcs] .= 1

    sizehint!(H, nvg)

    for v in srcs
        heappush!(H, DijkstraHeapEntry{T}(v, dists[v]))
    end

    while !isempty(H)
        hentry = heappop!(H)
        # info("Popped H - got $(hentry.vertex)")
        u = hentry.vertex
        for e in out_edges(g,u)
            v = dst(e)
            alt = (dists[u] == typemax(T)) ? typemax(T) : dists[u] + distmx[e]

            if !visited[v]
                dists[v] = alt
                parents[v] = u
                pathcounts[v] += pathcounts[u]
                visited[v] = true
                if allpaths
                    preds[v] = [u;]
                end
                heappush!(H, DijkstraHeapEntry{T}(v, alt))
                # info("Pushed $v")
            else
                if alt < dists[v]
                    dists[v] = alt
                    parents[v] = u
                    heappush!(H, DijkstraHeapEntry{T}(v, alt))
                end
                if alt == dists[v]
                    pathcounts[v] += pathcounts[u]
                    if allpaths
                        push!(preds[v], u)
                    end
                end
            end
        end
    end

    dists[srcs] .= 0
    pathcounts[srcs] .= 1
    parents[srcs] .= 0
    for src in srcs
        preds[src] = []
    end

    return DijkstraState{T}(parents, dists, preds, pathcounts)
end

dijkstra_shortest_paths(g::AGraphOrDiGraph, src::Int, distmx::AEdgeMap=weights(g); allpaths=false) =
    dijkstra_shortest_paths(g, [src;], distmx; allpaths=allpaths)