text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(data::Data, vars=nothing)
  # valid variables
  validvars = Tables.schema(values(data)).names

  # plot all variables by default
  isnothing(vars) && (vars = validvars)
  @assert vars ⊆ validvars "invalid variable names"

  # underlying domain where variables live
  dom = domain(data)

  # shared plot specs
  layout --> length(vars)

  for (i, var) in enumerate(vars)
    # data for variable
    vals = data[var]

    # handle categorical values
    v = eltype(vals) <: CategoricalValue ? levelcode.(vals) : copy(vals)

    # handle missing values
    v = replace(v, missing => NaN)

    @series begin
      subplot := i
      title --> string(var)
      legend --> false
      dom, v
    end
  end
end


text:
using Test

randop() = rand(["-->", "→",
                 "||",
                 "&&",
                 "<", "==", "<:", ">:",
                 "<|", "|>",
                 ":",
                 "+", "-",
                 ">>", "<<",
                 "*", "/",
                 "//",
                 "^", "↑",
                 "::",
                 ".", "->"])

test_expr_broken(str) = test_expr(str, false)

function test_expr(str, show_data=true)
    x, ps = CSTParser.parse(ParseState(str))

    x0 = Expr(x)
    x1 = remlineinfo!(Meta.parse(str))
    if CSTParser.has_error(ps) || x0 != x1
        if show_data
            println("Mismatch between flisp and CSTParser when parsing string $str")
            println("ParserState:\n $ps\n")
            println("CSTParser Expr:\n $x\n")
            println("Converted CSTParser Expr:\n $x0\n")
            println("Base EXPR:\n $x1\n")
        end
        return false
    end
    return true
end

