text:
using BinaryProvider
using Compat
using Compat.Test
using Compat.Libdl
using Compat.Pkg
using SHA

# The platform we're running on
const platform = platform_key()

# Useful command to launch `sh` on any platform
const sh = gen_sh_cmd

# Output of a few scripts we are going to run
const simple_out = "1\n2\n3\n4\n"
const long_out = join(["$(idx)\n" for idx in 1:100], "")
const newlines_out = join(["marco$(d)polo$(d)" for d in ("\n","\r","\r\n")], "")

# Explicitly probe platform engines in verbose mode to get coverage and make
# CI debugging easier
BinaryProvider.probe_platform_engines!(;verbose=true)

# Helper function to strip out color codes from strings to make it easier to
# compare output within tests that has been colorized
function strip_colorization(s)
    return replace(s, r"(\e\[\d+m)"m => "")
end

# Helper function to strip out log timestamps from strings
function strip_timestamps(s)
    return replace(s, r"^(\[\d\d:\d\d:\d\d\] )"m => "")
end

@testset "OutputCollector" begin
    cd("output_tests") do
        # Collect the output of `simple.sh``
        oc = OutputCollector(sh(`./simple.sh`))

        # Ensure we can wait on it and it exited properly
        @test wait(oc)

        # Ensure further waits are fast and still return 0
        let
            tstart = time()
            @test wait(oc)
            @test time() - tstart < 0.1
        end

        # Test that we can merge properly
        @test merge(oc) == simple_out

        # Test that merging twice works
        @test merge(oc) == simple_out

        # Test that `tail()` gives the same output as well
        @test tail(oc) == simple_out

        # Test that colorization works
        let
            red = Base.text_colors[:red]
            def = Base.text_colors[:default]
            gt = "1\n$(red)2\n$(def)3\n4\n"
            @test merge(oc; colored=true) == gt
            @test tail(oc; colored=true) == gt
        end

        # Test that we can grab stdout and stderr separately
        @test collect_stdout(oc) == "1\n3\n4\n"
        @test collect_stderr(oc) == "2\n"
    end

    # Next test a much longer output program
    cd("output_tests") do
        oc = OutputCollector(sh(`./long.sh`))

        # Test that it worked, we can read it, and tail() works
        @test wait(oc)
        @test merge(oc) == long_out
        @test tail(oc; len=10) == join(["$(idx)\n" for idx in 91:100], "")
    end

    # Next, test a command that fails
    cd("output_tests") do
        oc = OutputCollector(sh(`./fail.sh`))

        @test !wait(oc)
        @test merge(oc) == "1\n2\n"
    end

    # Next, test a command that kills itself (NOTE: This doesn't work on windows.  sigh.)
    @static if !Compat.Sys.iswindows()
        cd("output_tests") do
            oc = OutputCollector(sh(`./kill.sh`))

            @test !wait(oc)
            @test collect_stdout(oc) == "1\n2\n"
        end
    end

    # Next, test reading the output of a pipeline()
    grepline = pipeline(sh(`-c 'printf "Hello\nWorld\nJulia"'`), `grep ul`)
    oc = OutputCollector(grepline)

    @test wait(oc)
    @test merge(oc) == "Julia\n"

    # Next, test that \r and \r\n are treated like \n
    cd("output_tests") do
        oc = OutputCollector(sh(`./newlines.sh`))

        @test wait(oc)
        @test collect_stdout(oc) == newlines_out
    end

    # Next, test that tee'ing to a stream works
    cd("output_tests") do
        ios = IOBuffer()
        oc = OutputCollector(sh(`./simple.sh`); tee_stream=ios, verbose=true)
        @test wait(oc)
        @test merge(oc) == simple_out

        seekstart(ios)
        tee_out = String(read(ios))
        tee_out = strip_colorization(tee_out)
        tee_out = strip_timestamps(tee_out)
        @test tee_out == simple_out
    end

    # Also test that auto-tail'ing can be can be directed to a stream
    cd("output_tests") do
        ios = IOBuffer()
        oc = OutputCollector(sh(`./fail.sh`); tee_stream=ios)

        @test !wait(oc)
        @test merge(oc) == "1\n2\n"
        seekstart(ios)
        tee_out = String(read(ios))
        tee_out = strip_colorization(tee_out)
        tee_out = strip_timestamps(tee_out)
        @test tee_out == "1\n2\n"
    end

    # Also test that auto-tail'ing can be turned off
    cd("output_tests") do
        ios = IOBuffer()
        oc = OutputCollector(sh(`./fail.sh`); tee_stream=ios, tail_error=false)

        @test !wait(oc)
        @test merge(oc) == "1\n2\n"

        seekstart(ios)
        @test String(read(ios)) == ""
    end
end

@testset "PlatformNames" begin
    # Ensure the platform type constructors are well behaved
    @test_throws ArgumentError Linux(:not_a_platform)
    @test_throws ArgumentError Linux(:x86_64, :crazy_libc)
    @test_throws ArgumentError Linux(:x86_64, :glibc, :crazy_abi)
    @test_throws ArgumentError Linux(:x86_64, :glibc, :eabihf)
    @test_throws ArgumentError Linux(:armv7l, :glibc, :blank_abi)
    @test_throws ArgumentError MacOS(:i686)
    @test_throws ArgumentError MacOS(:x86_64, :glibc)
    @test_throws ArgumentError MacOS(:x86_64, :blank_libc, :eabihf)
    @test_throws ArgumentError Windows(:armv7l)
    @test_throws ArgumentError Windows(:x86_64, :glibc)
    @test_throws ArgumentError Windows(:x86_64, :blank_libc, :eabihf)
    @test_throws ArgumentError FreeBSD(:not_a_platform)
    @test_throws ArgumentError FreeBSD(:x86_64, :crazy_libc)
    @test_throws ArgumentError FreeBSD(:x86_64, :blank_libc, :crazy_abi)
    @test_throws ArgumentError FreeBSD(:x86_64, :blank_libc, :eabihf)
    @test_throws ArgumentError FreeBSD(:armv7l, :blank_libc, :blank_abi)

    # Test that we can get that arch of various platforms
    @test arch(Linux(:aarch64, :musl)) == :aarch64
    @test arch(Windows(:i686)) == :i686
    @test arch(UnknownPlatform()) == :unknown
    @test arch(FreeBSD(:amd64)) == :x86_64
    @test arch(FreeBSD(:i386)) == :i686

    # Test that our platform_dlext stuff works
    @test platform_dlext(Linux(:x86_64)) == platform_dlext(Linux(:i686))
    @test platform_dlext(Windows(:x86_64)) == platform_dlext(Windows(:i686))
    @test platform_dlext(MacOS()) != platform_dlext(Linux(:armv7l))
    @test platform_dlext(FreeBSD(:x86_64)) == platform_dlext(Linux(:x86_64))
    @test platform_dlext(UnknownPlatform()) == "unknown"
    @test platform_dlext() == platform_dlext(platform_key())

    # Test some valid dynamic library paths
    @test valid_dl_path("libfoo.so.1.2.3", Linux(:x86_64))
    @test valid_dl_path("libfoo-1.dll", Windows(:x86_64))
    @test valid_dl_path("libfoo.1.2.3.dylib", MacOS())
    @test !valid_dl_path("libfoo.dylib", Linux(:x86_64))
    @test !valid_dl_path("libfoo.so", Windows(:x86_64))

    # Make sure the platform_key() with explicit triplet works
    @test platform_key("x86_64-linux-gnu") == Linux(:x86_64)
    @test platform_key("i686-unknown-linux-gnu") == Linux(:i686)
    @test platform_key("x86_64-apple-darwin14") == MacOS()
    @test platform_key("armv7l-pc-linux-gnueabihf") == Linux(:armv7l)
    @test platform_key("arm-linux-gnueabihf") == Linux(:armv7l)
    @test platform_key("aarch64-unknown-linux-gnu") == Linux(:aarch64)
    @test platform_key("powerpc64le-linux-gnu") == Linux(:powerpc64le)
    @test platform_key("ppc64le-linux-gnu") == Linux(:powerpc64le)
    @test platform_key("x86_64-w64-mingw32") == Windows(:x86_64)
    @test platform_key("i686-w64-mingw32") == Windows(:i686)
    @test platform_key("x86_64-unknown-freebsd11.1") == FreeBSD(:x86_64)
    @test platform_key("i686-unknown-freebsd11.1") == FreeBSD(:i686)
    @test platform_key("amd64-unknown-freebsd12.0") == FreeBSD(:x86_64)
    @test platform_key("i386-unknown-freebsd10.3") == FreeBSD(:i686)

    # Make sure some of these things are rejected
    @test platform_key("totally FREEFORM text!!1!!!1!") == UnknownPlatform()
    @test platform_key("invalid-triplet-here") == UnknownPlatform()
    @test platform_key("aarch64-linux-gnueabihf") == UnknownPlatform()
    @test platform_key("armv7l-linux-gnu") == UnknownPlatform()
    @test platform_key("x86_64-w32-mingw64") == UnknownPlatform()

    # Test that we can indeed ask if something is linux or windows, etc...
    @test Compat.Sys.islinux(Linux(:aarch64))
    @test !Compat.Sys.islinux(Windows(:x86_64))
    @test Compat.Sys.iswindows(Windows(:i686))
    @test !Compat.Sys.iswindows(Linux(:x86_64))
    @test Compat.Sys.isapple(MacOS())
    @test !Compat.Sys.isapple(Linux(:powerpc64le))
    @test Compat.Sys.isbsd(MacOS())
    @test Compat.Sys.isbsd(FreeBSD(:x86_64))
    @test !Compat.Sys.isbsd(Linux(:powerpc64le, :musl))

    @test wordsize(Linux(:i686)) == wordsize(Linux(:armv7l)) == 32
    @test wordsize(MacOS()) == wordsize(Linux(:aarch64)) == 64
    @test wordsize(FreeBSD(:x86_64)) == wordsize(Linux(:powerpc64le)) == 64
    @test wordsize(UnknownPlatform()) == 0

    @test triplet(Windows(:i686)) == "i686-w64-mingw32"
    @test triplet(Linux(:x86_64, :musl)) == "x86_64-linux-musl"
    @test triplet(Linux(:armv7l, :musl)) == "arm-linux-musleabihf"
    @test triplet(Linux(:x86_64)) == "x86_64-linux-gnu"
    @test triplet(Linux(:armv7l)) == "arm-linux-gnueabihf"
    @test triplet(MacOS()) == "x86_64-apple-darwin14"
    @test triplet(FreeBSD(:x86_64)) == "x86_64-unknown-freebsd11.1"
    @test triplet(FreeBSD(:i686)) == "i686-unknown-freebsd11.1"
    @test triplet(UnknownPlatform()) == "unknown-unknown-unknown"

    @test repr(Windows(:x86_64)) == "Windows(:x86_64)"
    @test repr(Linux(:x86_64, :glibc, :blank_abi)) == "Linux(:x86_64, :glibc)"
    @test repr(MacOS()) == "MacOS(:x86_64)"

    for p in [Windows(:i686), Linux(:armv7l, :musl), FreeBSD(:x86_64), MacOS()]
        fakepath = "/path/to/nowhere/thingo." * triplet(p) * ".tar.gz"
        @test extract_platform_key(fakepath) == p
    end
end

