text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Lexer of SceneLang

"""
    _update_pos!(stream::InputStream, ch::Union{Char, Nothing})

Update `stream.location` after having read `ch` from the stream.
"""
function _update_pos!(stream::InputStream, ch::Union{Char, Nothing})
    if isnothing(ch)
        return
    elseif ch == '\n'
        stream.location.line_num += 1
        stream.location.col_num = 1
    elseif ch == '\t'
        stream.location.col_num += stream.tabulations
    else
        stream.location.col_num += 1
    end
end

"""
    read_char!(stream::InputStream)

Read a new character from the stream.
"""
function read_char!(stream::InputStream)
    if !isnothing(stream.saved_char)
        # Recover the "unread" character and return it
        ch = stream.saved_char
        stream.saved_char = nothing
    else
        # Read a new character from the stream
        ch = eof(stream) ? nothing : read(stream.stream, Char)
    end

    stream.saved_location = copy(stream.location)
    _update_pos!(stream, ch)

    return ch
end

"""
    unread_char!(stream::InputStream, ch::Union{Char, Nothing})

Push a character back to the stream. Accepts a `nothing` value to deal with eof.
"""
function unread_char!(stream::InputStream, ch::Union{Char, Nothing})
    @assert isnothing(stream.saved_char)
    stream.saved_char = ch
    stream.location = copy(stream.saved_location)
end

"""
    skip_whitespaces_and_comments(stream::InputStream)

Keep reading characters until a non-whitespace character is found out of a commented line.
"""
function skip_whitespaces_and_comments(stream::InputStream)
    ch = read_char!(stream)
    isnothing(ch) && return
    while isspace(ch) || ch == '#'
        if ch == '#'
            dump = read_char!(stream)
            exit_condition::Function = if dump == '='
                dump::Char -> dump != '=' ?
                    false :
                    read_char!(stream) == '#'
            else
                isnewline
            end
            while (dump = read_char!(stream); !isnothing(dump) && !exit_condition(dump)) end
        end
        ch = read_char!(stream)
        isnothing(ch) && return
    end

    # Put the non-whitespace character back
    unread_char!(stream, ch)
end

"""
    _parse_string_token(stream::InputStream, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`LiteralString`](@ref) value.
"""
function _parse_string_token(stream::InputStream, token_location::SourceLocation)
    str = ""
    while true
        ch = read_char!(stream)

        ch == '"' && break

        (isnothing(ch) || isnewline(ch)) && throw(UnfinishedExpression(token_location, "Unterminated string", length(str) + 1))

        str *= ch
    end

    return Token(token_location, LiteralString(str), length(str))
end

"""
    _parse_number_token(stream::InputStream, first_char::Char, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`LiteralNumber`](@ref) value.
"""
function _parse_number_token(stream::InputStream, first_char::Char, token_location::SourceLocation)
    str = first_char |> string
    while true
        ch = read_char!(stream)

        if isnothing(ch) || isspace(ch) || issymbol(ch)
            unread_char!(stream, ch)
            break
        end

        str *= ch
    end

    value = try
        parse(Float32, str)
    catch e
        isa(e, ArgumentError) &&
            rethrow(InvalidNumber(token_location, "'$str' is an invalid number format\n"*
                    "Valid formats only include digits, '+' or '-', dots, and letters 'e' or 'E'\n"*
                    "Examples: +9e3, -300.5, 5e-7", length(str)))
        rethrow(e)
    end

    return Token(token_location, LiteralNumber(value), length(str))
end

"""
    _parse_keyword_token(stream::InputStream, first_char::Char, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`Keyword`](@ref) value.
"""
function _parse_keyword_token(stream::InputStream, token_location::SourceLocation)
    str = ""
    while true
        ch = read_char!(stream)
        # Note that here we do not call "isalpha" but "isalnum": digits are ok after the first character
        !isnothing(ch) && (isdigit(ch) || isletter(ch) || ch == '_') || (unread_char!(stream, ch); break)

        str *= ch
    end

    sym = Symbol(str)
    Token(token_location, Keyword(sym), length(str) + 1)
end

"""
    _parse_command_or_type_token(stream::InputStream, first_char::Char, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`Command`](@ref) or [`LiteralType`](@ref) value.
"""
function _parse_command_or_type_token(stream::InputStream, first_char::Char, token_location::SourceLocation)
    start_loc = copy(stream.saved_location)
    str = first_char |> string
    while true
        ch = read_char!(stream)
        # Note that here digits are ok after the first character
        !isnothing(ch) && (isdigit(ch) || isletter(ch) || ch == '_') || (unread_char!(stream, ch); break)

        str *= ch
    end

    class, sym = all(ch -> isuppercase(ch) || ch == '_', str) ?
        (Command,     Symbol(str)) :
        (LiteralType, Symbol(str * "Type"))

    keywords = instances(class) .|> Symbol
    if (index = findfirst(s -> s == sym, keywords)) |> isnothing
        valid_instances = replace.(
            last.(
                split.(
                    repr.(
                        instances(class)
                    ),
                ".")
            ),
            Ref(r"Type$" => "")
        )
        throw((class <: Command ? InvalidCommand : InvalidType)(start_loc,
            "'$str' is not a valid '$class'.\nValid $(class)s are:\n\t" * join(valid_instances, "\n\t"), length(str)))
    end
    Token(token_location, class(index - 1), length(str))
end

"""
    _parse_identifier_token(stream::InputStream, first_char::Char, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`Identifier`](@ref) value.
"""
function _parse_identifier_token(stream::InputStream, first_char::Char, token_location::SourceLocation)
    str = first_char |> string
    while true
        ch = read_char!(stream)
        # Note that here we do not call "isalpha" but "isalnum": digits are ok after the first character
        !isnothing(ch) && (isdigit(ch) || isletter(ch) || ch == '_') || (unread_char!(stream, ch); break)

        str *= ch
    end

    sym = Symbol(str)
    Token(token_location, Identifier(sym), length(str))
end

"""
    _parse_math_expression_token(stream::InputStream, token_location::SourceLocation)

Parse the stream into a [`Token`](@ref) with [`MathExpression`](@ref) value.
"""
function _parse_math_expression_token(stream::InputStream, token_location::SourceLocation)
    str = ""
    endcol = 0
    while true
        ch = read_char!(stream)

        if ch == '$'
            endcol = stream.saved_location.col_num
            break
        end

        (isnothing(ch) || isnewline(ch)) &&
            throw(UnfinishedExpression(token_location, "Unterminated mathematical expression", stream.saved_location.col_num - token_location.col_num + 1))

        str *= ch
    end

    len = endcol - token_location.col_num
    expr = try
        Meta.parse(str)
    catch e
        isa(e, Base.Meta.ParseError) || rethrow(e)
        throw(InvalidExpression(token_location,"Mathemathical expression cannot be parsed:\n$(e.msg)", len + 1))
    end
    isa(expr, Symbol) && (expr = :(+$expr))
    Raytracer.isvalid(expr, len, token_location)
    return Token(token_location, MathExpression(expr), len + 1)
end

"""
    read_token(stream::InputStream)

Read the next token in the stream.
"""
function read_token(stream::InputStream)
    if !isnothing(stream.saved_token)
        result = stream.saved_token
        stream.saved_token = nothing
        return result
    end

    skip_whitespaces_and_comments(stream)

    # At this point we're sure that ch does *not* contain a whitespace character
    ch = read_char!(stream)
    if isnothing(ch)
        # No more characters in the file, so return a StopToken
        return Token(stream.location, StopToken(), 1)
    end

    # Check if we got some non ASCII character
    isascii(ch) || throw(BadCharacter(stream.saved_location, "Invalid character $ch: only ASCII charachters are supported"))

    # At this point we must check what kind of token begins with the "ch" character
    # (which has been put back in the stream with stream.unread_char). First,
    # we save the position in the stream
    token_location = copy(stream.saved_location)

    if issymbol(ch)
        # One-character symbol, like '(' or ','
        return Token(token_location, LiteralSymbol(Symbol(ch)), 1)
    elseif ch == '"'
        # A literal string (used for file names)
        return _parse_string_token(stream, token_location)
    elseif ch == '$'
        # A math expression
        return _parse_math_expression_token(stream, token_location)
    elseif ch == '.'
        return _parse_keyword_token(stream, token_location)
    elseif isdigit(ch) || ch ∈ ('+', '-')
        # A floating-point number
        return _parse_number_token(stream, ch, token_location)
    elseif isuppercase(ch)
        # Since it begins with an uppercase letter, it must be either a command or a type
        return _parse_command_or_type_token(stream, ch, token_location)
    elseif islowercase(ch) || ch == '_'
        # Since it begins with a lowercase letter, it must be an identifier
        return _parse_identifier_token(stream, ch, token_location)
    else
        # We got some weird ASCII character, like '@` or `&`
        throw(BadCharacter(stream.location, "Invalid character $ch"))
    end
end

"""
    unread_token(stream::InputStream, ch::Token)

Push a token back to the stream.
"""
function unread_token(stream::InputStream, token::Token)
    @assert isnothing(stream.saved_token)
    stream.saved_token = token
end


text:
using PackageCompiler: PackageCompiler, create_sysimage, create_app, create_library
using Test
using Libdl
using Pkg

ENV["JULIA_DEBUG"] = "PackageCompiler"

# Make a new depot
new_depot = mktempdir()
mkpath(joinpath(new_depot, "registries"))
cp(joinpath(DEPOT_PATH[1], "registries", "General"), joinpath(new_depot, "registries", "General"))
ENV["JULIA_DEPOT_PATH"] = new_depot
Base.init_depot_path()

is_slow_ci = haskey(ENV, "CI") && Sys.ARCH == :aarch64

if haskey(ENV, "CI")
    @show Sys.ARCH
end

