text:
"""
    TestPaths

This module is intended to be used for testing new path types to
ensure that they are adhering to the AbstractPath API.

# Example

```
# Create a PathSet
ps = PathSet(; symlink=true)

# Select the subset of tests to run
# Inspect TestPaths.TESTALL to see full list
testsets = [
    test_constructor,
    test_registration,
    test_show,
    test_parse,
    test_convert,
    test_components,
    test_parents,
    test_join,
    test_splitext,
    test_basename,
    test_filename,
    test_extensions,
    test_isempty,
    test_norm,
    test_real,
    test_relative,
    test_abs,
    test_isdir,
    test_isfile,
    test_stat,
    test_size,
    test_modified,
    test_created,
    test_cd,
    test_readpath,
    test_walkpath,
    test_read,
    test_write,
    test_mkdir,
    test_cp,
    test_mv,
    test_sync,
    test_symlink,
    test_touch,
    test_tmpname,
    test_tmpdir,
    test_mktmp,
    test_mktmpdir,
    test_download,
]

# Run all the tests
test(ps, testsets)
```
"""
module TestPaths
    using Dates
    using FilePathsBase
    using LinearAlgebra: norm
    using Test

    export PathSet,
        TESTALL,
        test,
        test_constructor,
        test_registration,
        test_show,
        test_parse,
        test_convert,
        test_components,
        test_parents,
        test_join,
        test_splitext,
        test_basename,
        test_filename,
        test_extensions,
        test_isempty,
        test_norm,
        test_real,
        test_relative,
        test_abs,
        test_isdir,
        test_isfile,
        test_stat,
        test_size,
        test_modified,
        test_created,
        test_issocket,
        test_isfifo,
        test_ischardev,
        test_isblockdev,
        test_ismount,
        test_isexecutable,
        test_isreadable,
        test_iswritable,
        test_cd,
        test_readpath,
        test_walkpath,
        test_read,
        test_write,
        test_mkdir,
        test_cp,
        test_mv,
        test_sync,
        test_symlink,
        test_touch,
        test_tmpname,
        test_tmpdir,
        test_mktmp,
        test_mktmpdir,
        test_chown,
        test_chmod,
        test_download

    # foo, bar, baz, qux, quux, quuz, corge, grault, garply, waldo, fred, plugh, xyzzy, and thud
    #=
    To Test:

    - [X] Construction (T, Path, p"...")
    - [X] show, parse, convert
    - [X] anchor, drive, root
    - [X] hasparent
    - [X] parent
    - [X] parents
    - [X] *, /, join
    - [X] basename
    - [X] filename
    - [X] isempty
    - [X] norm
    - [X] abs?
    - [X] relative?
    - [X] exists
    - [X] real
    - [X] stat/lstat
    - [X] size (could this default to using a FileBuffer?)
    - [X] modified
    - [X] created
    - [X] modified
    - [X] created
    - [X] issocket
    - [X] isfifo
    - [X] ischardev
    - [X] isblockdev
    - [X] ismount
    - [X] isdir
    - [X] isfile
    - [X] isexecutable
    - [X] iswritable
    - [X] isreadable
    - [X] cd?
    - [X] mkdir
    - [X] rm
    - [X] read
    - [X] write
    - [X] cp
    - [X] mv
    - [X] symlink
    - [X] touch
    - [X] tmpname
    - [X] tmpdir
    - [X] mktmp
    - [X] mktmpdir
    - [X] readpath
    - [X] walkpath
    - [X] open
    - [X] chmod
    - [X] chown
    - [X] download
    =#
    """
        PathSet(root::AbstractPath=tmpdir(); symlink=false)

    Constructs a common test path hierarchy to running shared API tests.

    Hierarchy:

    ```
    root
    |-- foo
    |   |-- baz.txt
    |-- bar
    |   |-- qux
    |       |-- quux.tar.gz
    |-- fred
    |   |-- plugh
    ```
    """
    struct PathSet{P<:AbstractPath}
        root::P
        foo::P
        baz::P
        bar::P
        qux::P
        quux::P
        fred::P
        plugh::P
        link::Bool
    end

    function PathSet(root=tmpdir() / "pathset_root"; symlink=false)
        root = abs(root)

        PathSet(
            root,
            root / "foo",
            root / "foo" / "baz.txt",
            root / "bar",
            root / "bar" / "qux",
            root / "bar" / "qux" / "quux.tar.gz",
            root / "fred",
            root / "fred" / "plugh",
            symlink,
        )
    end

    function initialize(ps::PathSet)
        @info "Initializing $(typeof(ps))"
        mkdir.([ps.foo, ps.qux, ps.fred]; recursive=true, exist_ok=true)
        write(ps.baz, "Hello World!")
        write(ps.quux, "Hello Again!")

        # If link is true then plugh is a symlink to foo
        if ps.link
            symlink(ps.foo, ps.plugh)
        else
            touch(ps.plugh)
        end
    end

    function test_constructor(ps::PathSet{P}) where P <: AbstractPath
        @testset "Constructor" begin
            str = string(ps.root)
            @test P(str) == ps.root
        end
    end

    function test_registration(ps::PathSet{P}) where P <: AbstractPath
        @testset "Path constructor" begin
            str = string(ps.root)
            @test FilePathsBase.ispathtype(P, str)
            @test Path(str) == ps.root
            @test p"foo/bar" == Path("foo/bar")
        end
    end

    function test_show(ps::PathSet)
        @testset "show" begin
            str = string(ps.root)
            # For windows paths
            str = replace(str, "\\" => "/")
            @test sprint(show, ps.root; context=:compat => true) == "p\"$str\""
            # TODO: Figure out why this is broken.
            @test_broken sprint(show, ps.root; context=:compat => false) == str
        end
    end

    function test_parse(ps::PathSet{P}) where P <: AbstractPath
        @testset "parsing" begin
            str = string(ps.root)
            @test parse(P, str) == ps.root
        end
    end

    function test_convert(ps::PathSet{P}) where P <: AbstractPath
        @testset "convert" begin
            str = string(ps.root)
            @test convert(P, str) == ps.root
            @test convert(String, ps.root) == str
        end
    end

    function test_components(ps::PathSet)
        @testset "components" begin
            str = string(ps.root)
            @test ps.root.anchor == ps.root.drive * ps.root.root
            @test ps.quux.segments[end-2:end] == ("bar", "qux", "quux.tar.gz")

            # Check that isless on the path segments works
            @test ps.bar < ps.foo
            @test sort([ps.foo, ps.bar, ps.fred]) == [ps.bar, ps.foo, ps.fred]
        end
    end

    function test_parents(ps::PathSet)
        @testset "parents" begin
            @test parent(ps.foo) == ps.root
            @test parent(ps.qux) == ps.bar

            @test hasparent(ps.qux)
            _parents = parents(ps.qux)
            @test _parents[end] == ps.bar
            @test _parents[end - 1] == ps.root

            # Test that relative paths with no parents return p"."
            @test parent(Path(basename(ps.foo))) == p"."

            # Test that parent on p"." should be ===
            path = p"."
            @test parent(path) === path
        end
    end

    function test_join(ps::PathSet)
        @testset "join" begin
            @test join(ps.root, "bar") == ps.bar
            @test ps.root / "foo" / "baz.txt" == ps.baz
            @test ps.root / "foobaz.txt" == ps.root / "foo" * "baz.txt"
        end
    end

    function test_splitext(ps::PathSet)
        @testset "splitext" begin
            @test splitext(ps.foo) == (ps.foo, "")
            @test splitext(ps.baz) == (ps.root / "baz", ".txt")
            @test splitext(ps.quux) == (ps.qux / "quux.tar", ".gz")
        end
    end
    function test_basename(ps::PathSet)
        @testset "basename" begin
            @test basename(ps.foo) == "foo"
            @test basename(ps.baz) == "baz.txt"
            @test basename(ps.quux) == "quux.tar.gz"
        end
    end

    function test_filename(ps::PathSet)
        @testset "filename" begin
            @test filename(ps.foo) == "foo"
            @test filename(ps.baz) == "baz"
            @test filename(ps.quux) == "quux"
        end
    end

    function test_extensions(ps::PathSet)
        @testset "extensions" begin
            @test extension(ps.foo) == ""
            @test extension(ps.baz) == "txt"
            @test extension(ps.quux) == "gz"
            @test extensions(ps.foo) == []
            @test extensions(ps.baz) == ["txt"]
            @test extensions(ps.quux) == ["tar", "gz"]
        end
    end

    function test_isempty(ps::PathSet{P}) where P <: AbstractPath
        @testset "isempty" begin
            @test !isempty(ps.foo)
            @test isempty(P())
        end
    end

    function test_norm(ps::PathSet)
        @testset "norm" begin
            @test norm(ps.bar / ".." / "foo") == ps.foo
            @test norm(ps.bar / ".") == ps.bar
        end
    end

    function test_real(ps::PathSet)
        @testset "real" begin
            # NOTE: We call `real` on ps.bar in the `norm` case because on
            # macOS the temp directory may include a symlink.
            @test real(ps.bar / ".." / "foo") == norm(real(ps.bar) / ".." / "foo")
            @test real(ps.bar / ".") == norm(real(ps.bar) / ".")

            if ps.plugh !== nothing
                if isa(ps.plugh, WindowsPath)
                    @test_broken real(ps.plugh) == real(ps.foo)
                else
                    @test real(ps.plugh) == real(ps.foo)
                end
            end
        end
    end

    function test_relative(ps::PathSet)
        @testset "relative" begin
            @test relative(ps.foo, ps.qux).segments == ("..", "..", "foo")
        end
    end

    function test_abs(ps::PathSet)
        @testset "abs" begin
            @test isabs(ps.root) || isabs(abs(ps.root))
        end
    end

    function test_isdir(ps::PathSet)
        @testset "isdir" begin
            @test isdir(ps.foo)
            @test !isdir(ps.baz)
        end
    end

    function test_isfile(ps::PathSet)
        @testset "isfile" begin
            @test isfile(ps.baz)
            @test !isfile(ps.foo)
        end
    end

    function test_stat(ps::PathSet)
        @testset "stat" begin
            s = stat(ps.root)
            fields = fieldnames(typeof(s))

            @test :size in fields
            @test :ctime in fields
            @test :mtime in fields
            @test :user in fields
            @test :mode in fields

            if ps.link
                @test lstat(ps.plugh) != stat(ps.plugh)
            end

            str = sprint(show, s)
        end
    end

    # Minimal testing of issocket, isfifo, ischardev, isblockdev and ismount which
    # people won't typically include.
    function test_issocket(ps::PathSet)
        @test !issocket(ps.root)
    end

    function test_isfifo(ps::PathSet)
        @test !isfifo(ps.root)
    end

    function test_ischardev(ps::PathSet)
        @test !ischardev(ps.root)
    end

    function test_isblockdev(ps::PathSet)
        @test !isblockdev(ps.root)
    end

    function test_ismount(ps::PathSet)
        @test !ismount(ps.root)
    end

    function test_size(ps::PathSet)
        @testset "size" begin
            @test size(ps.baz) > 0
        end
    end

    function test_modified(ps::PathSet)
        @testset "modified" begin
            @test isa(modified(ps.baz), Dates.AbstractDateTime)
            @test modified(ps.baz) >= modified(ps.root)
        end
    end

    function test_created(ps::PathSet)
        @testset "created" begin
            @test isa(created(ps.baz), Dates.AbstractDateTime)
            @test created(ps.baz) >= created(ps.root)
        end
    end

    function test_isexecutable(ps::PathSet{P}) where P <: AbstractPath
        @testset "isexecutable" begin
            # I'm not entirely sure how to test this generally
            @test !isexecutable(ps.baz)

            # Directories should be executable for system paths
            if P <: SystemPath
                @test isexecutable(ps.foo)
            end
        end
    end

    function test_isreadable(ps::PathSet)
        @testset "isreadable" begin
            # Our test files should be readable by default
            @test isreadable(ps.baz)
            @test isreadable(ps.quux)

            # Our test directories should also be readable by default
            @test isreadable(ps.foo)
            @test isreadable(ps.qux)
        end
    end

    function test_iswritable(ps::PathSet)
        @testset "iswritable" begin
            # Our test files should be writable by default
            @test iswritable(ps.baz)
            @test iswritable(ps.quux)

            # Our test directories should also be writable by default
            @test iswritable(ps.foo)
            @test iswritable(ps.qux)
        end
    end

    function test_cd(ps::PathSet{P}) where P <: AbstractPath
        if P <: SystemPath
            @testset "cd" begin
                init_path = cwd()

                cd(ps.foo) do
                    @test cwd() != real(init_path)
                    @test cwd() == real(ps.foo)
                end

                @test cwd() == real(init_path)

                cd(ps.qux)
                @test cwd() != real(init_path)
                @test cwd() == real(ps.qux)
                cd(init_path)
                @test cwd() == real(init_path)
            end
        end
    end

    function test_readpath(ps::PathSet)
        @testset "readpath" begin
            @test readdir(ps.root) == ["bar", "foo", "fred"]
            @test readdir(ps.qux) == ["quux.tar.gz"]
            @test readpath(ps.root) == [ps.bar, ps.foo, ps.fred]
            @test readpath(ps.qux) == [ps.quux]
        end
    end

    function test_walkpath(ps::PathSet)
        @testset "walkpath" begin
            topdown = [ps.bar, ps.qux, ps.quux, ps.foo, ps.baz, ps.fred, ps.plugh]
            bottomup = [ps.quux, ps.qux, ps.bar, ps.baz, ps.foo, ps.plugh, ps.fred]

            @test collect(walkpath(ps.root; topdown=true)) == topdown
            @test collect(walkpath(ps.root; topdown=false)) == bottomup
        end
    end

    function test_read(ps::PathSet)
        @testset "read" begin
            @test read(ps.baz, String) == "Hello World!"
            open(ps.quux, "r") do io
                @test read(io, String) == "Hello Again!"
            end
        end
    end

    function test_write(ps::PathSet)
        @testset "write" begin
            write(ps.baz, "Goodbye World!")
            @test read(ps.baz, String) == "Goodbye World!"

            @testset "truncate" begin
                open(ps.quux, "w") do io
                    write(io, "Hello?")
                end
                @test read(ps.quux, String) == "Hello?"
            end

            @testset "append" begin
                open(ps.quux, "a") do io
                    write(io, " Did you need something?")
                end
                @test read(ps.quux, String) == "Hello? Did you need something?"
            end

            @testset "read/write" begin
                open(ps.quux, "w+") do io
                    write(io, "Goodnight!")
                    seekstart(io)
                    @test read(io, String) == "Goodnight!"
                end

                open(ps.quux, "a+") do io
                    write(io, " Zzzz")
                    seekstart(io)
                    @test read(io, String) == "Goodnight! Zzzz"
                end
            end
        end
    end

    function test_mkdir(ps::PathSet)
        @testset "mkdir" begin
            garply = ps.root / "corge" / "grault" / "garply"
            @test_throws ErrorException mkdir(garply)
            mkdir(garply; recursive=true)
            @test exists(garply)
            @test_throws ErrorException mkdir(garply; recursive=true)
            mkdir(garply; recursive=true, exist_ok=true)
            rm(ps.root / "corge"; recursive=true)
            @test !exists(garply)
        end
    end

    function test_cp(ps::PathSet)
        @testset "cp" begin
            cp(ps.foo, ps.qux / "foo"; force=true)
            @test exists(ps.qux / "foo" / "baz.txt")
            @test_throws ArgumentError cp(ps.foo, ps.qux / "foo")
            cp(ps.foo, ps.qux / "foo"; force=true)
            rm(ps.qux / "foo"; recursive=true)
        end
    end

    function test_mv(ps::PathSet)
        @testset "mv" begin
            garply = ps.root / "corge" / "grault" / "garply"
            mkdir(garply; recursive=true, exist_ok=true)
            @test exists(garply)
            mv(ps.root / "corge", ps.foo / "corge"; force=true)
            @test exists(ps.foo / "corge" / "grault" / "garply")
            rm(ps.foo / "corge"; recursive=true)
        end
    end

    function test_sync(ps::PathSet)
        @testset "sync" begin
            @testset "empty destination" begin
                sync(ps.foo, ps.qux / "foo")
                @test exists(ps.qux / "foo" / "baz.txt")

                # Test that the copied baz file has a newer modified time
                baz_t = modified(ps.qux / "foo" / "baz.txt")
                @test modified(ps.baz) < baz_t
            end

            @testset "empty source" begin
                @test_throws ArgumentError sync(ps.root / "quux", ps.foo)
            end

            @testset "new source" begin
                # Don't cp unchanged files when a new file is added
                # NOTE: sleep before we make a new file, so it's clear that the
                # modified time has changed.
                baz_t = modified(ps.qux / "foo" / "baz.txt")
                sleep(1)
                write(ps.foo / "test.txt", "New src")
                sync(ps.foo, ps.qux / "foo")
                @test exists(ps.qux / "foo" / "test.txt")
                @test read(ps.qux / "foo" / "test.txt", String) == "New src"
                @test modified(ps.qux / "foo" / "baz.txt") == baz_t
                @test modified(ps.qux / "foo" / "test.txt") > baz_t
            end

            @testset "new destination" begin
                # Newer file of the same size is likely the result of an upload which
                # will always have a newer last modified time.
                test_t = modified(ps.foo / "test.txt")
                sleep(1)
                write(ps.qux / "foo" / "test.txt", "New dst")
                @test modified(ps.qux / "foo" / "test.txt") > test_t
                sync(ps.foo, ps.qux / "foo")
                @test read(ps.qux / "foo" / "test.txt", String) == "New dst"
                @test modified(ps.qux / "foo" / "test.txt") > test_t
            end

            @testset "no delete" begin
                # Test not deleting a file on sync
                rm(ps.foo / "test.txt")
                sync(ps.foo, ps.qux / "foo")
                @test exists(ps.qux / "foo" / "test.txt")
            end

            @testset "delete" begin
                # Test passing delete flag
                sync(ps.foo, ps.qux / "foo"; delete=true)
                @test !exists(ps.qux / "foo" / "test.txt")
                rm(ps.qux / "foo"; recursive=true)
            end

            @testset "mixed types" begin
                @testset "directory -> file" begin
                    @test_throws ArgumentError sync(ps.foo, ps.quux)
                end

                @testset "file -> directory" begin
                    @test_throws ArgumentError sync(ps.quux, ps.foo)
                end
            end

            @testset "walkpath order" begin
                # Test a condtion where the index could reorder the walkpath order.
                tmp_src = ps.root / "tmp-src"
                mkdir(tmp_src)
                src_file = tmp_src / "file1"
                write(src_file, "Hello World!")

                src_folder = tmp_src / "folder1"
                mkdir(src_folder)
                src_folder_file = src_folder / "file2"
                write(src_folder_file, "") # empty file

                src_folder2 = src_folder / "folder2"  # nested folders
                mkdir(src_folder2)
                src_folder2_file = src_folder2 / "file3"
                write(src_folder2_file, "Test")

                tmp_dst = ps.root / "tmp_dst"
                mkdir(tmp_dst)
                sync(tmp_src, tmp_dst)
                @test exists(tmp_dst / "folder1" / "folder2" / "file3")
                rm(tmp_src; recursive=true)
                rm(tmp_dst; recursive=true)
            end
        end
    end

    function test_symlink(ps::PathSet)
        if ps.link
            @testset "symlink" begin
                @test_throws ErrorException symlink(ps.foo, ps.plugh)
                symlink(ps.foo, ps.plugh; exist_ok=true, overwrite=true)
                symlink(ps.foo, ps.plugh; exist_ok=true)
                @test_throws ErrorException symlink(ps.foo / "thud", ps.plugh; exist_ok=true, overwrite=true)
            end
        end
    end

    function test_touch(ps::PathSet)
        @testset "touch" begin
            newfile = ps.root / "newfile"
            touch(newfile)
            @test exists(newfile)
            rm(newfile)
        end
    end

    function test_tmpname(ps::PathSet)
        @testset "tmpname" begin
            @test isa(tmpname(), AbstractPath)
            @test hasparent(tmpname())
            @test exists(parent(tmpname()))
        end
    end

    function test_tmpdir(ps::PathSet)
        @testset "tmpname" begin
            @test isa(tmpdir(), AbstractPath)
            @test exists(tmpdir())
            @test isdir(tmpdir())
        end
    end

    function test_mktmp(ps::PathSet)
        @testset "mktmp" begin
            mktmp(ps.root) do path, io
                @test exists(path)
                @test iswritable(io)
                write(io, "Foobar")
                seekstart(io)
                @test read(io, String) == "Foobar"
            end
        end
    end

    function test_mktmpdir(ps::PathSet)
        @testset "mktmpdir" begin
            mktmpdir(ps.root) do path
                @test exists(path)
                write(path / "test.txt", "Foobar")
                @test read(path / "test.txt", String) == "Foobar"
            end
        end
    end

    function test_chown(ps::PathSet)
        @testset "chown" begin
            newfile = ps.root / "newfile"
            touch(newfile)

            if haskey(ENV, "USER")
                if ENV["USER"] == "root"
                    chown(newfile, "nobody", "nogroup"; recursive=true)
                else
                    @test_throws ErrorException chown(newfile, "nobody", "nogroup"; recursive=true)
                end
            end

            rm(newfile)
        end
    end

    function test_chmod(ps::PathSet)
        @testset "chmod" begin
            newfile = ps.root / "newfile"
            newpath = ps.root / "thud"

            touch(newfile)
            mkdir(newpath)
            chmod(newfile, user=(READ+WRITE+EXEC), group=(READ+EXEC), other=READ)
            @test string(mode(newfile)) == "-rwxr-xr--"
            @test isexecutable(newfile)
            @test iswritable(newfile)
            @test isreadable(newfile)

            chmod(newfile, "-x")
            @test !isexecutable(newfile)

            @test string(mode(newfile)) == "-rw-r--r--"
            chmod(newfile, "+x")
            write(newfile, "foobar")
            @test read(newfile, String) == "foobar"
            chmod(newfile, "u=rwx")

            open(newfile, "r") do io
                @test read(io, String) == "foobar"
            end

            chmod(newpath, mode(newfile); recursive=true)
        end
    end

    function test_download(ps::PathSet)
        @testset "download" begin
            rm(ps.foo / "README.md"; force=true)
            download(
                "https://github.com/rofinn/FilePathsBase.jl/blob/master/README.md",
                ps.foo / "README.md"
            )
            @test exists(ps.foo / "README.md")

            # Test downloading from another path
            download(ps.foo / "README.md", ps.qux / "README.md")
            @test exists(ps.qux / "README.md")

             # Test downloading to a string
            download(ps.foo / "README.md", string(ps.fred / "README.md"))
            @test exists(ps.fred / "README.md")
            rm.([ps.foo / "README.md", ps.qux / "README.md", ps.fred / "README.md"])
        end
    end

    TESTALL = [
        test_constructor,
        test_registration,
        test_show,
        test_parse,
        test_convert,
        test_components,
        test_parents,
        test_join,
        test_basename,
        test_filename,
        test_extensions,
        test_isempty,
        test_norm,
        test_real,
        test_relative,
        test_abs,
        test_isdir,
        test_isfile,
        test_stat,
        test_size,
        test_modified,
        test_created,
        test_issocket,
        test_isfifo,
        test_ischardev,
        test_isblockdev,
        test_ismount,
        test_isexecutable,
        test_isreadable,
        test_iswritable,
        test_cd,
        test_readpath,
        test_walkpath,
        test_read,
        test_write,
        test_mkdir,
        test_cp,
        test_mv,
        test_symlink,
        test_touch,
        test_tmpname,
        test_tmpdir,
        test_mktmp,
        test_mktmpdir,
        test_chown,
        test_chmod,
        test_download,
    ]

    function test(ps::PathSet, test_sets=TESTALL)
        try
            initialize(ps)

            for ts in test_sets
                ts(ps)
            end
        finally
            rm(ps.root; recursive=true, force=true)
        end
    end