@testset "All tests" begin
    @test Meta.parse("(1,)") == Expr(:tuple, 1)
    @testset "Operators" begin
    # @testset "Binary Operators" begin
    #     for iter = 1:25
    #         println(iter)
    #         str = join([["x$(randop())" for i = 1:19];"x"])

    #         @test test_expr(str)
    #     end
    # end
        @testset "Conditional Operator" begin
            @test test_expr("a ? b : c")
            @test test_expr("a ? b : c : d")
            @test test_expr("a ? b : c : d : e")
            @test test_expr("a ? b : c : d : e")
        end

        @testset "Dot Operator" begin
            @test "a.b"  |> test_expr
            @test "a.b.c"  |> test_expr
            @test "(a(b)).c"  |> test_expr
            @test "(a).(b).(c)"  |> test_expr
            @test "(a).b.(c)"  |> test_expr
            @test "(a).b.(c+d)"  |> test_expr
        end

        @testset "Unary Operator" begin
            @test "+" |> test_expr
            @test "-" |> test_expr
            @test "!" |> test_expr
            @test "~" |> test_expr
            @test "&" |> test_expr
        # @test "::" |> test_expr
            @test "<:" |> test_expr
            @test ">:" |> test_expr
            @test "¬" |> test_expr
            @test "√" |> test_expr
            @test "∛" |> test_expr
            @test "∜" |> test_expr
        end

        @testset "Unary Operator" begin
            @test "a=b..." |> test_expr
            @test "a-->b..." |> test_expr
            @test "a&&b..." |> test_expr
            @test "a||b..." |> test_expr
            @test "a<b..." |> test_expr
            @test "a:b..." |> test_expr
            @test "a+b..." |> test_expr
            @test "a<<b..." |> test_expr
            @test "a*b..." |> test_expr
            @test "a//b..." |> test_expr
            @test "a^b..." |> test_expr
            @test "a::b..." |> test_expr
            @test "a where b..." |> test_expr
            @test "a.b..." |> test_expr
        end

        @testset "unary op calls" begin
            @test "+(a,b)" |> test_expr
            @test "-(a,b)" |> test_expr
            @test "!(a,b)" |> test_expr
            @test "¬(a,b)" |> test_expr
            @test "~(a,b)" |> test_expr
            @test "<:(a,b)" |> test_expr
            @test "√(a,b)" |> test_expr
            @test "\$(a,b)" |> test_expr
            @test ":(a,b)" |> test_expr
            @test "&a" |> test_expr
            @test "&(a,b)" |> test_expr
            @test "::a" |> test_expr
            @test "::(a,b)" |> test_expr
        end

        @testset "where precedence" begin
            @test "a = b where c = d" |> test_expr
            @test "a = b where c" |> test_expr
            @test "b where c = d" |> test_expr

            @test "a ? b where c : d" |> test_expr

            @test "a --> b where c --> d" |> test_expr
            @test "a --> b where c" |> test_expr
            @test "b where c --> d" |> test_expr

            @test "a || b where c || d" |> test_expr
            @test "a || b where c" |> test_expr
            @test "b where c || d" |> test_expr

            @test "a && b where c && d" |> test_expr
            @test "a && b where c" |> test_expr
            @test "b where c && d" |> test_expr

            @test "a <: b where c <: d" |> test_expr
            @test "a <: b where c" |> test_expr
            @test "b where c <: d" |> test_expr

            @test "a <| b where c <| d" |> test_expr
            @test "a <| b where c" |> test_expr
            @test "b where c <| d" |> test_expr

            @test "a : b where c : d" |> test_expr
            @test "a : b where c" |> test_expr
            @test "b where c : d" |> test_expr

            @test "a + b where c + d" |> test_expr
            @test "a + b where c" |> test_expr
            @test "b where c + d" |> test_expr

            @test "a << b where c << d" |> test_expr
            @test "a << b where c" |> test_expr
            @test "b where c << d" |> test_expr

            @test "a * b where c * d" |> test_expr
            @test "a * b where c" |> test_expr
            @test "b where c * d" |> test_expr

            @test "a // b where c // d" |> test_expr
            @test "a // b where c" |> test_expr
            @test "b where c // d" |> test_expr

            @test "a ^ b where c ^ d" |> test_expr
            @test "a ^ b where c" |> test_expr
            @test "b where c ^ d" |> test_expr

            @test "a :: b where c :: d" |> test_expr
            @test "a :: b where c" |> test_expr
            @test "b where c :: d" |> test_expr

            @test "a.b where c.d" |> test_expr
            @test "a.b where c" |> test_expr
            @test "b where c.d" |> test_expr

            @test "a where b where c" |> test_expr
        end
    end

    @testset "Type Annotations" begin
        @testset "Curly" begin
            @test "x{T}" |> test_expr
            @test "x{T,S}" |> test_expr
            @test "a.b{T}" |> test_expr
            @test "a(b){T}" |> test_expr
            @test "(a(b)){T}" |> test_expr
            @test "a{b}{T}" |> test_expr
            @test "a{b}(c){T}" |> test_expr
            @test "a{b}.c{T}" |> test_expr
            @test """x{T,
        S}""" |> test_expr
        end
    end

    @testset "Tuples" begin
        @static if VERSION > v"1.1-"
            @test typof(CSTParser.parse("1,")) === CSTParser.ErrorToken
        else
            @test "1," |> test_expr
        end
        @test "1,2" |> test_expr
        @test "1,2,3" |> test_expr
        @test "()" |> test_expr
        @test "(==)" |> test_expr
        @test "(1)" |> test_expr
        @test "(1,)" |> test_expr
        @test "(1,2)" |> test_expr
        @test "(a,b,c)" |> test_expr
        @test "(a...)" |> test_expr
        @test "((a,b)...)" |> test_expr
        @test "a,b = c,d" |> test_expr
        @test "(a,b) = (c,d)" |> test_expr
    end

    @testset "Function Calls" begin
        @testset "Simple Calls" begin
            @test "f(x)" |> test_expr
            @test "f(x,y)" |> test_expr
            @test "f(g(x))" |> test_expr
            @test "f((x,y))" |> test_expr
            @test "f((x,y), z)" |> test_expr
            @test "f(z, (x,y), z)" |> test_expr
            @test "f{a}(x)" |> test_expr
            @test "f{a<:T}(x::T)" |> test_expr
        end

        @testset "Keyword Arguments" begin
            @test "f(x=1)" |> test_expr
            @test "f(x=1,y::Int = 1)" |> test_expr
        end

        @testset "Compact Declaration" begin
            @test "f(x) = x" |> test_expr
            @test "f(x) = g(x)" |> test_expr
            @test "f(x) = (x)" |> test_expr
            @test "f(x) = (x;y)" |> test_expr
            @test "f(g(x)) = x" |> test_expr
            @test "f(g(x)) = h(x)" |> test_expr
        end

        @testset "Standard Declaration" begin
            @test "function f end" |> test_expr
            @test "function f(x) x end" |> test_expr
            @test "function f(x); x; end" |> test_expr
            @test "function f(x) x; end" |> test_expr
            @test "function f(x); x end" |> test_expr
            @test "function f(x) x;y end" |> test_expr
            @test """function f(x) x end""" |> test_expr
            @test """function f(x,y =1) x end""" |> test_expr
            @test """function f(x,y =1;z =2) x end""" |> test_expr
        end
        @testset "Anonymous" begin
            @test "x->y" |> test_expr
            @test "(x,y)->x*y" |> test_expr
            @test """function ()
            return
        end""" |> test_expr
        end
    end

    @testset "Types" begin
        @testset "Abstract" begin
            @test "abstract type t end" |> test_expr
            @test "abstract type t{T} end" |> test_expr
            @test "abstract type t <: S end" |> test_expr
            @test "abstract type t{T} <: S end" |> test_expr
        end

        @testset "primitive" begin
            @test "primitive type Int 64 end" |> test_expr
            @test "primitive type Int 4*16 end" |> test_expr
        end

        @testset "Structs" begin
            @test "struct a end" |> test_expr
            @test "struct a; end" |> test_expr
            @test "struct a; b;end" |> test_expr
            @test """struct a
            arg1
        end""" |> test_expr
            @test """struct a <: T
            arg1::Int
            arg2::Int
        end""" |> test_expr
            @test """struct a
            arg1::T
        end""" |> test_expr
            @test """struct a{T}
            arg1::T
            a(args) = new(args)
        end""" |> test_expr
            @test """struct a <: Int
            arg1::Vector{Int}
        end""" |> test_expr
            @test """mutable struct a <: Int
            arg1::Vector{Int}
        end""" |> test_expr
        end
    end

    @testset "Modules" begin
        @testset "Imports " begin
            @test "import ModA" |> test_expr
            @test "import .ModA" |> test_expr
            @test "import ..ModA.a" |> test_expr
            @test "import ModA.subModA" |> test_expr
            @test "import ModA.subModA: a" |> test_expr
            @test "import ModA.subModA: a, b" |> test_expr
            @test "import ModA.subModA: a, b.c" |> test_expr
            @test "import .ModA.subModA: a, b.c" |> test_expr
            @test "import ..ModA.subModA: a, b.c" |> test_expr
        end
        @testset "Export " begin
            @test "export ModA" |> test_expr
            @test "export a, b, c" |> test_expr
        end
    end

    @testset "Generators" begin
        @test "(y for y in X)" |> test_expr
        @test "((x,y) for x in X, y in Y)" |> test_expr
        @test "(y.x for y in X)" |> test_expr
        @test "((y) for y in X)" |> test_expr
        @test "(y,x for y in X)" |> test_expr
        @test "((y,x) for y in X)" |> test_expr
        @test "[y for y in X]" |> test_expr
        @test "[(y) for y in X]" |> test_expr
        @test "[(y,x) for y in X]" |> test_expr
        @test "Int[y for y in X]" |> test_expr
        @test "Int[(y) for y in X]" |> test_expr
        @test "Int[(y,x) for y in X]" |> test_expr
        @test """
    [a
    for a = 1:2]""" |> test_expr
        @test "[ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]" |> test_expr
        @test "all(d ≥ 0 for d in B.dims)" |> test_expr
        @test "(arg for x in X)" |> test_expr
        @test "(arg for x in X for y in Y)" |> test_expr
        @test "(arg for x in X for y in Y for z in Z)" |> test_expr
        @test "(arg for x in X if A)" |> test_expr
        @test "(arg for x in X if A for y in Y)" |> test_expr
        @test "(arg for x in X if A for y in Y if B)" |> test_expr
        @test "(arg for x in X if A for y in Y for z in Z)" |> test_expr
        @test "(arg for x in X if A for y in Y if B for z in Z)" |> test_expr
        @test "(arg for x in X if A for y in Y if B for z in Z if C)" |> test_expr
        @test "(arg for x in X, y in Y for z in Z)" |> test_expr
        @test "(arg for x in X, y in Y if A for z in Z)" |> test_expr
    end

    @testset "Macros " begin
        @test "macro m end" |> test_expr
        @test "macro m() end" |> test_expr
        @test "macro m() a end" |> test_expr
        @test "@mac" |> test_expr
        @test "@mac a b c" |> test_expr
        @test "@mac f(5)" |> test_expr
        @test "(@mac x)" |> test_expr
        @test "Mod.@mac a b c" |> test_expr
    # @test "[@mac a b]" |> test_expr
        @test "@inline get_chunks_id(i::Integer) = _div64(Int(i)-1)+1, _mod64(Int(i) -1)" |> test_expr
        @test "@inline f() = (), ()" |> test_expr
        @test "@sprintf(\"%08d\", id)" |> test_expr
        @test "[@m @n a for a in A]" |> test_expr # ensure closer.insquare propogates
    end

    @testset "Square " begin
        @testset "vect" begin
            @test "[x]" |> test_expr
            @test "[(1,2)]" |> test_expr
            @test "[x...]" |> test_expr
            @test "[1,2,3,4,5]" |> test_expr
        end

        @testset "ref" begin
            @test "t[i]" |> test_expr
            @test "t[i, j]" |> test_expr
        end

        @testset "vcat" begin
            @test "[x;]" |> test_expr
            @test "[x;y;z]" |> test_expr
            @test """[x
                  y
                  z]""" |> test_expr
            @test """[x
                  y;z]""" |> test_expr
            @test """[x;y
                  z]""" |> test_expr
            @test "[x,y;z]" |> test_expr
        end

        @testset "typed_vcat" begin
            @test "t[x;]" |> test_expr
            @test "t[x;y]" |> test_expr
            @test """t[x
                   y]""" |> test_expr
            @test "t[x;y]" |> test_expr
            @test "t[x y; z]" |> test_expr
            @test "t[x, y; z]" |> test_expr
        end

        @testset "hcat" begin
            @test "[x y]" |> test_expr
        end

        @testset "typed_hcat" begin
            @test "t[x y]" |> test_expr
        end

        @testset "Comprehension" begin
            @test "[i for i = 1:10]" |> test_expr
            @test "Int[i for i = 1:10]" |> test_expr
        end
    end

    @testset "Keyword Blocks" begin
        @testset "If" begin
            @test "if cond end" |> test_expr
            @test "if cond; a; end" |> test_expr
            @test "if cond a; end" |> test_expr
            @test "if cond; a end" |> test_expr
            @test """if cond
            1
            1
        end""" |> test_expr
            @test """if cond
        else
            2
            2
        end""" |> test_expr
            @test """if cond
            1
            1
        else
            2
            2
        end""" |> test_expr
            @test "if 1<2 end" |> test_expr
            @test """if 1<2
            f(1)
            f(2)
        end""" |> test_expr
            @test """if 1<2
            f(1)
        elseif 1<2
            f(2)
        end""" |> test_expr
            @test """if 1<2
            f(1)
        elseif 1<2
            f(2)
        else
            f(3)
        end""" |> test_expr
            @test "if cond a end" |> test_expr
        end

        @testset "Try" begin
        # @test "try f(1) end" |> test_expr
        # @test "try; f(1) end" |> test_expr
        # @test "try; f(1); end" |> test_expr
            @test "try; f(1); catch e; e; end" |> test_expr
            @test "try; f(1); catch e; e end" |> test_expr
            @test "try; f(1); catch e e; end" |> test_expr
            @test """try
            f(1)
        catch
        end""" |> test_expr
            @test """try
            f(1)
        catch
            error(err)
        end""" |> test_expr
            @test """try
            f(1)
        catch err
            error(err)
        end""" |> test_expr
            @test """try
            f(1)
        catch
            error(err)
        finally
            stop(f)
        end""" |> test_expr
            @test """try
            f(1)
        catch err
            error(err)
        finally
            stop(f)
        end""" |> test_expr
            @test """try
            f(1)
        finally
            stop(f)
        end""" |> test_expr
        end
        @testset "For" begin
            @test """for i = 1:10
            f(i)
        end""" |> test_expr
            @test """for i = 1:10, j = 1:20
            f(i)
        end""" |> test_expr
        end

        @testset "Let" begin
            @test """let x = 1
            f(x)
        end""" |> test_expr
            @test """let x = 1, y = 2
            f(x)
        end""" |> test_expr
            @test """let
            x
        end""" |> test_expr
        end

        @testset "Do" begin
            @test """f(X) do x
            return x
        end""" |> test_expr
            @test """f(X,Y) do x,y
            return x,y
        end""" |> test_expr
            @test "f() do x body end" |> test_expr
        end
    end

    @testset "Triple-quoted string" begin
        @test valof(CSTParser.parse("\"\"\" \" \"\"\"")) == " \" "
        @test valof(CSTParser.parse("\"\"\"a\"\"\"")) == "a"
        @test valof(CSTParser.parse("\"\"\"\"\"\"")) == ""
        @test valof(CSTParser.parse("\"\"\"\n\t \ta\n\n\t \tb\"\"\"")) == "a\n\nb"
        @test Expr(CSTParser.parse("\"\"\"\ta\n\tb \$c\n\td\n\"\"\"")) == Expr(:string, "\ta\n\tb ", :c, "\n\td\n")
        @test Expr(CSTParser.parse("\"\"\"\n\ta\n\tb \$c\n\td\n\"\"\"")) == Expr(:string, "\ta\n\tb ", :c, "\n\td\n")
        @test Expr(CSTParser.parse("\"\"\"\n\ta\n\tb \$c\n\td\n\t\"\"\"")) == Expr(:string, "a\nb ", :c, "\nd\n")
        @test Expr(CSTParser.parse("\"\"\"\n\t \ta\$(1+\n1)\n\t \tb\"\"\"")) == Expr(:string, "a", :(1 + 1), "\nb")
        ws = "                         "
        "\"\"\"\n$ws%rv = atomicrmw \$rmw \$lt* %0, \$lt %1 acq_rel\n$(ws)ret \$lt %rv\n$ws\"\"\"" |> test_expr
        ws1 = "        "
        ws2 = "    "
        "\"\"\"\n$(ws1)a\n$(ws1)b\n$(ws2)c\n$(ws2)d\n$(ws2)\"\"\"" |> test_expr
        "\"\"\"\n$(ws1)a\n\n$(ws1)b\n\n$(ws2)c\n\n$(ws2)d\n\n$(ws2)\"\"\"" |> test_expr
        @test "\"\"\"\n$(ws1)α\n$(ws1)β\n$(ws2)γ\n$(ws2)δ\n$(ws2)\"\"\"" |> test_expr
        @test "\"\"\"Float\$(bit)\"\"\"" |> test_expr
        @test kindof(CSTParser.parse("\"\"\"abc\$(de)fg\"\"\"")[3]) == CSTParser.Tokens.STRING
        @test kindof(CSTParser.parse("\"\"\"abc(de)fg\"\"\"")) == CSTParser.Tokens.TRIPLE_STRING
    end

    @testset "No longer broken things" begin
        @test "[ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]" |> test_expr
        @test "all(d ≥ 0 for d in B.dims)" |> test_expr
        @test ":(=)" |> test_expr
        @test ":(1)" |> test_expr
        @test ":(a)" |> test_expr
        @test "(@_inline_meta(); f(x))" |> test_expr
        @test "isa(a,b) != c" |> test_expr
        @test "isa(a,a) != isa(a,a)" |> test_expr
        @test "@mac return x" |> test_expr
        @static if VERSION > v"1.1-"
            @test typof(CSTParser.parse("a,b,").args[4]) === CSTParser.ErrorToken
        else
            @test "a,b," |> test_expr
        end
        @test "m!=m" |> test_expr
        @test "+(x...)" |> test_expr
        @test "+(promote(x,y)...)" |> test_expr
        @test "\$(x...)" |> test_expr #
        @test "ccall(:gethostname, stdcall, Int32, ())" |> test_expr
        @test "@inbounds @ncall a b c" |> test_expr
        @test "(a+b)``" |> test_expr
        @test "(-, ~)" |> test_expr
        @test """function +(x::Bool, y::T)::promote_type(Bool,T) where T<:AbstractFloat
                return ifelse(x, oneunit(y) + y, y)
            end""" |> test_expr
        @test """finalizer(x,x::GClosure->begin
                    ccall((:g_closure_unref,Gtk.GLib.libgobject),Void,(Ptr{GClosure},),x.handle)
                end)""" |> test_expr
        @test "function \$A end" |> test_expr
        @test "&ctx->exe_ctx_ref" |> test_expr
        @test ":(\$(docstr).\$(TEMP_SYM)[\$(key)])" |> test_expr
        @test "SpecialFunctions.\$(fsym)(n::Dual)" |> test_expr
        @test "(Base.@_pure_meta;)" |> test_expr
        @test "@M a b->(@N c = @O d e f->g)" |> test_expr
        @test "! = f" |> test_expr
        @test "[a=>1, b=>2]" |> test_expr
        @test "a.\$(b)" |> test_expr
        @test "a.\$f()" |> test_expr
        @test "4x/y" |> test_expr
        @test """
    ccall(:jl_finalize_th, Void, (Ptr{Void}, Any,),
                Core.getptls(), o)
    """ |> test_expr
        @test """
    A[if n == d
        i
    else
        (indices(A,n) for n = 1:nd)
    end...]
    """ |> test_expr
        @test """
    @spawnat(p,
        let m = a
            isa(m, Exception) ? m : nothing
        end)
    """ |> test_expr #
        @test "[@spawn f(R, first(c), last(c)) for c in splitrange(length(R), nworkers())]" |> test_expr
        @test "M.:(a)" |> test_expr
        @test """
            begin
                for i in I for j in J
                    if cond
                        a
                    end
                end end
            end""" |> test_expr
        @test "-f.(a.b + c)" |> test_expr
        @test ":(import Base: @doc)" |> test_expr
        @test "[a for a in A for b in B]" |> test_expr
        @test "+(a,b,c...)" |> test_expr
        @test """@testset a for t in T
        t
    end""" |> test_expr
        @test "import Base.==" |> test_expr
        @test "a`text`" |> test_expr
        @test "a``" |> test_expr
        @test "a`text`b" |> test_expr
        @test "[a; a 0]" |> test_expr
        @test "[a, b; c]" |> test_expr
        @test "t{a; b} " |> test_expr
        @test "a ~ b + c -d" |> test_expr
        @test "y[j=1:10,k=3:2:9; isodd(j+k) && k <= 8]" |> test_expr
        @test "(8=>32.0, 12=>33.1, 6=>18.2)" |> test_expr
        @test "(a,b = c,d)" |> test_expr
        @test "[ -1 -2;]" |> test_expr
        @test "-2y" |> test_expr # precedence
        @test "'''" |> test_expr # tokenize
        @test """
        if j+k <= deg +1
        end
        """ |> test_expr
        @test "function f() ::T end" |> test_expr # ws closer
        @test "import Base: +, -, .+, .-" |> test_expr
        @test "[a +   + l]" |> test_expr # ws closer
        @test "@inbounds C[i,j] = - α[i] * αjc" |> test_expr
        @test "@inbounds C[i,j] = - n * p[i] * pj" |> test_expr
        @test """
        if ! a
            b
        end
        """ |> test_expr # ws closer
        @test "[:-\n:+]" |> test_expr
        @test "::a::b" |> test_expr
        @test "-[1:nc]" |> test_expr
        @test "f() where {a} = b" |> test_expr
        @test "@assert .!(isna(res[2]))" |> test_expr # v0.6
        @test "-((attr.rise / PANGO_SCALE)pt).value" |> test_expr
        @test "!(a = b)" |> test_expr
        @test "-(1)a" |> test_expr
        @test "!(a)::T" |> test_expr
        @test "a::b where T<:S" |> test_expr
        @test "+(x::Bool, y::T)::promote_type(Bool,T) where T<:AbstractFloat" |> test_expr
        @test "T where V<:(T where T)" |> test_expr
        @test "function ^(z::Complex{T}, p::Complex{T})::Complex{T} where T<:AbstractFloat end" |> test_expr
        @test "function +(a) where T where S end" |> test_expr
        @test "function -(x::Rational{T}) where T<:Signed end" |> test_expr
        @test "\$(a)(b)" |> test_expr
        @test "if !(a) break end" |> test_expr
        @test "module a() end" |> test_expr
        @test "M.r\"str\" " |> test_expr
        @test "f(a for a in A if cond)" |> test_expr
        @test "\"dimension \$d is not 1 ≤ \$d ≤ \$nd\" " |> test_expr
        @test "-(-x)^1" |> test_expr
        @test """
        "\\\\\$ch"
        """ |> test_expr
        @test "µs" |> test_expr # normalize unicode
        @test """(x, o; p = 1) -> begin
    return o, p
    end""" |> test_expr # normalize unicode
        @test """(x, o...; p...) -> begin
    return o, p
    end""" |> test_expr # normalize unicode
        @test "function func() where {A where T} x + 1 end" |> test_expr # nested where
        @test "(;x)" |> test_expr # issue 39
        @test """let f = ((; a = 1, b = 2) -> ()),
    m = first(methods(f))
    @test DSE.keywords(f, m) == [:a, :b]
end""" |> test_expr
        @test "-1^a" |> test_expr
        @test "function(f, args...; kw...) end" |> test_expr
        @test "2a * b" |> test_expr
        @test "(g1090(x::T)::T) where {T} = x+1.0" |> test_expr
        @test "(:) = Colon()" |> test_expr
        @test "a + in[1]" |> test_expr
        @test "function f(ex) +a end" |> test_expr
        @test "x`\\\\`" |> test_expr
        @test "x\"\\\\\"" |> test_expr
        @test "x\"\\\\ \"" |> test_expr
        @test "a.{1}" |> test_expr
        @test "@~" |> test_expr
        @test "\$\$(x)" |> test_expr
        @test "\$\$(x)" |> test_expr
        @test CSTParser.typof(CSTParser.parse("=")) === CSTParser.ErrorToken
        @test CSTParser.typof(CSTParser.parse("~")) === CSTParser.OPERATOR
        @test "(1:\n2)" |> test_expr
        @test "a[: ]" |> test_expr
    end

    @testset "interpolation error catching" begin
        x = CSTParser.parse("\"a \$ b\"")
        @test x.fullspan == 7
        @test CSTParser.typof(x[2]) === CSTParser.ErrorToken
        x = CSTParser.parse("\"a \$# b\"")
        @test x.fullspan == 8
        @test CSTParser.typof(x[2]) === CSTParser.ErrorToken
    end

    @testset "Broken things" begin
        @test_broken "\$(a) * -\$(b)" |> test_expr_broken
    end

# test_fsig_decl(str) = (x->x.id).(CSTParser._get_fsig(CSTParser.parse(str)).defs)
# @testset "func-sig variable declarations" begin
#     @test test_fsig_decl("f(x) = x") == [:x]
#     @test test_fsig_decl("""function f(x)
#         x
#     end""") == [:x]