@testset "PackageCompiler.jl" begin
    new_project = mktempdir()
    old_project = Base.ACTIVE_PROJECT[]
    Base.ACTIVE_PROJECT[] = new_project
    try
        Pkg.add("Example")
    finally
        Base.ACTIVE_PROJECT[] = old_project
    end
    tmp = mktempdir()
    sysimage_path = joinpath(tmp, "sys." * Libdl.dlext)
    script = tempname()
    write(script, """
    script_func() = println(\"I am a script\")
    opt_during_sysimage = Base.JLOptions().opt_level
    print_opt() = println("opt: -O\$opt_during_sysimage")
    """)
    create_sysimage(; sysimage_path=sysimage_path,
                              project=new_project,
                              precompile_execution_file=joinpath(@__DIR__, "precompile_execution.jl"),
                              precompile_statements_file=joinpath.(@__DIR__, ["precompile_statements.jl",
                                                                              "precompile_statements2.jl"]),
                              script=script,
                              sysimage_build_args = `-O1`
                              )

    # Check we can load sysimage and that Example is available in Main
    str = read(`$(Base.julia_cmd()) -J $(sysimage_path) -e 'println(Example.hello("foo")); script_func(); print_opt()'`, String)
    @test occursin("Hello, foo", str)
    @test occursin("I am a script", str)
    @test occursin("opt: -O1", str)

    # Test creating an app
    app_source_dir = joinpath(@__DIR__, "..", "examples/MyApp/")
    app_compiled_dir = joinpath(tmp, "MyAppCompiled")
    for incremental in (is_slow_ci ? (false,) : (true, false))
        if incremental == false
            filter_stdlibs = (is_slow_ci ? (true, ) : (true, false))
        else
            filter_stdlibs = (false,)
        end
        for filter in filter_stdlibs
            tmp_app_source_dir = joinpath(tmp, "MyApp")
            cp(app_source_dir, tmp_app_source_dir)
            create_app(tmp_app_source_dir, app_compiled_dir; incremental=incremental, force=true, filter_stdlibs=filter,
                       precompile_execution_file=joinpath(app_source_dir, "precompile_app.jl"),
                       executables=["MyApp" => "julia_main",
                                    "SecondApp" => "second_main",
                                    "ReturnType" => "wrong_return_type",
                                    "Error" => "erroring",
                                    "Undefined" => "undefined",
                                    ])
            rm(tmp_app_source_dir; recursive=true)
            # Get rid of some local state
            rm(joinpath(new_depot, "packages"); recursive=true)
            rm(joinpath(new_depot, "compiled"); recursive=true)
            rm(joinpath(new_depot, "artifacts"); recursive=true)
            app_path(app_name) = abspath(app_compiled_dir, "bin", app_name * (Sys.iswindows() ? ".exe" : ""))
            app_output = read(`$(app_path("MyApp")) I get --args --julia-args --threads=3 --check-bounds=yes -O1`, String)

            # Check stdlib filtering
            if filter == true
                @test !(occursin("LinearAlgebra", app_output))
            else
                @test occursin("LinearAlgebra", app_output)
            end
            # Check dependency run
            @test occursin("Example.domath", app_output)
            # Check PROGRAM_FILE
            @test occursin("Base.PROGRAM_FILE = $(repr(app_path("MyApp")))", app_output)
            # Check jll package runs
            @test occursin("Hello, World!", app_output)
            # Check artifact runs
            @test occursin("The result of 2*5^2 - 10 == 40.000000", app_output)
            # Check artifact gets run from the correct place
            @test occursin("HelloWorld artifact at $(realpath(app_compiled_dir))", app_output)
            # Check ARGS
            @test occursin("""ARGS = ["I", "get", "--args"]""", app_output)
            # Check julia-args
            @test occursin("(Base.JLOptions()).opt_level = 1", app_output)
            @test occursin("(Base.JLOptions()).nthreads = 3", app_output)
            @test occursin("(Base.JLOptions()).check_bounds = 1", app_output)
            # Check transitive inclusion of dependencies
            @test occursin("is_crayons_loaded() = true", app_output)
            # Check app is precompiled in a normal process
            @test occursin("outputo: ok", app_output)
            @test occursin("myrand: ok", app_output)
            # Check distributed
            @test occursin("n = 20000000", app_output)
            @test occursin("From worker 2:\t8", app_output)
            @test occursin("From worker 3:\t8", app_output)
            @test occursin("From worker 4:\t8", app_output)
            @test occursin("From worker 5:\t8", app_output)

            # Test second app
            app_output = read(`$(app_path("SecondApp"))`, String)
            @test occursin("Hello from second main", app_output)

            io = IOBuffer()
            p = run(pipeline(ignorestatus(`$(app_path("ReturnType"))`), stderr=io;))
            @test occursin("ERROR: expected a Cint return value from function MyApp.wrong_return_type", String(take!(io)))
            @test p.exitcode == 1

            io = IOBuffer()
            p = run(pipeline(ignorestatus(`$(app_path("Error"))`), stderr=io;))
            @test occursin("MethodError: no method matching +(", String(take!(io)))
            @test p.exitcode == 1

            io = IOBuffer()
            p = run(pipeline(ignorestatus(`$(app_path("Undefined"))`), stderr=io;))
            @test occursin("UndefVarError: undefined not defined", String(take!(io)))
            @test p.exitcode == 1
        end
    end

    if !is_slow_ci
        # Test library creation
        lib_source_dir = joinpath(@__DIR__, "..", "examples/MyLib")
        lib_target_dir = joinpath(tmp, "MyLibCompiled")

        incremental = false
        filter = true
        lib_name = "inc"

        tmp_lib_src_dir = joinpath(tmp, "MyLib")
        cp(lib_source_dir, tmp_lib_src_dir)
        create_library(tmp_lib_src_dir, lib_target_dir; incremental=incremental, force=true, filter_stdlibs=filter,
                    precompile_execution_file=joinpath(lib_source_dir, "build", "generate_precompile.jl"),
                    precompile_statements_file=joinpath(lib_source_dir, "build", "additional_precompile.jl"),
                    lib_name=lib_name, version=v"1.0.0")
        rm(tmp_lib_src_dir; recursive=true)
    end

    # Test creating an empty sysimage
    if !is_slow_ci
        tmp = mktempdir()
        sysimage_path = joinpath(tmp, "empty." * Libdl.dlext)
        foreach(x -> touch(joinpath(tmp, x)), ["Project.toml", "Manifest.toml"])
        create_sysimage(String[]; sysimage_path=sysimage_path, incremental=false, filter_stdlibs=true, project=tmp)
        hello = read(`$(Base.julia_cmd()) -J $(sysimage_path) -e 'print("hello, world")'`, String)
        @test hello == "hello, world"
    end
end

ext: jl


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

using IntervalArithmetic
using Base.Test

setprecision(Interval, 128)
setprecision(Interval, Float64)

@testset "sin" begin
    @test sin(@interval(0.5)) == Interval(0.47942553860420295, 0.47942553860420301)
    @test sin(@interval(0.5, 1.67)) == Interval(4.7942553860420295e-01, 1.0)
    @test sin(@interval(1.67, 3.2)) == Interval(-5.8374143427580093e-02, 9.9508334981018021e-01)
    @test sin(@interval(2.1, 5.6)) == Interval(-1.0, 0.863209366648874)
    @test sin(@interval(0.5, 8.5)) == Interval(-1.0, 1.0)
    @test sin(@floatinterval(-4.5, 0.1)) == Interval(-1.0, 0.9775301176650971)
    @test sin(@floatinterval(1.3, 6.3)) == Interval(-1.0, 1.0)

    @test sin(@biginterval(0.5)) ⊆ sin(@interval(0.5))
    @test sin(@biginterval(0.5, 1.67)) ⊆ sin(@interval(0.5, 1.67))
    @test sin(@biginterval(1.67, 3.2)) ⊆ sin(@interval(1.67, 3.2))
    @test sin(@biginterval(2.1, 5.6)) ⊆ sin(@interval(2.1, 5.6))
    @test sin(@biginterval(0.5, 8.5)) ⊆ sin(@interval(0.5, 8.5))
    @test sin(@biginterval(-4.5, 0.1)) ⊆ sin(@interval(-4.5, 0.1))
    @test sin(@biginterval(1.3, 6.3)) ⊆ sin(@interval(1.3, 6.3))
end

@testset "cos" begin
    @test cos(@interval(0.5)) == Interval(0.87758256189037265, 0.87758256189037276)
    @test cos(@interval(0.5, 1.67)) == Interval(-0.09904103659872825, 0.8775825618903728)
    @test cos(@interval(2.1, 5.6)) == Interval(-1.0, 0.77556587851025016)
    @test cos(@interval(0.5, 8.5)) == Interval(-1.0, 1.0)
    @test cos(@interval(1.67, 3.2)) == Interval(-1.0, -0.09904103659872801)

    @test cos(@biginterval(0.5)) ⊆ cos(@interval(0.5))
    @test cos(@biginterval(0.5, 1.67)) ⊆ cos(@interval(0.5, 1.67))
    @test cos(@biginterval(1.67, 3.2)) ⊆ cos(@interval(1.67, 3.2))
    @test cos(@biginterval(2.1, 5.6)) ⊆ cos(@interval(2.1, 5.6))
    @test cos(@biginterval(0.5, 8.5)) ⊆ cos(@interval(0.5, 8.5))
    @test cos(@biginterval(-4.5, 0.1)) ⊆ cos(@interval(-4.5, 0.1))
    @test cos(@biginterval(1.3, 6.3)) ⊆ cos(@interval(1.3, 6.3))
end

@testset "tan" begin
    @test tan(@interval(0.5)) == Interval(0.54630248984379048, 0.5463024898437906)
    @test tan(@interval(0.5, 1.67)) == entireinterval()
    @test tan(@interval(1.67, 3.2)) == Interval(-10.047182299210307, 0.05847385445957865)
    @test tan(Interval(6.638314112824137, 8.38263151220128)) == entireinterval()  # https://github.com/JuliaIntervals/IntervalArithmetic.jl/pull/20

    @test tan(@biginterval(0.5)) ⊆ tan(@interval(0.5))
    @test tan(@biginterval(0.5, 1.67)) == entireinterval(BigFloat)
    @test tan(@biginterval(0.5, 1.67)) ⊆ tan(@interval(0.5, 1.67))
    @test tan(@biginterval(1.67, 3.2)) ⊆ tan(@interval(1.67, 3.2))
    @test tan(@biginterval(2.1, 5.6)) ⊆ tan(@interval(2.1, 5.6))
    @test tan(@biginterval(0.5, 8.5)) ⊆ tan(@interval(0.5, 8.5))
    @test tan(@biginterval(-4.5, 0.1)) ⊆ tan(@interval(-4.5, 0.1))
    @test tan(@biginterval(1.3, 6.3)) ⊆ tan(@interval(1.3, 6.3))

end

@testset "Inverse trig" begin
    @test asin(@interval(1)) == @interval(pi/2)#pi_interval(Float64)/2
    @test asin(@interval(0.9, 2)) == asin(@interval(0.9, 1))
    @test asin(@interval(3, 4)) == ∅

    @test asin(@biginterval(1)) ⊆ asin(@interval(1))
    @test asin(@biginterval(0.9, 2)) ⊆ asin(@interval(0.9, 2))
    @test asin(@biginterval(3, 4)) ⊆ asin(@interval(3, 4))

    @test acos(@interval(1)) == Interval(0., 0.)
    @test acos(@interval(-2, -0.9)) == acos(@interval(-1, -0.9))
    @test acos(@interval(3, 4)) == ∅

    @test acos(@biginterval(1)) ⊆ acos(@interval(1))
    @test acos(@biginterval(-2, -0.9)) ⊆ acos(@interval(-2, -0.9))
    @test acos(@biginterval(3, 4)) ⊆ acos(@interval(3, 4))

    @test atan(@interval(-1,1)) ==
        Interval(-pi_interval(Float64).hi/4, pi_interval(Float64).hi/4)
    @test atan(@interval(0)) == Interval(0.0, 0.0)
    @test atan(@biginterval(-1, 1)) ⊆ atan(@interval(-1, 1))
end

@testset "atan2" begin
    @test atan2(∅, entireinterval()) == ∅
    @test atan2(entireinterval(), ∅) == ∅
    @test atan2(@interval(0.0, 1.0), @biginterval(0.0)) == @biginterval(pi/2)
    @test atan2(@interval(0.0, 1.0), @interval(0.0)) == @interval(pi/2)
    @test atan2(@interval(-1.0, -0.1), @interval(0.0)) == @interval(-pi/2)
    @test atan2(@interval(-1.0, 1.0), @interval(0.0)) == @interval(-pi/2, pi/2)
    @test atan2(@interval(0.0), @interval(0.1, 1.0)) == @interval(0.0)
    @test atan2(@biginterval(0.0, 0.1), @biginterval(0.1, 1.0)) ⊆
        atan2(@interval(0.0, 0.1), @interval(0.1, 1.0))
    @test atan2(@interval(0.0, 0.1), @interval(0.1, 1.0)) ==
        Interval(0.0, 0.7853981633974484)
    @test atan2(@biginterval(-0.1, 0.0), @biginterval(0.1, 1.0)) ⊆
        atan2(@interval(-0.1, 0.0), @interval(0.1, 1.0))
    @test atan2(@interval(-0.1, 0.0), @interval(0.1, 1.0)) ==
        Interval(-0.7853981633974484, 0.0)
    @test atan2(@biginterval(-0.1, -0.1), @biginterval(0.1, Inf)) ⊆
        atan2(@interval(-0.1, -0.1), @interval(0.1, Inf))
    @test atan2(@interval(-0.1, 0.0), @interval(0.1, Inf)) ==
        Interval(-0.7853981633974484, 0.0)
    @test atan2(@biginterval(0.0, 0.1), @biginterval(-2.0, -0.1)) ⊆
        atan2(@interval(0.0, 0.1), @interval(-2.0, -0.1))
    @test atan2(@interval(0.0, 0.1), @interval(-2.0, -0.1)) ==
        Interval(2.356194490192345, 3.1415926535897936)
    @test atan2(@biginterval(-0.1, 0.0), @biginterval(-2.0, -0.1)) ⊆
        atan2(@interval(-0.1, 0.0), @interval(-2.0, -0.1))
    @test atan2(@interval(-0.1, 0.0), @interval(-2.0, -0.1)) ==
        @interval(-pi, pi)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(-Inf, -0.1)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(-Inf, -0.1))
    @test atan2(@interval(-0.1, 0.1), @interval(-Inf, -0.1)) ==
        @interval(-pi, pi)

    @test atan2(@biginterval(0.0, 0.0), @biginterval(-2.0, 0.0)) ⊆
        atan2(@interval(0.0, 0.0), @interval(-2.0, 0.0))
    @test atan2(@interval(-0.0, 0.0), @interval(-2.0, 0.0)) ==
        Interval(3.141592653589793, 3.1415926535897936)
    @test atan2(@biginterval(0.0, 0.1), @biginterval(-0.1, 0.0)) ⊆
        atan2(@interval(0.0, 0.1), @interval(-0.1, 0.0))
    @test atan2(@interval(-0.0, 0.1), @interval(-0.1, 0.0)) ==
        Interval(1.5707963267948966, 3.1415926535897936)
    @test atan2(@biginterval(-0.1, -0.1), @biginterval(-0.1, 0.0)) ⊆
        atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.0))
    @test atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.0)) ==
        Interval(-2.3561944901923453, -1.5707963267948966)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(-2.0, 0.0)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.0))
    @test atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.0)) ==
        @interval(-pi, pi)

    @test atan2(@biginterval(0.0, 0.0), @biginterval(-2.0, 0.0)) ⊆
        atan2(@interval(0.0, 0.0), @interval(-2.0, 0.0))
    @test atan2(@interval(-0.0, 0.0), @interval(-2.0, 0.0)) ==
        Interval(3.141592653589793, 3.1415926535897936)
    @test atan2(@biginterval(0.0, 0.1), @biginterval(-0.1, 0.0)) ⊆
        atan2(@interval(0.0, 0.1), @interval(-0.1, 0.0))
    @test atan2(@interval(-0.0, 0.1), @interval(-0.1, 0.0)) ==
        Interval(1.5707963267948966, 3.1415926535897936)
    @test atan2(@biginterval(-0.1, -0.1), @biginterval(-0.1, 0.0)) ⊆
        atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.0))
    @test atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.0)) ==
        Interval(-2.3561944901923453, -1.5707963267948966)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(-2.0, 0.0)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.0))
    @test atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.0)) ==
        @interval(-pi, pi)
    @test atan2(@biginterval(0.0, 0.1), @biginterval(-2.0, 0.1)) ⊆
        atan2(@interval(0.0, 0.1), @interval(-2.0, 0.1))
    @test atan2(@interval(-0.0, 0.1), @interval(-2.0, 0.1)) ==
        Interval(0.0, 3.1415926535897936)
    @test atan2(@biginterval(-0.1, -0.1), @biginterval(-0.1, 0.1)) ⊆
        atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.1))
    @test atan2(@interval(-0.1, -0.1), @interval(-0.1, 0.1)) ==
        Interval(-2.3561944901923453, -0.7853981633974482)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(-2.0, 0.1)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.1))
    @test atan2(@interval(-0.1, 0.1), @interval(-2.0, 0.1)) ==
        @interval(-pi, pi)

    @test atan2(@interval(-0.1, 0.1), @interval(0.1, 0.1)) ==
        Interval(-0.7853981633974484, 0.7853981633974484)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(0.1, 0.1)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(0.1, 0.1))
    @test atan2(@interval(0.0), @interval(-0.0, 0.1)) == @interval(0.0)
    @test atan2(@interval(0.0, 0.1), @interval(-0.0, 0.1)) ==
        Interval(0.0, 1.5707963267948968)
    @test atan2(@interval(-0.1, 0.0), @interval(0.0, 0.1)) ==
        Interval(-1.5707963267948968, 0.0)
    @test atan2(@interval(-0.1, 0.1), @interval(-0.0, 0.1)) ==
        Interval(-1.5707963267948968, 1.5707963267948968)
    @test atan2(@biginterval(-0.1, 0.1), @biginterval(-0.0, 0.1)) ⊆
        atan2(@interval(-0.1, 0.1), @interval(0.0, 0.1))
