text:
#nb # %% slideshow [markdown] {"slideshow": {"slide_type": "slide"}}
#md # # [Lecture 1 - Getting Started](@id lecture1)
#!md # # [Lecture 1 - Getting Started]
#
# **Essential Skills For Computational Biology**
#
# *2020-06-08*
#
# Kevin Bonham, PhD

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Learning Objectives
# After completing this lesson, students will be able to:
# - Navigate their computer's file system using a command line interface
# - Use a plain text editor (VS Code) to modify source code
# - Clone, commit to, and push from a git repository

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## What is a computer program?

#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
function hello(x)
    println("Hello, $(x)!")
end
#-
#nb %% a code block [code] {"slideshow": {"slide_type": "fragment"}}
hello("Students")

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Programs are just things and actions

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - "Things" in computer code are data

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - "Actions" are generally called "functions"

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - Real life is filled with programs

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Computer programs are just things and actions

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - Computer languages are procedural

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - Computer languages are literal

#nb %% a subslide [code] {"slideshow": {"slide_type": "fragment"}}
42 / 2
#-
#nb %% a subslide [code] {"slideshow": {"slide_type": "fragment"}}
## 5 + "2"

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## What are "essential" skills?

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - How do I think about writing a computer program?

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - When the code I've written has an error, what steps do I take to debug it?

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - How do I keep track of the code that I've written?

#nb # %% a subslide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - How do I get help when I'm stuck?

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Course components

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - Readings and exercises from [_Think Julia_](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html)
# ![think julia](https://imgur.com/yGSTdQe.png)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# - Reading and exercises from [BISC195 docs](https://wellesley-bisc195.github.io/BISC195.jl/dev/)
# ![BISC195 website](https://imgur.com/9AgXhjP.png)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# - Problem sets on github classroom
# ![Assignment 01](https://imgur.com/Tb7BLuy.png)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Grading
# ![Assignment 01](https://i.imgur.com/AM6tw55.jpg)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Assignments

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - There will be ~10 "Assignments," delivered on github classroom

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - All assignments will be worth the same amount

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - All assignments will have automated "unit tests" that will help you complete them

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# - All assignments will have "due dates" to help keep you up with course material,
#   but you will still get full credit for anything turned in by the end of the term.

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Course Schedule
#
# **Dates**: June 8 - July 17, 2020
#
# **Lecture**: Tu/Th 10-11am - Reserved for lecture
#
# **Office hours / lab**:
# - M/W 10am-12pm
# - Tu/Th 1-3pm

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# **All students are expected to attend lecture**
#
# **All students should be able to come to at least two lab periods**

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## Let's get started!
#
# - Course bible: https://wellesley-bisc195.github.io/BISC195.jl/stable/
# - [Lesson 1 - Getting Started](https://wellesley-bisc195.github.io/BISC195.jl/stable/Lessons/Lesson1/)


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

"""
Simple unit testing functionality:

* `@test`
* `@test_throws`

All tests belong to a *test set*. There is a default, task-level
test set that throws on the first failure. Users can choose to wrap
their tests in (possibly nested) test sets that will store results
and summarize them at the end of the test set with `@testset`.
"""
:Test # cf. #22288

module Test

export @test, @test_throws, @test_broken, @test_skip, @test_warn, @test_nowarn
export @testset
# Legacy approximate testing functions, yet to be included
export @inferred
export detect_ambiguities
export GenericString, GenericSet, GenericDict

#-----------------------------------------------------------------------

# Backtrace utility functions
function ip_matches_func_and_name(ip, func::Symbol, dir::String, file::String)
    for fr in StackTraces.lookup(ip)
        if fr === StackTraces.UNKNOWN || fr.from_c
            return false
        end
        path = string(fr.file)
        fr.func == func && dirname(path) == dir && basename(path) == file && return true
    end
    return false
end

function scrub_backtrace(bt)
    do_test_ind = findfirst(addr->ip_matches_func_and_name(addr, :do_test, ".", "test.jl"), bt)
    if do_test_ind != 0 && length(bt) > do_test_ind
        bt = bt[do_test_ind + 1:end]
    end
    name_ind = findfirst(addr->ip_matches_func_and_name(addr, Symbol("macro expansion"), ".", "test.jl"), bt)
    if name_ind != 0 && length(bt) != 0
        bt = bt[1:name_ind]
    end
    return bt
end

"""
    Result

All tests produce a result object. This object may or may not be
stored, depending on whether the test is part of a test set.
"""
abstract type Result end

"""
    Pass

The test condition was true, i.e. the expression evaluated to true or
the correct exception was thrown.
"""
struct Pass <: Result
    test_type::Symbol
    orig_expr
    data
    value
end
function Base.show(io::IO, t::Pass)
    print_with_color(:green, io, "Test Passed\n"; bold = true)
    if !(t.orig_expr === nothing)
        print(io, "  Expression: ", t.orig_expr)
    end
    if t.test_type == :test_throws
        # The correct type of exception was thrown
        print(io, "\n      Thrown: ", typeof(t.value))
    elseif t.test_type == :test && isa(t.data, Expr)
        # The test was an expression, so display the term-by-term
        # evaluated version as well
        print(io, "\n   Evaluated: ", t.data)
    end
end

"""
    Fail

The test condition was false, i.e. the expression evaluated to false or
the correct exception was not thrown.
"""
mutable struct Fail <: Result
    test_type::Symbol
    orig_expr
    data
    value
end
function Base.show(io::IO, t::Fail)
    print_with_color(Base.error_color(), io, "Test Failed\n"; bold = true)
    print(io, "  Expression: ", t.orig_expr)
    if t.test_type == :test_throws_wrong
        # An exception was thrown, but it was of the wrong type
        print(io, "\n    Expected: ", t.data)
        print(io, "\n      Thrown: ", isa(t.data, Type) ? typeof(t.value) : t.value)
    elseif t.test_type == :test_throws_nothing
        # An exception was expected, but no exception was thrown
        print(io, "\n    Expected: ", t.data)
        print(io, "\n  No exception thrown")
    elseif t.test_type == :test && isa(t.data, Expr)
        # The test was an expression, so display the term-by-term
        # evaluated version as well
        print(io, "\n   Evaluated: ", t.data)
    end
end

"""
    Error

The test condition couldn't be evaluated due to an exception, or
it evaluated to something other than a [`Bool`](@ref).
In the case of `@test_broken` it is used to indicate that an
unexpected `Pass` `Result` occurred.
"""
mutable struct Error <: Result
    test_type::Symbol
    orig_expr
    value
    backtrace
end
function Base.show(io::IO, t::Error)
    print_with_color(Base.error_color(), io, "Error During Test\n"; bold = true)
    if t.test_type == :test_nonbool
        println(io, "  Expression evaluated to non-Boolean")
        println(io, "  Expression: ", t.orig_expr)
        print(  io, "       Value: ", t.value)
    elseif t.test_type == :test_error
        println(io, "  Test threw an exception of type ", typeof(t.value))
        println(io, "  Expression: ", t.orig_expr)
        # Capture error message and indent to match
        errmsg = sprint(showerror, t.value, scrub_backtrace(t.backtrace))
        print(io, join(map(line->string("  ",line),
                            split(errmsg, "\n")), "\n"))
    elseif t.test_type == :test_unbroken
        # A test that was expected to fail did not
        println(io, " Unexpected Pass")
        println(io, " Expression: ", t.orig_expr)
        println(io, " Got correct result, please change to @test if no longer broken.")
    elseif t.test_type == :nontest_error
        # we had an error outside of a @test
        println(io, "  Got an exception of type $(typeof(t.value)) outside of a @test")
        # Capture error message and indent to match
        errmsg = sprint(showerror, t.value, t.backtrace)
        print(io, join(map(line->string("  ",line),
                            split(errmsg, "\n")), "\n"))
    end
end

"""
    Broken

The test condition is the expected (failed) result of a broken test,
or was explicitly skipped with `@test_skip`.
"""
mutable struct Broken <: Result
    test_type::Symbol
    orig_expr
end
function Base.show(io::IO, t::Broken)
    print_with_color(Base.warn_color(), io, "Test Broken\n"; bold = true)
    if t.test_type == :skipped && !(t.orig_expr === nothing)
        println(io, "  Skipped: ", t.orig_expr)
    elseif !(t.orig_expr === nothing)
        println(io, "Expression: ", t.orig_expr)
    end
end

#-----------------------------------------------------------------------

abstract type ExecutionResult end

struct Returned <: ExecutionResult
    value
    data
end

struct Threw <: ExecutionResult
    exception
    backtrace
end

function eval_test(evaluated::Expr, quoted::Expr)
    res = true
    i = 1
    args = evaluated.args
    quoted_args = quoted.args
    n = length(args)
    if evaluated.head == :comparison
        while i < n
            a, op, b = args[i], args[i+1], args[i+2]
            if res
                res = op(a, b) === true  # Keep `res` type stable
            end
            quoted_args[i] = a
            quoted_args[i+2] = b
            i += 2
        end

    elseif evaluated.head == :call
        op = args[1]
        res = op(args[2:n]...) === true
        for i in 2:n
            quoted_args[i] = args[i]
        end
    else
        throw(ArgumentError("Unhandled expression type: $(evaluated.head)"))
    end
    Returned(res, quoted)
end

const comparison_prec = Base.operator_precedence(:(==))

"""
    test_expr!(ex, kws...)

Preprocess test expressions of function calls with trailing keyword arguments
so that e.g. `@test a ≈ b atol=ε` means `@test ≈(a, b, atol=ε)`.
"""
test_expr!(m, ex) = ex

function test_expr!(m, ex, kws...)
    ex isa Expr && ex.head == :call || @goto fail
    for kw in kws
        kw isa Expr && kw.head == :(=) || @goto fail
        kw.head = :kw
        push!(ex.args, kw)
    end
    return ex
@label fail
    error("invalid test macro call: $m $ex $(join(kws," "))")
end

# @test - check if the expression evaluates to true
"""
    @test ex
    @test f(args...) key=val ...

Tests that the expression `ex` evaluates to `true`.
Returns a `Pass` `Result` if it does, a `Fail` `Result` if it is
`false`, and an `Error` `Result` if it could not be evaluated.

The `@test f(args...) key=val...` form is equivalent to writing
`@test f(args..., key=val...)` which can be useful when the expression
is a call using infix syntax such as approximate comparisons:

    @test a ≈ b atol=ε

This is equivalent to the uglier test `@test ≈(a, b, atol=ε)`.
It is an error to supply more than one expression unless the first
is a call expression and the rest are assignments (`k=v`).
"""
macro test(ex, kws...)
    test_expr!("@test", ex, kws...)
    orig_ex = Expr(:inert, ex)
    result = get_test_result(ex)
    :(do_test($result, $orig_ex))
end

"""
    @test_broken ex
    @test_broken f(args...) key=val ...

Indicates a test that should pass but currently consistently fails.
Tests that the expression `ex` evaluates to `false` or causes an
exception. Returns a `Broken` `Result` if it does, or an `Error` `Result`
if the expression evaluates to `true`.

The `@test_broken f(args...) key=val...` form works as for the `@test` macro.
"""
macro test_broken(ex, kws...)
    test_expr!("@test_broken", ex, kws...)
    orig_ex = Expr(:inert, ex)
    result = get_test_result(ex)
    # code to call do_test with execution result and original expr
    :(do_broken_test($result, $orig_ex))