end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

# SourceLocation print functions
function println(io::IO, location::SourceLocation)
     print(io, "location = [")
     if location.file_name == ""
          print(io, "nothing, ")
     else
          print(io, "$(location.file_name), ")
     end
     print(io, "(", location.line_num,",", location.col_num, ")")
end
print(io::IO, location::SourceLocation) = println(io, location)
println(location::SourceLocation) = println(stdout, location)
print(location::SourceLocation) = print(stdout,location)

# Token print functions
function println(io::IO, token::Token)
     print(io, "Token: \t")
     println(io, "token.location = ", token.location, " , token.value = ", token.value)
end
print(io::IO, token::Token) = println(io, token)
println(token::Token) = println(stdout, token)
print(token::Token) = print(stdout, token)

function println(inputstream::InputStream, scene::Scene)
     expect_symbol(inputstream, "(")
     token = read_token(inputstream)
     if typeof(token.value) == SymbolToken
          println("SymbolToken: ", token.location, ", ", token.value.symbol)
     elseif typeof(token.value) == LiteralNumberToken
          println("LiteralNumberToken: ", token.location, ",  ", token.value.number)
     elseif typeof(token.value) == StringToken
          println("StringToken: ", token.location, ",  ", token.value.string)
     elseif typeof(token.value) == KeywordToken
          println("KeywordToken: ", token.location, ",  ", token.value.keyword)
     elseif typeof(token.value) == IdentifierToken
          print("IdentifierToken: ", token.location, ", identifier = ", token.value.identifier, ", value =")
          if token.value.identifier ∈ keys(scene.float_variables)
               println(scene.float_variables[token.value.identifier])
          elseif token.value.identifier ∈ keys(scene.vector_variables)
               println(scene.vector_variables[token.value.identifier])
          elseif token.value.identifier ∈ keys(scene.color_variables)
               println(scene.color_variables[token.value.identifier])
          else
               println(typeof(token.value.identifier))
          end
     end
     expect_symbol(inputstream, ")")