end

@testset "Trig" begin
    for a in ( @interval(17, 19), @interval(0.5, 1.2) )
        @test tan(a) ⊆ sin(a)/cos(a)
    end

    @test sin(Interval(-pi/2, 3pi/2)) == Interval(-1, 1)
    @test cos(Interval(-pi/2, 3pi/2)) == Interval(-1, 1)
end

@testset "Trig with large arguments" begin
    x = Interval(2.)^1000   # this is a thin interval
    @test diam(x) == 0.0

    @test sin(x) == -1..1
    @test cos(x) == -1..1
    @test_skip tan(x) == Interval(-0.16125837995065806, -0.16125837995065803)

    x = Interval(prevfloat(∞), ∞)
    @test sin(x) == -1..1
    @test cos(x) == -1..1
    @test tan(x) == -∞..∞
end


text:
function free(cur_player::AbstractPlayer)
  cur_player.action = ""

  cur_player.clock.hold =
    Nullable{UUIDs.UUID}()

  return
end


text:
using Random, Test, TriangulatedWindingNumbers
using TriangulatedWindingNumbers: Vertex, Simplex, windingnumber, windingangle
using TriangulatedWindingNumbers: centroid, assessconvergence, position, value

@testset "TriangulatedWindingNumbers tests" begin

  Random.seed!(0)

  @testset "Simplex tests" begin

    irrelevant = [0.0, 0.0]
    @testset "Simplex encloses zero" begin
      v1 = Vertex(irrelevant, 1.0 - im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, -1.0 - im)
      encloseszero = Simplex{Float64, ComplexF64}([v1, v2, v3])
      @test isapprox(1, abs(windingangle(encloseszero)) / (2π))
      @test windingnumber(encloseszero) == 1
    end

    @testset "Simplex doesn't enclose zero" begin
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 2.0 + im)
      v3 = Vertex(irrelevant, 1.0 + im * 2)
      doesntenclosezero = Simplex{Float64, ComplexF64}([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 0.0 + 2*im)
      v3 = Vertex(irrelevant, -1.0 + im)
      doesntenclosezero = Simplex{Float64, ComplexF64}([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + 0im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, 1.0 + 0im)
      doesntenclosezero = Simplex{Float64, ComplexF64}([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
    end

    @testset "Simplices with identical vertices is converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      pos = T[1.0, 1.0]
      val = rand(U)
      v1 = Vertex(pos, val)
      v2 = Vertex(pos, val)
      v3 = Vertex(pos, val)
      s = Simplex{T,U}([v1, v2, v3])

      defaults = TriangulatedWindingNumbers.convergenceconfig(dim, T)
      isconverged, returncode = assessconvergence(s, defaults)
      @test isconverged
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplices in a simplex eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T) + eps(T), one(T)], rand(U))
      s = Simplex{T,U}([v1, v2, v3])
      defaults = TriangulatedWindingNumbers.convergenceconfig(dim, T)
      isconverged, returncode = assessconvergence(s, defaults)
      @test isconverged
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplices in a chain eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T), one(T) + eps(T) + eps(T)], rand(U))
      s = Simplex{T,U}([v1, v2, v3])
      defaults = TriangulatedWindingNumbers.convergenceconfig(dim, T)
      isconverged, returncode = assessconvergence(s, defaults)
      @assert all(isapprox(position(v1)[d], position(v2)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert all(isapprox(position(v2)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert !all(isapprox(position(v1)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @test isconverged
      @test returncode == :XTOL_REACHED
    end

  end

  @testset "End-to-end tests roots" begin
    @testset "Single root" begin
      for i in 1:10
        gridsize = [rand(1:10), rand(1:10)]
        xtol_abs=10.0^(-rand(3:15))
        function mock(x::Vector, root)
          return (x[1] + im * x[2]) - root
        end
        root = rand(ComplexF64)
        objective(x) = mock(x, root)
        lower = collect(reim(root)) .- rand(2)
        upper = collect(reim(root)) .+ rand(2)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, xtol_abs=xtol_abs)
        @test length(solutions) == 1
        for (s, reason) ∈ solutions
          @test isapprox(centroid(s)[1], real(root), atol=xtol_abs)
          @test isapprox(centroid(s)[2], imag(root), atol=xtol_abs)
        end
      end
    end

    @testset "Single pole" begin
      for i in 1:10
        gridsize = [rand(1:10), rand(1:10)]
        xtol_abs=10.0^(-rand(3:15))
        function mock(x::Vector, root)
          return 1 / ((x[1] + im * x[2]) - root)
        end
        root = rand(ComplexF64)
        objective(x) = mock(x, root)
        lower = collect(reim(root)) .- rand(2)
        upper = collect(reim(root)) .+ rand(2)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, xtol_abs=xtol_abs)
        @test length(solutions) == 1
        for (s, reason) ∈ solutions
          @test isapprox(centroid(s)[1], real(root), atol=xtol_abs)
          @test isapprox(centroid(s)[2], imag(root), atol=xtol_abs)
        end
      end
      for i in 1:10
        gridsize = [rand(1:10), rand(1:10)]
        function mock(x::Vector, root)
          return 1 / ((x[1] + im * x[2]) - root)
        end
        root = rand(ComplexF64)
        objective(x) = mock(x, root)
        lower = collect(reim(root)) .- rand(2)
        upper = collect(reim(root)) .+ rand(2)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, stopvalpole=1e15)
        @test length(solutions) == 1
        for (s, reason) ∈ solutions
          @test isapprox(centroid(s)[1], real(root))
          @test isapprox(centroid(s)[2], imag(root))
        end
      end

    end

    @testset "Multiple roots" begin
      for i in 1:100
        gridsize = [rand(5:20), rand(5:20)]
        xtol_abs=10.0^(-rand(3:15))
        function multimock(x::Vector, roots)
          return mapreduce(root->(x[1] + im * x[2]) - root, *, roots)
        end
        lower = [0.0, 0.0]
        upper = [1.0, 1.0]
        roots = [rand(ComplexF64) for i ∈ 1:rand(1:10)]
        objective(x) = multimock(x, roots)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, xtol_abs=xtol_abs, stopvalroot=1e-20)
        for (s, reason) ∈ solutions
          passed = false
          for root ∈ roots
            passed |= (isapprox(centroid(s)[1], real(root), atol=xtol_abs) &&
                       isapprox(centroid(s)[2], imag(root), atol=xtol_abs))
          end
          @test passed
        end
      end
    end

    @testset "Multiple poles" begin
      for i in 1:100
        gridsize = [rand(5:20), rand(5:20)]
        xtol_abs=10.0^(-rand(3:15))
        function multimock(x::Vector, roots)
          return mapreduce(root->1 / ((x[1] + im * x[2]) - root), *, roots)
        end
        lower = [0.0, 0.0]
        upper = [1.0, 1.0]
        roots = [rand(ComplexF64) for i ∈ 1:rand(1:10)]
        objective(x) = multimock(x, roots)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, xtol_abs=xtol_abs, stopvalpole=1e20)
        for (s, reason) ∈ solutions
          passed = false
          for root ∈ roots
            passed |= (isapprox(centroid(s)[1], real(root), atol=xtol_abs) &&
                       isapprox(centroid(s)[2], imag(root), atol=xtol_abs))
          end
          @test passed
        end
      end
    end

    @testset "xtol_rel is eps and xtol_abs zero stopvalroot is eps" begin
      for i in 1:10
        gridsize = [rand(1:10), rand(1:10)]
        function mock(x::Vector, root)
          return (x[1] + im * x[2]) - root
        end
        root = rand(ComplexF64)
        objective(x) = mock(x, root)
        lower = collect(reim(root)) .- rand(2)
        upper = collect(reim(root)) .+ rand(2)
        solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
          gridsize, xtol_abs=0.0, xtol_rel=eps(), stopvalroot=eps())
        @test !isempty(solutions)
        for (s, reason) ∈ solutions
          @test reason == :XTOL_REACHED || reason == :STOPVAL_ROOT_REACHED
          reason == :STOPVAL_ROOT_REACHED && continue
          difference = sqrt(sum(centroid(s) - [real(root), imag(root)]).^2)
          @test isapprox(difference, 0, atol=eps(), rtol=0)
          @test isapprox(centroid(s)[1], real(root), atol=0, rtol=eps())
          @test isapprox(centroid(s)[2], imag(root), atol=0, rtol=eps())
        end
      end
    end

  end

  @testset "test other API" begin
    for i in 1:1
      xtol_abs=10.0^(-rand(3:15))
      function mock(x::Vector, root)
        return 1 / ((x[1] + im * x[2]) - root)
      end
      root = rand(ComplexF64)
      objective(x) = mock(x, root)
      lower = collect(reim(root)) .- rand(2)
      upper = collect(reim(root)) .+ rand(2)
      solutions = TriangulatedWindingNumbers.solve(objective, lower, upper,
        1, xtol_abs=xtol_abs)
      @test length(solutions) == 1
      for (s, reason) ∈ solutions
        @test isapprox(centroid(s)[1], real(root), atol=xtol_abs)
        @test isapprox(centroid(s)[2], imag(root), atol=xtol_abs)
      end
    end

  end

  @testset "check errors are caught" begin
    @test_throws ArgumentError TriangulatedWindingNumbers.solve(x->im, [0.0], [1.0],
                                                                [1, 2])
    @test_throws ArgumentError TriangulatedWindingNumbers.solve(x->im, [0.0], [1.0],
                                                                [0])
    @test_throws ArgumentError TriangulatedWindingNumbers.solve(x->im, [0.0], [0.0],
                                                                [1])
    @test_throws ArgumentError TriangulatedWindingNumbers.solve(x->im, [0.0, 0.0],
      [1.0, 1.0], [2, 2], xtol_abs=0.0, xtol_rel=0.0)
  end