#     @test test_fsig_decl("f{T}(x::T) = x") == [:T, :x]
#     @test test_fsig_decl("""function f{T}(x::T)
#         x
#     end""") == [:T, :x]

#     @test test_fsig_decl("f(x::T) where T = x") == [:T, :x]
#     @test test_fsig_decl("""function f(x::T) where T
#         x
#     end""") == [:T, :x]

#     @test test_fsig_decl("f(x::T{S}) where T where S = x") == [:T, :S, :x]
#     @test test_fsig_decl("""function f(x::T{S}) where T where S
#         x
#     end""") == [:T, :S, :x]
# end

    @testset "Spans" begin
        CSTParser.parse(raw"""
    "ABC$(T)"
    """).fullspan >= 9
        CSTParser.parse("\"_\"").fullspan == 3
        CSTParser.parse("T.mutable && print(\"Ok\")").fullspan == 24
        CSTParser.parse("(\"\$T\")").fullspan == 6
        CSTParser.parse("\"\"\"\$T is not supported\"\"\"").fullspan == 25
        CSTParser.parse("using Compat: @compat\n").fullspan == 22
        CSTParser.parse("primitive = 1").fullspan == 13
    end

    @testset "Command or string with unicode" begin
        @test "```αhelloworldω```" |> test_expr
        @test "\"αhelloworldω\"" |> test_expr
    end

    @testset "conversion of floats with underscore" begin
        @test "30.424_876_125_859_513" |> test_expr
    end

    @testset "errors" begin
        @test typof(CSTParser.parse("1? b : c ")[1]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("1 ?b : c ")[2]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("1 ? b :c ")[4]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("1:\n2")[2]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("1.a")[1]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("f ()")) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("f{t} ()")) === CSTParser.ErrorToken
        @test typof(CSTParser.parse(": a")[1]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("const a")[2]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("const a = 1")[2]) === CSTParser.BinaryOpCall
        @test typof(CSTParser.parse("const global a")[2]) === CSTParser.ErrorToken
        @test typof(CSTParser.parse("const global a = 1")[2]) === CSTParser.Global
    end

    @testset "tuple params" begin
        @test "1,2,3" |> test_expr
        @test "1;2,3" |> test_expr
        @test "1,2;3" |> test_expr
        @test "(1,2,3)" |> test_expr
        @test "(1;2,3)" |> test_expr
        @test "(1,2;3)" |> test_expr
        @test "f(;)" |> test_expr
    end

    @testset "docs" begin
        @test "\"doc\"\nT" |> test_expr
        @test "@doc \"doc\" T" |> test_expr
        @test "@doc \"doc\"\nT" |> test_expr
        @test "@doc \"doc\n\n\n\"\nT" |> test_expr
        @test "begin\n@doc \"doc\"\n\nT\nend" |> test_expr
        @test "@doc \"I am a module\" ModuleMacroDoc" |> test_expr
    end

    @testset "braces" begin
        @test "{a}" |> test_expr
        @test "{a, b}" |> test_expr
        @test "{a, b; c}" |> test_expr
        @test "{a, b; c = 1}" |> test_expr
        @test "{a b}" |> test_expr
        @test "{a b; c}" |> test_expr
        @test "{a b; c = 1}" |> test_expr
    end

    @testset "import preceding dot whitespace" begin
        @test CSTParser.parse("using . M")[2].fullspan == 2
        @test CSTParser.parse("using .. M")[3].fullspan == 2
        @test CSTParser.parse("using ... M")[4].fullspan == 2
    end

    @testset "issue #116" begin
        @test """
    function foo() where {A <:B}
        body
    end""" |> test_expr

        @test """
    function foo() where {A <: B}
        body
    end""" |> test_expr
    end

    @testset "issue #165" begin
        x = CSTParser.parse("""
    a ? b
    function f end""")
        @test length(x) == 5 # make sure we always give out an EXPR of the right length
        @test typof(x[4]) === CSTParser.ErrorToken
        @test typof(x[4][1]) === CSTParser.OPERATOR
    end

    @testset "issue #182" begin
        x = CSTParser.parse("""
        quote
            \"\"\"
            txt
            \"\"\"
            sym
        end""")
        @test typof(x[2][1][1]) === CSTParser.GlobalRefDoc
    end
    if VERSION > v"1.3.0-"
        @testset "issue #198" begin
            @test test_expr(":var\"id\"")
        end
    end
    @testset "vscode issue #1632" begin
        @test test_expr("\"\$( a)\"")
        @test test_expr("\"\$(#=comment=# a)\"")
    end
end

ext: jl


text:
mutable struct BruteForceRiskRegion
    μ::Vector{Float64}
    Σ::Matrix{Float64}
    lattice::Array{Vector{Float64}}
    α::Float64
    function BruteForceRiskRegion(μ::Vector{Float64}, Σ::Matrix{Float64}, K::FiniteCone, α::Float64, lattice_width::Int64 = 100)
        checkRiskRegionArgs(μ, Σ, K, α)
        cone_coords = simplex(lattice_width, K.num_gen)
        num_points = length(cone_coords)
        lattice = Array{Vector{Float64}}(undef, num_points)
        for (i,p) in enumerate(cone_coords)
            lattice[i] = -K.A*p
        end
        return new(μ, Σ, lattice, α)
    end
end

function BruteForceRiskRegion(dist::MvNormal, K::FiniteCone, β::Float64, lattice_width::Int64 = 100)
    BruteForceRiskRegion(dist.μ, dist.Σ.mat, K, quantile(Normal(), β), lattice_width)
end

function BruteForceRiskRegion(dist::MvTDist, K::FiniteCone, β::Float64, lattice_width::Int64 = 100)
    BruteForceRiskRegion(dist.μ, dist.Σ.mat, K, quantile(TDist(dist.df), β), lattice_width)
end

function ∈(y::Vector{Float64}, Ω::BruteForceRiskRegion)
    for x in Ω.lattice
        if dot(x,y-Ω.μ) >= sqrt(dot(x,Ω.Σ*x)) * Ω.α
            return true
        end
    end
    false
end

##
# Returns a matrix whose columns are points of a simplex
# Arguments:
# N - the sum of the elements of each point in the simplex
# n -- the dimension of the space containing the simplex
##
function simplex(N::Int64, n::Int64)
    simplices = Array{Vector{Int64}}(undef, 0)
    if n > 1
        for i in 0:N
            for p in simplex(N-i, n-1)
                push!(p, i)
                push!(simplices, p)
            end
        end
        return simplices
    end
    if n == 1
        points = Array{Vector{Int64}}(undef, 1)
        points[1] = [N]
        return points
    end
end

dim(Ω::BruteForceRiskRegion) = length(Ω.μ)


text:
########################################
## File Name: sac_controller.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/18
## Description: SAC Controller for SACBP
########################################

using Convex
# using SCS
using ECOS
import Future
using Distributed

function getControlCoeffs(simModel::SimulatePosRangeLocalization2D,
                          s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                          UArray::Vector{<:PosControl},
                          rng::MersenneTwister)
    YArray,SArray,SArray_before_update = simulateForward(simModel,s_init,UArray,rng);
    CoSArray = simulateBackward(simModel,SArray,SArray_before_update,UArray,YArray);
    CoSParams = map(CoS -> vec(CoS.p), CoSArray);
    HArray = map(s -> trans_u_coeff(simModel.continuousTransModel,s.p), SArray);
    CoeffArray = similar(CoSParams);
    for ii = 1:length(CoeffArray)
        CoeffArray[ii] = HArray[ii]'*CoSParams[ii];
    end
    #if myid() == 5
    #    println(vec(CoSArray[1].p)[1])
    #end
    return CoeffArray;