end

"""
    @test_skip ex
    @test_skip f(args...) key=val ...

Marks a test that should not be executed but should be included in test
summary reporting as `Broken`. This can be useful for tests that intermittently
fail, or tests of not-yet-implemented functionality.

The `@test_skip f(args...) key=val...` form works as for the `@test` macro.
"""
macro test_skip(ex, kws...)
    test_expr!("@test_skip", ex, kws...)
    orig_ex = Expr(:inert, ex)
    testres = :(Broken(:skipped, $orig_ex))
    :(record(get_testset(), $testres))
end

# An internal function, called by the code generated by the @test
# macro to get results of the test expression.
# In the special case of a comparison, e.g. x == 5, generate code to
# evaluate each term in the comparison individually so the results
# can be displayed nicely.
function get_test_result(ex)
    orig_ex = Expr(:inert, ex)
    # Normalize non-dot comparison operator calls to :comparison expressions
    if isa(ex, Expr) && ex.head == :call && length(ex.args) == 3 &&
        first(string(ex.args[1])) != '.' &&
        (ex.args[1] === :(==) || Base.operator_precedence(ex.args[1]) == comparison_prec)
        ex = Expr(:comparison, ex.args[2], ex.args[1], ex.args[3])
    end
    if isa(ex, Expr) && ex.head == :comparison
        # pass all terms of the comparison to `eval_comparison`, as an Expr
        escaped_terms = [esc(arg) for arg in ex.args]
        quoted_terms = [QuoteNode(arg) for arg in ex.args]
        testret = :(eval_test(
            Expr(:comparison, $(escaped_terms...)),
            Expr(:comparison, $(quoted_terms...)),
        ))
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] in (:isequal, :isapprox)
        escaped_terms = [esc(arg) for arg in ex.args]
        quoted_terms = [QuoteNode(arg) for arg in ex.args]
        testret = :(eval_test(
            Expr(:call, $(escaped_terms...)),
            Expr(:call, $(quoted_terms...)),
        ))
    else
        testret = :(Returned($(esc(ex)), nothing))
    end
    result = quote
        try
            $testret
        catch _e
            Threw(_e, catch_backtrace())
        end
    end
    Base.remove_linenums!(result)
    result
end

# An internal function, called by the code generated by the @test
# macro to actually perform the evaluation and manage the result.
function do_test(result::ExecutionResult, orig_expr)
    # get_testset() returns the most recently added test set
    # We then call record() with this test set and the test result
    if isa(result, Returned)
        # expr, in the case of a comparison, will contain the
        # comparison with evaluated values of each term spliced in.
        # For anything else, just contains the test expression.
        # value is the evaluated value of the whole test expression.
        # Ideally it is true, but it may be false or non-Boolean.
        value = result.value
        testres = if isa(value, Bool)
            # a true value Passes
            value ? Pass(:test, nothing, nothing, value) :
                    Fail(:test, orig_expr, result.data, value)
        else
            # If the result is non-Boolean, this counts as an Error
            Error(:test_nonbool, orig_expr, value, nothing)
        end
    else
        # The predicate couldn't be evaluated without throwing an
        # exception, so that is an Error and not a Fail
        @assert isa(result, Threw)
        testres = Error(:test_error, orig_expr, result.exception, result.backtrace)
    end
    record(get_testset(), testres)
end

function do_broken_test(result::ExecutionResult, orig_expr)
    testres = Broken(:test, orig_expr)
    # Assume the test is broken and only change if the result is true
    if isa(result, Returned)
        value = result.value
        if isa(value, Bool) && value
            testres = Error(:test_unbroken, orig_expr, value, nothing)
        end
    end
    record(get_testset(), testres)
end

#-----------------------------------------------------------------------

"""
    @test_throws exception expr

Tests that the expression `expr` throws `exception`.
The exception may specify either a type,
or a value (which will be tested for equality by comparing fields).
Note that `@test_throws` does not support a trailing keyword form.
"""
macro test_throws(extype, ex)
    orig_ex = Expr(:inert, ex)
    result = quote
        try
            Returned($(esc(ex)), nothing)
        catch _e
            Threw(_e, nothing)
        end
    end
    Base.remove_linenums!(result)
    :(do_test_throws($result, $orig_ex, $(esc(extype))))
end

# An internal function, called by the code generated by @test_throws
# to evaluate and catch the thrown exception - if it exists
function do_test_throws(result::ExecutionResult, orig_expr::ANY, extype::ANY)
    if isa(result, Threw)
        # Check that the right type of exception was thrown
        success = false
        exc = result.exception
        if isa(extype, Type)
            success = isa(exc, extype)
        else
            if isa(exc, typeof(extype))
                success = true
                for fld in 1:nfields(extype)
                    if !(getfield(extype, fld) == getfield(exc, fld))
                        success = false
                        break
                    end
                end
            end
        end
        if success
            testres = Pass(:test_throws, nothing, nothing, exc)
        else
            testres = Fail(:test_throws_wrong, orig_expr, extype, exc)
        end
    else
        testres = Fail(:test_throws_nothing, orig_expr, extype, nothing)
    end
    record(get_testset(), testres)
end

#-----------------------------------------------------------------------
# Test for warning messages

ismatch_warn(s::AbstractString, output) = contains(output, s)
ismatch_warn(s::Regex, output) = ismatch(s, output)
ismatch_warn(s::Function, output) = s(output)
ismatch_warn(S::Union{AbstractArray,Tuple}, output) = all(s -> ismatch_warn(s, output), S)

"""
    @test_warn msg expr

Test whether evaluating `expr` results in [`STDERR`](@ref) output that contains
the `msg` string or matches the `msg` regular expression.  If `msg` is
a boolean function, tests whether `msg(output)` returns `true`.  If `msg` is a
tuple or array, checks that the error output contains/matches each item in `msg`.
Returns the result of evaluating `expr`.

See also [`@test_nowarn`](@ref) to check for the absence of error output.
"""
macro test_warn(msg, expr)
    quote
        let fname = tempname(), have_color = Base.have_color
            try
                @eval Base have_color = false
                ret = open(fname, "w") do f
                    redirect_stderr(f) do
                        $(esc(expr))
                    end
                end
                @test ismatch_warn($(esc(msg)), readstring(fname))
                ret
            finally
                eval(Base, Expr(:(=), :have_color, have_color))
                rm(fname, force=true)
            end
        end
    end
end

"""
    @test_nowarn expr

Test whether evaluating `expr` results in empty [`STDERR`](@ref) output
(no warnings or other messages).  Returns the result of evaluating `expr`.
"""
macro test_nowarn(expr)
    quote
        @test_warn r"^(?!.)"s $(esc(expr))
    end
end

#-----------------------------------------------------------------------

# The AbstractTestSet interface is defined by two methods:
# record(AbstractTestSet, Result)
#   Called by do_test after a test is evaluated
# finish(AbstractTestSet)
#   Called after the test set has been popped from the test set stack
abstract type AbstractTestSet end

"""
    record(ts::AbstractTestSet, res::Result)

Record a result to a testset. This function is called by the `@testset`
infrastructure each time a contained `@test` macro completes, and is given the
test result (which could be an `Error`). This will also be called with an `Error`
if an exception is thrown inside the test block but outside of a `@test` context.
"""
function record end

"""
    finish(ts::AbstractTestSet)

Do any final processing necessary for the given testset. This is called by the
`@testset` infrastructure after a test block executes. One common use for this
function is to record the testset to the parent's results list, using
`get_testset`.
"""
function finish end

"""
    TestSetException

Thrown when a test set finishes and not all tests passed.
"""
mutable struct TestSetException <: Exception
    pass::Int
    fail::Int
    error::Int
    broken::Int
    errors_and_fails::Vector{Union{Fail, Error}}
end

function Base.show(io::IO, ex::TestSetException)
    print(io, "Some tests did not pass: ")
    print(io, ex.pass,  " passed, ")
    print(io, ex.fail,  " failed, ")
    print(io, ex.error, " errored, ")
    print(io, ex.broken, " broken.")
end

function Base.showerror(io::IO, ex::TestSetException, bt; backtrace=true)
    print_with_color(Base.error_color(), io, string(ex))
end

#-----------------------------------------------------------------------

"""
    FallbackTestSet

A simple fallback test set that throws immediately on a failure.
"""
struct FallbackTestSet <: AbstractTestSet
end
fallback_testset = FallbackTestSet()

mutable struct FallbackTestSetException <: Exception
    msg::String
end

function Base.showerror(io::IO, ex::FallbackTestSetException, bt; backtrace=true)
    print_with_color(Base.error_color(), io, ex.msg)
end

# Records nothing, and throws an error immediately whenever a Fail or
# Error occurs. Takes no action in the event of a Pass or Broken result
record(ts::FallbackTestSet, t::Union{Pass,Broken}) = t
function record(ts::FallbackTestSet, t::Union{Fail,Error})
    println(t)
    throw(FallbackTestSetException("There was an error during testing"))
end
# We don't need to do anything as we don't record anything
finish(ts::FallbackTestSet) = ts

#-----------------------------------------------------------------------

"""
    DefaultTestSet

If using the DefaultTestSet, the test results will be recorded. If there
are any `Fail`s or `Error`s, an exception will be thrown only at the end,
along with a summary of the test results.
"""
mutable struct DefaultTestSet <: AbstractTestSet
    description::AbstractString
    results::Vector
    n_passed::Int
    anynonpass::Bool
end
DefaultTestSet(desc) = DefaultTestSet(desc, [], 0, false)

# For a broken result, simply store the result
record(ts::DefaultTestSet, t::Broken) = (push!(ts.results, t); t)
# For a passed result, do not store the result since it uses a lot of memory
record(ts::DefaultTestSet, t::Pass) = (ts.n_passed += 1; t)

# For the other result types, immediately print the error message
# but do not terminate. Print a backtrace.
function record(ts::DefaultTestSet, t::Union{Fail, Error})
    if myid() == 1
        print_with_color(:white, ts.description, ": ")
        print(t)
        # don't print the backtrace for Errors because it gets printed in the show
        # method
        isa(t, Error) || Base.show_backtrace(STDOUT, scrub_backtrace(backtrace()))
        println()
    end
    push!(ts.results, t)
    t, isa(t, Error) || backtrace()
end

# When a DefaultTestSet finishes, it records itself to its parent
# testset, if there is one. This allows for recursive printing of
# the results at the end of the tests
record(ts::DefaultTestSet, t::AbstractTestSet) = push!(ts.results, t)

function print_test_errors(ts::DefaultTestSet)
    for t in ts.results
        if (isa(t, Error) || isa(t, Fail)) && myid() == 1
            println("Error in testset $(ts.description):")
            Base.show(STDOUT,t)
            println()
        elseif isa(t, DefaultTestSet)
            print_test_errors(t)
        end
    end
end