end


text:
"""
This file is used to test the simulation. This is just the simulate_states code unwrapped
so will be a little less optimal than when directly run but should enable an easier
observation of correctness.
"""

using Revise
using Distributed
using Plots
using Chain
using ProfileView

include("simulate_states.jl")

# parameters to pass to the main function

file_date = "2022-04-12"

# set seed for consistent plots
rng = Random.Xoshiro(2022)

jurisdiction_assumptions = JurisdictionAssumptions()

# get the latest onset date
latest_start_date = Dates.Date(
    maximum(v for v in values(jurisdiction_assumptions.simulation_start_dates))
)

omicron_start_date = "2021-11-15"
omicron_dominant_date = "2021-12-15"
(dates, local_case_dict, import_case_dict) = read_in_cases(file_date, rng)
last_date_in_data = dates[end]
forecast_end_date = last_date_in_data + Dates.Day(35)
# create vector for dates
onset_dates = latest_start_date:Dates.Day(1):forecast_end_date
# add a small truncation to the simulations as we don't trust the most recent data
truncation_days = 7
# states to simulate
state = "QLD"
nsims = 1000
forecast_start_date = Dates.Date(jurisdiction_assumptions.simulation_start_dates[state])

# named tuple for initial conditions
D0 = jurisdiction_assumptions.initial_conditions[state]
N = jurisdiction_assumptions.pop_sizes[state]

# get the observed cases
local_cases = local_case_dict[state][dates .>= forecast_start_date]
# cutoff the last bit of the local cases
import_cases = import_case_dict[state]
local_cases = Int.(local_cases[begin:end-truncation_days+1])
import_cases = Int.(import_cases)

local_cases_input = CSV.read("data/local_cases_input_2022-04-13.csv", DataFrame)
local_cases_input = filter(:state => ==(state), local_cases_input)

plot(dates, local_case_dict[state])
plot!(local_cases_input.date_onset, local_cases_input.count)
xlims!((Dates.Date("2022-01-01"), Dates.Date("2022-04-15")))

include("simulate_states.jl")

(Z, D, U) = simulate_branching_process(
    D0,
    N,
    nsims,
    local_cases,
    import_cases,
    forecast_start_date,
    file_date,
    state,
    adjust_TP = false,
)

save_simulations(
    D,
    state,
    file_date,
    onset_dates,
    rng,
)

let
    forecast_start_date = Dates.Date(jurisdiction_assumptions.simulation_start_dates[state])
    local_cases = local_case_dict[state][dates .>= forecast_start_date]
    local_cases = local_cases[begin:end]
    D_local = D[:, 1, :] + D[:, 2, :]
    D_local_median = median(D_local, dims = 2)
    D_local_mean = mean(D_local, dims = 2)

    d_start = Dates.Date(forecast_start_date)
    Δd = Dates.Day(1)
    d_end = Dates.Date(forecast_start_date) + Δd * (size(D_local, 1) - 1)
    D_local_dates = d_start:Δd:d_end

    f = plot()
    # plot!(f, dates[dates .>= forecast_start_date], local_cases, legend = false, linealpha = 1, lc = 1)
    # plot!(f, Z_historical[36:end, :], legend = false, lc = 3, linealpha = 0.5)
    plot!(f, D_local_dates, D_local, legend = false, lc = 1, linealpha = 0.5)
    # plot!(f, sims_old[!, "onset date"], sims_old_mat, legend = false, lc = 2, linealpha = 0.5)
    # plot!(f, dates[dates .>= forecast_start_date][begin:end-truncation_days+1], local_cases[begin:end-truncation_days+1], legend = false, lc = :black, lw = 2)
    plot!(f, dates[dates .>= forecast_start_date], local_cases, legend = false, lc = :black, lw = 2)
    # plot!(f, local_cases2.date_onset, local_cases2.count, legend = false, linealpha = 1, lc = 1)
    plot!(f, D_local_dates, D_local_median, legend = false, lc = 2, lw = 2)
    # plot!(f, df_state[!, "onset date"], Matrix(Int.(df_state[!, 3:end-1])))
    vline!(f, [Dates.Date("2021-11-15")], ls = :dash, lc = :black)
    # xlims!(f, 0, length(local_cases) + 35)
    ylims!(f, 0, 40000)
    xlims!(f, Dates.value(Dates.Date("2021-12-15")), Dates.value(Dates.Date("2022-05-20")))
    # ylims!(f, 0, 500)
end

cases_adjusted = CSV.read("results/EpyReff/cases_adjusted_2022-04-12.csv", DataFrame)

cases_adjusted_state = filter(:STATE => ==(state), cases_adjusted)

plot(cases_adjusted_state.local_scaled_14)
plot!(cases_adjusted_state.local_scaled_28)
plot!(cases_adjusted_state.local_scaled_42)
plot!(cases_adjusted_state.local_scaled_72)
plot!(cases_adjusted_state.local_scaled_102)
plot!(cases_adjusted_state.local_scaled_132)

ϕ = read_in_susceptible_depletion(file_date)
ϕ_mean = mean(ϕ, dims = 1)

ψ = (
    ϕ_mean[6] * (cases_adjusted_state.local_scaled_14) / N
    + ϕ_mean[5] * (cases_adjusted_state.local_scaled_28) / N
    + ϕ_mean[4] * (cases_adjusted_state.local_scaled_42) / N
    + ϕ_mean[3] * (cases_adjusted_state.local_scaled_72) / N
    + ϕ_mean[2] * (cases_adjusted_state.local_scaled_102) / N
    + ϕ_mean[1] * (cases_adjusted_state.local_scaled_132) / N
)

scaling_factor = 1 .- ψ

scaling_factor_sim = zeros(size(Z, 1))
prop_inf_full = zeros(size(Z, 1), 6)

Z_historical = mean(Z, dims = 2)

for day in 1:length(Z_historical)

    ψ = 0.0

    τ_horizons_start = [0, 14, 28, 42, 72, 102]
    τ_horizons_end = [1, 14, 28, 42, 72, 102, 132] .- 1

    for i in 1:length(τ_horizons_start)

        # Get starting and ending index. Start must be at least the beginning of the array and end
        # must be at least the forecast length.
        τ_start = max(day - τ_horizons_end[i + 1], 1)
        τ_end = min(day - τ_horizons_start[i], size(Z_historical, 1))

        # sum infections over the time period of interest
        infections = sum(@view Z_historical[τ_start:τ_end])

        # we index from 4 - 1 as the largest reduction occurs as a result of recent infections
        prop_inf = min(1.0, infections / N)
        ψ += ϕ_mean[7 - i] * prop_inf

        prop_inf_full[day, i] = prop_inf

    end

    scale_factor = max(0.0, 1.0 - ψ)

    scaling_factor_sim[day] = scale_factor

end

d_start = Dates.Date(forecast_start_date) - Dates.Day(30)
Δd = Dates.Day(1)
d_end = Dates.Date(forecast_start_date) + Δd * (size(Z, 1) - 1) - Dates.Day(30)
D_local_dates = d_start:Δd:d_end

plot(cases_adjusted_state.date_inferred, scaling_factor)
plot!(D_local_dates, scaling_factor_sim)
xlims!((Dates.Date("2021-12-01"), D_local_dates[end]))
ylabel!("Scaling factor")

fig = plot(layout = (3, 2), legend = false)
plot!(
    fig,
    subplot = 1,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_14 ./ N,
    title = "14 days", titlefontsize=8
)
plot!(
    fig,
    subplot = 2,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_28 ./ N,
    title = "28-14 days", titlefontsize=8
)
plot!(
    fig,
    subplot = 3,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_42 ./ N,
    title = "42-28 days", titlefontsize=8
)
plot!(
    fig,
    subplot = 4,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_72 ./ N,
    title = "72-42 days", titlefontsize=8
)
plot!(
    fig,
    subplot = 5,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_102 ./ N,
    title = "102-72 days", titlefontsize=8)
plot!(
    fig,
    subplot = 6,
    cases_adjusted_state.date_inferred,
    cases_adjusted_state.local_scaled_132 ./ N,
    title = "132-102 days", titlefontsize=8)
plot!(
    fig,
    subplot = 1,
    D_local_dates,
    prop_inf_full[:, 1],
    ls = :dash
)
plot!(
    fig,
    subplot = 2,
    D_local_dates,
    prop_inf_full[:, 2],
    ls = :dash
)
plot!(
    fig,
    subplot = 3,
    D_local_dates,
    prop_inf_full[:, 3],
    ls = :dash
)
plot!(
    fig,
    subplot = 4,
    D_local_dates,
    prop_inf_full[:, 4],
    ls = :dash
)
plot!(
    fig,
    subplot = 5,
    D_local_dates,
    prop_inf_full[:, 5],
    ls = :dash
)
plot!(
    fig,
    subplot = 6,
    D_local_dates,
    prop_inf_full[:, 6],
    ls = :dash
)
xlims!((Dates.Date("2021-12-01"), D_local_dates[end]))
ylims!((0, 0.25))


text:
module q3

import Evolve

steps = 8000                # number of simulation steps
samplesteps = 50            # steps between measurements
L = 2                       # width of the box
xmin = -1                   # leftmost point in the box
N = 2000                    # number of subintervals
ϵ = L/N                     # width of each subinterval
ks = [200]                  # wavenumbers
σ = 0.05                    # wavepacket width
V0s = [5e3:1000:5e4]         # step potentials
x0 = -0.5                   # wavepacket starting position
Ts = Float64[]
Rs = Float64[]

# Run the trials
for k in ks
    for V0 in V0s
        print("Initializing for k = $k, V0 = $V0... ")
        V(xs) = map(x -> 0.5(1+sign(x))V0, xs)       # Heaviside potential
        AV = Evolve.A(V)            # evolution matrix constructor for square well
        ϕV = Evolve.Φ(k, σ, x0)     # initial state vector constructor for Gaussian packet

        A = AV(L, N, xmin)
        ϕ0 = ϕV(L, N, xmin)
        x = (L/N .* [0:N] + xmin)   # x positions

        # Evolve the system and record positions every so many steps
        print("evolving... ")
        test = x -> (abs(x[2]) > 1e-2 || abs(x[end-1]) > 1e-2)
        ϕ = Evolve.evolve(A, ϕ0, steps, test)
        phis = [ϕ0 ϕ[:,samplesteps:samplesteps:end]]

        # calculate transmission and reflection probabilities
        T = mean(abs(ϕ[N/2+1:end,end]).^2)
        R = mean(abs(ϕ[1:N/2,end]).^2)
        I = T + R
        T /= I
        R /= I
        push!(Ts, T)
        push!(Rs, R)
        print("T = $T, R = $R...")
        ### writedlm("trans_step_$V0", [T,R,I])

        println("done after $(size(ϕ)[2]) steps. Writing to file.")
        ### writedlm("phis_step_im_$k\_$V0.dat", imag(phis))
        ### writedlm("phis_step_re_$k\_$V0.dat", real(phis))
    end
end

writedlm("VTR.dat", [V0s Ts Rs])
end


text:
# Elias Kuepper, Holger Teichgraeber, 2019

 ######################
 # CapacityExpansion
 # Capacity Expansion Probelem formulation
 #
 #####################
#using ClustForOpt
include(joinpath("..","..","ClustForOpt","src","ClustForOpt_development.jl"))
using CSV
using DataFrames
using StatsBase
using JLD2
using FileIO
using JuMP

include(joinpath("utils","datastructs.jl"))
include(joinpath("utils","optvariable.jl"))
include(joinpath("utils","utils.jl"))
include(joinpath("utils","load_data.jl"))
include(joinpath("optim_problems","run_opt.jl"))
include(joinpath("optim_problems","opt_cep.jl"))


text:
"""
    Environment

Module pertinent to a sentential formula including mutable state (for new
variables and functors) and type information.
"""
module Environment

export BoolType, IntType, RationalType, RealType, Env, resetEnv!, userTypes,
  newVariable, newFunctor

import Base: length, show, haskey, setindex!, getindex, get
import FirstOrderLogic: Variable, Functor, FunctionTerm, PredicateTerm

abstract type PrimitiveType end

# the allowed logic types
struct InputType <: PrimitiveType end
show(io::IO, _::InputType) = print(io, "\$i")
struct BoolType <: PrimitiveType end
show(io::IO, _::BoolType) = print(io, "\$o")
struct IntType <: PrimitiveType  # we optionally allow an upper bound
  maxVal::Int
  IntType(d=0) = new(d)    # default to unbounded (indicated by 0)