@testset "Prefix" begin
    mktempdir() do temp_dir
        prefix = Prefix(temp_dir)

        # Test that it's taking the absolute path
        @test prefix.path == abspath(temp_dir)

        # Test that `bindir()`, `libdir()` and `includedir()` all work
        for dir in unique([bindir(prefix), libdir(prefix), includedir(prefix)])
            @test !isdir(dir)
            mkpath(dir)
        end

        # Create a little script within the bindir to ensure we can run it
        ppt_path = joinpath(bindir(prefix), "prefix_path_test.sh")
        open(ppt_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        chmod(ppt_path, 0o775)

        # Test that activation adds certain paths to our environment variables
        activate(prefix)

        # PATH[1] should be "<prefix>/bin" now
        @test BinaryProvider.split_PATH()[1] == bindir(prefix)
        @test Libdl.DL_LOAD_PATH[1] == libdir(prefix)

        # Test we can run the script we dropped within this prefix.  Once again,
        # something about Windows | busybox | Julia won't pick this up even though
        # the path clearly points to the file.  :(
        @static if !Compat.Sys.iswindows()
            @test success(sh(`$(ppt_path)`))
            @test success(sh(`prefix_path_test.sh`))
        end

        # Now deactivate and make sure that all traces are gone
        deactivate(prefix)
        @test BinaryProvider.split_PATH()[1] != bindir(prefix)
        @test Libdl.DL_LOAD_PATH[1] != libdir(prefix)

        # Test that we can control libdir() via platform arguments
        @test libdir(prefix, Linux(:x86_64)) == joinpath(prefix, "lib")
        @test libdir(prefix, Windows(:x86_64)) == joinpath(prefix, "bin")
    end
end

@testset "Products" begin
    temp_prefix() do prefix
        # Test that basic satisfication is not guaranteed
        e_path = joinpath(bindir(prefix), "fooifier")
        l_path = joinpath(libdir(prefix), "libfoo.$(Libdl.dlext)")
        e = ExecutableProduct(prefix, "fooifier", :fooifier)
        ef = FileProduct(prefix, joinpath("bin", "fooifier"), :fooifier)
        l = LibraryProduct(prefix, "libfoo", :libfoo)
        lf = FileProduct(l_path, :libfoo)

        @test !satisfied(e; verbose=true)
        @test !satisfied(ef; verbose=true)
        @test !satisfied(l, verbose=true)
        @test !satisfied(l, verbose=true, isolate=true)
        @test !satisfied(lf, verbose=true)

        # Test that simply creating a file that is not executable doesn't
        # satisfy an Executable Product (and say it's on Linux so it doesn't
        # complain about the lack of an .exe extension)
        mkpath(bindir(prefix))
        touch(e_path)
        @test satisfied(ef, verbose=true)
        @static if !Compat.Sys.iswindows()
            # Windows doesn't care about executable bit, grumble grumble
            @test !satisfied(e, verbose=true, platform=Linux(:x86_64))
        end

        # Make it executable and ensure this does satisfy the Executable
        chmod(e_path, 0o777)
        @test satisfied(e, verbose=true, platform=Linux(:x86_64))

        # Remove it and add a `$(path).exe` version to check again, this
        # time saying it's a Windows executable
        Base.rm(e_path; force=true)
        touch("$(e_path).exe")
        chmod("$(e_path).exe", 0o777)
        @test locate(e, platform=Windows(:x86_64)) == "$(e_path).exe"

        # Test that simply creating a library file doesn't satisfy it if we are
        # testing something that matches the current platform's dynamic library
        # naming scheme, because it must be `dlopen()`able.
        mkpath(libdir(prefix))
        touch(l_path)
        @test satisfied(lf, verbose=true)
        @test !satisfied(l, verbose=true)
        @test satisfied(lf, verbose=true, isolate=true)
        @test !satisfied(l, verbose=true, isolate=true)

        # But if it is from a different platform, simple existence will be
        # enough to satisfy a LibraryProduct
        @static if Compat.Sys.iswindows()
            p = Linux(:x86_64)
            mkpath(libdir(prefix, p))
            l_path = joinpath(libdir(prefix, p), "libfoo.so")
            touch(l_path)
            @test satisfied(l, verbose=true, platform=p)
            @test satisfied(l, verbose=true, platform=p, isolate=true)

            # Check LibraryProduct objects with explicit directory paths
            ld = LibraryProduct(libdir(prefix, p), "libfoo", :libfoo)
            @test satisfied(ld, verbose=true, platform=p)
            @test satisfied(ld, verbose=true, platform=p, isolate=true)
        else
            p = Windows(:x86_64)
            mkpath(libdir(prefix, p))
            l_path = joinpath(libdir(prefix, p), "libfoo.dll")
            touch(l_path)
            @test satisfied(l, verbose=true, platform=p)
            @test satisfied(l, verbose=true, platform=p, isolate=true)

            # Check LibraryProduct objects with explicit directory paths
            ld = LibraryProduct(libdir(prefix, p), "libfoo", :libfoo)
            @test satisfied(ld, verbose=true, platform=p)
            @test satisfied(ld, verbose=true, platform=p, isolate=true)
        end
    end

    # Ensure that the test suite thinks that these libraries are foreign
    # so that it doesn't try to `dlopen()` them:
    foreign_platform = @static if platform_key() == Linux(:aarch64)
        # Arbitrary architecture that is not dlopen()'able
        Linux(:powerpc64le)
    else
        # If we're not Linux(:aarch64), then say the libraries are
        Linux(:aarch64)
    end

    # Test for valid library name permutations
    for ext in ["1.so", "so", "so.1", "so.1.2", "so.1.2.3"]
        temp_prefix() do prefix
            l_path = joinpath(libdir(prefix, foreign_platform), "libfoo.$ext")
            l = LibraryProduct(prefix, "libfoo", :libfoo)
            mkdir(dirname(l_path))
            touch(l_path)
            @test satisfied(l; verbose=true, platform=foreign_platform)
        end
    end

    # Test for invalid library name permutations
    for ext in ["so.1.2.3a", "so.1.a"]
        temp_prefix() do prefix
            l_path = joinpath(libdir(prefix, foreign_platform), "libfoo.$ext")
            l = LibraryProduct(prefix, "libfoo", :libfoo)
            mkdir(dirname(l_path))
            touch(l_path)
            @test !satisfied(l; verbose=true, platform=foreign_platform)
        end
    end

    # Test for proper repr behavior
    temp_prefix() do prefix
        l = LibraryProduct(prefix, "libfoo", :libfoo)
        @test repr(l) == "LibraryProduct(prefix, $(repr(["libfoo"])), :libfoo)"
        l = LibraryProduct(libdir(prefix), ["libfoo", "libfoo2"], :libfoo)
        @test repr(l) == "LibraryProduct($(repr(libdir(prefix))), $(repr(["libfoo", "libfoo2"])), :libfoo)"

        e = ExecutableProduct(prefix, "fooifier", :fooifier)
        @test repr(e) == "ExecutableProduct(prefix, \"fooifier\", :fooifier)"
        e = ExecutableProduct(joinpath(bindir(prefix), "fooifier"), :fooifier)
        @test repr(e) == "ExecutableProduct($(repr(joinpath(bindir(prefix), "fooifier"))), :fooifier)"

        f = FileProduct(prefix, joinpath("etc", "fooifier"), :foo_conf)
        @test repr(f) == "FileProduct(prefix, $(repr(joinpath("etc", "fooifier"))), :foo_conf)"

        f = FileProduct(joinpath(prefix, "etc", "foo.conf"), :foo_conf)
        @test repr(f) == "FileProduct($(repr(joinpath(prefix, "etc", "foo.conf"))), :foo_conf)"
    end
end

@testset "Packaging" begin
    # Clear out previous build products
    for f in readdir(".")
        if !endswith(f, ".tar.gz") && !endswith(f, ".sha256")
            continue
        end
        Base.rm(f; force=true)
    end

    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    temp_prefix() do prefix
        # Create random files
        mkpath(bindir(prefix))
        mkpath(libdir(prefix))
        mkpath(joinpath(prefix, "etc"))
        bar_path = joinpath(bindir(prefix), "bar.sh")
        open(bar_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        baz_path = joinpath(libdir(prefix), "baz.so")
        open(baz_path, "w") do f
            write(f, "this is not an actual .so\n")
        end

        qux_path = joinpath(prefix, "etc", "qux.conf")
        open(qux_path, "w") do f
            write(f, "use_julia=true\n")
        end

        # Next, package it up as a .tar.gz file
        tarball_path, tarball_hash = package(prefix, "./libfoo", v"1.0.0"; verbose=true)
        @test isfile(tarball_path)

        # Check that we are calculating the hash properly
        tarball_hash_check = open(tarball_path, "r") do f
            bytes2hex(sha256(f))
        end
        @test tarball_hash_check == tarball_hash

        # Test that packaging into a file that already exists fails
        @test_throws ErrorException package(prefix, "./libfoo", v"1.0.0")
    end

    # Test that we can inspect the contents of the tarball
    contents = list_tarball_files(tarball_path)
    libdir_name = Compat.Sys.iswindows() ? "bin" : "lib"
    @test joinpath("bin", "bar.sh") in contents
    @test joinpath(libdir_name, "baz.so") in contents
    @test joinpath("etc", "qux.conf") in contents

    # Install it within a new Prefix
    temp_prefix() do prefix
        # Install the thing
        @test install(tarball_path, tarball_hash; prefix=prefix, verbose=true)

        # Ensure we can use it
        bar_path = joinpath(bindir(prefix), "bar.sh")
        baz_path = joinpath(libdir(prefix), "baz.so")
        qux_path = joinpath(prefix, "etc", "qux.conf")

        @test isfile(bar_path)
        @test isfile(baz_path)
        @test isfile(qux_path)

        # Ask for the manifest that contains these files to ensure it works
        manifest_path = manifest_for_file(bar_path; prefix=prefix)
        @test isfile(manifest_path)
        manifest_path = manifest_for_file(baz_path; prefix=prefix)
        @test isfile(manifest_path)
        manifest_path = manifest_for_file(qux_path; prefix=prefix)
        @test isfile(manifest_path)

        # Ensure that manifest_for_file doesn't work on nonexistent files
        @test_throws ErrorException manifest_for_file("nonexistent"; prefix=prefix)

        # Ensure that manifest_for_file doesn't work on orphan files
        orphan_path = joinpath(bindir(prefix), "orphan_file")
        touch(orphan_path)
        @test isfile(orphan_path)
        @test_throws ErrorException manifest_for_file(orphan_path; prefix=prefix)

        # Ensure that trying to install again over our existing files is an error
        @test_throws ErrorException install(tarball_path, tarball_hash; prefix=prefix)

        # Ensure we can uninstall this tarball
        @test isinstalled(tarball_path, tarball_hash; prefix=prefix)
        Base.rm(bar_path)
        @test !isinstalled(tarball_path, tarball_hash; prefix=prefix)
        @test uninstall(manifest_path; verbose=true)
        @test !isinstalled(tarball_path, tarball_hash; prefix=prefix)
        @test !isfile(bar_path)
        @test !isfile(baz_path)
        @test !isfile(qux_path)
        @test !isfile(manifest_path)

        # Also make sure that it removes the entire etc directory since it's now empty
        @test !isdir(dirname(qux_path))

        # Ensure that we don't want to install tarballs from other platforms
        other_path = "./libfoo.x86_64-juliaos-chartreuse.tar.gz"
        cp(tarball_path, other_path)
        @test_throws ArgumentError install(other_path, tarball_hash; prefix=prefix)
        Base.rm(other_path; force=true)

        # Ensure that hash mismatches throw errors
        fake_hash = reverse(tarball_hash)
        @test_throws ErrorException install(tarball_path, fake_hash; prefix=prefix)
    end

    # Get a valid platform tarball name that is not our native platform
    other_platform = Linux(:x86_64)
    if platform == Linux(:x86_64)
        other_platform = MacOS()
    end
    new_tarball_path = "libfoo.$(triplet(other_platform)).tar.gz"
    cp(tarball_path, new_tarball_path)

    # Also generate a totally bogus tarball pathname
    bogus_tarball_path = "libfoo.not-a-triplet.tar.gz"
    cp(tarball_path, bogus_tarball_path)

    # Check that installation fails with a valid but "incorrect" platform, but can be forced
    temp_prefix() do prefix
        @test_throws ArgumentError install(new_tarball_path, tarball_hash; prefix=prefix, verbose=true)
        @test install(new_tarball_path, tarball_hash; prefix=prefix, verbose=true, ignore_platform = true)
    end

    # Next, check installing with a bogus platform also fails, but can be forced
    temp_prefix() do prefix
        @test_throws ArgumentError install(bogus_tarball_path, tarball_hash; prefix=prefix, verbose=true)
        @test install(bogus_tarball_path, tarball_hash; prefix=prefix, verbose=true, ignore_platform = true)
    end

    # Cleanup after ourselves
    Base.rm(tarball_path; force=true)
    Base.rm("$(tarball_path).sha256"; force=true)
    Base.rm(bogus_tarball_path; force=true)
    Base.rm("$(bogus_tarball_path).sha256"; force=true)
    Base.rm(new_tarball_path; force=true)
    Base.rm("$(new_tarball_path).sha256"; force=true)
end

@testset "Verification" begin
    temp_prefix() do prefix
        foo_path = joinpath(prefix, "foo")
        open(foo_path, "w") do file
            write(file, "test")
        end
        foo_hash = bytes2hex(sha256("test"))

        # Check that verifying with the right hash works
        Compat.@info("This should say; no hash cache found")
        ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
        @test ret == true
        @test status == :hash_cache_missing

        # Check that it created a .sha256 file
        @test isfile("$(foo_path).sha256")

        # Check that it verifies the second time around properly
        Compat.@info("This should say; hash cache is consistent")
        ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
        @test ret == true
        @test status == :hash_cache_consistent

        # Sleep for imprecise filesystems
        sleep(2)

        # Get coverage of messing with different parts of the verification chain
        touch(foo_path)
        Compat.@info("This should say; file has been modified")
        ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
        @test ret == true
        @test status == :file_modified

        # Ensure that we throw an exception when we can't verify properly
        @test_throws ErrorException verify(foo_path, "0"^32; verbose=true)

        # Ensure that messing with the hash file works properly
        touch(foo_path)
        @test verify(foo_path, foo_hash; verbose=true)
        open("$(foo_path).sha256", "w") do file
            write(file, "this is not the right hash")
        end
        Compat.@info("This should say; hash has changed")
        ret, status = verify(foo_path, foo_hash; verbose=true, report_cache_status=true)
        @test ret == true
        @test status == :hash_cache_mismatch
    end
end

const socrates_urls = [
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.gz" =>
    "e65d2f13f2085f2c279830e863292312a72930fee5ba3c792b14c33ce5c5cc58",
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.bz2" =>
    "13fc17b97be41763b02cbb80e9d048302cec3bd3d446c2ed6e8210bddcd3ac76",
    "https://github.com/staticfloat/small_bin/raw/f1a92f5eafbd30a0c6a8efb6947485b0f6d1bec3/socrates.tar.xz" =>
    "61bcf109fcb749ee7b6a570a6057602c08c836b6f81091eab7aa5f5870ec6475",
]
const socrates_hash = "adcbcf15674eafe8905093183d9ab997cbfba9056fc7dde8bfa5a22dfcfb4967"

@testset "Downloading" begin
    for (url, hash) in socrates_urls
        temp_prefix() do prefix
            Base.rm(prefix.path; recursive=true, force=true)
            download_verify_unpack(url, hash, prefix.path; verbose=true)
            socrates_path = joinpath(prefix, "bin", "socrates")
            @test isfile(socrates_path)

            unpacked_hash = open(socrates_path) do f
                bytes2hex(sha256(f))
            end
            @test unpacked_hash == socrates_hash
        end
    end
end

# Use `build_libfoo_tarball.jl` in the BinaryBuilder.jl repository to generate more of these
const bin_prefix = "https://github.com/staticfloat/small_bin/raw/51b13b44feb2a262e2e04690bfa54d03167533f2/libfoo"
const libfoo_downloads = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/libfoo.aarch64-linux-gnu.tar.gz", "36886ac25cf5678c01fe20630b413f9354b7a3721c6a2c2043162f7ebd147ff5"),
    Linux(:armv7l, :glibc)  => ("$bin_prefix/libfoo.arm-linux-gnueabihf.tar.gz", "147ebaeb1a722da43ee08705689aed71ac87c3c2c907af047c6721c0025ba383"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/libfoo.powerpc64le-linux-gnu.tar.gz", "5c35295ac161272ada9a77d1f6b770e30ea864e521e31853258cbc36ad4c4468"),
    Linux(:i686, :glibc)    => ("$bin_prefix/libfoo.i686-linux-gnu.tar.gz", "97655b6a218d61284723b6923d7c96e6a256fa68b9419d723c588aa24404b102"),
    Linux(:x86_64, :glibc)  => ("$bin_prefix/libfoo.x86_64-linux-gnu.tar.gz", "5208c63a9d07e592c78f541fc13caa8cd191b11e7e77b31d407237c2b13ec391"),

    Linux(:aarch64, :musl)  => ("$bin_prefix/libfoo.aarch64-linux-musl.tar.gz", "81751477c1e3ee6c93e1c28ee7db2b99d1eed0d6ce86dc30d64c2e5dd4dfe88d"),
    Linux(:armv7l, :musl)   => ("$bin_prefix/libfoo.arm-linux-musleabihf.tar.gz", "bb65aad58f2e6fc39dc9688da1bca5e8103a3a3fa67dc589debbd2e98176f0e1"),
    Linux(:i686, :musl)     => ("$bin_prefix/libfoo.i686-linux-musl.tar.gz", "5f02fd1fe19f3a565fb128d3673b35c7b3214a101cef9dcbb202c0092438a87b"),
    Linux(:x86_64, :musl)   => ("$bin_prefix/libfoo.x86_64-linux-musl.tar.gz", "ea630600a12d2c1846bc93bcc8d9638a4991f63329205c534d93e0a3de5f641d"),

    FreeBSD(:x86_64)        => ("$bin_prefix/libfoo.x86_64-unknown-freebsd11.1.tar.gz", "5f6edd6247b3685fa5c42c98a53d2a3e1eef6242c2bb3cdbb5fe23f538703fe4"),
    MacOS(:x86_64)          => ("$bin_prefix/libfoo.x86_64-apple-darwin14.tar.gz", "fcc268772d6f21d65b45fcf3854a3142679b78e53c7673dac26c95d6ccc89a24"),

    Windows(:i686)          => ("$bin_prefix/libfoo.i686-w64-mingw32.tar.gz", "79181cf62ca8e0b2e0851fa0ace52f4ab335d0cad26fb7f9cd4ff356a9a96e70"),
    Windows(:x86_64)        => ("$bin_prefix/libfoo.x86_64-w64-mingw32.tar.gz", "7f8939e9529835b83810d3ae7e2556f6e002d571f619894e54ece42ea5262b7f"),
)

# Test manually downloading and using libfoo
@testset "Libfoo Downloading" begin
    temp_prefix() do prefix
        foo_path = joinpath(prefix,"foo")
        touch(foo_path)
        # Quick one-off tests for `safe_isfile()`:
        @test BinaryProvider.safe_isfile(foo_path)
        @test !BinaryProvider.safe_isfile("http://google.com")

        if !haskey(libfoo_downloads, platform)
            Compat.@warn("Platform $platform does not have a libfoo download, skipping download tests")
        else
            # Test a good download works
            url, hash = libfoo_downloads[platform]
            @test install(url, hash; prefix=prefix, verbose=true)

            fooifier = ExecutableProduct(prefix, "fooifier", :fooifier)
            libfoo = LibraryProduct(prefix, "libfoo", :libfoo)

            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)

            fooifier_path = locate(fooifier)
            libfoo_path = locate(libfoo)

            # We know that foo(a, b) returns 2*a^2 - b
            result = 2*2.2^2 - 1.1

            # Test that we can invoke fooifier
            @test !success(`$fooifier_path`)
            @test success(`$fooifier_path 1.5 2.0`)
            @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

            # Test that we can dlopen() libfoo and invoke it directly
            hdl = Libdl.dlopen_e(libfoo_path)
            @test hdl != C_NULL
            foo = Libdl.dlsym_e(hdl, :foo)
            @test foo != C_NULL
            @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
            Libdl.dlclose(hdl)

            # Test uninstallation
            @test uninstall(manifest_from_url(url; prefix=prefix); verbose=true)

            # Test that download_verify_unpack() works
            Base.rm(prefix.path; recursive=true, force=true)
            download_verify_unpack(url, hash, prefix.path)
            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)

            # Test that download_verify twice in a row works, and that mucking
            # with the file causes a redownload if `force` is true:
            tmpfile = joinpath(prefix, "libfoo.tar.gz")
            @test download_verify(url, hash, tmpfile; verbose=true)
            @test download_verify(url, hash, tmpfile; verbose=true)

            # We sleep for at least a second here so that filesystems with low
            # precision in their mtime implementations don't get confused
            sleep(2)

            open(tmpfile, "w") do f
                write(f, "not the correct contents")
            end

            @test_throws ErrorException download_verify(url, hash, tmpfile; verbose=true)

            # This should return `false`, signifying that the download had to erase
            # the previously downloaded file.
            @test !download_verify(url, hash, tmpfile; verbose=true, force=true)

            # Now let's test that install() works the same way; freaking out if
            # the local path has been messed with, unless `force` has been given:
            tarball_path = joinpath(prefix, "downloads", basename(url))
            try mkpath(dirname(tarball_path)) catch; end
            open(tarball_path, "w") do f
                write(f, "not the correct contents")
            end

            @test_throws ErrorException install(url, hash; prefix=prefix, verbose=true)
            @test install(url, hash; prefix=prefix, verbose=true, force=true)
            @test isinstalled(url, hash; prefix=prefix)
            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)

            # Test that installing with a custom tarball_path works:
            tarball_path = joinpath(prefix, "downloads2", "tarball.tar.gz")
            @test install(url, hash; prefix=prefix, tarball_path=tarball_path, verbose=true, force=true)

            # Check that the tarball exists and hashes properly
            @test isfile(tarball_path)
            hash_check = open(tarball_path, "r") do f
                bytes2hex(sha256(f))
            end
            @test hash_check == hash

            # Check that we're still satisfied
            @test isinstalled(url, hash; prefix=prefix)
            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)
        end

        # Test a bad download fails properly
        bad_url = "http://localhost:1/this_is_not_a_file.$(triplet(platform)).tar.gz"
        bad_hash = "0"^64
        @test_throws ErrorException install(bad_url, bad_hash; prefix=prefix, verbose=true)
    end