function print_test_results(ts::DefaultTestSet, depth_pad=0)
    # Calculate the overall number for each type so each of
    # the test result types are aligned
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    total_pass   = passes + c_passes
    total_fail   = fails  + c_fails
    total_error  = errors + c_errors
    total_broken = broken + c_broken
    dig_pass   = total_pass   > 0 ? ndigits(total_pass)   : 0
    dig_fail   = total_fail   > 0 ? ndigits(total_fail)   : 0
    dig_error  = total_error  > 0 ? ndigits(total_error)  : 0
    dig_broken = total_broken > 0 ? ndigits(total_broken) : 0
    total = total_pass + total_fail + total_error + total_broken
    dig_total = total > 0 ? ndigits(total) : 0
    # For each category, take max of digits and header width if there are
    # tests of that type
    pass_width   = dig_pass   > 0 ? max(length("Pass"),   dig_pass)   : 0
    fail_width   = dig_fail   > 0 ? max(length("Fail"),   dig_fail)   : 0
    error_width  = dig_error  > 0 ? max(length("Error"),  dig_error)  : 0
    broken_width = dig_broken > 0 ? max(length("Broken"), dig_broken) : 0
    total_width  = dig_total  > 0 ? max(length("Total"),  dig_total)  : 0
    # Calculate the alignment of the test result counts by
    # recursively walking the tree of test sets
    align = max(get_alignment(ts, 0), length("Test Summary:"))
    # Print the outer test set header once
    pad = total == 0 ? "" : " "
    print_with_color(:white, rpad("Test Summary:",align," "), " |", pad; bold = true)
    if pass_width > 0
        print_with_color(:green, lpad("Pass",pass_width," "), "  "; bold = true)
    end
    if fail_width > 0
        print_with_color(Base.error_color(), lpad("Fail",fail_width," "), "  "; bold = true)
    end
    if error_width > 0
        print_with_color(Base.error_color(), lpad("Error",error_width," "), "  "; bold = true)
    end
    if broken_width > 0
        print_with_color(Base.warn_color(), lpad("Broken",broken_width," "), "  "; bold = true)
    end
    if total_width > 0
        print_with_color(Base.info_color(), lpad("Total",total_width, " "); bold = true)
    end
    println()
    # Recursively print a summary at every level
    print_counts(ts, depth_pad, align, pass_width, fail_width, error_width, broken_width, total_width)
end

const TESTSET_PRINT_ENABLE = Ref(true)

# Called at the end of a @testset, behaviour depends on whether
# this is a child of another testset, or the "root" testset
function finish(ts::DefaultTestSet)
    # If we are a nested test set, do not print a full summary
    # now - let the parent test set do the printing
    if get_testset_depth() != 0
        # Attach this test set to the parent test set
        parent_ts = get_testset()
        record(parent_ts, ts)
        return ts
    end
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    total_pass   = passes + c_passes
    total_fail   = fails  + c_fails
    total_error  = errors + c_errors
    total_broken = broken + c_broken
    total = total_pass + total_fail + total_error + total_broken

    if TESTSET_PRINT_ENABLE[]
        print_test_results(ts)
    end

    # Finally throw an error as we are the outermost test set
    if total != total_pass + total_broken
        # Get all the error/failures and bring them along for the ride
        efs = filter_errors(ts)
        throw(TestSetException(total_pass,total_fail,total_error, total_broken, efs))
    end

    # return the testset so it is returned from the @testset macro
    ts
end

# Recursive function that finds the column that the result counts
# can begin at by taking into account the width of the descriptions
# and the amount of indentation. If a test set had no failures, and
# no failures in child test sets, there is no need to include those
# in calculating the alignment
function get_alignment(ts::DefaultTestSet, depth::Int)
    # The minimum width at this depth is
    ts_width = 2*depth + length(ts.description)
    # If all passing, no need to look at children
    !ts.anynonpass && return ts_width
    # Return the maximum of this width and the minimum width
    # for all children (if they exist)
    isempty(ts.results) && return ts_width
    child_widths = map(t->get_alignment(t, depth+1), ts.results)
    return max(ts_width, maximum(child_widths))
end
get_alignment(ts, depth::Int) = 0

# Recursive function that fetches backtraces for any and all errors
# or failures the testset and its children encountered
function filter_errors(ts::DefaultTestSet)
    efs = []
    for t in ts.results
        if isa(t, DefaultTestSet)
            append!(efs, filter_errors(t))
        elseif isa(t, Union{Fail, Error})
            append!(efs, [t])
        end
    end
    efs
end

# Recursive function that counts the number of test results of each
# type directly in the testset, and totals across the child testsets
function get_test_counts(ts::DefaultTestSet)
    passes, fails, errors, broken = ts.n_passed, 0, 0, 0
    c_passes, c_fails, c_errors, c_broken = 0, 0, 0, 0
    for t in ts.results
        isa(t, Fail)   && (fails  += 1)
        isa(t, Error)  && (errors += 1)
        isa(t, Broken) && (broken += 1)
        if isa(t, DefaultTestSet)
            np, nf, ne, nb, ncp, ncf, nce , ncb = get_test_counts(t)
            c_passes += np + ncp
            c_fails  += nf + ncf
            c_errors += ne + nce
            c_broken += nb + ncb
        end
    end
    ts.anynonpass = (fails + errors + c_fails + c_errors > 0)
    return passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken
end

# Recursive function that prints out the results at each level of
# the tree of test sets
function print_counts(ts::DefaultTestSet, depth, align,
                      pass_width, fail_width, error_width, broken_width, total_width)
    # Count results by each type at this level, and recursively
    # through any child test sets
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken = get_test_counts(ts)
    subtotal = passes + fails + errors + broken + c_passes + c_fails + c_errors + c_broken
    # Print test set header, with an alignment that ensures all
    # the test results appear above each other
    print(rpad(string("  "^depth, ts.description), align, " "), " | ")

    np = passes + c_passes
    if np > 0
        print_with_color(:green, lpad(string(np), pass_width, " "), "  ")
    elseif pass_width > 0
        # No passes at this level, but some at another level
        print(lpad(" ", pass_width), "  ")
    end

    nf = fails + c_fails
    if nf > 0
        print_with_color(Base.error_color(), lpad(string(nf), fail_width, " "), "  ")
    elseif fail_width > 0
        # No fails at this level, but some at another level
        print(lpad(" ", fail_width), "  ")
    end

    ne = errors + c_errors
    if ne > 0
        print_with_color(Base.error_color(), lpad(string(ne), error_width, " "), "  ")
    elseif error_width > 0
        # No errors at this level, but some at another level
        print(lpad(" ", error_width), "  ")
    end

    nb = broken + c_broken
    if nb > 0
        print_with_color(Base.warn_color(), lpad(string(nb), broken_width, " "), "  ")
    elseif broken_width > 0
        # None broken at this level, but some at another level
        print(lpad(" ", broken_width), "  ")
    end

    if np == 0 && nf == 0 && ne == 0 && nb == 0
        print_with_color(Base.info_color(), "No tests")
    else
        print_with_color(Base.info_color(), lpad(string(subtotal), total_width, " "))
    end
    println()

    # Only print results at lower levels if we had failures
    if np + nb != subtotal
        for t in ts.results
            if isa(t, DefaultTestSet)
                print_counts(t, depth + 1, align,
                    pass_width, fail_width, error_width, broken_width, total_width)
            end
        end
    end
end

#-----------------------------------------------------------------------

"""
    @testset [CustomTestSet] [option=val  ...] ["description"] begin ... end
    @testset [CustomTestSet] [option=val  ...] ["description \$v"] for v in (...) ... end
    @testset [CustomTestSet] [option=val  ...] ["description \$v, \$w"] for v in (...), w in (...) ... end

Starts a new test set, or multiple test sets if a `for` loop is provided.

If no custom testset type is given it defaults to creating a `DefaultTestSet`.
`DefaultTestSet` records all the results and, if there are any `Fail`s or
`Error`s, throws an exception at the end of the top-level (non-nested) test set,
along with a summary of the test results.

Any custom testset type (subtype of `AbstractTestSet`) can be given and it will
also be used for any nested `@testset` invocations. The given options are only
applied to the test set where they are given. The default test set type does
not take any options.

The description string accepts interpolation from the loop indices.
If no description is provided, one is constructed based on the variables.

By default the `@testset` macro will return the testset object itself, though
this behavior can be customized in other testset types. If a `for` loop is used
then the macro collects and returns a list of the return values of the `finish`
method, which by default will return a list of the testset objects used in
each iteration.
"""
macro testset(args...)
    isempty(args) && error("No arguments to @testset")

    tests = args[end]

    # Determine if a single block or for-loop style
    if !isa(tests,Expr) || (tests.head != :for && tests.head != :block)
        error("Expected begin/end block or for loop as argument to @testset")
    end

    if tests.head == :for
        return testset_forloop(args, tests)
    else
        return testset_beginend(args, tests)
    end
end

"""
Generate the code for a `@testset` with a `begin`/`end` argument
"""
function testset_beginend(args, tests)
    desc, testsettype, options = parse_testset_args(args[1:end-1])
    if desc === nothing
        desc = "test set"
    end
    # If we're at the top level we'll default to DefaultTestSet. Otherwise
    # default to the type of the parent testset
    if testsettype === nothing
        testsettype = :(get_testset_depth() == 0 ? DefaultTestSet : typeof(get_testset()))
    end

    # Generate a block of code that initializes a new testset, adds
    # it to the task local storage, evaluates the test(s), before
    # finally removing the testset and giving it a chance to take
    # action (such as reporting the results)
    quote
        ts = $(testsettype)($desc; $options...)
        # this empty loop is here to force the block to be compiled,
        # which is needed for backtrace scrubbing to work correctly.
        while false; end
        push_testset(ts)
        try
            $(esc(tests))
        catch err
            # something in the test block threw an error. Count that as an
            # error in this test set
            record(ts, Error(:nontest_error, :(), err, catch_backtrace()))
        end
        pop_testset()
        finish(ts)
    end
end

"""
Generate the code for a `@testset` with a `for` loop argument
"""
function testset_forloop(args, testloop)
    # Pull out the loop variables. We might need them for generating the
    # description and we'll definitely need them for generating the
    # comprehension expression at the end
    loopvars = Expr[]
    if testloop.args[1].head == :(=)
        push!(loopvars, testloop.args[1])
    elseif testloop.args[1].head == :block
        for loopvar in testloop.args[1].args
            push!(loopvars, loopvar)
        end
    else
        error("Unexpected argument to @testset")
    end

    desc, testsettype, options = parse_testset_args(args[1:end-1])

    if desc === nothing
        # No description provided. Generate from the loop variable names
        v = loopvars[1].args[1]
        desc = Expr(:string, "$v = ", esc(v)) # first variable
        for l = loopvars[2:end]
            v = l.args[1]
            push!(desc.args, ", $v = ")
            push!(desc.args, esc(v))
        end
    end

    if testsettype === nothing
        testsettype = :(get_testset_depth() == 0 ? DefaultTestSet : typeof(get_testset()))
    end

    # Uses a similar block as for `@testset`, except that it is
    # wrapped in the outer loop provided by the user
    tests = testloop.args[2]
    blk = quote
        # Trick to handle `break` and `continue` in the test code before
        # they can be handled properly by `finally` lowering.
        if !first_iteration
            pop_testset()
            push!(arr, finish(ts))
        end
        ts = $(testsettype)($desc; $options...)
        push_testset(ts)
        first_iteration = false
        try
            $(esc(tests))
        catch err
            # Something in the test block threw an error. Count that as an
            # error in this test set
            record(ts, Error(:nontest_error, :(), err, catch_backtrace()))
        end
    end
    quote
        arr = Array{Any,1}(0)
        first_iteration = true
        local ts
        try
            $(Expr(:for, Expr(:block, [esc(v) for v in loopvars]...), blk))
        finally
            # Handle `return` in test body
            if !first_iteration
                pop_testset()
                push!(arr, finish(ts))
            end
        end
        arr
    end