end
show(io::IO, t::IntType) = t.maxVal>0 ? print(io, "\$int($(t.maxVal))") : print(io, "\$int")
struct RationalType <: PrimitiveType end
show(io::IO, _::RationalType) = print(io, "\$rat")
struct RealType <: PrimitiveType end
# user-defined types are composites of primitive types; function types are
# represented with the last tuple element determining the output and the
# other elements representing the inputs
show(io::IO, _::RealType) = print(io, "\$real")
struct UserType <: PrimitiveType end
show(io::IO, _::UserType) = print(io, "\$tType")

"""
    PrimitiveType(str)

Build a primitive logic type (\$i,\$o,\$int,\$rat,\$real,\$tType) based on the
strings representing TPTP types: `str`. `str` may be i, o, int, rat, real.
"""
function PrimitiveType(s::AbstractString)
  if s == "\$i"
    InputType()
  elseif s == "\$o"
    BoolType()
  elseif s == "\$int"
    IntType()
  elseif s == "\$rat"
    RationalType()
  elseif s == "\$real"
    RealType()
  elseif s == "\$tType"
    UserType()
  else
    s
  end
end
const CompositeType = Tuple{Vararg{Union{String,PrimitiveType}}}

"""
    Environment

Struct representing and environment associated with a TPTP
description. Environments include type information and mutable state used to
generate fresh variables and functors.
"""
mutable struct Env
  varCount::Int                         # unique variable id
  funcCount::Int                        # unique functor id
  userTypes::Dict{String,CompositeType} # user-defined type definitions
end
Env() = Env(0, 0, Dict{String,CompositeType}())

"""
    length(e::Env)

The number of defined symbols in the environment `e`
"""
length(e::Env) = length(e.userTypes)

"""
    resetEnv!(e::Env)

Reset the state and types in an environment `e`.
"""
function resetEnv!(e::Env)
  e.varCount = e.funcCount = 0
  for k in keys(e.userTypes)
    delete!(e.userTypes, k)
  end
end

"""
    setindex!(e::Env, v, k::String)

Record the user type identified by key `k` and type `v` in the environment
`e`. `v` may be a `CompositeType`, a `String` or `Vector{String}` which are
converted to `CompositeType`.
"""
setindex!(e::Env, v::Union{PrimitiveType,String}, k) = setindex!(e.userTypes, (v,), k)
setindex!(e::Env, v::CompositeType, k) = setindex!(e.userTypes, v, k)
setindex!(e::Env, v::AbstractString, k) = setindex!(e.userTypes, (PrimitiveType(v),), k)
setindex!(e::Env, v::Vector, k) = setindex!(e.userTypes, tuple(v...), k)

"""
    getindex(e::Env, k::String)

Get the type associated with key `k`.
"""
getindex(e::Env, fv::Union{Functor,Variable}) = getindex(e.userTypes, fv.name)
getindex(e::Env, fp::Union{FunctionTerm,PredicateTerm}) = getindex(e, fp.name)
getindex(e::Env, s::AbstractString) = getindex(e.userTypes, s)

"""
    get(e::Env, k::String, default)

Get the type associated with key `k` and if not found return default.
"""
get(e::Env, fv::Union{Functor,Variable}, default) = get(e.userTypes, fv.name, default)
get(e::Env, fp::Union{FunctionTerm,PredicateTerm}, default) = get(e, fp.name, default)
get(e::Env, s::AbstractString, default) = get(e.userTypes, s, default)

"""
    haskey(e::Env, k)

Check whether there is a type associated with key `k`.
"""
haskey(e::Env, fv::Union{Functor,Variable}) = haskey(e.userTypes, fv.name)
haskey(e::Env, fp::Union{FunctionTerm,PredicateTerm}) = haskey(e, fp.name)
haskey(e::Env, s::AbstractString) = haskey(e.userTypes, s)

"""
    userTypes(e::Env)

Return the user types in the environment `e`.
"""
userTypes(e::Env) = e.userTypes

"""
    newVariable(vf, e::Env)

Create a fresh variable whose name is based on `vf` (either a `Variable` or
`Functor`) in the environment e.
"""
newVariable(v::Variable, e::Env) = Variable(v.name*"##"*string(e.varCount += 1))
newVariable(f::Functor, e::Env) = Functor(f.name * "#sk" * string(e.funcCount += 1))

"""
    newFunctor(v, e::Env), newFunctor(e::Env)

Create a fresh `Functor` whose name is based on `Variable` `vf` or using a
standard name in the environment e.
"""
newFunctor(v::Variable, e::Env) = Functor(lowercase(v.name)*"#sk"*string(e.funcCount += 1))
newFunctor(e::Env) = Functor("#f#"*string(e.funcCount += 1))

end # module Environment


text:
#!/usr/bin/env julia

using Luxor, Colors

function test_ellipse(fname)
  pagewidth, pageheight = 1200, 1400
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background("ivory")
  setopacity(0.5)
  setline(0.3)

  pagetiles = Tiler(pagewidth, pageheight, 4, 5, margin=50)
  for (pos, n) in pagetiles
    randomhue()
    ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :strokepreserve)
    clip()
    for i in 1:10
      ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :fill)
      randomhue()
      ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :stroke)
    end
    fill()
    sethue("black")
    ellipse(pos, 5, 5, :fill)
    clipreset()
  end

  finish()
  println("ellipse-test saved in $(fname)")
end

fname = "/tmp/ellipse-test1.pdf"
test_ellipse(fname)


text:
################################################################################
#
# Helmholtz-related linear algebra routines
#
################################################################################

export HelmholtzLinOp, computeNabc

function HelmholtzLinOp(h::NTuple{2, R}, m::Array{R, 2}, freq::R, nabc::NTuple{4, Int64})

    # Extension for ABC
    m_ext = extend(m, nabc)
    n_ext = size(m_ext)

    # Complex frequency for ABC
    abc_fact = 1f0
    ω = 2f0*pi*freq*(1f0.-abc_fact*1im*taper(size(m), nabc))

    # Laplacian
    Dxx = spdiagm(-1 => ones(R, n_ext[1]-1), 0 => -2f0*ones(R, n_ext[1]), 1 => ones(R, n_ext[1]-1))/h[1]^2
    Dxx[1, 1:2] = [1 -1]/h[1]^2 # Sommerfeld radiation condition
    Dxx[end, end-1:end] = [-1, 1]/h[1]^2 # Sommerfeld radiation condition
    Dzz = spdiagm(-1 => ones(R, n_ext[2]-1), 0 => -2f0*ones(R, n_ext[2]), 1 => ones(R, n_ext[2]-1))/h[2]^2
    Dzz[1, 1:2] = [1 -1]/h[2]^2 # Sommerfeld radiation condition
    Dzz[end, end-1:end] = [-1, 1]/h[2]^2 # Sommerfeld radiation condition
    Δ = kron(sparse(I, n_ext[2], n_ext[2]), Dxx)+kron(Dzz, sparse(I, n_ext[1], n_ext[1]))

    # Mass matrix
    w = ones(R, n_ext)
    w[[1, end], :] .= 0f0
    w[:, [1, end]] .= 0f0
    v = 1f0.-w
    v[:, [1, end]] = v[:, [1, end]]/h[2] # Sommerfeld radiation condition
    v[[1, end], :] = v[[1, end], :]/h[1] # Sommerfeld radiation condition
    # M = spdiagm(0 => vec(ω.^2f0.*w.*m_ext))+1im*spdiagm(0 => vec(ω.*v.*sqrt.(m_ext)))
    M = spdiagm(0 => vec(ω.^2f0.*w.*m_ext))+1im*spdiagm(0 => vec(ω.*v.*sqrt.(C.(m_ext))))

    # Assembly
    return -M-Δ

end

function taper(n::Tuple{Int64, Int64}, nabc::NTuple{4, Int64})

    taper_x = 1f0.-cat([range(-1f0, stop = 0f0, length = nabc[1])...].^2, zeros(R, n[1]), [range(0f0, stop = 1f0, length = nabc[2])...].^2, dims = 1)
    taper_z = 1f0.-cat([range(-1f0, stop = 0f0, length = nabc[3])...].^2, zeros(R, n[2]), [range(0f0, stop = 1f0, length = nabc[4])...].^2, dims = 1)
    return 1f0.-reshape(taper_x, :, 1)*reshape(taper_z, 1, :)

end

function computeNabc(h::NTuple{2, R}, m::Array{R, 2}, freq::R)

    return Int64.((round(1f0./(min(h...)*min(sqrt.(m[1, :])...)*freq)), round(1f0./(min(h...)*min(sqrt.(m[end, :])...)*freq)), round(1f0./(min(h...)*min(sqrt.(m[:, 1])...)*freq)), round(1f0./(min(h...)*min(sqrt.(m[:, end])...)*freq))))

end


text:
function hl_type(spec::SpecStructMember)
    @match s = spec begin
        if s.name == :pNext end => :Any
        GuardBy(is_version) => :VersionNumber
        GuardBy(is_arr) => begin
            T = hl_type(ptr_type(s.type))
            :(Vector{$T})
        end
        _ => hl_type(spec.type)
    end
end

function hl_type(type)
    @match t = type begin
        :Cstring => :String
        :VkBool32 => :Bool
        GuardBy(is_opaque_pointer) => t
        :(NTuple{$N,$T}) => begin
            @match T begin
                :UInt8 => :String
                # :(Ntuple{$M,$T}) => :(SMatrix{$N,$M,$T})
                # _ => :(SVector{$N,$T})
                _ => :(NTuple{$N,$(hl_type(T))})
            end
        end
        GuardBy(in([spec_structs.name; spec_unions.name])) => struct_name(t, true)
        GuardBy(in(spec_handles.name)) => remove_vk_prefix(t)
        GuardBy(is_fn_ptr) => :FunctionPtr
        :(Ptr{$T}) => hl_type(T)
        GuardBy(is_flag_bitmask) => bitmask_flag_type(t)
        GuardBy(in(spec_flags.name)) && if !isnothing(flag_by_name(t).bitmask)
        end => bitmask_flag_type(flag_by_name(t).bitmask)
        GuardBy(in(spec_constants.name)) => follow_constant(t)
        GuardBy(in(spec_enums.name)) => enum_type(t)
        GuardBy(is_vulkan_type) => remove_vk_prefix(t)
        _ => t
    end
end

function idiomatic_julia_type(type)
    @match t = type begin
        GuardBy(is_fn_ptr) => :FunctionPtr
        GuardBy(is_opaque_pointer) => t
        :(NTuple{$N,UInt8}) => :String
        :(NTuple{$N,$T}) => begin
            _N = @match N begin
                ::Symbol => :(Int($N))
                _ => N
            end
            :(NTuple{$_N,$(idiomatic_julia_type(T))})
        end
        :Cstring => :String
        :VkBool32 => :Bool
        :(Ptr{$pt}) => idiomatic_julia_type(pt)
        GuardBy(in([spec_structs.name; spec_unions.name])) => struct_name(t)
        GuardBy(in(spec_handles.name)) => remove_vk_prefix(t)
        GuardBy(in(spec_enums.name)) => enum_type(t)
        GuardBy(is_flag_bitmask) => bitmask_flag_type(t)
        GuardBy(in(spec_flags.name)) && if !isnothing(flag_by_name(t).bitmask)
        end => bitmask_flag_type(flag_by_name(t).bitmask)
        GuardBy(in(spec_constants.name)) => follow_constant(t)
        _ => t
    end
end

"""
Return a new type easier to deal with.
"""
function idiomatic_julia_type(spec::Spec)
    @match s = spec begin
        GuardBy(is_version) => :VersionNumber
        GuardBy(is_arr) => :(Vector{$(idiomatic_julia_type(ptr_type(s.type)))})
        GuardBy(is_data) => :(Ptr{Cvoid})
        _ => idiomatic_julia_type(s.type)
    end
end

function signature_type(type)
    @match type begin
        :UInt || :UInt8 || :UInt16 || :UInt32 || :UInt64 || :Int || :Int8 || :Int16 || :Int32 || :Int64 => :Integer
        :Float16 || :Float32 || :Float64 => :Real
        :String => :AbstractString
        :(Vector{$et}) => begin
            @match st = signature_type(et) begin
                :AbstractString || :Integer || :Real => :(AbstractArray{<:$st})
                _ => :(AbstractArray{$st})
            end
        end
        t => t
    end
end

