text:
using Weave
using Test

cleanup = true

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="tex")
result =  read("documents/pyplot_formats.tex", String)
ref =  read("documents/pyplot_formats_ref.tex", String)
result = replace(result, r"\s*PyObject.*\n", "\n") #Remove PyObjects, because they change
ref = replace(ref, r"\s*PyObject.*\n", "\n")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="github")
result =  read("documents/pyplot_formats.md", String)
ref =  read("documents/pyplot_formats_ref.md", String)
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

weave("documents/pyplot_formats.txt", plotlib="pyplot", doctype="rst", fig_ext=".svg")
result =  read("documents/pyplot_formats.rst", String)
ref =  read("documents/pyplot_formats_ref.rst", String)
result = replace(result, r"\s*PyObject.*\n", "")
ref = replace(ref, r"\s*PyObject.*\n", "")
@test result == ref

if cleanup
    rm("documents/pyplot_formats.tex")
    rm("documents/pyplot_formats.rst")
    rm("documents/pyplot_formats.md")
    rm("documents/figures", recursive = true)
end


text:
#
# Types
#

# BSONType mirrors C enum bson_type_t.
primitive type BSONType 8 end # 1 byte

Base.convert(::Type{T}, t::BSONType) where {T<:Number} = reinterpret(UInt8, t)
Base.convert(::Type{BSONType}, n::T) where {T<:Number} = reinterpret(BSONType, n)
BSONType(u::UInt8) = convert(BSONType, u)

#
# Constants for BSONType
#

const BSON_TYPE_EOD        = BSONType(0x00)
const BSON_TYPE_DOUBLE     = BSONType(0x01)
const BSON_TYPE_UTF8       = BSONType(0x02)
const BSON_TYPE_DOCUMENT   = BSONType(0x03)
const BSON_TYPE_ARRAY      = BSONType(0x04)
const BSON_TYPE_BINARY     = BSONType(0x05)
const BSON_TYPE_UNDEFINED  = BSONType(0x06)
const BSON_TYPE_OID        = BSONType(0x07)
const BSON_TYPE_BOOL       = BSONType(0x08)
const BSON_TYPE_DATE_TIME  = BSONType(0x09)
const BSON_TYPE_NULL       = BSONType(0x0A)
const BSON_TYPE_REGEX      = BSONType(0x0B)
const BSON_TYPE_DBPOINTER  = BSONType(0x0C)
const BSON_TYPE_CODE       = BSONType(0x0D)
const BSON_TYPE_SYMBOL     = BSONType(0x0E)
const BSON_TYPE_CODEWSCOPE = BSONType(0x0F)
const BSON_TYPE_INT32      = BSONType(0x10)
const BSON_TYPE_TIMESTAMP  = BSONType(0x11)
const BSON_TYPE_INT64      = BSONType(0x12)
const BSON_TYPE_DECIMAL128 = BSONType(0x13)
const BSON_TYPE_MAXKEY     = BSONType(0x7F)
const BSON_TYPE_MINKEY     = BSONType(0xFF)

# BSONSubType mirrors C enum bson_subtype_t.
primitive type BSONSubType 8 end

Base.convert(::Type{T}, t::BSONSubType) where {T<:Number} = reinterpret(UInt8, t)
Base.convert(::Type{BSONSubType}, n::T) where {T<:Number} = reinterpret(BSONSubType, n)
BSONSubType(u::UInt8) = convert(BSONSubType, u)

#
# Constants for BSONSubType
#

const BSON_SUBTYPE_BINARY            = BSONSubType(0x00)
const BSON_SUBTYPE_FUNCTION          = BSONSubType(0x01)
const BSON_SUBTYPE_BINARY_DEPRECATED = BSONSubType(0x02)
const BSON_SUBTYPE_UUID_DEPRECATED   = BSONSubType(0x03)
const BSON_SUBTYPE_UUID              = BSONSubType(0x04)
const BSON_SUBTYPE_MD5               = BSONSubType(0x05)
const BSON_SUBTYPE_USER              = BSONSubType(0x80)

#=
BSONIter mirrors C struct bson_iter_t and can be allocated in the stack.