end

"""
Parse the arguments to the `@testset` macro to pull out the description,
Testset Type, and options. Generally this should be called with all the macro
arguments except the last one, which is the test expression itself.
"""
function parse_testset_args(args)
    desc = nothing
    testsettype = nothing
    options = :(Dict{Symbol, Any}())
    for arg in args
        # a standalone symbol is assumed to be the test set we should use
        if isa(arg, Symbol)
            testsettype = esc(arg)
        # a string is the description
        elseif isa(arg, AbstractString) || (isa(arg, Expr) && arg.head == :string)
            desc = esc(arg)
        # an assignment is an option
        elseif isa(arg, Expr) && arg.head == :(=)
            # we're building up a Dict literal here
            key = Expr(:quote, arg.args[1])
            push!(options.args, Expr(:call, :(=>), key, arg.args[2]))
        else
            error("Unexpected argument $arg to @testset")
        end
    end

    (desc, testsettype, options)
end

#-----------------------------------------------------------------------
# Various helper methods for test sets

"""
    get_testset()

Retrieve the active test set from the task's local storage. If no
test set is active, use the fallback default test set.
"""
function get_testset()
    testsets = get(task_local_storage(), :__BASETESTNEXT__, AbstractTestSet[])
    return isempty(testsets) ? fallback_testset : testsets[end]
end

"""
    push_testset(ts::AbstractTestSet)

Adds the test set to the task_local_storage.
"""
function push_testset(ts::AbstractTestSet)
    testsets = get(task_local_storage(), :__BASETESTNEXT__, AbstractTestSet[])
    push!(testsets, ts)
    setindex!(task_local_storage(), testsets, :__BASETESTNEXT__)
end

"""
    pop_testset()

Pops the last test set added to the task_local_storage. If there are no
active test sets, returns the fallback default test set.
"""
function pop_testset()
    testsets = get(task_local_storage(), :__BASETESTNEXT__, AbstractTestSet[])
    ret = isempty(testsets) ? fallback_testset : pop!(testsets)
    setindex!(task_local_storage(), testsets, :__BASETESTNEXT__)
    return ret
end

"""
    get_testset_depth()

Returns the number of active test sets, not including the defaut test set
"""
function get_testset_depth()
    testsets = get(task_local_storage(), :__BASETESTNEXT__, AbstractTestSet[])
    return length(testsets)
end

_args_and_call(args...; kwargs...) = (args[1:end-1], kwargs, args[end](args[1:end-1]...; kwargs...))
"""
    @inferred f(x)

Tests that the call expression `f(x)` returns a value of the same type
inferred by the compiler. It is useful to check for type stability.

`f(x)` can be any call expression.
Returns the result of `f(x)` if the types match,
and an `Error` `Result` if it finds different types.

```jldoctest
julia> using Base.Test

julia> f(a,b,c) = b > 1 ? 1 : 1.0
f (generic function with 1 method)

julia> typeof(f(1,2,3))
Int64

julia> @code_warntype f(1,2,3)
Variables:
  #self#::#f
  a::Int64
  b::Int64
  c::Int64

Body:
  begin
      unless (Base.slt_int)(1, b::Int64)::Bool goto 3
      return 1
      3:
      return 1.0
  end::UNION{FLOAT64, INT64}

julia> @inferred f(1,2,3)
ERROR: return type Int64 does not match inferred return type Union{Float64, Int64}
Stacktrace:
 [1] error(::String) at ./error.jl:21

julia> @inferred max(1,2)
2
```
"""
macro inferred(ex)
    if Meta.isexpr(ex, :ref)
        ex = Expr(:call, :getindex, ex.args...)
    end
    Meta.isexpr(ex, :call)|| error("@inferred requires a call expression")

    Base.remove_linenums!(quote
        let
            $(if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex.args)
                # Has keywords
                args = gensym()
                kwargs = gensym()
                quote
                    $(esc(args)), $(esc(kwargs)), result = $(esc(Expr(:call, _args_and_call, ex.args[2:end]..., ex.args[1])))
                    inftypes = $(Base.gen_call_with_extracted_types(__module__, Base.return_types, :($(ex.args[1])($(args)...; $(kwargs)...))))
                end
            else
                # No keywords
                quote
                    args = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
                    result = $(esc(ex.args[1]))(args...)
                    inftypes = Base.return_types($(esc(ex.args[1])), Base.typesof(args...))
                end
            end)
            @assert length(inftypes) == 1
            rettype = isa(result, Type) ? Type{result} : typeof(result)
            rettype == inftypes[1] || error("return type $rettype does not match inferred return type $(inftypes[1])")
            result
        end
    end)
end

# Test approximate equality of vectors or columns of matrices modulo floating
# point roundoff and phase (sign) differences.
#
# This function is designed to test for equality between vectors of floating point
# numbers when the vectors are defined only up to a global phase or sign, such as
# normalized eigenvectors or singular vectors. The global phase is usually
# defined consistently, but may occasionally change due to small differences in
# floating point rounding noise or rounding modes, or through the use of
# different conventions in different algorithms. As a result, most tests checking
# such vectors have to detect and discard such overall phase differences.
#
# Inputs:
#     a, b:: StridedVecOrMat to be compared
#     err :: Default: m^3*(eps(S)+eps(T)), where m is the number of rows
#
# Raises an error if any columnwise vector norm exceeds err. Otherwise, returns
# nothing.
function test_approx_eq_modphase(a::StridedVecOrMat{S}, b::StridedVecOrMat{T},
                                 err = length(indices(a,1))^3*(eps(S)+eps(T))) where {S<:Real,T<:Real}
    @test indices(a,1) == indices(b,1) && indices(a,2) == indices(b,2)
    for i in indices(a,2)
        v1, v2 = a[:, i], b[:, i]
        @test min(abs(norm(v1-v2)),abs(norm(v1+v2))) ≈ 0.0 atol=err
    end
end

"""
    detect_ambiguities(mod1, mod2...; imported=false, ambiguous_bottom=false)

Returns a vector of `(Method,Method)` pairs of ambiguous methods
defined in the specified modules. Use `imported=true` if you wish to
also test functions that were imported into these modules from
elsewhere.

`ambiguous_bottom` controls whether ambiguities triggered only by
`Union{}` type parameters are included; in most cases you probably
want to set this to `false`. See [`Base.isambiguous`](@ref).
"""
function detect_ambiguities(mods...;
                            imported::Bool = false,
                            ambiguous_bottom::Bool = false,
                            allow_bottom::Union{Bool,Void} = nothing)
    if allow_bottom !== nothing
        Base.depwarn("the `allow_bottom` keyword to detect_ambiguities has been renamed to `ambiguous_bottom`", :detect_ambiguities)
        ambiguous_bottom = allow_bottom
    end
    function sortdefs(m1, m2)
        ord12 = m1.file < m2.file
        if !ord12 && (m1.file == m2.file)
            ord12 = m1.line < m2.line
        end
        ord12 ? (m1, m2) : (m2, m1)
    end
    ambs = Set{Tuple{Method,Method}}()
    for mod in mods
        for n in names(mod, true, imported)
            Base.isdeprecated(mod, n) && continue
            if !isdefined(mod, n)
                println("Skipping ", mod, '.', n)  # typically stale exports
                continue
            end
            f = getfield(mod, n)
            if isa(f, DataType) && isdefined(f.name, :mt)
                mt = Base.MethodList(f.name.mt)
                for m in mt
                    if m.ambig !== nothing
                        for m2 in m.ambig
                            if Base.isambiguous(m, m2, ambiguous_bottom=ambiguous_bottom)
                                push!(ambs, sortdefs(m, m2))
                            end
                        end
                    end
                end
            end
        end
    end
    return collect(ambs)
end

"""
The `GenericString` can be used to test generic string APIs that program to
the `AbstractString` interface, in order to ensure that functions can work
with string types besides the standard `String` type.
"""
struct GenericString <: AbstractString
    string::AbstractString
end
Base.convert(::Type{GenericString}, s::AbstractString) = GenericString(s)
Base.endof(s::GenericString) = endof(s.string)
Base.next(s::GenericString, i::Int) = next(s.string, i)

"""
The `GenericSet` can be used to test generic set APIs that program to
the `AbstractSet` interface, in order to ensure that functions can work
with set types besides the standard `Set` and `IntSet` types.
"""
struct GenericSet{T} <: AbstractSet{T}
    s::AbstractSet{T}
end

"""
The `GenericDict` can be used to test generic dict APIs that program to
the `Associative` interface, in order to ensure that functions can work
with associative types besides the standard `Dict` type.
"""
struct GenericDict{K,V} <: Associative{K,V}
    s::Associative{K,V}
end

for (G, A) in ((GenericSet, AbstractSet),
               (GenericDict, Associative))
    @eval begin
        Base.convert(::Type{$G}, s::$A) = $G(s)
        Base.done(s::$G, state) = done(s.s, state)
        Base.next(s::$G, state) = next(s.s, state)
    end
    for f in (:eltype, :isempty, :length, :start)
        @eval begin
            Base.$f(s::$G) = $f(s.s)
        end
    end
end

Base.get(s::GenericDict, x, y) = get(s.s, x, y)

end # module

ext: jl

max_stars_repo_name: nassarhuda/julia

max_issues_repo_name: nassarhuda/julia

max_forks_repo_name: nassarhuda/julia


text:
using WordTokenizers
using Zygote

# WIP


text:
const allocated = Dict{Tuple{ZeContext,ZePtr},Tuple{oneL0.DeviceBuffer,Int}}()

function allocate(ctx, dev, bytes::Int, alignment::Int)
    # 0-byte allocations shouldn't hit the pool
    bytes == 0 && return ZE_NULL

    buf = device_alloc(ctx, dev, bytes, alignment)
    make_resident(ctx, dev, buf)

    ptr = convert(ZePtr{Nothing}, buf)
    @assert !haskey(allocated, (ctx,ptr))
    allocated[(ctx,ptr)] = buf, 1

    return buf
end

function alias(ctx, dev, ptr)
    # 0-byte allocations shouldn't hit the pool
    ptr == ZE_NULL && return

    buf, refcount = allocated[(ctx,ptr)]
    allocated[(ctx,ptr)] = buf, refcount+1

    return
end

function release(ctx, dev, ptr)
    # 0-byte allocations shouldn't hit the pool
    ptr == ZE_NULL && return

    buf, refcount = allocated[(ctx,ptr)]
    if refcount == 1
        delete!(allocated, (ctx,ptr))
    else
        allocated[(ctx,ptr)] = buf, refcount-1
        return
    end

    evict(ctx, dev, buf)
    free(buf)

    return
end


text:
"""
Stores results data for one EmulationModel
"""
mutable struct EmulationModelStore <: AbstractModelStore
    data_container::DatasetContainer{DataFrameDataset}
    optimizer_stats::OrderedDict{Int, OptimizerStats}
end

get_data_field(store::EmulationModelStore, type::Symbol) =
    getfield(store.data_container, type)

function EmulationModelStore()
    return EmulationModelStore(
        DatasetContainer{DataFrameDataset}(),
        OrderedDict{Int, OptimizerStats}(),
    )
end