function relax_signature_type(type)
    @match type begin
        :(AbstractArray{$_}) || :(AbstractArray{<:$_}) => :AbstractArray
        t => t
    end
end

function relax_function_signature!(f, p::Dict)
    args = map(p[:args]) do arg
        if f(arg)
            relax_function_signature(arg)
        else
            arg
        end
    end
    p[:args] = args
end

function relax_function_signature(arg)
    @match arg begin
        :($identifier::$type) => :($identifier::$(relax_signature_type(type)))
        _ => arg
    end
end

bitmask_flag_type(type) = Symbol(replace(remove_vk_prefix(string(type)), "Bits" => ""))
bitmask_flag_type(spec::SpecBitmask) = bitmask_flag_type(spec.name)


text:
module NeuralStyles

using Metalhead
using Images
using Optim
using Flux
#using CuArrays
import Flux: @epochs
#CuArrays.allowscalar(false)

include("pre.jl")
include("styles.jl")

export train

end # module


text:
module Transforms

include("distancetransforms.jl")

include("statetransforms.jl")

# linear basis transform

include("compose.jl")

end


text:
"""
MPResults

This module has the code you need for the mixed precision example.
"""
module MPResults

export nsold01
export heqtest
export plotnsold
export fprintTeX
export data_populate
export data_harvest

using PyPlot
using LinearAlgebra
using JLD2
using Printf

include("Heq4nsold.jl")
include("nsoldv01.jl")
include("heqtest.jl")
include("plotnsold.jl")
include("fprintTeX.jl")
include("data_populate.jl")
include("data_harvest.jl")
#include("tabulatensold.jl")

using .Heq4nsold

end

max_stars_repo_name: ctkelley/MPResults


text:
import SimpleValueGraphs: typetuple, OneEdgeValGraph