end

# Test installation and failure modes of the bundled LibFoo.jl
@testset "LibFoo.jl" begin
    color="--color=$(Base.have_color ? "yes" : "no")"
    cd("LibFoo.jl") do
        Base.rm("./deps/deps.jl"; force=true)
        Base.rm("./deps/usr"; force=true, recursive=true)

        # Install `libfoo` and build the `deps.jl` file for `LibFoo.jl`
        coverage = "--code-coverage=$(Base.JLOptions().code_coverage != 0 ? "user" : "none")"
        run(`$(Base.julia_cmd()) $(coverage) $(color) deps/build.jl`)

        # Ensure `deps.jl` was actually created
        @test isfile("deps/deps.jl")
    end

    # Test that the generated deps.jl file gives us the important stuff
    cd("LibFoo.jl/deps") do
        dllist = Libdl.dllist()
        libjulia = filter(x -> occursin("libjulia", x), dllist)[1]
        julia_libdir = joinpath(dirname(libjulia), "julia")
        envvar_name = @static if Compat.Sys.isapple()
            "DYLD_LIBRARY_PATH"
        else Compat.Sys.islinux()
            "LD_LIBRARY_PATH"
        end

        original_libdirs = split(get(ENV, envvar_name, ""), ":")
        @static if !Compat.Sys.iswindows()
            original_libdirs = filter(x -> x != julia_libdir, original_libdirs)
            ENV[envvar_name] = join(original_libdirs, ":")
        end

        # Include deps.jl, run check_deps() and see if we get our products,
        # and also if the julia libdir gets added to the end of *_LIBRARY_PATH
        include("LibFoo.jl/deps/deps.jl")
        Base.invokelatest(check_deps)
        @test isfile(libfoo)
        @test isfile(fooifier)

        @static if !Compat.Sys.iswindows()
            @test julia_libdir in split(ENV[envvar_name], ":")
        end
    end

    cd("LibFoo.jl") do
        if VERSION < v"0.7-"
            # Now, run `LibFoo.jl`'s tests, adding `LibFoo.jl` to the LOAD_PATH
            # so that the tests can pick up the `LibFoo` module
            withenv("JULIA_LOAD_PATH"=>joinpath(pwd(),"src")) do
                run(`$(Base.julia_cmd()) $(color) test/runtests.jl`)
            end
        else
            # On julia 0.7, we can now rely on Project.toml to set the load path
            run(`$(Base.julia_cmd()) --project=$(pwd()) $(color) -e 'using Pkg; Pkg.test("LibFoo")'`)
        end
    end