end
function getControlCoeffs(simModel::SimulatePosRangeLocalization2D,
                          nominalPolicy::GradientMultiTargetLocalizationPolicy,
                          s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                          u_param_min::Vector{<:Real},
                          u_param_max::Vector{<:Real},
                          UArray::Vector{<:PosControl},
                          rng::MersenneTwister)
    YArray,SArray,SArray_before_update,UArrayClosed = simulateForward(simModel,nominalPolicy,s_init,u_param_min,u_param_max,UArray,rng);
    CoSArray = simulateBackward(simModel,SArray,SArray_before_update,UArrayClosed,YArray);
    CoSParams = map(CoS -> vec(CoS.p), CoSArray);
    HArray = map(s -> trans_u_coeff(simModel.continuousTransModel,s.p), SArray);
    CoeffArray = similar(CoSParams);
    NominalControlCostArray = Vector{Float64}(undef,length(CoeffArray)-1);
    for ii = 1:length(CoeffArray)
        CoeffArray[ii] = HArray[ii]'*CoSParams[ii];
        if ii != length(CoeffArray)
            NominalControlCostArray[ii] = dot(HArray[ii]'*CoSParams[ii],vec(UArrayClosed[ii]))+1/2*vec(UArrayClosed[ii])'*simModel.Cu*vec(UArrayClosed[ii]);
        end
    end
    #if myid() == 5
    #    println(vec(CoSArray[1].p)[1])
    #end
    return CoeffArray,NominalControlCostArray
end
function controlCoeffsExpected(simModel::SimulatePosRangeLocalization2D,  # This is slower than pmap() implementation.
                               s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                               UArray::Vector{<:PosControl},
                               numSamples::Int64,
                               rng::MersenneTwister)
    rngVec = [rng; accumulate(Future.randjump, fill(big(10)^20, nprocs()-1), init=rng)] # Use randjump to get around unresolved issue with parallel RNG.
    cosMat = @distributed (+) for ii = 1:numSamples hcat(getControlCoeffs(simModel,s_init,UArray,rngVec[myid()])...) end
    return cosMat./numSamples
end

function controlCoeffsExpected(simModel::SimulatePosRangeLocalization2D,
                               nominalPolicy::GradientMultiTargetLocalizationPolicy,
                               s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                               u_param_min::Vector{<:Real},
                               u_param_max::Vector{<:Real},
                               UArray::Vector{<:PosControl},
                               numSamples::Int64,
                               rng::MersenneTwister)
    rngVec = [rng; accumulate(Future.randjump, fill(big(10)^20, nprocs()-1), init=rng)] # Use randjump to get around unresolved issue with parallel RNG.
    hcat_coeffs(coeffs) = [hcat(coeffs[1]...),hcat(coeffs[2]...)];
    cosMat = @distributed (+) for ii = 1:numSamples hcat_coeffs(getControlCoeffs(simModel,nominalPolicy,s_init,u_param_min,u_param_max,UArray,rngVec[myid()])) end
    return cosMat[1]./numSamples,[cosMat[2]./numSamples...]
end

function optControlSchedule(simModel::SimulatePosRangeLocalization2D,
                            UArray::Vector{<:PosControl{T}},
                            coeffMat::Matrix{T}) where T <: Real
    UOptArray = similar(UArray);
    Cu_inv = inv(simModel.Cu);
    CostArray = Vector{T}(undef,length(UArray));
    for t = 1:length(UOptArray)
        coeffVec = coeffMat[:,t];
        u = -Cu_inv*coeffVec;
        CostArray[t] = -1/2*coeffVec'*Cu_inv*coeffVec;
        UOptArray[t] = PosControl(UArray[t].t,u);
    end
    return UOptArray,CostArray
end
function optControlSchedule(simModel::SimulatePosRangeLocalization2D,
                            UArray::Vector{<:PosControl{T}},
                            coeffMat::Matrix{T},
                            u_params_min::Vector{T},
                            u_params_max::Vector{T}) where T <: Real
    UOptArray = similar(UArray);
    Cu = simModel.Cu;
    Cu_inv = inv(Cu);
    CostArray = Vector{T}(undef,length(UArray));
    if isdiag(Cu)
        for t = 1:length(UOptArray)
            coeffVec = coeffMat[:,t];
            u_unconstrained = -Cu_inv*coeffVec;
            u_constrained = similar(u_unconstrained);
            for ii = 1:length(u_constrained)
                if u_unconstrained[ii] < u_params_min[ii]
                    u_constrained[ii] = u_params_min[ii];
                elseif u_unconstrained[ii] > u_params_max[ii]
                    u_constrained[ii] = u_params_max[ii];
                else
                    u_constrained[ii] = u_unconstrained[ii];
                end
            end
            UOptArray[t] = PosControl(UArray[t].t,vec(u_constrained));
            CostArray[t] = 1/2*u_constrained'*Cu*u_constrained + dot(u_constrained,coeffVec);
        end
        return UOptArray,CostArray
    else
        warn("Cost coefficient matrix not diagonal. Convex.jl is used for optimizing the control.")
        u = Variable(2);
        constraint = [u <= u_params_max, u >= u_params_min];
        for t = 1:length(UOptArray)
            coeffVec = coeffMat[:,t];
            cost = 1/2*quadform(u,Cu) + dot(u,coeffVec);
            prob = minimize(cost,constraint);
            solve!(prob,() -> ECOS.Optimizer(verbose=false))#,verbose=false);
            UOptArray[t] = PosControl(UArray[t].t,vec(evaluate(u)));
            CostArray[t] = prob.optval;
        end
        # Convex.clearmemory()
        return UOptArray,CostArray
    end
end

function determineControlTime(tcalc::Float64,
                              simModel::SimulatePosRangeLocalization2D,
                              UOptArray::Vector{<:PosControl},
                              CostArray::Vector{<:Real})
    start = time();
    t_current = UOptArray[1].t;
    t_allowed_min = round(t_current + tcalc*1.05 + simModel.dtexec,digits=5);  # Assuming the total t_calc is 1.05*t_calc so far.
    t_allowed_max = round(t_current + tcalc*1.05 + simModel.dto,digits=5); # This is required so that the control action is fully executed before the next control action takes over.
    index_allowed_min = length(filter(U -> U.t < t_allowed_min, UOptArray)) + 1;
    index_allowed_max = length(filter(U -> U.t <= t_allowed_max, UOptArray));
    index_chosen = findmin(CostArray[index_allowed_min:index_allowed_max])[2];
    tcalc += time() - start;
    return UOptArray[index_chosen+(index_allowed_min-1)],tcalc;
end
function determineControlTime(tcalc::Float64,
                              simModel::SimulatePosRangeLocalization2D,
                              UOptArray::Vector{<:PosControl},
                              CostArray::Vector{<:Real},
                              NominalControlCostArray::Vector{<:Real})
    start = time();
    t_current = UOptArray[1].t;
    t_allowed_min = round(t_current + tcalc*1.05 + simModel.dtexec,digits=5);  # Assuming the total t_calc is 1.05*t_calc so far.
    t_allowed_max = round(t_current + tcalc*1.05 + simModel.dto,digits=5); # This is required so that the control action is fully executed before the next control action takes over.
    index_allowed_min = length(filter(U -> U.t < t_allowed_min, UOptArray)) + 1;
    index_allowed_max = length(filter(U -> U.t <= t_allowed_max, UOptArray));
    index_chosen = findmin(CostArray[index_allowed_min:index_allowed_max]+NominalControlCostArray[index_allowed_min:index_allowed_max])[2];
    tcalc += time() - start;
    return UOptArray[index_chosen+(index_allowed_min-1)],tcalc;
end

function sacControlUpdate(simModel::SimulatePosRangeLocalization2D,
                          s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                          UArray::Vector{<:PosControl},
                          u_params_min::Vector{<:Real},
                          u_params_max::Vector{<:Real},
                          numSamples::Int64, # For forward-backward simulation.
                          rng::MersenneTwister; offline=true)
    tcalc = @elapsed coeffMat = controlCoeffsExpected(simModel,s_init,UArray,numSamples,rng);
    tcalc += @elapsed UOptSchedule,CostArray = optControlSchedule(simModel,UArray,coeffMat,u_params_min,u_params_max);
    if offline
        tcalc = 3/4*simModel.dto; # Assumption for offline simulation.
    end
    UOpt,tcalc = determineControlTime(tcalc,simModel,UOptSchedule,CostArray);
    UArrayNew = copy(UArray);
    ActTimes = Real[];
    for ii = 1:length(UArrayNew)
        if UArrayNew[ii].t >= round(UOpt.t - simModel.dtexec,digits=5) && UArrayNew[ii].t < round(UOpt.t,digits=5)
            push!(ActTimes,UArrayNew[ii].t);
            UArrayNew[ii] = PosControl(UArrayNew[ii].t,vec(UOpt));
        end
    end
    act_time_init,act_time_final = ActTimes[1],ActTimes[end];
    return UArrayNew,act_time_init,act_time_final,tcalc
end
function sacControlUpdate(simModel::SimulatePosRangeLocalization2D,
                          nominalPolicy::GradientMultiTargetLocalizationPolicy,
                          s_init::AugState{<:PhysPos,<:VecBelMvNormal},
                          UArray::Vector{<:PosControl},
                          u_param_min::Vector{<:Real},
                          u_param_max::Vector{<:Real},
                          numSamples::Int64, # For forward-backward simulation.
                          rng::MersenneTwister; offline=true)
    tcalc = @elapsed coeffMat,NominalControlCostArray = controlCoeffsExpected(simModel,nominalPolicy,s_init,u_param_min,u_param_max,UArray,numSamples,rng);
    tcalc += @elapsed UOptSchedule,CostArray = optControlSchedule(simModel,UArray,coeffMat,u_param_min,u_param_max);
    if offline
        tcalc = 3/4*simModel.dto; # Assumption for offline simulation.
    end
    UOpt,tcalc = determineControlTime(tcalc,simModel,UOptSchedule,CostArray,NominalControlCostArray);
    UArrayNew = copy(UArray);
    ActTimes = Real[];
    for ii = 1:length(UArrayNew)
        if UArrayNew[ii].t >= round(UOpt.t - simModel.dtexec,digits=5) && UArrayNew[ii].t < round(UOpt.t,digits=5)
            push!(ActTimes,UArrayNew[ii].t);
            UArrayNew[ii] = PosControl(UArrayNew[ii].t,vec(UOpt));
        end
    end
    act_time_init,act_time_final = ActTimes[1],ActTimes[end];
    return UArrayNew,act_time_init,act_time_final,tcalc
end

function getControlCoeffs(simModel::SimulateManipulate2D,
                          s_init::BelMvNormal,
                          UArray::Vector{<:MControl2D},
                          rng::MersenneTwister)
    YArray,SArray,SArray_before_update = simulateForward(simModel,s_init,UArray,rng);
    CoSArray = simulateBackward(simModel,SArray,SArray_before_update,UArray,YArray);
    CoSParams = map(CoS -> CoS.params, CoSArray);
    HArray = map(s -> trans_u_coeff(simModel.continuousTransModel,s,simModel.Q), SArray);
    CoeffArray = similar(CoSParams);
    for ii = 1:length(CoeffArray)
        CoeffArray[ii] = HArray[ii]'*CoSParams[ii]
    end
    #if myid() == 5
    #    println(CoSArray[1].params[1])
    #end
    return CoeffArray;
end
function getControlCoeffs(simModel::SimulateManipulate2D,
                          nominalPolicy::ManipulatePositionControlPolicy,
                          s_init::BelMvNormal,
                          u_param_min::Vector{<:Real},
                          u_param_max::Vector{<:Real},
                          posGain,
                          rotGain,
                          UArray::Vector{<:MControl2D},
                          rng::MersenneTwister)
    YArray,SArray,SArray_before_update,UArrayClosed = simulateForward(simModel,nominalPolicy,s_init,u_param_min,u_param_max,posGain,rotGain,UArray,rng);
    CoSArray = simulateBackward(simModel,SArray,SArray_before_update,UArrayClosed,YArray);
    CoSParams = map(CoS -> CoS.params, CoSArray);
    HArray = map(s -> trans_u_coeff(simModel.continuousTransModel,s,simModel.Q), SArray);
    CoeffArray = similar(CoSParams);
    NominalControlCostArray = Vector{Float64}(undef,length(CoeffArray)-1);
    for ii = 1:length(CoeffArray)
        CoeffArray[ii] = HArray[ii]'*CoSParams[ii]
        if ii != length(CoeffArray)
            NominalControlCostArray[ii] = dot(HArray[ii]'*CoSParams[ii],vec(UArrayClosed[ii]))+1/2*vec(UArrayClosed[ii])'*simModel.Cu*vec(UArrayClosed[ii]);
        end
    end
    #if myid() == 5
    #    println(CoSArray[1].params[1])
    #end
    return CoeffArray,NominalControlCostArray
end
function controlCoeffsExpected(simModel::SimulateManipulate2D,
                               s_init::BelMvNormal,
                               UArray::Vector{<:MControl2D},
                               numSamples::Int64,
                               rng::MersenneTwister)
    rngVec = [rng; accumulate(Future.randjump, fill(big(10)^20, nprocs()-1), init=rng)] # Use randjump to get around unresolved issue with parallel RNG.
    cosMat = @distributed (+) for ii = 1:numSamples hcat(getControlCoeffs(simModel,s_init,UArray,rngVec[myid()])...) end
    return cosMat./numSamples
end
function controlCoeffsExpected(simModel::SimulateManipulate2D,
                               nominalPolicy::ManipulatePositionControlPolicy,
                               s_init::BelMvNormal,
                               u_param_min::Vector{<:Real},
                               u_param_max::Vector{<:Real},
                               posGain,
                               rotGain,
                               UArray::Vector{<:MControl2D},
                               numSamples::Int64,
                               rng::MersenneTwister)
    rngVec = [rng; accumulate(Future.randjump, fill(big(10)^20, nprocs()-1), init=rng)] # Use randjump to get around unresolved issue with parallel RNG.
    hcat_coeffs(coeffs) = [hcat(coeffs[1]...),hcat(coeffs[2]...)];
    cosMat = @distributed (+) for ii = 1:numSamples hcat_coeffs(getControlCoeffs(simModel,nominalPolicy,s_init,u_param_min,u_param_max,posGain,rotGain,UArray,rngVec[myid()])) end
    return cosMat[1]./numSamples,[cosMat[2]./numSamples...]
end
function optControlSchedule(simModel::SimulateManipulate2D,
                            UArray::Vector{MControl2D{T}},
                            coeffMat::Matrix{T}) where T <: Real
    UOptArray = similar(UArray);
    Cu_inv = inv(simModel.Cu);
    CostArray = Vector{T}(undef,length(UArray));
    for t = 1:length(UOptArray)
        coeffVec = coeffMat[:,t];
        u = -Cu_inv*coeffVec;
        CostArray[t] = -1/2*coeffVec'*Cu_inv*coeffVec;
        UOptArray[t] = MControl2D(UArray[t].t,u)
    end
    return UOptArray,CostArray
end
function optControlSchedule(simModel::SimulateManipulate2D,
                            UArray::Vector{MControl2D{T}},
                            coeffMat::Matrix{T},
                            u_params_min::Vector{T},
                            u_params_max::Vector{T}) where T <: Real
    UOptArray = similar(UArray);
    Cu = simModel.Cu;
    Cu_inv = inv(Cu);
    CostArray = Vector{T}(undef, length(UArray));
    if isdiag(Cu)
        for t = 1:length(UOptArray)
            coeffVec = coeffMat[:,t];
            u_unconstrained = -Cu_inv*coeffVec;
            u_constrained = similar(u_unconstrained);
            for ii = 1:length(u_constrained)
                if u_unconstrained[ii] < u_params_min[ii]
                    u_constrained[ii] = u_params_min[ii];
                elseif u_unconstrained[ii] > u_params_max[ii]
                    u_constrained[ii] = u_params_max[ii];
                else
                    u_constrained[ii] = u_unconstrained[ii];
                end
            end
            UOptArray[t] = MControl2D(UArray[t].t,vec(u_constrained));
            CostArray[t] = 1/2*u_constrained'*Cu*u_constrained + dot(u_constrained,coeffVec);
        end
        return UOptArray,CostArray
    else
        warn("Cost coefficient matrix not diagonal. Convex.jl is used for optimizing the control.")
        u = Variable(3);
        constraint = [u <= u_params_max, u >= u_params_min];
        for t = 1:length(UOptArray)
            coeffVec = coeffMat[:,t];
            cost = 1/2*quadform(u,Cu) + dot(u,coeffVec);
            prob = minimize(cost,constraint);
            solve!(prob,() -> ECOS.Optimizer(verbose=false))#,verbose=false);
            UOptArray[t] = MControl2D(UArray[t].t,vec(evaluate(u)));
            CostArray[t] = prob.optval;
        end
        # Convex.clearmemory()
        return UOptArray,CostArray
    end
end

function determineControlTime(tcalc::Float64,
                              simModel::SimulateManipulate2D,
                              UOptArray::Vector{<:MControl2D},
                              CostArray::Vector{<:Real})
    start = time();
    t_current = UOptArray[1].t;
    t_allowed_min = round(t_current + tcalc*1.05 + simModel.dtexec,digits=5);  # Assuming the total t_calc is 1.05*t_calc so far.
    t_allowed_max = round(t_current + tcalc*1.05 + simModel.dto,digits=5); # This is required so that the control action is fully executed before the next control action takes over.
    index_allowed_min = length(filter(U -> U.t < t_allowed_min, UOptArray)) + 1;
    index_allowed_max = length(filter(U -> U.t <= t_allowed_max, UOptArray));
    index_chosen = findmin(CostArray[index_allowed_min:index_allowed_max])[2];
    tcalc += time() - start;
    return UOptArray[index_chosen+(index_allowed_min-1)],tcalc;
end
function determineControlTime(tcalc::Float64,
                              simModel::SimulateManipulate2D,
                              UOptArray::Vector{<:MControl2D},
                              CostArray::Vector{<:Real},
                              NominalControlCostArray::Vector{<:Real})
    start = time();
    t_current = UOptArray[1].t;
    t_allowed_min = round(t_current + tcalc*1.05 + simModel.dtexec,digits=5);  # Assuming the total t_calc is 1.05*t_calc so far.
    t_allowed_max = round(t_current + tcalc*1.05 + simModel.dto,digits=5); # This is required so that the control action is fully executed before the next control action takes over.
    index_allowed_min = length(filter(U -> U.t < t_allowed_min, UOptArray)) + 1;
    index_allowed_max = length(filter(U -> U.t <= t_allowed_max, UOptArray));
    index_chosen = findmin(CostArray[index_allowed_min:index_allowed_max]+NominalControlCostArray[index_allowed_min:index_allowed_max])[2];
    tcalc += time() - start;
    return UOptArray[index_chosen+(index_allowed_min-1)],tcalc;
end

function sacControlUpdate(simModel::SimulateManipulate2D,
                          s_init::BelMvNormal,
                          UArray::Vector{<:MControl2D},
                          u_params_min::Vector{<:Real},
                          u_params_max::Vector{<:Real},
                          numSamples::Int64, # For forward-backward simulation.
                          rng::MersenneTwister; offline=true)
    tcalc = @elapsed coeffMat = controlCoeffsExpected(simModel,s_init,UArray,numSamples,rng);
    tcalc += @elapsed UOptSchedule,CostArray = optControlSchedule(simModel,UArray,coeffMat,u_params_min,u_params_max);
    if offline
        tcalc = 3/4*simModel.dto; # For offline simulation.
    end
    UOpt,tcalc = determineControlTime(tcalc,simModel,UOptSchedule,CostArray);
    UArrayNew = copy(UArray);
    ActTimes = Real[];
    for ii = 1:length(UArrayNew)
        if UArrayNew[ii].t >= round(UOpt.t - simModel.dtexec,digits=5) && UArrayNew[ii].t < round(UOpt.t,digits=5)
            push!(ActTimes,UArrayNew[ii].t);
            UArrayNew[ii] = MControl2D(UArrayNew[ii].t,vec(UOpt));
        end
    end
    act_time_init,act_time_final = ActTimes[1],ActTimes[end];
    return UArrayNew,act_time_init,act_time_final,tcalc
end
function sacControlUpdate(simModel::SimulateManipulate2D,
                          nominalPolicy::ManipulatePositionControlPolicy,
                          s_init::BelMvNormal,
                          UArray::Vector{<:MControl2D},
                          u_param_min::Vector{<:Real},
                          u_param_max::Vector{<:Real},
                          posGain,
                          rotGain,
                          numSamples::Int64, # For forward-backward simulation.
                          rng::MersenneTwister; offline=true)
    tcalc = @elapsed coeffMat,NominalControlCostArray = controlCoeffsExpected(simModel,nominalPolicy,s_init,u_param_min,u_param_max,posGain,rotGain,UArray,numSamples,rng);
    tcalc += @elapsed UOptSchedule,CostArray = optControlSchedule(simModel,UArray,coeffMat,u_param_min,u_param_max);
    if offline
        tcalc = 3/4*simModel.dto; # For offline simulation.
    end
    UOpt,tcalc = determineControlTime(tcalc,simModel,UOptSchedule,CostArray,NominalControlCostArray);
    UArrayNew = copy(UArray);
    ActTimes = Real[];
    for ii = 1:length(UArrayNew)
        if UArrayNew[ii].t >= round(UOpt.t - simModel.dtexec,digits=5) && UArrayNew[ii].t < round(UOpt.t,digits=5)
            push!(ActTimes,UArrayNew[ii].t);
            UArrayNew[ii] = MControl2D(UArrayNew[ii].t,vec(UOpt));
        end
    end
    act_time_init,act_time_final = ActTimes[1],ActTimes[end];
    return UArrayNew,act_time_init,act_time_final,tcalc
end


text:
using Base: Slice, ScalarIndex

"""
    ContiguousSubArray{T,N,P,I,L}

Like `Base.FastContiguousSubArray` but without requirement for linear
indexing (i.e., type parameter `L` can be `false`).

# Examples

```
julia> A = view(ones(5, 5), :, [1,3]);

julia> A isa Base.FastContiguousSubArray
false

julia> A isa SIMD.ContiguousSubArray
true
```
"""
ContiguousSubArray{T,N,P,
                   I<:Union{Tuple{Union{Slice, AbstractUnitRange}, Vararg{Any}},
                            Tuple{Vararg{ScalarIndex}}},
                   L} = SubArray{T,N,P,I,L}

"""
    ContiguousArray{T,N}

Array types with contiguous first dimension.
"""
ContiguousArray{T,N} = Union{DenseArray{T,N}, ContiguousSubArray{T,N},
    Base.ReinterpretArray{T,N,T2,A} where {A <: Union{DenseArray{T2,N},
                                                      ContiguousSubArray{T2,N}}} where {T2}}

"""
    FastContiguousArray{T,N}

This is the type of arrays that `pointer(A, i)` works.
"""
FastContiguousArray{T,N} = Union{DenseArray{T,N}, Base.FastContiguousSubArray{T,N},
    Base.ReinterpretArray{T,N,T2,A} where {A <: Union{DenseArray{T2,N},
                                                      Base.FastContiguousSubArray{T2,N}}} where {T2}}
# https://github.com/eschnett/SIMD.jl/pull/40#discussion_r254131184
# https://github.com/JuliaArrays/MappedArrays.jl/pull/24#issuecomment-460568978

# vload
@propagate_inbounds function vload(::Type{Vec{N, T}}, ptr::Ptr{T}, mask::Union{Nothing, Vec{N, Bool}}=nothing,
                       ::Val{Aligned}=Val(false), ::Val{Nontemporal}=Val(false)) where {N, T, Aligned, Nontemporal}
    if mask === nothing
        Vec(Intrinsics.load(Intrinsics.LVec{N, T}, ptr, Val(Aligned), Val(Nontemporal)))
    else
        Vec(Intrinsics.maskedload(ptr, mask.data, Val(Aligned), Val(Nontemporal)))
    end
end

@propagate_inbounds function vload(::Type{Vec{N, T}}, a::FastContiguousArray{T,1}, i::Integer, mask=nothing,
                       ::Val{Aligned}=Val(false), ::Val{Nontemporal}=Val(false)) where {N, T, Aligned, Nontemporal}
    @boundscheck checkbounds(a, i:(i+N-1))
    GC.@preserve a begin
        ptr = pointer(a, i)
        vload(Vec{N, T}, ptr, mask, Val(Aligned), Val(Nontemporal))
    end
end
@propagate_inbounds vloada(::Type{Vec{N, T}}, ptr::Ptr{T}, mask=nothing) where {N, T} = vload(Vec{N, T}, ptr, mask, Val(true))
@propagate_inbounds vloadnt(::Type{Vec{N, T}}, ptr::Ptr{T}, mask=nothing) where {N, T} = vload(Vec{N, T}, ptr, mask, Val(true), Val(true))
@propagate_inbounds vloada(::Type{Vec{N, T}}, a::FastContiguousArray{T,1}, i::Integer, mask=nothing) where {N, T} = vload(Vec{N, T}, a, i, mask, Val(true))
@propagate_inbounds vloadnt(::Type{Vec{N, T}}, a::FastContiguousArray{T,1}, i::Integer, mask=nothing) where {N, T} = vload(Vec{N, T}, a, i, mask, Val(true), Val(true))

# vstore
@propagate_inbounds function vstore(x::Vec{N, T}, ptr::Ptr{T}, mask::Union{Nothing, Vec{N, Bool}}=nothing,
                       ::Val{Aligned}=Val(false), ::Val{Nontemporal}=Val(false)) where {N, T, Aligned, Nontemporal}
    if mask === nothing
        Intrinsics.store(x.data, ptr, Val(Aligned), Val(Nontemporal))
    else
        Intrinsics.maskedstore(x.data, ptr, mask.data, Val(Aligned), Val(Nontemporal))
    end
end
@propagate_inbounds function vstore(x::Vec{N, T}, a::FastContiguousArray{T,1}, i::Integer, mask=nothing,
               ::Val{Aligned}=Val(false), ::Val{Nontemporal}=Val(false)) where {N, T, Aligned, Nontemporal}
    @boundscheck checkbounds(a, i:(i+N-1))
    GC.@preserve a begin
        ptr = pointer(a, i)
        vstore(x, ptr, mask, Val(Aligned), Val(Nontemporal))
    end
    return a
end
@propagate_inbounds vstorea(x::Vec, ptr::Ptr, mask=nothing) = vstore(x, ptr, mask, Val(true))
@propagate_inbounds vstorent(x::Vec, ptr::Ptr, mask=nothing) = vstore(x, ptr, mask, Val(true), Val(true))
@propagate_inbounds vstorea(x::Vec, a, i, mask=nothing) = vstore(x, a, i, mask, Val(true))
@propagate_inbounds vstorent(x::Vec, a, i, mask=nothing) = vstore(x, a, i, mask, Val(true), Val(true))

@inline vloadx(ptr::Ptr, mask::Vec{<:Any, Bool}) =
    Vec(Intrinsics.maskedexpandload(ptr, mask.data))

@propagate_inbounds function vloadx(a::FastContiguousArray{T,1},
                                    i::Integer, mask::Vec{N, Bool}) where {N, T}
    @boundscheck checkbounds(a, i:i + N - 1)
    return GC.@preserve a begin
        ptr = pointer(a, i)
        vloadx(ptr, mask)
    end
end

@inline vstorec(x::Vec{N, T}, ptr::Ptr{T}, mask::Vec{N, Bool}) where {N, T} =
    Intrinsics.maskedcompressstore(x.data, ptr, mask.data)

@propagate_inbounds function vstorec(x::Vec{N, T}, a::FastContiguousArray{T,1},
                                     i::Integer, mask::Vec{N, Bool}) where {N, T}
    @boundscheck checkbounds(a, i:i + N - 1)
    GC.@preserve a begin
        ptr = pointer(a, i)
        vstorec(x, ptr, mask)
    end
    return a
end

function valloc(::Type{T}, N::Int, sz::Int) where T
    @assert N > 0
    @assert sz >= 0
    # We use padding to align the address of the first element, and
    # also to ensure that we can access past the last element up to
    # the next full vector width
    padding = N-1 + mod(-sz, N)
    mem = Vector{T}(undef, sz + padding)
    addr = Int(pointer(mem))
    off = mod(-addr, N * sizeof(T))
    @assert mod(off, sizeof(T)) == 0
    off = fld(off, sizeof(T))
    @assert 0 <= off <= padding
    res = view(mem, off+1 : off+sz)
    addr2 = Int(pointer(res))
    @assert mod(addr2, N * sizeof(T)) == 0
    res
end

function valloc(f, ::Type{T}, N::Int, sz::Int) where T
    mem = valloc(T, N, sz)
    @inbounds for i in 1:sz
        mem[i] = f(i)
    end
    mem
end

@inline function _get_vec_pointers(a, idx::Vec{N, Int}) where {N}
    ptrs = pointer(a) + (idx - 1) * sizeof(eltype(a))
end

# Have to be careful with optional arguments and @boundscheck,
# see https://github.com/JuliaLang/julia/issues/30411,
# therefore use @propagate_inbounds
@inline vgather(ptrs::Vec{N,Ptr{T}},
                 mask::Vec{N,Bool}=one(Vec{N,Bool}),
                 ::Val{Aligned}=Val(false)) where {N, T<:ScalarTypes, Aligned} =
    return Vec(Intrinsics.maskedgather(ptrs.data, mask.data))
@inline vgather(ptrs::Vec{<:Any,<:Ptr{<:ScalarTypes}}, ::Nothing, aligned::Val = Val(false)) =
    vgather(ptrs, one(Vec{length(ptrs),Bool}), aligned)
@propagate_inbounds function vgather(a::FastContiguousArray{T,1}, idx::Vec{N, Int},
                                     mask::Vec{N,Bool}=one(Vec{N,Bool}),
                                     ::Val{Aligned}=Val(false)) where {N, T<:ScalarTypes, Aligned}
    @boundscheck for i in 1:N
        checkbounds(a, @inbounds idx[i])
    end
    GC.@preserve a begin
        ptrs = _get_vec_pointers(a, idx)
        return vgather(ptrs, mask, Val(Aligned))
    end
end
@propagate_inbounds vgathera(a, idx, mask) = vgather(a, idx, mask, Val(true))
@propagate_inbounds vgathera(a, idx::Vec{N}) where {N} = vgather(a, idx, one(Vec{N,Bool}), Val(true))

@propagate_inbounds Base.getindex(a::FastContiguousArray{T,1}, idx::Vec{N,Int}) where {N,T} =
    vgather(a, idx)

@propagate_inbounds vscatter(x::Vec{N,T}, ptrs::Vec{N,Ptr{T}},
                             mask::Vec{N,Bool}, ::Val{Aligned}=Val(false)) where {N, T<:ScalarTypes, Aligned} =
    Intrinsics.maskedscatter(x.data, ptrs.data, mask.data)
@inline vscatter(x::Vec{N,T}, ptrs::Vec{N,Ptr{T}},
                 ::Nothing, aligned::Val=Val(false)) where {N, T<:ScalarTypes} =
    vscatter(x, ptrs, one(Vec{N, Bool}), aligned)
@propagate_inbounds function vscatter(x::Vec{N,T}, a::FastContiguousArray{T,1}, idx::Vec{N, Int},
                                      mask::Vec{N,Bool}=one(Vec{N, Bool}),
                                      ::Val{Aligned}=Val(false)) where {N, T<:ScalarTypes, Aligned}
    @boundscheck for i in 1:N
        checkbounds(a, @inbounds idx[i])
    end
    GC.@preserve a begin
        ptrs = _get_vec_pointers(a, idx)
        vscatter(x, ptrs, mask, Val(Aligned))
    end
    return
end
@propagate_inbounds vscattera(x, a, idx, mask) = vscatter(x, a, idx, mask, Val(true))
@propagate_inbounds vscattera(x, a, idx::Vec{N}) where {N}  = vscatter(x, a, idx, one(Vec{N,Bool}), Val(true))

@propagate_inbounds Base.setindex!(a::FastContiguousArray{T,1}, v::Vec{N,T}, idx::Vec{N,Int}) where {N, T} =
    vscatter(v, a, idx)

export VecRange

"""
    VecRange{N}(i::Int)
Analogous to `UnitRange` but for loading SIMD vector of width `N` at
index `i`.
# Examples
```jldoctest
julia> xs = ones(4);
julia> xs[VecRange{4}(1)]  # calls `vload(Vec{4,Float64}, xs, 1)`
<4 x Float64>[1.0, 1.0, 1.0, 1.0]
```
"""
struct VecRange{N}
    i::Int
end

@inline Base.length(idx::VecRange{N}) where {N} = N
@inline Base.first(idx::VecRange) = idx.i
@inline Base.last(idx::VecRange) = idx.i + length(idx) - 1

@inline Base.:+(idx::VecRange{N}, j::Integer) where N = VecRange{N}(idx.i + j)
@inline Base.:+(j::Integer, idx::VecRange{N}) where N = VecRange{N}(idx.i + j)
@inline Base.:-(idx::VecRange{N}, j::Integer) where N = VecRange{N}(idx.i - j)

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, idx::VecRange) =
    (first(inds) <= first(idx)) && (last(idx) <= last(inds))

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, idx::Vec) =
    all(first(inds) <= idx) && all(idx <= last(inds))