end

function parse_function(inputstream::InputStream, scene::Scene, open::Bool=false)
     token = read_token(inputstream)
     result = ""

     while true

          if typeof(token.value) == LiteralNumberToken
               result *= repr(token.value.number)

          elseif typeof(token.value) == StringToken
               result *= repr(token.value.string)

          elseif typeof(token.value) == SymbolToken && token.value.symbol ∈ OPEN_BRACKETS
               result *= token.value.symbol*parse_function(inputstream, scene, true)
               expect_symbol(inputstream, close_bracket(token.value.symbol))
               result *= close_bracket(token.value.symbol)

          elseif typeof(token.value) == SymbolToken && token.value.symbol ∉ CLOSED_BRACKETS && open==true
               result *= token.value.symbol

          elseif typeof(token.value) == KeywordToken
               keyword = token.value.keyword
               if keyword ∈ [DIFFUSE,  SPECULAR]
                    unread_token(inputstream, token)
                    result *= repr(parse_brdf(inputstream, scene))
               elseif keyword ∈ [UNIFORM,  CHECKERED,  IMAGE]
                    unread_token(inputstream, token)
                    result *= repr(parse_pigment(inputstream, scene))
               else
                    throw(GrammarError(token.location,
                    "keyword '$(keyword)' do not define anything that can be inside a function"))
               end

          elseif typeof(token.value) == IdentifierToken
               variable_name = token.value.identifier

               if (variable_name ∈ keys(SYM_NUM))
                    result *=  string(SYM_NUM[token.value.identifier])
               elseif variable_name ∈ keys(scene.float_variables)
                    result *= repr(scene.float_variables[variable_name])
               elseif variable_name ∈ keys(scene.string_variables)
                    result *= repr(scene.string_variables[variable_name])
               elseif variable_name ∈ keys(scene.bool_variables)
                    result *= repr(scene.bool_variables[variable_name])
               elseif variable_name ∈ keys(scene.vector_variables)
                    result *= repr(scene.vector_variables[variable_name])
               elseif variable_name ∈ keys(scene.color_variables)
                    result *= repr(scene.color_variables[variable_name])
               elseif variable_name ∈ keys(scene.pigment_variables)
                    result *= repr(scene.pigment_variables[variable_name])
               elseif variable_name ∈ keys(scene.brdf_variables)
                    result *= repr(scene.brdf_variables[variable_name])
               elseif variable_name ∈ keys(scene.transformation_variables)
                    result *= repr(scene.transformation_variables[variable_name])

               elseif isdefined(Raytracing, Symbol(variable_name)) || isdefined(Base, Symbol(variable_name))
                    result *= variable_name
                    var = eval(Symbol(variable_name))

                    if var isa Function
                         open_sym = expect_symbol(inputstream, ["{", "[", "("])
                         result *= open_sym*parse_function(inputstream, scene, true)
                         closed_sym = close_bracket(open_sym)
                         expect_symbol(inputstream, closed_sym)
                         result *= closed_sym
                    else
                         result *= var
                    end

               else
                    throw(GrammarError(token.location,
                    "identifier '$(variable_name)' do not define anything that can be compared"))
               end
          else
               unread_token(inputstream, token)
               break
          end

          token = read_token(inputstream)
     end

     return result
end

"""
     return_token_value(
          inputstream::InputStream,
          scene::Scene
          ) :: Union{Float64, Bool, String, Vec, RGB{Float32}, Pigment, BRDF}

Return the value of the token readed from `inputstream`.
If the token is an IdentifierToken, return the value associated with that
identifier insdie the given `scene`.
Throws `GrammarError` if an error occurs.
Call internally the following parsing functions:
- [`read_token`](@ref)
- [`unread_token`](@ref)
- [`parse_vector`](@ref)
- [`parse_color`](@ref)
- [`parse_pigment`](@ref)
- [`parse_brdf`](@ref)

See also: [`Vec`](@ref), [`Pigment`](@ref), [`BRDF`](@ref),
[`InputStream`](@ref), [`Scene`](@ref), [`GrammarError`](@ref)
"""
function return_token_value(inputstream::InputStream, scene::Scene, open::Bool = false)
     token = read_token(inputstream)
     result = ""

     while true

          if typeof(token.value) == LiteralNumberToken
               result *= repr(token.value.number)

          elseif typeof(token.value) == StringToken
               result *= repr(token.value.string)

          elseif typeof(token.value) == SymbolToken && token.value.symbol ∈ OPERATIONS
               result *= token.value.symbol

          elseif typeof(token.value) == SymbolToken && token.value.symbol=="["
               unread_token(inputstream, token)
               result *= repr(parse_vector(inputstream, scene))

          elseif typeof(token.value) == SymbolToken && token.value.symbol=="<"
               unread_token(inputstream, token)
               result *= repr(parse_color(inputstream, scene))

          elseif typeof(token.value) == SymbolToken && token.value.symbol=="("
               result *= "("*return_token_value(inputstream, scene, true)
               expect_symbol(inputstream, ")")
               result *= ")"

          elseif typeof(token.value) == KeywordToken
               keyword = token.value.keyword
               if keyword ∈ [DIFFUSE,  SPECULAR]
                    unread_token(inputstream, token)
                    result *= repr(parse_brdf(inputstream, scene))
               elseif keyword ∈ [UNIFORM,  CHECKERED,  IMAGE]
                    unread_token(inputstream, token)
                    result *= repr(parse_pigment(inputstream, scene))
               else
                    throw(GrammarError(token.location,
                    "keyword '$(keyword)' do not define anything that can be compared"))
               end

          elseif typeof(token.value) == IdentifierToken
               variable_name = token.value.identifier

               if (variable_name ∈ keys(SYM_NUM))
                    result *=  string(SYM_NUM[token.value.identifier])
               elseif variable_name ∈ keys(scene.float_variables)
                    result *= repr(scene.float_variables[variable_name])
               elseif variable_name ∈ keys(scene.string_variables)
                    result *= repr(scene.string_variables[variable_name])
               elseif variable_name ∈ keys(scene.bool_variables)
                    result *= repr(scene.bool_variables[variable_name])
               elseif variable_name ∈ keys(scene.vector_variables)
                    result *= repr(scene.vector_variables[variable_name])
               elseif variable_name ∈ keys(scene.color_variables)
                    result *= repr(scene.color_variables[variable_name])
               elseif variable_name ∈ keys(scene.pigment_variables)
                    result *= repr(scene.pigment_variables[variable_name])
               elseif variable_name ∈ keys(scene.brdf_variables)
                    result *= repr(scene.brdf_variables[variable_name])
               elseif variable_name ∈ keys(scene.transformation_variables)
                    result *= repr(scene.transformation_variables[variable_name])

               elseif isdefined(Raytracing, Symbol(variable_name))
                    result *= variable_name
                    var = eval(Symbol(variable_name))

                    if var isa Function
                         open_sym = expect_symbol(inputstream, ["{", "[", "("])
                         result *= open_sym*parse_function(inputstream, scene)
                         closed_sym = close_bracket(open_sym)
                         expect_symbol(inputstream, closed_sym)
                         result *= closed_sym
                    else
                         result *= var
                    end

               elseif isdefined(Base, Symbol(variable_name))
                    result *= variable_name
                    var = eval(Symbol(variable_name))

                    if var isa Function
                         open_sym = expect_symbol(inputstream, ["{", "[", "("])
                         result *= sym*parse_function(inputstream, scene)
                         closed_sym = close_bracket(open_sym)
                         expect_symbol(inputstream, closed_sym)
                         result *= closed_sym
                    else
                         result *= var
                    end
               else
                    throw(GrammarError(token.location,
                    "identifier '$(variable_name)' do not define anything that can be compared"))
               end
          else
               unread_token(inputstream, token)
               break
          end

          token = read_token(inputstream)
     end

     if open==true
          return result
     else
          return eval(Meta.parse(result))
     end
end

"""
     assert(inputstream::InputStream, scene::Scene)

Parse an assertion from the given `inputstream` and return it.
Throws `AssertionError` if the assertion is false, nothing otherwise.

Call internally the following parsing functions:
- [`expect_symbol`](@ref)
- [`expect_keywords`](@ref)
- [`expect_number`](@ref)
- [`expect_string`](@ref)
- [`return_token_value`](@ref)

## Examples

```text
assert(1, 1)             # Checks that 1==1
assert(1, 2, "<")        # Checks that 1<2
float var(1.0)           # Define var as a variable with value 1.0
assert(var, 2, "<=")     # Checks that var<=2
```

See also: [`InputStream`](@ref), [`Scene`](@ref)
"""
function assert(inputstream::InputStream, scene::Scene)
     expect_symbol(inputstream, "(")
     value1 =  return_token_value(inputstream, scene)
     expect_symbol(inputstream, ",")
     value2 =  return_token_value(inputstream, scene)

     token = read_token(inputstream)
     if typeof(token.value) == SymbolToken && token.value.symbol == ","
          unread_token(inputstream, token)
          expect_symbol(inputstream,",")
          operator = expect_string(inputstream, scene)
          expect_symbol(inputstream, ")")

          if operator ∈ ["=", "=="]
               @assert value1 == value2 "$(value1) and $(value2) are not equal!"
          elseif operator == "<"
               @assert value1 < value2 "$(value1) is not < $(value2)!"
          elseif operator == "<="
               @assert value1 <= value2 "$(value1) is not <= $(value2)!"
          elseif operator == ">"
               @assert value1 > value2 "$(value1) is not > $(value2)!"
          elseif operator == ">="
               @assert value1 >= value2 "$(value1) is not >= $(value2)!"
          elseif operator == "approx"
               @assert value1 ≈ value2 "$(value1) is not ≈ $(value2)!"
          else
               throw(GrammarError(token.location, "operator $(operator)  is not valid"))
          end
     else
          unread_token(inputstream, token)
          expect_symbol(inputstream, ")")
          @assert value1 == value2 "$(value1) and $(value2) are not equal!"
     end
end

ext: jl


text:
module Encoders

export AbstractEncoder, encode

abstract type AbstractEncoder end

function encode! end
function encode end