end


text:
import Base.isempty

export AbstractPolytope,
       vertices_list,
       singleton_list,
       isempty

"""
    AbstractPolytope{N<:Real} <: AbstractPolyhedron{N}

Abstract type for polytopic sets, i.e., bounded sets with finitely many flat
facets, or equivalently, bounded sets defined as an intersection of a finite
number of half-spaces, or equivalently, bounded sets with finitely many
vertices.

### Notes

Every concrete `AbstractPolytope` must define the following functions:
- `vertices_list(::AbstractPolytope{N})::Vector{Vector{N}}` -- return a list of
    all vertices

```jldoctest
julia> subtypes(AbstractPolytope)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope
```
"""
abstract type AbstractPolytope{N<:Real} <: AbstractPolyhedron{N} end

# =============================================
# Common AbstractPolytope functions
# =============================================

"""
    isbounded(P::AbstractPolytope)::Bool

Determine whether a polytopic set is bounded.

### Input

- `P` -- polytopic set

### Output

`true` (since a polytope must be bounded).
"""
function isbounded(::AbstractPolytope)::Bool
    return true
end

"""
    singleton_list(P::AbstractPolytope{N})::Vector{Singleton{N}} where {N<:Real}

Return the vertices of a polytopic set as a list of singletons.

### Input

- `P` -- polytopic set

### Output

List containing a singleton for each vertex.
"""
function singleton_list(P::AbstractPolytope{N}
                       )::Vector{Singleton{N}} where {N<:Real}
    return [Singleton(vi) for vi in vertices_list(P)]
end

"""
    isempty(P::AbstractPolytope)::Bool

Determine whether a polytope is empty.

### Input

- `P` -- abstract polytope

### Output

`true` if the given polytope contains no vertices, and `false` otherwise.

### Algorithm

This algorithm checks whether the `vertices_list` of the given polytope is empty
or not.
"""
function isempty(P::AbstractPolytope)::Bool
    return isempty(vertices_list(P))
end

function default_polyhedra_backend(P, N)
    @assert isdefined(@__MODULE__, :Polyhedra) "this function needs the package 'Polyhedra' to be loaded"
    error("no default backend for numeric type $N")
end

# given a polytope P, apply the linear map P to each vertex of P
# it is assumed that the interface function `vertices_list(P)` is available
@inline function _linear_map_vrep(M::AbstractMatrix{N}, P::AbstractPolytope{N}) where {N<:Real}
    vertices = broadcast(v -> M * v, vertices_list(P))
    m = size(M, 1) # output dimension
    if m == 1
        # TODO: substitute with Interval(convex_hull(vertices)...): convex hull 1D
        return Interval(minimum(vertices)[1], maximum(vertices)[1])
    elseif m == 2
        return VPolygon(vertices)
    else
        return VPolytope(vertices)
    end
end

@inline function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolytope{N},
                          use_inv::Bool) where {N<:Real}
    constraints = _linear_map_hrep_helper(M, P, use_inv)
    m = size(M, 1) # output dimension
    if m == 1
        return convert(Interval, HPolygon(constraints))
    elseif m == 2
        return HPolygon(constraints)
    else
        return HPolytope(constraints)
    end
end

# =============================================
# Functions that require Polyhedra
# =============================================

function load_polyhedra_abstractpolytope() # function to be loaded by Requires
return quote
import .Polyhedra: polyhedron
export polyhedron
using .Polyhedra: HRep, VRep,
                  removehredundancy!, removevredundancy!,
                  hreps, vreps,
                  intersect,
                  convexhull,
                  hcartesianproduct, vcartesianproduct,
                  points,
                  default_library

function default_polyhedra_backend(P, N::Type{<:AbstractFloat})
    return default_library(LazySets.dim(P), Float64)
end

function default_polyhedra_backend(P, N::Type{<:Rational})
    return default_library(LazySets.dim(P), Rational{Int})
end
end # quote
end # function load_polyhedra_hpolytope()

ext: jl


text:
#import GDAL

# # GDALRasterize
# function rasterize!(outar,shapefile;x0=-180.0,y0=-90.0,x1=180.0,y1=90.0)
#     sx,sy = size(outar)
#     ds_point = GDAL.openex(shapefile, GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)
#     ds_point == C_NULL && error("Could not open shapefile")
#     options = GDAL.rasterizeoptionsnew(["-of","MEM","-ts",string(sx),string(sy),"-te",string(x0),string(y0),string(x1),string(y1),"-a_srs","WGS84"], C_NULL)
#     ds_rasterize = GDAL.rasterize("data/point-rasterize.mem", Ptr{GDAL.GDALDatasetH}(C_NULL), ds_point, options, C_NULL)
#     GDAL.rasterizeoptionsfree(options)
#     GDAL.close(ds_point)
#
#     band = GDAL.getrasterband(ds_rasterize,1)
#
#     xsize = GDAL.getrasterbandxsize(band)
#     ysize = GDAL.getrasterbandysize(band)
#     dtype = GDAL.getdatatypename(GDAL.getrasterdatatype(band))
#     @assert (xsize,ysize)==(sx,sy)
#     GDAL.rasterio(band, GDAL.GF_Read, 0, 0, xsize, ysize,
#     outar, xsize, ysize, GDAL.GDT_Float64, 0, 0)
#     GDAL.close(ds_rasterize)
# end

import Shapefile
import GeoInterface: AbstractMultiPolygon, AbstractPoint
export cubefromshape
import ..Cubes.Axes: get_bb, axisfrombb
import WeightedOnlineStats: WeightedMean
import Dates: Day

function getlabeldict(shapepath,labelsym,T,labelsleft)
  t = Shapefile.Table(shapepath)
  labels = getproperty(t,labelsym)
  labeldict = Dict(T(i)=>stripc0x(labels[i]) for i in 1:length(labels) if T(i) in labelsleft)
  properties = Dict("labels"=>labeldict)
end
getlabeldict(shapepath, ::Nothing,T, labelsleft)=Dict{String,Any}()

function aggregate_out(allout, highmat, labelsleft,n)
  dsort = Dict(i[2]=>i[1] for i in enumerate(labelsleft))
  for jout in 1:size(allout,2)
    for iout in 1:size(allout,1)
      v = view(highmat,(iout*n-(n-1)):iout*n,(jout*n-(n-1)):jout*n,:)
      for value in skipmissing(v)
        allout[iout,jout,dsort[value]]+=1
      end
    end
  end
  map!(i->i/(n*n),allout,allout)
end

function cubefromshape_fraction(shapepath,lonaxis,lataxis;labelsym=nothing, T=Float64, nincrease=10)
  s = (length(lonaxis), length(lataxis))
  outmat = zeros(T,map(i->i*nincrease,s)...)
  lon1,lon2 = get_bb(lonaxis)
  lat1,lat2 = get_bb(lataxis)
  rasterize!(outmat, shapepath, bb = (left = lon1, right=lon2, top=lat1, bottom=lat2))

  outmat = replace(outmat,zero(T)=>missing)
  labelsleft = collect(skipmissing(unique(outmat)))
  pp = getlabeldict(shapepath,labelsym,T,Set(labelsleft))["labels"]

  allout = zeros(T,s...,length(labelsleft))
  aggregate_out(allout,outmat,labelsleft,nincrease)

  if labelsym===nothing
    newax = CategoricalAxis("Label", labelsleft)
  else
    newax = CategoricalAxis(labelsym,[pp[l] for l in labelsleft])
  end

  return CubeMem(CubeAxis[lonaxis, lataxis, newax], allout)

end
function cubefromshape_single(shapepath, lonaxis, lataxis; labelsym = nothing, T=Int32)
  s = (length(lonaxis), length(lataxis))
  outmat = zeros(T,s...)
  lon1,lon2 = get_bb(lonaxis)
  lat1,lat2 = get_bb(lataxis)
  rasterize!(outmat, shapepath, bb = (left = lon1, right=lon2, top=lat1, bottom=lat2))

  outmat = replace(outmat,zero(T)=>missing)
  labelsleft = collect(skipmissing(unique(outmat)))
  properties = getlabeldict(shapepath,labelsym,T,labelsleft)

  return CubeMem(CubeAxis[lonaxis, lataxis], outmat,properties)
end

"""
    cubefromshape_single(shapepath, refcube; samplefactor=nothing, labelsym = nothing, T=Int32)

This function will read the shapefile at `shapepath` (path to the file with .shp extension), which is required to contain
a list of polygons. The polygons will be *rasterized* to the same grid and bounding box as the reference data cube `refcube`.
If the shapefile comes with additional labels in a .dbf file, the labels can be mapped to the respective field codes by
providing the label to choose as a symbol. If a `samplefactor` is supplied, the polygons will be rasterized on an n-times higher
resolution and the fraction of area within each polygon will be returned.
"""
cubefromshape(shapepath, c::AbstractCubeData; kwargs...) = cubefromshape(shapepath, getAxis("Lon",c), getAxis("Lat",c);kwargs...)
function cubefromshape(args...; samplefactor=nothing, kwargs...)
  if samplefactor===nothing
    cubefromshape_single(args...; kwargs...)
  else
    cubefromshape_fraction(args...; nincrease = samplefactor, kwargs...)
  end
end

function prune_labels!(c::CubeMem)
  if haskey(c.properties,"labels")
    labelsleft = Set(skipmissing(unique(c.data)))
    dold = c.properties["labels"]
    dnew = Dict(k=>dold[k] for k in filter(i->in(i,labelsleft),keys(dold)))
    c.properties["labels"] = dnew
  end
  c
end
stripc0x(a) = replace(a, r"[^\x20-\x7e]"=> "")

function rasterize!(outar,shapefile;bb = (left = -180.0, right=180.0, top=90.0,bottom=-90.0),label=nothing)
  shapepath = shapefile
  t = Shapefile.Table(shapepath)
  p = Shapefile.shapes(t)
  if length(p)>1
    rasterizepoly!(outar,p,bb)
  else
    rasterizepoly!(outar,p[1],bb)
  end
end

function rasterizepoly!(outmat,poly::Vector{<:Union{Missing,AbstractMultiPolygon}},bb)
    foreach(1:length(poly)) do ipoly
        rasterizepoly!(outmat,poly[ipoly],bb,value=ipoly)
    end
    outmat
end

function rasterizepoly!(outmat,poly::Vector{T},bb;value=one(eltype(outmat)), wrap = (left = -180.0,right = 180.0)) where T<:AbstractPoint
nx,ny = size(outmat)
resx = (bb.right-bb.left)/nx
resy = (bb.top-bb.bottom)/ny
xr = range(bb.left+resx/2,bb.right-resx/2,length=nx)
yr = range(bb.top-resy/2,bb.bottom+resy/2,length=ny)
wrapwidth = wrap.right-wrap.left

for (iy,pixelY) in enumerate(yr)
    nodeX = Float64[]
    j = length(poly)
    for i = 1:length(poly)
        p1 = poly[i]
        p2 = poly[j]
        if wrap !==nothing && abs(p1.x-p2.x)>wrapwidth
            p1,p2 = p1.x < p2.x ? (T(p1.x+wrapwidth,p1.y),T(p2.x,p2.y)) : (T(p1.x,p1.y),T(p2.x+wrapwidth,p2.y))
        end
        if (p1.y < pixelY) && (p2.y >= pixelY) || (p2.y < pixelY) && (p1.y >= pixelY)
            push!(nodeX, p1.x + (pixelY-p1.y)/(p2.y-p1.y)*(p2.x-p1.x))
        end
        j = i
    end
    #Add intersect points at start and end for wrapped polygons
    if wrap!==nothing && any(i->i>wrap.right,nodeX)
        push!(nodeX,wrap.left)
        push!(nodeX,wrap.right)
        map!(i->i>wrap.right ? i-wrapwidth : i,nodeX,nodeX)
    end
    sort!(nodeX)
    @assert(iseven(length(nodeX)))
    for i = 1:2:length(nodeX)
        outmat[searchsortedfirst(xr,nodeX[i]):searchsortedlast(xr,nodeX[i+1]),iy].=value
    end