@inline _checkarity(::AbstractArray{<:Any,N}, ::Vararg{<:Any,N}) where {N} =
    nothing

@inline _checkarity(::T, ::Any) where {T <: AbstractArray} =
    if IndexStyle(T) isa IndexLinear
        nothing
    else
        throw(ArgumentError("""
        Array type $T does not support indexing with a single index.
        Exactly $(ndims(T)) (non-mask) indices have to be specified.
        """))
    end

_checkarity(::AbstractArray{<:Any,N}, ::Vararg{<:Any,M}) where {N,M} =
    throw(ArgumentError("""
    $M indices are given to $N-dimensional array.
    Exactly $N (non-mask) indices have to be specified when using SIMD.
    """))

# Combined with `_preprocessindices`, helper function `_extractmask`
# extracts `mask` in the tail position.  As slicing tuple is not
# type-stable, we use reverse-of-tail-of-reverse hack to extract
# `mask` at the end of `args`.
@inline _extractmask(mask::Vec{N,Bool}, R::Vararg{Integer}) where N =
    (reverse(R), mask)
@inline _extractmask(R::Vararg{Integer}) = (reverse(R), nothing)
@inline _extractmask(mask::Vec{N,Bool}) where {N} = ((), mask)
@inline _extractmask() = ((), nothing)