function Base.empty!(store::EmulationModelStore)
    stype = DatasetContainer
    for (name, _) in zip(fieldnames(stype), fieldtypes(stype))
        if name ∉ [:values, :timestamps]
            val = get_data_field(store, name)
            try
                empty!(val)
            catch
                @error "Base.empty! must be customized for type $stype or skipped"
                rethrow()
            end
        elseif name == :update_timestamp
            store.update_timestamp = UNSET_INI_TIME
        else
            setfield!(
                store.data_container,
                name,
                zero(fieldtype(store.data_container, name)),
            )
        end
    end
    empty!(store.optimizer_stats)
    return
end

function Base.isempty(store::EmulationModelStore)
    stype = DatasetContainer
    for (name, type) in zip(fieldnames(stype), fieldtypes(stype))
        if name ∉ [:values, :timestamps]
            val = get_data_field(store, name)
            try
                !isempty(val) && return false
            catch
                @error "Base.isempty must be customized for type $stype or skipped"
                rethrow()
            end
        elseif name == :update_timestamp
            store.update_timestamp != UNSET_INI_TIME && return false
        else
            val = get_data_fieldd(store, name)
            iszero(val) && return false
        end
    end
    return isempty(store.optimizer_stats)
end

function initialize_storage!(
    store::EmulationModelStore,
    container::OptimizationContainer,
    params::ModelStoreParams,
)
    num_of_executions = get_num_executions(params)
    for type in STORE_CONTAINERS
        field_containers = getfield(container, type)
        results_container = get_data_field(store, type)
        for (key, field_container) in field_containers
            @debug "Adding $(encode_key_as_string(key)) to EmulationModelStore" _group =
                LOG_GROUP_MODEL_STORE
            column_names = get_column_names(key, field_container)
            results_container[key] = DataFrameDataset(
                DataFrames.DataFrame(
                    OrderedDict(c => fill(NaN, num_of_executions) for c in column_names),
                ),
            )
        end
    end
    return
end

function write_result!(
    store::EmulationModelStore,
    name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    update_timestamp::Dates.DateTime,
    array::AbstractArray,
)
    df = axis_array_to_dataframe(array, key)
    write_result!(store, name, key, index, update_timestamp, df)
    return
end

function write_result!(
    store::EmulationModelStore,
    name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    update_timestamp::Dates.DateTime,
    df::DataFrames.DataFrame,
)
    @assert_op size(df)[1] == 1
    write_result!(store, name, key, index, update_timestamp, df[1, :])
    return
end

function write_result!(
    store::EmulationModelStore,
    ::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    update_timestamp::Dates.DateTime,
    df_row::DataFrames.DataFrameRow,
)
    container = get_data_field(store, get_store_container_type(key))
    set_value!(container[key], df_row, index)
    set_last_recorded_row!(container[key], index)
    set_update_timestamp!(container[key], update_timestamp)
    return
end

function read_results(
    store::EmulationModelStore,
    key::OptimizationContainerKey;
    index::Union{Int, Nothing}=nothing,
    len::Union{Int, Nothing}=nothing,
)
    container = get_data_field(store, get_store_container_type(key))
    df = container[key].values
    # Return a copy because callers may mutate it.
    if isnothing(index)
        @assert_op len === nothing
        return copy(df, copycols=true)
    elseif isnothing(len)
        return copy(df, copycols=true)[index:end, :]
    else
        return copy(df, copycols=true)[index:(index + len - 1), :]
    end
end

function get_dataset_size(store::EmulationModelStore, key::OptimizationContainerKey)
    container = get_data_field(store, get_store_container_type(key))
    return size(container[key].values)
end

function get_last_updated_timestamp(
    store::EmulationModelStore,
    key::OptimizationContainerKey,
)
    container = get_data_field(store, get_store_container_type(key))
    return get_update_timestamp(container[key])
end
function write_optimizer_stats!(
    store::EmulationModelStore,
    stats::OptimizerStats,
    index::EmulationModelIndexType,
)
    @assert !(index in keys(store.optimizer_stats))
    store.optimizer_stats[index] = stats
    return
end

function read_optimizer_stats(store::EmulationModelStore)
    return DataFrames.DataFrame([to_namedtuple(x) for x in values(store.optimizer_stats)])
end

function get_last_recorded_row(x::EmulationModelStore, key::OptimizationContainerKey)
    return get_last_recorded_row(x.data_container, key)
end