end
    outmat
end

function rasterizepoly!(outmat,pp::Shapefile.Polygon,bb;value=one(eltype(outmat)))
    points = map(1:length(pp.parts)) do ipart
        i0 = pp.parts[ipart]+1
        iend = ipart == length(pp.parts) ? length(pp.points) : pp.parts[ipart+1]
        pp.points[i0:iend]
    end
    foreach(1:length(points)) do ipoly
        rasterizepoly!(outmat,points[ipoly],bb,value=value)
    end
    outmat
end

function getboundingbox(data)
  defaultval = (CartesianIndex(size(data)),CartesianIndex{ndims(data)}())
  r1, r2 = mapreduce(
    i->ismissing(data[i]) ? defaultval : (i,i),
    (x,y)->(min(x[1],y[1]),max(x[2],y[2])),
    CartesianIndices(data),
    init = defaultval
  )
  CartesianIndices(map(Colon(),r1.I,r2.I))
end

"""
    aggregate_shapefile(shp,cube,refdate; npast=3, nfuture=10)

Calculates spatially aggregated statistics over a shapefile specified by `shp`.
Returns a time series cube for variables in `cube` that contains `npast` time steps
prior to and `nfuture` time steps after the event.
"""
function aggstats(shp,csubvar,refdate;npast=3,nfuture=10)

  craster = cubefromshape(shp, csubvar)
  bspace = getboundingbox(craster.data)
  blon,blat = map(i->CartesianIndices((i,)),bspace.indices)
  taxall = getAxis("Time",csubvar)
  itime = axVal2Index(taxall,refdate)
  refdate = taxall.values[itime]
  btime = LinearIndices((itime-npast:itime+nfuture,))
  csmallbox = csubvar[lon=blon,lat=blat,time=btime]
  maskcube = craster[lon=blon,lat=blat]
  if findAxis("Variable",csubvar)===nothing
    by = ()
    incax = ("lat","time")
  else
    by = (:var,)
    incax = ("lat","var","time")
  end
  ct = CubeTable(cdata = csmallbox, cmask = maskcube, include_axes=incax)

  by = (i->iszero(i.cmask),by...,i->Int((i.time-refdate)/Day(1)))
  fitres = cubefittable(ct,WeightedMean,:cdata,by=by,weight=i->cosd(i.lat),showprog=false)

  fr = fitres[Category1=false]
  renameaxis!(fr,"Category"=>RangeAxis("Time",(-npast:nfuture)*8))
  fr
end


text:
#=
  Bertrand's paradox

  From https://towardsdatascience.com/five-paradoxes-with-probabilities-that-will-puzzle-you-2f71201d6ee8
  """
  3. Bertrand's Box Paradox

  If you are familiar with the Montey Hall problem, this paradox is quite
  similar. In front of us there are three boxes:

    [image of three boxes]

  One box contains two silver coins, one box contains two gold coins and
  one box contains a gold- and a silver coin. We do not know which coins
  are in which box. Now, we pick a random box and blindly draw a coin from
  our box. It’s a gold coin!

  Now, the question is:

    What’s the probability that the second coin in our box is also a gold coin?

  My naive (and wrong) answer when encountering the problem for the first
  time was ½. I thought that because we drew a gold coin, our box is either
  the one with the two gold coins or the one with the mixed coins. In the
  first case, we would draw another gold coin and in the second case, we
  wouldn’t. Therefore, I presumed the probability should be ½.

    The real probability is ⅔.

  The reason for that is that the first gold coin we drew could either be
  the only gold coin in the mixed box, the first gold coin in the solely
  golden box, or the second gold coin in the solely golden box. And in
  two of these three possibilities, we will draw another gold coin.

  """

  Here we define the boxes as:
  - box1: 2 gold coins
  - box2: 1 gold and 1 silver coin
  - box3: 2 silver coins

  Distributions of variable coin1
  gold       =>    1000  (1.000000)
  Distributions of variable coin2
  gold       =>     668  (0.668000)
  silver     =>     332  (0.332000)
  Distributions of variable box
  b1         =>     668  (0.668000)
  b2         =>     332  (0.332000)

  I.e. probability that the second coin is also gold is 0.668000 (~ 2/3).

  Cf ~/blog/bertrands_paradox.blog
      ~/psi/bertrands_paradox.psi

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

# We observe that the first coin is a gold coin
@model function bertrands_paradox(coin1=1)
    # There are 3 boxes.
    b1 = 1
    b2 = 2
    b3 = 3

    # We select one of the boxes
    box ~ DiscreteUniform(1,3)

    # the coins
    gold = 1
    silver = 2

    # We draw one coin:
    coin1 ~ (box == b1) ? Categorical(simplex([2.0, 0.0])) :
            (box == b2) ? Categorical(simplex([1, 1])) :
            Categorical(simplex([0.0, 2.0]))

    # It's a gold coin.
    # true ~ Dirac(coin1 == gold)

    # What is the probability that the second coin
    # in the selected box is also a gold coin?
    coin2 ~ (box == b1) ? Categorical(simplex([1.0, 0.0])) :
            (box == b2) ? Categorical(simplex([0.0, 1.0])) :
            Categorical(simplex([0.0, 1.0]))

end

model = bertrands_paradox()
num_chns = 4

# chns = sample(model, Prior(), 10_000, num_chns)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)