@noinline _extractmask(rargs...) =
    throw(ArgumentError("""
    Using SIMD indexing `array[idx, i2, ..., iN, mask]` for `N`-dimensional
    array requires `i2` to `iN` to be all integers and `mask` to be optionally
    a SIMD vector `Vec` of `Bool`s.  Given `(i2, ..., iN, mask)` is
    $(summary(reverse(rargs)))
    """))

_maskedidx(idx, ::Nothing, ::Any) = idx
_maskedidx(idx::Vec, mask::Vec, fst) = vifelse(mask, idx, fst)
_maskedidx(idx::VecRange, mask::Vec, fst) =
    _maskedidx(Vec(ntuple(i -> i - 1 + idx.i, length(mask))), mask, fst)

Base.@propagate_inbounds function _preprocessindices(arr, idx, args)
    I, mask = _extractmask(reverse(args)...)
    _checkarity(arr, idx, I...)
    @boundscheck checkbounds(arr,
                             _maskedidx(idx, mask, first(axes(arr, 1))),
                             I...)
    return I, mask
end

"""
    _pointer(arr, i, I)
Pointer to the element `arr[i, I...]`.
"""
Base.@propagate_inbounds _pointer(arr::Union{Array, Base.ReinterpretArray}, i, I) =
    pointer(arr, LinearIndices(arr)[i, I...])
Base.@propagate_inbounds _pointer(arr::Base.FastContiguousSubArray, i, I) =
    pointer(arr, (i, I...))
Base.@propagate_inbounds _pointer(arr::Base.FastContiguousSubArray, i, I::Tuple{}) =
    pointer(arr, i)
# must be separate methods to resolve ambiguity
Base.@propagate_inbounds _pointer(arr::Base.ReinterpretArray, i, I::Tuple{}) =
    pointer(arr, i)
Base.@propagate_inbounds _pointer(arr::SubArray, i, I) =
    pointer(Base.unsafe_view(arr, 1, I...), i)

Base.@propagate_inbounds function Base.getindex(
        arr::ContiguousArray{T}, idx::VecRange{N},
        args::Vararg{Union{Integer,Vec{N,Bool}}}) where {N,T}
    I, mask = _preprocessindices(arr, idx, args)
    return vload(Vec{N,T}, _pointer(arr, idx.i, I), mask)
end

Base.@propagate_inbounds function Base.setindex!(
        arr::ContiguousArray{T}, v::Vec{N,T}, idx::VecRange{N},
        args::Vararg{Union{Integer,Vec{N,Bool}}}) where {N,T}
    I, mask = _preprocessindices(arr, idx, args)
    vstore(v, _pointer(arr, idx.i, I), mask)
    return arr
end

Base.@propagate_inbounds function Base.getindex(
        arr::ContiguousArray{T}, idx::Vec{N,<:Integer},
        args::Vararg{Union{Integer,Vec{N,Bool}}}) where {N,T}
    I, mask = _preprocessindices(arr, idx, args)
    ptrs = _pointer(arr, 1, I) - sizeof(T) + sizeof(T) * idx
    return vgather(ptrs, mask)
end

Base.@propagate_inbounds function Base.setindex!(
        arr::ContiguousArray{T}, v::Vec{N,T}, idx::Vec{N,<:Integer},
        args::Vararg{Union{Integer,Vec{N,Bool}}}) where {N,T}
    I, mask = _preprocessindices(arr, idx, args)
    ptrs = _pointer(arr, 1, I) - sizeof(T) + sizeof(T) * idx
    vscatter(v, ptrs, mask)
    return arr
end


text:
using Nextion
using Test

@testset "Constants" begin
    @testset "Return" begin
        @testset "Return Codes" begin
            @test Integer(Return.Code.CMD_FINISHED) == 0x01

            _code = Return.code(0x01)
            @test _code == Return.Code.CMD_FINISHED
            @test string(_code) == "CMD_FINISHED"
        end
        @testset "Return Modes" begin
            @test Integer(Return.Mode.NO_RETURN) == 0x00
        end
    end

    @testset "Alignement" begin
        @test Integer(Alignment.Vertical.CENTRE) == 1
        @test Integer(Alignment.Horizontal.CENTRE) == 1
    end

    @testset "Background" begin
        @test Integer(Background.SOLIDCOLOUR) == 1
    end

    @testset "Colour" begin
        @test Integer(Colour.WHITE) == Int32(65535)
    end

    @testset "Scroll" begin
        @testset "Scroll Direction" begin
            @test Integer(Scroll.Direction.LEFT_TO_RIGHT) == 0
        end
    end

    @testset "Baudrate" begin
        @test Baudrate.at(9600) == 9600
        @test Baudrate.id(2400) == 1  # is it necessary? ToDo
        @test_throws ErrorException Baudrate.at(9650) == 9600
    end

    @testset "GPIO" begin
        @testset "GPIO Modes" begin
            @test Integer(GPIO.Mode.PULL_UP_INPUT_MODE) == 0
        end
    end
end


text:
using Statistics: mean
using Oceananigans.Fields: compute!

struct WindowedSpatialAverage{F, S, D}
          field :: F
   field_slicer :: S
           dims :: D
end

"""
    WindowedSpatialAverage(field; dims, field_slicer=FieldSlicer())

Builds a `WindowedSpatialAverage` of `field` that returns the average of `field` over a

slice defined by `field_sicer`.

Example
=======

```jldoctest
julia> using Oceananigans

julia> using Oceananigans.Diagnostics: WindowedSpatialAverage, FieldSlicer

julia> using Oceananigans.OutputWriters: AveragedTimeInterval, NetCDFOutputWriter

julia> grid = RegularRectilinearGrid(size=(4, 6, 4), extent=(1,1,1));

julia> model = IncompressibleModel(grid=grid);

julia> set!(model.velocities.u, 1);

julia> slicer = FieldSlicer(j=3:6, k=1);

julia> U_wsa = WindowedSpatialAverage(model.velocities.u; dims=(1, 2), field_slicer=slicer);

julia> simulation = Simulation(model, Δt=10, stop_iteration=10);

julia> simulation.output_writers[:simple_output] = NetCDFOutputWriter(model, (U_wsa=U_wsa,),
                                                                      schedule = 10,
                                                                      filepath = "windowed_spatial_average_jldoctest.nc");
```
"""
WindowedSpatialAverage(field; dims, field_slicer=FieldSlicer()) = WindowedSpatialAverage(field, field_slicer, dims)

function (wsa::WindowedSpatialAverage)(model)
    compute!(wsa.field)
    window = slice_parent(wsa.field_slicer, wsa.field)
    return dropdims(mean(window, dims=wsa.dims), dims=wsa.dims)
end


text:
export stage, cached_stage, compute, debug_compute, free!, cleanup

###### Scheduler #######

compute(x; options=nothing) = compute(Context(), x; options=options)
compute(ctx, c::Chunk; options=nothing) = c

collect(ctx::Context, t::Thunk; options=nothing) =
    collect(ctx, compute(ctx, t; options=options); options=options)
collect(d::Union{Chunk,Thunk}; options=nothing) =
    collect(Context(), d; options=options)

abstract type Computation end

compute(ctx, c::Computation; options=nothing) =
    compute(ctx, stage(ctx, c); options=options)
collect(c::Computation; options=nothing) =
    collect(Context(), c; options=options)

"""
    compute(ctx::Context, d::Thunk; options=nothing) -> Chunk

Compute a Thunk - creates the DAG, assigns ranks to nodes for tie breaking and
runs the scheduler with the specified options. Returns a Chunk which references
the result.
"""
function compute(ctx::Context, d::Thunk; options=nothing)
    scheduler = get!(PLUGINS, :scheduler) do
        get_type(PLUGIN_CONFIGS[:scheduler])
    end
    res = scheduler.compute_dag(ctx, d; options=options)
    if ctx.log_file !== nothing
        if ctx.log_sink !== LocalEventLog
            logs = get_logs!(ctx.log_sink)
            open(ctx.log_file, "w") do io
                Dagger.show_plan(io, logs, d)
            end
        else
            @warn "Context log_sink not set to LocalEventLog, skipping"
        end
    end
    res
end

function debug_compute(ctx::Context, args...; profile=false, options=nothing)
    @time res = compute(ctx, args...; options=options)
    get_logs!(ctx.log_sink), res
end

function debug_compute(arg; profile=false, options=nothing)
    ctx = Context()
    dbgctx = Context(procs(ctx), LocalEventLog(), profile)
    debug_compute(dbgctx, arg; options=options)
end

Base.@deprecate gather(ctx, x) collect(ctx, x)
Base.@deprecate gather(x) collect(x)

cleanup() = cleanup(Context())
function cleanup(ctx::Context)
    if :scheduler in keys(PLUGINS)
        scheduler = PLUGINS[:scheduler]
        (scheduler).cleanup(ctx)
        delete!(PLUGINS, :scheduler)
    end
    nothing
end

function get_type(s::String)
    local T
    for t in split(s, ".")
        t = Symbol(t)
        if !@isdefined(T)
            T = Base.require(@__MODULE__, t)
        else
            T = Core.eval(T, t)
        end
    end
    T
end

##### Dag utilities #####

"""
    dependents(node::Thunk) -> Dict{Union{Thunk,Chunk}, Set{Thunk}}

Find the set of direct dependents for each task.
"""
function dependents(node::Thunk)
    deps = Dict{Union{Thunk,Chunk}, Set{Thunk}}()
    visited = Set{Thunk}()
    to_visit = Set{Thunk}()
    push!(to_visit, node)
    while !isempty(to_visit)
        next = pop!(to_visit)
        (next in visited) && continue
        if !haskey(deps, next)
            deps[next] = Set{Thunk}()
        end
        for inp in inputs(next)
            if istask(inp) || (inp isa Chunk)
                s = get!(()->Set{Thunk}(), deps, inp)
                push!(s, next)
                if istask(inp) && !(inp in visited)
                    push!(to_visit, inp)
                end
            end
        end
        push!(visited, next)
    end
    deps
end

"""
    noffspring(dpents::Dict{Union{Thunk,Chunk}, Set{Thunk}}) -> Dict{Thunk, Int}

Recursively find the number of tasks dependent on each task in the DAG.
Takes a Dict as returned by [`dependents`](@ref).
"""
function noffspring(dpents::Dict{Union{Thunk,Chunk}, Set{Thunk}})
    noff = Dict{Thunk,Int}()
    to_visit = collect(filter(istask, keys(dpents)))
    while !isempty(to_visit)
        next = popfirst!(to_visit)
        haskey(noff, next) && continue
        off = 0
        has_all = true
        for dep in dpents[next]
            if haskey(noff, dep)
                off += noff[dep] + 1
            else
                pushfirst!(to_visit, next)
                pushfirst!(to_visit, dep)
                has_all = false
                break
            end
        end
        has_all || continue
        noff[next] = off
    end
    noff
end

"""
    order(node::Thunk, ndeps) -> Dict{Thunk,Int}

Given a root node of the DAG, calculates a total order for tie-breaking.

  * Root node gets score 1,
  * rest of the nodes are explored in DFS fashion but chunks
    of each node are explored in order of `noffspring`,
    i.e. total number of tasks depending on the result of the said node.

Args:
- node: root node
- ndeps: result of [`noffspring`](@ref)
"""
function order(node::Thunk, ndeps)
    output = Dict{Thunk,Int}()
    to_visit = Thunk[]
    push!(to_visit, node)
    s = 0
    while !isempty(to_visit)
        next = popfirst!(to_visit)
        haskey(output, next) && continue
        s += 1
        output[next] = s
        parents = collect(Iterators.filter(istask, inputs(next)))
        sort!(parents, by=k->get(ndeps,k,0))
        append!(to_visit, parents)
    end
    output
end


text:
# Guide
# Pkg.dev("SixDOF")
# Pkg.add("SixDOF")
# Pkg.test("SixDOF")
# include("SixDOF.jl")