text:
export iauEors
"""
Equation of the origins, given the classical NPB matrix and the
quantity s.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rnpb  double[3][3]  classical nutation x precession x bias matrix
   s     double        the quantity s (the CIO locator)

Returned (function value):
         double        the equation of the origins in radians.

Notes:

   1.  The equation of the origins is the distance between the true
   equinox and the celestial intermediate origin and, equivalently,
   the difference between Earth rotation angle and Greenwich
   apparent sidereal time (ERA-GST).  It comprises the precession
   (since J2000.0) in right ascension plus the equation of the
   equinoxes (including the small correction terms).

   2.  The algorithm is from Wallace & Capitaine (2006).

References:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   Wallace, P. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauEors(double rnpb[3][3], double s)

function iauEors(rnbp::AbstractMatrix{<:Real}, s::Real)
   # Allocate return value
   rm = zeros(Float64, 3, 3)

   return ccall((:iauEors, libsofa_c), Cdouble,
               (Ptr{Cdouble}, Cdouble),
               convert(Array{Float64, 2}, rnbp'), convert(Float64, s))
end


text:
# Function for estimating a problem's time step
function update_a2a3adot(surf::TwoDSurf,dt)
    for ia = 2:3
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:3
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

function update_atermdot(surf::TwoDSurf,dt)
    for ia = 2:surf.naterm
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:surf.naterm
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

function update_adot(surf::TwoDSurf,dt)
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:3
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

# Function for updating the induced velocities
function update_indbound(surf::TwoDSurf, curfield::TwoDFlowField)
    surf.uind[1:surf.ndiv], surf.wind[1:surf.ndiv] = ind_vel([curfield.tev; curfield.lev; curfield.extv], surf.bnd_x, surf.bnd_z)
    return surf
end

function add_indbound_b(surf::TwoDSurf, surfj::TwoDSurf)
    uind, wind = ind_vel(surfj.bv, surf.bnd_x, surf.bnd_z)
    surf.uind[:] += uind
    surf.wind[:] += wind
    return surf

end

# Function for updating the downwash
function update_downwash(surf::TwoDSurf, vels::Vector{Float64})
    for ib = 1:surf.ndiv
        surf.downwash[ib] = -(surf.kinem.u + vels[1])*sin(surf.kinem.alpha) - surf.uind[ib]*sin(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*cos(surf.kinem.alpha) - surf.wind[ib]*cos(surf.kinem.alpha) - surf.kinem.alphadot*(surf.x[ib] - surf.pvt*surf.c) + surf.cam_slope[ib]*(surf.uind[ib]*cos(surf.kinem.alpha) + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))
    end
    return surf
end

# Function for a_0 and a_1 fourier coefficients
function update_a0anda1(surf::TwoDSurf)
    surf.a0[1] = simpleTrapz(surf.downwash,surf.theta)
    surf.aterm[1] = simpleTrapz(surf.downwash.*cos.(surf.theta),surf.theta)
    surf.a0[1] = -surf.a0[1]/(surf.uref*pi)
    surf.aterm[1] = 2. *surf.aterm[1]/(surf.uref*pi)
    return surf
end

# Function for calculating the fourier coefficients a_2 upwards to a_n
function update_a2toan(surf::TwoDSurf)
    for ia = 2:surf.naterm
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    return surf
end

# Function to update the external flowfield
function update_externalvel(curfield::TwoDFlowField, t)
    if (typeof(curfield.velX) == CosDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    elseif (typeof(curfield.velX) == SinDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    elseif (typeof(curfield.velX) == ConstDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    end
    if (typeof(curfield.velX) == StepGustDef)
        curfield.u[1] = curfield.velX(t)
    end
    if (typeof(curfield.velZ) == StepGustDef)
        curfield.w[1] = curfield.velZ(t)
    end
end

# Function updating the dimensional kinematic parameters
function update_kinem(surf::TwoDSurf, t)

    # Pitch kinematics
    if (typeof(surf.kindef.alpha) == EldUpDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldUptstartDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldRampReturnDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldRampReturntstartDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == ConstDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = 0.
    elseif (typeof(surf.kindef.alpha) == SinDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == CosDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == FileDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == VAWTalphaDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------

    # Plunge kinematics
    if (typeof(surf.kindef.h) == EldUpDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUptstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpIntDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpInttstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldRampReturnDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldRampReturntstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == ConstDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = 0.
    elseif (typeof(surf.kindef.h) == SinDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == CosDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == FileDef)
        surf.kinem.h = surf.kindef.h(t)
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == VAWThDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    end
    # ---------------------------------------------------------------------------------------------

    # Forward velocity
    if (typeof(surf.kindef.u) == EldUpDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == EldRampReturnDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == EldRampReturntstartDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == ConstDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = 0.
    elseif (typeof(surf.kindef.u) == SinDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == CosDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == LinearDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == FileDef)
        surf.kinem.u = surf.kindef.u(t)
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == VAWTuDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------
    return surf
end

# ---------------------------------------------------------------------------------------------
# Updates the bound vorticity distribution: eqn (2.1) in Ramesh et al. (2013)
# determines the strength of bound vortices
# determines the x, z components of the bound vortices
function update_bv(surf::TwoDSurf)
    gamma = zeros(surf.ndiv)
    for ib = 1:surf.ndiv
        gamma[ib] = (surf.a0[1]*(1 + cos(surf.theta[ib])))
        for ia = 1:surf.naterm
            gamma[ib] = gamma[ib] + surf.aterm[ia]*sin(ia*surf.theta[ib])*sin(surf.theta[ib])
        end
        gamma[ib] = gamma[ib]*surf.uref*surf.c
    end

    for ib = 2:surf.ndiv
        surf.bv[ib-1].s = (gamma[ib]+gamma[ib-1])*(surf.theta[2]-surf.theta[1])/2.
        surf.bv[ib-1].x = (surf.bnd_x[ib] + surf.bnd_x[ib-1])/2.
        surf.bv[ib-1].z = (surf.bnd_z[ib] + surf.bnd_z[ib-1])/2.
    end
end

# Function for calculating the wake rollup
function wakeroll(surf::TwoDSurf, curfield::TwoDFlowField, dt)

    nlev = length(curfield.lev)
    ntev = length(curfield.tev)
    nextv = length(curfield.extv)

    #Clean induced velocities
    for i = 1:ntev
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:nlev
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    for i = 1:nextv
        curfield.extv[i].vx = 0
        curfield.extv[i].vz = 0
    end

    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev; curfield.extv])

    #Add the influence of velocities induced by bound vortices
    utemp = zeros(ntev + nlev + nextv)
    wtemp = zeros(ntev + nlev + nextv)
    utemp, wtemp = ind_vel(surf.bv, [map(q -> q.x, curfield.tev); map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)], [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev); map(q -> q.z, curfield.extv) ])

    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Add the influence of freestream velocities
    for i = 1:ntev
        curfield.tev[i].vx += curfield.u[1]
        curfield.tev[i].vz += curfield.w[1]
    end
    for i = 1:nlev
        curfield.lev[i].vx += curfield.u[1]
        curfield.lev[i].vz += curfield.w[1]
    end
    for i = 1:nextv
        curfield.extv[i].vx += curfield.u[1]
        curfield.extv[i].vz += curfield.w[1]
    end

    #Convect free vortices with their induced velocities
    for i = 1:ntev
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:nlev
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:nextv
        curfield.extv[i].x += dt*curfield.extv[i].vx
        curfield.extv[i].z += dt*curfield.extv[i].vz
    end

    return curfield
end

function wakeroll(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, dt)

    nlev = length(curfield.lev)
    ntev = length(curfield.tev)
    nextv = length(curfield.extv)
    nsurf = length(surf)

    #Clean induced velocities
    for i = 1:ntev
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:nlev
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    for i = 1:nextv
        curfield.extv[i].vx = 0
        curfield.extv[i].vz = 0
    end

    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev; curfield.extv])

    #Add the influence of velocities induced by bound vortices
    utemp = zeros(ntev + nlev + nextv)
    wtemp = zeros(ntev + nlev + nextv)
    for is = 1:nsurf
        ut, wt = ind_vel(surf[is].bv, [map(q -> q.x, curfield.tev); map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)], [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev); map(q -> q.z, curfield.extv) ])
        utemp += ut
        wtemp += wt
    end
    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Add the influence of freestream velocities
    for i = 1:ntev
        curfield.tev[i].vx += curfield.u[1]
        curfield.tev[i].vz += curfield.w[1]
    end
    for i = 1:nlev
        curfield.lev[i].vx += curfield.u[1]
        curfield.lev[i].vz += curfield.w[1]
    end
    for i = 1:nextv
        curfield.extv[i].vx += curfield.u[1]
        curfield.extv[i].vz += curfield.w[1]
    end

    #Convect free vortices with their induced velocities
    for i = 1:ntev
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:nlev
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:nextv
        curfield.extv[i].x += dt*curfield.extv[i].vx
        curfield.extv[i].z += dt*curfield.extv[i].vz
    end

    return curfield
end

# Places a trailing edge vortex
function place_tev(surf::TwoDSurf,field::TwoDFlowField,dt)
    ntev = length(field.tev)
    if ntev == 0
        xloc = surf.bnd_x[surf.ndiv] + 0.5*surf.kinem.u*dt
        zloc = surf.bnd_z[surf.ndiv]
    else
        xloc = surf.bnd_x[surf.ndiv]+(1. /3.)*(field.tev[ntev].x - surf.bnd_x[surf.ndiv])

        zloc = surf.bnd_z[surf.ndiv]+(1. /3.)*(field.tev[ntev].z - surf.bnd_z[surf.ndiv])
    end
    push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))
    return field
end

function place_tev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt)
    nsurf = length(surf)
    ntev = length(field.tev)
    if ntev < nsurf
        for i = 1:nsurf
            xloc = surf[i].bnd_x[end] + 0.5*surf[i].kinem.u*dt
            zloc = surf[i].bnd_z[end]
            push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
        end
    else
        for i = 1:nsurf
            xloc = surf[i].bnd_x[end]+(1. /3.)*(field.tev[ntev-nsurf+i].x - surf[i].bnd_x[end])
            zloc = surf[i].bnd_z[end]+(1. /3.)*(field.tev[ntev-nsurf+i].z - surf[i].bnd_z[end])
            push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
        end
    end
    return field
end

# Places a leading edge vortex
function place_lev(surf::TwoDSurf,field::TwoDFlowField,dt)
    nlev = length(field.lev)

    le_vel_x = surf.kinem.u - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1]
    le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c- surf.kinem.hdot + surf.wind[1]

    if (surf.levflag[1] == 0)
        xloc = surf.bnd_x[1] + 0.5*le_vel_x*dt
        zloc = surf.bnd_z[1] + 0.5*le_vel_z*dt
    else
        xloc = surf.bnd_x[1]+(1. /3.)*(field.lev[nlev].x - surf.bnd_x[1])
        zloc = surf.bnd_z[1]+(1. /3.)*(field.lev[nlev].z - surf.bnd_z[1])
    end

    push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))

    return field
end

function place_lev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt,shed_ind::Vector{Int})
    nsurf = length(surf)
    nlev = length(field.lev)

    for i = 1:nsurf
        if i in shed_ind
            if (surf[i].levflag[1] == 0)
                le_vel_x = surf[i].kinem.u - surf[i].kinem.alphadot*sin(surf[i].kinem.alpha)*surf[i].pvt*surf[i].c + surf[i].uind[1]
                le_vel_z = -surf[i].kinem.alphadot*cos(surf[i].kinem.alpha)*surf[i].pvt*surf[i].c- surf[i].kinem.hdot + surf[i].wind[1]
                xloc = surf[i].bnd_x[1] + 0.5*le_vel_x*dt
                zloc = surf[i].bnd_z[1] + 0.5*le_vel_z*dt
                push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
            else
                xloc = surf[i].bnd_x[1]+(1. /3.)*(field.lev[nlev-nsurf+i].x - surf[i].bnd_x[1])
                zloc = surf[i].bnd_z[1]+(1. /3.)*(field.lev[nlev-nsurf+i].z - surf[i].bnd_z[1])
                push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
            end
        else
            push!(field.lev, TwoDVort(0., 0., 0., 0., 0., 0.))
        end
    end
    return field
end

# Function for updating the positions of the bound vortices
function update_boundpos(surf::TwoDSurf, dt::Float64)
    for i = 1:surf.ndiv
        surf.bnd_x[i] = surf.bnd_x[i] + dt*((surf.pvt*surf.c - surf.x[i])*sin(surf.kinem.alpha)*surf.kinem.alphadot - surf.kinem.u + surf.cam[i]*cos(surf.kinem.alpha)*surf.kinem.alphadot)
        surf.bnd_z[i] = surf.bnd_z[i] + dt*(surf.kinem.hdot + (surf.pvt*surf.c - surf.x[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot - surf.cam[i]*sin(surf.kinem.alpha)*surf.kinem.alphadot)
    end
    return surf
end

# Function to calculate induced velocities by a set of vortices at a target location
function ind_vel(src::Vector{TwoDVort},t_x,t_z)
    #'s' stands for source and 't' for target
    uind = zeros(length(t_x))
    wind = zeros(length(t_x))

    for itr = 1:length(t_x)
        for isr = 1:length(src)
            xdist = src[isr].x - t_x[itr]
            zdist = src[isr].z - t_z[itr]
            distsq = xdist*xdist + zdist*zdist
            uind[itr] = uind[itr] - src[isr].s*zdist/(2*pi*sqrt(src[isr].vc*src[isr].vc*src[isr].vc*src[isr].vc + distsq*distsq))
            wind[itr] = wind[itr] + src[isr].s*xdist/(2*pi*sqrt(src[isr].vc*src[isr].vc*src[isr].vc*src[isr].vc + distsq*distsq))
        end
    end
    return uind, wind
end
# Same as above except allows for only one vortex and multiple collocation points
#   Used for influence coefficients
function ind_vel(src::TwoDVort,t_x,t_z)
    #'s' stands for source and 't' for target
    xdist = src.x .- t_x
    zdist = src.z .- t_z
    distsq = xdist.*xdist + zdist.*zdist
    uind = -src.s.*zdist./(2*pi*sqrt.(src.vc*src.vc*src.vc*src.vc .+ distsq.*distsq))
    wind = src.s.*xdist./(2*pi*sqrt.(src.vc*src.vc*src.vc*src.vc .+ distsq.*distsq))
    return uind, wind
end

# Function determining the effects of interacting vortices - velocities induced on each other - classical n-body problem
function mutual_ind(vorts::Vector{TwoDVort})
    for i = 1:length(vorts)
        for j = i+1:length(vorts)
            dx = vorts[i].x - vorts[j].x
            dz = vorts[i].z - vorts[j].z
            #source- tar
            dsq = dx*dx + dz*dz

            magitr = 1. /(2*pi*sqrt(vorts[j].vc*vorts[j].vc*vorts[j].vc*vorts[j].vc + dsq*dsq))
            magjtr = 1. /(2*pi*sqrt(vorts[i].vc*vorts[i].vc*vorts[i].vc*vorts[i].vc + dsq*dsq))

            vorts[j].vx -= dz * vorts[i].s * magjtr
            vorts[j].vz += dx * vorts[i].s * magjtr

            vorts[i].vx += dz * vorts[j].s * magitr
            vorts[i].vz -= dx * vorts[j].s * magitr
        end
    end
    return vorts
end

"""
    controlVortCount(delvort, surf, curfield)

Performs merging or deletion operation on free vortices in order to
control computational cost according to algorithm specified.

Algorithms for parameter delvort

 - delNone

    Does nothing, no vortex count control.

 - delSpalart(limit=500, dist=10, tol=1e-6)

    Merges vortices according to algorithm given in Spalart,
    P. R. (1988). Vortex methods for separated flows.

     - limit: min number of vortices present for merging to occur

     - dist: small values encourage mergin near airfoil, large values in
         wake (see paper)

     - tol: tolerance for merging, merging is less likely to occur for low
        values (see paper)

There is no universal set of parameters that work for all problem. If
using vortex control, test and calibrate choice of parameters