show_var_dist_pct(chns, :coin1,["gold","silver"])
show_var_dist_pct(chns, :coin2,["gold","silver"])
show_var_dist_pct(chns, :box,["b1","b2","b3"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using CSV
using Dates
using DataFrames
using DelimitedFiles
using XLSX
using YAML
using InteractiveUtils
const IU = InteractiveUtils

using SLiDE  # see src/SLiDE.jl

READ_DIR = joinpath("data", "readfiles")

files_map = [
    XLSXInput("generate_yaml.xlsx", "map_parse",     "B1:Z150", "map_parse"),
    XLSXInput("generate_yaml.xlsx", "map_scale",     "B1:Z150", "map_scale"),
    XLSXInput("generate_yaml.xlsx", "map_bluenote",  "B1:Z150", "map_bluenote"),
    XLSXInput("generate_yaml.xlsx", "map_crosswalk", "B1:Z150", "map_crosswalk")
]

files_map = write_yaml(READ_DIR, files_map)
files_map = run_yaml(files_map)
include(joinpath(SLIDE_DIR, "dev", "datastream", "adjust_maps.jl"))

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using DelimitedFiles, Printf
using Interpolations, Plots

using Oceananigans
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.Utils

#####
##### Some useful constants
#####

const km = 1000
const Ω_Earth = 7.292115e-5  # [s⁻¹]
const φ = -75  # degrees latitude

#####
##### Model grid and domain size
#####

arch = CPU()
FT = Float64

Nx = 32
Ny = 32
Nz = 32

Lx = 10km
Ly = 10km
Lz = 1km

end_time = 7day

#####
##### Set up source of meltwater: We will implement a source of meltwater as
##### a relaxation term towards a reference T and S value at a single point.
##### This is in effect weakly imposing a Value/Dirchlet boundary condition.
#####

const source_type = :point
# const source_type = :line

λ = 1/(1minute)  # Relaxation timescale [s⁻¹].

# Temperature and salinity of the meltwater outflow.
T_source = -1
S_source = 33.95

# Index of the point source at the middle of the southern wall.
source_index = (Int(Nx/2), 1, Int(Nz/2))

# Point source
@inline T_point_source(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse((i, j, k) == p.source_index, -p.λ * (C.T[i, j, k] - p.T_source), 0)

@inline S_point_source(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse((i, j, k) == p.source_index, -p.λ * (C.S[i, j, k] - p.S_source), 0)

# Line source
@inline T_line_source(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse((j, k) == (p.source_index[2], p.source_index[3]), -p.λ * (C.T[i, j, k] - p.T_source), 0)

@inline S_line_source(i, j, k, grid, time, U, C, p) =
    @inbounds ifelse((j, k) == (p.source_index[2], p.source_index[3]), -p.λ * (C.S[i, j, k] - p.S_source), 0)

params = (source_index=source_index, T_source=T_source, S_source=S_source, λ=λ)

if source_type == :point
    forcing = ModelForcing(T = T_point_source, S = S_point_source)
elseif source_type == :line
    forcing = ModelForcing(T = T_line_source, S = S_line_source)
end

#####
##### Set up model
#####

# eos = LinearEquationOfState()
eos = RoquetIdealizedNonlinearEquationOfState(:freezing)

model = Model(
           architecture = arch,
             float_type = FT,
                   grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(-Lx/2, Lx/2), y=(0, Ly), z=(-Lz, 0)),
                tracers = (:T, :S, :meltwater),
               coriolis = FPlane(rotation_rate=Ω_Earth, latitude=φ),
               buoyancy = SeawaterBuoyancy(equation_of_state=eos),
                closure = AnisotropicMinimumDissipation(),
    boundary_conditions = ChannelSolutionBCs(),
                forcing = forcing,
             parameters = params
)

#####
##### Read reference profiles from disk
##### As these profiles are derived from observations, we will have to do some
##### post-processing to be able to use them as initial conditions.
#####
##### We will get rid of all NaN values and use the remaining data to linearly
##### interpolate the T and S profiles to the model's grid.
#####

# The pressure is given in dbar so we will convert to depth (meters) assuming
# 1 dbar = 1 meter (this is approximately true).
z = readdlm("reference_pressure.txt")[:]

# We also flatten the arrays by indexing with a Colon [:] to convert the arrays
# from N×1 arrays to 1D arrays of length N.
T = readdlm("reference_temperature.txt")[:]
S = readdlm("reference_salinity.txt")[:]

# Get the indices of all the non-NaN values.
T_good_inds = findall(!isnan, T)
S_good_inds = findall(!isnan, S)

# Create T and S arrays that do not contain NaNs, along with corresponding
# z values.
T_good = T[T_good_inds]
S_good = S[S_good_inds]

z_T = z[T_good_inds]
z_S = z[S_good_inds]

# Linearly interpolate T and S profiles to model grid.
Ti = LinearInterpolation(z_T, T_good, extrapolation_bc=Interpolations.Flat())
Si = LinearInterpolation(z_S, S_good, extrapolation_bc=Interpolations.Flat())

zC = model.grid.zC
T₀ = Ti.(-zC)
S₀ = Si.(-zC)

# Plot and save figures of reference and interpolated profiles.
T_fpath = "temperature_profiles.png"
S_fpath = "salinity_profiles.png"

T_plot = plot(T_good, -z_T, grid=false, dpi=300, label="Reference",
              xlabel="Temperature (C)", ylabel="Depth (m)")
plot!(T_plot, T₀, zC, label="Interpolation")
@info "Saving temperature profiles to $T_fpath..."
savefig(T_plot, T_fpath)

S_plot = plot(S_good, -z_S, grid=false, dpi=300, label="Reference",
              xlabel="Salinity (ppt)", ylabel="Depth (m)", )
plot!(S_plot, S₀, zC, label="Interpolation")
@info "Saving temperature profiles to $S_fpath..."
savefig(S_plot, S_fpath)

#####
##### Setting up initial conditions
#####

T₀_3D = repeat(reshape(T₀, 1, 1, Nz), Nx, Ny, 1)
S₀_3D = repeat(reshape(S₀, 1, 1, Nz), Nx, Ny, 1)

set!(model.tracers.T, T₀_3D)
set!(model.tracers.S, S₀_3D)

# Set meltwater concentration to 1 at the source.
if source_type == :point
    model.tracers.meltwater.data[source_index...] = 1
elseif source_type == :line
    model.tracers.meltwater.data[:, source_index[2], source_index[3]] .= 1  # Line source
end

#####
##### Write out 3D fields and slices to NetCDF files.
#####

fields = Dict(
        "u" => model.velocities.u,
        "v" => model.velocities.v,
        "w" => model.velocities.w,
        "T" => model.tracers.T,
        "S" => model.tracers.S,
"meltwater" => model.tracers.meltwater,
       "nu" => model.diffusivities.νₑ,
   "kappaT" => model.diffusivities.κₑ.T,
   "kappaS" => model.diffusivities.κₑ.S
)

output_attributes = Dict(
        "u" => Dict("longname" => "Velocity in the x-direction", "units" => "m/s"),
        "v" => Dict("longname" => "Velocity in the y-direction", "units" => "m/s"),
        "w" => Dict("longname" => "Velocity in the z-direction", "units" => "m/s"),
        "T" => Dict("longname" => "Temperature", "units" => "C"),
        "S" => Dict("longname" => "Salinity", "units" => "g/kg"),
"meltwater" => Dict("longname" => "Meltwater concentration"),
       "nu" => Dict("longname" => "Nonlinear LES viscosity", "units" => "m^2/s"),
   "kappaT" => Dict("longname" => "Nonlinear LES diffusivity for temperature", "units" => "m^2/s"),
   "kappaS" => Dict("longname" => "Nonlinear LES diffusivity for salinity", "units" => "m^2/s")
)

eos_name(::LinearEquationOfState) = "LinearEOS"
eos_name(::RoquetIdealizedNonlinearEquationOfState) = "RoquetEOS"
prefix = "ice_shelf_meltwater_outflow_$(source_type)_$(eos_name(eos))_"

model.output_writers[:fields] =
    NetCDFOutputWriter(model, fields, filename = prefix * "fields.nc",
                       interval = 6hour, output_attributes = output_attributes)

model.output_writers[:depth_slice] =
    NetCDFOutputWriter(model, fields, filename = prefix * "middepth_xy_slice.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       zC = source_index[3], zF = source_index[3])

model.output_writers[:surface_slice] =
    NetCDFOutputWriter(model, fields, filename = prefix * "surface_xy_slice.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       zC = Nz, zF = Nz)

model.output_writers[:calving_front_slice] =
    NetCDFOutputWriter(model, fields, filename = prefix * "calving_front_xz_slice.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       yC = 1, yF = 2)

model.output_writers[:along_channel_slice] =
    NetCDFOutputWriter(model, fields, filename = prefix * "along_channel_yz_slice.nc",
                       interval = 5minute, output_attributes = output_attributes,
                       xC = source_index[1], xF = source_index[1])

#####
##### Print banner
#####

@printf("""

    Simulating ocean dynamics of meltwater outflow from beneath Antarctic ice shelves
        N : %d, %d, %d
        L : %.3g, %.3g, %.3g [km]
        Δ : %.3g, %.3g, %.3g [m]
        φ : %.3g [latitude]
        f : %.3e [s⁻¹]
     days : %d
   source : %s
 T_source : %.2f [°C]
 S_source : %.2f [g/kg]
  closure : %s
      EoS : %s

""", model.grid.Nx, model.grid.Ny, model.grid.Nz,
     model.grid.Lx / km, model.grid.Ly / km, model.grid.Lz / km,
     model.grid.Δx, model.grid.Δy, model.grid.Δz,
     φ, model.coriolis.f, end_time / day,
     source_type, T_source, S_source,
     typeof(model.closure), typeof(model.buoyancy.equation_of_state))

#####
##### Time step!
#####

# Wizard utility that calculates safe adaptive time steps.
wizard = TimeStepWizard(cfl=0.3, Δt=1second, max_change=1.2, max_Δt=30second)

# CFL utilities for reporting stability criterions.
cfl = AdvectiveCFL(wizard)
dcfl = DiffusiveCFL(wizard)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 50

# Convinient alias
C_mw = model.tracers.meltwater

while model.clock.time < end_time
    walltime = @elapsed begin
        time_step!(model; Nt=Ni, Δt=wizard.Δt)

        if source_type == :point
            C_mw.data[source_index...] = 1
        elseif source_type == :line
            C_mw.data[:, source_index[2], source_index[3]] .= 1
        end

        # Normalize meltwater concentration to be 0 <= C_mw <= 1.
        C_mw.data .= max.(0, C_mw.data)
        C_mw.data .= C_mw.data ./ maximum(C_mw.data)
    end

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Find maximum velocities.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    # Find maximum ν and κ.
    νmax = maximum(model.diffusivities.νₑ.data.parent)
    κmax = maximum(model.diffusivities.κₑ.T.data.parent)

    # Calculate a new adaptive time step.
    update_Δt!(wizard, model)

    # Print progress statement.
    i, t = model.clock.iteration, model.clock.time
    @printf("[%06.2f%%] i: %d, t: %5.2f days, umax: (%6.3g, %6.3g, %6.3g) m/s, CFL: %6.4g, νκmax: (%6.3g, %6.3g), νκCFL: %6.4g, next Δt: %8.5g s, ⟨wall time⟩: %s\n",
            progress, i, t / day, umax, vmax, wmax, cfl(model), νmax, κmax, dcfl(model), wizard.Δt, prettytime(walltime / Ni))
end

for ow in model.output_writers
	ow isa NetCDFOutputWriter && close(ow)
end


text:
using Compose, Viznet

set_default_graphic_size(12cm, 12cm)

dot = nodestyle(:circle;r = 0.005)
diamond = nodestyle(:diamond)
tri1 = nodestyle(:triangle; θ=-π/2)
tri2 = nodestyle(:triangle; θ=π/2)
sq = nodestyle(:square; r=0.015)
e1 = bondstyle(:default, stroke("black"))
lt = SquareLattice(10, 10)

function vnode(node)
    canvas() do
        lt = SquareLattice(6, 6)
        for v in vertices(lt)
            node >> lt[v]
        end
    end
end

function circuit_body(lt)
    for i in vertices(lt)
        dot >> lt[i]
    end
    for (i, j) in bonds(lt)
        center = (lt[i] .+ lt[j]) ./ 2
        (lt[i][1]!=lt[j][1] ? sq : diamond) >> center
    end
    for i=1:lt.Nx
        e1 >> lt[(i,1);(i,lt.Ny)]
    end
    for j=1:lt.Ny
        e1 >> lt[(1,j);(lt.Nx,j)]
    end
end

function circuit_ends(lt)
    offset = 0.7
    for j=1:lt.Ny
        e1 >> lt[(1-offset,j);(lt.Nx+offset,j)]
        tri1 >> lt[1-offset,j]
        tri2 >> lt[lt.Nx+offset,j]
    end
end

function circuit_body2(lt)
    offset = 0.15
    for j=1:lt.Ny-1
        i = 1
        _offset = j % 2 == 0 ? offset : -offset
        e1 >> lt[(i+_offset,j); (i+_offset, j+1)]
        diamond >> (lt[i+_offset,j]./2 .+ lt[i+_offset, j+1]./2)
    end

    for i=2:lt.Nx
        for j=1:lt.Ny
            _offset = j % 2 == 0 ? offset : -offset
            sq >> (lt[i-1,j]./2 .+ lt[i, j]./2)
        end
        for j=1:lt.Ny-1
            _offset = j % 2 == 0 ? offset : -offset
            e1 >> lt[(i+_offset,j); (i+_offset, j+1)]
            diamond >> (lt[i+_offset,j]./2 .+ lt[i+_offset, j+1]./2)
        end
    end
end

compose(context(0.1,0.1,0.8,0.8), canvas() do
    circuit_body2(lt)
    circuit_ends(lt)
end)


text:
methodid(args...) = hash(args)
methodid(ir::IR) = methodid(IRTools.argtypes(ir)...)

const MethodId = UInt

"Location in trace"
struct Loc{V}
  mthdid::MethodId
  statement::V
end


text:
###
### Particle Filtering and Particle MCMC Samplers.
###

#######################
# Particle Transition #
#######################

"""
    ParticleTransition{T, F<:AbstractFloat} <: AbstractTransition

Fields:
- `θ`: The parameters for any given sample.
- `lp`: The log pdf for the sample's parameters.
- `le`: The log evidence retrieved from the particle.
- `weight`: The weight of the particle the sample was retrieved from.
"""
struct ParticleTransition{T, F<:AbstractFloat} <: AbstractTransition
    θ::T
    lp::F
    le::F
    weight::F
end

transition_type(spl::Sampler{<:ParticleInference}) = ParticleTransition

function additional_parameters(::Type{<:ParticleTransition})
    return [:lp,:le, :weight]
end

####
#### Generic Sequential Monte Carlo sampler.
####

"""
    SMC()

Sequential Monte Carlo sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
SMC()
```
"""
struct SMC{space, R} <: ParticleInference
    resampler::R
end

function SMC(resampler = Turing.Core.ResampleWithESSThreshold(), space::Tuple = ())
    SMC{space, typeof(resampler)}(resampler)
end
SMC(::Tuple{}) = SMC()
SMC(space::Symbol...) = SMC(space)

mutable struct SMCState{V<:VarInfo, F<:AbstractFloat} <: AbstractSamplerState
    vi                   ::   V
    # The logevidence after aggregating all samples together.
    average_logevidence  ::   F
    particles            ::   ParticleContainer
end

function SMCState(model::Model)
    vi = VarInfo(model)
    particles = ParticleContainer(model, Trace[])

    return SMCState(vi, 0.0, particles)
end

function Sampler(alg::SMC, model::Model, s::Selector)
    dict = Dict{Symbol, Any}()
    state = SMCState(model)
    return Sampler(alg, dict, s, state)
end

function sample_init!(
    ::AbstractRNG,
    model::Turing.Model,
    spl::Sampler{<:SMC},
    N::Integer;
    kwargs...
)
    # set the parameters to a starting value
    initialize_parameters!(spl; kwargs...)

    # reset the VarInfo
    vi = spl.state.vi
    vi.num_produce = 0
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)
    empty!(vi)

    # create a new set of particles
    T = Trace{typeof(spl),typeof(vi),typeof(model)}
    particles = T[Trace(model, spl, vi) for _ in 1:N]

    # create a new particle container
    spl.state.particles = pc = ParticleContainer(model, particles)

    while consume(pc) !== Val{:done}
        resample!(pc, spl.alg.resampler)
    end
end

function step!(
    ::AbstractRNG,
    model::Turing.Model,
    spl::Sampler{<:SMC},
    ::Integer;
    iteration=-1,
    kwargs...
)
    # check that we received a real iteration number
    @assert iteration >= 1 "step! needs to be called with an 'iteration' keyword argument."

    # grab the weights
    pc = spl.state.particles
    Ws = getweights(pc)

    # update the master vi
    particle = pc.vals[iteration]
    params = tonamedtuple(particle.vi)
    lp = getlogp(particle.vi)

    return ParticleTransition(params, lp, pc.logE, Ws[iteration])
end

####
#### Particle Gibbs sampler.
####

"""
    PG(n_particles::Int)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
PG(100, 100)
```
"""
struct PG{space} <: ParticleInference
  n_particles           ::    Int         # number of particles used
  resampler             ::    Function    # function to resample
end
function PG(n_particles::Int, resampler::Function, space::Tuple)
    return PG{space}(n_particles, resampler)
end
PG(n1::Int, ::Tuple{}) = PG(n1)
function PG(n1::Int, space::Symbol...)
    PG(n1, resample_systematic, space)
end

mutable struct PGState{V<:VarInfo, F<:AbstractFloat} <: AbstractSamplerState
    vi                   ::   V
    # The logevidence after aggregating all samples together.
    average_logevidence  ::   F
end

function PGState(model::Model)
    vi = VarInfo(model)
    return PGState(vi, 0.0)
end

const CSMC = PG # type alias of PG as Conditional SMC

"""
    Sampler(alg::PG, model::Model, s::Selector)

Return a `Sampler` object for the PG algorithm.
"""
function Sampler(alg::PG, model::Model, s::Selector)
    info = Dict{Symbol, Any}()
    state = PGState(model)
    return Sampler(alg, info, s, state)
end

function step!(
    ::AbstractRNG,
    model::Turing.Model,
    spl::Sampler{<:PG},
    ::Integer;
    kwargs...
)
    # obtain or create reference particle
    vi = spl.state.vi
    ref_particle = isempty(vi) ? nothing : forkr(Trace(model, spl, vi))

    # reset the VarInfo before new sweep
    vi.num_produce = 0
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)

    # create a new set of particles
    num_particles = spl.alg.n_particles
    T = Trace{typeof(spl),typeof(vi),typeof(model)}
    if ref_particle === nothing
        particles = T[Trace(model, spl, vi) for _ in 1:num_particles]
    else
        particles = Vector{T}(undef, num_particles)
        @inbounds for i in 1:(num_particles - 1)
            particles[i] = Trace(model, spl, vi)
        end
        @inbounds particles[num_particles] = ref_particle
    end

    # create a new particle container
    pc = ParticleContainer(model, particles)

    # run the particle filter
    while consume(pc) !== Val{:done}
        resample!(pc, spl.alg.resampler, ref_particle)
    end

    # pick a particle to be retained.
    Ws = getweights(pc)
    indx = randcat(Ws)

    # extract the VarInfo from the retained particle.
    params = tonamedtuple(vi)
    spl.state.vi = pc.vals[indx].vi
    lp = getlogp(spl.state.vi)

    # update the master vi.
    return ParticleTransition(params, lp, pc.logE, 1.0)
end

function sample_end!(
    ::AbstractRNG,
    ::Model,
    spl::Sampler{<:ParticleInference},
    N::Integer,
    ts::Vector{ParticleTransition};
    kwargs...
)
    # Set the default for resuming the sampler.
    resume_from = get(kwargs, :resume_from, nothing)

    # Exponentiate the average log evidence.
    # loge = exp(mean([t.le for t in ts]))
    loge = mean(t.le for t in ts)

    # If we already had a chain, grab the logevidence.
    if resume_from !== nothing   # concat samples
        @assert resume_from isa Chains "resume_from needs to be a Chains object."
        # pushfirst!(samples, resume_from.info[:samples]...)
        pre_loge = resume_from.logevidence
        # Calculate new log-evidence
        pre_n = length(resume_from)
        loge = (pre_loge * pre_n + loge * N) / (pre_n + N)
    end

    # Store the logevidence.
    spl.state.average_logevidence = loge
end

function assume(spl::Sampler{<:Union{PG,SMC}}, dist::Distribution, vn::VarName, ::VarInfo)
    vi = current_trace().vi
    if vn in getspace(spl)
        if ~haskey(vi, vn)
            r = rand(dist)
            push!(vi, vn, r, dist, spl)
        elseif is_flagged(vi, vn, "del")
            unset_flag!(vi, vn, "del")
            r = rand(dist)
            vi[vn] = vectorize(dist, r)
            setgid!(vi, spl.selector, vn)
            setorder!(vi, vn, vi.num_produce)
        else
            updategid!(vi, vn, spl)
            r = vi[vn]
        end
    else # vn belongs to other sampler <=> conditionning on vn
        if haskey(vi, vn)
            r = vi[vn]
        else
            r = rand(dist)
            push!(vi, vn, r, dist, Selector(:invalid))
        end
        acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    end
    return r, 0
end

function observe(spl::Sampler{<:Union{PG,SMC}}, dist::Distribution, value, vi)
    produce(logpdf(dist, value))
    return 0
end

####
#### Resampling schemes for particle filters
####

# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# Default resampling scheme
function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
    return resample_systematic(w, num_particles)
end

# More stable, faster version of rand(Categorical)
function randcat(p::AbstractVector{<:Real})
    T = eltype(p)
    r = rand(T)
    s = 1
    for j in eachindex(p)
        r -= p[j]
        if r <= zero(T)
            s = j
            break
        end
    end
    return s
end

function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
    return rand(Distributions.sampler(Categorical(w)), num_particles)
end

function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)

    M = length(w)

    # "Repetition counts" (plus the random part, later on):
    Ns = floor.(length(w) .* w)

    # The "remainder" or "residual" count:
    R = Int(sum(Ns))

    # The number of particles which will be drawn stocastically:
    M_rdn = num_particles - R

    # The modified weights:
    Ws = (M .* w - floor.(M .* w)) / M_rdn

    # Draw the deterministic part:
    indx1, i = Array{Int}(undef, R), 1
    for j in 1:M
        for k in 1:Ns[j]
            indx1[i] = j
            i += 1
        end
    end

    # And now draw the stocastic (Multinomial) part:
    return append!(indx1, rand(Distributions.sampler(Categorical(w)), M_rdn))
end

"""
    resample_stratified(weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by stratified resampling.

In stratified resampling `n` ordered random numbers `u₁`, ..., `uₙ` are generated, where
``uₖ \\sim U[(k - 1) / n, k / n)``. Based on these numbers the samples `x₁`, ..., `xₙ`
are selected according to the multinomial distribution defined by the normalized `weights`,
i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_stratified(weights::AbstractVector{<:Real}, n::Integer)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]

    # generate all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # sample next `u` (scaled by `n`)
        u = oftype(v, i - 1 + rand())

        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample
    end

    return samples
end

"""
    resample_systematic(weights, n)

Return a vector of `n` samples `x₁`, ..., `xₙ` from the numbers 1, ..., `length(weights)`,
generated by systematic resampling.

In systematic resampling a random number ``u \\sim U[0, 1)`` is used to generate `n` ordered
numbers `u₁`, ..., `uₙ` where ``uₖ = (u + k − 1) / n``. Based on these numbers the samples
`x₁`, ..., `xₙ` are selected according to the multinomial distribution defined by the
normalized `weights`, i.e., `xᵢ = j` if and only if
``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``.
"""
function resample_systematic(weights::AbstractVector{<:Real}, n::Integer)
    # check input
    m = length(weights)
    m > 0 || error("weight vector is empty")

    # pre-calculations
    @inbounds v = n * weights[1]
    u = oftype(v, rand())

    # find all samples
    samples = Array{Int}(undef, n)
    sample = 1
    @inbounds for i in 1:n
        # as long as we have not found the next sample
        while v < u
            # increase and check the sample
            sample += 1
            sample > m &&
                error("sample could not be selected (are the weights normalized?)")

            # update the cumulative sum of weights (scaled by `n`)
            v += n * weights[sample]
        end

        # save the next sample
        samples[i] = sample

        # update `u`
        u += one(u)
    end

    return samples
end


text:
# We put the definition of the Cassette context in its own file to
# prevent Revise from re-evaluating it
Cassette.@context CassetteCtx


text:
@block SimonDanisch [layout] begin

    @cell "Layouting" [scatter, lines, surface, heatmap, vbox] begin
        p1 = scatter(rand(10), markersize = 1)
        p2 = lines(rand(10), rand(10))
        p3 = surface(0..1, 0..1, rand(100, 100))
        p4 = heatmap(rand(100, 100))
        x = 0:0.1:10
        p5 = lines(0:0.1:10, sin.(x))
        pscene = vbox(
            hbox(p1, p2),
            p3,
            hbox(p4, p5, sizes = [0.7, 0.3]),
            sizes = [0.2, 0.6, 0.2]
        )
    end

    @cell "Comparing contours, image, surfaces and heatmaps" [image, contour, surface, heatmap, vbox] begin
        N = 20
        x = LinRange(-0.3, 1, N)
        y = LinRange(-1, 0.5, N)
        z = x .* y'
        hbox(
            vbox(
                contour(x, y, z, levels = 20, linewidth =3),
                contour(x, y, z, levels = 0, linewidth = 0, fillrange = true),
                heatmap(x, y, z),
            ),
            vbox(
                image(x, y, z, colormap = :viridis),
                surface(x, y, fill(0f0, N, N), color = z, shading = false),
                image(-0.3..1, -1..0.5, AbstractPlotting.logo())
            )
        )
    end
end

@block JuliusKrumbiegel ["layout"] begin
    using MakieLayout

    @cell "Faceting" [faceting, grid] begin
        # layoutscene is a convenience function that creates a Scene and a GridLayout
        # that are already connected correctly and with Outside alignment
        scene, layout = layoutscene(30, resolution = (1200, 900))

        using ColorSchemes

        ncols = 4
        nrows = 4

        # create a grid of LAxis objects
        axes = [LAxis(scene) for i in 1:nrows, j in 1:ncols]
        # and place them into the layout
        layout[1:nrows, 1:ncols] = axes

        # link x and y axes of all LAxis objects
        linkxaxes!(axes...)
        linkyaxes!(axes...)

        lineplots = [lines!(axes[i, j], (1:0.1:8pi) .+ i, sin.(1:0.1:8pi) .+ j,
                color = get(ColorSchemes.rainbow, ((i - 1) * nrows + j) / (nrows * ncols)), linewidth = 4)
            for i in 1:nrows, j in 1:ncols]

        for i in 1:nrows, j in 1:ncols
            # remove unnecessary decorations in some of the facets, this will have an
            # effect on the layout as the freed up space will be used to make the axes
            # bigger
            i > 1 && (axes[i, j].titlevisible = false)
            j > 1 && (axes[i, j].ylabelvisible = false)
            j > 1 && (axes[i, j].yticklabelsvisible = false)
            j > 1 && (axes[i, j].yticksvisible = false)
            i < nrows && (axes[i, j].xticklabelsvisible = false)
            i < nrows && (axes[i, j].xticksvisible = false)
            i < nrows && (axes[i, j].xlabelvisible = false)
        end

        autolimits!(axes[1]) # hide

        legend = LLegend(scene, permutedims(lineplots, (2, 1)), ["Line $i" for i in 1:length(lineplots)],
            ncols = 2)
        # place a legend on the side by indexing into one column after the current last
        layout[:, end+1] = legend

        # index into the 0th row, thereby adding a new row into the layout and place
        # a text object across the first four columns as a super title
        layout[0, 1:4] = LText(scene, text="MakieLayout Facets", textsize=50)

        scene

    end

    @cell "ManualTicks" [layout, ticks] begin
        scene = Scene(resolution = (1000, 1000));
        campixel!(scene);

        maingl = GridLayout(scene, alignmode = Outside(30))

        la = maingl[1, 1] = LAxis(scene)
        la.attributes.yautolimitmargin = (0f0, 0.05f0)

        poly!(la, BBox(0, 1, 4, 0), color=:blue)
        poly!(la, BBox(1, 2, 7, 0), color=:red)
        poly!(la, BBox(2, 3, 1, 0), color=:green)

        la.attributes.xticks = ManualTicks([0.5, 1.5, 2.5], ["blue", "red", "green"])
        la.attributes.xlabel = "Color"
        la.attributes.ylabel = "Value"

        scene
    end

    @cell "Protrusion changes" [protrusions] begin
        using MakieLayout.Animations

        scene, layout = layoutscene(resolution = (600, 600))

        axes = [LAxis(scene) for i in 1:2, j in 1:2]
        layout[1:2, 1:2] = axes

        a_title = Animation([0, 2], [30.0, 50.0], sineio(n=2, yoyo=true, prewait=0.2))
        a_xlabel = Animation([2, 4], [20.0, 40.0], sineio(n=2, yoyo=true, prewait=0.2))
        a_ylabel = Animation([4, 6], [20.0, 40.0], sineio(n=2, yoyo=true, prewait=0.2))

        record(scene, @replace_with_a_path(mp4), 0:1/60:6, framerate = 60) do t
            axes[1, 1].titlesize = a_title(t)
            axes[1, 1].xlabelsize = a_xlabel(t)
            axes[1, 1].ylabelsize = a_ylabel(t)
        end
    end

    @cell "Hiding decorations" [decorations] begin
        scene = Scene(resolution = (600, 600), camera=campixel!)

        layout = GridLayout(
            scene, 2, 2, # we need to specify rows and columns so the gap sizes don't get lost
            addedcolgaps = Fixed(0),
            addedrowgaps = Fixed(0),
            alignmode = Outside(30)
        )

        axes = [LAxis(scene) for j in 1:2, i in 1:2]
        layout[1:2, 1:2] = axes

        re = record(scene, @replace_with_a_path(mp4), framerate=3) do io
            recordframe!(io)
            for ax in axes
                ax.titlevisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.xlabelvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.ylabelvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.xticklabelsvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.yticklabelsvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.xticksvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.yticksvisible = false
                recordframe!(io)
            end
            for ax in axes
                ax.bottomspinevisible = false
                ax.leftspinevisible = false
                ax.topspinevisible = false
                ax.rightspinevisible = false
                recordframe!(io)
            end
        end
        return re
    end

    @cell "Axis aspects" [axis, aspect] begin
        using FileIO
        scene, layout = layoutscene(30, resolution = (1200, 900))

        axes = [LAxis(scene) for i in 1:2, j in 1:3]
        tightlimits!.(axes)
        layout[1:2, 1:3] = axes

        img = rotr90(MakieGallery.loadasset("cow.png"))

        for ax in axes
            image!(ax, img)
        end

        axes[1, 1].title = "Default"

        axes[1, 2].title = "DataAspect"
        axes[1, 2].aspect = DataAspect()

        axes[1, 3].title = "AxisAspect(418/348)"
        axes[1, 3].aspect = AxisAspect(418/348)

        axes[2, 1].title = "AxisAspect(1)"
        axes[2, 1].aspect = AxisAspect(1)

        axes[2, 2].title = "AxisAspect(2)"
        axes[2, 2].aspect = AxisAspect(2)

        axes[2, 3].title = "AxisAspect(0.5)"
        axes[2, 3].aspect = AxisAspect(0.5)

        scene
    end

    @cell "Linked axes" [axis, link] begin
        scene, layout = layoutscene()

        layout[1, 1:3] = axs = [LAxis(scene) for i in 1:3]
        linkxaxes!(axs[1:2]...)
        linkyaxes!(axs[2:3]...)

        axs[1].title = "x linked"
        axs[2].title = "x & y linked"
        axs[3].title = "y linked"

        for i in 1:3
            lines!(axs[i], 1:10, 1:10, color = "green")
            if i != 1
                lines!(axs[i], 1:10, 11:20, color = "blue")
            end
            if i != 3
                lines!(axs[i], 11:20, 1:10, color = "red")
            end
        end

        scene
    end

    @cell "Nested grids" [grid] begin
        scene, layout = layoutscene(30, resolution = (1200, 900))

        subgl_left = GridLayout()
        subgl_left[1:2, 1:2] = [LAxis(scene) for i in 1:2, j in 1:2]

        subgl_right = GridLayout()
        subgl_right[1:3, 1] = [LAxis(scene) for i in 1:3]

        layout[1, 1] = subgl_left
        layout[1, 2] = subgl_right

        scene
    end

    @cell "Grid alignment" [grid] begin
        scene, layout = layoutscene(30, resolution = (1200, 1200))

        layout[1, 1] = LAxis(scene, title="No grid layout")
        layout[2, 1] = LAxis(scene, title="No grid layout")
        layout[3, 1] = LAxis(scene, title="No grid layout")

        subgl_1 = layout[1, 2] = GridLayout(alignmode=Inside())
        subgl_2 = layout[2, 2] = GridLayout(alignmode=Outside())
        subgl_3 = layout[3, 2] = GridLayout(alignmode=Outside(50))

        subgl_1[1, 1] = LAxis(scene, title="Inside")
        subgl_2[1, 1] = LAxis(scene, title="Outside")
        subgl_3[1, 1] = LAxis(scene, title="Outside(50)")

        layout[1:3, 2] = [LRect(scene, color = :transparent, strokecolor = :red) for i in 1:3]

        scene
    end

    @cell "Spanned grid content" [grid] begin

        scene, layout = layoutscene(4, 4, 30, resolution = (1200, 1200))

        layout[1, 1:2] = LAxis(scene, title="[1, 1:2]")
        layout[2:4, 1:2] = LAxis(scene, title="[2:4, 1:2]")
        layout[:, 3] = LAxis(scene, title="[:, 3]")
        layout[1:3, end] = LAxis(scene, title="[1:3, end]")
        layout[end, end] = LAxis(scene, title="[end, end]")

        scene
    end

    @cell "Indexing outside grid" [grid] begin
        scene, layout = layoutscene(30, resolution = (1200, 1200))

        layout[1, 1] = LAxis(scene)
        for i in 1:3
            layout[:, end+1] = LAxis(scene)
            layout[end+1, :] = LAxis(scene)
        end

        layout[0, :] = LText(scene, text="Super Title", textsize=50)
        layout[end+1, :] = LText(scene, text="Sub Title", textsize=50)
        layout[2:end-1, 0] = LText(scene, text="Left Text", textsize=50,
        rotation=pi/2)
        layout[2:end-1, end+1] = LText(scene, text="Right Text", textsize=50,
        rotation=-pi/2)

        scene
    end

    @cell "Column and row sizes" [grid] begin
        scene = Scene(resolution = (1200, 900), camera=campixel!)

        layout = GridLayout(
        scene, 6, 6,
        colsizes = [Fixed(200), Relative(0.25), Auto(), Auto(), Auto(2), Auto()],
        rowsizes = [Auto(), Fixed(100), Relative(0.25), Aspect(2, 1), Auto(), Auto()],
        alignmode = Outside(30, 30, 30, 30))

        for i in 2:6, j in 1:5
            if i == 6 && j == 3
                layout[i, j] = LText(scene, text="My Size is Inferred")
            else
                layout[i, j] = LRect(scene)
            end
        end

        for (j, label) in enumerate(["Fixed(200)", "Relative(0.25)", "Auto()", "Auto()", "Auto(2)"])
            layout[1, j] = LText(scene, width = Auto(false), text = label)
        end

        for (i, label) in enumerate(["Fixed(100)", "Relative(0.25)", "Aspect(2, 1)", "Auto()", "Auto()"])
            layout[i + 1, 6] = LText(scene, height = Auto(false), text = label)
        end

        scene
    end

    @cell "Trimming grids" [grid] begin

        scene, layout = layoutscene(resolution = (600, 600))

        record(scene, @replace_with_a_path(mp4), framerate=1) do io

            ax1 = layout[1, 1] = LAxis(scene, title = "Axis 1")
            recordframe!(io)

            ax2 = layout[1, 2] = LAxis(scene, title = "Axis 2")
            recordframe!(io)

            layout[2, 1] = ax2
            recordframe!(io)

            trim!(layout)
            recordframe!(io)

            layout[2, 3:4] = ax1
            recordframe!(io)

            trim!(layout)
            recordframe!(io)
        end
    end

    @cell "3D scenes" ["3d"] begin
        scene, layout = layoutscene()
        makescene() = LScene(scene, camera = cam3d!, raw = false,
            scenekw = (backgroundcolor = RGBf0(0.9, 0.9, 0.9), clear = true))
        layout[1, 1]   = makescene()
        layout[1, 2]   = makescene()
        layout[2, 1:2] = makescene()
        layout[1:2, 3] = makescene()
        foreach(LScene, layout) do s
            scatter!(s, rand(100, 3));
        end
        scene
    end

    @cell "Window resizing" [record] begin
        using MakieLayout.Animations

        container_scene = Scene(camera = campixel!, resolution = (1200, 1200))

        t = Node(0.0)

        a_width = Animation([1.0, 7.0], [1200.0, 800.0], sineio(n=2, yoyo=true, postwait=0.5))
        a_height = Animation([2.5, 8.5], [1200.0, 800.0], sineio(n=2, yoyo=true, postwait=0.5))
        scene_area = lift(t) do t
            IRect(0, 0, a_width(t), a_height(t))
        end

        scene = Scene(container_scene, scene_area, camera = campixel!)

        rect = poly!(scene, scene_area,
            raw=true, color=RGBf0(0.97, 0.97, 0.97), strokecolor=:transparent, strokewidth=0)[end]

        outer_gl = GridLayout(scene, alignmode = Outside(30))

        inner_gl = outer_gl[1, 1] = GridLayout()

        ax1 = inner_gl[1, 1] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))
        ax2 = inner_gl[1, 2] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))
        ax3 = inner_gl[2, 1:2] = LAxis(scene, xautolimitmargin=(0, 0), yautolimitmargin=(0, 0))

        guigl = inner_gl[3, 1:2] = GridLayout()
        b1 = guigl[1, 1] = LButton(scene, label = "prev", width = Auto())
        sl = guigl[1, 2] = LSlider(scene, startvalue = 6, height = 40)
        b2 = guigl[1, 3] = LButton(scene, label = "next", width = Auto())

        data = randn(200, 200) .+ 3 .* sin.((1:200) ./ 20) .* sin.((1:200)' ./ 20)
        h1 = heatmap!(ax1, data)
        h2 = heatmap!(ax2, data, colormap = :blues)
        h3 = heatmap!(ax3, data, colormap = :heat)

        agl1 = gridnest!(inner_gl, 1, 1)
        agl1[1, 2] = LColorbar(scene, h1, width = 30, label = "normal bar")
        agl1[2, 1:2] = LSlider(scene, height = 20, startvalue = 4)
        agl1[3, 1:2] = LSlider(scene, height = 20, startvalue = 5)
        agl1[4, 1:2] = LSlider(scene, height = 20, startvalue = 6)

        agl2 = gridnest!(inner_gl, 1, 2)
        agl2[1, 2] = LColorbar(scene, h2, width = 30, height = Relative(0.66), label = "two thirds bar")
        agl2gl = agl2[2, :] = GridLayout()
        agl2gl[1, 1] = LButton(scene, label = "Run", height = Auto())
        agl2gl[1, 2] = LButton(scene, label = "Start")

        agl3 = gridnest!(inner_gl, 2, 1:2)
        agl3[:, 3] = LColorbar(scene, h3, width = 30, height=200, label = "fixed height bar")
        rowsize!(agl3, 1, Auto(false, 1.0))

        inner_gl[0, :] = LText(scene, text = "MakieLayout", textsize = 50)

        record(container_scene, @replace_with_a_path(mp4), 0:1/30:9; framerate=30) do ti
        t[] = ti
        end

    end

    @cell "Aspect ratios stretching circles" [layout] begin
        using Animations

        # scene setup for animation

        container_scene = Scene(camera = campixel!, resolution = (1200, 1200))

        t = Node(0.0)

        a_width = Animation([1, 7], [1200.0, 800], sineio(n=2, yoyo=true, postwait=0.5))
        a_height = Animation([2.5, 8.5], [1200.0, 800], sineio(n=2, yoyo=true, postwait=0.5))

        scene_area = lift(t) do t
        IRect(0, 0, round(Int, a_width(t)), round(Int, a_height(t)))
        end

        scene = Scene(container_scene, scene_area, camera = campixel!)

        rect = poly!(scene, scene_area,
        raw=true, color=RGBf0(0.97, 0.97, 0.97), strokecolor=:transparent, strokewidth=0)[end]

        outer_layout = GridLayout(scene, alignmode = Outside(30))

        # example begins here

        layout = outer_layout[1, 1] = GridLayout()

        titles = ["aspect via layout", "axis aspect", "no aspect", "data aspect"]
        axs = layout[1:2, 1:2] = [LAxis(scene, title = t) for t in titles]

        # plot to all axes using an array comprehension
        [lines!(a, Circle(Point2f0(0, 0), 100f0)) for a in axs]

        rowsize!(layout, 1, Fixed(400))
        # force the layout cell [1, 1] to be square
        colsize!(layout, 1, Aspect(1, 1))

        axs[2].aspect = 1
        axs[4].autolimitaspect = 1

        rects = layout[1:2, 1:2] = [LRect(scene, color = (:black, 0.05),
        strokecolor = :transparent) for _ in 1:4]

        record(container_scene, @replace_with_a_path(mp4), 0:1/60:9; framerate=60) do ti
            t[] = ti
        end

    end

end


text:
module Interface

export simulate!
export RunTask, visualize, simulate, test, runTask, run, runModels

import Modia3D

include("run.jl")

end