using Plots
```

First, we import the module.

```
using SixDOF
```

There are several structs we need to define.  The inputs used in this example correspond to the Zagi flying wing in Appendix E of Small Unmanned Aircraft: Theory and Practice by Beard and McLain.  We specify the mass properties:

```
```

For this example:
```
m = 1.56
Ixx = 0.1147
Iyy = 0.0576
Izz = 0.1712
Ixz = 0.0015
mp = MassProp(m, Ixx, Iyy, Izz, Ixz)
nothing # hide
```

Next, we specify reference areas and lengths:

```
Sref = 0.2589
bref = 1.4224
cref = 0.3302
ref = Reference(Sref, bref, cref)
nothing # hide

controller = ConstantController(0.0, 0.0, 0.0, 0.0, 0.8)
nothing # hide

Wi = [0.0, 0.0, 0.0]
Wb = [0.0, 0.0, 0.0]
rho = 1.2682
asound = 300.0
g = 9.81
atm = ConstantAtmosphere(Wi, Wb, rho, asound, g)
nothing # hide

CL0 = 0.09167 # Zero-alpha lift
CLalpha = 3.5016  # lift curve slope
CLq = 2.8932 # Pitch rate derivative
CLM = 0.0 # Mach derivative
CLdf = 0.0  # flaps derivative
CLde = 0.2724  # elevator derivative
CLmax = 1.4  # max CL (stall)
CLmin = -0.9  # min CL (neg stall)
alphas = 20*pi/180

CD0 = 0.01631  # zero-lift drag coerff
U0 = 10.0  # velocity corresponding to Reynolds number of CD0
exp_Re = -0.2  # exponent in Reynolds number calc
e = 0.8  # Oswald efficiency
Mcc = 0.7  # crest critcal Mach number
CDdf = 0.0  # flaps
CDde = 0.3045  # elevator
CDda = 0.0  # aileron
CDdr = 0.0  # rudder

CYbeta = -0.07359 # Sideslip derivative
CYp = 0.0  # Roll rate derivative
CYr = 0.0 # Yaw rate derivative
CYda = 0.0 # Roll control (aileron) derivative
CYdr = 0.0 # Yaw control (rudder) derivative

Clbeta = -0.02854  # Sideslip derivative
Clp = -0.3209  # Roll rate derivative
Clr = 0.03066  # Yaw rate derivative
Clda = 0.1682  # Roll (aileron) control derivative
Cldr = 0.0  #Yaw (rudder) control derivative

Cm0 = -0.02338 # Zero-alpha pitch
Cmalpha = -0.5675 # Alpha derivative
Cmq = -1.3990 # Pitch rate derivative
CmM = 0.0
Cmdf = 0.0
Cmde = -0.3254 # Pitch control derivative

Cnbeta = -0.00040  # Slideslip derivative
Cnp = -0.01297  # Roll rate derivative
Cnr = -0.00434  # Yaw rate derivative
Cnda = -0.00328  # Roll (aileron) control derivative
Cndr = 0.0  # Yaw (rudder) control derivative

sd = StabilityDeriv(CL0, CLalpha, CLq, CLM, CLdf, CLde, alphas,
    CD0, U0, exp_Re, e, Mcc, CDdf, CDde, CDda, CDdr,
    CYbeta, CYp, CYr, CYda, CYdr,
    Clbeta, Clp, Clr, Clda, Cldr,
    Cm0, Cmalpha, Cmq, CmM, Cmdf, Cmde,
    Cnbeta, Cnp, Cnr, Cnda, Cndr)
nothing # hide

CT0 = 0.11221
CT1 = -0.13803
CT2 = -0.047394
CQ0 = 0.0062
CQ1 = 0.00314
CQ2 = -0.015729
D = 10*0.0254
num = 2
type = COCOUNTER
R = 0.5
Kv = 2500.0 * pi/30
i0 = 0.3
voltage = 8.0
propulsion = MotorPropBatteryDataFit(CT2, CT1, CT0, CQ2, CQ1, CQ0, D, num, type, R, Kv, i0, voltage)
nothing # hide

inertial2 = UniformGravitationalField()
nothing # hide

Vinf = U0
alpha = 3.0*pi/180
s0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, Vinf*cos(alpha), 0.0, Vinf*sin(alpha), 0.0, 0.0, 0.0]
tspan = (0.0, 4.0)
p = mp, ref, sd, propulsion, inertial2, atm, controller

# using DifferentialEquations
# prob = DifferentialEquations.ODEProblem(sixdof!, s0, tspan, p)
# sol = DifferentialEquations.solve(prob)
# nothing # hide

#Tyson testing his stuff
JACOBIAN=finite_jacobian(s0,p)
s0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, Vinf*cos(alpha), 0.0, Vinf*sin(alpha), 0.0, 0.0, 0.0]
p = mp, ref, sd, propulsion, inertial2, atm, controller
NEWJACOBIAN=getjacobian(s0,p)

#=
struct wing{TV, TF}
    pos::TV
    dir::TV
    cr::TF
    ct::TF
    tilt::TF
    twist::TF
end
=#

leftWing = Wing([0.0 0.0 0.0],[-.1 -1 -.05], .3, .1, 0.0, 3.0)
rightWing = Wing([0.0 0.0 0.0],[-.1  1 -.05], .3, .1, 0.0, 3.0)

planeModel = (leftWing,rightWing)

animvtk(JACOBIAN, planeModel,"Zagi")
animvtk(NEWJACOBIAN, planeModel,"Zagi2")


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @save

# flags -
has_TM_flag = 0

# load training set -
full_training_data_frame = load_fibrinolysis_training_data()

# filter out the data -
experimental_data_table = filter(:visitid => x -> (x == 2 || x == 3 || x == 1), full_training_data_frame)
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:16]))
validation_data_array = convert.(Float32, Matrix(experimental_data_table[!, 17:20]))

# what is the size of the system?
(P, D) = size(experimental_data_table)

# setup sampling -
δ = 0.9
number_of_samples = 100
number_of_input_types = 14
input_range_array = Array{Float32,2}(undef, number_of_input_types, 2)
for i ∈ 1:number_of_input_types
    L = (1 - δ) * minimum(input_data[:, i])
    U = (1 + δ) * maximum(input_data[:, i])
    input_range_array[i, 1] = L
    input_range_array[i, 2] = U
end

# archive -
ensemble_archive = Array{Vector{Float32},1}()

# build an input array -
sample_input_array_random = Array{Float32,2}(undef, number_of_samples, number_of_input_types)
for i ∈ 1:number_of_samples
    for j ∈ 1:number_of_input_types
        L = input_range_array[j, 1]
        U = input_range_array[j, 2]
        r = rand()
        v = L + r * (U - L)
        sample_input_array_random[i, j] = v
    end
end

# main training loop -
for i ∈ 1:39

    # load a model -
    model_name = "deep_fibrinolysis_model-L$(i)O-TF-TM-$(has_TM_flag)-ALL.bson"
    model_file_path = joinpath(_PATH_TO_MODELS, "network-2-all-data-epoch-12k", model_name)
    @load model_file_path deep_fibrinolysis_model

    # sample -
    simulated_output_data = deep_fibrinolysis_model(input_data[i,:])
    push!(ensemble_archive, simulated_output_data)
end