"""
function controlVortCount(delvort :: delVortDef, surf_locx :: Float64, surf_locz :: Float64, curfield :: TwoDFlowField)

    if typeof(delvort) == delNone

    elseif typeof(delvort) == delSpalart
        D0 = delvort.dist
        V0 = delvort.tol
        if length(curfield.tev) > delvort.limit
            #Check possibility of merging the last vortex with the closest 20 vortices
            gamma_j = curfield.tev[1].s
            d_j = sqrt((curfield.tev[1].x- surf_locx)^2 + (curfield.tev[1].z - surf_locz)^2)
            z_j = sqrt(curfield.tev[1].x^2 + curfield.tev[1].z^2)

            for i = 2:20
                gamma_k = curfield.tev[i].s
                d_k = sqrt((curfield.tev[i].x - surf_locx)^2 + (curfield.tev[i].z - surf_locz)^2)
                z_k = sqrt(curfield.tev[i].x^2 + curfield.tev[i].z^2)

                fact = abs(gamma_j*gamma_k)*abs(z_j - z_k)/(abs(gamma_j + gamma_k)*(D0 + d_j)^1.5*(D0 + d_k)^1.5)

                if fact < V0
                    #Merge the 2 vortices into the one at k
                    curfield.tev[i].x = (abs(gamma_j)*curfield.tev[1].x + abs(gamma_k)*curfield.tev[i].x)/(abs(gamma_j + gamma_k))
                    curfield.tev[i].z = (abs(gamma_j)*curfield.tev[1].z + abs(gamma_k)*curfield.tev[i].z)/(abs(gamma_j + gamma_k))
                    curfield.tev[i].s += curfield.tev[1].s

                    popfirst!(curfield.tev)

                    break
                end
            end

            if length(curfield.lev) > delvort.limit
                #Check possibility of merging the last vortex with the closest 20 vortices

                gamma_j = curfield.lev[1].s
                d_j = sqrt((curfield.lev[1].x - surf_locx)^2 + (curfield.lev[1].z - surf_locz)^2)
                z_j = sqrt(curfield.lev[1].x^2 + curfield.lev[1].z^2)

                for i = 2:20
                    gamma_k = curfield.lev[i].s
                    d_k = sqrt((curfield.lev[i].x - surf_locx)^2 + (curfield.lev[i].z - surf_locz)^2)
                    z_k = sqrt(curfield.lev[i].x^2 + curfield.lev[i].z^2)

                    fact = abs(gamma_j*gamma_k)*abs(z_j - z_k)/(abs(gamma_j + gamma_k)*(D0 + d_j)^1.5*(D0 + d_k)^1.5)

                    if fact < V0
                        #Merge the 2 vortices into the one at k
                        curfield.lev[i].x = (abs(gamma_j)*curfield.lev[1].x + abs(gamma_k)*curfield.lev[i].x)/(abs(gamma_j + gamma_k))
                        curfield.lev[i].z = (abs(gamma_j)*curfield.lev[1].z + abs(gamma_k)*curfield.lev[i].z)/(abs(gamma_j + gamma_k))
                        curfield.lev[i].s += curfield.lev[1].s

                        popfirst!(curfield.lev)

                        break
                    end
                end
            end
        end
    end
end

function update_kinem2DOF(surf::TwoDSurf, strpar :: TwoDOFPar, kinem :: KinemPar2DOF, dt, cl, cm)
    #Update previous terms
    kinem.alpha_pr = kinem.alpha
    kinem.h_pr = kinem.h

    kinem.alphadot_pr3 = kinem.alphadot_pr2
    kinem.alphadot_pr2 = kinem.alphadot_pr
    kinem.alphadot_pr = kinem.alphadot

    kinem.hdot_pr3 = kinem.hdot_pr2
    kinem.hdot_pr2 = kinem.hdot_pr
    kinem.hdot_pr = kinem.hdot

    kinem.alphaddot_pr3 = kinem.alphaddot_pr2
    kinem.alphaddot_pr2 = kinem.alphaddot_pr

    kinem.hddot_pr3 = kinem.hddot_pr2
    kinem.hddot_pr2 = kinem.hddot_pr

    #Calculate hddot and alphaddot from forces based on 2DOF structural system
    m11 = 2. /surf.c
    m12 = -strpar.x_alpha*cos(kinem.alpha)
    m21 = -2. *strpar.x_alpha*cos(kinem.alpha)/surf.c
    m22 = strpar.r_alpha*strpar.r_alpha

    R1 = 4*strpar.kappa*surf.uref*surf.uref*cl/(pi*surf.c*surf.c) - 2*strpar.w_h*strpar.w_h*(strpar.cubic_h_1*kinem.h + strpar.cubic_h_3*kinem.h^3)/surf.c - strpar.x_alpha*sin(kinem.alpha)*kinem.alphadot*kinem.alphadot

    R2 = 8*strpar.kappa*surf.uref*surf.uref*cm/(pi*surf.c*surf.c) - strpar.w_alpha*strpar.w_alpha*strpar.r_alpha*strpar.r_alpha*(strpar.cubic_alpha_1*kinem.alpha + strpar.cubic_alpha_3*kinem.alpha^3)

    kinem.hddot_pr = (1/(m11*m22-m21*m12))*(m22*R1-m12*R2)
    kinem.alphaddot_pr = (1/(m11*m22-m21*m12))*(-m21*R1+m11*R2)

    #Kinematics are updated according to the 2DOF response
    kinem.alphadot = kinem.alphadot_pr + (dt/12.)*(23*kinem.alphaddot_pr - 16*kinem.alphaddot_pr2 + 5*kinem.alphaddot_pr3)
    kinem.hdot = kinem.hdot_pr + (dt/12)*(23*kinem.hddot_pr-16*kinem.hddot_pr2 + 5*kinem.hddot_pr3)

    kinem.alpha = kinem.alpha_pr + (dt/12)*(23*kinem.alphadot_pr-16*kinem.alphadot_pr2 + 5*kinem.alphadot_pr3)
    kinem.h = kinem.h_pr + (dt/12)*(23*kinem.hdot_pr - 16*kinem.hdot_pr2 + 5*kinem.hdot_pr3)

    #Copy over these values to the TwoDSurf
    surf.kinem.alphadot = kinem.alphadot
    surf.kinem.hdot = kinem.hdot
    surf.kinem.alpha = kinem.alpha
    surf.kinem.h = kinem.h

    return surf, kinem
end

function IFR(surf::TwoDSurf,x,z,t::Float64)
    # Given body frame and kinematics, find global positions
    alpha = surf.kindef.alpha(t)
    X0 = -surf.kindef.u(t)*t
    Z0 = surf.kindef.h(t)
    pvt = surf.pvt

    X = (x .- pvt).*cos(alpha) + z.*sin(alpha) .+ X0 .+ pvt
    Z = -(x .- pvt).*sin(alpha) + z.*cos(alpha) .+ Z0

    return X,Z
end

function BFR(surf::TwoDSurf,X,Z,t::Float64)
    # Given inertial frame and kinematics, find body positions
    alpha = surf.kindef.alpha(t)
    X0 = -surf.kindef.u(t)*t
    Z0 = surf.kindef.h(t)
    pvt = surf.pvt

    x = (X .- X0 .- pvt).*cos(alpha) - (Z .- Z0).*sin(alpha) .+ pvt
    z = (X .- X0 .- pvt).*sin(alpha) + (Z .- Z0).*cos(alpha)

    return x,z
end

function refresh_vortex(surf::TwoDSurf,vor_loc)
    # Updates vortex locations to vor_loc
    for i=1:length(surf.bv)
        surf.bv[i].x = vor_loc[i,1]
        surf.bv[i].z = vor_loc[i,2]
    end
    return surf
end

function place_tev2(surf::TwoDSurf,field::TwoDFlowField,dt,t)
    TE_x, TE_z = IFR(surf,surf.x[end],surf.cam[end],t)
    #=
    if size(field.tev,1) == 0
        xloc = surf.bv[end].x + 0.5*surf.kinem.u*dt
        zloc = surf.bv[end].z
    else
        xloc = surf.bv[end].x + (1. /3.)*(field.tev[end].x - surf.bv[end].x)
        zloc = surf.bv[end].z + (1. /3.)*(field.tev[end].z - surf.bv[end].z)
    end
    =#
    if size(field.tev,1) == 0
        xloc = TE_x + 0.5*surf.kinem.u*dt
        zloc = TE_z + 0.5*surf.kinem.hdot*dt
    else
        xloc = TE_x + (1. /3.)*(field.tev[end].x - TE_x)
        zloc = TE_z + (1. /3.)*(field.tev[end].z - TE_z)
    end
    push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))
    return field
end

# Places a leading edge vortex
function place_lev2(surf::TwoDSurf,field::TwoDFlowField,dt,t)
    LE_x, LE_z = IFR(surf,surf.x[1],surf.cam[1],t)

    le_vel_x = surf.kinem.u - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1]
    le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c- surf.kinem.hdot + surf.wind[1]

    if (surf.levflag[1] == 0)
        xloc = LE_x + 0.5*le_vel_x*dt
        zloc = LE_z + 0.5*le_vel_z*dt
    else
        xloc = LE_x + (1. /3.)*(field.lev[end].x - LE_x)
        zloc = LE_z + (1. /3.)*(field.lev[end].z - LE_z)
    end

    push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))

    return field
end

function vor_BFR(surf::TwoDSurf,t,IFR_field::TwoDFlowField,curfield::TwoDFlowField)
    # finds and updates the BFR positions of the vortexes in IFR_field to
    #   curfield
    ntev = size(IFR_field.tev,1)
    nlev = size(IFR_field.lev,1)
    vorX = [map(q -> q.x,IFR_field.tev); map(q -> q.x,IFR_field.lev)]
    vorZ = [map(q -> q.z,IFR_field.tev); map(q -> q.z,IFR_field.lev)]

    # Body frame conversion
    vorx, vorz = BFR(surf,vorX,vorZ,t)

    # Update bfr TEV
    for i = 1:ntev
        curfield.tev[i].x = vorx[i]
        curfield.tev[i].z = vorz[i]
    end
    # Update bfr LEV
    for i = 1:nlev
        curfield.lev[i].x = vorx[i+ntev]
        curfield.lev[i].z = vorz[i+ntev]
    end
    return curfield
end

function influence_coeff(surf::TwoDSurf,curfield::TwoDFlowField,coll_loc,n,dt,x_w,z_w)
    # With the surface and flowfield, determine the influence matrix "a"
    a = zeros(surf.ndiv,surf.ndiv)
    a[end,:] .= 1. # for wake portion
    vc = surf.bv[1].vc

    for i = 1:surf.ndiv-1, j = 1:surf.ndiv
        # i is test location, j is vortex source
        t_x = coll_loc[i,1]
        t_z = coll_loc[i,2]

        if j < surf.ndiv # Bound vortices (a_ij)
            src = surf.bv[j]
            xdist = src.x .- t_x
            zdist = src.z .- t_z

            distsq = xdist.*xdist + zdist.*zdist # dist_type 1
            u = -zdist./(2*pi*distsq)
            w = xdist./(2*pi*distsq)
        else # Wake vorticy (a_iw)
            xdist = x_w .- t_x
            zdist = z_w .- t_z

            distsq = xdist.*xdist + zdist.*zdist # dist type 2
            u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
            w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
        end

        a[i,j] = u*n[i,1] + w*n[i,2]
    end
    return a
end

function mod_influence_coeff(surf::TwoDSurf,curfield::TwoDFlowField,coll_loc,n,dt,x_w,z_w,x_lev,z_lev)
    # With the surface and flowfield, determine the influence matrix "a" including the
    # modifications needed to calculate the LEV strength
    a = zeros(surf.ndiv,surf.ndiv)
    a[end,:] .= 1. # for wake portion
    vc = surf.bv[1].vc

    for i = 1:surf.ndiv-1, j = 1:surf.ndiv
        # i is test location, j is vortex source
        t_x = coll_loc[i,1]
        t_z = coll_loc[i,2]

        if j < surf.ndiv # Bound vortices (a_ij)
            src = surf.bv[j]
            xdist = src.x .- t_x
            zdist = src.z .- t_z

            distsq = xdist.*xdist + zdist.*zdist # dist_type 1
            u = -zdist./(2*pi*distsq)
            w = xdist./(2*pi*distsq)
        else # Wake vorticy (a_iw)
            xdist = x_w .- t_x
            zdist = z_w .- t_z

            distsq = xdist.*xdist + zdist.*zdist # dist type 2
            u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
            w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
        end

        a[i,j] = u*n[i,1] + w*n[i,2]
    end
    a1 = a[:,1]
    a[:,1:end-1] = a[:,2:end] # shift over a matrix

    for i = 1:surf.ndiv-1 # calc LEV terms
        t_x = coll_loc[i,1]
        t_z = coll_loc[i,2]

        xdist = x_lev .- t_x
        zdist = z_lev .- t_z

        distsq = xdist.*xdist + zdist.*zdist # dist type 2
        u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
        w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))

        a[i,end] = u*n[i,1] + w*n[i,2]
    end

    return a,a1
end


text:
using Symata
using Test

## FIXME: needed ?
import Test: @test

@Symata.sym TimeOff()   # don't print hundreds of diagnostic lines

# For debugging
# @Symata.ex VersionInfo()

# Run tests for Julia-side interface
include("juliainterface_test.jl")
include("mxpr_test.jl")

# Run Symata language tests
function runtests()
    eval(Meta.parse("@sym Tests()"))
end
@test (runtests() ; true)


text:
# # RNN Character level

# In this example, we create a character-level recurrent neural network.
# A recurrent neural network (RNN) outputs a prediction and a hidden state at each step
# of the computation. The hidden state captures historical information of a sequence
# (i.e. the neural network has memory) and the output is the final prediction of the model.
# We use this type of neural network to model sequences such as text or time series.

# ![char-rnn](../char-rnn/docs/rnn-train.png)

# Source: https://d2l.ai/chapter_recurrent-neural-networks/rnn.html#rnn-based-character-level-language-models

# This example demonstrates the use of Flux’s implementation of the
# [Long Short Term Memory recurrent layer](https://www.researchgate.net/publication/13853244_Long_Short-term_Memory)(LSTM)
# which is an RNN that generally exhibits a longer memory span over sequences as well as
# [Flux utility functions](https://fluxml.ai/Flux.jl/stable/utilities/).

# If you need more information about how RNNs work and related technical concepts,
# check out the following resources:

# * [The Unreasonable Effectiveness of Recurrent Neural Networks](http://karpathy.github.io/2015/05/21/rnn-effectiveness/)
# * [Understanding LSTM Networks](https://colah.github.io/posts/2015-08-Understanding-LSTMs/)
# * [Illustrated Guide to Recurrent Neural Networks: Understanding the Intuition](https://www.youtube.com/watch?v=LHXXI4-IEns)

# To run this example, we need the following packages:

using Flux
using Flux: onehot, chunk, batchseq, throttle, logitcrossentropy
using StatsBase: wsample
using Base.Iterators: partition
using Parameters: @with_kw

# We set default values for the hyperparameters:

@with_kw mutable struct Args
    lr::Float64 = 1e-2	# Learning rate
    seqlen::Int = 50	# Length of batch sequences
    nbatch::Int = 50	# Number of batches text is divided into
    throttle::Int = 30	# Throttle timeout
    epochs::Int = 2     # Number of Epochs
end

# ## Data

# We create the function `getdata` to download the training data and create arrays of batches
# for training the model:

function getdata(args)
    ## Download the data if not downloaded as 'input.txt'
    isfile("input.txt") ||
        download("https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt","input.txt")

    text = collect(String(read("input.txt")))

    ## an array of all unique characters
    alphabet = [unique(text)..., '_']

    text = map(ch -> onehot(ch, alphabet), text)
    stop = onehot('_', alphabet)

    N = length(alphabet)

    ## Partitioning the data as sequence of batches, which are then collected as array of batches
    Xs = collect(partition(batchseq(chunk(text, args.nbatch), stop), args.seqlen))
    Ys = collect(partition(batchseq(chunk(text[2:end], args.nbatch), stop), args.seqlen))

    return Xs, Ys, N, alphabet
end

# The function `getdata` performs the following tasks:

# * Downloads a dataset of [all of Shakespeare's works (concatenated)](https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt)
# if not previously downloaded. This function loads the data as a vector of characters with the function `collect`.
# * Gets the alphabet. It consists of the unique characters of the data and the stop character ‘_’.
# * One-hot encodes the alphabet and the stop character.
# * Gets the size of the alphabet N.
# * Partitions the data as an array of batches. Note that the `Xs` array contains the sequence of characters in the text whereas the `Ys` array contains the next character of the sequence.

# ## Model

# We create the RNN with two Flux’s LSTM layers and an output layer of the size of the alphabet:

function build_model(N)
    return Chain(
            LSTM(N, 128),
            LSTM(128, 128),
            Dense(128, N))
end

# The size of the input and output layers is the same as the size of the alphabet.
# Also, we set the size of the hidden layers to 128.

# ## Train the model

# Now, we define the function `train` that creates the model and the loss function as well as the training loop:

function train(; kws...)
    ## Initialize the parameters
    args = Args(; kws...)

    ## Get Data
    Xs, Ys, N, alphabet = getdata(args)

    ## Constructing Model
    m = build_model(N)

    function loss(xs, ys)
        Flux.reset!(m)
        return sum(logitcrossentropy.([m(x) for x in xs], ys))
    end

    ## Training
    opt = ADAM(args.lr)
    tx, ty = (Xs[5], Ys[5])
    evalcb = () -> @show loss(tx, ty)

    @info "Start Training, total $(args.epochs) epochs"
    for epoch = 1:args.epochs
        @info "Epoch $(epoch) / $(args.epochs)"
        Flux.train!(loss, params(m), zip(Xs, Ys), opt, cb = throttle(evalcb, args.throttle))
    end
    return m, alphabet
end

# The function `train` performs the following tasks:

# * Calls the function `getdata` to obtain the train and test data as well as the alphabet and its size.
# * Calls the function `build_model` to create the RNN.
# * Defines the loss function. For this type of neural network, we use the [logitcrossentropy](https://fluxml.ai/Flux.jl/stable/models/losses/#Flux.Losses.logitcrossentropy)
# loss function. Notice that it is important that we call the function [reset!](https://fluxml.ai/Flux.jl/stable/models/layers/#Flux.reset!)
# before computing the loss so that it resets the hidden state of a recurrent layer back to its original value
# * Sets the [ADAM optimiser](https://fluxml.ai/Flux.jl/stable/training/optimisers/#Flux.Optimise.RADAM) with the learning rate *lr* we defined above.
# * Creates a [callback](https://fluxml.ai/Flux.jl/stable/training/training/#Callbacks) *evalcb* so that you can observe the training process (print the loss value).
# * Runs the training loop using [Flux’s train!](https://fluxml.ai/Flux.jl/stable/training/training/#Flux.Optimise.train!).
# It uses the function [throttle](https://fluxml.ai/Flux.jl/stable/utilities/#Flux.throttle) so that the callback *evalcb*
# can only be triggered at most once during timeout seconds (as defined above).

# ## Test the model

# We define the function `sample_data` to test the model.
# It generates samples of text with the alphabet that the function `getdata` computed.
# Notice that it obtains the model’s prediction by calling the
# [softmax function](https://fluxml.ai/Flux.jl/stable/models/nnlib/#Softmax)
# to get the probability distribution of the output and then it chooses randomly the prediction.

function sample_data(m, alphabet, len; seed="")
    m = cpu(m)
    Flux.reset!(m)
    buf = IOBuffer()
    if seed == ""
        seed = string(rand(alphabet))
    end
    write(buf, seed)
    c = wsample(alphabet, softmax([m(onehot(c, alphabet)) for c in collect(seed)][end]))
    for i = 1:len
        write(buf, c)
        c = wsample(alphabet, softmax(m(onehot(c, alphabet))))
    end
    return String(take!(buf))
end

# Finally, to run this example we call the functions `train` and `sample_data`:

cd(@__DIR__)
m, alphabet = train()
sample_data(m, alphabet, 1000) |> println


text:
using Distributions
using Test
using DataFrames
using CSVFiles
using DataStructures

import GlobalSensitivityAnalysis: ishigami, split_output, sample

################################################################################
## JULIA
################################################################################

# define the (uncertain) parameters of the problem and their distributions
data = SobolData(
    params = OrderedDict(:x1 => Uniform(-3.14159265359, 3.14159265359),
        :x2 => Uniform(-3.14159265359, 3.14159265359),
        :x3 => Uniform(-3.14159265359, 3.14159265359)),
    N = 1000
)

N = data.N
D = length(data.params)

# sampling
julia_samples = sample(data)
julia_ishigami = ishigami(julia_samples)

# analysis
julia_A, julia_B, julia_AB, julia_BA = split_output(julia_ishigami, N, D, data.calc_second_order)
julia_results = analyze(data, julia_ishigami; num_resamples = 10_000)

################################################################################
## Python
################################################################################

# sampling
py_samples = load("data/sobol/py_ishigami/py_samples.csv", header_exists=false, colnames = ["x1", "x2", "x3"]) |> DataFrame
py_ishigami = load("data/sobol/py_ishigami/py_ishigami.csv", header_exists=false) |> DataFrame

# analysis
py_A = load("data/sobol/py_ishigami/py_A.csv", header_exists=false) |> DataFrame
py_B = load("data/sobol/py_ishigami/py_B.csv", header_exists=false) |> DataFrame
py_AB = load("data/sobol/py_ishigami/py_AB.csv", header_exists=false) |> DataFrame
py_BA = load("data/sobol/py_ishigami/py_BA.csv", header_exists=false) |> DataFrame

py_firstorder = load("data/sobol/py_ishigami/py_firstorder.csv", header_exists=false) |> DataFrame
py_secondorder = load("data/sobol/py_ishigami/py_secondorder.csv", header_exists=false) |> DataFrame
py_totalorder = load("data/sobol/py_ishigami/py_totalorder.csv", header_exists=false) |> DataFrame

py_firstorder_conf = load("data/sobol/py_ishigami/py_firstorder_conf.csv", header_exists=false) |> DataFrame
py_secondorder_conf = load("data/sobol/py_ishigami/py_secondorder_conf.csv", header_exists=false) |> DataFrame
py_totalorder_conf = load("data/sobol/py_ishigami/py_totalorder_conf.csv", header_exists=false) |> DataFrame

################################################################################
## Testing
################################################################################

@testset "Ishigami Function - Sampling" begin
    @test julia_samples ≈ Matrix(py_samples) atol = ATOL_sample
end

@testset "Ishigami Function - Analysis" begin
    @test julia_ishigami ≈ Matrix(py_ishigami) atol = ATOL_sobol
    @test julia_A ≈ Matrix(py_A) atol = ATOL_sobol
    @test julia_B ≈ Matrix(py_B) atol = ATOL_sobol
    @test julia_AB ≈ Matrix(py_AB) atol = ATOL_sobol
    @test julia_BA ≈ Matrix(py_BA) atol = ATOL_sobol

    @test julia_results[:firstorder] ≈ Matrix(py_firstorder) atol = ATOL_sobol
    @test julia_results[:totalorder] ≈ Matrix(py_totalorder) atol = ATOL_sobol

    for i = 1:D
        for j = i+1:D
            @test julia_results[:secondorder][i,j] ≈ Matrix(py_secondorder)[i,j] atol = ATOL_sobol
        end
    end

    # test confidence intervals
    @test julia_results[:firstorder_conf] ≈ Matrix(py_firstorder_conf) atol = ATOL_CI
    @test julia_results[:totalorder_conf] ≈ Matrix(py_totalorder_conf) atol = ATOL_CI

    for i = 1:D
        for j = i+1:D
            @test julia_results[:secondorder_conf][i,j] ≈ Matrix(py_secondorder_conf)[i,j] atol = ATOL_CI
        end
    end
end

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
using Distributions

nworkers = 3
nwait = 2
comm_distributions = [Uniform(0, 1), Uniform(0, 1), Uniform(0, 1)]
comp_distributions = [Uniform(0, 1), Uniform(10, 11), Uniform(20, 21)]
sim = EventDrivenSimulator(;nwait, nworkers, comm_distributions, comp_distributions)
step!(sim)
@test sim.epoch == 1
@test 10 <= sim.time <= 12
@test sim.isidle == [true, true, false]
@test sim.nfresh == [1, 1, 0]

step!(sim, 3)
@test sim.epoch == 4
@test sim.nfresh == [4, 4, 0]

sim = EventDrivenSimulator(sim)
step!(sim)
@test sim.epoch == 1
@test 10 <= sim.time <= 12
@test sim.isidle == [true, true, false]
@test sim.nfresh == [1, 1, 0]

empty!(sim)
step!(sim)
@test sim.epoch == 1
@test 10 <= sim.time <= 12
@test sim.isidle == [true, true, false]
@test sim.nfresh == [1, 1, 0]