########################### Generic encoding methods ###########################
"""
    encode(encoder::E, x::T) where {E <: AbstractEncoder,T <: Real}

Encode `x` as a binary vector according to `encoder`.

# Examples
```jldoctest
julia> encode(Thermometer(Float64, 0, 1, 5), 0.5)
5-element Array{Bool,1}:
 1
 1
 1
 0
 0

 julia> encode(CircularThermometer(Float64, 0, 1, 5), 0.5)
 5-element Array{Bool,1}:
 0
 0
 1
 1
 0
```
"""
function encode(encoder::E, x::T) where {E <: AbstractEncoder,T <: Real}
    pattern = Vector{Bool}(undef, encoder.resolution)

    encode!(encoder, x, pattern)
end

"""
    encode!(encoder::E, x::AbstractVector{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}

Encode `x` as a binary matrix, where each component of `x` is mapped into a
column, according to `encoder` and place the result in `pattern`.

# Examples
```jldoctest
julia> encode!(Thermometer(Float64, 0, 1, 5), [0.1, 0.5], Matrix{Bool}(undef, 5, 2))
5×2 Array{Bool,2}:
 1  1
 0  1
 0  1
 0  0
 0  0

julia> encode!(CircularThermometer(Float64, 0, 1, 5), [0.1, 0.5], Matrix{Bool}(undef, 5, 2))
5×2 Array{Bool,2}:
 1  0
 1  0
 0  1
 0  1
 0  0
```
"""
function encode!(encoder::E, x::AbstractVector{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}
    validate_input_output(encoder, x, pattern)

    if encoder.min isa AbstractVector
        for (i, col) in Iterators.enumerate(eachcol(pattern))
            encode!(encoder, encoder.min[i], encoder.max[i], x[i], col)
        end
    else
        for (i, col) in Iterators.enumerate(eachcol(pattern))
            encode!(encoder, x[i], col)
        end
    end

    pattern
end

"""
    encode!(encoder::E, x::AbstractVector{T}, pattern::AbstractVector{Bool}) where {E <: AbstractEncoder,T <: Real}

Encode `x` as a binary vector, where each component of `x` is mapped into a
partition of the binary vector, according to `encoder` and place the result in
`pattern`.

# Examples
```jldoctest
julia> encode!(Thermometer(Float64, 0, 1, 5), [0.1, 0.5], Vector{Bool}(undef, 10))
10-element Array{Bool,1}:
 1
 0
 0
 0
 0
 1
 1
 1
 0
 0

julia> encode!(CircularThermometer(Float64, 0, 1, 5), [0.1, 0.5], Vector{Bool}(undef, 10))
10-element Array{Bool,1}:
 1
 1
 0
 0
 0
 0
 0
 1
 1
 0
```
"""
function encode!(encoder::E, x::AbstractVector{T}, pattern::AbstractVector{Bool}) where {E <: AbstractEncoder,T <: Real}
    validate_input_output(encoder, x, pattern)

    if encoder.min isa AbstractVector
        for (i, slice) in Iterators.enumerate(Iterators.partition(pattern, encoder.resolution))
            encode!(encoder, encoder.min[i], encoder.max[i], x[i], slice)
        end
    else
        for (i, slice) in Iterators.enumerate(Iterators.partition(pattern, encoder.resolution))
            encode!(encoder, x[i], slice)
        end
    end

    pattern
end

"""
    encode(encoder::E, x::AbstractVector{T}; flat=true) where {E <: AbstractEncoder,T <: Real}

Encode `x` according to `encoder`.

If keyword argument `flat` is `true`, the result is as a binary vector, where
each component of `x` is mapped into a partition of the output vector. If `flat`
is false, the result is a binary matrix, where each component is mapped to a
column of the output matrix.

# Examples
```jldoctest
julia> encode(Thermometer(Float64, 0, 1, 5), [0.1, 0.5]; flat=true)
10-element Array{Bool,1}:
 1
 0
 0
 0
 0
 1
 1
 1
 0
 0

 julia> encode(Thermometer(Float64, 0, 1, 5), [0.1, 0.5]; flat=false)
 5×2 Array{Bool,2}:
  1  1
  0  1
  0  1
  0  0
  0  0
```
"""
function encode(encoder::E, x::AbstractVector{T}; flat=true) where {E <: AbstractEncoder,T <: Real}
    if flat
        pattern = Vector{Bool}(undef, encoder.resolution * length(x))
    else
        pattern = Matrix{Bool}(undef, encoder.resolution, length(x))
    end

    encode!(encoder, x, pattern)
end

"""
    encode!(encoder::E, X::AbstractMatrix{T}, pattern::AbstractArray{Bool,3}) where {E <: AbstractEncoder,T <: Real}

Encode each row of `X` according to `encoder`.

Each row of `X` is encoded as a binary matrix and mapped to a slice of `pattern`.

# Examples
```jldoctest
julia> encode!(Thermometer(Float64, 0, 1, 5), [0 0.1; 0.3 0.7], Array{Bool,3}(undef, 5, 2, 2))
5×2×2 Array{Bool,3}:
[:, :, 1] =
 0  1
 0  0
 0  0
 0  0
 0  0

[:, :, 2] =
 1  1
 1  1
 0  1
 0  1
 0  0

```
"""
function encode!(encoder::E, X::AbstractMatrix{T}, pattern::AbstractArray{Bool,3}) where {E <: AbstractEncoder,T <: Real}
    validate_input_output(encoder, X, pattern)

    for (slice, row) in Iterators.zip(eachslice(pattern; dims=3), eachrow(X))
        encode!(encoder, row, slice)
    end

    pattern
end

# Flat version. All values in each row of the input matrix are encoded and placed in a single row of the output matrix
"""
    encode!(encoder::E, X::AbstractMatrix{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}

Encode `X` according to `encoder`.

Each row of `X` is encoded as a binary vector and mapped to a row of `pattern`.

# Examples
```jldoctest
julia> encode!(Thermometer(Float64, 0, 1, 5), [0 0.1; 0.3 0.7], Matrix{Bool}(undef, 2, 10))
2×10 Array{Bool,2}:
 0  0  0  0  0  1  0  0  0  0
 1  1  0  0  0  1  1  1  1  0

```
"""
function encode!(encoder::E, X::AbstractMatrix{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}
    validate_input_output(encoder, X, pattern)

    for (slice, row) in Iterators.zip(eachrow(pattern), eachrow(X))
        encode!(encoder, row, slice)
    end

    pattern
end

"""
    encode(encoder::E, X::AbstractMatrix{T}; flat=true) where {E <: AbstractEncoder,T <: Real}

Encode `x` according to `encoder`.

If keyword argument `flat` is `true`, the result is as a binary matrix, where
each component of row `X` is mapped into a row of the output matrix. If `flat`
is false, the result is a 3-dimensional binary array, where each row is mapped
to a slice of the output array.

# Examples
```jldoctest
julia> encode(Thermometer(Float64, 0, 1, 5), [0 0.1; 0.3 0.7]; flat=true)
2×10 Array{Bool,2}:
 0  0  0  0  0  1  0  0  0  0
 1  1  0  0  0  1  1  1  1  0

julia> encode(Thermometer(Float64, 0, 1, 5), [0 0.1; 0.3 0.7]; flat=false)
5×2×2 Array{Bool,3}:
[:, :, 1] =
 0  1
 0  0
 0  0
 0  0
 0  0

[:, :, 2] =
 1  1
 1  1
 0  1
 0  1
 0  0
```
"""
function encode(encoder::E, X::AbstractMatrix{T}; flat=true) where {E <: AbstractEncoder,T <: Real}
    if flat
        pattern = Matrix{Bool}(undef, size(X, 1), encoder.resolution * size(X, 2))
    else
        pattern = Array{Bool,3}(undef, encoder.resolution, size(X, 2), size(X, 1))
    end

    encode!(encoder, X, pattern)
end

########################## Generic validation methods ##########################
function validate_input_output(encoder::E, x::T, pattern::AbstractVector{Bool}) where {E <: AbstractEncoder,T <: Real}
    # If more then one minimum/maximum were especified, it doesn't make sense to try to encode a scalar
    length(encoder.min) != 1 && throw(DomainError(x, "expected vector of length $(length(encoder.min))"))

    length(pattern) != encoder.resolution && throw(
      DimensionMismatch(
        "output pattern must have the length of the thermometer's resolution"))

    nothing
end

function validate_input_output(encoder::E, x::AbstractVector{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}
    size(pattern, 1) != encoder.resolution && throw(
      DimensionMismatch(
        "the number of rows of the output pattern must equal the thermometer's resolution"))

    size(pattern, 2) != length(x) && throw(
      DimensionMismatch(
          "the number of columns of the output pattern must equal x's length"))

    encoder.min isa AbstractVector && length(x) != length(encoder.min) && throw(
      DimensionMismatch(
        "the length of x must match the lengths of the min and max vectors"))

    nothing
end

function validate_input_output(encoder::E, x::AbstractVector{T}, pattern::AbstractVector{Bool}) where {E <: AbstractEncoder,T <: Real}
    length(pattern) != length(x) * encoder.resolution && throw(
      DimensionMismatch(
        "the length of the output pattern must equal the thermometer's resolution times the length of x"))

    encoder.min isa AbstractVector && length(x) != length(encoder.min) && throw(
      DimensionMismatch(
        "the length of x must match the lengths of the min and max vectors"))
    nothing
end

function validate_input_output(encoder::E, X::AbstractMatrix{T}, pattern::AbstractArray{Bool,3}) where {E <: AbstractEncoder,T <: Real}
    size(pattern, 1) != encoder.resolution && throw(
      DimensionMismatch(
        "the number of rows of the output pattern must equal the thermometer's resolution"))

    size(pattern, 2) != size(X, 2) && throw(
      DimensionMismatch(
          "the number of columns of the output pattern must equal X's number of columns"))

    size(pattern, 3) != size(X, 1) && throw(
      DimensionMismatch(
          "the depth of the output pattern must equal X's number of rows"))

    encoder.min isa AbstractVector && size(X, 2) != length(encoder.min) && throw(
      DimensionMismatch(
        "the number of columns of X must match the lengths of the min and max vectors"))

    nothing
end

function validate_input_output(encoder::E, X::AbstractMatrix{T}, pattern::AbstractMatrix{Bool}) where {E <: AbstractEncoder,T <: Real}
    size(pattern, 1) != size(X, 1) && throw(
        DimensionMismatch(
          "the number of rows of the output pattern must equal the the number of rows of X"))

    size(pattern, 2) != size(X, 2) * encoder.resolution && throw(
          DimensionMismatch(
              "the number of columns of the output pattern must equal the thermometer's resolution times X's number of columns"))

    encoder.min isa AbstractVector && size(X, 2) != length(encoder.min) && throw(
      DimensionMismatch(
        "the number of columns of X must match the lengths of the min and max vectors"))

    nothing
end

include("Thermometer.jl")
export Thermometer

include("CircularThermometer.jl")
export CircularThermometer

end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
"""
`Basics.integrate()`
    ... integrates a one- or two-dimensional function by different numerical methods, either on a given grid or for a general function.

  + `("function: on radial grid, Newton-Cotes", F::Array{Float64,1}, grid::Radial.Grid)`
    ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r) by using a 5-point Newton-Cotes
        integration formula; a value::Float64 is returned.

  + `("function: on radial grid, Simpson rule", F::Array{Float64,1}, grid::Radial.Grid)`
    ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r) by using Simpson's rule;
        a value::Float64 is returned.

  + `("function: on radial grid, trapez rule", F::Array{Float64,1}, grid::Radial.Grid)`
    ... to integrate the (radial) function F over the given (radial) grid, int_0^infinity dr F(r), by using a simple trapez rule;
        a value::Float64 is returned.
"""
function Basics.integrate(sa::String, F::Array{Float64,1}, grid::Radial.Grid)
    if       sa == "function: on radial grid, Newton-Cotes"     wa = integrateOnGridNewtonCotes(F, grid)
    elseif   sa == "function: on radial grid, Simpson rule"     wa = integrateOnGridSimpsonRule(F, grid)
    elseif   sa == "function: on radial grid, trapez rule"      wa = integrateOnGridTrapezRule( F, grid)
    else     error("Unsupported keystring = $sa.")
    end

end

"""
`Basics.integrateOnGridNewtonCotes()`  ... integrates by using a 5-point Newton-Cotes integration formula; see Basic.integrate().
"""
function Basics.integrateOnGridNewtonCotes(F::Array{Float64,1}, grid::Radial.Grid)
    coefficients = Array{Float64}([2 * 7, 32, 12, 32, 7]) * 2 / 45 * grid.h
    n = size(coefficients, 1)

    result = 0.;    i0     = 1

    while  abs(F[i0]) == 0
        i0 += 1
    end
    gamma  = log(F[i0]/F[i0+1] * grid.rp[i0+1]/grid.rp[i0]) / log(grid.r[i0]/grid.r[i0+1])
    result = 1/(gamma + 1) * grid.r[i0] * F[i0] / grid.rp[i0]

    result += F[i0] * coefficients[n]
    for i = i0+1:size(F, 1)
        result += coefficients[(i - i0) % (n - 1) + 1] * F[i]
    end

    if (size(F, 1) - i0) % (n - 1) + 1 == 1
        result -= coefficients[n] * F[size(F, 1)]
    end

    return( result )
end

"""
`Basics.integrateOnGridSimpsonRule()`  ... integrates by using Simpson's rule; see JAC.Basics.integrate().
"""
function Basics.integrateOnGridSimpsonRule(F::Array{Float64,1}, grid::Radial.Grid)
    coefficients = Array{Float64}([2 * 1, 4, 1]) / 3 * grid.h
    n = size(coefficients, 1)

    result = 0.
    result += F[1] * coefficients[n]
    for i = 2:size(F, 1)
        result += coefficients[(i - 1) % (n - 1) + 1] * F[i]
    end

    if (size(F, 1) - 1) % (n - 1) + 1 == 1
        result -= coefficients[n] * F[size(F, 1)]
    end

    return( result )
end

"""
`Basics.integrateOnGridTrapezRule()`  ... integrates by using a simple trapez rule; see JAC.Basics.integrate().
"""
function Basics.integrateOnGridTrapezRule(F::Array{Float64,1}, grid::Radial.Grid)
    coefficients = Array{Float64}([2 * 1, 1]) / 2 * grid.h
    n = size(coefficients, 1)

    result = 0.

    result += F[1] * coefficients[n]
    for i = 2:size(F, 1)
        result += coefficients[(i - 1) % (n - 1) + 1] * F[i]
    end

    if (size(F, 1) - 1) % (n - 1) + 1 == 1
        result -= coefficients[n] * F[size(F, 1)]
    end

    return( result )
end


text:
@test "add_test_files skips . and .." begin
  @assert !in(".", Testfast.files)
  @assert !in("..", Testfast.files)
end

@test "add_test_files scans through all directories and adds all files" begin
  @assert in("./folder/nested_folder/testfast_nested_folder_test.jl", Testfast.files)
  @assert in("./folder/testfast_folder_test.jl", Testfast.files)
  @assert in("./testfast_test.jl", Testfast.files)
end

@test "before_all is added properly" begin
  @assert before_all_test == 1
end

@test "setup is added properly" begin
  @assert setup_test == 4 # This is the fourth run test
end

@test "teardown is added properly" begin
  @assert teardown_test == 4 # This is the fifth test, so teardowns have run four times
end

@test "after_all is added properly" begin
  @assert after_all_test == 0 # This hasn't run at all yet
end

@test "no leaking from setups" begin
  error = false

  try
    this_should_not_leak == true
  catch
    error = true
  end

  @assert error

end

@test "no leaking from teardowns" begin
  error = false

  try
    this_also_should_not_leak == true
  catch
    error = true
  end

  @assert error

end


text:
using Pipe
using DataFrames
using Statistics
using JLD2
using CategoricalArrays
using DocStringExtensions

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct Dataset
    X::Array{Union{Missing,Float64},2}
    means::Array{Float64,1}
    stdevs::Array{Float64,1}
    value_columns::Array{String,1}
    xmask::BitArray{2}
    mv::Bool
    mvs::Array{Float64,1}
    varvalues::Array{Int64,1}
    ranges::Array{Float64,1}
end

abstract type MultivariateModel end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct PCA <: MultivariateModel
    T::DataFrame
    P::DataFrame
    r2x_cum::Vector{Float64}
    r2x::Vector{Float64}
    eigenvalues::Vector{Float64}
end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct PLS <: MultivariateModel
    T::DataFrame
    P::DataFrame
    C::DataFrame
    W::DataFrame
    U::DataFrame
end

function copydataset(dataset)::Dataset
    Dataset(
        copy(dataset.X),
        copy(dataset.means),
        copy(dataset.stdevs),
        copy(dataset.value_columns),
        copy(dataset.xmask),
        dataset.mv,
        copy(dataset.mvs),
        copy(dataset.varvalues),
        copy(dataset.ranges)
    )
end

function copydataset(dataset,obsmask::BitArray{1}, varmask::BitArray{1})::Dataset
    Dataset(
        dataset.X[obsmask,varmask],
        dataset.means[varmask],
        dataset.stdevs[varmask],
        dataset.value_columns[varmask],
        dataset.xmask[obsmask,varmask],
        dataset.mv,
        dataset.mvs[varmask],
        dataset.varvalues[varmask],
        dataset.ranges[varmask]
    )
end

function fillmissings(dataset::Dataset; fillvalue = 0.0)
    dataset.X[.~dataset.xmask] .= fillvalue
end

function hasVariation(x)
    !isnan(x) && x > 0
end

function norm(v)
    sqrt(sum(v.^2))
end

function valueRange(array)
    valrng = NaN

    try
        valrng = maximum(skipmissing(array))-minimum(skipmissing(array))
    catch
        valrng = NaN
    end

    valrng
end

isnumcol(type)= type <: Union{Missing,Number}

function parseMatrix(X::Array{Union{Missing, Float64},2},value_columns::Array{String,1})

    ranges = [valueRange(skipmissing(col)) for col in eachcol(X) ]

    var_means::Array{Float64,1} = [mean(skipmissing(col)) for col in eachcol(X) ]
    mean_mask = (!isnan).(var_means)

    var_stdevs::Array{Float64,1} = [std(skipmissing(col)) for col in eachcol(X) ]
    std_mask = hasVariation.(var_stdevs)

    inc_mask = mean_mask .& std_mask

    #Xtr = X[:,inc_mask]
    Xtr = X[:,:]

    xmask = (!ismissing).(Xtr)

    varvalues = sum(xmask,dims=1)[:]

    mvs = 1 .- (varvalues ./ size(X,1))

    Dataset(Xtr, var_means, var_stdevs, value_columns, xmask, sum(.~xmask) > 0, mvs, varvalues,ranges)
end

function parseDataFrame(X::Matrix{Union{Missing,Float64}},value_columns::Vector{String}; filters=[dataset-> .~isnan.(dataset.means)])
    return parseMatrix(X,value_columns) |> ds -> filterDataset(ds,filters=filters)
end

function parseDataFrame(X::Matrix{Union{Missing,Float64}})

    value_columns = collect(range(1,length=size(X,2))) .|> string .|> i -> "var"*i

    return parseDataFrame(X,value_columns)

end

function parseDataFrame(df::AbstractDataFrame; filters=[dataset-> .~isnan.(dataset.means)])::Dataset

    value_columns = names(df[:,eltype.(eachcol(df)) |> coltypes -> isnumcol.(coltypes)])

    X = Matrix{Union{Missing,Float64}}(df[:,value_columns])

    return parseMatrix(X,value_columns) |> ds -> filterDataset(ds,filters=filters)
end

function calcObsMissingValues(dataset::Dataset)
    calcMissingValues(dataset.xmask,2)
end

function calcObsMissingValues(mask::BitArray{2})
    calcMissingValues(mask,2)
end

function calcVarMissingValues(mask::BitArray{2})
    calcMissingValues(mask,1)
end

function calcMissingValues(mask::BitArray{2}, dim)
    sum(.~mask,dims=dim) ./ size(mask,dim)
end

function filterDataset(dataset::Dataset; filters=[], obsfilters=[])::Dataset
    mask = filterVariables(dataset,filters=filters)

    obsmask = filterObservations(dataset,filters=obsfilters)

    Dataset(
        dataset.X[obsmask,mask],
        dataset.means[mask],
        dataset.stdevs[mask],
        dataset.value_columns[mask],
        dataset.xmask[obsmask,mask],
        sum(.~(dataset.xmask[obsmask,mask])) > 0,
        dataset.mvs[mask],
        dataset.varvalues[mask],
        dataset.ranges[mask]
    )
end

function filterVariables(dataset::Dataset; filters=[])

    allvariables = allvars(dataset)

    varmask = @pipe [filter(dataset) for filter in filters] |> [allvariables,_...] |> tomatrix |> m->all(m,dims=2) |> vec
end

function filterObservations(dataset::Dataset; filters=[])

    allobservations = allobs(dataset)

    varmask = @pipe [filter(dataset) for filter in filters] |> [allobservations,_...] |> tomatrix |> m->all(m,dims=2) |> vec
end

function nobs(dataset::Dataset)
    size(dataset.X,1)
end

function nvars(dataset::Dataset)
    size(dataset.X,2)
end

function tomatrix(arrayOfArrays)
    hcat(arrayOfArrays...)
end

function allvars(dataset::Dataset)::BitArray{1}
    BitArray{1}(ones(nvars(dataset)))
end

function allobs(dataset::Dataset)::BitArray{1}
    BitArray{1}(ones(nobs(dataset)))
end

function normalize!(dataset::Dataset; doscale::Bool=false, stdevs=dataset.stdevs, means=dataset.means, transformations=[])

    #TODO:double check that it works
    mean_mask = (!isnan).(means)
    std_mask = hasVariation.(stdevs)

    inc_mask = mean_mask .& std_mask

    # has variables to be excluded
    excludevars = sum(inc_mask) < nvars(dataset)

    if !isempty(transformations)
        [transform(dataset) transform in transformations]
    else
        dataset.X[:,mean_mask] .-= dataset.means[mean_mask]'

        if doscale
            dataset.X[:,inc_mask] ./= dataset.stdevs[inc_mask]'
        end
    end

    dataset.X[.~dataset.xmask] .= 0

    dataset
end

function normalize(dataset::Dataset; doscale::Bool=false, stdevs=dataset.stdevs, means=dataset.means)

    dscopy = copydataset(dataset)

    normalize!(dscopy,doscale=doscale)
end

function normalizedata(X::Array{Union{Missing,Float64},2};normalize::Bool=false)
    var_means::Array{Float64,1} = [mean(skipmissing(col)) for col in eachcol(X) ]
    mean_mask = (!isnan).(var_means)

    var_stdevs::Array{Float64,1} = [std(skipmissing(col)) for col in eachcol(X) ]
    std_mask = hasVariation.(var_stdevs)

    inc_mask = mean_mask .& std_mask

    Xtr = X[:,inc_mask]

    xmask = (!ismissing).(Xtr)

    Xtr .-= var_means[inc_mask]'

    if normalize
        Xtr ./= var_stdevs[inc_mask]'
    end

    Xtr[.~xmask] .= 0
end

"""
$(FUNCTIONNAME)(model::T, dataset::Dataset, name::String) where T <: MultivariateModel

    Save PCA or PLS model as JLD2 file
"""
function savemodel(model::T, dataset::Dataset, path::String, transformations::Array{String,1}) where T <: MultivariateModel

    values = fieldnames(T) |> fns -> getfield.(Ref(model), fns)

    zipped = zip(string.(fieldnames(T)), values)

    jldopen(path, "w") do file

        file["modeltype"] = string(T)

        foreach(fv -> file[fv[1]] = fv[2], zipped)

        file["means"] = DataFrame(zip(dataset.value_columns, dataset.means)) |> df -> rename!(df, [:var,:mean])

        file["stdevs"] = DataFrame(zip(dataset.value_columns, dataset.stdevs)) |> df -> rename!(df, [:var,:stdev])

        file["transformations"] = transformations

        close(file)
    end
end

"""
$(FUNCTIONNAME)(path::String)::Tuple{MultivariateModel,Array{Float64,1},Array{Float64,1}}

    Load PCA or PLS model from JLD2 file into a tuple containing the model, variable standard deviations from variable means from model calibration

"""
function loadmodel(modelfile::String)::Tuple{MultivariateModel,Array{Float64,1},Array{Float64,1},Array{String,1},Array{String,1}}

    jldfile = jldopen(modelfile, "r")

    modeltype = jldfile["modeltype"]
    stdevs = jldfile["stdevs"] |> df -> convert(Array{Float64,2}, df[:,[:stdev]])[:]
    means = jldfile["means"] |> df -> convert(Array{Float64,2}, df[:,[:mean]])[:]

    transformations = jldfile["transformations"]

    variables = jldfile["means"][:,:var]

    type::DataType =  modeltype == "PCA" ? PCA : PLS

    values = map(fname -> jldfile[string(fname)], fieldnames(type))

    close(jldfile)

    type(values...),stdevs,means,variables,transformations
end

function selectNumerical(df)

    isnumcol(type) = type <: Union{Missing,Number}

    select(df, eltype.(eachcol(df)) |> (coltypes -> isnumcol.(coltypes)) |> findall)
end

function selectColumns(df,typesel)

    select(df, eltype.(eachcol(df)) |> typesel |> findall)
end

```
Get a onehot representation of a CategoricalArray
```
function onehot(carray::CategoricalArray)
    lvls = carray |> levels

    hots = [(carray .== lv) |> bvals -> convert(Array{Union{Missing,Float64},1},bvals) for lv in lvls]

    DataFrame(hots) |> v -> rename!(v,Symbol.(lvls))
end

function predictLevel(carray::CategoricalArray,predMatrix::Array{Union{Missing, Float64},2})::Array{String,1}
    predicted_classes = [findmax(row) |> last for row in eachrow(predMatrix)] |> maxidxs -> getindex(levels(carray),maxidxs)
end

function getTransformations(model::T) where T <: MultivariateModel
    return []
end


text:
###############################################################################
##### Predict functions
###############################################################################

function getμ(tree::Tree)
  [leaf.μ for leaf in leafnodes(tree.root)]
end

function StatsBase.predict(bt::BartTree)
  bt.S * getμ(bt.tree)
end

function StatsBase.predict(bs::BartState, bm::BartModel)
  fhat = zeros(bm.td.n)
  for bt in bs.ensemble.trees
    fhat += predict(bt)
  end
  fhat
end

function StatsBase.predict(trees::Vector{Tree}, X::Matrix{Float64})
  yhat = zeros(size(X, 1))
  for tree in trees
    yhat += leafprob(X, tree) * getμ(tree)
  end
  yhat
end

function StatsBase.predict(bc::RegBartChain, X::Matrix{Float64})
  treedraws = reshape(bc.treedraws, size(bc.treedraws, 1)*size(bc.treedraws, 3))
  reduce(hcat, pmap(t -> predict(t, X), treedraws)) .+ bc.bm.td.ybar
end

function StatsBase.predict(bc::ProbitBartChain, X::Matrix{Float64})
  treedraws = reshape(bc.treedraws, size(bc.treedraws, 1)*size(bc.treedraws, 3))
  reduce(hcat, pmap(t -> cdf.(Normal(), predict(t, X)), treedraws))
end


text:
#=
using BenchmarkTools
using ConstrainedRootSolvers

function benchmark_ConstrainedRootSolvers(FT)
    # A S-shape function, e.g. Weibull function
    @inline function _r_func(x::FT) where {FT<:AbstractFloat}
        return exp( -1 * (x/2) ^ FT(0.8) ) - FT(0.5)
    end

    # A R-shape function
    @inline function _s_func(x::FT) where {FT<:AbstractFloat}
        return exp( -1 * (x/2) ^ 5 ) - FT(0.5)
    end

    # segmented function with solution 1
    @inline function _f0(x::FT) where {FT<:AbstractFloat}
        return x>1 ? FT(1) : x
    end

    # segmented function with solution 1
    @inline function _f1(x::FT) where {FT<:AbstractFloat}
        return x>1 ? 2-x : x
    end

    # segmented function with solution of upper end
    @inline function _f2(x::FT) where {FT<:AbstractFloat}
        return x>1 ? (1+x)/2 : x
    end

    # segmented function with solution 1
    @inline function _f3(x::FT) where {FT<:AbstractFloat}
        return x>1 ? 2-x : FT(0)
    end

    # segmented function with solution of lower bound
    @inline function _f4(x::FT) where {FT<:AbstractFloat}
        return x>1 ? 2-x : FT(1)
    end

    # segmented function with solution 1
    @inline function _f5(x::FT) where {FT<:AbstractFloat}
        return x>1 ? FT(0.5) : x
    end

    # surface function
    @inline function _surf_func(x::Array{FT,1}) where {FT<:AbstractFloat}
        return -(x[1]-2)^2 - (x[2]-3)^2
    end

    println("\nBenchmarking find_peak with Bisection method...");
    ms = BisectionMethod{FT}(x_min=0, x_max=5);
    st = SolutionTolerance{FT}(1e-3, 50);
    lr = [1, 1, 5, 1, 0, 1]
    lf = [_f0, _f1, _f2, _f3, _f4, _f5]
    for f in lf
        @btime find_peak($f, $ms, $st);
    end

    println("\nBenchmarking find_peak with NelderMeadMethod method...");
    ms = NelderMeadMethod{FT}(N=2, x_inis=FT[rand(FT)+1, rand(FT)+2, 0]);
    rt = ResidualTolerance{FT}(1e-5, 50);
    st = SolutionToleranceND{FT}(FT[1e-3, 1e-3], 50);
    @btime find_peak($_surf_func, $ms, $rt);
    @btime find_peak($_surf_func, $ms, $st);

    println("\nBenchmarking find_peak with ReduceStepMethod method...");
    ms = ReduceStepMethod{FT}(0, 5, 1+rand(FT), 1);
    st = SolutionTolerance{FT}(1e-3, 50);
    lr = [1, 1, 5, 1, 0, 1]
    lf = [_f0, _f1, _f2, _f3, _f4, _f5]
    for f in lf
        @btime find_peak($f, $ms, $st);
    end

    println("\nBenchmarking find_peak with ReduceStepMethodND method...");
    ms = ReduceStepMethodND{FT}(x_mins = FT[0,0],
                                x_maxs = FT[5,5],
                                x_inis = [1+rand(FT), 2+rand(FT)],
                                Δ_inis = FT[1,1]);
    st = SolutionToleranceND{FT}(FT[1e-3, 1e-3], 50);
    @btime find_peak($_surf_func, $ms, $st);

    println("\nBenchmarking find_zero with BisectionMethod method...");
    ms = BisectionMethod{FT}(x_min=0, x_max=5);
    rt = ResidualTolerance{FT}(1e-3, 50);
    st = SolutionTolerance{FT}(1e-3, 50);
    @btime find_zero($_r_func, $ms, $rt);
    @btime find_zero($_r_func, $ms, $st);
    @btime find_zero($_s_func, $ms, $rt);
    @btime find_zero($_s_func, $ms, $st);

    println("\nBenchmarking find_zero with NewtonBisectionMethod method...");
    ms = NewtonBisectionMethod{FT}(0, 5, 1+rand(FT));
    rt = ResidualTolerance{FT}(1e-3, 50);
    st = SolutionTolerance{FT}(1e-3, 50);
    @btime find_zero($_r_func, $ms, $rt);
    @btime find_zero($_r_func, $ms, $st);
    @btime find_zero($_s_func, $ms, $rt);
    @btime find_zero($_s_func, $ms, $st);

    println("\nBenchmarking find_zero with NewtonRaphsonMethod method...");
    ms = NewtonRaphsonMethod{FT}(2);
    rt = ResidualTolerance{FT}(1e-3, 50);
    st = SolutionTolerance{FT}(1e-3, 50);
    @btime find_zero($_r_func, $ms, $rt);
    @btime find_zero($_r_func, $ms, $st);
    @btime find_zero($_s_func, $ms, $rt);
    @btime find_zero($_s_func, $ms, $st);

    println("\nBenchmarking find_zero with ReduceStepMethod method...");
    ms = ReduceStepMethod{FT}(0, 5, 1+rand(FT), 1);
    rt = ResidualTolerance{FT}(1e-3, 50);
    st = SolutionTolerance{FT}(1e-3, 50);
    @btime find_zero($_r_func, $ms, $rt);
    @btime find_zero($_r_func, $ms, $st);
    @btime find_zero($_s_func, $ms, $rt);
    @btime find_zero($_s_func, $ms, $st);

    return nothing
end

benchmark_ConstrainedRootSolvers(Float32);
benchmark_ConstrainedRootSolvers(Float64);
=#


text:
module SimpleDrawing
using Plots, LinearAlgebra

export newdraw, draw_circle, draw_arc, draw_segment, draw_point, draw_vector
export draw_rectangle, find_center, non_colinear_check
export finish, draw

"""
`newdraw()` clears the screen and sets up a blank drawing canvas with
no axis, grid, or legend, and sets the aspect ratio to 1.
"""
function newdraw()
    plot(aspectratio = 1, legend = false, axis = false, grid = false, ticks = false)
end

function draw() end

"""
    draw_circle(x,y,r)

Draw a circle centered
at `(x,y)` with radius `r`.

Also: `draw_circle(z,r)` where `z` is complex.
"""
function draw_circle(x::Real, y::Real, r::Real; opts...)
    f(t) = r * cos(t) + x
    g(t) = r * sin(t) + y
    plot!(f, g, 0, 2pi; opts...)
end

draw_circle(z::Complex, r::Real; opts...) = draw_circle(real(z), imag(z), r; opts...)

"""
    draw_circle(a,b,c)

Draw a circle through the points given by
the three complex arguments.
"""
function draw_circle(a::Complex, b::Complex, c::Complex; opts...)
    z = find_center(a, b, c)
    r = abs(a - z)
    x, y = reim(z)
    draw_circle(x, y, r; opts...)
end

export draw_disc

"""
    draw_disc(x,y,r)

Draw a disc centered at `(x,y)` and with radius `r`. May also be invoked as:
* `draw_disc(z,r)` where `z` is complex.
* `draw_disc(a,b,c)` where `a,b,c` are complex which gives a disc passing through `a`, `b`, and `c`.

For example: `draw_disc(0,1,3; color=:yellow, linecolor=:red)`
"""
function draw_disc(x::Real, y::Real, r::Real; opts...)
    draw_circle(x, y, r, seriestype = [:shape], fillalpha = 1; opts...)
end

function draw_disc(z::Complex, r::Real; opts...)
    draw_disc(real(z), imag(z), r; opts...)
end

function draw_disc(a::Complex, b::Complex, c::Complex; opts...)
    z = find_center(a, b, c)
    r = abs(a - z)
    x, y = reim(z)
    draw_disc(x, y, r; opts...)
end

"""
`draw_arc(x,y,r,t1,t2)` draws the shorter arc centered at `(x,y)`
with radius `r` between angles `t1` and `t2`.
"""
function draw_arc(x::Real, y::Real, r::Real, t1::Real, t2::Real; opts...)
    f(t) = r * cos(t) + x
    g(t) = r * sin(t) + y

    plot!(f, g, t1, t2; opts...)
end

"""
`draw_arc(a,b,c)` where the arguments are `Complex`: Draw the
arc from `a` through `b` to `c`.
"""
function draw_arc(a::Complex, b::Complex, c::Complex; opts...)
    m2pi(x::Real) = mod(x, 2π)

    z = find_center(a, b, c)
    x, y = reim(z)
    r = abs(b - z)
    ta = m2pi(angle(a - z))
    tb = m2pi(angle(b - z))
    tc = m2pi(angle(c - z))

    if ta > tc
        ta, tc = tc, ta
    end

    if ta < tb < tc
        return draw_arc(x, y, r, ta, tc; opts...)
    end

    return draw_arc(x, y, r, tc - 2pi, ta; opts...)
end

"""
`draw_segment(x,y,xx,yy)` draws a line segment from `(x,y)` to
`(xx,yy)`. Also `draw_segment(z,zz)` for `Complex` arguments.
"""
function draw_segment(a::Real, b::Real, c::Real, d::Real; opts...)
    plot!([a, c], [b, d]; opts...)
end

function draw_segment(a::Complex, b::Complex; opts...)
    x, y = reim(a)
    xx, yy = reim(b)
    draw_segment(x, y, xx, yy; opts...)
end

"""
`draw_point(x,y)` [or `draw_point(z)`] draws a point at coordinates `(x,y)`
[or complex location `z`].

`draw_point(list)` plots all the points in the one-dimensional `list` of
complex values.
"""
function draw_point(x::Real, y::Real; opts...)
    plot!([x], [y]; marker = 1, opts...)
end

draw_point(z::Complex; opts...) = draw_point(real(z), imag(z); opts...)

function draw_point(pts::Array{Complex{T},1}; opts...) where {T}
    for p in pts
        draw_point(p; opts...)
    end
    plot!()
end

"""
`draw_rectangle(x,y,xx,yy)` draws a rectangle with opposite corners
`(x,y)` and `(xx,yy)`.  May also be called with complex arguments
`draw_rectangle(w,z)`.
"""
function draw_rectangle(x::Real, y::Real, xx::Real, yy::Real; opts...)
    xlist = [x, xx, xx, x, x]
    ylist = [y, y, yy, yy, y]
    plot!(xlist, ylist; opts...)
end

function draw_rectangle(w::Complex, z::Complex; opts...)
    draw_rectangle(real(w), imag(w), real(z), imag(z); opts...)
end

"""
`draw_vector` is used to draw vectors (line segments with an arrow at one end).
The variations are:
+ `draw_vector(x,y)` draws a vector from the origin to `(x,y)`.
+ `draw_vector(x,y,basex,basey)` draws a vector from `(basex,basey)` to
`(basex+x,basey+y)`.
+ `draw_vector(z)` draws a vector from the origin to the complex value `z`.
+ `draw_vector(z,basez)` draws a vector from the complex location `basez` to
`z+basez`.
"""
draw_vector(x::Real, y::Real; opts...) = draw_segment(0, 0, x, y; arrow = :arrow, opts...)

draw_vector(z::Complex; opts...) = draw_segment(0 + 0im, z; arrow = :arrow, opts...)

draw_vector(x::Real, y::Real, basex::Real, basey::Real; opts...) =
    draw_segment(basex, basey, basex + x, basey + y; arrow = :arrow, opts...)

draw_vector(z::Complex, basez::Complex; opts...) =
    draw_segment(basez, basez + z; arrow = :arrow, opts...)

"""
`finish()` is used to clean up a drawing after various calls to `draw`.
It removes the axes and the grid lines, and sets the aspect ratio to one.
"""
function finish()
    plot!(aspectratio = 1, legend = false, axis = false, grid = false, ticks = false)
end

export find_center, non_colinear_check

"""
`find_center(a,b,c)`: Given three points in the complex plane,
find the point `z` that is equidistant from all three. If the three
points are collinear then return `Inf + Inf*im`.
"""
function find_center(a::Complex, b::Complex, c::Complex)::Complex
    if !non_colinear_check(a, b, c)
        return Inf + im * Inf
    end
    A = collect(reim(a))
    B = collect(reim(b))
    C = collect(reim(c))
    AB = 0.5 * (A + B)
    BC = 0.5 * (B + C)
    M = [(A - B)'; (B - C)']
    rhs = [dot(AB, A - B), dot(BC, B - C)]

    Z = M \ rhs
    return Z[1] + im * Z[2]
end

"""
`non_colinear_check(a,b,c)`: test if the complex numbers are distinct
and noncollinear.
"""
function non_colinear_check(a::Complex, b::Complex, c::Complex)::Bool
    if a == b || b == c || a == c
        return false
    end
    z = (b - a) / (c - a)
    return imag(z) != 0
end

"""
    resize_gr_window(wide=800, tall=600)

Change the size of the `gr` drawing window. This is particularly
useful in VS code.
"""
function resize_gr_window(wide::Int = 800, tall::Int = 600)
    gr(size = (wide, tall))
end
export resize_gr_window

include("cubic.jl")
include("curve.jl")
include("my_spy.jl")
include("axes.jl")
end


text:
# Misc.jl

function checks_weights(v::AbstractVector{<:Real}, w::AbstractVector{<:Real})
    length(v) == length(w) ? nothing : throw(ArgumentError("`v` and `w` vectors must be the same size, got $(length(v)) and $(length(w))"))
    any([isnan(x) for x in w]) ? throw(ArgumentError("`w` vector cannot contain NaN values")) : nothing
    all(w .>= 0) ? nothing : throw(ArgumentError("`w` vector cannot contain negative entries"))
 end

 function checks_weights(v::AbstractVector{<:Real}, w::AbstractWeights)
    length(v) == length(w) ? nothing : throw(ArgumentError("`v` and `w` vectors must be the same size, got $(length(v)) and $(length(w))"))
    any([isnan(x) for x in w]) ? throw(ArgumentError("`w` vector cannot contain NaN values")) : nothing
    all(w .>= 0) ? nothing : throw(ArgumentError("`w` vector cannot contain negative entries"))
 end


text:
export local_field, inertia_degree, absolute_inertia_degree, absolute_ramification_index,
        eisenstein_extension, unramified_extension

################################################################################
#
#  Show function
#
################################################################################

function show(io::IO, a::LocalField{S, EisensteinLocalField}) where S
  print(io, "Eisenstein extension with defining polynomial ", defining_polynomial(a))
  print(io, " over ", base_field(a))
end

function show(io::IO, a::LocalField{S, UnramifiedLocalField}) where S
  print(io, "Unramified extension with defining polynomial ", defining_polynomial(a))
  print(io, " over ", base_field(a))
end

function show(io::IO, a::LocalField{S, GenericLocalField}) where S
  print(io, "Extension with defining polynomial ", defining_polynomial(a))
  print(io, " over ", base_field(a))
end

################################################################################
#
#  Characteristic
#
################################################################################

characteristic(K::LocalField) = characteristic(base_field(K))
prime(K::LocalField) = prime(base_field(K))

################################################################################
#
#  Type derivation
#
################################################################################

base_field_type(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter} = parent_type(S)
base_field_type(::Type{LocalField{S, T}}) where {S <: FieldElem, T <: LocalFieldParameter} = parent_type(S)

elem_type(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter} = LocalFieldElem{S, T}
elem_type(::Type{LocalField{S, T}}) where {S <: FieldElem, T <: LocalFieldParameter} = LocalFieldElem{S, T}

dense_matrix_type(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter} =  Generic.MatSpaceElem{LocalFieldElem{S, T}}
dense_matrix_type(::Type{LocalField{S, T}}) where {S <: FieldElem, T <: LocalFieldParameter} =  Generic.MatSpaceElem{LocalFieldElem{S, T}}

dense_poly_type(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter} = Generic.Poly{LocalFieldElem{S, T}}
dense_poly_type(::Type{LocalField{S, T}}) where {S <: FieldElem, T <: LocalFieldParameter} = Generic.Poly{LocalFieldElem{S, T}}

################################################################################
#
#  Basic predicates
#
################################################################################

is_domain_type(::Type{S}) where S <: LocalField = true
is_exact_type(::Type{S}) where S <: LocalField = false
isfinite(K::LocalField) = isfinite(base_field(K))
is_infinite(K::LocalField) = is_infinite(base_field(K))

################################################################################
#
#  Generating polynomials properties
#
################################################################################

function is_eisenstein_polynomial(f::PolyElem{S}) where S <: Union{padic, qadic, LocalFieldElem}
  if !iszero(valuation(leading_coefficient(f)))
    return false
  end
  if !isone(absolute_ramification_index(base_ring(f))*valuation(constant_coefficient(f)))
    return false
  end
  for i = 1:degree(f)-1
    c = coeff(f, i)
    if !iszero(c) && valuation(c) <= 0
      return false
    end
  end
  return true
end

function is_eisenstein_polynomial(f::T, p::S) where {T <: Union{fmpq_poly, fmpz_poly}, S<: Union{fmpz, Int}}
  @assert is_prime(p)
  if !iszero(valuation(leading_coefficient(f), p))
    return false
  end
  if !isone(valuation(constant_coefficient(f), p))
    return false
  end
  for i = 1:degree(f)-1
    c = coeff(f, i)
    if !iszero(c) && valuation(c, p) <= 0
      return false
    end
  end
  return true
end

function is_eisenstein_polynomial(f::PolyElem{<:NumFieldElem}, p::NumFieldOrdIdl)
  @assert is_prime(p)
  if !iszero(valuation(leading_coefficient(f), p))
    return false
  end
  if !isone(valuation(constant_coefficient(f), p))
    return false
  end
  for i = 1:degree(f)-1
    c = coeff(f, i)
    if !iszero(c) && valuation(c, p) <= 0
      return false
    end
  end
  return true
end

function _generates_unramified_extension(f::PolyElem{S}) where S <: Union{padic, qadic, LocalFieldElem}
  K = base_ring(f)
  F, mF = ResidueField(K)
  g = map_coefficients(mF, f)
  return is_irreducible(g)
end

var(a::LocalField) = a.S

function gen(K::LocalField)
  g = gen(parent(defining_polynomial(K)))
  el = K(g)
  return setprecision!(el, precision(K))
end

################################################################################
#
#  Subfields
#
################################################################################

function prime_field(L::Union{FlintQadicField, LocalField})
  L = base_ring(defining_polynomial(L))
  while typeof(L) != FlintPadicField
    L = base_ring(defining_polynomial(L))
  end
  return L
end

function base_field(L::LocalField)
  return base_ring(L.defining_polynomial)
end

function absolute_base_field(L::LocalField)
  return absolute_base_field(base_field(L))
end

absolute_base_field(L::FlintPadicField) = L
absolute_base_field(L::FlintQadicField) = base_field(L)

################################################################################
#
#  Degree
#
################################################################################

function degree(K::LocalField)
  return degree(defining_polynomial(K))
end

function absolute_degree(::FlintPadicField)
  return 1
end

function absolute_degree(K::FlintQadicField)
  return degree(K)
end
function absolute_degree(K::LocalField)
  return degree(K)*absolute_degree(base_field(K))
end

################################################################################
#
#  Traces
#
################################################################################

function assure_traces(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  if isdefined(K, :traces_basis)
    return nothing
  end
  res = S[base_field(K)(degree(K))]
  append!(res, polynomial_to_power_sums(defining_polynomial(K), degree(K)-1))
  K.traces_basis = res
  return nothing
end

################################################################################
#
#  Ramification index
#
################################################################################

function ramification_index(K::FlintPadicField)
  return 1
end

function ramification_index(K::FlintQadicField)
  return 1
end

function ramification_index(K::LocalField{S, EisensteinLocalField}) where S <: FieldElem
  return degree(K)
end

function ramification_index(K::LocalField{S, UnramifiedLocalField}) where S <: FieldElem
  return 1
end

function ramification_index(K::LocalField{S, GenericLocalField}) where S <: FieldElem
  error("Not yet implemented")
end

absolute_ramification_index(K::PadicField) = 1
absolute_ramification_index(K::QadicField) = 1

function absolute_ramification_index(K::LocalField{S, T}) where {S <: FieldElem, T <: LocalFieldParameter}
  return ramification_index(K)*absolute_ramification_index(base_field(K))
end

function ramification_index(L::LocalField, K::LocalField)
  if base_field(L) === K
    return ramification_index(L)
  else
    return ramification_index(L)*ramification_index(base_field(L), K)
  end
end

################################################################################
#
#  Inertia degree
#
################################################################################

function inertia_degree(K::FlintPadicField)
  return 1
end

function inertia_degree(K::FlintQadicField)
  return degree(K)
end

function inertia_degree(K::LocalField{S, EisensteinLocalField}) where S
  return 1
end

function inertia_degree(K::LocalField{S, UnramifiedLocalField}) where S
  return degree(K)
end

function inertia_degree(L::LocalField, K::LocalField)
  if base_field(L) === K
    return inertia_degree(L)
  else
    return inertia_degree(L)*inertia_degree(base_field(L), K)
  end
end

function absolute_inertia_degree(L::LocalField)
  return inertia_degree(L)*absolute_inertia_degree(base_field(L))
end

absolute_inertia_degree(::PadicField) = 1
absolute_inertia_degree(K::QadicField) = degree(K)

################################################################################
#
#  Basis
#
################################################################################

function basis(K::Union{FlintQadicField, LocalField})
  return powers(gen(K), degree(K)-1)
end

function absolute_basis(K::LocalField)
  Bk = absolute_basis(base_field(K))
  BKr = basis(K)
  BK = Vector{elem_type(K)}()
  for i = 1:length(BKr)
    for j = 1:length(Bk)
      push!(BK, BKr[i]*K(Bk[j]))
    end
  end
  return BK
end

absolute_basis(K::FlintQadicField) = basis(K)
absolute_basis(K::FlintPadicField) = padic[one(K)]

################################################################################
#
#  Constructors
#
################################################################################

#=
function find_irreducible_polynomial(K, n::Int)
  Zx, x = PolynomialRing(FlintZZ, "x", cached = false)
  f = cyclotomic(ppio(degree(K), n)*n, x)
  lf = factor(f, K)
  return first(keys(lf[1]))
end

function unramified_extension(L::LocalField, n::Int, prec::Int, s::String = "z")
  K, mK = ResidueField(L)
  f = find_irreducible_polynomial(K, n)
  coeffs =
  return local
end
=#

function eisenstein_extension(f::Generic.Poly{S}, s::String = "a"; check::Bool = true, cached::Bool = true) where S
  return local_field(f, s, EisensteinLocalField, check = check, cached = cached)
end

function unramified_extension(f::Generic.Poly{S}, s::String = "a"; check::Bool = true, cached::Bool = true) where S
  return local_field(f, s, UnramifiedLocalField, check = check, cached = cached)
end

function local_field(f::Generic.Poly{S},::Type{T}; check::Bool = true, cached::Bool = true) where {S <: FieldElem, T <: LocalFieldParameter}
  return local_field(f, "a", T, check = check, cached = cached)
end

function local_field(f::Generic.Poly{S}, s::String, ::Type{EisensteinLocalField}; check::Bool = true, cached::Bool = true) where {S <: FieldElem}
  symb = Symbol(s)
  if check && !is_eisenstein_polynomial(f)
    error("Defining polynomial is not Eisenstein")
  end
  K = LocalField{S, EisensteinLocalField}(f, symb)
  return K, gen(K)
end

function local_field(f::Generic.Poly{S}, s::String, ::Type{UnramifiedLocalField}; check::Bool = true, cached::Bool = true) where {S <: FieldElem}
  symb = Symbol(s)
  if check && !_generates_unramified_extension(f)
    error("Defining polynomial is not irreducible over the residue field!")
  end
  K = LocalField{S, UnramifiedLocalField}(f, symb)
  return K, gen(K)
end

function local_field(f::Generic.Poly{S}, s::String, ::Type{T} = GenericLocalField; check::Bool = true, cached::Bool = true) where {S <: FieldElem, T <: LocalFieldParameter}
  symb = Symbol(s)
  if check && !is_irreducible(f)
    error("Defining polynomial is not irreducible")
  end
  K = LocalField{S, T}(f, symb)
  return K, gen(K)
end

function local_field(f::fmpq_poly, p::Int, precision::Int, s::String, ::Type{T} = GenericLocalField; check::Bool = true, cached::Bool = true) where T <: LocalFieldParameter
  @assert is_prime(p)
  K = PadicField(p, precision)
  fK = map_coefficients(K, f)
  return local_field(fK, s, T, cached = cached, check = check)
end

function defining_polynomial(K::LocalField)
  return K.defining_polynomial
end

function precision(K::LocalField)
  return precision(defining_polynomial(K))*ramification_index(K)
end

function setprecision!(K::LocalField, n::Int)
  K.defining_polynomial = setprecision(defining_polynomial(K), n)
  return nothing
end

################################################################################
#
#  Uniformizer
#
################################################################################

function uniformizer(K::LocalField{S, EisensteinLocalField}) where S <: FieldElem
  return gen(K)
end

function uniformizer(K::LocalField{S, UnramifiedLocalField}) where S <: FieldElem
  return K(uniformizer(base_field(K)))
end

function uniformizer(K::LocalField{S, GenericLocalField}) where S <: FieldElem
  error("Not yet implemented")
end

################################################################################
#
#  ResidueField
#
################################################################################

function ResidueField(K::LocalField{S, EisensteinLocalField}) where {S <: FieldElem}
  if isdefined(K, :residue_field_map)
    mp = K.residue_field_map
    return codomain(mp), mp
  end
  k = base_field(K)
  ks, mks = ResidueField(k)

  function proj(a::LocalFieldElem)
    @assert parent(a) === K
    for i = 1:degree(a.data)
      if valuation(coeff(a, i)) < 0
        error("The projection is not well defined!")
      end
    end
    return mks(coeff(a, 0))
  end

  function lift(a)
    @assert parent(a) === ks
    return setprecision(K(mks\(a)), precision(K))
  end
  mp = MapFromFunc(proj, lift, K, ks)

  K.residue_field_map = mp

  return ks, mp
end

 ########### Residue field of unramified local field ext ################
function ResidueField(K::LocalField{ S, UnramifiedLocalField}) where {S <: FieldElem}
   if isdefined(K, :residue_field_map)
     mp = K.residue_field_map
     return codomain(mp), mp
   end
   k = base_field(K)
   ks, mks = ResidueField(k)
   Fpt = PolynomialRing(ks, cached = false)[1]
   g = defining_polynomial(K)
   f = Fpt([ks(mks(coeff(g, i))) for i=0:degree(K)])
   kk = FiniteField(f)[1]
   bas = basis(K)
   u = gen(kk)
   function proj(a:: Hecke.LocalFieldElem)
     col = typeof(kk(1))[]
     for i = 0:degree(K)-1
       push!(col, mks(coeff(a,i)) * u^i )
     end
     return sum(col)
   end
   function lift(b:: Hecke.RelFinFieldElem)
     col = typeof(K(1))[]
     for i = 0:degree(kk)-1
       push!(col, K(mks\(coeff(b,i))) * bas[i+1] )
     end
     return sum(col)
   end
   mp = MapFromFunc(proj, lift, K, kk)
   K.residue_field_map = mp
  return kk, mp
end

 ################### unramified extension over local field L of a given degree n ####################

 function unramified_extension(L::Union{FlintPadicField, FlintQadicField, LocalField}, n::Int)
   R, mR = ResidueField(L)
   f = polynomial(R, push!([rand(R) for i = 0:n-1], one(R)))
   while !is_irreducible(f)
     f = polynomial(R, push!([rand(R) for i = 0:n-1], one(R)))
   end
   f_L = polynomial(L, [mR\(coeff(f, i)) for i = 0:degree(f)])
   return unramified_extension(f_L)
 end


text:
# ------------------ #
# Exception handling #
# ------------------ #

struct VariableNotAllowedError <: Exception
    bad_var::Symbol
    eq::Expr
    shifts::Set{Int}
end

function Base.showerror(io::IO, v::VariableNotAllowedError)
    bad_exprs = [:($(v.bad_var)($shift)) for shift in v.shifts]
    bad_str = join(bad_exprs, ", ", " and ")
    print(io, "Invalid symbol(s) $(bad_str) found in equation $(v.eq). ")
    print(io, "Try adding $(v.bad_var) to the list of args")
end

# # ---------- #
# # Base Types #
# # ---------- #

const FlatArgs = AbstractVector
const GroupedArgs = AbstractDict{Symbol,<:Any}
const ArgType = Union{FlatArgs,GroupedArgs}

const FlatParams = AbstractVector
const GroupedParams = AbstractDict{Symbol,<:Any}
const ParamType = Union{FlatParams,GroupedParams}

struct Der{T} end
const TDer{n} = Type{Der{n}}

struct SkipArg end

_to_flat(f::FlatArgs) = f
_to_flat(g::GroupedArgs) = vcat(values(g)...)

# ---------------- #
# helper functions #
# ---------------- #
"""
    build_definition_map(defs::AbstractDict, incidence::IncidenceTable,
                         dynvars::Vector{Symbol})

For each key `k` in `defs`, and each time shift in `incidence.by_eq[k]`,
construct a mapping between the time shifted stringified symbol for `k` to the
stringified expression that it should be replaced with.

## Example

```jlcon
julia> defs = Dict(:x=>:(a(0)/(1-c(1))));

julia> incidence = Dolang.IncidenceTable(:(foo(0) = log(a(0))+b(0)/x(-1)));

julia> Dolang.build_definition_map(defs, incidence)
Dict{Symbol,Union{Expr, Number, Symbol}} with 1 entry:
  :_x_m1_ => :(_a_m1_ / (1 - _c__0_))
```
"""
function build_definition_map(defs::AbstractDict, incidence::IncidenceTable)
    out = Dict{Symbol,Union{Symbol,Expr,Number}}()

    # NOTE: we need to make the keys of definitions have the form
    #       (def, shift)::Tuple{Symbol,Int} in order for csubs to work
    #       properly below
    norm_defs = OrderedDict{Tuple{Symbol,Int},Any}()
    for (_def, _ex) in defs
        norm_defs[(_def, 0)] = _ex
    end

    for def_var in keys(defs)
        if haskey(incidence.by_var, def_var)
            _ex = csubs(defs[def_var], norm_defs)
            for time in incidence.by_var[def_var]
                new_key = stringify((def_var, time))
                out[new_key] = stringify(
                    time_shift(_ex, time)
                )
            end
        end
    end
    out
end

# --------------- #
# FunctionFactory #
# --------------- #

SymExpr = Union{Expr,Symbol,Number}

struct FunctionFactory
        # stringified equations
        equations::OrderedDict{Symbol,SymExpr}
        # list of group of (stringified) variables
        arguments::OrderedDict{Symbol, Vector{Symbol}}
        # list of assigned variables
        # preamble: definitions
        preamble::OrderedDict{Symbol, SymExpr}
        # name of function
        funname::Symbol
end

const FlatFunctionFactory = FunctionFactory

# function FunctionFactory(equations::OrderedDict, arguments::OrderedDict, definitions::AbstractDict; eliminate_preamble=false, funname=:anon)

#     # eqs: OrderedDict (targets are keys)
#     # args: OrderedDict
#     # defs:
#     all_vars = union([list_variables(eq) for eq in values(equations)]...)
#     present_vars = union(union(values(arguments)...), keys(equations))
#     filter!(x->(x isa Tuple), present_vars)

#     # variables appearing in equations
#     needed_defs = setdiff(all_vars, present_vars)

#     defs = solve_definitions(definitions, needed_defs)
#     defs_stringified = OrderedDict{Symbol,SymExpr}([stringify(k)=>stringify(v) for (k,v) in defs])
#     args_stringified = OrderedDict{Symbol,Vector{Symbol}}([k=>stringify.(v) for (k,v) in arguments])
#     eqs_stringified = OrderedDict{Symbol, SymExpr}([stringify(k)=>stringify(v) for (k,v) in equations])
#     targets = [keys(eqs_stringified)...]
#     return FlatFunctionFactory( eqs_stringified, args_stringified, targets, defs_stringified, funname)

# end

# function FunctionFactory(equations::Vector, arguments, definitions; kwargs...)
#     eqs = OrderedDict(Symbol("eq_",i)=>eq for (i,eq) in enumerate(equations))
#     FlatFunctionFactory(eqs, arguments, definitions; kwargs...)
# end

# function FunctionFactory(equations, arguments::Vector, definitions; kwargs...)
#     args = OrderedDict(Symbol("arg_",i)=>arg for (i,arg) in enumerate(arguments))
#     FunctionFactory(equations, args, definitions; kwargs...)
# end

# function FunctionFactory(equations::Vector, arguments::Vector, definitions; kwargs...)
#     eqs = OrderedDict(Symbol("eq_",i)=>eq for (i,eq) in enumerate(equations))
#     args = OrderedDict(Symbol("arg_",i)=>arg for (i,arg) in enumerate(arguments))
#     FunctionFactory(eqs, args, definitions; kwargs...)
# end

max_issues_repo_name: EconForge/Dolang