S = transpose(hcat(ensemble_archive...))

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
type BSONIter
    _wrap_::Vector{UInt8}
    done::Bool

    function BSONIter(bsonObject::BSONObject)
        bsonIter = new(Array(UInt8, 128), false)
        ccall(
            (:bson_iter_init, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            bsonObject._wrap_
            ) || error("BSONIter(): failure")
        bsonIter.done = !ccall(
            (:bson_iter_next, libbson),
            Bool, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
        return bsonIter
    end

    function BSONIter(bsonObject::BSONObject, key::AbstractString)
        bsonIter = new(Array(UInt8, 128), false)
        ccall(
            (:bson_iter_init, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            bsonObject._wrap_
            ) || error("BSONIter(): failure")
        keyCStr = string(key)
        bsonIter.done = !ccall(
            (:bson_iter_find, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            keyCStr
            )
        return bsonIter
    end

    function BSONIter(bsonArray::BSONArray)
        bsonIter = new(Array(UInt8, 128), false)
        ccall(
            (:bson_iter_init, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            bsonArray._wrap_
            ) || error("BSONIter(): failure")
        bsonIter.done = !ccall(
            (:bson_iter_next, libbson),
            Bool, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
        return bsonIter
    end

    function BSONIter(bsonArray::BSONArray, key::Integer)
        bsonIter = new(Array(UInt8, 128), false)
        ccall(
            (:bson_iter_init, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            bsonArray._wrap_
            ) || error("BSONIter(): failure")
        keyCStr = string(key)
        bsonIter.done = !ccall(
            (:bson_iter_find, libbson),
            Bool, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            keyCStr
            )
        return bsonIter
    end
end
export BSONIter

# Index

Base.getindex(bsonObject::BSONObject, key::AbstractString) = begin
    bsonIter = BSONIter(bsonObject, key)
    bsonIter.done && error("key not found: $(repr(key))")
    value(bsonIter)
end

Base.getindex(bsonArray::BSONArray, key::Integer) = begin
    bsonIter = BSONIter(bsonArray, key - 1)
    bsonIter.done && error("key not found: $(repr(key))")
    value(bsonIter)
end

# Associative collection

Base.get!(bsonObject::BSONObject, key::AbstractString, default) = begin
    iter = BSONIter(bsonObject, key)
    iter.done || return value(iter)
    append(bsonObject, key, default)
    return default
end

# Iterator

start(bsonObject::BSONObject) = begin
    return BSONIter(bsonObject)
end
export start

next(bsonObject::BSONObject, bsonIter::BSONIter) = begin
    ((key(bsonIter), value(bsonIter)), next!(bsonIter))
end
export next

done(bsonObject::BSONObject, bsonIter::BSONIter) = begin
    bsonIter.done
end
export done

start(bsonArray::BSONArray) = begin
    return BSONIter(bsonArray)
end
export start

next(bsonArray::BSONArray, bsonIter::BSONIter) = begin
    (value(bsonIter), next!(bsonIter))
end
export next

done(bsonArray::BSONArray, bsonIter::BSONIter) = begin
    bsonIter.done
end
export done

# Private

function key(bsonIter::BSONIter)
    bsonIter.done && error("alredy done iteration")
    cstr = ccall(
        (:bson_iter_key, libbson),
        Ptr{UInt8}, (Ptr{UInt8}, ),
        bsonIter._wrap_
        )
    cstr == C_NULL && error("bson_iter_key: failure")
    String(unsafe_string(cstr))
end

function value_type(bsonIter::BSONIter)
    bsonIter.done && error("alredy done iteration")
    return ccall(
        (:bson_iter_type, libbson),
        BSONType, (Ptr{UInt8}, ),
        bsonIter._wrap_
        )
end

done(bsonIter::BSONIter) = bsonIter.done
export done

function value(bsonIter::BSONIter)
    ty = value_type(bsonIter)
    if ty == BSON_TYPE_DOUBLE
        return ccall(
            (:bson_iter_double, libbson),
            Cdouble, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
    elseif ty == BSON_TYPE_INT32
        return ccall(
            (:bson_iter_int32, libbson),
            Int32, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
    elseif ty == BSON_TYPE_INT64
        return ccall(
            (:bson_iter_int64, libbson),
            Int64, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
    elseif ty == BSON_TYPE_BOOL
        return ccall(
            (:bson_iter_bool, libbson),
            Bool, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
    elseif ty == BSON_TYPE_DATE_TIME
        ts = ccall(
            (:bson_iter_date_time, libbson),
            Int64, (Ptr{UInt8}, ),
            bsonIter._wrap_
            )
        return Dates.unix2datetime(ts / 1000)
    elseif ty == BSON_TYPE_NULL
        return nothing
    elseif ty == BSON_TYPE_MINKEY
        return :minkey
    elseif ty == BSON_TYPE_MAXKEY
        return :maxkey
    elseif ty == BSON_TYPE_UTF8
        return String(unsafe_string(ccall(
            (:bson_iter_utf8, libbson),
            Ptr{UInt8}, (Ptr{UInt8}, Ptr{UInt8}),
            bsonIter._wrap_,
            C_NULL
            )))
    elseif ty == BSON_TYPE_OID
        data = Array(UInt8, 12)
        ptr = ccall(
                (:bson_iter_oid, libbson),
                Ptr{UInt8}, (Ptr{UInt8},),
                bsonIter._wrap_
                )
        unsafe_copy!(pointer(data), ptr, 12)
        return BSONOID(data)
    elseif ty == BSON_TYPE_DOCUMENT
        length = Array(UInt32, 1)
        data = Array(Ptr{UInt8}, 1)
        ccall(
            (:bson_iter_document, libbson),
            Ptr{Void}, (Ptr{UInt8}, Ptr{UInt32}, Ptr{Ptr{UInt8}}),
            bsonIter._wrap_, length, data
            )
        return BSONObject(data[1], length[1], bsonIter)
    elseif ty == BSON_TYPE_ARRAY
        length = Array(UInt32, 1)
        data = Array(Ptr{UInt8}, 1)
        ccall(
            (:bson_iter_array, libbson),
            Ptr{Void}, (Ptr{UInt8}, Ptr{UInt32}, Ptr{Ptr{UInt8}}),
            bsonIter._wrap_, length, data
            )
        return BSONArray(data[1], length[1], bsonIter)
    else
        error("unhandled BSONType $ty")
    end
end
export value

function next!(bsonIter::BSONIter)
    bsonIter.done = !ccall(
        (:bson_iter_next, libbson),
        Bool, (Ptr{UInt8}, ),
        bsonIter._wrap_
        )
    bsonIter
end


text: auc_val(labels, scores) = auc_trapezoidal(roccurve(labels, scores)...)

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
module Draw

using PyPlot
using PyCall
patch = pyimport("matplotlib.patches")

#using ..Voronoi.Geometry
#using ..Voronoi.Diagram
#using ..Voronoi.BeachLine
#using ..Voronoi.EventQueue

import Voronoi

plt = PyPlot
ax = nothing
WIDTH = nothing
HEIGHT = nothing
FRAME = 0

#colors = ["white"]
colors = ["xkcd:grey", "xkcd:crimson", "xkcd:gold", "xkcd:green", "xkcd:azure", "xkcd:beige", "xkcd:silver", "xkcd:lavender", "xkcd:lightgreen", "xkcd:magenta", "xkcd:ivory", "xkcd:maroon", "xkcd:orange", "xkcd:orangered", "xkcd:orchid", "xkcd:pink", "xkcd:plum", "xkcd:gold", "xkcd:salmon", "xkcd:sienna", "xkcd:lime", "xkcd:tan", "xkcd:tomato", "xkcd:violet", "xkcd:wheat", "xkcd:indigo", "xkcd:yellowgreen", "xkcd:chocolate", "xkcd:coral", "xkcd:brown"]

function init(w, h)
   global plt, ax
   global WIDTH, HEIGHT

   plt.pygui(true)
   plt.ion()
   plt.clf()
   #plt.xkcd() # uncomment for generalized wobbliness
   ax = plt.gca() # get current axes

   WIDTH = w
   HEIGHT = h
end

function clear(title::String)
   global plt, ax
   global WIDTH, HEIGHT, FRAME

   plt.cla()
   plt.title(title)
   if WIDTH > 0
      ax.set_aspect("equal")
      ax.set_xlim([0 - FRAME, WIDTH + FRAME])
      ax.set_ylim([0 - FRAME, HEIGHT + FRAME])
      ax.grid(false)
      ax.get_xaxis().set_visible(false)
      ax.get_yaxis().set_visible(false)
   else
      ax.grid(true)
      ax.get_xaxis().set_visible(true)
      ax.get_yaxis().set_visible(true)
   end
end

function legend(xlabel, ylabel)
   global plt

   plt.xlabel(xlabel)
   plt.ylabel(ylabel)
   plt.legend()
end

function commit()
   global plt

   plt.draw()
end

function circle_Internal(p::Tuple{Real, Real}, color, fill, r, l)
   global ax

   ax.add_artist(patch.Circle((p[1], p[2]), color=color, radius=r, fill=fill, zorder=3, linewidth=l))
end

function point(p::Tuple{Real, Real}, color, fill)
   circle_Internal(p, color, fill, 0.5, 1)
end

function circle(p::Tuple{Real, Real}, color, r)
   circle_Internal(p, color, false, r, 1)
end

function line(p1::Tuple{Real, Real}, p2::Tuple{Real, Real}, color)
   global plt

   plt.plot([p1[1], p2[1]], [p1[2], p2[2]], color=color, linestyle="-", linewidth=3, zorder=1)
end

function thinLine(p1::Tuple{Real, Real}, p2::Tuple{Real, Real}, color)
   global plt

   plt.plot([p1[1], p2[1]], [p1[2], p2[2]], color=color, linestyle="-", linewidth=1, zorder=1)
end

function plot(f, color, start, finish)
   global plt

   x = range(start, stop=finish, length=1000)
   y = map(f, x)

   plt.plot(x, y, color=color, linewidth=3, zorder=2)
end

function plot(x, y, color, label, style)
   global plt

   plt.plot(x, y, color=color, label=label, linestyle=style, linewidth=3, zorder=2)
end

function fortuneIteration(V::Voronoi.Diagram.DCEL, T::Voronoi.BeachLine.BST, Q::Voronoi.EventQueue.Heap, points::Array{Tuple{Real, Real}}, ly::Real)
   #clear("Computing Voronoi Voronoi.Diagram using Fortune's Algorithm")
   clear("")

   # draw points
   for p in points
      point(p, "xkcd:black", (ly <= p[2]))
      if ly < p[2]
         f = Voronoi.Geometry.parabola(p, ly)
         plot(f, "xkcd:silver", 0, WIDTH)
      end
   end

   # calculate parabolas and breakpoints
   beachLineFoci = Voronoi.BeachLine.beachLine(T, ly)

   ### DEBUGGING SECTION
   #Voronoi.BeachLine.printTree(T)
   #println(beachLineFoci)

   ## print x foci of beach line, without breakpoints
   #i = 1
   #while i <= size(beachLineFoci)[1]
   #   print(beachLineFoci[i][1], " ")
   #   i += 2
   #end
   #println()

   ## print forward list of beachline arcs x's
   #a = T.root
   #if a != nothing
   #   while !isa(a, Voronoi.BeachLine.Arc)
   #      a = a.leftChild
   #   end
   #   while a != nothing
   #      print(a.focus[1], " ")
   #      a = a.next
   #   end
   #end
   #println()

   ## print backwards list of beachline arcs x's
   #a = T.root
   #if a != nothing
   #   while !isa(a, Voronoi.BeachLine.Arc)
   #      a = a.rightChild
   #   end
   #   while a != nothing
   #      print(a.focus[1], " ")
   #      a = a.prev
   #   end
   #end
   #println()
   ### DEBUGGING SECTION END

   # draw beachline
   start = (0, HEIGHT)
   i = 2
   while i <= size(beachLineFoci)[1] + 2
      if i <= size(beachLineFoci)[1]
         finish = beachLineFoci[i]
      else
         finish = (WIDTH, start[2])
      end

      p = beachLineFoci[i-1]
		f = Voronoi.Geometry.parabola(p, ly)

		if f == nothing # point is over the sweep line
         if start == nothing # special case where there the first couple of points are on the same y coordinate
         else
            Draw.line((p[1], ly), (p[1], start[2]), "xkcd:azure")
            #Draw.line((p[1], ly), (p[1], start[2]), "xkcd:black")
         end
		else
         if 0 <= start[1]
            st = start[1]
         else
            st = 0
         end

         if finish[1] <= WIDTH
            fn = finish[1]
         else
            fn = WIDTH
         end

			plot(f, "xkcd:azure", st, fn)
			#plot(f, "xkcd:black", st, fn)
		end

      start = finish
      i += 2
	end

   # draw circumcircles
   i = 1
   while i < Q.pos
      if isa(Q.data[i], Voronoi.EventQueue.CircleEvent) && !(Q.data[i].removed)
         b = Q.data[i].disappearingArc
         a = b.prev
         c = b.next
         O = Q.data[i].center
         r = O[2] - Q.data[i].coordinates[2]
         circle(O, "xkcd:orangered", r)
         thinLine(O, b.region.generator, "xkcd:orangered")
         point(O, "xkcd:magenta", true)
         point((O[1], O[2] - r), "xkcd:orangered", true)
      end

      i += 1
   end

   # draw diagram edges
   for he in V.halfEdges
      if he.origin != nothing && he.twin.origin != nothing
         Draw.line(he.origin, he.twin.origin, "xkcd:black")
      end
   end

   # draw sweepline
	line((0, ly), (WIDTH, ly), "xkcd:gold")
	#line((0, ly), (WIDTH, ly), "xkcd:black")
	commit()
end

function voronoiDiagram(V::Voronoi.Diagram.DCEL)
    global WIDTH, HEIGHT

    #clear("Computing Centroidal Voronoi Tesselations using Lloyd's Algorithm")
	clear("")

    regions = Voronoi.Diagram.regionBorders(V)

    # draw regions
    i = 1
    for region in regions
        if size(region[1])[1] > 0
            #color = (V.regions[i].generator[1]/WIDTH, V.regions[i].generator[2]/HEIGHT, 0.5)
            #ax[:fill](region[1], region[2], color=color)
            ax.fill(region[1], region[2], colors[(i % size(colors)[1])+ 1])
        end
        i += 1
    end

    # draw diagram edges
    for he in V.halfEdges
        if he.origin != nothing
            if he.twin != nothing
                Draw.line(he.origin, he.twin.origin, "xkcd:black")
            elseif he.next != nothing
                Draw.line(he.origin, he.next.origin, "xkcd:black")
            end
        end
    end

    # draw generators
    for r in V.regions
        point(r.generator, "xkcd:black", true)
    end

    line((0, 0), (0, HEIGHT), "cyan")
    line((0, 0), (WIDTH, 0), "cyan")
    line((WIDTH, 0), (WIDTH, HEIGHT), "cyan")
    line((0, HEIGHT), (WIDTH, HEIGHT), "cyan")

    commit()
end

end # module

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
"""
    niw_hyperparams(κ::Float32, m::AbstractArray{Float32}, ν::Float32, ψ::AbstractArray{Float32})

[Normal Inverse Wishart](https://en.wikipedia.org/wiki/Normal-inverse-Wishart_distribution)
"""
struct niw_hyperparams <: distribution_hyper_params
    κ::Float32
    m::AbstractArray{Float32}
    ν::Float32
    ψ::AbstractArray{Float32}
end

mutable struct niw_sufficient_statistics <: sufficient_statistics
    N::Float32
    points_sum::AbstractArray{Float32,1}
    S::AbstractArray{Float32,2}
end

function calc_posterior(prior:: niw_hyperparams, suff_statistics::niw_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    κ = prior.κ + suff_statistics.N
    ν = prior.ν + suff_statistics.N
    m = (prior.m.*prior.κ + suff_statistics.points_sum) / κ
    ψ = (prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m'+ suff_statistics.S) / ν
    ψ = Matrix(Hermitian(ψ))
    return niw_hyperparams(κ,m,ν,ψ)
end

function sample_distribution(hyperparams::niw_hyperparams)
    Σ = rand(Distributions.InverseWishart(hyperparams.ν, hyperparams.ν* hyperparams.ψ))
    μ = rand(Distributions.MvNormal(hyperparams.m, Σ/hyperparams.κ))
    invΣ = inv(Σ)
    chol = cholesky(Hermitian(invΣ))
    return mv_gaussian(μ,Σ,invΣ,logdet(Σ),chol.U)
end

function create_sufficient_statistics(hyper::niw_hyperparams,posterior::niw_hyperparams,points::AbstractArray{Float32,2}, pts_to_group = 0)
    if size(points,2) == 0
        return niw_sufficient_statistics(size(points,2),zeros(Float32,length(hyper.m)),zeros(Float32,length(hyper.m),length(hyper.m)))
    end
    pts = copy(points)
    points_sum = sum(pts, dims = 2)[:]
    S = pts * pts'
    return niw_sufficient_statistics(size(points,2),points_sum,S)
end

function log_marginal_likelihood(hyper::niw_hyperparams, posterior_hyper::niw_hyperparams, suff_stats::niw_sufficient_statistics)
    D = length(suff_stats.points_sum)
    logpi = log(pi)
    return -suff_stats.N*D*0.5*logpi +
        log_multivariate_gamma(posterior_hyper.ν/2, D)-
        log_multivariate_gamma(hyper.ν/2, D) +
         (hyper.ν/2)*(D*log(hyper.ν)+logdet(hyper.ψ))-
         (posterior_hyper.ν/2)*(D*log(posterior_hyper.ν) + logdet(posterior_hyper.ψ)) +
         (D/2)*(log(hyper.κ/posterior_hyper.κ))
end

function aggregate_suff_stats(suff_l::niw_sufficient_statistics, suff_r::niw_sufficient_statistics)
    return niw_sufficient_statistics(suff_l.N+suff_r.N, suff_l.points_sum + suff_r.points_sum, suff_l.S+suff_r.S)
end


text:
macro _const(ex)
    ex = esc(ex)
    if VERSION < v"1.8.0-DEV.1148"
        return ex
    else
        return Expr(:const, ex)
    end
end

const var"@const" = var"@_const"

#=
const Historic = if USE_HISTORIC
    Base.require(Base.PkgId(Base.UUID(0xe6ec0b50ef98488aa141280ae2eaf113), "Historic"))
else
    nothing
end

@static if Historic === nothing
    macro record(_...)
        nothing
    end
else
    Historic.@define Debug
    using .Debug: @record
end

"""
    @yield_unsafe expression

Document that `expression` must not yield to the Julia scheduler.

This macro is used purely for documentation.  The `expression` is evaluated
as-is.
"""
macro yield_unsafe(ex)
    esc(ex)
end
=#

function ConcurrentUtils.spinloop()
    GC.safepoint()
    ccall(:jl_cpu_pause, Cvoid, ())
end

function ConcurrentUtils.spinfor(nspins)
    for _ in oneto(nspins)
        spinloop()
    end
end

oneto(::Nothing) = ()
oneto(n) = Base.OneTo(n)

struct Infinity end
const ∞ = Infinity()
Base.:<(::Int, ::Infinity) = true

struct Forever end
Base.iterate(::Forever, _state = nothing) = (nothing, nothing)
Base.eltype(::Type{Forever}) = Nothing
Base.IteratorSize(::Type{Forever}) = Base.IsInfinite()

oneto(::Infinity) = Forever()

struct NegativeInfinity end
Base.:-(::Infinity) = NegativeInfinity()
Base.:<(::Int, ::NegativeInfinity) = false

fieldoffset_by_name(T, name) = fieldoffset(T, findfirst(==(name), fieldnames(T)))

_typeof(x) = typeof(x)
_typeof(::Type{T}) where {T} = Type{T}

pad7() = ntuple(_ -> 0, Val(7))


text:
using Combo
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
srand(1234)
a = rand(1:2, 10)
seq, val = cakewalk(seq -> sum((seq - a).^2), 2, length(a), verbose = false);
@test val[1] == Combo.ZERO