According to [libbson documentation](http://mongoc.org/libbson/current/bson_iter_t.html),
it is meant to be used on the stack and can be discarded at any time
as it contains no external allocation.
The contents of the structure should be considered private
and may change between releases, however the structure size will not change.

Inspecting its size in C, we get:

```c
sizeof(bson_iter_t) == 80
```
=#
primitive type BSONIter 80 * 8 end # 80 bytes

"""
A `BSONObjectId` represents a unique identifier
for a BSON document.

# Example

The following generates a new `BSONObjectId`.

```julia
julia> Mongoc.BSONObjectId()
```

# C API

`BSONObjectId` instances addresses are passed
to libbson/libmongoc API using `Ref(oid)`,
and are owned by the Julia process.

Mirrors C struct `bson_oid_t`:

```c
typedef struct {
   uint8_t bytes[12];
} bson_oid_t;
```
"""
struct BSONObjectId
    bytes::NTuple{12, UInt8}
end

function BSONObjectId()
    oid_ref = Ref{BSONObjectId}()
    bson_oid_init(oid_ref, C_NULL)
    return oid_ref[]
end

function BSONObjectId(oid_string::String)
    if !bson_oid_is_valid(oid_string)
        error("'$oid_string' is not a valid ObjectId.")
    end
    oid_ref = Ref{BSONObjectId}()
    bson_oid_init_from_string(oid_ref, oid_string)
    return oid_ref[]
end

function BSONObjectId(oid::BSONObjectId)
    return BSONObjectId(oid.bytes)
end

Base.copy(oid::BSONObjectId) = BSONObjectId(oid)

"""
`BSONError` is the default `Exception`
for BSON/MongoDB function call errors.

# C API

Mirrors C struct `bson_error_t`.

`BSONError` instances addresses are passed
to libbson/libmongoc API using `Ref(error)`,
and are owned by the Julia process.

```c
typedef struct {
   uint32_t domain;
   uint32_t code;
   char message[504];
} bson_error_t;
```
"""
struct BSONError <: Exception
    domain::UInt32
    code::UInt32
    message::NTuple{504, UInt8}
end

"""
`BSONCode` is a BSON element
with JavaScript source code.

# Example

```julia
julia> bson = Mongoc.BSON("source" => Mongoc.BSONCode("function() = 1"))
BSON("{ "source" : { "\$code" : "function() = 1" } }")
```
"""
struct BSONCode
    code::String
end

"""
A `BSON` represents a document in *Binary JSON* format,
defined at http://bsonspec.org/.

In Julia, you can manipulate a `BSON` instance
just like a `Dict`.

# Example

```julia
bson = Mongoc.BSON()
bson["name"] = "my name"
bson["num"] = 10.0
```

# C API

`BSON` is a wrapper for C struct `bson_t`.
"""
mutable struct BSON <: AbstractDict{String, Any}
    handle::Ptr{Cvoid}

    function BSON(handle::Ptr{Cvoid}; enable_finalizer::Bool=true)
        new_bson = new(handle)
        if enable_finalizer
            finalizer(destroy!, new_bson)
        end
        return new_bson
    end
end

function destroy!(bson::BSON)
    if bson.handle != C_NULL
        bson_destroy(bson.handle)
        bson.handle = C_NULL
    end
    nothing
end

"""
Wrapper for bson_value_t.

See [`Mongoc.get_as_bson_value`](@ref).
"""
mutable struct BSONValue
    handle::Ptr{Cvoid}
end

function Base.deepcopy(bson::BSON) :: BSON
    return BSON(bson_copy(bson.handle))
end

abstract type AbstractBSONReader end

mutable struct BSONReader <: AbstractBSONReader
    handle::Ptr{Cvoid}
    data::Vector{UInt8}

    function BSONReader(handle::Ptr{Cvoid}, data::Vector{UInt8})
        new_reader = new(handle, data)
        finalizer(destroy!, new_reader)
        return new_reader
    end
end

BSONReader(handle::Ptr{Cvoid}) = BSONReader(handle, Vector{UInt8}())

function destroy!(reader::BSONReader)
    if reader.handle != C_NULL
        bson_reader_destroy(reader.handle)
        reader.handle = C_NULL
    end
    nothing
end

const DEFAULT_BSON_WRITER_BUFFER_CAPACITY = 2^10

# buffer_handle must be a valid pointer to a Buffer object.
# A reference to this buffer must be kept from outside this function
# to avoid GC on it.
function unsafe_buffer_realloc(buffer_ptr::Ptr{UInt8}, num_bytes::Csize_t, writer_objref::Ptr{Cvoid})
    local writer::BSONWriter = unsafe_pointer_to_objref(writer_objref)
    @assert buffer_ptr == pointer(writer.buffer) "Is this the same BSONWriter?"

    current_len = length(writer.buffer)
    inc = num_bytes - current_len
    if inc > 0
        append!(writer.buffer, [ UInt8(0) for i in 1:inc ])
    end

    @assert length(writer.buffer) == num_bytes
    writer.buffer_length_ref[] = num_bytes

    return pointer(writer.buffer)
end

mutable struct BSONWriter
    handle::Ptr{Cvoid}
    buffer::Vector{UInt8}
    buffer_handle_ref::Ref{Ptr{UInt8}}
    buffer_length_ref::Ref{Csize_t}

    function BSONWriter(initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY,
                        buffer_offset::Integer=0)

        buffer = zeros(UInt8, initial_buffer_capacity)
        new_writer = new(C_NULL, buffer, Ref(pointer(buffer)), Ref(Csize_t(initial_buffer_capacity)))
        finalizer(destroy!, new_writer)

        realloc_func = @cfunction(unsafe_buffer_realloc, Ptr{UInt8}, (Ptr{UInt8}, Csize_t, Ptr{Cvoid}))

        handle = bson_writer_new(new_writer.buffer_handle_ref,
                                 new_writer.buffer_length_ref,
                                 Csize_t(buffer_offset),
                                 realloc_func,
                                 pointer_from_objref(new_writer))
        if handle == C_NULL
            error("Failed to create a BSONWriter.")
        end
        new_writer.handle = handle

        return new_writer
    end
end

function destroy!(writer::BSONWriter)
    if writer.handle != C_NULL
        bson_writer_destroy(writer.handle)
        writer.handle = C_NULL
    end
    nothing
end

#
# API
#

Base.convert(::Type{String}, oid::BSONObjectId) = bson_oid_to_string(oid)
Base.string(oid::BSONObjectId) = convert(String, oid)
Base.convert(::Type{BSONObjectId}, oid_string::String) = BSONObjectId(oid_string)

Base.convert(::Type{String}, code::BSONCode) = code.code
Base.string(code::BSONCode) = convert(String, code)
Base.convert(::Type{BSONCode}, code_string::String) = BSONCode(code_string)

Base.show(io::IO, oid::BSONObjectId) = print(io, "BSONObjectId(\"", string(oid), "\")")
Base.show(io::IO, bson::BSON) = print(io, "BSON(\"", as_json(bson), "\")")
Base.show(io::IO, code::BSONCode) = print(io::IO, "BSONCode(\"$(code.code)\")")

function Base.show(io::IO, err::BSONError)
    print(io, "BSONError: domain=$(Int(err.domain)), code=$(Int(err.code)), message=")
    for c in err.message
        c_char = Char(c)
        if c_char == '\0'
            break
        else
            print(io, c_char)
        end
    end
end

Base.showerror(io::IO, err::BSONError) = show(io, err)

function get_time(oid::BSONObjectId)
    return Dates.unix2datetime(bson_oid_get_time_t(oid))
end

function BSON()
    handle = bson_new()
    if handle == C_NULL
        error("Failed to create a new BSON document.")
    end
    return BSON(handle)
end

function BSON(json_string::String)
    handle = bson_new_from_json(json_string)
    if handle == C_NULL
        error("Failed parsing JSON to BSON. $json_string.")
    end
    return BSON(handle)
end

function BSON(vector::Vector)
    result = BSON()

    for (i, v) in enumerate(vector)
        result[string(i-1)] = v
    end

    return result
end

function BSON(args::Pair...)
    result = BSON()

    for (k, v) in args
        result[k] = v
    end

    return result
end

BSON(dict::Dict) = BSON(dict...)

macro bson_str(json_string)
    BSON(json_string)
end

"""
    as_json(bson::BSON; canonical::Bool=false) :: String

Converts a `bson` object to a JSON string.

# Example

```julia
julia> document = Mongoc.BSON("{ \"hey\" : 1 }")
BSON("{ "hey" : 1 }")

julia> Mongoc.as_json(document)
"{ \"hey\" : 1 }"

julia> Mongoc.as_json(document, canonical=true)
"{ \"hey\" : { \"\$numberInt\" : \"1\" } }"
```

# C API

* [`bson_as_canonical_extended_json`](http://mongoc.org/libbson/current/bson_as_canonical_extended_json.html)

* [`bson_as_relaxed_extended_json`](http://mongoc.org/libbson/current/bson_as_relaxed_extended_json.html)
"""
function as_json(bson::BSON; canonical::Bool=false) :: String
    local bson_cstring::Cstring

    if canonical
        bson_cstring = bson_as_canonical_extended_json(bson.handle)
    else
        bson_cstring = bson_as_relaxed_extended_json(bson.handle)
    end

    if bson_cstring == C_NULL
        error("Couldn't convert bson to json.")
    end

    result = unsafe_string(bson_cstring)
    bson_free(convert(Ptr{Cvoid}, bson_cstring))

    return result
end

Base.length(document::BSON) = bson_count_keys(document.handle)

#
# Read values from BSON
#

has_field(bson::BSON, key::AbstractString) = bson_has_field(bson.handle, key)
Base.haskey(bson::BSON, key::AbstractString) = has_field(bson, key)

function bson_iter_init(document::BSON) :: Ref{BSONIter}
    iter_ref = Ref{BSONIter}()
    ok = bson_iter_init(iter_ref, document.handle)
    if !ok
        error("Couldn't create iterator for BSON document.")
    end
    return iter_ref
end

abstract type BSONIteratorMode end

struct IterateKeyValuePairs <: BSONIteratorMode
end

struct IterateKeys <: BSONIteratorMode
end

struct IterateValues <: BSONIteratorMode
end

struct BSONIterator{M<:BSONIteratorMode}
    bson_iter_ref::Ref{BSONIter}
    document::BSON

    function BSONIterator(document::BSON, ::Type{M}) where {M<:BSONIteratorMode}
        iter_ref = bson_iter_init(document)
        return new{M}(iter_ref, document)
    end
end

Base.iterate(document::BSON) = iterate(document, BSONIterator(document, IterateKeyValuePairs))
Base.keys(document::BSON) = BSONIterator(document, IterateKeys)
Base.values(document::BSON) = BSONIterator(document, IterateValues)

function Base.iterate(::BSON, itr::BSONIterator{IterateKeyValuePairs})

    if bson_iter_next(itr.bson_iter_ref)
        key = unsafe_string(bson_iter_key(itr.bson_iter_ref))
        value = get_value(itr.bson_iter_ref)

        return key => value, itr
    end

    return nothing
end

function Base.iterate(itr::BSONIterator{IterateKeys}, ::Nothing=nothing)

    if bson_iter_next(itr.bson_iter_ref)
        key = unsafe_string(bson_iter_key(itr.bson_iter_ref))
        return key, nothing
    end

    return nothing
end

function Base.iterate(itr::BSONIterator{IterateValues}, ::Nothing=nothing)

    if bson_iter_next(itr.bson_iter_ref)
        value = get_value(itr.bson_iter_ref)
        return value, nothing
    end

    return nothing
end

"""
    as_dict(document::BSON) :: Dict

Converts a BSON document to a Julia `Dict`.
"""
function as_dict(document::BSON) :: Dict
    result = Dict()
    for (k, v) in document
        result[k] = v
    end
    return result
end

function as_dict(iter_ref::Ref{BSONIter}) :: Dict
    result = Dict()
    while bson_iter_next(iter_ref)
        result[unsafe_string(bson_iter_key(iter_ref))] = get_value(iter_ref)
    end
    return result
end

function get_value(iter_ref::Ref{BSONIter})
    local bson_type::BSONType = bson_iter_type(iter_ref)

    if bson_type == BSON_TYPE_UTF8
        return unsafe_string(bson_iter_utf8(iter_ref))
    elseif bson_type == BSON_TYPE_INT64
        return bson_iter_int64(iter_ref)
    elseif bson_type == BSON_TYPE_INT32
        return bson_iter_int32(iter_ref)
    elseif bson_type == BSON_TYPE_DOUBLE
        return bson_iter_double(iter_ref)
    elseif bson_type == BSON_TYPE_OID
        # converts Ptr{BSONObjectId} to BSONObjectId
        return unsafe_load(bson_iter_oid(iter_ref))
    elseif bson_type == BSON_TYPE_BOOL
        return bson_iter_bool(iter_ref)
    elseif bson_type == BSON_TYPE_DATE_TIME
        millis = Int64(bson_iter_date_time(iter_ref))
        return isodate2datetime(millis)

    elseif bson_type == BSON_TYPE_ARRAY || bson_type == BSON_TYPE_DOCUMENT

        child_iter_ref = Ref{BSONIter}()
        ok = bson_iter_recurse(iter_ref, child_iter_ref)
        if !ok
            error("Couldn't iterate array inside BSON.")
        end

        if bson_type == BSON_TYPE_ARRAY
            result_vector = Vector()
            while bson_iter_next(child_iter_ref)
                push!(result_vector, get_value(child_iter_ref))
            end
            return result_vector
        else
            @assert bson_type == BSON_TYPE_DOCUMENT
            return as_dict(child_iter_ref)
        end
    elseif bson_type == BSON_TYPE_BINARY

        length_ref = Ref{UInt32}()
        buffer_ref = Ref{Ptr{UInt8}}()
        bson_iter_binary(iter_ref, length_ref, buffer_ref)

        result_data = Vector{UInt8}(undef, length_ref[])
        unsafe_copyto!(pointer(result_data), buffer_ref[], length_ref[])

        return result_data

    elseif bson_type == BSON_TYPE_CODE
        return BSONCode(unsafe_string(bson_iter_code(iter_ref)))
    elseif bson_type == BSON_TYPE_NULL
        return nothing
    else
        error("BSON Type not supported: $bson_type.")
    end
end

function Base.getindex(document::BSON, key::AbstractString)
    iter_ref = Ref{BSONIter}()
    ok = bson_iter_init_find(iter_ref, document.handle, key)
    if !ok
        error("Key $key not found.")
    end
    return get_value(iter_ref)
end

"""
    get_as_bson_value(doc, key) :: BSONValue

Returns a value stored in a bson document `doc`
as a `BSONValue`.

See also [Mongoc.BSONValue](@ref).
"""
function get_as_bson_value(document::BSON, key::AbstractString) :: BSONValue
   iter_ref = Ref{BSONIter}()
    ok = bson_iter_init_find(iter_ref, document.handle, key)
    if !ok
        error("Key $key not found.")
    end
    return get_as_bson_value(iter_ref)
end

function get_as_bson_value(iter_ref::Ref{BSONIter})
    # this BSONValue may be valid only during
    # the lifetime of the BSON document
    return BSONValue(bson_iter_value(iter_ref))
end

#
# Write values to BSON
#

function Base.setindex!(document::BSON, value::BSONObjectId, key::AbstractString)
    ok = bson_append_oid(document.handle, key, -1, value)
    if !ok
        error("Couldn't append oid to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Int64, key::AbstractString)
    ok = bson_append_int64(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Int64 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Int32, key::AbstractString)
    ok = bson_append_int32(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Int32 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::AbstractString, key::AbstractString)
    ok = bson_append_utf8(document.handle, key, -1, value, -1)
    if !ok
        error("Couldn't append String to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Bool, key::AbstractString)
    ok = bson_append_bool(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Bool to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Float64, key::AbstractString)
    ok = bson_append_double(document.handle, key, -1, value)
    if !ok
        error("Couldn't append Float64 to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::DateTime, key::AbstractString)
    ok = bson_append_date_time(document.handle, key, -1, datetime2isodate(value))
    if !ok
        error("Couldn't append DateTime to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::BSON, key::AbstractString)
    ok = bson_append_document(document.handle, key, -1, value.handle)
    if !ok
        error("Couldn't append Sub-Document BSON to BSON document.")
    end
    nothing
end

Base.setindex!(document::BSON, value::Dict, key::AbstractString) = setindex!(document, BSON(value), key)

function Base.setindex!(document::BSON, value::Vector{T}, key::AbstractString) where T
    sub_document = BSON(value)
    ok = bson_append_array(document.handle, key, -1, sub_document.handle)
    if !ok
        error("Couldn't append array to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::BSONCode, key::AbstractString)
    ok = bson_append_code(document.handle, key, -1, value.code)
    if !ok
        error("Couldn't append String to BSON document.")
    end
    nothing
end

function Base.setindex!(document::BSON, value::Date, key::AbstractString)
    error("BSON format does not support `Date` type. Use `DateTime` instead.")
end

function Base.setindex!(document::BSON, value::Vector{UInt8}, key::AbstractString)
  ok = bson_append_binary(document.handle, key, -1, BSON_SUBTYPE_BINARY, value, UInt32(length(value)))
  if !ok
      error("Couldn't append array to BSON document.")
  end
  nothing
end

function Base.setindex!(document::BSON, ::Nothing, key::AbstractString)
    ok = bson_append_null(document.handle, key, -1)
    if !ok
        error("Couldn't append missing value to BSON document.")
    end
    nothing
end

#
# Write/Read BSON to/from IO
#

function _get_number_of_bytes_written_to_buffer(buffer::Vector{UInt8}) :: Int64
    isempty(buffer) && return 0

    # the first 4 bytes in data is a size int32
    doc_size = reinterpret(Int32, buffer[1:4])[1]
    local total_size::Int64 = 0
    while doc_size != 0
        total_size += doc_size

        if length(buffer) < total_size + 4
            break
        end

        doc_size = reinterpret(Int32, buffer[total_size+1:total_size+4])[1]
    end
    return total_size
end

function _get_number_of_bytes_written_to_buffer(writer::BSONWriter)
    _get_number_of_bytes_written_to_buffer(writer.buffer)
end

function bson_writer(f::Function, io::IO;
                     initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY)

    writer = BSONWriter(initial_buffer_capacity)

    try
        f(writer)

        @assert bson_writer_get_length(writer.handle) == _get_number_of_bytes_written_to_buffer(writer)

        for byte_index in 1:bson_writer_get_length(writer.handle)
            write(io, writer.buffer[byte_index])
        end

        flush(io)
    finally
        destroy!(writer)
    end

    nothing
end

function write_bson(f::Function, writer::BSONWriter)
    bson_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
    ok = bson_writer_begin(writer.handle, bson_handle_ref)
    if !ok
        error("Failed to write bson document to IO: there was not enough space in buffer. Increase the buffer initial capacity.")
    end

    bson = BSON(bson_handle_ref[], enable_finalizer=false)
    f(bson)
    bson_writer_end(writer.handle)
end

"""
    write_bson(io::IO, bson::BSON;
        initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY)

Writes a single BSON document to `io` in binary format.
"""
function write_bson(io::IO, bson::BSON;
                    initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY)

    bson_writer(io, initial_buffer_capacity=initial_buffer_capacity) do writer
        write_bson(writer) do dest
            bson_copy_to_excluding_noinit(bson.handle, dest.handle)
        end
    end

    nothing
end

struct BufferedBSON
    bson_data::Vector{UInt8}
end

function BufferedBSON(bson::BSON)
    io = IOBuffer()
    write_bson(io, bson)
    return BufferedBSON(take!(io))
end

BSON(buff::BufferedBSON) :: BSON = read_bson(buff.bson_data)[1]

function Serialization.serialize(s::AbstractSerializer, bson::BSON)
    Serialization.serialize_type(s, BSON)
    Serialization.serialize(s.io, BufferedBSON(bson))
end

function Serialization.deserialize(s::AbstractSerializer, ::Type{BSON})
    BSON(Serialization.deserialize(s.io))
end

Base.write(io::IO, bson::BSON) = serialize(io, bson)
Base.read(io::IO, ::Type{BSON}) = deserialize(io)::BSON

"""
    write_bson(io::IO, bson_list::Vector{BSON};
        initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY)

Writes a vector of BSON documents to `io` in binary format.

# Example

```julia
list = Vector{Mongoc.BSON}()

let
    src = Mongoc.BSON()
    src["id"] = 1
    src["name"] = "1st"
    push!(list, src)
end

let
    src = Mongoc.BSON()
    src["id"] = 2
    src["name"] = "2nd"
    push!(list, src)
end

open("documents.bson", "w") do io
    Mongoc.write_bson(io, list)
end
```
"""
function write_bson(io::IO, bson_list::Vector{BSON};
                    initial_buffer_capacity::Integer=DEFAULT_BSON_WRITER_BUFFER_CAPACITY)

    bson_writer(io, initial_buffer_capacity=initial_buffer_capacity) do writer
        for src_bson in bson_list
            write_bson(writer) do dest
                bson_copy_to_excluding_noinit(src_bson.handle, dest.handle)
            end
        end
    end

    nothing
end

"""
    read_bson(io::IO) :: Vector{BSON}

Reads all BSON documents from `io`.
This method will continue to read from `io` until
it reaches eof.
"""
read_bson(io::IO) :: Vector{BSON} = read_bson(read(io))

"""
    read_bson(data::Vector{UInt8}) :: Vector{BSON}

Parses a vector of bytes to a vector of BSON documents.
Useful when reading BSON as binary from a stream.
"""
function read_bson(data::Vector{UInt8}) :: Vector{BSON}
    reader = BSONReader(data)

    try
        return read_bson(reader)
    finally
        destroy!(reader)
    end
end

function BSONReader(data::Vector{UInt8})
    if isempty(data)
        return result
    end

    data_shrinked = data[1:_get_number_of_bytes_written_to_buffer(data)]

    reader_handle = bson_reader_new_from_data(pointer(data_shrinked), length(data_shrinked))
    if reader_handle == C_NULL
        error("Failed to create a BSONReader.")
    end

    return BSONReader(reader_handle, data_shrinked)
end

"""
    read_bson(filepath::AbstractString) :: Vector{BSON}

Reads all BSON documents from a file located at `filepath`.

This will open a `Mongoc.BSONReader` pointing to the file
and will parse file contents to BSON documents.
"""
function read_bson(filepath::AbstractString) :: Vector{BSON}
    reader = BSONReader(filepath)

    try
        return read_bson(reader)
    finally
        destroy!(reader)
    end
end

function BSONReader(filepath::AbstractString)
    @assert isfile(filepath) "$filepath not found."
    err_ref = Ref{BSONError}()
    reader_handle = bson_reader_new_from_file(filepath, err_ref)
    if reader_handle == C_NULL
        throw(err_ref[])
    end
    return BSONReader(reader_handle)
end

"""
    read_next_bson(reader::BSONReader) :: Union{Nothing, BSON}

Reads the next BSON document available in the stream pointed by `reader`.
Returns `nothing` if reached the end of the stream.
"""
function read_next_bson(reader::BSONReader) :: Union{Nothing, BSON}
    reached_eof_ref = Ref(false)
    bson_handle = bson_reader_read(reader.handle, reached_eof_ref)

    if bson_handle == C_NULL
        if !reached_eof_ref[]
            @warn("Finished reading BSON from stream, but didn't reach stream's eof.")
        end
        return nothing
    end

    bson_copy_handle = bson_copy(bson_handle) # creates a copy with its own lifetime
    return BSON(bson_copy_handle)
end

mutable struct BSONJSONReader <: AbstractBSONReader
    handle::Ptr{Cvoid}

    function BSONJSONReader(handle::Ptr{Cvoid})
        new_reader = new(handle)
        finalizer(destroy!, new_reader)
        return new_reader
    end
end

function destroy!(reader::BSONJSONReader)
    if reader.handle != C_NULL
        bson_json_reader_destroy(reader.handle)
        reader.handle = C_NULL
    end
    nothing
end

function BSONJSONReader(filepath::AbstractString)
    @assert isfile(filepath) "$filepath not found."
    err_ref = Ref{BSONError}()
    reader_handle = bson_json_reader_new_from_file(filepath, err_ref)
    if reader_handle == C_NULL
        throw(err_ref[])
    end
    return BSONJSONReader(reader_handle)
end

function read_next_bson(reader::BSONJSONReader, buffer::BSON=BSON()) :: Union{Nothing, BSON}
    err_ref = Ref{BSONError}()
    ok = bson_json_reader_read(reader.handle, buffer.handle, err_ref)

    if ok == 1 # successful and data was read
        bson_copy_handle = bson_copy(buffer.handle)
        return BSON(bson_copy_handle)
    elseif ok == 0 # successful and no data was read
        return nothing
    elseif ok == -1 # there was an error
        throw(err_ref[])
    else # unexpected
        error("Unexpected result from bson_json_reader_read: $ok.")
    end
end

function Base.iterate(reader::AbstractBSONReader, state=nothing)
    next_bson = read_next_bson(reader)
    if next_bson == nothing
        return nothing
    else
        return (next_bson, nothing)
    end
end

Base.collect(reader::AbstractBSONReader) = read_bson(reader)

"""
    read_bson(reader::BSONReader) :: Vector{BSON}

Reads all BSON documents from a `reader`.
"""
function read_bson(reader::AbstractBSONReader) :: Vector{BSON}
    result = Vector{BSON}()
    for bson in reader
        push!(result, bson)
    end
    return result
end

"""
    read_bson_from_json(filepath::AbstractString) :: Vector{BSON}

Reads a JSON file into a list of BSON documents.
The file should contain a sequence of valid JSON documents.

# Example of a valid JSON file

```json
{ "num" : 1, "str" : "two" }
{ "num" : 2, "str" : "three" }
```
"""
function read_bson_from_json(filepath::AbstractString) :: Vector{BSON}
    reader = BSONJSONReader(filepath)
    try
        return read_bson(reader)
    finally
        destroy!(reader)
    end
end

ext: jl


text:
#=
Test the diffusion function generate the matrix which is square and with same length with the dimension of ODE.
=#

m_diff = LeastAction.DiffusionMatrix(length(sp.u0))

@testset "Size of m_diff should be square and equal to the length of input vetor" for i in size(m_diff)
    @test i == length(sp.u0)
end;


text:
# Base Interface
import Base:
    size, getindex, *
import Base.LinAlg:
    A_mul_B!, At_mul_B, A_ldiv_B!

# TODO: make sure this implementation is stable

"""
    DenseBlockAngular{Tv<:Real, Ti<:Integer}

# Attributes
-`m::Ti`: Number of linking constraints
-`n::Ti`: Total number of columns
-`R::Ti`: Number of blocks
-`colptr::Vector{Ti}`: Index of first column of each block. The columns of
    block `r` are columns `colptr[k] ... colptr[k+1]-1`.
-`cols::Matrix{Tv}`: A `m x n` matrix containing the columns.
"""
struct DenseBlockAngular{Tv<:Real, Ti<:Integer} <: AbstractMatrix{Tv}
    m::Ti
    n::Ti
    R::Ti

    colptr::Vector{Ti}
    cols::Matrix{Tv}

    DenseBlockAngular(m::Ti, n::Ti, R::Ti, colptr::Vector{Ti}, cols::Matrix{Tv}) where {Tv<:Real, Ti<:Integer} =
        new{Tv, Ti}(m, n, R, colptr, cols)
end

function DenseBlockAngular(col_list::Vector{Matrix{T}}) where T<:Real

    R = size(col_list, 1)
    if R == 0
        return DenseBlockAngular(0, 0, 0, Vector{Int}(0), Matrix{T}(0, 0))
    end
    m = size(col_list[1], 1)

    colptr = Vector{Int}(R+1)
    colptr[1] = 1

    # Dimension check
    for r in 1:R
        m == size(col_list[r], 1) || throw(DimensionMismatch("Block 1 has dimensions $(size(col_list[1])) but block $(r) has dimensions $(size(col_list[r]))"))
        colptr[r+1] = colptr[r] + size(col_list[r], 2)
    end

    cols = hcat(col_list...)

    return DenseBlockAngular(m, colptr[end]-1, R, colptr, cols)
end

size(M::DenseBlockAngular) = (M.m+M.R, M.n)
function getindex(M::DenseBlockAngular{Tv, Ti}, i::Integer, j::Integer) where {Tv<:Real, Ti<:Integer}
    if !(1 <= i <= (M.m+M.R) && 1 <= j <= (M.n)); throw(BoundsError()); end

    if i > M.R
        return M.cols[i - M.R, j]
    else
        # find if column j belongs to block i
        return (M.colptr[i] <= j < M.colptr[i+1]) ? oneunit(Tv) : zero(Tv)
    end

    # never reached, for type stability
    return zero(Tv)
end

# Matrix-vector Multiplication
function *(A::DenseBlockAngular{Tv, Ti}, x::AbstractVector{Tv}) where{Tv<:Real, Ti<:Integer}

    n = size(x, 1)
    n == A.n || throw(DimensionMismatch("A has dimensions $(size(A)) but x has dimension $(size(x))"))

    y = zeros(Tv, size(A, 1))
    r = 1
    for i in 1:n
        if i == A.colptr[r+1]
            r += 1
        end
        y[r] += x[i]
    end

    Base.BLAS.gemv!('N', 1.0, A.cols, x, 0.0, view(y, (A.R+1):size(A, 1)))

    return y
end

function A_mul_B!(y::AbstractVector{Tv}, A::DenseBlockAngular{Tv, Ti}, x::AbstractVector{Tv}) where{Tv<:Real, Ti<:Integer}

    m, n = size(A)
    n == size(x, 1) || throw(DimensionMismatch("A has dimensions $(size(A)) but x has dimension $(size(x))"))
    m == size(y, 1) || throw(DimensionMismatch("A has dimensions $(size(A)) but y has dimension $(size(y))"))

    @inbounds for r in 1:A.R
        @views y[r] = sum(x[A.colptr[r]:(A.colptr[r+1]-1)])
    end

    Base.BLAS.gemv!('N', 1.0, A.cols, x, 0.0, view(y, (A.R+1):size(A, 1)))

    return y
end

function At_mul_B(A::DenseBlockAngular{Tv, Ti}, y::AbstractVector{Tv}) where{Tv<:Real, Ti<:Integer}

    m, n = size(A)
    m == size(y, 1) || throw(DimensionMismatch("A has dimension $(size(A)) but y has dimension $(size(y))"))

    x = ones(n)
    @inbounds for r in 1:A.R
        for i in A.colptr[r]:(A.colptr[r+1]-1)
            x[i] = y[r]
        end
    end
    @views Base.BLAS.gemv!('T', 1.0, A.cols, y[(A.R+1):end], 1.0, x)

    return x
end

"""
    FactorBlockAngular

# Attributes
-`m::Integer`: number of linking constraints
-`n::Integer`: total number of columns
-`R::Integer`: number of blocks
-`colptr::Vector{Int}`: index of first column of each block. The columns of
    block `r` are columns `colptr[k] ... colptr[k+1]-1`.
-`d::Vector{T}`: `R`-dimensional vector containing the diagonal coefficients
    of the implicit Cholesky factor.
-`η::Matrix{T}`: An `m x R` matrix containing the lower blocks of the
    implicit Cholesky factor.
-`Fc::Base.LinAlg.Cholesky{T, Matrix{T}}`: Cholesky factor of the dense
    Schur complement.
"""
mutable struct FactorBlockAngular{Tv<:Real, Ti<:Integer} <: Factorization{Tv}
    m::Ti
    n::Ti
    R::Ti

    colptr::Vector{Ti}

    d::Vector{Tv}
    η::Matrix{Tv}

    Fc::Base.LinAlg.Cholesky{Tv, Matrix{Tv}}

    FactorBlockAngular(
        m::Ti, n::Ti, R::Ti,
        colptr::Vector{Ti}, d::Vector{Tv}, η::Matrix{Tv},
        Fc::Base.LinAlg.Cholesky{Tv, Matrix{Tv}}
    ) where{Tv<:Real, Ti<:Integer} = new{Tv, Ti}(m, n, R, colptr, d, η, Fc)
end

size(F::FactorBlockAngular) = (F.m+F.R, F.n)

"""
    A_ldiv_B!

"""
function A_ldiv_B!(F::FactorBlockAngular{Tv, Ti}, b::AbstractVector{Tv}) where{Tv<:Real, Ti<:Integer}

    # Dimension check
    (F.m + F.R) == size(b, 1) || throw(DimensionMismatch("F has dimensions $(size(F)) but b has dimensions $(size(b))"))

    @views b[1:F.R] .*= F.d

    # right-hand side for global solve (b0)
    @views Base.BLAS.gemv!('N', -1.0, F.η, b[1:F.R], 1.0, b[(F.R+1):end])

    # global solve
    @views Base.LinAlg.A_ldiv_B!(F.Fc, b[(F.R+1):end])

    # local backward pass
    @views b[1:F.R] ./= F.d
    @views Base.BLAS.gemv!('T', -1.0, F.η, b[(F.R+1):end], 1.0, b[1:F.R])
    @views b[1:F.R] .*= F.d

    return b
end

function A_ldiv_B!(y::AbstractVector{Tv}, F::FactorBlockAngular{Tv, Ti}, b::AbstractVector{Tv}) where{Tv<:Real, Ti<:Integer}

    # Dimension check
    (F.m + F.R) == size(b, 1) || throw(DimensionMismatch("F has dimensions $(size(F)) but b has dimensions $(size(b))"))
    (F.m + F.R) == size(y, 1) || throw(DimensionMismatch("F has dimensions $(size(F)) but y has dimensions $(size(y))"))

    # over-write y with b
    copy!(y, b)
    @views y[1:F.R] .*= F.d

    # right-hand side for global solve
    @views Base.BLAS.gemv!('N', -1.0, F.η, y[1:F.R], 0.0, y[(F.R+1):end])
    @views Base.BLAS.axpy!(1.0, b[(F.R+1):end], y[(F.R+1):end])

    # global solve
    @views Base.LinAlg.A_ldiv_B!(F.Fc, y[(F.R+1):end])

    # local backward pass
    @views Base.BLAS.gemv!('T', -1.0, F.η, y[(F.R+1):end], 0.0, y[1:F.R])
    @views Base.BLAS.axpy!(1.0, b[1:F.R], y[1:F.R])
    @views y[1:F.R] .*= F.d

    return y
end

"""
    CpBDBt!(C, B, d)

Compute C += B*D*B'
    where B' is the transpose of B, and D = diag(d1, ..., dn)
    and B is mxn dense, C is mxm (dense), d is nx1 (dense).

C is modified in-place (no memory allocated),
only its upper-triangular part is modified.

# Arguments
-`C::StridedMatrix{T}`: An `m`-by-`m` dense matrix, modified in place
-`B::StridedMatrix{T}`: An `m`-by-`k` dense matrix
-`d::StridedVector{T}`: A vector of length `k`
"""
function CpBDBt!(C::StridedMatrix{T}, B::StridedMatrix{T}, d::StridedVector{T}) where {T<:Real}
    # TODO: allow user to specify which triangular part should be update
    # TODO: add parameter for computing C = α C + B*D*B'

    # Dimension checks
    (m, k) = size(B)
    (m, m) == size(C) || throw(DimensionMismatch("C has dimensions $(size(C)) but B has dimensions $(size(B))"))
    k == size(d, 1) || throw(DimensionMismatch("d has dimensions $(size(d)) but B has dimensions $(size(B))"))

    # compute symmetric rank-k update
    # only the upper-triangular part of C is modified
    temp = zero(T)
    for j=1:m
        for l = 1:k
            @inbounds temp = d[l] * B[j, l]
            if temp == zero(T)
                # skip loop if temp is zero
                continue
            end
            for i=1:j
                @inbounds C[i, j] = C[i, j] + temp * B[i, l]
            end
        end
    end

    return C
end

# TODO: blocked version
# Compare to BLAS `syrk!`, there is a performance factor of ~2-3 to be gained

function cholesky!(
    A::DenseBlockAngular{Tv, Ti},
    θ::AbstractVector{Tv},
    F::FactorBlockAngular{Tv, Ti}
    ) where{Tv<:Real, Ti<:Integer}

    # Schur complement
    C = zeros(Tv, A.m, A.m)

    for r in 1:A.R
        θ_ = view(θ, A.colptr[r]:(A.colptr[r+1]-1))
        A_ = view(A.cols, :, A.colptr[r]:(A.colptr[r+1]-1))
        η_ = view(F.η, :, r)

        # diagonal elements
        F.d[r] = oneunit(Tv) / sum(θ_)

        # lower factors
        Base.BLAS.gemv!('N', 1.0, A_, θ_, 0.0, η_)

        # Schur complement
        CpBDBt!(C, A_, θ_)
        Base.BLAS.syr!('U', -F.d[r], η_, C)
    end

    F.Fc = cholfact(Symmetric(C))

    return F
end

function cholesky(
    A::DenseBlockAngular{Tv, Ti},
    θ::AbstractVector{Tv}
    ) where{Tv<:Real, Ti<:Integer}

    # Schur complement
    C = zeros(Tv, A.m, A.m)
    d = zeros(Tv, A.R)
    η = zeros(Tv, A.m, A.R)

    for r in 1:A.R
        θ_ = view(θ, A.colptr[r]:(A.colptr[r+1]-1))
        A_ = view(A.cols, :, A.colptr[r]:(A.colptr[r+1]-1))
        η_ = view(η, :, r)

        # diagonal elements
        d[r] = oneunit(Tv) / sum(θ_)

        # lower factors
        Base.BLAS.gemv!('N', 1.0, A_, θ_, 0.0, η_)

        # Schur complement
        CpBDBt!(C, A_, θ_)
        Base.BLAS.syr!('U', -d[r], η_, C)
    end

    # compute Cholesky factor of C
    Fc = cholfact(Symmetric(C))

    F = FactorBlockAngular(A.m, A.n, A.R, A.colptr, d,η, Fc)

    return F
end


text:
# ---
# title: 357. Count Numbers with Unique Digits
# id: problem357
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/count-numbers-with-unique-digits/description/>
# hidden: true
# ---
#
# Given a **non-negative** integer n, count all numbers with unique digits, x,
# where 0 ≤ x < 10n.
#
# **Example:**
#
#
#
#     Input: 2
#     Output: 91
#     Explanation: The answer should be the total numbers in the range of 0 ≤ x < 100,
#                  excluding 11,22,33,44,55,66,77,88,99
#
#
#
#
# **Constraints:**
#
#   * `0 <= n <= 8`
#
#
## @lc code=start
using LeetCode

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


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

"""
A data source is supposed to be imported, exported or both.
Subtypes representing importable data source must implement method `import_data`.
Subtypes representing exportable data source must implement method `export_data`.
The method `pathof` should return the path to the database file or directory.
"""
abstract type AbstractDataSource end
abstract type AbstractDataFile <: AbstractDataSource end

const DataSourceDesignator = Union{AbstractDataSource,Symbol}

"""
    pathof(source::AbstractDataFile)

Return the path to the file or directory of `source`. Defaults to `source.path`.
"""
function pathof(source::AbstractDataFile)
    hasproperty(source, :path) ? source.path : nothing
end

"""
    import_data(source; kwargs...)

Import data from `source` and converts it either into a table or a named collection
thereof. `source` can be a `Symbol` which will be searched in `DATA_SOURCES`.
"""
function import_data end

"""
    import_data(sources::AbstractVector, data; kwargs...)

Import each source in `sources` and return a `Dict`:
```
Dict(source => import_data(source, data; kwargs...) for source ∈ sources)
```
"""
function import_data(
    sources::AbstractVector = option(:sources_to_update);
    cache_data_type = option(:cache_data_type),
    force_update::Bool = option(:force_update),
    kwargs...,
)
    @argcheck !isempty(sources)
    @debug "Importing data." sources
    source_keys = make_symbols(sources)
    try
        if !force_update && cache_data_type != nothing
            data = import_data(
                cache_data_type;
                cache_data_type = nothing,
                force_update = true,
                kwargs...,
            )
            if !(data isa AbstractDict)
                return data
            elseif source_keys ⊆ keys(data)
                return subdict(data, source_keys)
            end
        end
    catch
    end

    result = DataDict()
    for (source, sourcesym) ∈ zip(sources, source_keys)
        data =
            import_data(source; cache_data_type = nothing, force_update = true, kwargs...)
        result[sourcesym] = data
    end
    if cache_data_type != nothing
        export_data(cache_data_type, result; kwargs...)
    end
    @debug "Data imported." sources _debuginfo(result)
    result
end

function import_data(
    sourcesym::Symbol;
    path = nothing,
    cache_data_type = option(:cache_data_type),
    force_update::Bool = option(:force_update),
    kwargs...,
)
    try
        if !force_update && cache_data_type != nothing
            data = import_data(
                cache_data_type;
                cache_data_type = nothing,
                force_update = false,
                kwargs...,
            )
            if !(data isa AbstractDict)
                return data
            elseif haskey(data, sourcesym)
                return data[sourcesym]
            end
        end
    catch
    end

    source = DATA_SOURCES[sourcesym](nothing; path = path, kwargs...)
    data = import_data(source; cache_data_type = nothing, force_update = true, kwargs...)
    if cache_data_type != nothing
        export_data(cache_data_type, Dict(sourcesym => data); kwargs...)
    end
    data
end

"""
    export_data(source, data; kwargs...)

Export `data` to `source`; `data` must be a table or a named collection thereof.
If `source` can be a `Symbol`, search it in `DATA_SOURCES`.
"""
function export_data end

"""
    export_data(sources::AbstractVector, data; kwargs...)

Export each source in `sources` and return an `Dict`:
```
Dict(source => export_data(source, data; kwargs...) for source ∈ sources)
```
"""
function export_data(
    sources::AbstractVector,
    data;
    source_keys::AbstractVector{Symbol} = make_symbols(sources),
    kwargs...,
)
    Dict(
        sourcesym => export_data(source; kwargs...)
        for (source, sourcesym) ∈ zip(sources, source_keys)
    )
end

function export_data(sourcesym::Symbol, data; output = nothing, kwargs...)
    source = DATA_SOURCES[sourcesym](nothing; path = output, kwargs...)
    export_data(source, data; kwargs...)
end

function export_data(data; cache_data_type = option(:cache_data_type), kwargs...)
    export_data(cache_data_type, data; kwargs...)
end

"""
    cache_path(
        ext::Union{Nothing, AbstractString} = databaseExtension;
        directory = cache_directory,
        sub = cache_filename,
        nsamples::Integer = cache_nsamples,
        kwargs...,
    )

Returns the default cache file path with the given `ext`ension, or no extension if `ext`
is `nothing`.
"""
function cache_path(
    ext::Union{Nothing,AbstractString} = data_extension(option(:cache_data_type));
    directory = option(:cache_directory),
    update_period::Period = option(:update_period),
    subpath = option(:cache_filename),
    nsamples::Integer = option(:cache_nsamples),
    kwargs...,
)
    @argcheck nsamples >= 1
    (period_directory, lastdir) = _cache_subpath(Path(directory), update_period)
    subpath = isnothing(ext) ? subpath : subpath * "." * ext
    result = joinpath(period_directory, lastdir, Path(subpath))
    curdir = Path(dirname(result))
    if !exists(curdir)
        mkdir(curdir, recursive = true)
        gc_cache_directory(period_directory, nsamples)
    end
    result
end

function _cache_subpath(directory::AbstractPath, update_period::Period)
    try
        update_period = update_period < Second(1) ? Second(1) : update_period
    catch
    end
    n = Dates.value(update_period)
    timeformat = _dateformat_aux(update_period)
    period_directory = joinpath(directory, Path("Update every " * string(update_period)))
    (perfloor, perceil) = Dates.floorceil(now(), update_period)
    if n == 1
        lastdir = Path(Dates.format(perfloor, timeformat))
    else
        lastdir = Path(
            Dates.format(perfloor, timeformat) * " to " * Dates.format(perceil, timeformat),
        )
    end
    (period_directory, lastdir)
end

function cache_path(sourcetype::Type{T}; kwargs...) where {T<:AbstractDataFile}
    cache_path(data_extension(T); kwargs...)
end
function cache_path(sourcesym::Symbol; kwargs...)
    source = DATA_SOURCES[sourcesym]
    cache_path(data_extension(source); kwargs...)
end

_dateformat_aux(::Year) = dateformat"yyyy"
_dateformat_aux(per::Month) =
    Dates.value(per) >= 12 ? dateformat"yyyy" : dateformat"yyyy-mm"
_dateformat_aux(per::Week) =
    Dates.value(per) >= 5 ? dateformat"yyyy-mm" : dateformat"yyyy-mm-dd"
_dateformat_aux(per::Day) =
    Dates.value(per) >= 31 ? dateformat"yyyy-mm" : dateformat"yyyy-mm-dd"
_dateformat_aux(per::Hour) =
    Dates.value(per) >= 24 ? dateformat"yyyy-mm-dd" : dateformat"yyyy-mm-dd HHh"
_dateformat_aux(per::Minute) =
    Dates.value(per) >= 60 ? dateformat"yyyy-mm-dd HHh" : dateformat"yyyy-mm-dd HHhMM\m"
_dateformat_aux(per::Second) =
    Dates.value(per) >= 60 ? dateformat"yyyy-mm-dd HHhMM\m" :
    dateformat"yyyy-mm-dd HHhMM\mSS\s"

_dateformat_aux(period::Period) = dateformat"yyyy-mm-dd HHhMM\mSS\s"

"""
    gc_cache_directory(directory::AbstractPath, nsamples::Integer)

Garbage collect `directory`. Each subdirectory in cache `directory` is named by
date and/or time it was saved. This function keeps only the oldest `nsamples`
subdirectories.
"""
function gc_cache_directory(directory::AbstractPath, nsamples::Integer)
    @argcheck nsamples >= 1
    dirlist = readpath(directory)
    ndel = length(dirlist) - nsamples
    if ndel > 0
        todelete = sort!(dirlist; by = basename)[1:ndel]
        foreach(d -> rm(d; recursive = true), todelete)
    end
end

tableKeys = [:country, :state, :date]
tableColumns =
    [:country, :state, :latitude, :longitude, :date, :confirmed, :recovered, :deaths]

"""
Auxiliar structure to handle downloads.
"""
struct DownloadSource <: AbstractDataSource
    url::AbstractString
    function DownloadSource(url; kwargs...)
        new(string(url))
    end
end

function import_data(data::DownloadSource; kwargs...)
    url = data.url
    ext = split(url, '.')[end]
    ext = all(isletter, ext) ? ext : nothing
    subpath = joinpath(p"_DOWNLOADS_/", Path(string(hash(url))))
    path = cache_path(ext; subpath = subpath, kwargs...)
    @debug "Downloading data." url path
    download(url, path)
end

# CSV file and directory of files
"""
Path to a CSV file or directory
"""
struct CsvPath <: AbstractDataFile
    path::AbstractPath
end
function CsvPath(::Nothing; path = nothing, kwargs...)
    if path == nothing
        path = cache_path("csv"; kwargs...)
    end
    CsvPath(Path(path))
end

DATA_SOURCES[:csv] = CsvPath
data_extension(::Type{CsvPath}) = "csv"

function import_data(source::CsvPath; simple::Bool = false, kwargs...)
    path = pathof(source)
    @argcheck exists(path)
    if isfile(path)
        @debug "Importing DataFrame from CSV file." path
        data = csv_read(path; copycols = true)
    else
        @debug "Importing Dict from CSV directory." path
        data = DataDict(
            # Parameter `simple` don't need to be passed on
            Symbol(filename(elt)) => import_data(CsvPath(elt); kwargs...)
            for elt ∈ readpath(path)
        )
    end
    if simple
        data = simplify(data; kwargs...)
    end
    data
end

function export_data(
    source::CsvPath,
    data::AbstractDataFrame;
    pretty::Bool = false,
    kwargs...,
)
    destiny = pathof(source)
    @debug "Exporting DataFrame to CSV file." destiny _debuginfo(data)
    if pretty
        data = prettify(data; kwargs...)
    end
    exists(destiny) && rm(destiny; recursive = true)
    csv_write(destiny, data; transform = (col, val) -> something(val, missing))
    destiny
end

function export_data(source::CsvPath, data::AbstractDict; pretty::Bool = false, kwargs...)
    destiny = pathof(source)
    @debug "Exporting Dict to CSV directory." destiny _debuginfo(data)
    if pretty
        data = prettify(data; kwargs...)
    end
    !exists(destiny) && mkdir(destiny; recursive = true)
    for (fname, fdata) ∈ data
        fname = replace(string(fname), '/' => '_')
        fpath = joinpath(destiny, Path(fname * ".csv"))
        export_data(CsvPath(fpath), fdata; pretty = pretty, kwargs...)
    end
    destiny
end

struct OdsPath <: AbstractDataFile
    path::AbstractPath
end
function OdsPath(::Nothing; path = nothing, kwargs...)
    if path == nothing
        path = cache_path("ods"; kwargs...)
    end
    OdsPath(Path(path))
end

DATA_SOURCES[:ods] = OdsPath
data_extension(::Type{OdsPath}) = "ods"

function import_data(source::OdsPath; simple::Bool = false, kwargs...)
    path = pathof(source)
    @argcheck exists(path)
    if isfile(path)
        data = DataDict(Symbol(k) => v for (k, v) ∈ ods_readall(path))
    else
        data = DataDict(
            Symbol(filename(elt)) => import_data(OdsPath(elt); kwargs...)
            for elt ∈ readpath(path)
        )
        dfdict = get(data, :_DEFAULT_, nothing)
        if dfdict != nothing
            delete!(dfdict::Dict, :_DEFAULT_)
            merge!(data, dfdict)
        end
        data
    end
    if simple
        data = simplify(data; kwargs...)
    end
end

function export_data(
    source::OdsPath,
    data::AbstractDataFrame;
    pretty::Bool = false,
    table_name = option(:cache_table_name),
    kwargs...,
)
    if pretty
        data = prettify(data; kwargs...)
    end
    export_data(
        source,
        DataFrameDict(Symbol(table_name) => data);
        pretty = pretty,
        kwargs...,
    )
end

function export_data(source::OdsPath, data::DataFrameDict; pretty::Bool = false, kwargs...)
    destiny = pathof(source)
    rm(destiny; recursive = true, force = true)
    sheetdict = SortedDict{Tuple,DataFrame}()
    for (tname, tbl) ∈ data
        tname = string(tname)
        if pretty
            tname = prettify(tname; kwargs...)
            tbl = prettify(tbl; kwargs...)
        end
        push!(sheetdict, (tname, 1, 1) => tbl)
    end
    ods_write(destiny, sheetdict)
    destiny
end

function export_data(
    source::OdsPath,
    data::AbstractDataDict;
    pretty::Bool = false,
    kwargs...,
)
    if pretty
        data = prettify(data; kwargs...)
    end
    destiny = pathof(source)
    #rm(destiny; recursive = true, force = true)
    subdata = filter(p -> p[2] isa AbstractDict, data)
    dataframes = filter(p -> !(p[2] isa AbstractDict), data)
    dataframes = Dict{Symbol,DataFrame}(k => DataFrame(v) for (k, v) ∈ dataframes)
    if isempty(subdata)
        dfsource = source
    else
        !exists(destiny) && mkdir(destiny; recursive = true)
        dfsource = OdsPath(join(destiny, p"_DEFAULT_.ods"))

        for (fname, fdata) ∈ subdata
            fname = replace(string(fname), '/' => '_')
            fpath = joinpath(destiny, Path(fname * ".ods"))
            export_data(OdsPath(fpath), fdata; pretty = pretty, kwargs...)
        end
    end
    export_data(dfsource::OdsPath, dataframes; kwargs...)
    destiny
end


text:
using LaTeXStrings

export tex

@api tex => (TeX <: Tile) begin
    doc("Create TeX/LaTeX tiles from `LaTeXString` object.")
    arg(source::AbstractString, doc="The source TeX object.")
    kwarg(
        block::Bool=false,
        doc="""If set to true, the resulting tile will be a block. It is inline
             by default"""
    )
end

render(l::TeX, state) = Elem("ka-tex", attributes=@d(:source=>replace(l.source, r"(^\$)|(\$$)", ""), :block=>boolattr(l.block)))


text:
@testset "Arithmethic" begin
    @testset "basic arithmetic: $T" for T in [Arb, Acb]
        @test T(3) + 1 == 4
        @test T(3) + 1 isa T
        @test T(2) + T(3) == 5
        @test T(2) + T(3) isa T
        @test T(3) - 1 == 2
        @test T(3) - 1 isa T
        @test T(2) - T(3) == -1
        @test T(2) - T(3) isa T
        @test T(3) * 1 == 3
        @test T(3) * 1 isa T
        @test Bool(Arblib.contains(T(3) / 1, T(3)))
        @test T(3) / 1 isa T
        @test Bool(Arblib.contains(T(3) / T(1), T(3)))
        @test T(3) / T(1) isa T
    end

    @testset "promote: $T" for (T, TMat) in [(Arb, ArbMatrix), (Acb, AcbMatrix)]
        A = TMat(3, 3)
        # promotion of TRef to T
        a = A[1, 1]
        promote(a, 1) isa Tuple{T,T}
        promote(a, T(1)) isa Tuple{T,T}
        @test a + 3 == 3
        @test a + 3 isa T
    end

    @testset "real/imag" begin
        x, y = Arb(rand()), Arb(rand())
        z = Acb(x, y)

        @test Arblib.realref(z) isa ArbRef
        @test Arblib.realref(z) == x
        @test real(z) == x
        @test Arblib.imagref(z) isa ArbRef
        @test Arblib.imagref(z) == y
        @test imag(z) == y
    end

    @testset "midref" begin
        x = Arb(0.25)
        @test Arblib.midref(x) isa ArfRef
        @test startswith(sprint(show, x), "0.250")
        @test Float64(Arblib.midref(x)) isa Float64
        @test Float64(Arblib.midref(x)) == 0.25
        @test Float64(x) == 0.25
        @test sprint(show, x) == sprint(show, x[])
    end

    @testset "radref" begin
        x = Arb(0.25)
        m = Arblib.radref(x)
        @test m isa MagRef
        m[] = 1.0
        @test Float64(m) ≥ 1.0
        @test sprint(show, m) == sprint(show, m[])
    end

    @testset "convert to Float64/ComplexF64" begin
        x = Arb(0.25)
        @test Float64(x) isa Float64
        @test Float64(x) == 0.25
        z = Acb(2.0 + 0.125im)
        @test ComplexF64(z) isa ComplexF64
        @test ComplexF64(z) == 2.0 + 0.125im
    end
end


text:
module TestCBSOData3

    using Test
    using CBSOData3

    include("01_tables_interface.jl")
    include("02_get_meta.jl")
    include("03_get_tables.jl")
end


text:
using StringDistances
using Test

include("distances.jl")
include("modifiers.jl")
include("pairwise.jl")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export card

"""
    card(;kwargs...)
    card(children::Any;kwargs...)
    card(children_maker::Function;kwargs...)

A Card component.
Card to display content
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The children of the element
- `style` (Dict; optional): The inline styles
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `className` (String; optional): The class of the element
- `title` (String; optional): Title of the card
- `info` (String; optional): Additional information about the content of this card.
- `actions` (optional): Actions available on the side menu, button clicks will be outputted to the actionPropName prop of this card. actions has the following type: Array of lists containing elements 'displayName', 'actionPropName'.
Those elements have the following types:
  - `displayName` (String; optional)
  - `actionPropName` (String; optional)s
- `action_click` (String; optional): The action click value
"""
function card(; kwargs...)
        available_props = Symbol[:children, :style, :id, :className, :title, :info, :actions, :action_click]
        wild_props = Symbol[]
        return Component("card", "Card", "dash_carbon_components", available_props, wild_props; kwargs...)
end

card(children::Any; kwargs...) = card(;kwargs..., children = children)
card(children_maker::Function; kwargs...) = card(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related orbit elements.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A., McClain, W. D (2013). Fundamentals of Astrodynamics
#       and Applications. Microcosm Press.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Get the current EOP Data.
eop = read_iers_eop("./eop_IAU1980_2.txt", Val(:IAU1980))

# File: ./src/transformations/orbit_elements.jl
# =============================================

# Function: change_oe_frame
# -------------------------

################################################################################
#                                 TEST RESULTS
################################################################################
#
# Scenario 01
# ===========
#
# Using AGI STK, it was verified that:
#
#   ╔══════════════════════╦═════════════════════════════╦═══════════════╗
#   ║   Reference Frame    ║            TEME             ║     J2000     ║
#   ╠══════════════════════╬═════════════════════════════╬═══════════════╣
#   ║ Epoch                ║ 1-Jun-2016 11:00:00.000 UTC ║ -             ║
#   ║ Mean motion          ║ 14.4 revs/day               ║ 14.4 revs/day ║
#   ║ Eccentricity         ║ 0.001111                    ║ 0.001111      ║
#   ║ Inclination          ║ 98.405˚                     ║ 98.3365˚      ║
#   ║ Argument of Perigee  ║ 90˚                         ║ 90.0604˚      ║
#   ║ RAAN                 ║ 227.336˚                    ║ 227.134˚      ║
#   ║ Mean Anomaly         ║ 320˚                        ║ 320˚          ║
#   ╚══════════════════════╩═════════════════════════════╩═══════════════╝
#
################################################################################

@testset "Function change_oe_frame" begin
    # This is the Amazonia-1 orbit, which has a = 7130.982 km.
    a = 7130.982e3
    e = 0.001111
    i = 98.405 * π / 180
    Ω = 227.336 * π / 180
    ω = 90 * π / 180
    f = M_to_f(e, 320 * π / 180)
    epoch = date_to_jd(2016, 6, 1, 11, 0, 0)

    oe_TEME  = KeplerianElements(0, a, e, i, Ω, ω, f)
    oe_J2000 = change_oe_frame(oe_TEME, TEME(), J2000(), epoch, eop)
    m = get_M(oe_J2000)

    @test oe_J2000.a / 1000    ≈ 7130.982 atol=1e-8
    @test oe_J2000.e           ≈ e        atol=1e-8
    @test oe_J2000.i * 180 / π ≈ 98.3365  atol=1e-4
    @test oe_J2000.Ω * 180 / π ≈ 227.134  atol=1e-3
    @test oe_J2000.ω * 180 / π ≈ 90.0604  atol=1e-4
    @test m * 180/π            ≈ 320      atol=1e-8

    # Test the other signature of `change_oe_frame`.
    oe_J2000_2 = change_oe_frame(
        a,
        e,
        i,
        Ω,
        ω,
        f,
        TEME(),
        J2000(),
        epoch,
        eop
    )

    @test oe_J2000 === oe_J2000_2
end


text:
# Compute a few eigenvalues of the dielectric matrix (q=0,ω=0) iteratively

using DFTK
using Plots
using KrylovKit
using Printf

# Calculation parameters
kgrid = [1, 1, 1]
Ecut = 5

# Silicon lattice
a = 10.26
lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

# change the symmetry to compute the dielectric operator with and without symmetries
model = model_LDA(lattice, atoms, symmetries=false)
basis = PlaneWaveBasis(model; Ecut, kgrid)
scfres = self_consistent_field(basis, tol=1e-14)

# Apply ε† = 1 - χ0 (vc + fxc)
function eps_fun(δρ)
    δV = apply_kernel(basis, δρ; ρ=scfres.ρ)
    χ0δV = apply_χ0(scfres.ham, scfres.ψ, scfres.εF, scfres.eigenvalues, δV)
    δρ - χ0δV
end

# eager diagonalizes the subspace matrix at each iteration
eigsolve(eps_fun, randn(size(scfres.ρ)), 5, :LM; eager=true, verbosity=3)


text:
abstract type AbstractProjector end

struct L2Projector <: AbstractProjector
    fe_values::CellValues
    M_cholesky # ::SuiteSparse.CHOLMOD.Factor{Float64}
    dh::DofHandler
    set::Vector{Integer}
    node2dof_map
end

function L2Projector(fe_values::JuAFEM.Values, interp::Interpolation, grid::JuAFEM.AbstractGrid, set=1:getncells(grid))

    dim, T, shape = typeof(fe_values).parameters

    # Create an internal scalar valued field. This is enough since the projection is done on a component basis, hence a scalar field.
    dh = DofHandler(grid)
    push!(dh, :_, 1, interp)
    _, vertex_dict, _, _ = close!(dh, true)

    M = _assemble_L2_matrix(fe_values, set, dh)  # the "mass" matrix
    M_cholesky = cholesky(M)  # TODO maybe have a lazy eval instead of precomputing? / JB
    return L2Projector(fe_values, M_cholesky, dh, set, vertex_dict[1])
end

function _assemble_L2_matrix(fe_values, set, dh)

    n = JuAFEM.getn_scalarbasefunctions(fe_values)
    M = create_symmetric_sparsity_pattern(dh)
    assembler = start_assemble(M)

    Me = zeros(n, n)
    cell_dofs = zeros(Int, n)

    function symmetrize_to_lower!(K)
       for i in 1:size(K, 1)
           for j in i+1:size(K, 1)
               K[j, i] = K[i, j]
           end
       end
    end

    ## Assemble contributions from each cell
    for cellnum in set
        celldofs!(cell_dofs, dh, cellnum)

        fill!(Me, 0)
        Xe = getcoordinates(dh.grid, cellnum)
        reinit!(fe_values, Xe)

        ## ∭( v ⋅ u )dΩ
        for q_point = 1:getnquadpoints(fe_values)
            dΩ = getdetJdV(fe_values, q_point)
            for j = 1:n
                v = shape_value(fe_values, q_point, j)
                for i = 1:j
                    u = shape_value(fe_values, q_point, i)
                    Me[i, j] += v ⋅ u * dΩ
                end
            end
        end
        symmetrize_to_lower!(Me)
        assemble!(assembler, cell_dofs, Me)
    end
    return M
end

"""
    project(Vector{Vector{<:Tensor}}}, L2Projector, project_to_nodes=true)

Makes a L2 projection of tensor values to the nodes of the grid. This is commonly used for turning values, computed at integration points, to nodal values, so they can be visualized easier. It is also useful for error estimation and recovering more accurate estimations of secondary unknowns.

If the parameter `project_to_nodes` is true, then the projection returns the values in the order of the mesh nodes. If false, it returns the values corresponding to the degrees of freedom for a scalar field over the domain, which is useful if one wants to interpolate the projected values.
"""
function project(
    vars::Array{Array{Tensor{order,dim,T,M},1},1},
    proj::L2Projector;
    project_to_nodes=true) where {order,dim,T,M}

    projected_vals = _project(vars, proj, M)
    proj.node2dof_map
    if project_to_nodes
        # NOTE we may have more projected values than verticies in the mesh => not all values are returned
        num_nodes = length(proj.node2dof_map)
        reordered_vals = [projected_vals[proj.node2dof_map[node], :] for node in 1:num_nodes]
        return Tensor{order,dim,T}.(reordered_vals)
    else
        # convert back to the original tensor type
        return Tensor{order,dim,T}.(eachrow(projected_vals))
    end

    return Tensor{order,dim,T}.(eachrow(projected_vals))
end

function project(
    vars::Array{Array{SymmetricTensor{order,dim,T,M},1},1},
    proj::L2Projector
    ) where {order,dim,T,M}

    projected_vals = _project(vars, proj, M)
    # convert back to the original tensor type
    return SymmetricTensor{order,dim,T}.(eachrow(projected_vals))
end

function _project(vars, proj::L2Projector, M::Integer)
    # Assemble the multi-column rhs, f = ∭( v ⋅ x̂ )dΩ
    # The number of columns corresponds to the length of the data-tuple in the tensor x̂.

    f = zeros(ndofs(proj.dh), M)
    fe_values = proj.fe_values
    n = getnbasefunctions(proj.fe_values)
    fe = zeros(n, M)

    cell_dofs = zeros(Int, n)
    nqp = getnquadpoints(fe_values)

    ## Assemble contributions from each cell
    for cellnum in proj.set
        celldofs!(cell_dofs, proj.dh, cellnum)
        fill!(fe, 0)
        Xe = getcoordinates(proj.dh.grid, cellnum)
        cell_vars = vars[cellnum]
        reinit!(fe_values, Xe)

        for q_point = 1:nqp
            dΩ = getdetJdV(fe_values, q_point)
            qp_vars = cell_vars[q_point]
            for i = 1:n
                v = shape_value(fe_values, q_point, i)
                fe[i, :] += v * qp_vars[1:M] * dΩ
            end
        end

        # Assemble cell contribution
        for (num, dof) in enumerate(cell_dofs)
            f[dof, :] += fe[num, :]
        end
    end

    # solve for the projected nodal values
    return proj.M_cholesky\f

end


text:
using LiterateOrg
tangle_package(joinpath(dirname(@__FILE__), "..", "src", "SphericalOperators.org"), "SphericalOperators")


text:
"""
Simple minded density for weighted data

Compute cdf
Take its derivative

OUT
    densityV  ::  Float64
        mass / interval width
    midV  ::  Float64
        interval mid points
"""
function density_weighted(xV :: Array{T,1}, wtV :: Array{T,1}, dbg :: Bool) where T <: Real

	## Input check
	n = length(xV);
	if dbg
		@assert (all(wtV .>= 0))
		@assert (size(wtV) == size(xV))
	end

	## Main

	# Construct cdf
	# Drop too small weights
	idxV = findall(wtV .> 1e-6);
	n = length(idxV);
	if n < 2
		error("No observations");
	end

	sortM = sortslices([xV[idxV] wtV[idxV]],  dims = 1);
	cumWtV = cumsum(sortM[:,2]) ./ sum(sortM[:,2]);
	xSortV = sortM[:,1];

	# Density = d(cumulative weight) / dx
	densityV = (cumWtV[2:n] - cumWtV[1 : (n-1)]) ./ (xSortV[2:n] - xSortV[1 : (n-1)]);
	xMidV = 0.5 .* (xSortV[2:n] + xSortV[1 : (n-1)]);

	## Output check
	if dbg > 10
		@assert (all(densityV .> 0));
	end
	return densityV, xMidV
end


text:
using Flux
using Flux.Tracker
using Flux: @epochs

using JuliaRL

import JuliaRL.start!, JuliaRL.step!

TileCoder = JuliaRL.FeatureCreators.TileCoder
# using JuliaRL.Environments
using Random
using LinearAlgebra
using ProgressMeter

function test_flux()

    W = param(rand(2, 5))
    b = param(rand(2))

    predict(x) = W*x .+ b
    loss(x, y) = Flux.mse(predict(x), y)

    W_opt = rand(2,5)
    b_opt = rand(2)

    get_opt(x) = W_opt*x .+ b_opt .+ 0.001*randn(length(size(x)) == 1 ? 2 : (2, size(x)[2]))

    X = [rand(5) for i = 1:1000]
    Y = [get_opt(X[i]) for i = 1:1000]
    data = zip(X, Y)

    train_X = rand(5, 100)
    train_Y = get_opt(train_X)

    evalcb() = @show(loss(train_X, train_Y))

    @epochs 10 Flux.train!(loss, data, ADAM([W, b], 0.001), cb = Flux.throttle(evalcb, 5))

end

function mountain_car_test_flux(α=0.5/8, ϵ=0.1, tilings=32, tiles=2)

    ## Inefficient version of tilecoded q-learning

    ϵ = 0.1
    α = 0.5/tilings

    feat_size = (tilings*(tiles+1)^2)*3
    weights = param(zeros((tilings*(tiles+1)^2)*3))
    iht = TileCoder.IHT(tilings*(tiles+1)^2)

    Q(ϕ) = weights'*ϕ
    function make_feats(t, action)
        feats = zeros(feat_size)
        feats[t .+ ((tilings*(tiles+1)^2)*action + 1)] .= 1
        return feats
    end

    watkins_q_target(t, r) = r + maximum([Q(make_feats(t, a)) for a = 0:2])
    get_action(t) = findmax([Q(make_feats(t, a)) for a = 0:2])[2] - 1
    loss(ϕ, y) = Flux.mse(Q(ϕ), y)

    optimizer = Descent(α)
    env = MountainCar(Random.GLOBAL_RNG)
    # env_ns = MountainCar
    cumulative_reward_array = zeros(Int64, 1000)
    @showprogress 0.1 "Episode: " for episode = 1:1000
        terminal = false
        num_steps = 0
        cumulative_reward = 0
        start!(env)
        action = 0

        t = TileCoder.tiles!(iht, 8, normalized_state(env).*4)
        # println(weights)
        while !terminal

            action = get_action(t)
            # println(action)
            ϕ = make_feats(t, action)
            # println([Q(make_feats(t, a)).data for a = 0:2])
            if rand() < ϵ
                action = rand(0:2)
            end

            state, reward, terminal = step!(env, action)
            t_prime = TileCoder.tiles!(iht, tilings, normalized_state(env).*tiles)
            target = watkins_q_target(t_prime, reward).data

            Flux.train!(loss, [weights], [(ϕ, target)], optimizer)

            num_steps += 1
            cumulative_reward += reward

            t = copy(t_prime)
        end
        cumulative_reward_array[episode] = cumulative_reward
        # println("Episode: $episode, Steps: $num_steps, Reward: $cumulative_reward")
    end
    return cumulative_reward_array
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
#############################################################################
# real_imag.jl
# Handles real and imaginary part of the variables, constants
# and expressions.
#############################################################################

import Base.real, Base.imag
export real, imag
export sign, monotonicity, curvature, evaluate

### Real
struct RealAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function RealAtom(x::AbstractExpr)
        children = (x,)
        return new(:real, hash(children), children, x.size)
    end
end

function sign(x::RealAtom)
    if sign(x.children[1]) == ComplexSign()
        return NoSign()
    else
        return sign(x.children[1])
    end
end

function monotonicity(x::RealAtom)
    return (Nondecreasing(),)
end

function curvature(x::RealAtom)
    return ConstVexity()
end

function evaluate(x::RealAtom)
    return real.(evaluate(x.children[1]))
end

function conic_form!(x::RealAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, x)
        new_objective = ConicObj()
        objective = conic_form!(x.children[1], unique_conic_forms)

        for var in keys(objective)
            re = real.(objective[var][1])
            im = real.(objective[var][2])
            new_objective[var] = (re,im)
        end

        cache_conic_form!(unique_conic_forms, x, new_objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

real(x::AbstractExpr) = RealAtom(x)
real(x::Value) = RealAtom(Constant(x))

### Imaginary
struct ImaginaryAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function ImaginaryAtom(x::AbstractExpr)
        children = (x,)
        return new(:imag, hash(children), children, x.size)
    end
end

function sign(x::ImaginaryAtom)
    sign(x.children[1]) == ComplexSign()
    return NoSign()
end

function monotonicity(x::ImaginaryAtom)
    return (Nondecreasing(),)
end

function curvature(x::ImaginaryAtom)
    return ConstVexity()
end

function evaluate(x::ImaginaryAtom)
    return imag.(evaluate(x.children[1]))
end

function conic_form!(x::ImaginaryAtom, unique_conic_forms::UniqueConicForms=UniqueConicForms())
    if !has_conic_form(unique_conic_forms, x)
        new_objective = ConicObj()
        objective = conic_form!(x.children[1], unique_conic_forms)

        for var in keys(objective)
            re = imag.(objective[var][1])
            im = imag.(objective[var][2])
            new_objective[var] = (re,im)
        end
        cache_conic_form!(unique_conic_forms, x, new_objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

imag(x::AbstractExpr) = ImaginaryAtom(x)
imag(x::Value) = ImaginaryAtom(Constant(x))


text:
mutable struct ConstantMean{T<:Real} <: PriorMean{T}
    C::T
    opt::Optimizer
end

"""
**ConstantMean**
```julia
ConstantMean(c::T=1.0;opt::Optimizer=Adam(α=0.01))
```

Construct a prior mean with constant `c`
Optionally set an optimizer `opt` (`Adam(α=0.01)` by default)
"""
function ConstantMean(c::T=1.0;opt::Optimizer=Adam(α=0.01)) where {T<:Real}
    ConstantMean{T}(c,opt)
end

function update!(μ::ConstantMean{T},grad::AbstractVector{T}) where {T<:Real}
    μ.C += update(μ.opt,sum(grad))
end

(μ::ConstantMean{T})(x::Real) where {T<:Real} = μ.C
(μ::ConstantMean{T})(x::AbstractArray) where {T<:Real} = μ.C*ones(T,size(X))


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

__precompile__(true)

module UUIDs

using Random

export UUID, uuid1, uuid4, uuid_version

import Base: UUID

"""
    uuid_version(u::UUID) -> Int

Inspects the given UUID and returns its version
(see [RFC 4122](https://www.ietf.org/rfc/rfc4122)).

# Examples
```jldoctest
julia> uuid_version(uuid4())
4
```
"""
uuid_version(u::UUID) = Int((u.value >> 76) & 0xf)

"""
    uuid1([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 1 (time-based) universally unique identifier (UUID), as specified
by RFC 4122. Note that the Node ID is randomly generated (does not identify the host)
according to section 4.5 of the RFC.

# Examples
```jldoctest; filter = r"[a-z0-9]{8}-([a-z0-9]{4}-){3}[a-z0-9]{12}"
julia> rng = MersenneTwister(1234);

julia> uuid1(rng)
fcbd9b64-1bc2-11e8-1f13-43a2532b2fa8
```
"""
function uuid1(rng::AbstractRNG=Random.GLOBAL_RNG)
    u = rand(rng, UInt128)

    # mask off clock sequence and node
    u &= 0x00000000000000003fffffffffffffff

    # set the unicast/multicast bit and version
    u |= 0x00000000000010000000010000000000

    # 0x01b21dd213814000 is the number of 100 nanosecond intervals
    # between the UUID epoch and Unix epoch
    timestamp = round(UInt64, time() * 1e7) + 0x01b21dd213814000
    ts_low = timestamp & typemax(UInt32)
    ts_mid = (timestamp >> 32) & typemax(UInt16)
    ts_hi = (timestamp >> 48) & 0x0fff

    u |= UInt128(ts_low) << 96
    u |= UInt128(ts_mid) << 80
    u |= UInt128(ts_hi) << 64

    UUID(u)
end

"""
    uuid4([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 4 (random or pseudo-random) universally unique identifier (UUID),
as specified by RFC 4122.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> uuid4(rng)
196f2941-2d58-45ba-9f13-43a2532b2fa8
```
"""
function uuid4(rng::AbstractRNG=Random.GLOBAL_RNG)
    u = rand(rng, UInt128)
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000040008000000000000000
    UUID(u)
end

end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function rankdense{T<:Real}(a::Array{T,1})
    indexin(a, unique(a))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Napier
using Base.Test
using FactCheck

facts("Testing basics") do
  for (i, (n, e)) in enumerate([
    (0.0001, -9.210340),
    (0.001, -6.907755),
    (0.01, -4.605169),
    (0.1, -2.302584),
    (1, 0),
    (2, 0.693147),
    (e, 1), # e=EULER CONSTANT
    (3, 1.098612),
    (4, 1.386293),
    (5, 1.609437),
    (6, 1.791759),
    (7, 1.945910),
    (8, 2.079441),
    (9, 2.197224),
    (10, 2.302584),
    (20, 2.995732),
    (30, 3.401197),
    (40, 3.688879),
    (50, 3.912023),
    (60, 4.094344),
    (70, 4.248495),
    (80, 4.382026),
    (90, 4.499809),
    (100, 4.605169),
    (200, 5.298317),
    (300, 5.703782),
    (400, 5.991464),
    (500, 6.214607),
    (600, 6.396929),
    (700, 6.551080),
    (800, 6.684611),
    (900, 6.802395),
    (1000, 6.907755),
    (10000, 9.210340)
  ])
    r = ln(n)
    r2 = trunc(r, 6)
    @fact r2 --> e "Mismatch at [$i] r => ln($n)=$r, r2 => trunc(r, 6))=$r2"
  end

  @fact_throws DomainError ln(0)
  @fact_throws DomainError ln(-0.5)
  @fact_throws DomainError ln(-5)
end

@test 1 == 1


text:
using Compose
using Colors

tomato_bisque =
           compose(context(),
                   (context(), circle(), fill("bisque")),
                   (context(), rectangle(), fill("tomato")))
img = SVG("introspect.svg", 4inch, 4(sqrt(3)/2)inch)

introspect(tomato_bisque) |> img


text:
"""
    StandardDriver(csampler::AbstractContextSampler, policy::AbstractPolicy,
                   rsampler::AbstractRewardSampler[,
                   metrics::Tuple{Vararg{<:AbstractMetric}}])

A simple driver that samples contexts, passes them to the policy to generate actions,
then observes the rewards.

# Arguments
- `csampler::AbstractContextSampler`: Context sampler.
- `policy::AbstractPolicy`: Policy to generate actions given contexts.
- `rsampler::AbstractRewardSampler`: Sampler for rewards, given the contexts and actions.
- `metrics::Tuple{Vararg{<:AbstractMetric}}`: A tuple of metrics that will each be
    called as `metric(X, a, r)`.
"""
mutable struct StandardDriver{
    T1<:AbstractContextSampler,
    T2<:AbstractRewardSampler,
    T3<:AbstractPolicy,
    T4<:Tuple{Vararg{<:AbstractMetric}},
} <: AbstractDriver
    csampler::T1
    policy::T3
    rsampler::T2
    metrics::T4
end

function StandardDriver(
    csampler::AbstractContextSampler,
    policy::AbstractPolicy,
    rsampler::AbstractRewardSampler,
)
    return StandardDriver(csampler, policy, rsampler, ())
end

function (driver::StandardDriver)(batch_size::Int64)
    X = driver.csampler(batch_size)
    a = driver.policy(X)
    r = driver.rsampler(X, a)
    for met in driver.metrics
        met(X, a, r)
    end
    return X, a, r
end

mutable struct LatentDriver{
    T1<:AbstractContextSampler,
    T2<:AbstractRewardSampler,
    T3<:AbstractPolicy,
    T4<:Tuple{Vararg{<:AbstractMetric}},
    T5<:Function,
} <: AbstractDriver
    csampler::T1
    policy::T3
    rsampler::T2
    tform::T5
    metrics::T4
end

function LatentDriver(
    csampler::AbstractContextSampler,
    policy::AbstractPolicy,
    rsampler::AbstractRewardSampler,
    tform::Function,
)
    return LatentDriver(csampler, policy, rsampler, tform, ())
end

function (driver::LatentDriver)(batch_size::Int64)
    Z = driver.csampler(batch_size)
    X = mapslices(driver.tform, Z; dims=1)
    a = driver.policy(X)
    r = driver.rsampler(Z, a)
    for met in driver.metrics
        met(Z, a, r)
    end
    return X, a, r
end

"""
    run!(num_batches::Int64, batch_size::Int64, driver::AbstractDriver; <keyword arguments>)

Run a `driver` for `num_batches` batches.
"""
function run!(
    num_batches::Int64, batch_size::Int64, driver::AbstractDriver; verbose::Bool=true
)
    if num_batches <= 0 || batch_size <= 0
        throw(ArgumentError("num_batches and batch_size must be positive"))
    end
    for i in 1:num_batches
        if verbose
            print("\rStep $i/$num_batches")
        end
        X, a, r = driver(batch_size)
        update!(driver.policy, X, a, r)
    end
end


text:
abstract type AbstractRegularization end
abstract type AbstractRegularizer end
# Regularization API #
# ------------------------------------------------------------
initialize_regularization!(lshaped::AbstractLShaped) = initialize_regularization!(lshaped, lshaped.regularization)
restore_regularized_master!(lshaped::AbstractLShaped) = restore_regularized_master!(lshaped, lshaped.regularization)
filter_variables!(lshaped::AbstractLShaped, list) = filter_variables!(lshaped.regularization, list)
filter_constraints!(lshaped::AbstractLShaped, list) = filter_constraints!(lshaped.regularization, list)
log_regularization!(lshaped::AbstractLShaped) = log_regularization!(lshaped, lshaped.regularization)
log_regularization!(lshaped::AbstractLShaped, t::Integer) = log_regularization!(lshaped, t, lshaped.regularization)
take_step!(lshaped::AbstractLShaped) = take_step!(lshaped, lshaped.regularization)
decision(lshaped::AbstractLShaped) = decision(lshaped, lshaped.regularization)
objective_value(lshaped::AbstractLShaped) = objective_value(lshaped, lshaped.regularization)
gap(lshaped::AbstractLShaped) = gap(lshaped, lshaped.regularization)
process_cut!(lshaped::AbstractLShaped, cut::AbstractHyperPlane) = process_cut!(lshaped, cut, lshaped.regularization)
# ------------------------------------------------------------
# Attributes #
# ------------------------------------------------------------
"""
    RawRegularizationParameter

An optimizer attribute used for raw parameters of the regularizer. Defers to `RawParameter`.
"""
struct RawRegularizationParameter <: RegularizationParameter
    name::Any
end
"""
    RegularizationPenaltyterm

An optimizer attribute used to set the proximal term in regulariztion procedures (RD/LV). Options are:

- [`Quadratic`](@ref) (default)
- [`Linearized`](@ref)
- [`InfNorm`](@ref)
- [`ManhattanNorm`](@ref)
"""
struct RegularizationPenaltyterm <: RegularizationParameter end

include("common.jl")
include("no_regularization.jl")
include("rd.jl")
include("tr.jl")
include("lv.jl")


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Gaston test framework
#
# There are two kinds of tests:
#   - Tests that should error -- things that Gaston should not allow, such as
#     building a figure with an invalid p = plotstyle.
#   - Tests that should succeed -- things that should work without producing
#     an error.

using Gaston, Test

@testset "Figure and set commands" begin
    closeall()
    set(reset=true)
    set(mode="null")
    # figures
    @test figure() == 1
    @test figure() == 2
    @test figure(4) == 4
    @test figure() == 3
    @test closefigure(4) == 3
    @test closefigure() == 2
    @test closeall() == 2
    @test begin
        closeall()
        figure()
        figure()
        figure(4);
        closefigure(1,2)
    end == 4
    @test begin
        closeall()
        figure()
        figure()
        figure(4);
        closefigure(4)
    end == 2
    @test begin
        closeall()
        figure()
        figure()
        figure(4);
        closefigure(1,2,3,4,5)
    end == nothing
    @test begin
        closeall()
        figure(3)
        figure(3)
    end == 3
    if !Sys.iswindows()
        @test set(term="x11") === nothing # This test does not pass in Windows
    end
    @test set(termopts="noenhanced") === nothing
    @test set(debug=true) === nothing
    @test set(debug=false) === nothing
end

@testset "2-D plots" begin
    closeall()
    set(reset=true)
    set(mode="null")
    @test closeall() == 0
    @test begin
        plot(1:10)
    end isa Gaston.Figure
    @test begin
        plot(1:10,handle=2)
    end isa Gaston.Figure
    @test begin
        plot(1:10,handle=4)
    end isa Gaston.Figure
    @test closeall() == 3
    @test begin
        plot(1:10,cos)
        Gaston.gnuplot_state.gp_error
    end == false
    @test begin
        plot!(1:10,cos)
        Gaston.gnuplot_state.gp_error
    end == false
    @test begin
        stem(1:10,cos)
        Gaston.gnuplot_state.gp_error
    end == false
    @test begin
        set(termopts="feed noenhanced ansi")
        plot(sin.(-3:0.01:3),
             plotstyle = "linespoints",
             legend = :sine,
             linecolor = "'blue'",
             linewidth = 2,
             pointtype = "ecircle",
             pointsize = "1.1",
             fillstyle = "",
             Axes(
                  title = :test_plot_1,
                  style = "line 20 dt '-.-'",
                  xlabel = :x,
                  ylabel = :y,
                  key = "inside horizontal left top",
                  axis ="loglog",
                  xrange = "[2:3]",
                  yrange = "[2:3]",
                  xzeroaxis = :off,
                  yzeroaxis = :on,
                  font = "'Arial, 12'",
                  size = "79,24"
                 )
            )
    end isa Gaston.Figure
    @test begin
        plot(range(-10,10,length=100),sin,
             w = :lp,
             ls = 15,
             pi = -10,
             Axes(style = "line 15 dt '-...-' lc 'red' lw 3 pt 6"))
    end isa Gaston.Figure
    set(termopts="")
    @test begin
        plot!(cos.(-3:0.01:3),
              legend = "'cos'",
              plotstyle = "linespoints",
              linecolor = "'red'",
              linewidth = 2,
              dashtype = "'.'",
              pointtype = "fcircle",
              pointsize = "1"
             )
    end isa Gaston.Figure
    for ptype in ('λ', 'a', '%', "⋅", "plus")
        @test begin
            plot!(cos.(-3:0.01:3),
                  legend = :cos,
                  w = :linespoints,
                  lc = :red,
                  lw = 2,
                  dt = :.,
                  pt = ptype,
                  ps = 1
                 )
        end isa Gaston.Figure
    end
    @test plot(1:10,cos) isa Gaston.Figure
    @test plot!(1:10,cos) isa Gaston.Figure
    @test stem(1:10,cos) isa Gaston.Figure
    @test plot(1:10, Axes(xrange = "[:3.4]")) isa Gaston.Figure
    @test plot(1:10, Axes(xrange = "[3.4:*]")) isa Gaston.Figure
    @test begin
        plot(1:3, 4:6, plotstyle="points", pointsize=3,
             Axes(xrange="[2.95:3.05]", yrange="[3.95:4.045]"))
    end isa Gaston.Figure
    @test begin
        plot(rand(10) .+ im.*rand(10),
             w = :l,
             leg = :complex_vector,
             Axes(title = :complex_vector,
                  xtics = 1:2:10)
            )
    end isa Gaston.Figure
    @test begin
        y=rand(10)
        ylow=y.-rand(10)
        yhigh=y.+3rand(10)
        plot(y, supp=[ylow yhigh], ps=:errorlines)
    end isa Gaston.Figure
    @test begin
        x=6rand(10).-3
        y=6rand(10).-3
        xdelta=0.5rand(10)
        ydelta=rand(10).-1
        plot(x, y, supp=[xdelta ydelta], ps="vectors head filled", lt=2)
    end isa Gaston.Figure
    @test_skip begin
        x=250randn(10)
        y=10randn(10)
        c1=randn(10)
        c2=50randn(10)
        p = plot(x,y,supp=[c1 c2], ps=:parallelaxes)
    end isa Gaston.Figure
    @test_skip begin
        x=1:5
        y=1:5
        l=["1" "ss" "3" "yy" "x5"]
        p = plot(x, y, supp=l, ps=:labels)
    end isa Gaston.Figure
    @test begin
        stem(0:0.1:3pi,sin)
    end isa Gaston.Figure
    @test begin
        stem(0:0.1:3pi,sin,onlyimpulses=true)
    end isa Gaston.Figure
    @test begin
        t=0:0.05:3
        y1=exp.(-t)+0.1rand(length(t))
        y2=exp.(-t)-0.1rand(length(t))
        p = plot(t,y1,supp=y2,ps=:filledcurves,lc="brown")
        p = plot!(t,exp.(-t),lc="black")
    end isa Gaston.Figure
    # This test is cool, but it fails way too often due to slight
    # OS, configuration, and gnuplot version differences. Disabled for now.
    @test_skip begin
        set(reset=true)
        set(terminal="dumb", size="27,13")
        x="\f                           \n  11 +-----------------+   \n  10 |-+ + + + + + + +*|   \n   9 |-+          ****-|   \n   7 |-+        **   +-|   \n   6 |-+      **     +-|   \n   5 |-+    **       +-|   \n   4 |-+  **         +-|   \n   3 |****           +-|   \n   1 |-+ + + + + + + +-|   \n   0 +-----------------+   \n     1 2 3 4 5 6 7 8 9 10  \n                           \n"
        a=repr("text/plain", p = plot(1:10))
        a == x
    end == true
    @test begin
        set(reset=true)
        stem(rand(10))
    end isa Gaston.Figure
    @test begin
        stem(rand(10),onlyimpulses=true)
    end isa Gaston.Figure
    @test begin
        scatter(rand(10),rand(10))
    end isa Gaston.Figure
    @test begin
        scatter(randn(10), randn(10), pointtype="λ")
    end isa Gaston.Figure
    @test begin
        scatter(rand(10),rand(10), lc=:red)
    end isa Gaston.Figure
    @test begin
        scatter(rand(10),rand(10), lc=:red, Axes(title=:test))
    end isa Gaston.Figure
    @test begin
        scatter(complex.(rand(10), rand(10)), lc=:red)
    end isa Gaston.Figure
    @test begin
        scatter(complex.(rand(10), rand(10)), lc=:red, Axes(title="'test'"))
    end isa Gaston.Figure
    @test begin
        scatter(complex.(rand(10), rand(10)),
                linecolor="'green'",
                Axes(title = :test))
        scatter!(complex.(rand(10), rand(10)),linecolor="'blue'")
    end isa Gaston.Figure
    @test begin
        stem(exp.(0:0.01:1))
    end isa Gaston.Figure
    @test begin
        stem(exp.(0:0.01:1),onlyimpulses=true)
    end isa Gaston.Figure
    @test begin
        t=0:0.01:1
        stem(t,exp)
    end isa Gaston.Figure
    @test begin
        t=0:0.01:1
        stem(t,exp,Axes(xlabel=:x),lc=:red)
    end isa Gaston.Figure
    @test begin
        t=0:0.01:1
        stem(exp.(t),Axes(xlabel=:x))
    end isa Gaston.Figure
    @test begin
        bar(1:10)
    end isa Gaston.Figure
    @test begin
        bar(rand(10),Axes(xlabel=:x),lc=:red)
    end isa Gaston.Figure
    @test begin
        bar(11:20,rand(10),Axes(xlabel=:x),lc=:red)
    end isa Gaston.Figure
    if Gaston.gnuplot_state.gnuplot_available
        @test begin
            set(reset=true)
            s = MIME"image/svg+xml"()
            p = plot(1:10)
            a = repr(s, p)
            a[1:29] == "<?xml version=\"1.0\" encoding="
        end == true
        @test begin
            set(reset=true)
            s = MIME"image/svg+xml"()
            p = plot(1:10)
            a = repr(s, p)
            idx = 0
            while a[end-idx] == '\n' || a[end-idx] == '\r'
                idx = idx + 1
                idx > 20 && return false
            end
            a[end-5-idx:end-idx] == "</svg>"
        end == true
    end
end

@testset "Histograms" begin
    closeall()
    set(reset=true)
    set(mode="null")
    @test begin
        histogram(rand(1000))
    end isa Gaston.Figure
    @test begin
        histogram(randn(1000),
                  bins = 100,
                  norm = 1,
                  linecolor = :blue,
                  linewidth = 2,
                  Axes(
                       title = "'test histogram'",
                       xlabel = "'x'",
                       ylabel = "'y'",
                       style = "fill solid",
                       key = "inside horizontal left top",
                       xrange = "[*:*]",
                       yrange = "[*:*]",
                       font = "'Arial, 12'",
                       size = "0.9,0.9")
                 )
    end isa Gaston.Figure
end

@testset "Images" begin
    closeall()
    set(reset=true)
    set(mode="null")
    z = rand(5,6)
    @test begin
        imagesc(z,
                Axes(title = :test_imagesc_1,
                     xlabel = "'xx'",
                     ylabel = "'yy'",
                     font = "'Arial, 12'")
               )
    end isa Gaston.Figure
    @test begin
        imagesc(5:10, 21:25, z,
                Axes(title=:test_imagesc_2,
                     xlabel=:xx,
                     ylabel=:yy)
               )
    end isa Gaston.Figure
    @test begin
        R = [x+y for x=0:5:120, y=0:5:120]
        G = [x+y for x=0:5:120, y=120:-5:0]
        B = [x+y for x=120:-5:0, y=0:5:120]
        Z = zeros(3,25,25)
        Z[1,:,:] = R
        Z[2,:,:] = G
        Z[3,:,:] = B
        imagesc(Z, Axes(title="'RGB Image'"))
    end isa Gaston.Figure
end

@testset "3-D plots" begin
    closeall()
    set(reset=true)
    set(mode="null")
    @test begin
        surf(rand(10,10))
    end isa Gaston.Figure
    @test begin
        x=[0,1,2]; y=[0,1,2,3]; Z=[10 10 10; 10 5 10; 10 1 10; 10 0 10]
        surf(x, y, Z,
             legend = :test,
             w = :lines,
             lc = :black,
             Axes(
                  title = :test,
                  xlabel = :x,
                  ylabel = :y,
                  zlabel = :z,
                  key = "inside horizontal left top",
                  xzeroaxis = :on,
                  yzeroaxis = :on,
                  zzeroaxis = :on,
                  font = "'Arial, 12'",
                  size = "0.9,0.9",
                  view = (90,60))
            )
    end isa Gaston.Figure
    @test begin
        surf(0:9,2:11,(x,y)->x*y)
    end isa Gaston.Figure
    @test begin
        surf(0:9,2:11,(x,y)->x*y,
             legend = "'test'",
             plotstyle="pm3d",
             Axes(title="'test'")
            )
    end isa Gaston.Figure
    @test begin
        surf(0:9,2:11,(x,y)->x*y)
        surf!(0:9,2:11,(x,y)->x/y)
    end isa Gaston.Figure
    @test begin
        scatter3(rand(10),rand(10),rand(10),pointtype="fcircle",linecolor=:red)
    end isa Gaston.Figure
    @test begin
        scatter3(rand(10),rand(10),rand(10),
                 Axes(title=:Scatter3, xlabel=:x),
                 lc=:red)
        scatter3!(rand(8),rand(8),rand(8),lc=:black)
    end isa Gaston.Figure
    @test begin
        scatter3(0:0.1:10pi, cos, sin, pt=7, lc="palette",
                 Axes(title=:Spiral,xlabel=:x,ylabel=:y,zlabel=:z))
    end isa Gaston.Figure
    @test begin
        x = 0:0.1:10pi
        scatter3(x, cos, sin, pt=7, lc="palette", ps="variable",
                 supp = x./20,
                 Axes(title=:Spiral,xlabel=:x,ylabel=:y,zlabel=:z))
    end isa Gaston.Figure
    @test begin
        x = -1:0.05:1
        y = -1.5:0.05:2
        egg(x,y) = x^2 + y^2/(1.4 + y/5)^2
        segg = [egg(x,y) for x in x, y in y]
        a = Axes(auto="fix",
                 size="ratio -1",
                 cntrparam="levels incremental 0,0.02,1",
                 palette=:cool)
        contour(x,y,segg',a,w=:lines,lc="palette",labels=false)
    end isa Gaston.Figure
    @test begin
        x = y = -5:0.1:5
        heatmap(x,y,(x,y)->cos.(x/2).*sin.(y/2))
    end isa Gaston.Figure
    @test begin
        x = y = -5:0.1:5
        heatmap(x,y,(x,y)->cos.(x/2).*sin.(y/2), Axes(title=:heatmap))
    end isa Gaston.Figure
end

@testset "Multiplot" begin
    @test begin
        closeall()
        x = y = -15:0.33:15
        z = rand(5,6)
        p1 = scatter(rand(10), rand(10), Axes(title=:p1), handle=1);
        p2 = imagesc(z,Axes(title = :p2), handle=2);
        p3 = surf(x, y, (x,y)->sin.(sqrt.(x.*x+y.*y))./sqrt.(x.*x+y.*y),
                  Axes(title=:p3), plotstyle="pm3d",handle=3);
        plot([p1 p2 ; nothing p3])
    end isa Gaston.Figure
end

@testset "Saving plots" begin
    closeall()
    set(reset=true)
    set(mode="null")
    p = plot(1:10)
    filename = tempname()
    @test begin
        save(output=filename,term="pdf")
    end == nothing
    @test begin
        save(handle=1,term="png",output=filename)
    end == nothing
    @test begin
        save(term="eps",output=filename)
    end == nothing
    @test begin
        save(term="pdf",
             output=filename,
             font = "Arial, 12",
             size = "5,3",
             linewidth = 3)
    end == nothing
    @test begin
        save(term="svg",output=filename,linewidth=3)
    end == nothing
    @test begin
        save(term="gif",output=filename,size="640,480")
    end == nothing
end

@testset "Tests that should fail" begin
    closeall()
    set(reset=true)
    set(mode="null")
    # figure-related
    @test_throws MethodError figure("invalid")
    @test_throws MethodError figure(1.0)
    @test_throws MethodError figure(1:2)
    @test_throws DomainError closefigure(-1)
    @test_throws TypeError closefigure("invalid")
    @test_throws TypeError closefigure(1.0)
    @test_throws TypeError closefigure(1:2)
    # plot
    @test_throws ArgumentError p = plot(0:10,0:11)
    #@test_throws DimensionMismatch surf([1,2],[3,4],[5,6,7])
    #@test_throws DimensionMismatch surf([1,2,3],[3,4],[5,6,7])
    #plot!
    closeall()
    @test_throws ErrorException p = plot!(0:10)
    # imagesc
    z = rand(5,6)
    @test_throws ArgumentError imagesc(1:5,1:7,z)
    # histogram
    #@test_throws MethodError histogram(0:10+im*0:10)
    # save
    @test_throws DomainError begin
        save(handle=2,term="png",output="test")
        Gaston.gnup = plot_state.gp_error
    end
    @test_throws ErrorException p = plot!(rand(10),handle=10)
    # set
    @test_throws MethodError set(color=3)
    @test_throws TypeError set(debug="oo")
    closeall()
end


text:
using DrWatson
@quickactivate :Swalbe
using CSV
# Take default values
sys = Swalbe.SysConst_1D()
# Kapital gamma, mobility M, diffusion coefficient D and dumping interval
GG = 0.5
Mob = 1.0
Dif = 0.1
time_dump = 100

"""
    run_active_thin_film()

Simulation with `active` colloids given by a field rho.
"""
function run_active_thin_film(
    sys::Swalbe.SysConst_1D;
    h₀=1,
    ϵ=0.001,
    θ₀=1/9,
    verbos=true,
    dump = time_dump,
    rho = zeros(sys.L),
    differentials = zeros(sys.L, 4),
    fluid=zeros(sys.Tmax÷dump, sys.L),
    rho_evo=zeros(sys.Tmax÷dump, sys.L),
    Gam = GG,
    mobl = Mob,
    difu = Dif,
    T=Float64
)
    println("Simulating a active thin film")
    fout, ftemp, feq, height, vel, pressure, dgrad, F, slip, h∇p, fluc = Swalbe.Sys(sys, true, T)
    Swalbe.randinterface!(height, h₀, ϵ)
    # Swalbe.sinewave1D!(height, h₀, 1, ϵ, 1)
    rho .= 0.1
    Swalbe.equilibrium!(feq, height, vel)
    ftemp .= feq
    println("Starting the lattice Boltzmann time loop")
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(height)
            if verbos
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end
        # Time update of rho
        Swalbe.update_rho!(rho, fluc, height, dgrad, differentials, M=mobl, D=difu)
        # Film pressure with rho field
        Swalbe.filmpressure!(pressure, height, dgrad, rho, sys.γ, θ₀, sys.n, sys.m, sys.hmin, sys.hcrit, Gamma=Gam)
        Swalbe.∇f!(h∇p, pressure, dgrad, height)
        Swalbe.slippage!(slip, height, vel, sys.δ, sys.μ)
        # HereSwalbe.thermal!(fluc, height, sys.kbt, sys.μ, sys.δ)
        # Here we a force that is like pull of an inclined plane
        F .= h∇p .+ slip # .+ fluc
        Swalbe.equilibrium!(feq, height, vel)
        Swalbe.BGKandStream!(fout, feq, ftemp, -F)
        Swalbe.moments!(height, vel, fout)
        # Make a snaphot of the configuration
        Swalbe.snapshot!(fluid, height, t, dumping = dump)
        Swalbe.snapshot!(rho_evo, rho, t, dumping = dump)
    end

    return fluid, rho_evo# height, vel, rho

end
# Now run the actual simulation
fluid, rho_evo = run_active_thin_film(sys, Gam = 1.0, mobl = 0.1, difu = 0.1)
# Push the results into a dict for writting to CSV
df_fluid = Dict()
df_rho = Dict()
# Fill the dict with the time steps recorded
for t in 1:sys.Tmax÷sys.tdump
    df_fluid["h_$(t*time_dump)"] = fluid[t,:]
    df_rho["rho_$(t*time_dump)"] = rho_evo[t,:]
end
# Write the desired CSV files
CSV.write("height_run_Gamma_$(GG)_Dif_$(Dif)_Mob_$(Mob).csv", df_fluid)
CSV.write("rho_run_Gamma_$(GG)_Dif_$(Dif)_Mob_$(Mob).csv", df_rho)
# Simulation is done and output is save to files height_...csv and rho_...csv
println("Sim done, results in height_{...}.csv and rho_{...}.csv")