text:
module Test

using Test

#include("../src/bracket.jl")
include("../src/bracket_point.jl")

function test(func,result)
    vm = Vm()
    #reset!(vm)
    # load prelude
    vm.bra = load_file(vm,"../src/prelude.clj")
    vm.ket = NIL
    eval_bra(vm)

    #vm.depth = 1  # normal programs start in depth 1
    vm.ket = NIL
    vm.bra = make_bra(vm,func)
    res, b =reverse_list(vm, make_bra(vm,result))
    eval_bra(vm)
    #println("ket is ")
    #printKet(vm,vm.ket)
    #println("res is")
    #printKet(vm,res)
    isEqual(vm,res,vm.ket)
end

function tests()

   println("Bracket Tests")

   @testset "Stack shuffling and list manipulation" begin

   @test test("1 2 3",     "1 2 3")  # values of bra are shifted on ket
   @test test("1 2 3; this is a comment",  "1 2 3")  # values of bra are shifted on ket
   @test test("1 [ ]",     "1 [ ]")  # values of bra are shifted on ket
   @test test("[1 2 3]",   "[1 2 3]")  # values of bra are shifted on ket
   @test test("1 [a ; comment \n [+ 2 3 ]]",  "1 [a [+ 2 3 ]]")

   @test test("x esc 4",   "x 4")      # escape a symbol
   @test test("x' 4",      "x 4")      # short notation
   #@test test("[1 2 3]'",  "[3 2 1]")  # escaping a list -> reverse
   @test test("[1 2 3]'",  "[1 2 3]")   # list and numbs: ..
   @test test("3'",        "3")         # .. just move to the ket

   # Stack shuffling
   @test test("dup 2",     "2 2")
   @test test("dup 2 3",   "2 2 3")
   @test test("dup [2 3]", "[2 3] [2 3]")
   @test test("drop 2 3",  "3")
   @test test("drop 2",    "")

   @test test("swap 2 3",  "3 2")
   @test test("swap f' g'",  "g f")
   @test test("swap [] [1 2]",  "[1 2] []")

   @test test("rot 1 2 3",  "3 1 2")
   @test test("rot 1",         "")
   @test test("rot [x y] f' 3", "3 [x y] f")

   # car, cdr, cons
   #@test test("car [1 2 3] 4", "3 [1 2] 4")
   @test test("car [1 2 3] 4", "3 4")
   #@test test("car [x]",      "x []")   # bring a symbol on ket
   @test test("car [x]",      "x")   # bring a symbol on ket
   @test test("car [] 10", "10")        # car empty list
   @test test("car 1", "")              # car atom
   @test test("car x' def x' [1 2]",  "2")          # car a symbol ..
   @test test("car x' car x' def x' [1 2]",  "2 2") # .. leaves the symbol intact
   @test test("car x' def x' 1",  "")

   @test test("cdr [1 2 3] 4",   "[1 2] 4")
   @test test("cdr []",        "[]")     # cdr empty list
   @test test("cdr 1",         "[]")     # cdr atom
   @test test("cdr x' def x' [1 2 3]",  "[1 2]")           # cdr a symbol ..
   @test test("x` cdr x' def x' [1 2 3]", "[1 2 3] [1 2]") # .. leaves symbol intact
   @test test("cdr foo' def foo' bar'",  "[]")

   @test test("cons 1 []",          "[1]")
   @test test("cons 4 [1 2 3]",     "[1 2 3 4]")
   @test test("cons [4 5] [1 2 3]", "[1 2 3 [4 5]]")
   #@test test("cons car [1 2 3 4]",  "[1 2 3 4]") # identity
   @test test("cons car swap cdr dup [1 2 3]", "[1 2 3]")  # identity
   @test test("cons car x' cdr x' def x' [1 2 3 4]",  "[1 2 3 4]") # identity
   # cons into variable
   #@test test("cons 4 x' def x' [1 2 3]",  "[x . 4]")     # dotted list ..

   # no dotted list
   @test test("cons 1 2", "[2 1]") # consing two atoms automatically creates a list
   @test test("cons 1 x' ", "[x 1]")

   # def, val
   @test test("def x' 2 10",  "10")   # x bound to 2, def consumes value on ket
   @test test("x",            "[]")   # unbound variable evaluates to NIL
   @test test("x def x' 2",   "2")
   @test test("x y x def y' 3 def x' 2",  "2 3 2")

   @test test("foo-bar def foo-bar' 2", "2")         # dash in symbol name
   @test test("foo_bar def foo_bar' 2", "2")         # underscore in symbol name
   @test test("foo1 def foo1' 2", "2")               # digit in symbol name
   @test test("Fo-13g def Fo-13g' 2", "2")           # digit in symbol name
   @test test("f123456789  def f123456789' 2", "2")  # symbol name with 10 characters

   @test test("val x' def x' 2",       "2")
   @test test("x vesc def x' 2",       "2")  # vesc, escape value of symbol to ket
   @test test("x` def x' 2",           "2")  # backtick = short for vesc
   @test test("val [1 2]",             "[1 2]")
   @test test("[1 2]`",                "[1 2]")
   @test test("x` def x' [1 2 3]",     "[1 2 3]")
   @test test("val x' def x' [1 2 3]", "[1 2 3]")
   @test test("x def x' [1 2 3]",      "1 2 3")
   @test test("val x' def x' 3 x ` def x' 2",  "3 2")

   @test test("f 2 3 def f' [add]",      "5")
   @test test("foo def foo' [add 1] 2",  "3")
   @test test("foo 1 2 def foo' add'", "add 1 2") # def a symbol that is a builtin
   @test test("foo 1 2 def foo' [add]", "3")      # [add] is quotation

   @test test("a def [a] 2",             "2")    # def: bind a list of keys
   @test test("a b def [a b] 1 2",       "1 2")
   @test test("a b c def [a b c] 1 2 3", "1 2 3")
   @test test("a b c def [a b c] 1 2",   "[] 1 2")
   @test test("a b` def [a b] 1 [2 3]",  "1 [2 3]")
   @test test("def [] 1 2",              "2")

   @test test("a` b` def [[a b]] 2", "2 2")       # pattern matching
   @test test("a` b` c` d` def [[a b [c d]]] 2", "2 2 2 2")
   @test test("b` def [[b]] [1 2 3]", "[1 2 3]")
   @test test("a` b` def [[a b]] [1 2 3]", "[1 2] 3")
   @test test("a` b` c` def [[a b c]] [1 2 3]", "[1] 2 3")

   @test test("a def [a`] 2",             " 2")    # backtick interpreted as set
   @test test("a b def [a` b] 2 3",       " 2 3")
   @test test("b a b def [a` b`] 2 3",    "3 2 3")
   @test test("f 2 3 def [f`] [add]",     "5")
   @test test("eval [x def [x`] 2]",      "2")

   #test("x set x' 2",  "2")                  # set as a primary operator, obsolete now
   #test("x y x set y' 3 set x' 2",  "2 3 2")
   #test("f 2 3 set f' [add]",     "5")
   #test("eval [x set x' 2]",  "2")

  # numbers as symbols
  # @test test("eval 10 def 10 5",        "5")
  # @test test("eval 10 def [10] 5",      "5")
  # @test test("eval 10 def 10 5",        "5")
  # @test test("eval 10 1 2 def 10 [add]", "3")
  # @test test("val 10 1 2 def 10 [add]", "[add] 1 2")
  # @test test("10` 1 2 def 10 [add]",    "[add] 1 2")
  # @test test("val 10 1 2 def 10 add'",  "add 1 2")

   @test test("eval [x def x' 2]",  "2")
   @test test("eval [x def x' 2] def x' 3",  "2")          # local scope
   @test test("eval [x] def x' 2]",  "2")    # inner scope can use value defined outside
   #@test test("eval [x] set x' 2]",  "2")    # inner scope can use value defined outside
   @test test("eval [x] def [x`] 2]", "2")    # inner scope can use value defined outside
   @test test("x eval [x def x' 2] x def x' 3",  "3 2 3")  # def changes only within scope
   #@test test("x eval [x set x' 2] x set x' 3",  "2 2 3")  # set changes also outside
   @test test("x eval [x def [x`] 2] x def [x`] 3",  "2 2 3")  # set changes also outside
   @test test("x eval [x def [x`] 2] x def x' 3",  "2 2 3")  # set changes also outside

   @test test("eval [add 1] 2",          "3")
   @test test("eval [add 1 2]",          "3")
   @test test("eval [1 2 3] 4",          "1 2 3 4")
   @test test("eval [] 1",               "1") # eval empty list
   @test test("eval add' 1 2",           "3")
   @test test("eval foo' def foo' 5",    "5")
   @test test("eval foo' def foo' bar'", "bar")

   # lambda and lexical scoping
   @test test("lambda x' [+ x 1]","[+ x 1 def x']")
   @test test("f def f' lambda x' [+ x 1] 10","11")
   @test test("eval lambda x' [+ x 1] 10","11")
   # alternative use interpreted strings and "\\""
   @test test("eval \\x' [+ x 1] 10","11")

   @test test("eval \\[x] [+ x 1] 10","11")      # lambda with list of arguments
   @test test("\\[x y] [- x y]","[- x y def [x y]]")
   @test test("eval \\[x y] [- x y] 10 2","8")

   @test test("eval foo` def foo' \\[x] [+ x 1] 10","11")      # backquote or val ..
   @test test("eval val foo' def foo' \\[x] [+ x 1] 10","11")  # .. puts quotqtion of closure on ket
   #@test test("eval val foo' def foo' \\[x] [+ x 1] 10","11")  # val on closure puts quotation on stack
   @test test("eval eval [       [x def x' + 1 x`] def x' 10] def x' 1", "2")
   @test test("eval eval [\\ [] [x def x' + 1 x`] def x' 10] def x' 1", "11")
   @test test("eval swap eval dup eval [\\[] [x def [x`] + 1 x`]] def x' 10",         "12 11")
   @test test("x foo foo def foo' eval [\\[] [x def [x`] + 1 x`]] def x' 10",         "12 12 11")
   #@test test("x foo foo def foo' eval [\\[] [x set x' + 1 x`]]  def x' 10",          "12 12 11")
   @test test("x foo foo def foo' eval [\\[] [x def x' + 1 x`]]  def x' 10",          "10 11 11")
   @test test("x foo foo def foo' eval [\\[] [x def [x`] + 1 x`] def x' 10] def x' 1", "1 12 11")
   @test test("x foo foo def foo' eval [\\[] [x def x' + 1 x`]   def x' 10] def x' 1", "1 11 11")
   @test test("x foo foo def foo' [x def x' + 1 x`]   def x' 10",                     "10 11 11")
   @test test("x foo foo def foo' [x def [x`] + 1 x`] def x' 10",                     "12 12 11")

   @test test("100 eval [f def x' 20] def f' \\[] [x] def x' 10","100 10")
   @test test("100 eval [f def x' 20] def f'     [x] def x' 10","100 20")
   @test test("100 g def g' [f def x' 20] def f' \\[] [x] def x' 10","100 10")
   @test test(" g g def g' [f def x' 20] def f' \\[] [x] def x' 10","20 10")
   # the 2nd evaluation of g is the last command and due to tail elimination
   #  f is evaluated in the global scope

   # closure
   @test test("a 2 a 3 def a' [make-adder 4] a 2 a 3 def a' [make-adder 5]
       def make-adder' [addx def x']
       def addx' [+ x z def z']" , "6 7 7 8")

   @test test("add1 10 add1 20 def add1' eval \\[x] [\\[][+ x]] 1", "11 21")

   @test test("add1 add1 add1 def add1' eval [\\[][x def [x`] + x 1] def x']
   0", "3 2 1")

   end # testset

# ------------------------------------
   @testset "Math and logic" begin
   # Math
   @test test("add 2 2",            "4")
   @test test("+ 2 2",             "4")
   @test test("sub 5 2",           "3")
   @test test("- 5 2",             "3")
   @test test("mul 2 3",           "6")
   @test test("* 2 3",             "6")
   @test test("div 9 2",           "4")
   @test test("/ 8 2",             "4")
   @test test("/ 9 2",             "4")
   @test test("/ 9 0",             "0")   # division by zero returns 0
   @test test("- * 2 4 10",        "-2")
   #@test test("+ 2.0 2.0",         "4.0")
   #@test test("+ 2.0 2",           "4.0")
   #@test test("- 2.0 2",           "0.0")
   #@test test("* 3.14 -2.0",       "-6.28")
   #@test test("/ 3.14 2.0",        "1.57")
   #@test test("/ 3.14 2",          "1.57")
   #@test test("/ 3.14 0.0",        "0.0")  # division by zero
   #@test test("/ 3.14 0",          "0.0")  # division by zero
   #@test test("/ 3 0.0",           "0.0")  # division by zero

   @test test("+ 2 x' def x' 3",     "5")
   @test test("- 5 x' def x' 2",     "3")
   @test test("- y' 2 def y' 5",     "3")
   @test test("- y' x' def y' 5 def x' 2",  "3")

   @test test("+ [2 3] 1", "[3 4]")  # add to a list
   @test test("- [2 3] 1", "[1 2]")
   @test test("+ 1 [2 3]", "[3 4]")  # add to a list
   @test test("- 4 [2 3]", "[2 1]")
   @test test("+ [10 20][1 2]", "[11 22]")
   @test test("+ [10 20][3 1 2]", "[11 22]") # list of different length
   @test test("+ x' 2 def [x] [5 6]", "[7 8]")
   @test test("+ x' 2 def x' 3", "5")
   @test test("+ 2 x' def x' 3", "5")
   @test test("- x' y' def [x y] 5 3", "2")
   @test test("- x' y' def [x y] [5 6] 3", "[2 3]")
   @test test("- [5 6] x' def [x] 3", "[2 3]")
   @test test("foo foo def foo' [+ 1] 2", "4")

   @test test("lt 4 10",               "1")
   @test test("lt 10 4",               "0")
   @test test("lt 4 4",                "0")
   @test test("gt 10 4",               "1")
   @test test("gt 10 10",              "0")
   @test test("gt 4 10",               "0")
   @test test("lt 5 [4 5 6 7]" , "[0 0 1 1]")
   #@test test("lt 4.0 10",             "1")
   #@test test("lt 4.0 10.0",           "1")
   #@test test("lt 4 10.0",             "1")
   #@test test("lt 5 [4 5 6 7]",        "[0 0 1 1]")
   #@test test("gt 4 [-2.0 3 4.0 5]",   "[1 1 0 0]")
   #@test test("lt [1 1.0] [0.2 0.2 0.2]",  "[0 0]")

   # eq
   @test test("eq 2 2",                "1")
   @test test("eq 2 3",                "0")
   @test test("eq [] []",              "1")
   @test test("eq [1 2 3] [1 2 3] ",   "1")
   @test test("eq [1 2 3] [1 2 4] ",   "0")
   @test test("eq 2 []",               "0")
   @test test("eq 2 x'",               "0")
   @test test("eq x' x'",              "1")
   @test test("eq y' x'",              "0")
   @test test("eq y' x' def y' x'",    "0")
   #@test test("eq 2 car [eq]",         "0")
   #@test test("eq car [eq] car [eq]",  "1")
   #@test test("eq car [car] car [eq]", "0")
   #@test test("eq [car [1 4] 2] [car [1 4] 2]", "1")

   # logic, if
   @test test("if 1 20 30",            "20")
   @test test("if 0 20 30",            "30")
   @test test("if [7] 20 30",          "20") # nonempty list is true
   @test test("if [] 20 30",           "30") # empty list is false
   @test test("if foo' 20 30",         "20") # any symbol not Nil is true
   @test test("if 20 30",              "")
   @test test("if 20",                 "")
   @test test("if 1 foo' bar'",        "foo")

   # dip
   @test test("dip [+ 1] 5 2", "5 3")
   @test test("dip [+ 1] [+ 10] 2", "[+ 10] 3")
   @test test("dip [1 2 3] 4", "4 1 2 3")

   # cond
   #@test test("cond [[2]]",       "2")
   #@test test("cond [[10] [11] [0]]",   "10")
   #@test test("cond [[10] [11] [1]]",   "11")
   #@test test("cond [[3] [1] [eq 4] [2 drop] [lt 4 dup]] 3", "3")
   #@test test("cond [[3] [1] [eq 4] [2 drop] [lt 4 dup]] 5", "2")
   #@test test("cond foo' 5 def foo' [[3] [1] [eq 4] [2 drop] [lt 4 dup]]", "2")

   @test test("typ 20",     "1")
   @test test("typ add'",   "2")
   @test test("typ foo'",   "3")
   @test test("typ [2 3]",  "4")
   @test test("typ []",     "2")   # == NIL
   @test test("typ \\[x][2 3]",  "5")

   end # testset

# ------------------------------------
   @testset "prelude" begin

   # @test test("splt [1 2 3 4]", "4 [1 2 3]")
   @test test("over 1 2", "2 1 2")
   @test test("over 1"  , "")
   @test test("over"    , "")
   @test test("rot1 1 2 3", "2 3 1")
   @test test("rot1 1 2 [+ 4]", "2 [+ 4] 1")
   @test test("rot1 1",         "")
   @test test("drop2 10 11 12", "12")
   @test test("drop3 10 11 12", "")
   @test test("dup2 2 3", "2 3 2 3")
   @test test("dupd 2 3", "2 3 3")
   @test test("nip 2 3 4", "2 4")
   @test test("nip2 2 3 4", "2")
   @test test("swapd 2 3 4", "2 4 3")
   @test test("dupd 2 3", "2 3 3")
   @test test("dup2 2 3", "2 3 2 3")
   @test test("rot4 1 2 3 4", "4 1 2 3")
   @test test("rot14 1 2 3 4", "2 3 4 1")
   @test test("splt [1 2 3]", "3 [1 2]")
   @test test("cons splt [1 2 3]", "[1 2 3]") # identity

   @test test("dip2 [+ 1] 1 2 3", "1 2 4")
   @test test("keep [+] 2 3", "2 5")
   #@test test("keep2 [+] 2 3", "2 3 5")

   #@test test("+ 2 3 meta stop 10", "10")   // stop execution
   #@test test("empty_ket 1 2 3","")         // set ket to empty list
   #@test test("print_list [x y] def x' 2 def y' 3", "")  // print list

   @test test("curry [eq] foo'", "[eq foo']")
   @test test("eqfoo bar' eqfoo foo' def eqfoo' curry [eq] foo'", "0 1")

   @test test("not 5", "0")
   @test test("not 0", "1")
   @test test("and 1 0", "0")
   @test test("and 0 1", "0")
   @test test("and 2 5", "5")
   @test test("and [] 1", "[]")
   @test test("or 0 1", "1")
   @test test("or 1 0", "1")
   @test test("or 0 0", "0")
   @test test("or 2 5", "2")
   @test test("when 1 [+ 10] 20","30")
   @test test("when 0 [+ 10] 20","20")
   @test test("when1 6 [+ 10]","16")  # retains the argument of the logical decision
   @test test("when1 0 [+ 10]","")
   @test test("unless 1 [+ 10] 20","20")
   @test test("unless 0 [+ 10] 20","30")
   #@test test("reverse [1 2 3]","[3 2 1]")

   @test test("caar [6 [4 5][1 2 3]", "3")
   @test test("cadr [6 [4 5][1 2 3]", "[4 5]")
   @test test("cdar [6 [4 5][1 2 3]", "[1 2]")
   @test test("cddr [6 [4 5][1 2 3]", "[6]")

   @test test("keep [+ 1] 2","2 3")
   @test test("keep2 [+] 2 3","2 3 5")
   @test test("bi [* dup][+ 1] 2","4 3")
   @test test("bi2 [*][+] 3 4","12 7")
   @test test("tri [- 1] [* dup][+ 1] 2","-1 4 3")
   @test test("tri2 [-][*][+] 3 4","-1 12 7")
   @test test("bistar [* dup][+ 1] 3 2","9 3")
   @test test("bi2star [*][+] 4 3 2 1","12 3")

   @test test("cleave [[* dup][+ 1]] 2","4 3")
   @test test("cleave [[- 1] [* dup][+ 1]] 2","-1 4 3")
   @test test("cleave2 [[-][*][+]] 3 4","-1 12 7")

   @test test("each [* dup] [4 3 2 1]", "16 9 4 1")
   #@test test("map [* dup] [4 3 2 1]", "[16 9 4 1]")
   @test test("map [* dup] [4 3 2 1]", "[1 4 9 16]")    # reverse still missing
   #@test test("map1 [* dup] [4 3 2 1]", "[1 4 9 16]")
   @test test("unstack [4 3 2 1]", "4 3 2 1")

   @test test("sum [2 5 10]", "17")
   @test test("prod [2 5 10]", "100")
   @test test("size [2 5 foo [3 4] 10]", "5")
   @test test("repeat 4 [+ 2] 0", "8")
   @test test("filter [gt swap 0] [2 -1 5] ", "[5 2]")  # reverse still missing
   @test test("filter [gt swap 0] [-2 -1 -10] ", "[]")  # reverse still missing
   @test test("drop drop loop [lt 0 dup - swap 1 keep [*]] 4 1", "24")

   end

# ------------------------------------
   @testset "Recur and small examples" begin

   # recur
   @test test("eval [rec gt 0 dup add 1 dup] -5", "0 -1 -2 -3 -4 -5")   #simple loop
   @test test("foo def foo' [rec gt 0 dup add 1 dup] -5", "0 -1 -2 -3 -4 -5")  #simple loop

   # simple closure

   @test test("eval f 3 2
       a 5 a 6
       def a' f 1
       def f' \\[x][
         \\[y][+ x y]]","5 6 7")

   # simple closure for bank account
   @test test("
    acc withdraw' 60
    acc deposit' 100 acc withdraw' 60 acc withdraw' 60 acc deposit' 40
    def acc' make-acc 50
    def make-acc' [
        \\[m][eval if eq m withdraw'
                [withdraw]
            [eval if eq m deposit'
                [deposit]
                [unknown']]]
        def withdraw' [
            eval if gt balance rot
            [balance def [balance`] - balance]
            [insuff' drop] dup]
        def deposit' [balance def [balance`] + balance]
        def balance' ]"    , "70 130 insuff 30 90")

   #more fun with closures
   @test test("
    show-bal
    eval deposit1 acc' 3
    def deposit1' \\[ac]
      [\\[def [bal`] + bal] ac`]
    show-bal deposit 5
    def deposit' \\[def [bal`] + bal] acc'
    show-bal
    def show-bal' \\[bal] acc'
    acc
    def acc' make-acc 10
    def make-acc' [
       \\[][do-stuff']
       def bal'
    ]", "18 15 10 do-stuff")

    # Factorial
    #simple recursive
    #@test test("fac 4 def fac' [cond [[* fac - swap 1 dup] [1 drop] [eq 1 dup]]]", "24")
    @test test("fac 4 def fac' [eval if eq 1 rot [1 drop] [* fac - swap 1 dup] dup]", "24")
    #@test test("fac 4 def fac' \\[n] [cond [[* fac - n 1 n] 1 [eq 1 n]]]", "24")
    @test test("fac 4 def fac' \\[n] [eval if eq 1 n 1 [* fac - n 1 n]] " , "24")

    # tail recursive
    @test test("fac 4 def fac' \\[n]
        [drop swap eval \\[acc cnt] [rec lt cnt n * acc cnt + cnt 1] 1 1]", "24")

    @test test("fac 4 def fac' \\[n]
        [drop swap eval [rec lt cnt n * acc cnt + cnt 1 def [acc cnt]] 1 1]", "24")

   # factorial with loop
   @test test("drop drop loop [lt 0 dup - swap 1 keep [*]] 4 1", "24")
   @test test("drop eval [rec lt 0 dup - swap 1 keep [*]] 4 1", "24")

    # Fibonacci numbers with simple recursion
    # simple recursion
    #@test test("fib 6
    # def fib' \\[n] [cond [ [+ fib - n 1 fib - n 2] n [lt n 2]]]", "8")

    @test test("fib 6
     def fib' \\[n] [eval if lt n 2 [n] [+ fib - n 1 fib - n 2]]" , "8")

    @test test("fib 6 def fib' [eval if gt 2 rot [] [+ fib - swap 1 swap fib - swap 2 dup] dup]", "8")

    # tail recursive
    @test test("fib 6 def fib' \\[n][
               loop 1 1 n def loop' \\[a b n][
               eval if rot [a] [loop b + a b - n 1]
             eq n 0 ] ]", "13")

    @test test("fib 6 def fib' \\[max] [
        drop drop fib-iter 1 0 1
        def fib-iter' [rec lt n max + n 1 j + i j def [n i j]] ]", "13")

     @test test("fib 6 def fib' \\[max] [
        drop drop fib-iter 1 0 1
        def fib-iter' \\[n i j] [rec lt n max + n 1 j + i j]] ]", "13")

    # Ackermann function
    #@test test("ack 3 4 def ack' \\[m n]
    #  [cond
    #    [ [ack - m 1 ack m - n 1]
    #      [ack - m 1 1]  [eq 0 n]
    #      [+ n 1]  [eq 0 m]] ]",  "125")

    @test test("ack 3 4 def ack' \\[m n]
     [eval if eq 0 m
        [+ n 1]
     [eval if eq 0 n
         [ack - m 1 1]
     [ack - m 1 ack m - n 1] ]] ", "125")

   end # testset

end

tests()

end #module


text:
"""
"""
struct DomainContribution <: GridapType
  dict::IdDict{Triangulation,AbstractArray}
end

DomainContribution() = DomainContribution(IdDict{Triangulation,AbstractArray}())

num_domains(a::DomainContribution) = length(a.dict)

get_domains(a::DomainContribution) = keys(a.dict)

function get_contribution(a::DomainContribution,trian::Triangulation)
  if haskey(a.dict,trian)
     return a.dict[trian]
  else
    @unreachable """\n
    There is not contribution associated with the given mesh in this DomainContribution object.
    """
  end
end

Base.getindex(a::DomainContribution,trian::Triangulation) = get_contribution(a,trian)

function add_contribution!(a::DomainContribution,trian::Triangulation,b::AbstractArray,op=+)

  S = eltype(b)
  if !(S<:AbstractMatrix || S<:AbstractVector || S<:Number)
    @unreachable """\n
    You are trying to add a contribution with eltype $(S).
    Only cell-wise matrices, vectors, or numbers are accepted.

    Make sure that you are defining the terms in your weak form correclty.
    """
  end

  if length(a.dict) > 0
    T = eltype(first(values(a.dict)))
    if T <: AbstractMatrix
      @assert S<:AbstractMatrix """\n
      You are trying to add a contribution with eltype $(S) to a DomainContribution that
      stores cell-wise matrices.

      Make sure that you are defining the terms in your weak form correclty.
      """
    elseif T <: AbstractVector
      @assert S<:AbstractVector """\n
      You are trying to add a contribution with eltype $(S) to a DomainContribution that
      stores cell-wise vectors.

      Make sure that you are defining the terms in your weak form correclty.
      """
    elseif T <: Number
      @assert S<:Number """\n
      You are trying to add a contribution with eltype $(S) to a DomainContribution that
      stores cell-wise numbers.

      Make sure that you are defining the terms in your weak form correclty.
      """
    end
  end

  if haskey(a.dict,trian)
    a.dict[trian] = lazy_map(Broadcasting(op),a.dict[trian],b)
  else
    if op == +
     a.dict[trian] = b
    else
     a.dict[trian] = lazy_map(Broadcasting(op),b)
    end
  end
  a
end

Base.sum(a::DomainContribution)= sum(map(sum,values(a.dict)))

Base.copy(a::DomainContribution) = DomainContribution(copy(a.dict))

function (+)(a::DomainContribution,b::DomainContribution)
  c = copy(a)
  for (trian,array) in b.dict
    add_contribution!(c,trian,array)
  end
  c
end

function (-)(a::DomainContribution,b::DomainContribution)
  c = copy(a)
  for (trian,array) in b.dict
    add_contribution!(c,trian,array,-)
  end
  c
end

function (*)(a::Number,b::DomainContribution)
  c = DomainContribution()
  for (trian,array_old) in b.dict
    s = size(get_cell_map(trian))
    array_new = lazy_map(Broadcasting(*),Fill(a,s),array_old)
    add_contribution!(c,trian,array_new)
  end
  c
end

(*)(a::DomainContribution,b::Number) = b*a

function get_array(a::DomainContribution)
  @assert num_domains(a) == 1 """\n
  Method get_array(a::DomainContribution) can be called only
  when the DomainContribution object involves just one domain.
  """
  a.dict[first(keys(a.dict))]
end

struct Measure <: GridapType
  quad::CellQuadrature
end

get_cell_points(a::Measure) = get_cell_points(a.quad)

Measure(args...) = Measure(CellQuadrature(args...))

function integrate(f,b::Measure)
  c = integrate(f,b.quad)
  cont = DomainContribution()
  add_contribution!(cont,b.quad.trian,c)
  cont
end

function (*)(a::Integrand,b::Measure)
  integrate(a.object,b)
end

(*)(b::Measure,a::Integrand) = a*b

ext: jl


text:
module KnetLayers

import Knet # seed!, atype, dir
using LinearAlgebra, CUDA
using AutoGrad: value, grad, cat1d, params, Param, @diff, @zerograd
using Knet.Ops20: relu, sigm, elu, invx, mat, bmm, logp, logsumexp, nll, bce, logistic, accuracy, zeroone, dropout, softmax, rnnforw, rnninit, rnnparam, RNN, BNMoments, conv4, deconv4, pool, unpool
using Knet.Train20: Data, minibatch, param, param0, gaussian, xavier, bilinear, SGD, Momentum, Nesterov, Adam, Adagrad, Adadelta, Rmsprop, update!, optimizers, train!, converge!, adam!, sgd!, nesterov!, rmsprop!, progress!, adam, converge, sgd, nesterov, rmsprop, progress, hyperband, goldensection
using Knet.KnetArrays: KnetArray, setseed, cpucopy, gpucopy, save, load, gc, JLDMODE
using Knet.LibKnet8: gpu # deprecated
import Knet.KnetArrays: _ser

export gpu, KnetArray,
       relu, sigm, elu, invx, mat, bmm, conv4, deconv4, pool, unpool,
       logp, logsumexp, nll, bce, logistic, accuracy, zeroone, dropout, softmax,
       Data, minibatch,
       param, param0, params,
       gaussian, xavier, bilinear, setseed,
       SGD, Momentum, Nesterov, Adam, Adagrad, Adadelta, Rmsprop, update!, optimizers,
       train!, converge!, adam!, sgd!, nesterov!, rmsprop!, progress!,
       adam, converge, sgd, nesterov, rmsprop, progress,
       hyperband, goldensection, cpucopy, gpucopy,
       value, grad, cat1d, Param, @diff, @zerograd

"""
    KnetLayers.dir(path...)
Construct a path relative to KnetLayers root.
# Example
```julia
julia> KnetLayers.dir("src")
"/Users/ekin/git/KnetLayers/src"
```
"""
dir(path...) = joinpath(dirname(@__DIR__),path...)
seed! = Knet.seed!
#Setters for atype
arrtype = Array{Float32}

"""
Used for setting default underlying array type for layer parameters.

    settype!(t::T) where T<:Type{KnetArray{V}} where V <: AbstractFloat = CUDA.functional() ? (global arrtype=t) : error("No GPU available")
    settype!(t::T) where T<:Type{Array{V}} where V <: AbstractFloat = (global arrtype=t)
    settype!(t::T) where T<:Union{Type{CuArray{V}},Type{KnetArray{V}}} where V <: AbstractFloat = CUDA.functional() ? (global arrtype=t) : error("No GPU available")
    settype!(t::Union{Type{CuArray},Type{KnetArray},Type{Array}}) = settype!(t{Float32})

# Example
```julia
julia> KnetLayers.settype!(KnetArray) # on a GPU machine
KnetArray{Float32}
```
"""
settype!(t::T) where T<:Type{KnetArray{V}} where V <: AbstractFloat = CUDA.functional() ? (global arrtype=t) : error("No GPU available")
settype!(t::T) where T<:Type{Array{V}} where V <: AbstractFloat = (global arrtype=t)
settype!(t::T) where T<:Union{Type{CuArray{V}},Type{KnetArray{V}}} where V <: AbstractFloat = CUDA.functional() ? (global arrtype=t) : error("No GPU available")
settype!(t::Union{Type{CuArray},Type{KnetArray},Type{Array}}) = settype!(t{Float32})

include("core.jl");
include("primitive.jl");   export Bias, Multiply, Embed, Linear, Dense, BatchNorm, Diagonal, LayerNorm
include("nonlinear.jl");   export NonAct, ReLU,Sigm,Tanh,LeakyReLU,ELU, Dropout, LogSoftMax, SoftMax, LogSumExp, GeLU
include("loss.jl");        export CrossEntropyLoss, BCELoss, LogisticLoss, SigmoidCrossEntropyLoss
include("cnn.jl");         export Pool,UnPool,DeConv,Conv
include("special.jl");     export MLP
include("rnn.jl");         export RNN,SRNN,LSTM,GRU,RNNOutput,PadRNNOutput,PadSequenceArray
include("chain.jl");       export Chain
include("attention.jl");   export MultiheadAttention
include("transformer.jl"); export Transformer, TransformerDecoder, PositionEmbedding, TransformerModel
include("datasets/Datasets.jl"); export Datasets
include("../data/IndexedDict.jl");
include("../data/one_hot.jl");

function __init__()
    global arrtype = CUDA.functional() ? KnetArray{Float32} : Array{Float32}
end

end # module


text:
using Test, NeuroAnalysis, BenchmarkTools, DataFrames, Plots, FileIO, FFTW

@testset "NeuroAnalysis" begin

    include("spiketest.jl")
    include("imagetest.jl")
    include("conditiontest.jl")
    include("regressiontest.jl")
    include("visualizationtest.jl")

## 1D grating
# @manipulate for f in 0:0.1:1, p in 0:0.001:1
#     plot([x->gratingf(x,f=f,phase=p),x->cas(x,f=f,phase=sin2cas(p))],-2,2,labels=["sin" "cas"])
# end
## 2D grating
# @manipulate for θ in 0:0.01:π, f in 0:0.05:1, sinp in 0:0.001:1, kx in -1:0.05:1,ky in -1:0.05:1, casp in 0:0.001:1
#     x=y=-2:0.05:2
#     sg = [gratingf(i,j,θ=θ,f=f,phase=sinp) for j in reverse(y),i in x]
#     scg = [cas(i,j;sin2cas(θ,f,sinp)...) for j in reverse(y),i in x]
#     cg = [cas(i,j,kx=kx,ky=ky,phase=casp) for j in reverse(y),i in x]
#     csg = [gratingf(i,j;cas2sin(kx,ky,casp)...) for j in reverse(y),i in x]
#
#     p = plot(layout=(2,3),leg=false,size=(600,400),clims=(-1,1),frame=:none,aspect_ratio=:equal,yflip=true)
#     heatmap!(p,subplot=1,sg,color=:grays,title="sin")
#     heatmap!(p,subplot=2,scg,color=:grays,title="sin2cas")
#     heatmap!(p,subplot=3,sg.-scg,color=:grays,title="sin - sin2cas")
#     heatmap!(p,subplot=4,cg,color=:grays,title="cas")
#     heatmap!(p,subplot=5,csg,color=:grays,title="cas2sin")
#     heatmap!(p,subplot=6,cg.-csg,color=:grays,title="cas - cas2sin")
#     p
# end

## functions
plot(vmf,-π,π)
plot(gvmf,-π,π)
plot(gratingf,-3,3)
plot(gaussianf,-3,3)
plot(gaborf,-3,3)
plot(dogf,-3,3)

funnames = ["grating","gaussian","gabor","dog"]
plot([gratingf,gaussianf,gaborf,dogf],-3,3,labels=permutedims(funnames),lw=[2 2 3 3],xlabel="y′")

x=y=-3:0.05:3;z=[]
push!(z,[gratingf(i,j,θ=0.25π,f=0.5) for j in y,i in x])
push!(z,[gaussianf(i,j,σ₁=0.8,σ₂=0.5,θ=0.25π) for j in y,i in x])
push!(z,[gaborf(i,j,σ₁=0.8,σ₂=0.5,θ=0.25π,f=0.5) for j in y,i in x])
push!(z,[dogf(i,j,σₑ₁=0.8,σₑ₂=0.5,σᵢ₁=1,σᵢ₂=0.7,θₑ=0.25π,θᵢ=0.25π) for j in y,i in x])

p = plot(layout=(2,2),legend=false,size=(600,600))
foreach(i->heatmap!(p,z[i],subplot=i,aspect_ratio=:equal,frame=:none,color=:coolwarm,clims=(-1,1),title=funnames[i]),1:4)
p

fs = 20
f = 2
t = 0:1/fs:10
x = cos.(2π*(f*t))
plot(t,x,xlabel="Second")
Fs = dft(x,fs,0,f)
FF = rfft(x)
plot(abs.(FF))

i = round(Int,f*length(x)/fs) + 1
@test FF[i] ≈ Fs[2]
@test FF[1] ≈ Fs[1]

end


text:
import MPI

# Abstract type for all possible bases that can be used in DFTK. Right now this is just
# one, but this type helps to resolve method ambiguities while avoiding an uninformative ::Any.
abstract type AbstractBasis{T <: Real} end

# There are two kinds of plane-wave basis sets used in DFTK.
# The k-dependent orbitals are discretized on spherical basis sets {G, 1/2 |k+G|^2 ≤ Ecut}.
# Potentials and densities are expressed on cubic basis sets large enough to contain
# products of orbitals. This also defines the real-space grid
# (as the dual of the cubic basis set).

"""
Discretization information for ``k``-point-dependent quantities such as orbitals.
More generally, a ``k``-point is a block of the Hamiltonian;
eg collinear spin is treated by doubling the number of kpoints.
"""
struct Kpoint{T <: Real}
    spin::Int                     # Spin component can be 1 or 2 as index into what is
                                  # returned by the `spin_components` function
    coordinate::Vec3{T}           # Fractional coordinate of k-point
    mapping::Vector{Int}          # Index of G_vectors[i] on the FFT grid:
                                  # G_vectors(basis)[kpt.mapping[i]] == G_vectors(basis, kpt)[i]
    mapping_inv::Dict{Int, Int}   # Inverse of `mapping`:
                                  # G_vectors(basis)[i] == G_vectors(basis, kpt)[mapping_inv[i]]
    G_vectors::Vector{Vec3{Int}}  # Wave vectors in integer coordinates:
                                  # ({G, 1/2 |k+G|^2 ≤ Ecut})
end

@doc raw"""
A plane-wave discretized `Model`.
Normalization conventions:
- Things that are expressed in the G basis are normalized so that if ``x`` is the vector,
  then the actual function is ``\sum_G x_G e_G`` with
  ``e_G(x) = e^{iG x} / \sqrt(\Omega)``, where ``\Omega`` is the unit cell volume.
  This is so that, eg ``norm(ψ) = 1`` gives the correct normalization.
  This also holds for the density and the potentials.
- Quantities expressed on the real-space grid are in actual values.

`G_to_r` and `r_to_G` convert between these representations.
"""
struct PlaneWaveBasis{T} <: AbstractBasis{T}
    model::Model{T}

    ## Global grid information
    # fft_size defines both the G basis on which densities and
    # potentials are expanded, and the real-space grid
    fft_size::Tuple{Int, Int, Int}
    # factor for integrals in real space: sum(ρ) * dvol ~ ∫ρ
    dvol::T  # = model.unit_cell_volume ./ prod(fft_size)
    # Information used to construct the k-point-specific basis
    # (not used directly after that)
    Ecut::T  # The basis set is defined by {e_{G}, 1/2|k+G|^2 ≤ Ecut}
    variational::Bool  # Is the k-point specific basis variationally consistent with
    #                    the basis used for the density / potential?

    ## Plans for forward and backward FFT
    # All these plans are *completely unnormalized* (eg FFT * BFFT != I)
    # The normalizations are performed in G_to_r/r_to_G according to
    # the DFTK conventions (see above)
    opFFT   # out-of-place FFT plan
    ipFFT   # in-place FFT plan
    opBFFT  # inverse plans (unnormalized plan; backward in FFTW terminology)
    ipBFFT
    r_to_G_normalization::T  # r_to_G = r_to_G_normalization * FFT
    G_to_r_normalization::T  # G_to_r = G_to_r_normalization * BFFT

    ## MPI-local information of the kpoints this processor treats
    # Irreducible kpoints. In the case of collinear spin,
    # this lists all the spin up, then all the spin down
    kpoints::Vector{Kpoint{T}}
    # BZ integration weights, summing up to model.n_spin_components
    kweights::Vector{T}

    ## (MPI-global) information on the k-point grid
    ## These fields are not actually used in computation, but can be used to reconstruct a basis
    # Monkhorst-Pack grid used to generate the k-points, or nothing for custom k-points
    kgrid::Union{Nothing,Vec3{Int}}
    kshift::Union{Nothing,Vec3{T}}
    # full list of (non spin doubled) k-point coordinates in the irreducible BZ
    kcoords_global::Vector{Vec3{T}}
    kweights_global::Vector{T}

    ## Setup for MPI-distributed processing over k-points
    comm_kpts::MPI.Comm           # communicator for the kpoints distribution
    krange_thisproc::Vector{Int}  # indices of kpoints treated explicitly by this
    #                               processor in the global kcoords array
    krange_allprocs::Vector{Vector{Int}}  # indices of kpoints treated by the
    #                                       respective rank in comm_kpts

    ## Symmetry operations that leave the discretized model (k and r grids) invariant.
    # Subset of model.symmetries.
    symmetries::Vector{SymOp}
    # Whether the symmetry operations leave the rgrid invariant
    # If this is true, the symmetries are a property of the complete discretized model.
    # Therefore, all quantities should be symmetric to machine precision
    symmetries_respect_rgrid::Bool

    ## Instantiated terms (<: Term). See Hamiltonian for high-level usage
    terms::Vector{Any}
end

# prevent broadcast
import Base.Broadcast.broadcastable
Base.Broadcast.broadcastable(basis::PlaneWaveBasis) = Ref(basis)

Base.eltype(::PlaneWaveBasis{T}) where {T} = T

@timing function build_kpoints(model::Model{T}, fft_size, kcoords, Ecut;
                               variational=true) where T
    kpoints_per_spin = [Kpoint[] for _ in 1:model.n_spin_components]
    for k in kcoords
        k = Vec3{T}(k)  # rationals are sloooow
        mapping = Int[]
        Gvecs_k = Vec3{Int}[]
        # provide a rough hint so that the arrays don't have to be resized so much
        n_guess = div(prod(fft_size), 8)
        sizehint!(mapping, n_guess)
        sizehint!(Gvecs_k, n_guess)
        for (i, G) in enumerate(G_vectors(fft_size))
            if !variational || sum(abs2, model.recip_lattice * (G + k)) / 2 ≤ Ecut
                push!(mapping, i)
                push!(Gvecs_k, G)
            end
        end
        mapping_inv = Dict(ifull => iball for (iball, ifull) in enumerate(mapping))
        for iσ = 1:model.n_spin_components
            push!(kpoints_per_spin[iσ],
                  Kpoint(iσ, k, mapping, mapping_inv, Gvecs_k))
        end
    end

    vcat(kpoints_per_spin...)  # put all spin up first, then all spin down
end
function build_kpoints(basis::PlaneWaveBasis, kcoords)
    build_kpoints(basis.model, basis.fft_size, kcoords, basis.Ecut;
                  variational=basis.variational)
end

# Lowest-level constructor, should not be called directly.
# All given parameters must be the same on all processors
# and are stored in PlaneWaveBasis for easy reconstruction.
function PlaneWaveBasis(model::Model{T}, Ecut::Number, fft_size, variational,
                        kcoords, kweights, kgrid, kshift,
                        symmetries_respect_rgrid, comm_kpts) where {T <: Real}
    # Validate fft_size
    if variational
        max_E = sum(abs2, model.recip_lattice * floor.(Int, Vec3(fft_size) ./ 2)) / 2
        Ecut > max_E && @warn(
            "For a variational method, Ecut should be less than the maximal kinetic " *
            "energy the grid supports ($max_E)"
        )
    end
    if !(all(fft_size .== next_working_fft_size(T, fft_size)))
        @show fft_size next_working_fft_size(T, fft_size)
        error("Selected fft_size will not work for the buggy generic " *
              "FFT routines; use next_working_fft_size")
    end
    fft_size = Tuple{Int, Int, Int}(fft_size)  # explicit conversion in case passed as array

    # filter out the symmetries that don't preserve the real-space grid
    symmetries = model.symmetries
    if symmetries_respect_rgrid
        symmetries = symmetries_preserving_rgrid(symmetries, fft_size)
    end

    # build or validate the kgrid, and get symmetries preserving the kgrid
    if isnothing(kcoords)
        # MP grid based on kgrid/kshift
        @assert !isnothing(kgrid)
        @assert !isnothing(kshift)
        @assert isnothing(kweights)
        kcoords, kweights, symmetries = bzmesh_ir_wedge(kgrid, symmetries; kshift)
    else
        # Manual kpoint set based on kcoords/kweights
        @assert length(kcoords) == length(kweights)
        all_kcoords = unfold_kcoords(kcoords, symmetries)
        symmetries = symmetries_preserving_kgrid(symmetries, all_kcoords)
    end

    # Init MPI, and store MPI-global values for reference
    MPI.Init()
    kcoords_global  = kcoords
    kweights_global = kweights

    # Setup FFT plans
    (ipFFT, opFFT, ipBFFT, opBFFT) = build_fft_plans(T, fft_size)

    # Normalization constants
    # r_to_G = r_to_G_normalization * FFT
    # The convention we want is
    # ψ(r) = sum_G c_G e^iGr / sqrt(Ω)
    # so that the G_to_r has to normalized by 1/sqrt(Ω).
    # The other constant is chosen because FFT * BFFT = N
    G_to_r_normalization = 1/sqrt(model.unit_cell_volume)
    r_to_G_normalization = sqrt(model.unit_cell_volume) / length(ipFFT)

    # Compute k-point information and spread them across processors
    # Right now we split only the kcoords: both spin channels have to be handled
    # by the same process
    n_kpt   = length(kcoords_global)
    n_procs = mpi_nprocs(comm_kpts)
    if n_procs > n_kpt
        # XXX Supporting more processors than kpoints would require
        # fixing a bunch of "reducing over empty collections" errors
        # In the unit tests it is really annoying that this fails so we hack around it, but
        # generally it leads to duplicated work that is not in the users interest.
        if parse(Bool, get(ENV, "CI", "false"))
            comm_kpts = MPI.COMM_SELF
            krange_thisproc = 1:n_kpt
            krange_allprocs = fill(1:n_kpt, n_procs)
        else
            error("No point in trying to parallelize $n_kpt kpoints over $n_procs " *
                  "processes; reduce the number of MPI processes.")
        end
    else
        # get the slice of 1:n_kpt to be handled by this process
        # Note: MPI ranks are 0-based
        krange_allprocs = split_evenly(1:n_kpt, n_procs)
        krange_thisproc = krange_allprocs[1 + MPI.Comm_rank(comm_kpts)]
        @assert mpi_sum(length(krange_thisproc), comm_kpts) == n_kpt
        @assert !isempty(krange_thisproc)
    end
    kweights_thisproc = kweights[krange_thisproc]

    # Setup k-point basis sets
    !variational && @warn(
        "Non-variational calculations are experimental. " *
        "Not all features of DFTK may be supported or work as intended."
    )
    kpoints = build_kpoints(model, fft_size, kcoords_global[krange_thisproc], Ecut; variational)
    # kpoints is now possibly twice the size of krange. Make things consistent
    if model.n_spin_components == 2
        krange_thisproc   = vcat(krange_thisproc, n_kpt .+ krange_thisproc)
        krange_allprocs   = [vcat(range, n_kpt .+ range) for range in krange_allprocs]
        kweights_thisproc = vcat(kweights_thisproc, kweights_thisproc)
    end
    @assert mpi_sum(sum(kweights_thisproc), comm_kpts) ≈ model.n_spin_components
    @assert length(kpoints) == length(kweights_thisproc)

    dvol  = model.unit_cell_volume ./ prod(fft_size)
    terms = Vector{Any}(undef, length(model.term_types))  # Dummy terms array, filled below
    basis = PlaneWaveBasis{T}(
        model, fft_size, dvol,
        Ecut, variational,
        opFFT, ipFFT, opBFFT, ipBFFT,
        r_to_G_normalization, G_to_r_normalization,
        kpoints, kweights_thisproc, kgrid, kshift,
        kcoords_global, kweights_global, comm_kpts, krange_thisproc, krange_allprocs,
        symmetries, symmetries_respect_rgrid, terms)

    # Instantiate the terms with the basis
    for (it, t) in enumerate(model.term_types)
        term_name = string(nameof(typeof(t)))
        @timing "Instantiation $term_name" basis.terms[it] = t(basis)
    end
    basis
end

# This is an intermediate-level constructor, which allows for the
# custom specification of k points and G grids.
# For regular usage, the higher-level one below should be preferred
@timing function PlaneWaveBasis(model::Model{T}, Ecut::Number,
                                kcoords ::Union{Nothing, AbstractVector},
                                kweights::Union{Nothing, AbstractVector};
                                variational=true, fft_size=nothing,
                                kgrid=nothing, kshift=nothing,
                                symmetries_respect_rgrid=isnothing(fft_size),
                                comm_kpts=MPI.COMM_WORLD) where {T <: Real}
    if isnothing(fft_size)
        @assert variational
        if symmetries_respect_rgrid
            # ensure that the FFT grid is compatible with the "reasonable" symmetries
            # (those with fractional translations with denominators 2, 3, 4, 6,
            #  this set being more or less arbitrary) by forcing the FFT size to be
            # a multiple of the denominators.
            # See https://github.com/JuliaMolSim/DFTK.jl/pull/642 for discussion
            denominators = [denominator(rationalize(sym.w[i]; tol=SYMMETRY_TOLERANCE))
                            for sym in model.symmetries for i = 1:3]
            factors = intersect((2, 3, 4, 6), denominators)
        else
            factors = (1,)
        end
        fft_size = compute_fft_size(model, Ecut, kcoords; factors=factors)
    end
    PlaneWaveBasis(model, Ecut, fft_size, variational, kcoords, kweights,
                   kgrid, kshift, symmetries_respect_rgrid, comm_kpts)
end

@doc raw"""
Creates a `PlaneWaveBasis` using the kinetic energy cutoff `Ecut` and a Monkhorst-Pack
``k``-point grid. The MP grid can either be specified directly with `kgrid` providing the
number of points in each dimension and `kshift` the shift (0 or 1/2 in each direction).
If not specified a grid is generated using `kgrid_from_minimal_spacing` with
a minimal spacing of `2π * 0.022` per Bohr.
"""
function PlaneWaveBasis(model::Model; Ecut,
                        kgrid=kgrid_from_minimal_spacing(model, 2π * 0.022),
                        kshift=zeros(3),
                        kwargs...)
    PlaneWaveBasis(model, austrip(Ecut), nothing, nothing;
                   kgrid, kshift, kwargs...)
end

"""
Creates a new basis identical to `basis`, but with a custom set of kpoints
"""
@timing function PlaneWaveBasis(basis::PlaneWaveBasis, kcoords::AbstractVector,
                                kweights::AbstractVector)
    kgrid = kshift = nothing
    PlaneWaveBasis(basis.model, basis.Ecut,
                   basis.fft_size, basis.variational,
                   kcoords, kweights, kgrid, kshift,
                   basis.symmetries_respect_rgrid, basis.comm_kpts)
end

"""
    G_vectors(fft_size::Tuple)

The wave vectors `G` in reduced (integer) coordinates for a cubic basis set
of given sizes.
"""
function G_vectors(fft_size::Union{Tuple,AbstractVector})
    # Note that a collect(G_vectors_generator(fft_size)) is 100-fold slower
    # than this implementation, hence the code duplication.
    start = .- cld.(fft_size .- 1, 2)
    stop  = fld.(fft_size .- 1, 2)
    axes  = [[collect(0:stop[i]); collect(start[i]:-1)] for i in 1:3]
    [Vec3{Int}(i, j, k) for i in axes[1], j in axes[2], k in axes[3]]
end
function G_vectors_generator(fft_size::Union{Tuple,AbstractVector})
    # The generator version is used mainly in symmetry.jl for lowpass_for_symmetry! and
    # accumulate_over_symmetries!, which are 100-fold slower with G_vector(fft_size).
    start = .- cld.(fft_size .- 1, 2)
    stop  = fld.(fft_size .- 1, 2)
    axes = [[collect(0:stop[i]); collect(start[i]:-1)] for i in 1:3]
    (Vec3{Int}(i, j, k) for i in axes[1], j in axes[2], k in axes[3])
end

@doc raw"""
    G_vectors(basis::PlaneWaveBasis)
    G_vectors(basis::PlaneWaveBasis, kpt::Kpoint)

The list of wave vectors ``G`` in reduced (integer) coordinates of a `basis`
or a ``k``-point `kpt`.
"""
G_vectors(basis::PlaneWaveBasis) = G_vectors(basis.fft_size)
G_vectors(::PlaneWaveBasis, kpt::Kpoint) = kpt.G_vectors

@doc raw"""
    G_vectors_cart(basis::PlaneWaveBasis)
    G_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)

The list of ``G`` vectors of a given `basis` or `kpt`, in cartesian coordinates.
"""
G_vectors_cart(basis::PlaneWaveBasis) = recip_vector_red_to_cart.(basis.model, G_vectors(basis))
function G_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)
    recip_vector_red_to_cart.(basis.model, G_vectors(basis, kpt))
end

@doc raw"""
    Gplusk_vectors(basis::PlaneWaveBasis, kpt::Kpoint)

The list of ``G + k`` vectors, in reduced coordinates.
"""
function Gplusk_vectors(basis::PlaneWaveBasis, kpt::Kpoint)
    map(G -> G + kpt.coordinate, G_vectors(basis, kpt))
end

@doc raw"""
    Gplusk_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)

The list of ``G + k`` vectors, in cartesian coordinates.
"""
function Gplusk_vectors_cart(basis::PlaneWaveBasis, kpt::Kpoint)
    recip_vector_red_to_cart.(basis.model, Gplusk_vectors(basis, kpt))
end

@doc raw"""
    r_vectors(basis::PlaneWaveBasis)

The list of ``r`` vectors, in reduced coordinates. By convention, this is in [0,1)^3.
"""
function r_vectors(basis::PlaneWaveBasis{T}) where T
    N1, N2, N3 = basis.fft_size
    [Vec3{T}(T(i-1) / N1, T(j-1) / N2, T(k-1) / N3) for i = 1:N1, j = 1:N2, k = 1:N3]
end

@doc raw"""
    r_vectors_cart(basis::PlaneWaveBasis)

The list of ``r`` vectors, in cartesian coordinates.
"""
r_vectors_cart(basis::PlaneWaveBasis) = vector_red_to_cart.(basis.model, r_vectors(basis))

"""
Return the index tuple `I` such that `G_vectors(basis)[I] == G`
or the index `i` such that `G_vectors(basis, kpoint)[i] == G`.
Returns nothing if outside the range of valid wave vectors.
"""
function index_G_vectors(basis::PlaneWaveBasis, G::AbstractVector{T}) where {T <: Integer}
    start = .- cld.(basis.fft_size .- 1, 2)
    stop  = fld.(basis.fft_size .- 1, 2)
    lengths = stop .- start .+ 1

    function mapaxis(lengthi, Gi)
        Gi >= 0 && return 1 + Gi
        return 1 + lengthi + Gi
    end
    if all(start .<= G .<= stop)
        CartesianIndex(Tuple(mapaxis.(lengths, G)))
    else
        nothing  # Outside range of valid indices
    end
end

function index_G_vectors(basis::PlaneWaveBasis, kpoint::Kpoint,
                         G::AbstractVector{T}) where {T <: Integer}
    fft_size = basis.fft_size
    idx = index_G_vectors(basis, G)
    isnothing(idx) && return nothing
    idx_linear = LinearIndices(fft_size)[idx]
    get(kpoint.mapping_inv, idx_linear, nothing)
end

"""
Return the index range of ``k``-points that have a particular spin component.
"""
function krange_spin(basis::PlaneWaveBasis, spin::Integer)
    n_spin = basis.model.n_spin_components
    @assert 1 ≤ spin ≤ n_spin
    spinlength = div(length(basis.kpoints), n_spin)
    (1 + (spin - 1) * spinlength):(spin * spinlength)
end

"""
Sum an array over kpoints, taking weights into account
"""
function weighted_ksum(basis::PlaneWaveBasis, array)
    res = sum(basis.kweights .* array)
    mpi_sum(res, basis.comm_kpts)
end

"""
Gather the distributed ``k``-point data on the master process and return
it as a `PlaneWaveBasis`. On the other (non-master) processes `nothing` is returned.
The returned object should not be used for computations and only to extract data
for post-processing and serialisation to disk.
"""
function gather_kpts(basis::PlaneWaveBasis)
    # No need to allocate and setup a new basis object
    mpi_nprocs(basis.comm_kpts) == 1 && return basis

    # Gather k-point info on master
    kcoords  = getproperty.(basis.kpoints, :coordinate)
    kcoords  = gather_kpts(kcoords, basis)
    kweights = gather_kpts(basis.kweights, basis)

    # Number of distinct k-point coordinates is number of k-points with spin 1
    n_spinup_thisproc = count(kpt.spin == 1 for kpt in basis.kpoints)
    n_kcoords = mpi_sum(n_spinup_thisproc, basis.comm_kpts)

    if isnothing(kcoords)  # i.e. master process
        nothing
    else
        PlaneWaveBasis(basis.model,
                       basis.Ecut,
                       kcoords[1:n_kcoords],
                       kweights[1:n_kcoords];
                       basis.variational,
                       basis.kgrid,
                       basis.kshift,
                       basis.symmetries_respect_rgrid,
                       comm_kpts=MPI.COMM_SELF,
                      )
    end
end

"""
Gather the distributed data of a quantity depending on `k`-Points on the master process
and return it. On the other (non-master) processes `nothing` is returned.
"""
function gather_kpts(data::AbstractArray, basis::PlaneWaveBasis)
    master = tag = 0
    n_kpts = sum(length, basis.krange_allprocs)

    if MPI.Comm_rank(basis.comm_kpts) == master
        allk_data = similar(data, n_kpts)
        allk_data[basis.krange_allprocs[1]] = data
        for rank in 1:mpi_nprocs(basis.comm_kpts) - 1  # Note: MPI ranks are 0-based
            rk_data, status = MPI.recv(rank, tag, basis.comm_kpts)
            @assert MPI.Get_error(status) == 0  # all went well
            allk_data[basis.krange_allprocs[rank + 1]] = rk_data
        end
        allk_data
    else
        MPI.send(data, master, tag, basis.comm_kpts)
        nothing
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct InboundNatRulePropertiesFormat <: SwaggerModel
    frontendIPConfiguration::Any # spec type: Union{ Nothing, SubResource } # spec name: frontendIPConfiguration
    backendIPConfiguration::Any # spec type: Union{ Nothing, NetworkInterfaceIPConfiguration } # spec name: backendIPConfiguration
    protocol::Any # spec type: Union{ Nothing, TransportProtocol } # spec name: protocol
    frontendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: frontendPort
    backendPort::Any # spec type: Union{ Nothing, Int32 } # spec name: backendPort
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    enableFloatingIP::Any # spec type: Union{ Nothing, Bool } # spec name: enableFloatingIP
    enableTcpReset::Any # spec type: Union{ Nothing, Bool } # spec name: enableTcpReset
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function InboundNatRulePropertiesFormat(;frontendIPConfiguration=nothing, backendIPConfiguration=nothing, protocol=nothing, frontendPort=nothing, backendPort=nothing, idleTimeoutInMinutes=nothing, enableFloatingIP=nothing, enableTcpReset=nothing, provisioningState=nothing)
        o = new()
        validate_property(InboundNatRulePropertiesFormat, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        setfield!(o, Symbol("frontendIPConfiguration"), frontendIPConfiguration)
        validate_property(InboundNatRulePropertiesFormat, Symbol("backendIPConfiguration"), backendIPConfiguration)
        setfield!(o, Symbol("backendIPConfiguration"), backendIPConfiguration)
        validate_property(InboundNatRulePropertiesFormat, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(InboundNatRulePropertiesFormat, Symbol("frontendPort"), frontendPort)
        setfield!(o, Symbol("frontendPort"), frontendPort)
        validate_property(InboundNatRulePropertiesFormat, Symbol("backendPort"), backendPort)
        setfield!(o, Symbol("backendPort"), backendPort)
        validate_property(InboundNatRulePropertiesFormat, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(InboundNatRulePropertiesFormat, Symbol("enableFloatingIP"), enableFloatingIP)
        setfield!(o, Symbol("enableFloatingIP"), enableFloatingIP)
        validate_property(InboundNatRulePropertiesFormat, Symbol("enableTcpReset"), enableTcpReset)
        setfield!(o, Symbol("enableTcpReset"), enableTcpReset)
        validate_property(InboundNatRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type InboundNatRulePropertiesFormat

const _property_map_InboundNatRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("frontendIPConfiguration")=>Symbol("frontendIPConfiguration"), Symbol("backendIPConfiguration")=>Symbol("backendIPConfiguration"), Symbol("protocol")=>Symbol("protocol"), Symbol("frontendPort")=>Symbol("frontendPort"), Symbol("backendPort")=>Symbol("backendPort"), Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("enableFloatingIP")=>Symbol("enableFloatingIP"), Symbol("enableTcpReset")=>Symbol("enableTcpReset"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_InboundNatRulePropertiesFormat = Dict{Symbol,String}(Symbol("frontendIPConfiguration")=>"SubResource", Symbol("backendIPConfiguration")=>"NetworkInterfaceIPConfiguration", Symbol("protocol")=>"TransportProtocol", Symbol("frontendPort")=>"Int32", Symbol("backendPort")=>"Int32", Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("enableFloatingIP")=>"Bool", Symbol("enableTcpReset")=>"Bool", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ InboundNatRulePropertiesFormat }) = collect(keys(_property_map_InboundNatRulePropertiesFormat))
Swagger.property_type(::Type{ InboundNatRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_InboundNatRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ InboundNatRulePropertiesFormat }, property_name::Symbol) =  _property_map_InboundNatRulePropertiesFormat[property_name]

function check_required(o::InboundNatRulePropertiesFormat)
    true
end

function validate_property(::Type{ InboundNatRulePropertiesFormat }, name::Symbol, val)
end


text:
function solve()
    xs = [parse(Int, x) for x in split(readline())]
    15 - sum(xs)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
# This file was generated, do not modify it. # hide
r_mpi = range(mdl, :(tree.max_depth), lower=1, upper=10)
r_msl = range(mdl, :(tree.min_samples_leaf), lower=1, upper=50)

tm = TunedModel(model=mdl, ranges=[r_mpi, r_msl], tuning=Grid(resolution=8),
                resampling=CV(nfolds=5, rng=112),
                operation=predict_mode, measure=misclassification_rate)
mtm = machine(tm, X, y)
fit!(mtm, rows=train)

ypred = predict_mode(mtm, rows=test)
misclassification_rate(ypred, y[test])


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ManagedByTenant <: SwaggerModel
    tenantId::Any # spec type: Union{ Nothing, String } # spec name: tenantId

    function ManagedByTenant(;tenantId=nothing)
        o = new()
        validate_property(ManagedByTenant, Symbol("tenantId"), tenantId)
        setfield!(o, Symbol("tenantId"), tenantId)
        o
    end
end # type ManagedByTenant

const _property_map_ManagedByTenant = Dict{Symbol,Symbol}(Symbol("tenantId")=>Symbol("tenantId"))
const _property_types_ManagedByTenant = Dict{Symbol,String}(Symbol("tenantId")=>"String")
Base.propertynames(::Type{ ManagedByTenant }) = collect(keys(_property_map_ManagedByTenant))
Swagger.property_type(::Type{ ManagedByTenant }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagedByTenant[name]))}
Swagger.field_name(::Type{ ManagedByTenant }, property_name::Symbol) =  _property_map_ManagedByTenant[property_name]

function check_required(o::ManagedByTenant)
    true
end

function validate_property(::Type{ ManagedByTenant }, name::Symbol, val)
end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------


text:
include("linear_approximator.jl")
include("dyna_agents/dyna_agents.jl")
include("monte_carlo_learner.jl")
include("td_learner.jl")
include("gradient_bandit_learner.jl")
include("double_learner.jl")
include("policy_iteration.jl")
include("tabular_policy.jl")
include("off_policy.jl")


text:
using QuantumInfo, RandomQuantum
using Base.Test

trnorm(A) = sum(svdvals(A))

#@testloop "Randomized GinUE tests"
for _ in 1:10
    rows = rand(2:10)
    cols = rand(1:10)
    X = rand(GinibreEnsemble(rows,cols))
    @test var(vec(real(X))) > 0
    @test var(vec(imag(X))) > 0
end

#@testloop "Randomized FubiniStudyEnsemble tests"
for _ in 1:10
    d = rand(2:10)
    ψ = rand(FubiniStudyPureState(d))
    @test var(real(ψ)) > 0
    @test var(imag(ψ)) > 0
    @test_approx_eq norm(ψ) 1.0
end

#@testloop "Randomized FSOE tests"
for _ in 1:10
    d = rand(2:10)
    db = rand(2:10)
    ρ = rand(FubiniStudyMixedState(d,db))
    @test_approx_eq trnorm(ρ) 1
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end

#@testloop "Randomized HSOE tests"
for _ in 1:10
    d = rand(2:10)
    ρ = rand(HilbertSchmidtMixedState(d))
    @test_approx_eq trnorm(ρ) 1
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end

#BuresMixedState
for _ in 1:10
    d = rand(2:10)
    ρ = rand(BuresMixedState(d))
    @test_approx_eq trnorm(ρ) 1
    @test ishermitian(ρ)
    @test ispossemidef(ρ)
    @test var(vec(real(ρ))) > 0
    @test var(vec(imag(ρ))) > 0
end

#ClosedHaarEnsemble
for _ in 1:10
    d = rand(2:10)
    U = rand(ClosedHaarEnsemble(d))
    @test_approx_eq_eps norm(U*U'-eye(d)) 0.0 1e-12
    @test_approx_eq_eps norm(U'*U-eye(d)) 0.0 1e-12
    @test_approx_eq_eps abs(eigvals(U)) ones(d) 1e-12
end

#OpenHaarEnsemble
for _ in 1:10
    d = 2 # rand(2:10)
    db = rand(2:10)
    E = rand(OpenHaarEnsemble(d,db))
    @test istp(E,tol=1e-14)
    @test iscp(E,tol=1e-14)
    @test ischannel(E,tol=1e-14)
end

#ClosedEvolution
for _ in 1:10
    α = rand()
    d = rand(2:10)
    U = rand(RandomClosedEvolution(d,α))
    @test_approx_eq_eps norm(U*U'-eye(d)) 0.0 1e-12
    @test_approx_eq_eps norm(U'*U-eye(d)) 0.0 1e-12
    @test_approx_eq_eps abs(eigvals(U)) ones(d) 1e-12
end

#OpenEvolution
for _ in 1:10
    α = rand()
    d = rand(2:10)
    db = rand(2:10)
    E = rand(RandomOpenEvolution(d,db,α))
    @test ischannel(E,tol=1e-12)
end


text:
module drbasisj

  export drbasis,holajulia

  include("drbasis-function.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LIBMINCd_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LIBMINCd")
JLLWrappers.@generate_main_file("LIBMINCd", UUID("7b82580a-3892-5ced-a724-ff39eef97d7e"))
end  # module LIBMINCd_jll


text:
export simulate, TrajectoryData, visualize

import Plots
using Makie
using AbstractPlotting.MakieLayout
using StaticArrays
using Unrolled
using Base.Iterators: drop
using ColorSchemes

"""
You can use a component name as the index to retrive the corresponding component.
The result is a matrix indexed by (time, agent).
"""
@kwdef struct TrajectoryData
    times::Vector{𝕋}
    "values are indexed by (time, agent, component)"
    values::Array{ℝ,3}
    components::Vector{String}
    function TrajectoryData(times::AbstractVector{𝕋}, num_agents::ℕ, components::Vector{String})
        values = Array{ℝ}(undef, length(times), num_agents, length(components))
        new(times, values, components)
    end
end

function Base.getindex(data::TrajectoryData, comp::String)
    c_id = indexin([comp], data.components)[1]
    @assert c_id !== nothing "component '$comp' not found."
    view(data.values, :,:, (c_id::Int))
end

"""
A generic result visualization function that draws multiple line plots,
one for each component of the state vector.
"""
function visualize(
    result::TrajectoryData; delta_t = nothing, loss = nothing,
)::Plots.Plot
    if delta_t === nothing
        times = result.times
        xlabel = "time step"
    else
        times = result.times .* delta_t
        xlabel = "time (s)"
    end
    ## currently assume all agents have the same type of states and observations
    N = size(result.values, 2)
    labels = reshape(["agent $i" for i in 1:N], (1, :))
    n = length(result.components)
    colors = [get(ColorSchemes.rainbow, j / N) for j in 1:N]
    ps = map(enumerate(result.components)) do (c_id, comp)
        ys = result.values[:,:,c_id]
        Plots.plot(times, ys, label=labels, title=comp)
    end

    if loss !== nothing
        let ts = collect(AxisArrays.axes(loss)[1])
            (delta_t !== nothing) && (ts = ts .* delta_t)
            ps2 = map([(:u, "Action"), (:x, "State")]) do (label, name)
                ys = collect(loss[comp=label])
                Plots.plot(ts, ys, label=name, title=name)
            end
            append!(ps, ps2)
        end
    end

    Plots.plotlyjs()
    Plots.plot(ps..., layout=(length(ps), 1), size=(500, length(ps) * 200), legend=false)
end

function shorten_queue(q, len)
    @assert(len ≤ length(q),
        "Actual communication delay ($len) should ≤ the modeled delay ($(length(q))).")
    to_drop = length(q) - len
    FixedQueue(collect(drop(collect(q), to_drop)))
end

"""
Simulate multiple distributed agents with delayed communication.

Returns `logs`.

# Arguments
- `callback`: a function of the form `callback(xs,zs,us,loss,t) -> nothing` that runs at
every time step. `loss` will be `nothing` if loss_model is `nothing`.
"""
function simulate(
    world_dynamics::WorldDynamics{N},
    delay_model::DelayModel,
    framework::ControllerFramework{X,Z,U,Msg,Log},
    init_status::Each{Tuple{X,Z,U}},
    callback::Function,
    (t0, tf)::Tuple{𝕋, 𝕋};
    loss_model::Union{RegretLossModel, Nothing} = nothing,
    xs_observer = (xs, t) -> xs,
    zs_observer = (zs, t) -> zs,
) where {X,Z,U,Msg,Log,N}
    @assert tf ≥ t0
    @unpack Tx, Tu, Tc = short_delay_names(delay_model)
    @assert Tc > 0

    if loss_model !== nothing
        # the *ideal* state and observation at the current time step
        xs_idl, _, us_idl = @unzip(MVector{N}(init_status), MVector{N}{Tuple{X,Z,U}})
        @unpack central, x_weights, u_weights = loss_model
        modeled_dynamics = loss_model.world_model
        s_c = init_state(central, t0)
        loss_cache = AxisArray(-ones(ℝ, N, 2), id=1:N, comp=[:u, :x])
    else
        loss_cache = nothing
    end

    # the *actual* state, observation, and actuation at the current time step
    xs, zs, us = @unzip(MVector{N}(init_status), MVector{N}{Tuple{X,Z,U}})

    (controllers, msg_qs) = make_controllers(framework, init_status, t0)
    # the head of these queues corresponding to values that are currently taking effect
    msg_qs = shorten_queue.(msg_qs, Tc)
    state_qs = SVector{N}(constant_queue.(xs, Tx))
    obs_qs = SVector{N}(constant_queue.(zs, Tx))
    act_qs = SVector{N}(constant_queue.(us, Tu))

    # we store the newest messages/actions into these caches and wait until
    # the next control/actuation step to push them into the queues to take effect.
    msg_cache = SizedMatrix{N,N,Msg}(hcat(first.(msg_qs)...)) # [sender, receiver]
    act_cache = MVector{N, U}(us)
    for t in t0:tf
        xs_read = pushpop!.(state_qs, xs_observer(xs, t))
        zs_read = pushpop!.(obs_qs, zs_observer(zs, t))
        # run controllers
        for i in SOneTo(N)
            x, z = xs_read[i], zs_read[i]
            ms = first(msg_qs[i])
            u, ms′ = control!(controllers[i], x, z, ms)
            u = limit_control(world_dynamics.dynamics[i], u, x, t)
            msg_cache[i, :] = ms′
            act_cache[i] = u
        end
        us .= pushpop!.(act_qs, act_cache)
        for i in SOneTo(N)
            pushpop!(msg_qs[i], convert(Vector{Msg}, msg_cache[:, i]))
        end

        # record loss
        if loss_model !== nothing
            us_idl .= control_all(central, s_c, xs_idl, zs, t, SOneTo(N))
            us_idl .= limit_control.(modeled_dynamics.dynamics, us_idl, xs_idl, t)
            for i in 1:N
                x_loss = sum((xs[i] .- xs_idl[i]).^2 .* x_weights[i])
                u_loss = sum((collect(us[i]) .- us_idl[i]).^2 .* u_weights[i])
                loss_cache[id = i, comp=:u] = u_loss
                loss_cache[id = i, comp=:x] = x_loss
            end
        end

        # record results
        callback(xs, zs, us, loss_cache, t)

        # update physics
        xs1 = sys_forward.(world_dynamics.dynamics, xs, us, t)
        zs1 = obs_forward.(world_dynamics.obs_dynamics, xs, zs, t)
        if loss_model !== nothing
            xs2 = sys_forward.(modeled_dynamics.dynamics, xs, us, t)
            xs_idl .= sys_forward.(modeled_dynamics.dynamics, xs_idl, us_idl, t) .+ (xs1.-xs2)
            xs .= xs1
        end
        xs .= xs1
        zs .= zs1
    end
    logs = Dict(i => write_logs(controllers[i]) for i in 1:N)
    return logs
end

"""
Simulate multiple distributed agents with the given initial states and controllers and
return `(result::TrajectoryData, (logs [, loss_history]))`.

`loss_history` is an `AxisArray` indexed by `(time, id, comp)`, where `comp` could be
either `:u` (action loss) or `:x` (state loss).

# Arguments
- `times::Vector{𝕋}`: the time instances at which the simulation should save the states in
the result. The first element in this sequence should be the time of the given initial states.
- `recorder`: a tuple of the shape `(components::Vector{String}, record_f)`, where
   `record_f` should takes in `(xs, zs, us)` and returns a real matrix of the size
   `(num agents, length(components))`.
"""
function simulate(
    world_dynamics::WorldDynamics{N},
    delay_model::DelayModel,
    framework::ControllerFramework{X,Z,U,Msg,Log},
    init_status::Each{Tuple{X,Z,U}},
    recorder::Tuple{Vector{String},Function},
    times::AbstractVector{𝕋};
    loss_model::Union{RegretLossModel, Nothing} = nothing,
    kwargs...
) where {X,Z,U,Msg,Log,N}
    @assert !isempty(times)

    comps, record_f = recorder
    result = TrajectoryData(times, N, comps)
    data_idx = 1
    t0, tf = times[1], times[end]
    if loss_model !== nothing
        loss_history = AxisArray(-ones(ℝ, tf-t0+1,  N, 2), time=t0:tf, id=1:N, comp=[:u, :x])
    end

    function callback(xs,zs,us,loss,t)
        if t == times[data_idx]
            data = record_f(xs, zs, us)
            @assert size(data) == (N, length(comps)) ("The recorder should "
               * "return a matrix of size (N * Num curves), got size $(size(data))")
            result.values[data_idx, :, :] = data
            data_idx += 1
        end
        if loss_model !== nothing
            @assert loss !== nothing
            loss_history[atvalue(t), :, :] = loss
        end
    end

    logs = simulate(world_dynamics, delay_model, framework, init_status,
        callback, (t0, tf); loss_model, kwargs...)
    sim_out = (loss_model !== nothing) ? (logs, loss_history) : logs

    result, sim_out
end


text:
module test_Plot4

import ModiaMath

#  Desired:
#    using Test
#    using Unitful
#
#  In order that these packages need not to be defined in the user environment, they are included via ModiaMath:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
    t = linspace(0.0, 10.0, 100)
else
    using ModiaMath.Test
    t = range(0.0, stop=10.0, length=100)
end

using ModiaMath.Unitful

result = Dict{AbstractString,Any}()
result["time"] = t * u"s"
result["phi"]  = sin.(t)u"rad"
result["phi2"] = 0.5 * sin.(t)u"rad"
result["w"]    = cos.(t)u"rad/s"
result["w2"]   = 0.6 * cos.(t)u"rad/s"
result["r"]    = hcat(0.4 * cos.(t), 0.5 * sin.(t), 0.3*cos.(t))
result["r2"]   = hcat(0.1*cos.(t), 0.2*cos.(t), 0.3*cos.(t), 0.4*cos.(t), 0.5*cos.(t),
                      0.6*cos.(t), 0.7*cos.(t), 0.8*cos.(t), 0.9*cos.(t), 1.0*cos.(t), 1.1*cos.(t))

ModiaMath.plot(result, [ (:phi), (:phi, :phi2, :w), (:w, :w2) ], heading="Vector of plots")

ModiaMath.plot(result, [ "r", "r2[3]", "r2"], heading="Vector of plots", figure=2)

end


text:
# State
# =====

mutable struct State
    state::Symbol
    State() = new(:init)
end

function start!(state::State)
    state.state = :running
end

function finish!(state::State)
    state.state = :finished
end

function is_running(state::State)
    return state.state == :running
end


text:
# struct LinkResult{T<:FirmIdentifier, D1<:Union{Nothing, Date}, D2<:Union{Nothing, Date}}
#     id::T
#     date_start::D1
#     date_end::D2
# end

# LinkResult(id, d1::Union{Nothing, Date}=nothing, d2::Union{Nothing, Date}=nothing) = LinkResult(id, d1, d2)
# LinkResult(id, d1::Date, ::Missing) = LinkResult(id, d1)
# LinkResult(id, ::Missing, ::Missing) = LinkResult(id)

# struct IDLinkDict{T<:FirmIdentifier, W, D1, D2}
#     links::Dict{T, Vector{LinkResult{W, D1, D2}}}
# end

# get_id(x::LinkResult) = x.id

# function Base.getindex(links::IDLinkDict{T}, id::T) where {T <: FirmIdentifier}
#     get_id.(links.links[id]) |> unique
# end

# function Base.getindex(links::IDLinkDict{T, W, D1, D2}, id::T, date::Date) where {T, W <: FirmIdentifier, D1, D2 <: Nothing}
#     links[id]
# end

# function Base.getindex(links::IDLinkDict{T, W, D1, D2}, id::T, date::Date) where {T, W <: FirmIdentifier, D1, D2 <: Date}
#     res = links.links[id]
#     out = W[]
#     for v in res
#         if v.date_start <= date <= v.date_end
#             push!(out, get_id(v))
#         end
#     end
#     return out |> unique
# end

# function Base.getindex(links::IDLinkDict{T, W, D1, D2}, id::T, date::Date) where {T, W <: FirmIdentifier, D1 <: Date, D2 <: Nothing}
#     res = links.links[id]
#     # assumes the vector is sorted by dates
#     for v in res
#         if v.date_start <= date <= v.date_end
#             return [get_id(v)]
#         end
#     end
#     return W[]
# end

mutable struct LinkTable
    table::String # table name
    id_cols::Vector{String} # id cols as they exist in the database
    date_col_min::Union{Missing, String} # date column as it is in the database, missing for
    # cases where the link is 1:1
    date_col_max::Union{Missing, String} # date column in the database, if missing and min
    # is not missing, then date_col_min is assumed to be the minimum date until a new
    # date is specified
    filters::Dict{String, <:Any} # filters that will be applied when downloading data
    type_translations::Vector{Pair{String, Type{<:FirmIdentifier}}} # translate the column
    # as downloaded to the correct type
end

const LINKTABLE_CACHE = Dict{LinkTable, DataFrame}()

Base.hash(x::LinkTable) = Base.hash("$(x.table) $(join(x.id_cols, ","))")
Base.hash(x::LinkTable, h::UInt) = Base.hash("$(x.table) $(join(x.id_cols, ","))", h)
Base.:(==)(x::LinkTable, y::LinkTable) = Base.hash(x) == Base.hash(y)

LinkTable(::Type{Permno}, ::Type{T}) where {T<:CusipAll} = LinkTable(
    default_tables["crsp_stocknames"],
    ["permno", lowercase(string(T))],
    "namedt",
    "nameenddt",
    Dict{String, Any}("ncusip" => missing),
    ["permno" => Permno, lowercase(string(T)) => T]
)
LinkTable(::Type{T}, ::Type{Permno}) where {T<:CusipAll} = LinkTable(Permno, T)

LinkTable(::Type{K}, ::Type{T}) where {K<:CusipAll,T<:CusipAll} = LinkTable(
    default_tables["crsp_stocknames"],
    [lowercase(string(K)), lowercase(string(T))],
    "namedt",
    "nameenddt",
    Dict{String, Any}("ncusip" => missing),
    [lowercase(string(T)) => T, lowercase(string(K)) => K]
)

LinkTable(::Type{Ticker}, ::Type{T}) where {T<:CusipAll} = LinkTable(
    default_tables["crsp_stocknames"],
    ["ticker", lowercase(string(T))],
    "namedt",
    "nameenddt",
    Dict{String, Any}("ncusip" => missing),
    ["ticker" => Ticker, lowercase(string(T)) => T]
)
LinkTable(::Type{T}, ::Type{Ticker}) where {T<:CusipAll} = LinkTable(Ticker, T)

LinkTable(::Type{Permno}, ::Type{Ticker}) = LinkTable(
    default_tables["crsp_stocknames"],
    ["permno", "ticker"],
    "namedt",
    "nameenddt",
    Dict{String, Any}("ncusip" => missing),
    ["permno" => Permno, "ticker" => Ticker]
)
LinkTable(::Type{Ticker}, ::Type{Permno}) = LinkTable(Permno, Ticker)

LinkTable(::Type{Permno}, ::Type{GVKey}) = LinkTable(
    default_tables["crsp_a_ccm_ccmxpf_lnkhist"],
    ["lpermno", "gvkey"],
    "linkdt",
    "linkenddt",
    Dict{String, Any}(
        "linktype" => ["LU", "LC"],
        "linkprim" => ["P", "C"],
        "lpermno" => missing
    ),
    ["lpermno" => Permno, "gvkey" => GVKey]
)
LinkTable(::Type{GVKey}, ::Type{Permno}) = LinkTable(Permno, GVKey)

LinkTable(::Type{Permno}, ::Type{IbesTicker}) = LinkTable(
    default_tables["wrdsapps_ibcrsphist"],
    ["ticker", "permno"],
    "sdate",
    "edate",
    Dict{String, Any}(
        "score" => 1:4,
        "permno" => missing
    ),
    ["permno" => Permno, "ticker" => IbesTicker]
)
LinkTable(::Type{IbesTicker}, ::Type{Permno}) = LinkTable(Permno, IbesTicker)

LinkTable(::Type{NCusip}, ::Type{IbesTicker}) = LinkTable(
    default_tables["wrdsapps_ibcrsphist"],
    ["ticker", "ncusip"],
    "sdate",
    "edate",
    Dict{String, Any}(
        "score" => 1:4,
        "ncusip" => missing
    ),
    ["ncusip" => NCusip, "ticker" => IbesTicker]
)
LinkTable(::Type{IbesTicker}, ::Type{NCusip}) = LinkTable(NCusip, IbesTicker)

LinkTable(::Type{GVKey}, ::Type{CIK}) = LinkTable(
    default_tables["comp_company"],
    ["gvkey", "cik"],
    missing,
    missing,
    Dict{String, Any}(
        "cik" => missing
    ),
    ["gvkey" => GVKey, "cik" => CIK]
)
LinkTable(::Type{CIK}, ::Type{GVKey}) = LinkTable(GVKey, CIK)

# generic function to convert a pair into the appropriate table
LinkTable(x::Pair{Type{<:FirmIdentifier}, Type{<:FirmIdentifier}}) = LinkTable(x[1], x[2])

function Base.merge(
    t1::LinkTable,
    t2::LinkTable
)
    @assert(t1.table == t2.table)
    LinkTable(
        t1.table,
        vcat(t1.id_cols, t2.id_cols) |> unique,
        t1.date_col_min,
        t1.date_col_max,
        merge(t1.filters, t2.filters),
        vcat(t1.type_translations, t2.type_translations) |> unique,
    )
end

function date_cols(table::LinkTable)
    x = String[]
    if !ismissing(table.date_col_min)
        push!(x, table.date_col_min)
    end
    if !ismissing(table.date_col_max)
        push!(x, table.date_col_max)
    end
    x
end

function date_subset(
    x::Tuple{Date, Date},
    y::Tuple{Date, Date}
)
    x[1] >= y[1] && x[2] <= y[2]
end
function date_checks(
    d_vec1,
    d_vec2
)
    if length(d_vec1) == 1
        return [true]
    end
    out = ones(Bool, length(d_vec1))
    for (i, x) in enumerate(zip(d_vec1, d_vec2))
        for y in zip(d_vec1, d_vec2)
            if x != y && date_subset(x, y)
                out[i] = false
                break
            end
        end
    end
    out
end

"""
function link_table(
    conn,
    table::LinkTable,
    fil_type::Vector{T}=FirmIdentifier[]
) where {T<:FirmIdentifier}

Generic function to download data from a linking table
"""
function link_table(
    conn,
    table::LinkTable,
    fil_type::Vector{T}=FirmIdentifier[]
) where {T<:FirmIdentifier}
    if 0 < length(fil_type) <= 1000
        temp_filter = Dict{String, Any}()
        for (key, val) in table.filters
            temp_filter[key] = val
        end
        table.filters = temp_filter
        col_str_temp = ""
        for (col, t) in table.type_translations
            if t == T
                col_str_temp = col
            end
        end
        table.filters[col_str_temp] = value.(fil_type)
    end
    fil = create_filter(table.filters)
    col_str = join(vcat(table.id_cols, date_cols(table)), ", ")
    query = """
    SELECT DISTINCT $col_str FROM $(table.table)
    $fil
    """
    df = WRDSMerger.run_sql_query(conn, query) |> DataFrame
    for (col, t) in table.type_translations
        df[!, col] = t.(df[:, col])
        rename!(df, col => string(t))
    end
    return df
end

"""
Takes a built tree and converts it into a list of pairs
"""
function find_item(T::Type{<:FirmIdentifier}, node::FirmIdentifierNode)
    if node.data == T
        return node
    end
    out = 0
    for x in node

        if x.data == T
            out = x
            break
        end
        out = find_item(T, x)
        if out != 0
            break
        end
    end

    out
end

function parent_list(node::FirmIdentifierNode; out=Type{<:FirmIdentifier}[])
    push!(out, node.data)
    if parent(node) !== nothing
        return parent_list(parent(node); out)
    end
    out
end

function build_list(T::Type{<:FirmIdentifier}, tree::FirmIdentifierNode)
    bot_node = find_item(T, tree)
    if bot_node == 0
        error("The identifier $T is not in the tree, check that all the links work.")
    end
    out = parent_list(bot_node) |> reverse
    out2 = Pair{Type{<:FirmIdentifier}, Type{<:FirmIdentifier}}[]
    for i in 1:length(out)-1
        push!(out2, out[i] => out[i+1])
    end
    out2
end

function adjust_date_cols(df::DataFrame, table::LinkTable, date_min::Date, date_max::Date)
    if ismissing(table.date_col_max) && !ismissing(table.date_col_min)
        col = Dict(table.type_translations...)[table.id_cols[1]] |> string
        df[!, table.date_col_min] = coalesce.(df[:, table.date_col_min], date_min)
        sort!(df, [col, table.date_col_min])
        gdf = groupby(df, [col])
        df = transform(gdf, table.date_col_min => lead => "date_max")
        df[!, "date_max"] = coalesce.(df[:, "date_max"] .- Day(1), date_max)# I subtract a day since use <= later
        table.date_col_max = "date_max"
    elseif ismissing(table.date_col_min) && ismissing(table.date_col_max)
        table.date_col_min = "date_min"
        table.date_col_max = "date_max"
        df[!, "date_min"] .= date_min
        df[!, "date_max"] .= date_max
    else
        df[!, table.date_col_min] = coalesce.(df[:, table.date_col_min], date_min)
        df[!, table.date_col_max] = coalesce.(df[:, table.date_col_max], date_max)

        """
        CRSP seems to have updated the main table used (crsp.stocknames) so that
        there are a bunch of repeated values where a linked pair is a complete
        subset of another set of dates, this tries to correct for that.

        This is only necessary in the case that all data comes from WRDS,
        if both date columns are missing then this does not matter and if
        one is missing the adjust_date_cols function handles it as well.
        """
        d_cols = [table.date_col_min, table.date_col_max]
        g_cols = [string(t) for (i, t) in table.type_translations]
        subset!(df, d_cols => ByRow((x, y) -> y >= x))
        transform!(df, d_cols => ByRow((x, y) -> x:Day(1):y) => :date_range)
        df = combine(groupby(df, g_cols), :date_range => merge_date_ranges => :date_range)
        transform!(df, :date_range => ByRow(x -> (f=x[1], g=x[end])) => d_cols)
        select!(df, Not(:date_range))
        df = subset(
            groupby(df, g_cols),
            d_cols => (x, y) -> date_checks(x, y)
        )
    end
    df[!, table.date_col_min] = Date.(df[:, table.date_col_min])
    df[!, table.date_col_max] = Date.(df[:, table.date_col_max])
    return df
end

function unique_tables(
    list,
    tables
)
    table_names = String[]
    out_l = Pair{Type{<:FirmIdentifier}, Type{<:FirmIdentifier}}[]
    out_tables = LinkTable[]
    for (l, table) in zip(list, tables)
        if table.table ∉ table_names
            push!(out_l, l)
            push!(out_tables, table)
            push!(table_names, table.table)
        else
            i = findfirst(table.table .== table_names)
            out_tables[i] = merge(out_tables[i], table)
        end
    end
    (out_l, out_tables)
end

"""
    function link_identifiers(
        conn,
        cur_ids::Vector{T},
        dates::Vector{Date},
        new_types::Type{<:FirmIdentifier}...;
        convert_to_values::Bool=true
    ) where {T<:FirmIdentifier}

Provides links between a firm identifier (on a specific date) and
other firm identifiers provided. Generally, these can then be joined into
a DataFrame.

This relies on the ability to build a tree from the provided type to other types,
that function in turn relies on `LinkTable` existing for pairs of functions.

Returns a DataFrame with a column of the ID provided, date, and a column of each
of the requested identifiers. The identifiers have capitalization that follows the
type (i.e., the GVKey column will be titled "GVKey").

## Example

```julia
df = DataFrame(
    cik=["0001341439", "0000004447", "0000723254"],
    date=[Date(2020), Date(2019), Date(2020)]
)

leftjoin(
    df,
    link_identifiers(
        CIK.(df.cik),
        df.date,
        Permno,
        Ticker
    ),
    on=["cik" => "CIK", "date"]
)
```

"""
function link_identifiers(
    conn,
    cur_ids::Vector{T},
    dates::Vector{Date},
    new_types::Type{<:FirmIdentifier}...;
    convert_to_values::Bool=true,
    validate::Bool=true,
    show_tree::Bool=false,
    colnames::AbstractVector{<:Pair}=Pair{String, String}[],
    use_cache::Bool=true
) where {T<:FirmIdentifier}
    df = DataFrame(
        ids=cur_ids,
        date=dates
    ) |> unique
    rename!(df, "ids" => string(T))
    tree = build_tree(T) # tree starts from provided type and goes to all available types
    show_tree && print_tree(tree)
    list = vcat([
        build_list(K, tree) for K in new_types
    ]...) |> unique # pairs of types from new_types to the provided type
    tables = LinkTable.(list) # get the list of tables, already in order
    list, tables = unique_tables(list, tables) # make the tables unique, prevents downloading unnecessary data
    for (l, table) in zip(list, tables)
        if use_cache && haskey(LINKTABLE_CACHE, table)
            new_table = LINKTABLE_CACHE[table]
        else
            new_table = link_table(conn, table)
            if use_cache
                LINKTABLE_CACHE[table] = new_table
            end
        end
        if nrow(new_table) == 0
            continue
        end
        new_table = adjust_date_cols(new_table, table, minimum(df.date), maximum(df.date))
        df = range_join(
            df,
            new_table,
            [string(l[1])],
            [
                Conditions("date", >=, table.date_col_min),
                Conditions("date", <=, table.date_col_max)
            ],
            validate=(validate, false),
            jointype=:left
        )
        select!(df, Not([table.date_col_min, table.date_col_max]))
    end
    select!(df, unique(string.(vcat([string(T), "date"], [x for x in new_types]))))

    if convert_to_values
        for col in names(df)
            if col != "date"
                df[!, col] = value.(df[:, col])
            end
        end
    end
    for col in colnames
        if col[1] ∈ names(df)
            rename!(df, col)
        end
    end
    df

end


text:
using Pkg
# import Pkg;
# Pkg.activate(".")
# ]add CUDA

##
# Pkg.add("Test")
Pkg.add("Zygote")
Pkg.add("BSON")
Pkg.add("ForwardDiff")
Pkg.add("Flux")
Pkg.add("CUDA")
Pkg.add("Functors")
# Pkg.add("Rotations")
# Pkg.add("JLD2")
Pkg.add("Documenter")
Pkg.add("Example")
Pkg.add("DSP")
Pkg.add("Plots")
Pkg.add("FileIO")
Pkg.add("JLD2")
Pkg.add("Lazy")
# Pkg.add(name="GLMakie",version="0.2.9")
Pkg.add("GLMakie")
# Pkg.add("CairoMakie")
# Pkg.add("FFTW")

max_stars_repo_name: aced-differentiate/equiv-docs

max_issues_repo_name: aced-differentiate/equiv-docs

max_forks_repo_name: aced-differentiate/equiv-docs


text:
const BBox = Rect2D{Float32}

const Optional{T} = Union{Nothing, T}

struct RectSides{T}
    left::T
    right::T
    bottom::T
    top::T
end

abstract type Side end

struct Left <: Side end
struct Right <: Side end
struct Top <: Side end
struct Bottom <: Side end
# for protrusion content:
struct TopLeft <: Side end
struct TopRight <: Side end
struct BottomLeft <: Side end
struct BottomRight <: Side end

struct Inner <: Side end
struct Outer <: Side end

abstract type GridDir end
struct Col <: GridDir end
struct Row <: GridDir end

struct RowCols{T <: Union{Number, Vector{Float64}}}
    lefts::T
    rights::T
    tops::T
    bottoms::T
end

"""
    struct Span

Used to specify space that is occupied in a grid. Like 1:1|1:1 for the first square,
or 2:3|1:4 for a rect over the 2nd and 3rd row and the first four columns.
"""
struct Span
    rows::UnitRange{Int64}
    cols::UnitRange{Int64}
end

"""
    mutable struct GridContent{G, T}

Wraps content elements of a `GridLayout`. It keeps track of the `parent`, the `content` and its position in the grid via `span` and `side`.
"""
mutable struct GridContent{G, T} # G should be GridLayout but can't be used before definition
    parent::Optional{G}
    content::T
    span::Span
    side::Side
    needs_update::Node{Bool}
    protrusions_handle::Optional{Function}
    computedsize_handle::Optional{Function}
end

abstract type AlignMode end

"AlignMode that excludes the protrusions from the bounding box."
struct Inside <: AlignMode end

"AlignMode that includes the protrusions within the bounding box, plus paddings."
struct Outside <: AlignMode
    padding::RectSides{Float32}
end
Outside() = Outside(0f0)
Outside(padding::Real) = Outside(RectSides{Float32}(padding, padding, padding, padding))
Outside(left::Real, right::Real, bottom::Real, top::Real) =
    Outside(RectSides{Float32}(left, right, bottom, top))

"AlignMode that is Inside where padding is Nothing and Outside where it is Real."
struct Mixed <: AlignMode
    padding::RectSides{Union{Nothing, Float32}}
end
function Mixed(; left = nothing, right = nothing, bottom = nothing, top = nothing)
    paddings = map((left, right, bottom, top)) do side
        isnothing(side) ? side : Float32(side)
    end
    Mixed(RectSides{Union{Nothing, Float32}}(paddings...))
end

abstract type ContentSize end
abstract type GapSize <: ContentSize end

"""
    struct Auto <: ContentSize

If used as a `GridLayout`'s row / column size and `trydetermine == true`, signals to the `GridLayout` that the row / column should shrink to match the largest determinable element inside.
If no size of a content element can be determined, the remaining space is split between all `Auto` rows / columns according to their `ratio`.

If used as width / height of a layoutable element and `trydetermine == true`, the element's computed width / height will report the auto width / height if it can be determined.
This enables a parent `GridLayout` to adjust its column / rowsize to the element's width / height.
If `trydetermine == false`, the element's computed width / height will report `nothing` even if an auto width / height can be determined, which will prohibit a parent `GridLayout` from adjusting a row / column to the element's width / height.
This is useful to, e.g., prohibit a `GridLayout` from shrinking a column's width to the width of a super title, even though the title's width can be auto-determined.

The `ratio` is ignored if `Auto` is used as an element size.
"""
struct Auto <: ContentSize
    trydetermine::Bool # false for determinable size content that should be ignored
    ratio::Float64 # float ratio in case it's not determinable

    Auto(trydetermine::Bool = true, ratio::Real = 1.0) = new(trydetermine, ratio)
end
Auto(ratio::Real) = Auto(true, ratio)

struct Fixed <: GapSize
    x::Float64
end
struct Relative <: GapSize
    x::Float64
end
struct Aspect <: ContentSize
    index::Int
    ratio::Float64
end

"""
    mutable struct LayoutNodes{T, G}

`T` is the same type parameter of contained `GridContent`, `G` is `GridLayout` which is defined only after `LayoutNodes`.

A collection of `Node`s and an optional `GridContent` that are needed to interface with the MakieLayout layouting system.

- `suggestedbbox::Node{BBox}`: The bounding box that an element should place itself in. Depending on the element's `width` and `height` attributes, this is not necessarily equal to the computedbbox.
- `protrusions::Node{RectSides{Float32}}`: The sizes of content "sticking out" of the main element into the `GridLayout` gaps.
- `computedsize::Node{NTuple{2, Optional{Float32}}}`: The width and height that the element computes for itself if possible (else `nothing`).
- `autosize::Node{NTuple{2, Optional{Float32}}}`: The width and height that the element reports to its parent `GridLayout`. If the element doesn't want to cause the parent to adjust to its size, autosize can hide the computedsize from it by being set to `nothing`.
- `computedbbox::Node{BBox}`: The bounding box that the element computes for itself after it has received a suggestedbbox.
- `gridcontent::Optional{GridContent{G, T}}`: A reference of a `GridContent` if the element is currently placed in a `GridLayout`. This can be used to retrieve the parent layout, remove the element from it or change its position, and assign it to a different layout.
"""
mutable struct LayoutNodes{T, G} # G again GridLayout
    suggestedbbox::Node{BBox}
    protrusions::Node{RectSides{Float32}}
    computedsize::Node{NTuple{2, Optional{Float32}}}
    autosize::Node{NTuple{2, Optional{Float32}}}
    computedbbox::Node{BBox}
    gridcontent::Optional{GridContent{G, T}} # the connecting link to the gridlayout
end

mutable struct GridLayout
    content::Vector{GridContent}
    nrows::Int
    ncols::Int
    rowsizes::Vector{ContentSize}
    colsizes::Vector{ContentSize}
    addedrowgaps::Vector{GapSize}
    addedcolgaps::Vector{GapSize}
    alignmode::AlignMode
    equalprotrusiongaps::Tuple{Bool, Bool}
    needs_update::Node{Bool}
    block_updates::Bool
    layoutnodes::LayoutNodes
    attributes::Attributes
    parentscene::Optional{AbstractPlotting.Scene}
    _update_func_handle::Optional{Function} # stores a reference to the result of on(obs)

    function GridLayout(
        content, nrows, ncols, rowsizes, colsizes,
        addedrowgaps, addedcolgaps, alignmode, equalprotrusiongaps, needs_update,
        layoutnodes, attributes, parentscene)

        gl = new(content, nrows, ncols, rowsizes, colsizes,
            addedrowgaps, addedcolgaps, alignmode, equalprotrusiongaps,
            needs_update, false, layoutnodes, attributes, parentscene, nothing)

        validategridlayout(gl)

        # attach_parent!(gl, parent)

        # on(needs_update) do update
        #     request_update(gl)
        # end

        gl
    end
end

struct AxisAspect
    aspect::Float32
end

struct DataAspect end

abstract type Ticks end

struct AutoLinearTicks <: Ticks
    idealtickdistance::Float32
end

struct ManualTicks <: Ticks
    values::Vector{Float32}
    labels::Vector{String}
end

"""
    struct CustomTicks{F1<:Function, F2<:Function} <: Ticks

For the use of custom functions that compute tick values and labels in an `LAxis`.

- `f_tickvalues::F1`: A function that takes minimum_value, maximum_value, and pixelwidth as arguments and returns a `Float` array of tick values.
- `f_ticklabels::F2`: A function that takes in an array of `Float`s and returns an array of `String` labels.
"""
struct CustomTicks{F1<:Function, F2<:Function} <: Ticks
    f_tickvalues::F1
    f_ticklabels::F2
end

mutable struct LineAxis
    parent::Scene
    protrusion::Node{Float32}
    attributes::Attributes
    decorations::Dict{Symbol, Any}
    tickpositions::Node{Vector{Point2f0}}
    tickvalues::Node{Vector{Float32}}
    ticklabels::Node{Vector{String}}
end

abstract type LObject end

mutable struct LAxis <: AbstractPlotting.AbstractScene
    parent::Scene
    scene::Scene
    xaxislinks::Vector{LAxis}
    yaxislinks::Vector{LAxis}
    limits::Node{BBox}
    layoutnodes::LayoutNodes
    needs_update::Node{Bool}
    attributes::Attributes
    block_limit_linking::Node{Bool}
    decorations::Dict{Symbol, Any}
end

mutable struct LColorbar <: LObject
    parent::Scene
    scene::Scene
    layoutnodes::LayoutNodes
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

mutable struct LText <: LObject
    parent::Scene
    layoutnodes::LayoutNodes
    textobject::AbstractPlotting.Text
    attributes::Attributes
end

mutable struct LRect <: LObject
    parent::Scene
    layoutnodes::LayoutNodes
    rect::AbstractPlotting.Poly
    attributes::Attributes
end

struct LSlider <: LObject
    scene::Scene
    layoutnodes::LayoutNodes
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

struct LButton <: LObject
    scene::Scene
    layoutnodes::LayoutNodes
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

struct LToggle <: LObject
    scene::Scene
    layoutnodes::LayoutNodes
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

abstract type LegendElement end

struct LineElement <: LegendElement
    attributes::Attributes
end

struct MarkerElement <: LegendElement
    attributes::Attributes
end

struct PolyElement <: LegendElement
    attributes::Attributes
end

struct LegendEntry
    elements::Vector{LegendElement}
    attributes::Attributes
end

struct LLegend <: LObject
    scene::Scene
    entries::Node{Vector{LegendEntry}}
    layoutnodes::LayoutNodes
    attributes::Attributes
    decorations::Dict{Symbol, Any}
    entrytexts::Vector{LText}
    entryplots::Vector{Vector{AbstractPlot}}
end

const Indexables = Union{UnitRange, Int, Colon}

struct GridPosition
    layout::GridLayout
    rows::Indexables # this doesn't warrant type parameters I think
    cols::Indexables # as these objects will only be used briefly
end

struct LScene <: AbstractPlotting.AbstractScene
    scene::Scene
    attributes::Attributes
    layoutnodes::MakieLayout.LayoutNodes
end


text:
function save(str::Stream{format"OFF"}, msh::AbstractMesh)
    # writes an OFF geometry file, with colors
    #  see http://people.sc.fsu.edu/~jburkardt/data/off/off.html
    #  for format description
    io = stream(str)
    vts = coordinates(msh)
    fcs = faces(msh)

    cs  = hasproperty(msh, :color) ? msh.color : RGBA{Float32}(0,0,0,1)

    nV = size(vts, 1)
    nF = size(fcs, 1)
    nE = nF*3

    # write the header
    println(io, "OFF")
    println(io, "$nV $nF $nE")
    # write the data
    for v in vts
        println(io, join(Vec{3, Float32}(v), " "))
    end

    for i = 1:nF
        f = fcs[i]
        c = isa(cs, Array) ? RGBA{Float32}(cs[i]) : cs
        facelen = length(f)
        println(io,
            facelen, " ", join(raw.(ZeroIndex.(f)), " "), " ",
            join((red(c), green(c), blue(c), alpha(c)), " ")
        )
    end
    close(io)
end

function load(st::Stream{format"OFF"}; facetype=GLTriangleFace, pointtype=Point3f0)
    io = stream(st)
    points = pointtype[]
    faces = facetype[] # faces might be triangulated, so we can't assume count
    n_points = 0
    n_faces = 0
    found_counts = false
    read_verts = 0

    while !eof(io)
        txt = readline(io)
        if startswith(txt, "#") || isempty(txt) || all(iscntrl, txt) #comment or others
            continue
        elseif found_counts && read_verts < n_points # read verts
            vert = pointtype(parse.(eltype(pointtype), split(txt)))
            if length(vert) == 3
                read_verts += 1
                points[read_verts] = vert
            end
            continue
        elseif found_counts # read faces
            splitted = split(txt)
            facelen = parse(Int, popfirst!(splitted))
            if facelen == 3
                push!(faces, GLTriangleFace(reinterpret(ZeroIndex{Cuint}, parse.(Cuint, splitted[1:3]))))
            elseif facelen == 4
                push!(faces, convert_simplex(facetype, QuadFace{Cuint}(reinterpret(ZeroIndex{Cuint}, parse.(Cuint, splitted[1:4]))))...)
            end
            continue
        elseif !found_counts && all(isdigit, split(txt)[1]) # vertex and face counts
            counts = Int[parse(Int, s) for s in split(txt)]
            n_points = counts[1]
            n_faces = counts[2]
            resize!(points, n_points)
            found_counts = true
        end
    end
    return Mesh(points, faces)
end


text:
using ArgCheck
using LinearAlgebra
using Polynomials

function rand_gram_matrix_rank(n, r)

	"""random gram matrix with rank at most r, and with great likelihood r
	"""

	function normalized_random_vector(r)
		v = rand(ComplexF64, r)
		1/norm(v) .* v
	end

	generating_vectors = 0.
	generating_vectors = hcat([normalized_random_vector(r) for i = 1:n]...)

	generating_vectors' * generating_vectors

end

function incomplete_rank(A::Matrix; atol = 1000eps())

	"""computes the permanent of a general matrix n*n A of rank r < n
	following  Two Algorithmic Results for the Traveling Salesman Problem
	Alexander I. Barvinok

	"""

	n = size(A,1)
	r = rank(A, atol = atol)

	@argcheck r < n

	qr_dec = qr(A)
	G = qr_dec.Q
	B = qr_dec.R

end

atol = 1000eps()
A = rand_gram_matrix_rank(4,3)

n = size(A,1)
r = rank(A, atol = atol)

@argcheck r < n

qr_dec = qr(A)
G = qr_dec.Q
B = qr_dec.R

j = 1

G[j,:][1:r]

L = prod([Polynomial(B[:,j][1:r]) for j in 1:n]).coeffs
R = prod([Polynomial(G[i,:][1:r]) for i in 1:n]).coeffs

[Polynomial(B[:,j][1:r]) for j in 1:n]

n = 35
t = @elapsed Permanents.ryser(ones(n,n))