@testset "edge iterator $G" for
    G in (ValGraph, ValOutDiGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
              edgeval_types=typetuple(E_VALS),
              edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        iter = edges(g, :)

        E_should_be = (is_directed(g) ? ValDiEdge : ValEdge){V, E_VALS}

        @testset "eltype" begin
            @test eltype(iter) <: E_should_be
        end

        @testset "length" begin
            @test length(iter) == ne(g)
        end

        @testset "correct edges" begin
            @test all(iter) do e
                s, d, vs = src(e), dst(e), get_edgeval(e, :)
                e isa E_should_be &&
                has_edge(g, s, d) &&
                get_edgeval(g, s, d, :) == vs
            end
        end

        @testset "lexically strictly increasing" begin
            s_prev, d_prev = 0, 0
            @test all(iter) do (e)
                s, d = src(e), dst(e)
                if (s < s_prev ||  (s == s_prev && d <= d_prev))
                    return false
                end
                s_prev, d_prev = s, d
                return true
            end
        end

    end
end

@testset "inneighbors $G" for
    G in (ValGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
            edgeval_types=typetuple(E_VALS),
            edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        @testset "inneighbors are unique" begin
            @test all(vertices(g)) do u
                allunique(inneighbors(g, u))
            end
        end

        @testset "inneighbors same as simple graph" begin
            @test all(vertices(g)) do u
                issetequal(inneighbors(g, u), inneighbors(gs, u))
            end
        end

    end
end

@testset "outneighbors $G" for
    G in (ValGraph, ValOutDiGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
            edgeval_types=typetuple(E_VALS),
            edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        @testset "outneighbors are unique" begin
            @test all(vertices(g)) do u
                allunique(outneighbors(g, u))
            end
        end

        @testset "outneighbors same as simple graph" begin
            @test all(vertices(g)) do u
                issetequal(outneighbors(g, u), outneighbors(gs, u))
            end
        end

    end
end

@testset "outedgevals $G" for
    G in (ValGraph, ValOutDiGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
            edgeval_types=typetuple(E_VALS),
            edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        @testset "outedgevals are correct for key $key" for key in allkeys_for_E_VALS(E_VALS)

            @test all(vertices(g)) do u
                outedgevals(g, u, key) == [get_edgeval(g, u, v, key) for v in outneighbors(g, u)]
            end
        end

        if g isa OneEdgeValGraph
            @testset "outedgevals for graphs with single key" begin
                @test all(vertices(g)) do u
                    outedgevals(g, u) == outedgevals(g, u, 1)
                end
            end
        end
    end
end

@testset "inedgevals $G" for
    G in (ValGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
            edgeval_types=typetuple(E_VALS),
            edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        @testset "inedgevals are correct for key $key" for key in allkeys_for_E_VALS(E_VALS)

            @test all(vertices(g)) do u
                inedgevals(g, u, key) == [get_edgeval(g, v, u, key) for v in inneighbors(g, u)]
            end
        end

        if g isa OneEdgeValGraph
            @testset "inedgevals for graphs with single key" begin
                @test all(vertices(g)) do u
                    inedgevals(g, u) == inedgevals(g, u, 1)
                end
            end
        end
    end
end

@testset "edges $G" for
    G in (ValGraph, ValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL,
    (gs, info) in make_testgraphs(is_directed(G) ? SimpleDiGraph{V} : SimpleGraph{V})

    @testset "Params: $G{$V, $E_VALS}, gs = $info" begin

        g = G(gs;
            edgeval_types=typetuple(E_VALS),
            edgeval_init=(s,d) -> rand_sample(E_VALS)
        )

        @testset "edges are lexicographical sorted, unique" begin
            lexicographical_sorted = true
            s = 0
            d = 0
            for e in edges(g, :)
                s_new = src(e)
                d_new = dst(e)

                if s_new < s || (s_new == s && d_new <= d)
                    lexicographical_sorted = false
                    break
                end
            end
            @test lexicographical_sorted
        end

        @testset "correct src, dst" begin
            @test all(edges(g, :)) do e
                has_edge(g, src(e), dst(e))
            end
        end

        @testset "correct values" begin
            @test all(edges(g, :)) do e
                get_edgeval(e, :) == get_edgeval(g, src(e), dst(e), :)
            end
        end

        @testset "correct length" begin
            @test length(edges(g, :)) == ne(g)
        end
    end
end


text:
function solve_stokes_threaded(
    color_list,
    KS,
    GS,
    MS,
    FS,
    TT,
    DoF_U,
    DoF_P,
    coordinates,
    ScratchStokes::ScratchThermalGlobal,
    SF_Stokes::ShapeFunctionsThermal,
    U,
    Upolar,
    Ucartesian,
    P,
    g,
    ρ,
    η,
    𝓒,
    UBC,
    to,
)

    # unpack boundary conditions
    ∂Ωu = UBC.Ω
    ufix = UBC.vfix
    ufree =  UBC.ifree
    # unpack fixed element coordinates
    θ, r = coordinates.θ, coordinates.r

    @timeit to "Thermal diffusion threaded" begin
        # Reset Matrices
        fill!(KS.nzval, 0.0)
        fill!(GS.nzval, 0.0)
        fill!(MS.nzval, 0.0)
        fill!(FS, 0.0)

        ## Element assembly
        @timeit to "Assembly" stokes_assembly_threaded!(
            color_list,
            KS,
            MS,
            GS,
            FS,
            DoF_U,
            DoF_P,
            coordinates,
            ScratchStokes,
            SF_Stokes,
            g,
            ρ,
            η,
            𝓒,
        )

        @timeit to "BCs" begin
            KS, GS, FS = _prepare_matrices(KS, GS, FS, TT)
            U, FS = _apply_bcs(U, KS, FS, ∂Ωu, ufix)
        end

        @timeit to "PCG solver" U,P = StokesPcCG(U, P, KS, MS, GS, FS, ufree)

        @timeit to "Remove net rotation" U, Ucart, Upolar, Ucartesian = updatevelocity2(U, Ucartesian, Upolar, ρ, TT, coordinates, gr)

    end
    return T, to
end

"""
assembly_threaded!(color_list, K, M, F, dof, coordinates::ElementCoordinates,
VarT::ThermalParameters, A::ScratchThermalGlobal,
B::ShapeFunctionsThermal, valA, ρ, Δt)

Threaded assembly of FEM matrices for thermal diffusion. Based on a colored mesh.
"""
function stokes_assembly_threaded!(
    color_list,
    KS,
    MS,
    GS,
    FS,
    DoF_U,
    DoF_P,
    coordinates::ElementCoordinates,
    ScratchStokes::ScratchThermalGlobal,
    SF_Stokes::ShapeFunctionsStokes,
    g,
    ρ,
    η,
    𝓒,
)

    for color in color_list
        Threads.@threads for element in color
            stokes_assemble_element!(
                element,
                KS,
                MS,
                GS,
                FS,
                DoF_U,
                DoF_P,
                coordinates.θ,
                coordinates.r,
                ScratchStokes,
                SF_Stokes,
                g,
                ρ,
                η,
                𝓒,
            )
        end
    end

    return KS, GS, MS, FS
end

mydet(A::SMatrix{2, 2, T, 4}) where {T} = A[1] * A[4] - A[2] * A[3]

function stokes_assemble_element!(
    iel,
    KS,
    MS,
    GS,
    FS,
    DoF_U,
    DoF_P,
    θ,
    r,
    A::ScratchThermalGlobal,
    SF_Stokes::ShapeFunctionsStokes,
    g,
    ρ,
    η,
    𝓒,
)

    N, ∇N, NP, dN3ds, w_ip, N3 =
        SF_Stokes.N, SF_Stokes.∇N, SF_Stokes.NP, SF_Stokes.dN3ds, SF_Stokes.w_ip, SF_Stokes.N3

    valU, valP = Val(12), Val(3)

    # Polar coordinates of element nodes
    θ_el = @SVector [θ[i, iel] for i in 1:3]
    r_el = @SVector [r[i, iel] for i in 1:3]
    coords = SMatrix{3,2}([θ_el r_el])

    # Jacobian n. 1 (p:=polar, l:=local): reference element --> current element
    J_pl = dN3ds * coords
    detJ_pl = mydet(J_pl)

    R_21 = r_el[2] - r_el[1]  # = -detJa_PL*deta_dth
    R_31 = r_el[3] - r_el[1]  # =  detJa_PL*dxi_dth
    Th_31 = θ_el[3] - θ_el[1] # = -detJa_PL*dxi_dr
    Th_21 = θ_el[2] - θ_el[1] # =  detJa_PL*deta_dr

    # the Jacobian ∂ξ∂θ to transform local (ξ, η) into global (θ,r) derivatives
    #     ∂ξ∂θ = [ R_31    -R_21
    #             -Th_31    Th_21] / detJa_PL
    # ∂ξ∂θ = @SMatrix [(r_el[3]-r_el[1]) (r_el[2]-r_el[1])
    #                  (θ_el[3]-θ_el[1])  (θ_el[2]-θ_el[1])]

    Ke = empty_element_matrices(valU)
    Ge = empty_element_matrices(valU, valP)
    Me = empty_element_matrices(valP)
    fe = empty_element_force_vector(valU)

    Pdofs = DoF_P[iel]
    Udofs = DoF_U[iel]
    ρ_el = @SVector [ρ[Pdofs[i]] for i in 1:3]

    # INTEGRATION LOOP
    @inbounds for ip in 1:(A.nip)
        # Unpack shape functions
        N_ip = N[ip]
        NP_ip = NP[ip]
        N3_ip = N3[ip]
        ∇N_ip = ∇N[ip]

        # ρ at ith integration point
        ρ_ip = mydot(ρ_el, N3_ip)
        η_ip = _element_viscosity(η, gr.e2n, PhaseID, iel, N3_ip)

        # Polar coordinates of the integration points
        θ_ip = mydot(θ_el, N3_ip)
        r_ip = mydot(r_el, N3_ip)
        cos_ip = cos(θ_ip)
        sin_ip = sin(θ_ip)
        cos_ip_r_ip = cos_ip / r_ip
        sin_ip_r_ip = sin_ip / r_ip
        # Build inverse of the 2nd Jacobian
        # invJ_double = @SMatrix [(∂ξ∂θ[1,1]*cos_ip_r_ip-∂ξ∂θ[2,1]*sin_ip) (-∂ξ∂θ[1,2]*cos_ip_r_ip+∂ξ∂θ[2,2]*sin_ip);
        #                        -(∂ξ∂θ[1,1]*sin_ip_r_ip+∂ξ∂θ[2,1]*cos_ip) (∂ξ∂θ[1,2]*sin_ip_r_ip+∂ξ∂θ[2,2]*cos_ip) ]
        invJ_double = @SMatrix [
             R_31*cos_ip_r_ip-Th_31*sin_ip   -R_21*cos_ip_r_ip+Th_21*sin_ip
            -R_31*sin_ip_r_ip-Th_31*cos_ip    R_21*sin_ip_r_ip+Th_21*cos_ip
        ]

        # Partial derivatives
        ∂N∂x = invJ_double * ∇N_ip / detJ_pl
        # Integration weight
        ω = r_ip * detJ_pl * w_ip[ip]

        # Update elemental matrices
        D = DMatrix(𝓒, iel, ip, Val(η))
        B = @SMatrix [
            ∂N∂x[1,1]   0           ∂N∂x[1,2]   0           ∂N∂x[1,3]   0           ∂N∂x[1,4]   0           ∂N∂x[1,5]   0           ∂N∂x[1,6]   0
            0           ∂N∂x[2,1]   0           ∂N∂x[2,2]   0           ∂N∂x[2,3]   0           ∂N∂x[2,4]   0           ∂N∂x[2,5]   0           ∂N∂x[2,6]
            ∂N∂x[2,1]   ∂N∂x[1,1]   ∂N∂x[2,2]   ∂N∂x[1,2]   ∂N∂x[2,3]   ∂N∂x[1,3]   ∂N∂x[2,4]   ∂N∂x[1,4]   ∂N∂x[2,5]   ∂N∂x[1,5]   ∂N∂x[2,6]   ∂N∂x[1,6]
        ]

        # Element matrices
        Ke += (B' * D * B) * ω
        Ge += (vec(∂N∂x).*NP_ip) * ω
        Me += (NP_ip.*NP_ip') * (ω/η_ip)

        # Force vector -- right hand side
        dummy1 = @SVector [sin_ip, cos_ip, sin_ip, cos_ip, sin_ip, cos_ip, sin_ip, cos_ip, sin_ip, cos_ip, sin_ip, cos_ip]
        dummy2 = @SVector [N_ip[1], N_ip[1], N_ip[2], N_ip[2], N_ip[3], N_ip[3], N_ip[4], N_ip[4], N_ip[5], N_ip[5], N_ip[6], N_ip[6]]
        fe +=  (g * ω * ρ_ip * cos_ip) * (dummy2 .* dummy1)

    end

    # Update global matrices
    update_KS!(KS, Udofs, Ke)
    update_KS!(MS, Pdofs, Me)
    update_GS!(GS, Udofs, Pdofs, Ge)
    update_force_vector!(FS, fe, Udofs)

end

function update_KS!(
    K::SparseMatrixCSC{Float64,Int64}, Udofs, Ke
)
    for (_j, j) in enumerate(Udofs), (_i, i) in enumerate(Udofs)
         @inbounds @fastmath K[i, j] += Ke[_i, _j]
    end
end

function update_GS!(
    G::SparseMatrixCSC{Float64,Int64}, Udofs, Pdofs, Ge
)
    for (_j, j) in enumerate(Pdofs), (_i, i) in enumerate(Udofs)
         @inbounds @fastmath G[i, j] += Ge[_i, _j]
    end
end

function update_MS!(
    M::SparseMatrixCSC{Float64,Int64}, Pdofs, Me
)
    for (_j, j) in enumerate(Pdofs), (_i, i) in enumerate(Pdofs)
         @inbounds @fastmath M[i, j] += Me[_i, _j]
    end
end

function stokes_immutables(gr, nnod, ndim, nvert, nnodel, nel, nip)
    # Constant mesh values
    ScratchStokes = ScratchThermalGlobal(ndim, # := model dimension
                                         nvert, # := number of vertices per element
                                         nnodel, # := number of nodes per element
                                         nel, # := number of elements
                                         nip) # := number of vertices per element

    ni, nn, nnP, nn3 = Val(nip), Val(nnodel), Val(3), Val(3)
    N, ∇N, _, w_ip = _get_SF(ni,nn)
    NP,_,_,_ = _get_SF(ni,nnP)
    N3,_, dN3ds,_ = _get_SF(ni,nn3)

    # Shape function and their derivatives
    SF_Stokes = ShapeFunctionsStokes( N, ∇N, NP, dN3ds, w_ip, N3)

    # Degrees of freedom per element
    DoF_U = DoFs_Velocity(gr.e2n)
    DoF_P = DoFs_Pressure(gr.e2nP)
    # Diffusion sparse matrices and force vector
    _, _, _, _, KS, GS, MS, _ = sparsitystokes(gr)
    FS = Vector{Float64}(undef, 2*nnod)

    return  KS, GS, MS, FS, DoF_U, DoF_P, nn, SF_Stokes, ScratchStokes
end

function stress_shape_functions(; nip = 6)
    ni, nn, nnP, nn3 = Val(nip), Val(6), Val(3), Val(3)
    N, ∇N, _, w_ip = _get_SF(ni,nn)
    NP,_,_,_ = _get_SF(ni,nnP)
    N3,_, dN3ds,_ = _get_SF(ni,nn3)

    # Shape function and their derivatives
    ShapeFunctionsStress( N, ∇N, NP, dN3ds, w_ip, N3)

end

DMatrix(𝓒, iel, ip, ::Val{Anisotropic}) =
    @SMatrix [𝓒.η11[iel,ip] 𝓒.η13[iel,ip] 𝓒.η15[iel,ip]
              𝓒.η13[iel,ip] 𝓒.η33[iel,ip] 𝓒.η35[iel,ip]
              𝓒.η15[iel,ip] 𝓒.η35[iel,ip] 𝓒.η55[iel,ip]]

DMatrix(𝓒, iel, ip, ::Val{Isotropic}) =
    @SMatrix [4/3  -2/3 0
             -2/3   4/3 0
              0     0   1]

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
using ALifeBenchmark
using Printf
using Serialization
using Dates
import Random
using Profile
using PProf

const NUM_ITERATIONS = 1e6
const NUM_TRIALS = 1

create_model(logger) = GebModel(logger)
function print_function(trial, step, model)
    if step % 500 == 0
        @info "$step \t $(step / NUM_ITERATIONS) \t $(length(model.organisms))"
    end
end

collect_distribution(create_model, execute!, print_function, NUM_ITERATIONS, NUM_TRIALS, 1)

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
@testset "Cholesky" begin
    # Check sensitivities for lower-triangular version.
    let rng = MersenneTwister(123456), N = 10
        for _ in 1:10
            B, VB = randn.(Ref(rng), [N, N], [N, N])
            Ū = UpperTriangular(randn(rng, N, N))
            @test check_errs(B->cholesky(B'B + I).U, Ū, B, 1e-2 .* VB)
        end
    end

    let
        X_ = Matrix{Float64}(I, 5, 5)
        X = Leaf(Tape(), X_)
        C = cholesky(X)
        @test C isa Branch{<:Cholesky}
        @test getfield(C, :f) == LinearAlgebra.cholesky
        U = C.U
        @test U isa Branch{<:UpperTriangular}
        @test getfield(U, :f) == Base.getproperty
        @test unbox(U) ≈ cholesky(X_).U

        @test_throws MethodError ∇(X->cholesky(X).info)(X_)
    end

    let
        X_ = Matrix{Float64}(I, 5, 5)
        X = Leaf(Tape(), X_)
        U = cholesky(X).U
        C = Cholesky(U, 'U', 0)
        @test C isa Branch{<:Cholesky}
        @test getfield(C, :f) == LinearAlgebra.Cholesky
        @test unbox(C) == Cholesky(UpperTriangular(X_), 'U', 0)
        # Ensure we can still directly extract the .U field
        UU = C.U
        @test UU isa Branch{<:UpperTriangular}
        # And access .L as well
        LL = C.L
        @test LL isa Branch{<:LowerTriangular}
        # Make sure that computing the Cholesky and already having the Cholesky
        # produce the same results
        expected = Matrix(0.5I, 5, 5)
        @test ∇(X->det(cholesky(X).U))(X_)[1] ≈ expected
        @test ∇(X->det(cholesky(X).L))(X_)[1] ≈ expected
        @test ∇(X->det(Cholesky(cholesky(X).U, :U, 0).U))(X_)[1] ≈ expected
        @test ∇(X->det(Cholesky(cholesky(X).L, 'L', 0).U))(X_)[1] ≈ expected
    end
end


text:
using EzXML

# Tests for TiffData adapted from
# https://docs.openmicroscopy.org/ome-model/5.5.4/ome-tiff/specification.html

function get_result(ifd_index, ifd_files, dimlist)
    io = IOBuffer()
    for (idx, ifd) in enumerate(ifd_index)
        (ifd == nothing) && continue
        # get dimensions corresponding to the right image
        # the last index in ifd is position
        dims = dimlist[ifd[end]]
        println(io, idx-1, "\tZ", ifd[loc(dims, :Z)]-1, "-T", ifd[loc(dims, :T)]-1, "-C", ifd[loc(dims, :C)]-1)
    end
    String(take!(io));
end

loc(dims, dimname) = findfirst(isequal(dimname), keys(dims))-2

function get_ifds(fragment)
    omexml = root(parsexml(wrap(fragment)))
    containers = findall("//*[@DimensionOrder]", omexml)
    n_ifds = 0
    dimlist = []
    for container in containers
        dims, _ = OMETIFF.build_axes(container)
        push!(dimlist, dims)
        n_ifds += dims[:Z]*dims[:C]*dims[:T]
    end
    ifd_index = Array{Union{NTuple{4, Int}, Nothing}}(nothing, n_ifds)
    ifd_files = Array{Union{Tuple{String, String}, Nothing}}(nothing, n_ifds)
    obs_filepaths = Set{String}()
    for (idx, container) in enumerate(containers)
        OMETIFF.ifdindex!(ifd_index, ifd_files, obs_filepaths, container, dimlist[idx], "", idx)
    end
    ifd_index, ifd_files, dimlist
end

@testset "Multiple planes per TiffData" begin
    @testset "Empty TiffData" begin
        fragment1 = """
        <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows325"
                Type="uint8" DimensionOrder="XYZTC"
                SizeX="512" SizeY="512" SizeZ="3" SizeT="2" SizeC="2">
            <TiffData/>
        </Pixels>
        """

        expected = """
        0	Z0-T0-C0
        1	Z1-T0-C0
        2	Z2-T0-C0
        3	Z0-T1-C0
        4	Z1-T1-C0
        5	Z2-T1-C0
        6	Z0-T0-C1
        7	Z1-T0-C1
        8	Z2-T0-C1
        9	Z0-T1-C1
        10	Z1-T1-C1
        11	Z2-T1-C1
        """

        @test get_result(get_ifds(fragment1)...) == expected
    end

    @testset "Multiple planes" begin
        fragment2 = """
        <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows462"
                Type="uint8" DimensionOrder="XYCTZ"
                SizeX="512" SizeY="512" SizeZ="4" SizeT="3" SizeC="2">
            <TiffData PlaneCount="10"/>
        </Pixels>
        """

        expected2 = """
        0	Z0-T0-C0
        1	Z0-T0-C1
        2	Z0-T1-C0
        3	Z0-T1-C1
        4	Z0-T2-C0
        5	Z0-T2-C1
        6	Z1-T0-C0
        7	Z1-T0-C1
        8	Z1-T1-C0
        9	Z1-T1-C1
        """

        @test get_result(get_ifds(fragment2)...) == expected2
    end

    @testset "Offset IFD, multiple planes" begin
        fragment3 = """
        <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows197"
                Type="uint8" DimensionOrder="XYZTC"
                SizeX="512" SizeY="512" SizeZ="4" SizeC="3" SizeT="2">
            <TiffData IFD="3" PlaneCount="5"/>
        </Pixels>
        """

        expected3 = """
        3	Z0-T0-C0
        4	Z1-T0-C0
        5	Z2-T0-C0
        6	Z3-T0-C0
        7	Z0-T1-C0
        """

        @test get_result(get_ifds(fragment3)...) == expected3
    end
end

@testset "Reverse temporal order" begin
    fragment4 = """
    <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows789"
            Type="uint8" DimensionOrder="XYZTC"
            SizeX="512" SizeY="512" SizeZ="1" SizeC="1" SizeT="6">
        <TiffData IFD="0" FirstT="5"/>
        <TiffData IFD="1" FirstT="4"/>
        <TiffData IFD="2" FirstT="3"/>
        <TiffData IFD="3" FirstT="2"/>
        <TiffData IFD="4" FirstT="1"/>
        <TiffData IFD="5" FirstT="0"/>
    </Pixels>
    """

    expected4 = """
    0	Z0-T5-C0
    1	Z0-T4-C0
    2	Z0-T3-C0
    3	Z0-T2-C0
    4	Z0-T1-C0
    5	Z0-T0-C0
    """

    @test get_result(get_ifds(fragment4)...) == expected4
end

@testset "Different pixel dimensions" begin
    fragment5 = """
    <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows789"
            Type="uint8" DimensionOrder="XYZTC"
            SizeX="512" SizeY="512" SizeZ="1" SizeC="1" SizeT="6">
        <TiffData IFD="0" FirstT="5"/>
        <TiffData IFD="2" FirstT="4"/>
        <TiffData IFD="4" FirstT="3"/>
        <TiffData IFD="6" FirstT="2"/>
        <TiffData IFD="7" FirstT="1"/>
        <TiffData IFD="8" FirstT="0"/>
    </Pixels>
    <Pixels ID="urn:lsid:loci.wisc.edu:Pixels:ows789"
            Type="uint8" DimensionOrder="XYZTC"
            SizeX="512" SizeY="512" SizeZ="1" SizeC="1" SizeT="3">
        <TiffData IFD="1" FirstT="1"/>
        <TiffData IFD="3" FirstT="2"/>
        <TiffData IFD="5" FirstT="3"/>
    </Pixels>
    """

    expected5 = """
    0	Z0-T5-C0
    1	Z0-T1-C0
    2	Z0-T4-C0
    3	Z0-T2-C0
    4	Z0-T3-C0
    5	Z0-T3-C0
    6	Z0-T2-C0
    7	Z0-T1-C0
    8	Z0-T0-C0
    """
    @test get_result(get_ifds(fragment5)...) == expected5
end

@testset "Missing IFDs" begin

    fragment6 = """
    <Pixels BigEndian="false" DimensionOrder="XYCZT" ID="Pixels:0" PhysicalSizeX="5.2304" PhysicalSizeXUnit="µm" PhysicalSizeY="5.2304" PhysicalSizeYUnit="µm" SizeC="1" SizeT="10" SizeX="256" SizeY="256" SizeZ="1" TimeIncrement="5000.0" TimeIncrementUnit="ms" Type="uint16">
    <Channel ID="Channel:0:0" Name="Default" SamplesPerPixel="1">
        <LightPath/>
    </Channel>
    <TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        <UUID FileName="181003_multi_pos_time_course_1_MMStack.ome.tif">urn:uuid:3872d873-e5a6-421b-a3a2-4d6539f77442</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="1" FirstZ="0" IFD="1" PlaneCount="1">
        <UUID FileName="181003_multi_pos_time_course_1_MMStack.ome.tif">urn:uuid:3872d873-e5a6-421b-a3a2-4d6539f77442</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="2" FirstZ="0" IFD="2" PlaneCount="1">
        <UUID FileName="181003_multi_pos_time_course_1_MMStack.ome.tif">urn:uuid:3872d873-e5a6-421b-a3a2-4d6539f77442</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="3" FirstZ="0" IFD="3" PlaneCount="1">
        <UUID FileName="181003_multi_pos_time_course_1_MMStack.ome.tif">urn:uuid:3872d873-e5a6-421b-a3a2-4d6539f77442</UUID>
    </TiffData>
    </Pixels>
    """

    ifd_index, ifd_files, dimlist = get_ifds(fragment6)

    @test all(ifd_index[1:4] .== [(1, 1, 1, 1), (1, 1, 2, 1), (1, 1, 3, 1), (1, 1, 4, 1)])
    @test all(ifd_index[5:end] .== nothing)
end

@testset "IFD Index Sharing" begin

    fragment7 = """
    <Pixels BigEndian="false" DimensionOrder="XYCTZ" ID="Pixels:0:0" Type="uint8" SizeC="1" SizeT="1" SizeX="18" SizeY="24" SizeZ="5">
    <Channel Color="-2147483648" ID="Channel:0"/>
    <TiffData FirstC="0" FirstT="0" FirstZ="0" IFD="0" PlaneCount="1">
        <UUID FileName="multifile-Z1.ome.tiff">urn:uuid:d25bfea2-2708-4a4f-bcd9-8ab1ac478041</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="0" FirstZ="1" IFD="0" PlaneCount="1">
        <UUID FileName="multifile-Z2.ome.tiff">urn:uuid:7f283a5c-828e-4dec-955f-1e5ed26773b6</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="0" FirstZ="2" IFD="0" PlaneCount="1">
        <UUID FileName="multifile-Z3.ome.tiff">urn:uuid:e3bdd73e-d2a8-4999-9135-b8422265ba18</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="0" FirstZ="3" IFD="0" PlaneCount="1">
        <UUID FileName="multifile-Z4.ome.tiff">urn:uuid:be142aa6-c8d8-40cb-9494-7ef5123308ff</UUID>
    </TiffData>
    <TiffData FirstC="0" FirstT="0" FirstZ="4" IFD="0" PlaneCount="1">
        <UUID FileName="multifile-Z5.ome.tiff">urn:uuid:38731a88-9908-4aec-8df3-18d63c0d24dd</UUID>
    </TiffData>
    </Pixels>
    """

    ifd_index, ifd_files, dimlist = get_ifds(fragment7)

    @test all(ifd_index .== [(1, 1, 1, 1), (1, 1, 2, 1), (1, 1, 3, 1), (1, 1, 4, 1), (1, 1, 5, 1)])
end


text:
@testset "1110.delete-nodes-and-return-forest.jl" begin
    @test del_nodes1110(TreeNode{Int}([1, 2, 3, 4, 5, 6, 7]), [3, 5]) ==
          [TreeNode{Int}([1, 2, nothing, 4]), TreeNode(6), TreeNode(7)]
end


text:
using Plots, LaTeXStrings

Plots.gr()

function example_orbit_propagation_kepler(plot_dir)
    # Simulate Keplerian Orbit
    plot_name = "keplerian_orbit"
    println("Generating: $plot_name.svg")

    epc0 = Epoch(2019, 1, 1, 12, 0, 0, 0.0)
    oe0  = [R_EARTH + 500e3, 0.01, 75.0, 45.0, 30.0, 0.0]
    eci0 = sOSCtoCART(oe0, use_degrees=true)

    # Simulate orbit for one orbit
    T    = orbit_period(oe0[1])
    epcf = epc0 + T
    orb  = EarthInertialState(epc0, eci0, dt=1.0,
            mass=1.0, n_grav=0, m_grav=0,
            drag=false, srp=false,
            moon=false, sun=false,
            relativity=false
    )

    t, epc, eci = sim!(orb, epcf)

    # Create plot with one series
    AX_LIM = R_EARTH + 1000e3
    plt = plot3d(1, xlim=(-AX_LIM,AX_LIM),
                    ylim=(-AX_LIM,AX_LIM),
                    zlim=(-AX_LIM, AX_LIM),
                    title = "Keplerian Orbit",
                    xformatter=:scientific,
                    yformatter=:scientific,
                    zformatter=:scientific,
                    linewidth=2,
                    legend=false)

    for i in 1:length(t)
        push!(plt, eci[1, i], eci[2, i], eci[3, i])
    end
    Plots.savefig(plot_dir * "/$plot_name.svg")

    # Plot Keplerian Elements
    plot_name = "keplerian_elements"
    println("Generating: $plot_name.svg")

    # Compute orbital Elements
    elements = zeros(Float64, 6, length(t))
    plots = []
    for (i, e) in enumerate(epc)
        elements[:, i] = sCARTtoOSC(eci[:, i], use_degrees=true)
    end

    ylabels = [L"a \; [km]", L"e", L"i \; [deg]", L"\Omega \; [deg]", L"\omega \; [deg]", L"M \; [deg]"]
    for i in 1:6
        y = elements[i, 1]
        if i == 6
            push!(plots, plot(t/T, elements[i, :],
                        xlim=(0.0, 1.0),
                        ylim=(0.0, 360.0),
                        ylabel=ylabels[i],
                        legend=false,
                        linewidth=2))
        else
            push!(plots, plot(t/T, elements[i, :],
                        xlim=(0.0, 1.0),
                        ylim=(y-0.1*y, y+0.1*y),
                        ylabel=ylabels[i],
                        legend=false,
                        linewidth=2))
        end
    end

    plot(plots..., layout=(3,2), xlabel=L"T \; [Orbits]")

    Plots.savefig(plot_dir * "/$plot_name.svg")
end

function example_orbit_propagation_fullforce(plot_dir)
    # Simulate Keplerian Orbit
    plot_name = "fullforce_orbit"
    println("Generating: $plot_name.svg")

    epc0 = Epoch(2019, 1, 1, 12, 0, 0, 0.0)
    oe0  = [R_EARTH + 500e3, 0.01, 75.0, 45.0, 30.0, 0.0]
    eci0 = sOSCtoCART(oe0, use_degrees=true)

    # Simulate orbit for one orbit
    T    = orbit_period(oe0[1])
    epcf = epc0 + T
    orb  = EarthInertialState(epc0, eci0, dt=1.0,
            mass=100.0, n_grav=20, m_grav=20,
            drag=true, srp=true,
            moon=true, sun=true,
            relativity=true
    )

    t, epc, eci = sim!(orb, epcf)

    # Create plot with one series
    AX_LIM = R_EARTH + 1000e3
    plt = plot3d(1, xlim=(-AX_LIM,AX_LIM),
                    ylim=(-AX_LIM,AX_LIM),
                    zlim=(-AX_LIM, AX_LIM),
                    title = "Keplerian Orbit",
                    xformatter=:scientific,
                    yformatter=:scientific,
                    zformatter=:scientific,
                    linewidth=2,
                    legend=false)

    for i in 1:length(t)
        push!(plt, eci[1, i], eci[2, i], eci[3, i])
    end
    Plots.savefig(plot_dir * "/$plot_name.svg")

    # Plot Keplerian Elements
    plot_name = "fullforce_elements"
    println("Generating: $plot_name.svg")

    # Compute orbital Elements
    elements = zeros(Float64, 6, length(t))
    plots = []
    for (i, e) in enumerate(epc)
        elements[:, i] = sCARTtoOSC(eci[:, i], use_degrees=true)
    end

    # ylabels = [L"a \; [km]", L"e", L"i \; [deg]", L"\Omega \; [deg]", L"\omega \; [deg]", L"M \; [deg]"]
    ylabels = ["a [km]", "e", "i [deg]", "\\Omega [deg]", "\\omega [deg]", "M [deg]"]
    for i in 1:6
        y_min = minimum(elements[i, :])
        y_max = maximum(elements[i, :])
        if i == 6
            push!(plots, plot(t/T, elements[i, :],
                        xlim=(0.0, 1.0),
                        ylim=(0.0, 360.0),
                        ylabel=ylabels[i],
                        legend=false,
                        linewidth=2))
        else
            push!(plots, plot(t/T, elements[i, :],
                        xlim=(0.0, 1.0),
                        ylim=(y_min-0.01*y_min, y_max+0.01*y_max),
                        ylabel=ylabels[i],
                        legend=false,
                        linewidth=2))
        end
    end

    # plot(plots..., layout=(3,2), xlabel=L"T \; [Orbits]")
    plot(plots..., layout=(3,2), xlabel="T [Orbits]")

    Plots.savefig(plot_dir * "/$plot_name.svg")
end

# Generate all plots
function makeplots()
    # Start building plots
    println("Generating plots")

    # Make output directory
    plot_dir = (pwd()[end-3:end] == "docs") ? "build/plots" : "docs/build/plots"
    if isdir(plot_dir)
        rm(plot_dir, recursive=true)
    end
    mkdir(plot_dir)

    example_orbit_propagation_kepler(plot_dir)
    example_orbit_propagation_fullforce(plot_dir)
end


text:
@testset "bits type" begin
    x = TimeDag.Maybe{Float64}()
    @test !TimeDag.valid(x)
    @test_throws ArgumentError TimeDag.value(x)
    # Checking that this *doesn't* throw, although it could return anything.
    TimeDag.unsafe_value(x)

    y = TimeDag.Maybe(3.0)
    @test TimeDag.valid(y)
    @test TimeDag.value(y) == 3.0
    @test TimeDag.unsafe_value(y) == 3.0
end

@testset "non bits type" begin
    x = TimeDag.Maybe{Vector{Float64}}()
    @test !TimeDag.valid(x)
    @test_throws ArgumentError TimeDag.value(x)
    # This accesses an undefined reference, which is illegal.
    @test_throws UndefRefError TimeDag.unsafe_value(x)

    y = TimeDag.Maybe([3.0])
    @test TimeDag.valid(y)
    @test TimeDag.value(y) == [3.0]
    @test TimeDag.unsafe_value(y) == [3.0]
end