text:
"""
    @abstractmethod

Macro used in generic functions that must be overloaded by derived types.
"""
macro abstractmethod(message="This function belongs to an interface definition and cannot be used.")
  quote
    error($(esc(message)))
  end
end

"""
    @notimplemented
    @notimplemented "Error message"

Macro used to raise an error, when something is not implemented.
"""
macro notimplemented(message="This function is not yet implemented")
  quote
    error($(esc(message)))
  end
end

"""
    @notimplementedif condition
    @notimplementedif condition "Error message"

Macro used to raise an error if the `condition` is true
"""
macro notimplementedif(condition,message="This function is not yet implemented")
  quote
    if $(esc(condition))
      @notimplemented $(esc(message))
    end
  end
end

"""
    @unreachable
    @unreachable "Error message"

Macro used to make sure that a line of code is never reached.
"""
macro unreachable(message="This line of code cannot be reached")
  quote
    error($(esc(message)))
  end
end

"""
  @check condition
  @check condition "Error message"

Macro used to make sure that condition is fulfilled, like `@assert`
but the check gets deactivated when running Julia with --boundscheck=no
"""
macro check(test,msg="A check failed")
  quote
    @boundscheck @assert $(esc(test)) $(esc(msg))
  end
end

struct Table{T} <: AbstractVector{Vector{T}}
  data::Vector{T}
  ptrs::Vector{Int32}
end

get_data(a::Table) = a.data
get_ptrs(a::Table) = a.ptrs

Base.size(a::Table) = (length(a.ptrs)-1,)
Base.IndexStyle(::Type{<:Table}) = IndexLinear()
function Base.getindex(a::Table{T},i::Integer) where T
  pini = a.ptrs[i]
  pend = a.ptrs[i+1]-1
  ps = pini:pend
  v = zeros(T,length(ps))
  for (i,p) in enumerate(ps)
    v[i] = a.data[p]
  end
  v
end
Base.copy(a::Table) = Table(copy(a.data),copy(a.ptrs))

function Table(a::AbstractArray{<:AbstractArray})
  data, ptrs = generate_data_and_ptrs(a)
  Table(data,ptrs)
end

function empty_table(a::Table)
  ptrs = similar(a.ptrs,eltype(a.ptrs),0+1)
  data = similar(a.data,eltype(a.data),0)
  Table(data,ptrs)
end

function generate_data_and_ptrs(vv::AbstractArray{<:AbstractArray{T}}) where T
  ptrs = Vector{Int32}(undef,length(vv)+1)
  _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = Vector{T}(undef,ndata)
  _generate_data_and_ptrs_fill_data!(data,vv)
  (data, ptrs)
end

function _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  k = 1
  for i in eachindex(vv)
    v = vv[i]
    ptrs[k+1] = length(v)
    k += 1
  end
end

function _generate_data_and_ptrs_fill_data!(data,vv)
  k = 1
  for i in eachindex(vv)
    v = vv[i]
    for vi in v
      data[k] = vi
      k += 1
    end
  end
end

function length_to_ptrs!(ptrs::AbstractArray{<:Integer})
  ptrs[1] = 1
  @inbounds for i in 1:(length(ptrs)-1)
    ptrs[i+1] += ptrs[i]
  end
end

function counts_to_ptrs(counts)
  n = length(counts)
  ptrs = Vector{Int32}(undef,n+1)
  @inbounds for i in 1:n
    ptrs[i+1] = counts[i]
  end
  length_to_ptrs!(ptrs)
  ptrs
end

function ptrs_to_counts(ptrs)
  counts = similar(ptrs,eltype(ptrs),length(ptrs)-1)
  @inbounds for i in 1:length(counts)
    counts[i] = ptrs[i+1]-ptrs[i]
  end
  counts
end

function rewind_ptrs!(ptrs::AbstractVector{<:Integer})
  @inbounds for i in (length(ptrs)-1):-1:1
    ptrs[i+1] = ptrs[i]
  end
  ptrs[1] = 1
end


text:
#!/usr/bin/env julia

import Celeste.ParallelRun: BoundingBox, infer_box
import Celeste.Log

const usage_info =
"""
Usage:
  infer-box.jl <ramin> <ramax> <decmin> <decmax>
"""

const stagedir = ENV["CELESTE_STAGE_DIR"]

if length(ARGS) != 4
    println(usage_info)
else
    box = BoundingBox(ARGS...)
    # Output gets written to the top level of the staging directory.
    # We may want to modify that in the future by changing the third argument.
    tic()
    infer_box(box, stagedir, stagedir)
    Log.info("infer_box for $box took $(toq()) seconds")
end

ext: jl


text:
using BinaryBuilder.Auditor
using BinaryBuilder.Auditor: compatible_marchs, valid_library_path

# Tests for our auditing infrastructure

@testset "Auditor - cppfilt" begin
    # We take some known x86_64-linux-gnu symbols and pass them through c++filt
    mangled_symbol_names = [
        "_ZNKSt7__cxx1110_List_baseIiSaIiEE13_M_node_countEv",
        "_ZNKSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEE6lengthEv@@GLIBCXX_3.4.21",
        "_Z10my_listlenNSt7__cxx114listIiSaIiEEE",
        "_ZNKSt7__cxx114listIiSaIiEE4sizeEv",
    ]
    unmangled_symbol_names = Auditor.cppfilt(mangled_symbol_names, Platform("x86_64", "linux"))
    @test all(unmangled_symbol_names .== [
        "std::__cxx11::_List_base<int, std::allocator<int> >::_M_node_count() const",
        "std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >::length() const@@GLIBCXX_3.4.21",
        "my_listlen(std::__cxx11::list<int, std::allocator<int> >)",
        "std::__cxx11::list<int, std::allocator<int> >::size() const",
    ])
end

@testset "Auditor - ISA tests" begin
    @test compatible_marchs(Platform("x86_64", "linux")) == ["x86_64"]
    @test compatible_marchs(Platform("x86_64", "linux"; march="x86_64")) == ["x86_64"]
    @test compatible_marchs(Platform("x86_64", "linux"; march="avx")) == ["x86_64", "avx"]
    @test compatible_marchs(Platform("x86_64", "linux"; march="avx2")) == ["x86_64", "avx", "avx2"]
    @test compatible_marchs(Platform("x86_64", "linux"; march="avx512")) == ["x86_64", "avx", "avx2", "avx512"]
    @test compatible_marchs(Platform("armv7l", "linux")) == ["armv7l"]
    @test compatible_marchs(Platform("i686", "linux"; march="prescott")) == ["pentium4", "prescott"]
    @test compatible_marchs(Platform("aarch64", "linux"; march="armv8_1")) == ["armv8_0", "armv8_1"]

    product = ExecutableProduct("main", :main)

    # The microarchitecture of the product doesn't match the target architecture: complain!
    mktempdir() do build_path
        platform = Platform("x86_64", "linux"; march="avx")
        build_output_meta = nothing
        @test_logs (:info, "Building for x86_64-linux-gnu-march+avx") (:warn, r"is avx512, not avx as desired.$") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "isa_tests",
                v"1.0.0",
                [DirectorySource(build_tests_dir)],
                # Build the test suite, install the binaries into our prefix's `bin`
                raw"""
                cd ${WORKSPACE}/srcdir/isa_tests
                make -j${nproc} CFLAGS="-march=skylake-avx512 -mtune=skylake-avx512" install
                install_license /usr/include/ltdl.h
                """,
                # Build for our platform
                [platform],
                # Ensure our executable products are built
                [product],
                # No dependencies
                Dependency[];
                preferred_gcc_version=v"6",
                lock_microarchitecture=false,
            )

            # Extract our platform's build
            @test haskey(build_output_meta, platform)
            tarball_path, tarball_hash = build_output_meta[platform][1:2]
            @test isfile(tarball_path)

            # Unpack it somewhere else
            @test verify(tarball_path, tarball_hash)
            testdir = joinpath(build_path, "testdir")
            mkdir(testdir)
            unpack(tarball_path, testdir)
            prefix = Prefix(testdir)

            # Run ISA test
            readmeta(locate(product, prefix)) do oh
                detected_isa = Auditor.analyze_instruction_set(oh, platform; verbose=true)
                @test detected_isa == "avx512"
            end
        end
    end

    # The instruction set of the product is compatible with the target
    # architecture, but it's lower than desired: issue a gentle warning
    mktempdir() do build_path
        platform = Platform("x86_64", "linux"; march="avx512")
        build_output_meta = nothing
        @test_logs (:info, "Building for x86_64-linux-gnu-march+avx512") (:warn, r"is avx, not avx512 as desired. You may be missing some optimization flags during compilation.$") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "isa_tests",
                v"1.0.0",
                [DirectorySource(build_tests_dir)],
                # Build the test suite, install the binaries into our prefix's `bin`
                raw"""
                cd ${WORKSPACE}/srcdir/isa_tests
                make -j${nproc} CFLAGS="-march=sandybridge -mtune=sandybridge" install
                install_license /usr/include/ltdl.h
                """,
                # Build for our platform
                [platform],
                # Ensure our executable products are built
                [product],
                # No dependencies
                Dependency[];
                preferred_gcc_version=v"6",
                lock_microarchitecture=false,
            )

            # Extract our platform's build
            @test haskey(build_output_meta, platform)
            tarball_path, tarball_hash = build_output_meta[platform][1:2]
            @test isfile(tarball_path)

            # Unpack it somewhere else
            @test verify(tarball_path, tarball_hash)
            testdir = joinpath(build_path, "testdir")
            mkdir(testdir)
            unpack(tarball_path, testdir)
            prefix = Prefix(testdir)

            # Run ISA test
            readmeta(locate(product, prefix)) do oh
                detected_march = Auditor.analyze_instruction_set(oh, platform; verbose=true)
                @test detected_march == "avx"
            end
        end
    end

    # The microarchitecture of the product matches the target architecture: no warnings!
    for march in ("x86_64", "avx", "avx2", "avx512")
        mktempdir() do build_path
            platform = Platform("x86_64", "linux"; march=march)
            build_output_meta = nothing
            @test_logs (:info, "Building for x86_64-linux-gnu-march+$(march)") match_mode=:any begin
                build_output_meta = autobuild(
                    build_path,
                    "isa_tests",
                    v"1.0.0",
                    [DirectorySource(build_tests_dir)],
                    # Build the test suite, install the binaries into our prefix's `bin`
                    raw"""
                    cd ${WORKSPACE}/srcdir/isa_tests
                    make -j${nproc} install
                    install_license /usr/include/ltdl.h
                    """,
                    # Build for our platform
                    [platform],
                    # Ensure our executable products are built
                    [product],
                    # No dependencies
                    Dependency[];
                    # Use a recent version of GCC to make sure we can detect the
                    # ISA accurately even with new optimizations
                    preferred_gcc_version=v"8"
                )
            end

            # Extract our platform's build
            @test haskey(build_output_meta, platform)
            tarball_path, tarball_hash = build_output_meta[platform][1:2]
            @test isfile(tarball_path)

            # Unpack it somewhere else
            @test verify(tarball_path, tarball_hash)
            testdir = joinpath(build_path, "testdir")
            mkdir(testdir)
            unpack(tarball_path, testdir)
            prefix = Prefix(testdir)

            # Run ISA test
            readmeta(locate(product, prefix)) do oh
                detected_march = Auditor.analyze_instruction_set(oh, platform; verbose=true)
                if march == "avx2"
                    # Detecting the ISA isn't 100% reliable and it's even less
                    # accurate when looking for AVX2 features
                    @test_broken march == detected_march
                else
                    @test march == detected_march
                end
            end
        end
    end
end

@testset "Auditor - cxxabi selection" begin
    libcxxstringabi_test = LibraryProduct("libcxxstringabi_test", :libcxxstringabi_test)

    # Factor the autobuild() out
    function do_build(build_path, script, platform, gcc_version)
        autobuild(
            build_path,
            "libcxxstringabi_test",
            v"1.0.0",
            # Copy in the build_tests sources
            [DirectorySource(build_tests_dir)],
            script,
            # Build for this platform
            [platform],
            # The products we expect to be build
            [libcxxstringabi_test],
            # No depenedencies
            Dependency[];
            preferred_gcc_version=gcc_version
        )
    end

    for platform in (Platform("x86_64", "linux"; cxxstring_abi="cxx03"),
                     Platform("x86_64", "linux"; cxxstring_abi="cxx11"))
        # Look across multiple gcc versions; there can be tricksy interactions here
        for gcc_version in (v"4", v"6", v"9")
            # Do each build within a separate temporary directory
            mktempdir() do build_path
                script = raw"""
                    cd ${WORKSPACE}/srcdir/cxxstringabi_tests
                    make install
                    install_license /usr/share/licenses/libuv/LICENSE
                """
                build_output_meta = do_build(build_path, script, platform, gcc_version)
                # Extract our platform's build
                @test haskey(build_output_meta, platform)
                tarball_path, tarball_hash = build_output_meta[platform][1:2]
                @test isfile(tarball_path)

                # Unpack it somewhere else
                @test verify(tarball_path, tarball_hash)
                testdir = joinpath(build_path, "testdir")
                mkdir(testdir)
                unpack(tarball_path, testdir)
                prefix = Prefix(testdir)

                # Ensure that the library detects as the correct cxxstring_abi:
                readmeta(locate(libcxxstringabi_test, prefix)) do oh
                    detected_cxxstring_abi = Auditor.detect_cxxstring_abi(oh, platform)
                    @test detected_cxxstring_abi == cxxstring_abi(platform)
                end

                # Explicitly test cxx string abi mismatches
                if gcc_version > v"4"
                    script = """
                        mkdir -p \${libdir}
                        /opt/\${target}/bin/\${target}-g++ -fPIC \\
                            -D_GLIBCXX_USE_CXX11_ABI=$(cxxstring_abi(platform) == "cxx03" ? "1" : "0") \\
                            -o \${libdir}/libcxxstringabi_test.\${dlext} \\
                            -shared \${WORKSPACE}/srcdir/cxxstringabi_tests/lib.cc
                        install_license /usr/share/licenses/libuv/LICENSE
                    """
                    @test_logs (:warn, r"ignoring our choice of compiler") match_mode=:any begin
                        do_build(build_path, script, platform, gcc_version)
                    end
                end
            end
        end
    end

    # Explicitly test not setting a cxx string abi at all
    script = raw"""
        cd ${WORKSPACE}/srcdir/cxxstringabi_tests
        make install
        install_license /usr/share/licenses/libuv/LICENSE
    """
    platform = Platform("x86_64", "linux")
    mktempdir() do build_path
        @test_logs (:warn, r"contains std::string values") match_mode=:any begin
            do_build(build_path, script, platform, v"6")
        end
    end
end

@testset "Auditor - .dll moving" begin
    for platform in [Platform("x86_64", "windows")]
        mktempdir() do build_path
            build_output_meta = nothing
            @test_logs (:warn, r"lib/libfoo.dll should be in `bin`") (:warn, r"Simple buildsystem detected") match_mode=:any begin
                build_output_meta = autobuild(
                    build_path,
                    "dll_moving",
                    v"1.0.0",
                    GitSource[],
                    # Install a .dll into lib
                    raw"""
                    mkdir -p ${prefix}/lib
                    cc -o ${prefix}/lib/libfoo.${dlext} -shared /usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
                    install_license /usr/include/ltdl.h
                    """,
                    # Build for our platform
                    [platform],
                    # Ensure our library product is built
                    [LibraryProduct("libfoo", :libfoo)],
                    # No dependencies
                    Dependency[]
                )
            end

            @test haskey(build_output_meta, platform)
            tarball_path, tarball_hash = build_output_meta[platform][1:2]
            @test isfile(tarball_path)

            # Test that `libfoo.dll` gets moved to `bin` if it's a windows
            contents = list_tarball_files(tarball_path)
            @test "bin/libfoo.$(platform_dlext(platform))" in contents
        end
    end
end

@testset "Auditor - .dylib identity mismatch" begin
    mktempdir() do build_path
        no_id = LibraryProduct("no_id", :no_id)
        abs_id = LibraryProduct("abs_id", :wrong_id)
        wrong_id = LibraryProduct("wrong_id", :wrong_id)
        right_id = LibraryProduct("right_id", :wrong_id)
        platform = Platform("x86_64", "macos")

        build_output_meta = autobuild(
            build_path,
            "dll_moving",
            v"1.0.0",
            FileSource[],
            # Intsall a .dll into lib
            raw"""
            mkdir -p "${libdir}"
            SRC=/usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
            cc -o ${libdir}/no_id.${dlext} -shared $SRC
            cc -o ${libdir}/abs_id.${dlext} -Wl,-install_name,${libdir}/abs_id.${dlext} -shared $SRC
            cc -o ${libdir}/wrong_id.${dlext} -Wl,-install_name,@rpath/totally_different.${dlext} -shared $SRC
            cc -o ${libdir}/right_id.${dlext} -Wl,-install_name,@rpath/right_id.${dlext} -shared $SRC
            install_license /usr/include/ltdl.h
            """,
            # Build for MacOS
            [platform],
            # Ensure our executable products are built
            Product[no_id, abs_id, wrong_id, right_id],
            # No dependencies
            Dependency[],
        )

        # Extract our platform's build
        @test haskey(build_output_meta, platform)
        tarball_path, tarball_hash = build_output_meta[platform][1:2]
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        prefix = Prefix(testdir)

        # Helper to extract the dylib id of a path
        function get_dylib_id(path)
            return readmeta(path) do oh
                dylib_id_lcs = [lc for lc in MachOLoadCmds(oh) if isa(lc, MachOIdDylibCmd)]
                @test !isempty(dylib_id_lcs)
                return dylib_name(first(dylib_id_lcs))
            end
        end

        # Locate the build products within the prefix, ensure that all the dylib ID's
        # now match the pattern `@rpath/$(basename(p))`
        no_id_path = locate(no_id, prefix; platform=platform)
        abs_id_path = locate(abs_id, prefix; platform=platform)
        right_id_path = locate(right_id, prefix; platform=platform)
        for p in (no_id_path, abs_id_path, right_id_path)
            @test any(startswith.(p, libdirs(prefix)))
            @test get_dylib_id(p) == "@rpath/$(basename(p))"
        end

        # Only if it already has an `@rpath/`-ified ID, it doesn't get touched.
        wrong_id_path = locate(wrong_id, prefix; platform=platform)
        @test any(startswith.(wrong_id_path, libdirs(prefix)))
        @test get_dylib_id(wrong_id_path) == "@rpath/totally_different.dylib"

        # Ensure that this bianry is codesigned
        @test BinaryBuilder.Auditor.check_codesigned(right_id_path, platform)
    end
end

@testset "Auditor - absolute paths" begin
    mktempdir() do build_path
        sharedir = joinpath(realpath(build_path), "share")
        mkpath(sharedir)
        open(joinpath(sharedir, "foo.conf"), "w") do f
            write(f, "share_dir = \"$sharedir\"")
        end

        # Test that `audit()` warns about an absolute path within the prefix
        @test_logs (:warn, r"share/foo.conf contains an absolute path") match_mode=:any begin
            Auditor.audit(Prefix(build_path); verbose=true)
        end
    end
end

@testset "Auditor - broken symlinks" begin
    mktempdir() do build_path
        bindir = joinpath(realpath(build_path), "bin")
        mkpath(bindir)
        # Test both broken and working (but external) symlinks
        symlink("../../artifacts/1a2b3/lib/libzmq.dll.a", joinpath(bindir, "libzmq.dll.a"))
        # The following symlinks shouldn't raise a warning
        symlink("/bin/bash", joinpath(bindir, "bash.exe"))
        symlink("libfoo.so.1.2.3", joinpath(bindir, "libfoo.so"))

        # Test that `audit()` warns about broken symlinks
        @test_logs (:warn, r"Broken symlink: bin/libzmq.dll.a") match_mode=:any begin
            Auditor.warn_deadlinks(build_path)
        end
    end
end

@testset "Auditor - gcc version" begin
    # These tests assume our gcc version is concrete (e.g. that Julia is linked against libgfortran)
    our_libgfortran_version = libgfortran_version(platform)
    @test our_libgfortran_version != nothing

    mktempdir() do build_path
        hello_world = ExecutableProduct("hello_world_fortran", :hello_world_fortran)
        build_output_meta = @test_logs (:warn, r"CompilerSupportLibraries_jll") match_mode=:any begin
            autobuild(
                build_path,
                "hello_fortran",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build the test suite, install the binaries into our prefix's `bin`
                raw"""
                    # Build fortran hello world
                    make -j${nproc} -sC /usr/share/testsuite/fortran/hello_world install
                    # Install fake license just to silence the warning
                    install_license /usr/share/licenses/libuv/LICENSE
                    """,
                # Build for our platform
                [platform],
                #
                Product[hello_world],
                # No dependencies
                Dependency[];
            )
        end

        # Extract our platform's build, run the hello_world tests:
        output_meta = select_platform(build_output_meta, platform)
        @test output_meta != nothing
        tarball_path, tarball_hash = output_meta[1:2]

        # Ensure the build products were created
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)

        # Attempt to run the executable, we expect it to work since it's our platform:
        hello_world_path = locate(hello_world, Prefix(testdir); platform=platform)
        with_libgfortran() do
            @test readchomp(`$hello_world_path`) == "Hello, World!"
        end

        # If we audit the testdir, pretending that we're trying to build an ABI-agnostic
        # tarball, make sure it warns us about it.
        @test_logs (:warn, r"links to libgfortran!") match_mode=:any begin
            @test !Auditor.audit(Prefix(testdir); platform=BinaryBuilderBase.abi_agnostic(platform), autofix=false)
            # Make sure audit is otherwise happy with the executable
            @test Auditor.audit(Prefix(testdir); platform=platform, autofix=false)
        end

        # Let's pretend that we're building for a different libgfortran version:
        # audit should warn us.
        libgfortran_versions = (3, 4, 5)
        other_libgfortran_version = libgfortran_versions[findfirst(v -> v != our_libgfortran_version.major, libgfortran_versions)]
        @test_logs (:warn, Regex("but we are supposedly building for libgfortran$(other_libgfortran_version)")) readmeta(hello_world_path) do oh
            p = deepcopy(platform)
            p["libgfortran_version"] = "$(other_libgfortran_version).0.0"
            @test !Auditor.audit(Prefix(testdir); platform=p, autofix=false)
        end
    end
end

@testset "Auditor - soname matching" begin
    mktempdir() do build_path
        build_output_meta = nothing
        linux_platform = Platform("x86_64", "linux")
        @test_logs (:info, r"creating link to libfoo\.so\.1\.0\.0") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "soname_matching",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build the library only with the versioned name
                raw"""
                mkdir -p "${libdir}"
                cc -o "${libdir}/libfoo.${dlext}.1.0.0" -fPIC -shared /usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
                # Set the soname to a non-existing file
                patchelf --set-soname libfoo.so "${libdir}/libfoo.${dlext}.1.0.0"
                """,
                # Build for Linux
                [linux_platform],
                # Ensure our library product is built
                [LibraryProduct("libfoo", :libfoo)],
                # No dependencies
                Dependency[];
                autofix = true,
                verbose = true,
                require_license = false
            )
        end
        # Extract our platform's build
        @test haskey(build_output_meta, linux_platform)
        tarball_path, tarball_hash = build_output_meta[linux_platform][1:2]
        # Ensure the build products were created
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        @test readlink(joinpath(testdir, "lib", "libfoo.so")) == "libfoo.so.1.0.0"
    end
end

@testset "Auditor - execution permission" begin
    mktempdir() do build_path
        build_output_meta = nothing
        product = LibraryProduct("libfoo", :libfoo)
        @test_logs (:info, r"Making .*libfoo.* executable") match_mode=:any begin
            build_output_meta = autobuild(
                build_path,
                "exec",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build a library without execution permissions
                raw"""
                mkdir -p "${libdir}"
                cc -o "${libdir}/libfoo.${dlext}" -fPIC -shared /usr/share/testsuite/c/dyn_link/libfoo/libfoo.c
                chmod 640 "${libdir}/libfoo.${dlext}"
                """,
                # Build for our platform
                [platform],
                # Ensure our library product is built
                [product],
                # No dependencies
                Dependency[];
                verbose = true,
                require_license = false
            )
        end

        # Extract our platform's build
        @test haskey(build_output_meta, platform)
        tarball_path, tarball_hash = build_output_meta[platform][1:2]
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        libfoo_path = joinpath(testdir, build_output_meta[platform][4][product]["path"])
        # Tar.jl normalizes permissions of executable files to 0o755, instead of
        # recording exact original permissions:
        # https://github.com/JuliaIO/Tar.jl/blob/37766a22f5a6ac9f07022d83debd5db7d7a4b896/README.md#permissions
        @test_broken filemode(libfoo_path) & 0o777 == 0o750
    end
end

@testset "Auditor - other checks" begin
    platform = Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc")
    mktempdir() do build_path
        build_output_meta = @test_logs (:error, r"libsoft.so does not match the hard-float ABI") match_mode=:any begin
            autobuild(
                build_path,
                "hard_float_ABI",
                v"1.0.0",
                # No sources
                FileSource[],
                # Build a library which doesn't link to the standard library and
                # forces the soft-float ABI
                raw"""
                mkdir -p "${libdir}" "${bindir}"
                # This library has hard-float ABI
                echo 'int test() { return 0; }' | cc -shared -fPIC -o "${libdir}/libhard.${dlext}" -x c -
                # This library has soft-float ABI
                echo 'int _start() { return 0; }' | /opt/${target}/bin/${target}-gcc -nostdlib -shared -mfloat-abi=soft -o "${libdir}/libsoft.${dlext}" -x c -
                # hello_world built by Go doesn't specify any float ABI
                make -C /usr/share/testsuite/go/hello_world/
                cp "/tmp/testsuite/${target}/go/hello_world/hello_world" "${bindir}/hello_world"
                """,
                # Build for Linux armv7l hard-float
                [platform],
                # Ensure our library product is built
                [
                    LibraryProduct("libhard", :libhard),
                    LibraryProduct("libsoft", :libsoft),
                    ExecutableProduct("hello_world", :hello_world),
                ],
                # No dependencies
                Dependency[];
                compilers = [:c, :go],
                verbose = true,
                require_license = false
            )
        end

        @test haskey(build_output_meta, platform)
        tarball_path, tarball_hash = build_output_meta[platform][1:2]
        @test isfile(tarball_path)

        # Unpack it somewhere else
        @test verify(tarball_path, tarball_hash)
        testdir = joinpath(build_path, "testdir")
        mkdir(testdir)
        unpack(tarball_path, testdir)
        # Remove libsoft.so, we want to run audit only on the other products
        rm(joinpath(testdir, "lib", "libsoft.so"))
        # Make sure `hello_world` passes the float ABI check even if it doesn't
        # set `EF_ARM_ABI_FLOAT_HARD`.
        @test Auditor.audit(Prefix(testdir); platform=platform, require_license=false)
    end
end

@testset "valid_library_path" begin
    linux = Platform("x86_64", "linux")
    macos = Platform("x86_64", "macos")
    windows = Platform("x86_64", "windows")

    @test valid_library_path("/usr/libc.dylib", macos)
    @test !valid_library_path("/usr/libc.dylib.", macos)
    @test !valid_library_path("/usr/libc.dylib.1", macos)
    @test !valid_library_path("/usr/libc.dylib", linux)
    @test !valid_library_path("/usr/libc.dylib", windows)

    @test valid_library_path("libc.dll", windows)
    @test !valid_library_path("libc.dll.1", windows)
    @test !valid_library_path("libc.dll", linux)
    @test !valid_library_path("libc.dll", macos)

    @test valid_library_path("/usr/libc.so", linux)
    @test valid_library_path("/usr/libc.so.1", linux)
    @test valid_library_path("/usr/libc.so.1.2", linux)
    @test !valid_library_path("/usr/libc.so.", linux)
    @test !valid_library_path("/usr/libc.sot", linux)
    @test !valid_library_path("/usr/libc.so", macos)
    @test !valid_library_path("/usr/libc.so", windows)
end


text:
include("file1.jl")
include("file2.jl")


text:
using ..SCIP: libscip

const uint8_t = UInt8
const LLONG_MAX = typemax(Clonglong)
const LLONG_MIN = typemin(Clonglong)
const DBL_MAX = typemax(Cfloat)
const PRIx64 = "llx"
const SIZE_MAX = typemax(Csize_t)
const UINT64_C = UInt64
const UINT32_C = UInt32


text:
function InterpolatingMicrosphere(arg0::jint, arg1::jint, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::UnitSphereRandomVectorGenerator)
    return InterpolatingMicrosphere((jint, jint, jdouble, jdouble, jdouble, UnitSphereRandomVectorGenerator), arg0, arg1, arg2, arg3, arg4, arg5)
end

function copy(obj::InterpolatingMicrosphere)
    return jcall(obj, "copy", InterpolatingMicrosphere, ())
end

function get_dimension(obj::InterpolatingMicrosphere)
    return jcall(obj, "getDimension", jint, ())
end

function get_size(obj::InterpolatingMicrosphere)
    return jcall(obj, "getSize", jint, ())
end

function value(obj::InterpolatingMicrosphere, arg0::Vector{jdouble}, arg1::Vector{Vector{jdouble}}, arg2::Vector{jdouble}, arg3::jdouble, arg4::jdouble)
    return jcall(obj, "value", jdouble, (Vector{jdouble}, Vector{Vector{jdouble}}, Vector{jdouble}, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end


text:
# operation_type: Forward Modeling, RTM, FWI, TWRI
function H5Modeling(;
  # models
  h5vel,
  h5density=nothing,
  h5qualityFactor=nothing,
  nb::Real=40,  # number of absorbing boundaries points on each side
  h5epsilon=nothing,
  h5delta=nothing,
  h5tetha=nothing,
  h5phi=nothing,
  model_xkey::String="CDP_X",
  model_ykey::String="CDP_Y",

  # geometry
  h5geom=nothing,
  geom_src_pkey::String="SP",
  geom_src_xkey::String="SRCX",
  geom_src_ykey::String="SRCY",
  geom_src_zkey::String="SES",
  geom_rec_xkey::String="GRPX",
  geom_rec_ykey::String="GRPY",
  geom_rec_zkey::String="RGE",

  # judi
  opt::JUDI.JUDIOptions,

  # common settings
  computation_type::String,
  # LSRTM
  lsrtm_niter::Int=nothing,
  lsrtm_batchsize::Int=nothing,
  # FWI
  fwi_niter::Int=nothing,
  fwi_batchsize::Int=nothing,
  fwi_vmin::Number,
  fwi_vmax::Number,
  # TWRI
  twri_opt::TWRIOptions=nothing
  )

  # pylogging = pyimport("logging")
  h5gt = pyimport("h5gtpy._h5gt")
  h5geo = pyimport("h5geopy._h5geo")

  @info "H5Modeling started\n"

  if isnothing(h5vel)
    @error "No velocity model\n"
    return
  end

  if isnothing(opt)
    @error "No JUDI Options for modeling\n"
    return
  end

  if !isnothing(opt.dt_comp) && opt.dt_comp <= 0
    @error "dt_comp can't be less or equal to zero\n"
    return
  end

  # prepare physical parameters
  phpvel, model_orientation = H5ReadPhysicalParameter(
    h5vel, phptype=VELOCITY, xkey=model_xkey, ykey=model_ykey)
  if isnothing(phpvel)
    @error "Unable to read PhysicalParameter: VELOCITY\n"
    return
  end

  # declare PHP vars to be able use thm as arguments (or unknown var error will appear)
  phpdensity = nothing
  phpqualityFactor = nothing
  phpepsilon = nothing
  phpdelta = nothing
  phptheta = nothing
  phpphi = nothing

  if !isnothing(h5density)
    phpdensity, model_orientation = H5ReadPhysicalParameter(
      h5density, phptype=DENSITY, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phpdensity)
      @error "Unable to read PhysicalParameter: DENSITY\n"
      return
    end
  end

  if !isnothing(h5qualityFactor)
    phpqualityFactor, model_orientation = H5ReadPhysicalParameter(
      h5qualityFactor, phptype=QUALITYFACTOR, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phpqualityFactor)
      @error "Unable to read PhysicalParameter: QUALITYFACTOR\n"
      return
    end
  end

  if !isnothing(h5epsilon)
    phpepsilon, model_orientation = H5ReadPhysicalParameter(
      h5epsilon, phptype=EPSILON, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phpepsilon)
      @error "Unable to read PhysicalParameter: EPSILON\n"
      return
    end
  end

  if !isnothing(h5delta)
    phpdelta, model_orientation = H5ReadPhysicalParameter(
      h5delta, phptype=DELTA, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phpdelta)
      @error "Unable to read PhysicalParameter: DELTA\n"
      return
    end
  end

  if !isnothing(h5tetha)
    phptheta, model_orientation = H5ReadPhysicalParameter(
      h5tetha, phptype=THETA, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phptheta)
      @error "Unable to read PhysicalParameter: THETA\n"
      return
    end
  end

  if !isnothing(h5phi)
    phpphi, model_orientation = H5ReadPhysicalParameter(
      h5phi, phptype=PHI, xkey=model_xkey, ykey=model_ykey)
    if isnothing(phpphi)
      @error "Unable to read PhysicalParameter: PHI\n"
      return
    end
  end

  model = JUDI.Model(phpvel.n, phpvel.d, phpvel.o, phpvel.data,
          rho=isnothing(phpdensity) ? nothing : phpdensity.data,
          qp=isnothing(phpqualityFactor) ? nothing : phpqualityFactor.data,
          epsilon=isnothing(phpepsilon) ? nothing : phpepsilon.data,
          delta=isnothing(phpdelta) ? nothing : phpdelta.data,
          theta=isnothing(phptheta) ? nothing : phptheta.data,
          phi=isnothing(phpphi) ? nothing : phpphi.data,
          nb=nb)

  # Geometry
  if isnothing(h5geom)
    @error "Geometry H5Seis is NULL\n"
    return
  end

  model_origin_x = model.o[1]
  model_origin_y = 0
  if length(model.o) == 3
    model_origin_y = model.o[2]
  end

  con = H5SeisCon(seis=h5geom, pkey=geom_src_pkey)
  recGeometry = H5GeometryOOC(
    h5geo=h5geo,
    container=con,
    key="receiver",
    xkey=geom_rec_xkey,
    ykey=geom_rec_ykey,
    zkey=geom_rec_zkey,
    do_coord_transform=true,
    model_origin_x=model_origin_x,
    model_origin_y=model_origin_y,
    model_orientation=model_orientation)
  srcGeometry = H5GeometryOOC(
    h5geo=h5geo,
    container=con,
    key="source",
    xkey=geom_src_xkey,
    ykey=geom_src_ykey,
    zkey=geom_src_zkey,
    do_coord_transform=true,
    model_origin_x=model_origin_x,
    model_origin_y=model_origin_y,
    model_orientation=model_orientation)

  # recGeometry = Geometry(
  #   container=con,
  #   xkey=geom_rec_xkey,
  #   ykey=geom_rec_ykey,
  #   zkey=geom_rec_zkey,
  #   do_coord_transform=true,
  #   model_origin_x=model_origin_x,
  #   model_origin_y=model_origin_y,
  #   model_orientation=model_orientation)
  # srcGeometry = Geometry(
  #   container=con,
  #   xkey=geom_src_xkey,
  #   ykey=geom_src_ykey,
  #   zkey=geom_src_zkey,
  #   do_coord_transform=true,
  #   model_origin_x=model_origin_x,
  #   model_origin_y=model_origin_y,
  #   model_orientation=model_orientation)

  if isnothing(srcGeometry) || isnothing(srcGeometry)
    @error "Unable to prepare Source and Receiver JUDI Geometry objects from H5Seis\n"
    return
  end

  # setup wavelet
  dt = Float32(abs(h5geom.getSampRate("ms")))
  nt = Integer(h5geom.getNSamp())
  t = Float32((nt-1)*dt)
  wavelet = ricker_wavelet(t, dt, opt.f0)
  q = judiVector(srcGeometry, wavelet)

  # JUDI (Devito) doesn't create a dir to save data
  if opt.save_data_to_disk
    mkpath(opt.file_path)
  end

  # simply to print the value
  ntComp = get_computational_nt(srcGeometry, recGeometry, model, dt=opt.dt_comp)
  dtComp = calculate_dt(model)

  # print modeling info
  @info "computation_type: $computation_type\n"
  @info "ntComp: $ntComp\n"
  @info "dtComp: $dtComp\n"
  @info "model settings:\n"
  @info "model.o: $(model.o)\n"
  @info "model.n: $(model.n)\n"
  @info "model.d: $(model.d)\n"
  @info "model_xkey: $model_xkey\n"
  @info "model_ykey: $model_ykey\n"
  @info "model_origin_x: $model_origin_x\n"
  @info "model_origin_y: $model_origin_y\n"
  @info "model_orientation: $model_orientation\n"
  @info "geometry settings:\n"
  @info "geom_src_xkey: $geom_src_xkey\n"
  @info "geom_src_ykey: $geom_src_ykey\n"
  @info "geom_src_zkey: $geom_src_zkey\n"
  @info "geom_rec_xkey: $geom_rec_xkey\n"
  @info "geom_rec_ykey: $geom_rec_ykey\n"
  @info "geom_rec_zkey: $geom_rec_zkey\n"
  @info "dt: $dt\n"
  @info "t: $t\n"
  @info "src_frq (kHz): $(opt.f0)"
  @info "nsrc: $(length(con))\n"
  @info "pkey: $(con.pkey)\n"
  @info "pkeyvals: $(con.pkeyvals)\n"
  @info "options: $opt\n"

  if computation_type == "Forward Modeling"
    H5ForwardModeling(model=model, opt=opt, q=q, recGeometry=recGeometry)
  elseif computation_type == "RTM"
    dobs = judiVector(recGeometry, con)
    mkpath(opt.file_path)
    H5RTM(model=model, opt=opt, q=q, dobs=dobs)
  elseif computation_type == "LSRTM"
    @error "LSRTM not ready yet\n"
  elseif computation_type == "FWI"
    dobs = judiVector(recGeometry, con)
    mkpath(opt.file_path)
    H5FWI(model=model, opt=opt, q=q, dobs=dobs,
          niterations=fwi_niter, batchsize=fwi_batchsize,
          vmin=fwi_vmin, vmax=fwi_vmax)
  elseif computation_type == "TWRI"
    @error "TWRI not ready yet\n"
  end

end


text:
n, p = 50, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=12341222)

@testset "Tools" begin
    lr  = LogisticRegression(1.0, 2.0; fit_intercept=false,
                             scale_penalty_with_samples = false)
    obj = objective(lr, X, y)
    J   = LogisticLoss() + L2Penalty() + 2L1Penalty()
    @test obj(θ) ≈ J(y, X*θ, θ)
    lr  = LogisticRegression(1.0, 2.0;
                             fit_intercept=true, penalize_intercept=true,
                             scale_penalty_with_samples = false)
    obj = objective(lr, X, y)
    J   = LogisticLoss() + L2Penalty() + 2L1Penalty()
    @test obj(θ1) ≈ J(y, X1*θ1, θ1)
end

@testset "Objectives" begin
    λ = 0.5
    γ = 0.3
    δ = 0.7

    # L2 / L2
    r = RidgeRegression(λ; fit_intercept=false, scale_penalty_with_samples = false)
    J = objective(r, X, y)
    @test J(θ) ≈ sum(abs2.(X*θ - y))/2 + λ * sum(abs2.(θ))/2
    r = RidgeRegression(λ; fit_intercept=true, penalize_intercept=true,
                        scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ sum(abs2.(X1*θ1 - y1))/2 + λ * sum(abs2.(θ1))/2
    r = RidgeRegression(λ; fit_intercept=true, penalize_intercept=false,
                           scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ sum(abs2.(X1*θ1 - y1))/2 + λ * sum(abs2.(θ1 .* vcat(ones(p), 0)))/2

    # Scale with samples
    r = RidgeRegression(λ; fit_intercept=false, scale_penalty_with_samples = true)
    J = objective(r, X, y)
    @test J(θ) ≈ sum(abs2.(X*θ - y))/2 + λ * length(y) * sum(abs2.(θ))/2

    # L2 / L2+L1
    r = ElasticNetRegression(λ, γ; fit_intercept=true, penalize_intercept=true,
                                   scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ sum(abs2.(X1*θ1 - y1))/2 +
                    λ * sum(abs2.(θ1))/2 +
                    γ * sum(abs.(θ1))
    r = ElasticNetRegression(λ, γ; fit_intercept=true, penalize_intercept=false,
                                   scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ sum(abs2.(X1*θ1 - y1))/2 +
                    λ * sum(abs2.(θ1 .* vcat(ones(p), 0)))/2 +
                    γ * sum(abs.(θ1 .* vcat(ones(p), 0)))

    # Logistic / L2 + L1
    r = LogisticRegression(λ, γ; fit_intercept=true, penalize_intercept=true,
                                 scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ LogisticLoss()(y1, X1*θ1) +
                    λ * sum(abs2.(θ1))/2 +
                    γ * sum(abs.(θ1))
    r = LogisticRegression(λ, γ; fit_intercept=true, penalize_intercept=false,
                                 scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ LogisticLoss()(y1, X1*θ1) +
                    λ * sum(abs2.(θ1 .* vcat(ones(p), 0)))/2 +
                    γ * sum(abs.(θ1 .* vcat(ones(p), 0)))

    # Robust / L2
    r = HuberRegression(δ, λ, γ; fit_intercept=true, penalize_intercept=true,
                                 scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ RobustLoss(HuberRho(δ))(y1, X1*θ1) +
                    λ * sum(abs2.(θ1))/2 +
                    γ * sum(abs.(θ1))
    r = HuberRegression(δ, λ, γ; fit_intercept=true, penalize_intercept=false,
                                 scale_penalty_with_samples = false)
    J = objective(r, X, y1)
    @test J(θ1) ≈ RobustLoss(HuberRho(δ))(y1, X1*θ1) +
                    λ * sum(abs2.(θ1 .* vcat(ones(p), 0)))/2 +
                    γ * sum(abs.(θ1 .* vcat(ones(p), 0)))
end


text:
"""
The [centroid of a triangle](https://en.wikipedia.org/wiki/Centroid#Of_a_triangle) is the point of
intersection of its medians (the lines joining each vertex with the midpoint of the opposite side).
The [Centroids Exists Theorem](https://sites.math.rutgers.edu/~zeilberg/PG/CentroidExists.html)
states that the three medians indeed intersect at single point.
"""
module Centroid

using Plots
using SymPy
using PlaneGeometry
using PlaneGeometry.GeoPlots

export
    centroid_proof, centroid_draw, centroid_rand, centroid

"""
    centroid_draw(📐::Triangle)

Verify Centroid Exists Theorem for the triangle `📐`.
"""
function centroid_draw(tri::Triangle)
    trishape=shape(tri)

    cpt, medians = centroid(tri)

    if cpt == nothing
        return (nothing, false)
    end

    plt = plot(tri, fill=(0, :green), aspect_ratio=:equal, fillalpha= 0.2, leg=false)

    for e in medians
        plot!(e, leg=false, color=:orange)
    end

    plt = scatter!(shape([vertices(tri)..., cpt]),
            leg=false,
            aspect_ratio=:equal,
            color=[:red, :red, :red, :green],
            series_annotations = text.(["A", "B", "C", "cpt"], :bottom))

    plt, true
end

"""
    centroid_draw(A, B, C)

Verify Centroid Exists theorem.
"""
function centroid_draw(A, B, C)
    tri = Triangle(A, B, C)

    centroid_draw(tri)
end

"""
    centroid_draw(xA, yA, xB, yB, xC, yC)

Verify Centroid Exists theorem.
"""
function centroid_draw(xA, yA, xB, yB, xC, yC)
    A = Point(xA, yA);
    B = Point(xB, yB);
    C = Point(xC, yC);

    centroid_draw(A, B, C)
end

"""
    centroid_rand()

Verify Centroid Exists Theorem for a random triangle.
"""
function centroid_rand()
    pts = rand(0:1//10:1,6);
    plt, hold = centroid_draw(pts...)
end

"""
    centroid_proof()

Prove Centroid Exists Theorem.
"""
function centroid_proof()
    @vars by cx positive=true;
    @vars cy;

    A = Point(0, 0); B = Point(0, by); C = Point(cx, cy);
    tri = Triangle(A, B, C)

    cpt, medians = centroid(tri)
    cpt != nothing
end

"""
    centroid(📐)

Find the centroid and the three medians of `📐`.
"""
function centroid(tri)
    medians = map(i->median(circshift(vertices(tri), i)...), 0:2)
    cpt = concurrent(medians)
    cpt, medians
end

end


text:
"""
Executes a elastic flag benchmark proposed by Turek et al in "Proposal for
numerical benchmarking of fluid-structure interaction between an elastic
object and laminar incompressible flow"

The default parameters correspond to FSI2 test case.
FS1 test can be run modifying the default values of the following parameters:
Um=0.2,
Re=20,
rho_s=1.0e-3
FS3 test can be run modifying the default values of the following parameters:
Um=2.0,
Re=200,
rho_s=1.0e-3
E_s=5.6e6

It supports the following mesh motion strategies:
 - Laplacian
 - Linear elasticity
 - Noehookean elasticity (not tested)
 - Biharmonic (Laplacian squared)

The fluid and structure can be strongly coupled, i.e. same variational space for
fluid and solid displacements and velocities, or weakly coupled with different
variational spaces and coupling enforced through Nitche's method.
"""
function execute(problem::FSIProblem{:elasticFlag}; kwargs...)

  # Problem setting (Default FSI-2)
  println("Setting Elastic flag problem parameters")
  Um = _get_kwarg(:Um,kwargs,1.0)
  H = _get_kwarg(:H,kwargs,0.41)
  ⌀ = _get_kwarg(:D,kwargs,0.1)
  # Solid properties
  E_s = _get_kwarg(:E_s,kwargs,1.4e6)
  ν_s = _get_kwarg(:nu_s,kwargs,0.4)
  ρ_s = _get_kwarg(:rho_s,kwargs,1.0e4)
  # Fluid properties
  ρ_f = _get_kwarg(:rho_f,kwargs,1.0e3)
  Re = _get_kwarg(:Re,kwargs, 100.0)
  μ_f = ρ_f * Um * ⌀ / Re
  γ_f = _get_kwarg(:gamma_f,kwargs,1.0)
  # Mesh properties
  E_m = _get_kwarg(:E_m,kwargs,1.0)
  ν_m = _get_kwarg(:nu_m,kwargs,-0.1)
  α_m = _get_kwarg(:alpha_m,kwargs,1.0e-5)
  weight_strategy = _get_kwarg(:alpha_m_weight,kwargs,"constant")
  # Time stepping
  t0 = _get_kwarg(:t0,kwargs,0.0)
  tf = _get_kwarg(:tf,kwargs,0.1)
  dt = _get_kwarg(:dt,kwargs,0.1)
  θ  = _get_kwarg(:theta,kwargs,0.5)
  # Post-process
  is_vtk = _get_kwarg(:is_vtk,kwargs,false)

  # Mesh strategy
  strategyName = _get_kwarg(:strategy,kwargs,"laplacian")
  strategy = WeakForms.MeshStrategy{Symbol(strategyName)}()

  # Fluid-Solid coupling type
  couplingName = _get_kwarg(:coupling,kwargs,"strong")
  coupling = WeakForms.Coupling{Symbol(couplingName)}()

  # Define BC functions
  println("Defining Boundary conditions")
  u_in(x, t) = VectorValue(1.5 * Um * x[2] * (H - x[2]) / ((H / 2)^2), 0.0)
  u_noSlip(x, t) = VectorValue(0.0, 0.0)
  u_in(t::Real) = x -> u_in(x, t)
  u_noSlip(t::Real) = x -> u_noSlip(x, t)
  ∂tu_in(t) = x -> VectorValue(0.0, 0.0)
  ∂tu_in(x, t) = ∂tu_in(t)(x)
  T_in = typeof(u_in)
  T_noSlip = typeof(u_noSlip)
  @eval ∂t(::$T_in) = $∂tu_in
  @eval ∂t(::$T_noSlip) = $∂tu_in
  bconds = get_boundary_conditions(problem,strategy,coupling,u_in,u_noSlip)

  # Forcing terms
  f(t) = x -> VectorValue(0.0,0.0)

  # Discrete model
  println("Defining discrete model")
  modelName = _get_kwarg(:model,kwargs,"../models/elasticFlag.json")
  model = DiscreteModelFromFile(modelName)
  model_solid = DiscreteModel(model,tags="solid")
  model_fluid = DiscreteModel(model,tags="fluid")
  models = Dict(:Ω => model, :Ωf => model_fluid, :Ωs => model_solid)

  # Triangulations
  println("Defining triangulations")
  Tₕ = get_FSI_triangulations(models,coupling)

  # Quadratures
  println("Defining quadratures")
  order = _get_kwarg(:order,kwargs,2)
  dTₕ = get_FSI_measures(Tₕ,order)

  # Test FE Spaces
  println("Defining FE spaces")
  Y_ST, X_ST, Y_FSI, X_FSI = get_FE_spaces(strategy,coupling,models,order,bconds)

  # Stokes problem for initial solution
  println("Defining Stokes operator")
  op_ST = get_Stokes_operator(X_ST,Y_ST,strategy,dTₕ[:Ωf],μ_f,f(0.0))

  # Setup equation parameters
  mesh_params = Dict{Symbol,Any}(
    :w_strategy=>weight_strategy,
    :α=>α_m,
    :E=>E_m,
    :ν=>ν_m,
  )
  fluid_params = Dict{Symbol,Any}(
    :μ=>μ_f,
    :ρ=>ρ_f,
    :fu=>f,
    :fv=>f,
  )
  solid_params = Dict{Symbol,Any}(
    :ρ=>ρ_s,
    :E=>E_s,
    :ν=>ν_s,
    :fu=>f,
    :fv=>f,
  )
  Γi_params = Dict{Symbol,Any}(
    :μ=>μ_f,
    :γ=>γ_f,
    :dt=>dt
  )
  params = (mesh_params, fluid_params, solid_params, Γi_params)

  # FSI problem
  println("Defining FSI operator")
  op_FSI = get_FSI_operator(X_FSI,Y_FSI,coupling,strategy,Tₕ,dTₕ,params)

  # Setup output files
  folderName = "fsi-results"
  fileName = "fields"
  if !isdir(folderName)
    mkdir(folderName)
  end
  filePath = join([folderName, fileName], "/")

  # Solve Stokes problem
  @timeit "ST problem" begin
  println("Defining Stokes solver")
  xh = solve(op_ST)
  if(is_vtk)
    writePVD(filePath, Tₕ[:Ωf], [(xh, 0.0)])
  end
end

# Solve FSI problem
@timeit "FSI problem" begin
println("Defining FSI solver")
xh0  = interpolate(xh,X_FSI(0.0))
nls = NLSolver(
  show_trace = true,
  method = :newton,
  linesearch = BackTracking(),
  ftol = 1.0e-6,
  iterations = 50
  )
odes =  ThetaMethod(nls, dt, θ)
solver = TransientFESolver(odes)
xht = solve(solver, op_FSI, xh0, t0, tf)
end

# Compute outputs
out_params = Dict{Symbol,Any}(
  :μ=>μ_f,
  :Um=>Um,
  :⌀=>⌀,
  :ρ=>ρ_f,
  :θ=>θ,
  :bdegree=>2*order,
  :filePath=>filePath,
  :is_vtk=>is_vtk,
  )
output = computeOutputs(xh0,xht,coupling,strategy,models,Tₕ,dTₕ,out_params)

end

function get_boundary_conditions(
  problem::FSIProblem{:elasticFlag},
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:strong},
  u_in,
  u_noSlip
  )
  boundary_conditions = (
    # Tags
    FSI_Vu_tags = ["inlet", "noSlip", "cylinder","fixed","outlet"],
    FSI_Vv_tags = ["inlet", "noSlip", "cylinder","fixed"],
    ST_Vu_tags = ["inlet", "noSlip", "cylinder","interface","outlet"],
    ST_Vv_tags = ["inlet", "noSlip", "cylinder","interface"],
    # Values,
    FSI_Vu_values = [u_noSlip, u_noSlip, u_noSlip, u_noSlip, u_noSlip],
    FSI_Vv_values = [u_in, u_noSlip, u_noSlip, u_noSlip],
    ST_Vu_values = [u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)],
    ST_Vv_values = [u_in(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)],
  )
end

function get_boundary_conditions(
  problem::FSIProblem{:elasticFlag},
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:weak},
  u_in,
  u_noSlip
  )
  boundary_conditions = (
    # Tags
    FSI_Vw_f_tags = ["inlet", "noSlip", "cylinder","outlet"],
    FSI_Vu_f_tags = ["inlet", "noSlip", "cylinder","outlet"],
    FSI_Vv_f_tags = ["inlet", "noSlip", "cylinder"],
    FSI_Vu_s_tags = ["fixed"],
    FSI_Vv_s_tags = ["fixed"],
    ST_Vu_tags = ["inlet", "noSlip", "cylinder","interface","outlet"],
    ST_Vv_tags = ["inlet", "noSlip", "cylinder","interface"],
    # Values,
    FSI_Vw_f_values = [u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)],
    FSI_Vu_f_values = [u_noSlip, u_noSlip, u_noSlip, u_noSlip],
    FSI_Vv_f_values = [u_in, u_noSlip, u_noSlip],
    FSI_Vu_s_values = [u_noSlip],
    FSI_Vv_s_values = [u_noSlip],
    ST_Vu_values = [u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)],
    ST_Vv_values = [u_in(0.0), u_noSlip(0.0), u_noSlip(0.0), u_noSlip(0.0)],
  )
end

function computeOutputs(xh0,xht,coupling::WeakForms.Coupling,strategy::WeakForms.MeshStrategy,models,Tₕ,quads,params)

  # Unpack parameters
  bdegree = params[:bdegree]
  μ = params[:μ]
  Um = params[:Um]
  ⌀ = params[:⌀]
  ρ = params[:ρ]
  θ = params[:θ]
  filePath = params[:filePath]
  is_vtk = params[:is_vtk]
  if( typeof(strategy) == WeakForms.MeshStrategy{:biharmonic} )
    uvpindex = [2,3,4]
  else
    uvpindex = [1,2,3]
  end

  ## Surface triangulation
  Γc = BoundaryTriangulation(models[:Ω], tags="cylinder")
  dΓc = Measure(Γc, bdegree)
  n_Γc = get_normal_vector(Γc)

  # Aux function
  traction_boundary(n,u,v,p) = n ⋅ WeakForms.Pᵥ_Ωf(μ,u,v)  + WeakForms.Pₚ_Ωf(u,p) * n
  function traction_closure(coupling,n,u,v,p)
    if typeof(coupling) == WeakForms.Coupling{:weak}
      traction_boundary(n,u,v,p).⁺
    else
      traction_boundary(n,u,v,p)
    end
  end
  traction_interface(n,u,v,p) = traction_closure(coupling,n,u,v,p)

  # Get interface normal Vector
  n_Γi = get_normal_vector(Tₕ[:Γi])

  ## Initialize arrays
  tpl = Real[]
  FDpl = Real[]
  FLpl = Real[]
  uhn = xh0[uvpindex[1]]
  vhn = xh0[uvpindex[2]]
  phn = xh0[uvpindex[3]]

  ## Loop over steps
  outfiles = paraview_collection(filePath, append=true
  ) do pvd
  for (i,(xh, t)) in enumerate(xht)
    println("STEP: $i, TIME: $t")
    println("============================")

    ## Get the solution at n+θ (where velocity and pressure are balanced)
    uh = xh[uvpindex[1]]
    vh = xh[uvpindex[2]]
    ph = xh[uvpindex[3]]
    xθ(xₙ₊₁,xₙ) = θ*xₙ₊₁ + (1-θ)*xₙ

    # Integrate on the cylinder
    FDc, FLc = ∑(∫(traction_boundary(n_Γc,xθ(uh,uhn),xθ(vh,vhn),xθ(ph,phn)) )dΓc )

    # Integrate on the interface
    FDi, FLi = ∑(∫(traction_boundary(n_Γc,xθ(uh,uhn),xθ(vh,vhn),xθ(ph,phn)) )dΓc )

    FD = FDc + FDi
    FL = FLc + FLi

    ## Drag and lift Forces
    push!(tpl, t)
    push!(FDpl, -FD)
    push!(FLpl, -FL)

    ## store step n
    uhn = uh
    vhn = vh
    phn = ph

    # Write to PVD
    if(is_vtk)
      uh = xh[uvpindex[1]]
      vh = xh[uvpindex[2]]
      ph = xh[uvpindex[3]]
      pvd[t] = createvtk(
      Tₕ[:Ω],
      filePath * "_$t.vtu",
      cellfields = ["uh" => uh, "vh" => vh, "ph" => ph]
      )
    end
  end
end

return (tpl, FDpl, FLpl)
end


text:
# to run the script, you should have these packages installed
using Plots
using StatsPlots
using BenchmarkTools
using Statistics
using Indicators
using RollingFunctions
using ImageFiltering
using MaxMinFilters

y = randn(1000000)
windows = [5, 15, 50, 150]
t_mf1 = zeros(size(windows))
t_mf2 = zeros(size(windows))
t_img = zeros(size(windows))
t_rol = zeros(size(windows))
t_ind = zeros(size(windows))

for i in eachindex(windows)
    w = windows[i]

    t = @benchmark MaxMinFilters.movmax($y, $w) samples=3
    t_mf1[i] = median(t).time / 10^6

    t = @benchmark MaxMinFilters.movmaxmin($y, $w) samples=3
    t_mf2[i] = median(t).time / 10^6

    t = @benchmark ImageFiltering.mapwindow($extrema, $y, $w) samples=3
    t_img[i] = median(t).time / 10^6

    t = @benchmark RollingFunctions.runmax($y, $w) samples=3
    t_rol[i] = median(t).time / 10^6

    t = @benchmark Indicators.runmax($y, n = $w, cumulative = $false) samples=3
    t_ind[i] = median(t).time / 10^6
end

ctg = repeat(["1 - MaxMinFilters . movmax",
              "2 - MaxMinFilters . movmaxmin",
              "3 - ImageFiltering . mapwindow extrema",
              "4 - RollingFunctions . runmax",
              "5 - Indicators . runmax",
              ], inner = length(windows))

#bars = log10.(hcat(t_mf1, t_mf2, t_img, t_rol, t_ind)) # log scale
bars = hcat(t_mf1, t_mf2, t_img, t_rol, t_ind)
xnames = repeat(collect(1:length(windows)), outer = size(bars,2))

groupedbar(xnames, bars, group = ctg, xlabel = "Window lengths", ylabel = "Time, ms",
        title = "Moving maximum performance, array length = 1000000.", bar_width = 0.7,
        lw = 0, framestyle = :box, legend=:topleft, ylims = (0, 530))

xticks!(collect(1:4), string.(windows))
#yticks!(collect(0:1:3), "10^{" .* string.(collect(0:1:3)) .* "}") # log scale

#png("plot.png")


text:
# Write a macro that runs a block of code backward.
macro backward(b)
  quote
    map(eval,reverse($b.args))[1]
  end
end

max_stars_repo_name: hemslo/seven-in-seven

max_issues_repo_name: hemslo/seven-in-seven

max_forks_repo_name: hemslo/seven-in-seven


text:
"""
 Holds the tableau of a stochastic implicit Runge-Kutta method.
 qdrift holds the RK coefficients for the drift part,
 and qdiff holds the RK coefficients for the diffusion part of the SDE.
"""
struct TableauSIRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift and qdiff are understood as the classical orders of these methods.

    function TableauSIRK{T}(name, s, qdrift, qdiff) where {T}
        # THE COMMENTED OUT PART WAS FOR TableauFIRK. MAY IMPLEMENT SOMETHING
        # SIMILAR FOR TableauSIRK LATER.

        # if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
        #     warn("Initializing TableauFIRK with explicit tableau ", q.name, ".\n",
        #          "You might want to use TableauERK instead.")
        # elseif q.s > 1 && istril(q.a)
        #     warn("Initializing TableauFIRK with diagonally implicit tableau ", q.name, ".\n",
        #          "You might want to use TableauDIRK instead.")
        # end

        @assert s == qdrift.s == qdiff.s

        new(name, s, qdrift, qdiff)
    end
end

function TableauSIRK(name::Symbol, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}) where {T}
    TableauSIRK{T}(name, qdrift.s, qdrift, qdiff)
end

function TableauSIRK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T}, order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T}) where {T}
    TableauSIRK{T}(name, length(c_drift), CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff))
end

# TODO function readTableauSIRKFromFile(dir::AbstractString, name::AbstractString)

"""
Parameters for right-hand side function of implicit Runge-Kutta methods.
   A - if positive, the upper bound of the Wiener process increments; if A=0.0, no truncation
"""
mutable struct ParametersSIRK{DT, TT, ET <: SDE{DT,TT}, D, M, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauSIRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    A::DT

    t::TT
    q::Vector{DT}
end

function ParametersSIRK(equ::ET, tab::TableauSIRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}, A::DT) where {DT, TT, ET <: SDE{DT,TT}}
    @assert equ.m == length(ΔW) == length(ΔZ)
    ParametersSIRK{DT, TT, ET, equ.d, equ.m, tab.s}(equ, tab, Δt, ΔW, ΔZ, A, 0, zeros(DT, equ.d))
end

"""
Structure for holding the internal stages Q, the values of the drift vector
and the diffusion matrix evaluated at the internal stages V=v(Q), B=B(Q),
and the increments Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW
"""
struct NonlinearFunctionCacheSIRK{DT}
    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    B::Vector{Matrix{DT}}
    Y::Vector{Vector{DT}}

    v::Vector{DT}
    b::Matrix{DT}
    y::Vector{DT}

    function NonlinearFunctionCacheSIRK{DT}(d, m, s) where {DT}

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, d, s)
        V = create_internal_stage_vector(DT, d, s)
        B = create_internal_stage_vector(DT, d, m, s)
        Y = create_internal_stage_vector(DT, d, s)

        # create velocity and update vector
        v = zeros(DT,d)
        b = zeros(DT,d,m)
        y = zeros(DT,d)

        new(Q, V, B, Y, v, b, y)
    end
end

"""
Unpacks the data stored in x = (Y[1][1], Y[1][2], ... Y[1][D], Y[2][1], ...)
into Y::Vector{Vector}, calculates the internal stages Q, the values of the RHS
of the SDE ( v(Q) and B(Q) ), and assigns them to V and B.
Unlike for FIRK, here Y = Δt a v(Q) + ̃a B(Q) ΔW
"""
function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, B::Vector{Matrix{ST}}, Y::Vector{Vector{ST}},
                         params::ParametersSIRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    local tᵢ::TT

    @assert S == length(Q) == length(V) == length(B)

    # copy x to Y and calculate Q
    for i in eachindex(Q)
        @assert D == size(B[i],1) == length(Q[i]) == length(V[i])
        @assert M == size(B[i],2)
        for k in eachindex(Q[i])
            Y[i][k] = x[D*(i-1)+k]
        end
        Q[i] .= params.q .+ Y[i]
    end

    # compute VQ = v(Q) and BQ=B(Q)
    for i in eachindex(Q,V,B)
        tᵢ = params.t + params.Δt * params.tab.qdrift.c[i]
        # calculates v(t,Q[i]) and assigns to the i-th column of V
        params.equ.v(tᵢ, Q[i], V[i])
        # calculates B(t,Q[i]) and assigns to the matrix B[i]
        params.equ.B(tᵢ, Q[i], B[i])
    end
end

"Compute stages of stochastic implicit Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersSIRK{DT,TT,ET,D,M,S}) where {ST,DT,TT,ET,D,M,S}

    cache = NonlinearFunctionCacheSIRK{ST}(D, M, S)

    quote
        compute_stages!(x, $cache.Q, $cache.V, $cache.B, $cache.Y, params)

        local y1::ST
        local y2::ST

        # compute b = - (Y-AV)
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                for j in 1:S
                    y1 += params.tab.qdrift.a[i,j] * $cache.V[j][k] * params.Δt + params.tab.qdiff.a[i,j] * dot($cache.B[j][k,:], params.ΔW)
                    y2 += params.tab.qdrift.â[i,j] * $cache.V[j][k] * params.Δt + params.tab.qdiff.â[i,j] * dot($cache.B[j][k,:], params.ΔW)
                end
                b[D*(i-1)+k] = - $cache.Y[i][k] + (y1 + y2)
            end
        end
    end
end

"Stochastic implicit Runge-Kutta integrator."
struct IntegratorSIRK{DT, TT, PT <: ParametersSIRK{DT,TT},
                              ST <: NonlinearSolver{DT}, N} <: StochasticIntegrator{DT,TT}
    params::PT
    solver::ST
    # InitialGuessSDE not implemented for SIRK
    #iguess::IT
    fcache::NonlinearFunctionCacheSIRK{DT}

    q::Matrix{Vector{TwicePrecision{DT}}}
end

# K - the integer in the bound A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; K=0 no truncation
function IntegratorSIRK(equation::SDE{DT,TT,VT,BT,N}, tableau::TableauSIRK{TT}, Δt::TT; K::Int=0) where {DT,TT,VT,BT,N}
    D = equation.d
    M = equation.m
    NS= equation.ns
    NI= equation.n
    S = tableau.s

    # create params
    K==0 ? A = 0.0 : A = sqrt( 2*K*Δt*abs(log(Δt)) )
    params = ParametersSIRK(equation, tableau, Δt, zeros(DT,M), zeros(DT,M), A)

    # create solver
    solver = create_nonlinear_solver(DT, D*S, params)

    # Not implementing InitialGuessSDE
    # create initial guess
    #iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache for internal stage vectors and update vectors
    fcache = NonlinearFunctionCacheSIRK{DT}(D, M, S)

    # create solution vectors
    q = create_solution_vector(DT, D, NS, NI)

    # create integrator
    IntegratorSIRK{DT, TT, typeof(params), typeof(solver), N}(params, solver, fcache, q)
end

equation(integrator::IntegratorSIRK) = integrator.params.equ
timestep(integrator::IntegratorSIRK) = integrator.params.Δt
tableau(integrator::IntegratorSIRK) = integrator.params.tab
Base.eltype(integrator::IntegratorSIRK{DT, TT, PT, ST, N}) where {DT, TT, PT, ST, N} = DT

function initialize!(int::IntegratorSIRK{DT,TT}, sol::SolutionSDE, k::Int, m::Int) where {DT,TT}
    check_solution_dimension_asserts(sol, k, m)

    # copy the m-th initial condition for the k-th sample path
    get_initial_conditions!(sol, int.q[k,m], k, m)

    # Not implementing InitialGuessSDE
    # # initialise initial guess
    # initialize!(int.iguess, m, sol.t[0], int.q[m])
end

"""
This function computes initial guesses for Y and assigns them to int.solver.x
The prediction is calculated using an explicit integrator.
"""
function initial_guess!(int::IntegratorSIRK{DT,TT}) where {DT,TT}

    # NOT IMPLEMENTING InitialGuessSDE

    # SIMPLE SOLUTION
    # The simplest initial guess for Y is 0
    # int.solver.x .= zeros(eltype(int), int.params.tab.s*ndims(int))

    # USING AN EXPLICIT INTEGRATOR TO COMPUTE AN INITIAL GUESS
    # Below we use the R2 method of Burrage & Burrage to calculate
    # the internal stages at the times c[1]...c[s].
    # This approach seems to give very good approximations if the time step
    # and magnitude of noise are not too large. If the noise intensity is too big,
    # one may have to perform a few iterations of the explicit method with a smaller
    # time step, use a higher-order explicit method (e.g. CL or G5), or use
    # the simple solution above.

    local tV1 = zeros(DT,int.params.equ.d)
    local tV2 = zeros(DT,int.params.equ.d)
    local tB1 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local tB2 = zeros(DT,int.params.equ.d, int.params.equ.m)
    local Q   = zeros(DT,int.params.equ.d)
    local t2::TT
    local Δt_local::TT
    local ΔW_local = zeros(DT,int.params.equ.m)

    # When calling this function, int.params should contain the data:
    # int.params.q - the solution at the previous time step
    # int.params.t - the time of the previous step
    # int.params.ΔW- the increment of the Brownian motion for the current step

    #Evaluating the functions v and B at t,q - same for all stages
    int.params.equ.v(int.params.t, int.params.q, tV1)
    int.params.equ.B(int.params.t, int.params.q, tB1)

    for i in 1:int.params.tab.s

        Δt_local  = int.params.tab.qdrift.c[i]*int.params.Δt
        ΔW_local .= int.params.tab.qdrift.c[i]*int.params.ΔW

        Q = int.params.q + 2. / 3. * Δt_local * tV1 + 2. / 3. * tB1 * ΔW_local

        t2 = int.params.t + 2. / 3. *Δt_local

        int.params.equ.v(t2, Q, tV2)
        int.params.equ.B(t2, Q, tB2)

        #Calculating the Y's and assigning them to the array int.solver.x as initial guesses
        for j in 1:int.params.equ.d
            int.solver.x[(i-1)*int.params.equ.d+j] =  Δt_local*(1. / 4. * tV1[j] + 3. / 4. * tV2[j]) + dot( (1. / 4. * tB1[j,:] + 3. / 4. * tB2[j,:]), ΔW_local )
        end
    end

end

"""
Integrate SDE with a stochastic implicit Runge-Kutta integrator.
  Integrating the k-th sample path for the m-th initial condition
"""
function integrate_step!(int::IntegratorSIRK{DT,TT}, sol::SolutionSDE{DT,TT,NQ,NW}, k::Int, m::Int, n::Int) where {DT,TT,NQ,NW}
    check_solution_dimension_asserts(sol, k, m, n)

    # set time for nonlinear solver
    int.params.t  = sol.t[0] + (n-1)*int.params.Δt
    int.params.q .= int.q[k, m]

    # copy the increments of the Brownian Process
    if NW==1
        #1D Brownian motion, 1 sample path
        int.params.ΔW[1] = sol.W.ΔW[n-1]
        int.params.ΔZ[1] = sol.W.ΔZ[n-1]
    elseif NW==2
        #Multidimensional Brownian motion, 1 sample path
        int.params.ΔW .= sol.W.ΔW[n-1]
        int.params.ΔZ .= sol.W.ΔZ[n-1]
    elseif NW==3
        #1D or Multidimensional Brownian motion, k-th sample path
        int.params.ΔW .= sol.W.ΔW[n-1,k]
        int.params.ΔZ .= sol.W.ΔZ[n-1,k]
    end

    # truncate the increments ΔW with A
    if int.params.A>0
        for i in eachindex(int.params.ΔW)
            if int.params.ΔW[i]<-int.params.A
                int.params.ΔW[i] = -int.params.A
            elseif int.params.ΔW[i]>int.params.A
                int.params.ΔW[i] = int.params.A
            end
        end
    end

    # compute initial guess and assign to int.solver.x
    initial_guess!(int)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, int.fcache.Q, int.fcache.V, int.fcache.B, int.fcache.Y, int.params)

    # compute final update
    update_solution!(int.q[k,m], int.fcache.V, int.fcache.B, int.params.tab.qdrift.b, int.params.tab.qdrift.b̂, int.params.tab.qdiff.b, int.params.tab.qdiff.b̂, int.params.Δt, int.params.ΔW)

    # # NOT IMPLEMENTING InitialGuessSDE
    # # # copy solution to initial guess
    # # update!(int.iguess, m, sol.t[0] + n*int.params.Δt, int.q[m])

    # take care of periodic solutions
    cut_periodic_solution!(int.q[k,m], int.params.equ.periodicity)

    # # copy to solution
    set_solution!(sol, int.q[k,m], n, k, m)
end


text:
module ViscousFlow

#using DocStringExtensions
using Reexport
using UnPack
@reexport using ImmersedLayers
@reexport using GridUtilities

export ViscousIncompressibleFlowProblem
export setup_grid, viscousflow_system, setup_problem, surface_point_spacing

#= Supporting functions =#

setup_problem(g;kwargs...) =
    ViscousIncompressibleFlowProblem(g;timestep_func=DEFAULT_TIMESTEP_FUNC,
                                       kwargs...)

setup_problem(g,bl;bc=nothing,kwargs...) =
    ViscousIncompressibleFlowProblem(g,bl;timestep_func=DEFAULT_TIMESTEP_FUNC,
                                       bc=get_bc_func(bc),
                                       kwargs...)

function grid_spacing(phys_params)
    gridRe = get(phys_params,"grid Re",DEFAULT_GRID_RE)
    Δx = gridRe/get_Reynolds_number(phys_params)
end

function get_Reynolds_number(phys_params)
  haskey(phys_params,"Re") || error("No Reynolds number set")
  return phys_params["Re"]
end

function default_timestep(sys)
    @unpack phys_params, motions, forcing = sys
    g = get_grid(sys)
    Fo = get(phys_params,"Fourier",DEFAULT_FOURIER_NUMBER)
    Co = get(phys_params,"CFL",DEFAULT_CFL_NUMBER)
    Re = get_Reynolds_number(phys_params)

    Uscale = 1.0
    if !isnothing(motions)
      Uscale, _ = maxlistvelocity(sys)
    end

    Δt = min(Fo*Re*cellsize(g)^2,Co*cellsize(g)/Uscale)
    return Δt
end

function default_freestream(t,phys_params)
    Vinfmag = get(phys_params,"freestream speed",0.0)
    Vinf_angle = get(phys_params,"freestream angle",0.0)
    Uinf = Vinfmag*cos(Vinf_angle)
    Vinf = Vinfmag*sin(Vinf_angle)
    return Uinf, Vinf
end

function default_vsplus(t,base_cache,phys_params,motions)
  vsplus = zeros_surface(base_cache)
  return vsplus
end

function default_vsminus(t,base_cache,phys_params,motions)
    vsminus = zeros_surface(base_cache)
    return vsminus
end

const DEFAULT_GRID_RE = 2.0
const DEFAULT_FOURIER_NUMBER = 1.0
const DEFAULT_CFL_NUMBER = 0.5
const DEFAULT_DS_TO_DX_RATIO = 1.4
const DEFAULT_FREESTREAM_FUNC = default_freestream
const DEFAULT_TIMESTEP_FUNC = default_timestep
const DEFAULT_VSPLUS_FUNC = default_vsplus
const DEFAULT_VSMINUS_FUNC = default_vsminus

#=
Process keywords
=#

function get_freestream_func(forcing::Dict)
    return get(forcing,"freestream",DEFAULT_FREESTREAM_FUNC)
end

get_freestream_func(::Nothing) = get_freestream_func(Dict())

function get_forcing_models(forcing::Dict)
    return get(forcing,"forcing models",nothing)
end

get_forcing_models(::Nothing) = get_forcing_models(Dict())

function get_bc_func(bc_in::Dict)
    bc = Dict()
    bc["exterior"] = haskey(bc_in,"exterior") ? bc_in["exterior"] : DEFAULT_VSPLUS_FUNC
    bc["interior"] = haskey(bc_in,"interior") ? bc_in["interior"] : DEFAULT_VSMINUS_FUNC
    return bc
end

get_bc_func(::Nothing) = get_bc_func(Dict())

#=
Defining the extra cache and extending prob_cache
=#

@ilmproblem ViscousIncompressibleFlow vector

struct ViscousIncompressibleFlowCache{CDT,FRT,DVT,VFT,VORT,DILT,VELT,FCT} <: AbstractExtraILMCache
   cdcache :: CDT
   fcache :: FRT
   dvb :: DVT
   vb_tmp :: DVT
   v_tmp :: VFT
   dv :: VFT
   dv_tmp :: VFT
   w_tmp :: VORT
   divv_tmp :: DILT
   velcache :: VELT
   f :: FCT
end

function ImmersedLayers.prob_cache(prob::ViscousIncompressibleFlowProblem,
                                   base_cache::BasicILMCache{N,scaling}) where {N,scaling}
    @unpack phys_params, forcing = prob
    @unpack gdata_cache, gcurl_cache, g = base_cache

    dvb = zeros_surface(base_cache)
    vb_tmp = zeros_surface(base_cache)

    v_tmp = zeros_grid(base_cache)
    dv = zeros_grid(base_cache)
    dv_tmp = zeros_grid(base_cache)
    w_tmp = zeros_gridcurl(base_cache)
    divv_tmp = zeros_griddiv(base_cache)

    velcache = VectorFieldCache(base_cache)

    # Construct a Lapacian outfitted with the viscosity
    Re = get_Reynolds_number(phys_params)
    viscous_L = Laplacian(base_cache,gcurl_cache,1.0/Re)

    # Create cache for the convective derivative
    cdcache = ConvectiveDerivativeCache(base_cache)

    fcache = nothing

    # Create cache for the forcing regions
    fmods = get_forcing_models(forcing)
    fcache = ForcingModelAndRegion(fmods,base_cache)

    # The state here is vorticity, the constraint is the surface traction
    f = _get_ode_function_list(viscous_L,base_cache)

    ViscousIncompressibleFlowCache(cdcache,fcache,dvb,vb_tmp,v_tmp,dv,dv_tmp,w_tmp,divv_tmp,velcache,f)
end

_get_ode_function_list(viscous_L,base_cache::BasicILMCache{N}) where {N} =
                ODEFunctionList(state = zeros_gridcurl(base_cache),
                                constraint = zeros_surface(base_cache),
                                ode_rhs=viscousflow_vorticity_ode_rhs!,
                                lin_op=viscous_L,
                                bc_rhs=viscousflow_vorticity_bc_rhs!,
                                constraint_force = viscousflow_vorticity_constraint_force!,
                                bc_op = viscousflow_vorticity_bc_op!)

_get_ode_function_list(viscous_L,base_cache::BasicILMCache{0}) =
                 ODEFunctionList(state = zeros_gridcurl(base_cache),
                                 ode_rhs=viscousflow_vorticity_ode_rhs!,
                                 lin_op=viscous_L)

#= setup API =#

"""
    setup_grid(xlim::Tuple,ylim::Tuple,phys_params::Dict)

Construct a Cartesian grid with limits `xlim` and `ylim`
and spacing determined by the Reynolds number in the `phys_params`.
"""
function setup_grid(xlim::Tuple,ylim::Tuple,phys_params)
    PhysicalGrid(xlim,ylim,grid_spacing(phys_params))
end

"""
    surface_point_spacing(g::PhysicalGrid,phys_params)

Calculate the surface point spacing for a given grid, using
the specified parameter "point spacing ratio" in the physical parameters,
or the default value (1.4) if not specified.
"""
function surface_point_spacing(g::PhysicalGrid,phys_params)
    ds_to_dx = get(phys_params,"point spacing ratio",DEFAULT_DS_TO_DX_RATIO)
    return ds_to_dx*cellsize(g)
end

"""
    viscousflow_system(grid,[bodies];kwargs...)

Construct the operators and cache variables for a viscous flow problem.

The `kwargs` are the optional keyword aguments. There are several, some of
which are crucial for certain types of problems.

- `ddftype = ` to set the DDF type. The default is `CartesianGrids.Yang3`.
- `scaling = ` to set the scaling type, `GridScaling` (default) or `IndexScaling`.
- `dtype = ` to set the element type to `Float64` (default) or `ComplexF64`.
- `phys_params = ` A dictionary to pass in physical parameters
- `bc = ` A dictionary to pass in boundary condition data or functions, using "external"
          and "internal" keys to pass in functions that provide the
          corresponding surface data outside and inside the surface(s).
- `forcing = ` A dictionary to pass in forcing models (via the "forcing models" key),
                or to pass in an alternative model for
               the freestream velocity (with the "freestream" key)
- `motions = ` to provide function(s) that specify the velocity of the immersed surface(s). Note: if this keyword is used, it is assumed that surfaces will move.
- `timestep_func =` to pass in a function for time-dependent problems that provides the time-step size.
                  It is expected that this function takes in two arguments,
                  the `grid::PhysicalGrid` and `phys_params`, and returns the time step.
                  It defaults to the basic Fourier/CFL type function `default_timestep`
"""
function viscousflow_system(args...;kwargs...)
  prob = setup_problem(args...;kwargs...)
  return construct_system(prob)
end

#= ODE functions =#

function viscousflow_vorticity_ode_rhs!(dw,w,sys::ILMSystem,t)
  @unpack extra_cache, base_cache = sys
  @unpack v_tmp, dv = extra_cache

  velocity!(v_tmp,w,sys,t)
  viscousflow_velocity_ode_rhs!(dv,v_tmp,sys,t)
  curl!(dw,dv,base_cache)

  return dw
end

function viscousflow_velocity_ode_rhs!(dv,v,sys::ILMSystem,t)
    @unpack bc, forcing, phys_params, extra_cache, base_cache = sys
    @unpack dvb, dv_tmp, cdcache, fcache = extra_cache

    Re = get_Reynolds_number(phys_params)
    over_Re = 1.0/Re

    fill!(dv,0.0)

    # Calculate the convective derivative
    convective_derivative!(dv_tmp,v,base_cache,cdcache)
    dv .-= dv_tmp

    # Calculate the double-layer term
    prescribed_surface_jump!(dvb,t,sys)
    surface_divergence_symm!(dv_tmp,over_Re*dvb,base_cache)
    dv .-= dv_tmp

    # Apply forcing
    apply_forcing!(dv_tmp,v,t,fcache,phys_params)
    dv .+= dv_tmp

    return dv
end

function viscousflow_vorticity_bc_rhs!(vb,sys::ILMSystem,t)
    @unpack bc, forcing,extra_cache, base_cache, phys_params = sys
    @unpack dvb, vb_tmp, v_tmp, velcache, divv_tmp = extra_cache
    @unpack dcache, ϕtemp = velcache

    viscousflow_velocity_bc_rhs!(vb,sys,t)

    # Subtract influence of scalar potential field
    fill!(divv_tmp,0.0)
    prescribed_surface_jump!(dvb,t,sys)
    scalarpotential_from_masked_divv!(ϕtemp,divv_tmp,dvb,base_cache,dcache)
    vecfield_from_scalarpotential!(v_tmp,ϕtemp,base_cache)
    interpolate!(vb_tmp,v_tmp,base_cache)
    vb .-= vb_tmp

    # Subtract influence of free stream
    freestream_func = get_freestream_func(forcing)
    Uinf, Vinf = freestream_func(t,phys_params)
    vb.u .-= Uinf
    vb.v .-= Vinf

    return vb
end

function viscousflow_velocity_bc_rhs!(vb,sys::ILMSystem,t)
    prescribed_surface_average!(vb,t,sys)
    return vb
end

function viscousflow_vorticity_constraint_force!(dw,τ,sys)
    @unpack extra_cache, base_cache = sys
    @unpack dv = extra_cache

    viscousflow_velocity_constraint_force!(dv,τ,sys)
    curl!(dw,dv,base_cache)

    return dw
end

function viscousflow_velocity_constraint_force!(dv,τ,sys::ILMSystem{S,P,N}) where {S,P,N}
    @unpack base_cache = sys
    regularize!(dv,τ,base_cache)
    return dv
end

function viscousflow_velocity_constraint_force!(dv,τ,sys::ILMSystem{S,P,0}) where {S,P}
    return dv
end

function viscousflow_vorticity_bc_op!(vb,w,sys::ILMSystem)
    @unpack extra_cache, base_cache = sys
    @unpack velcache, v_tmp = extra_cache
    @unpack ψtemp = velcache

    vectorpotential_from_curlv!(ψtemp,w,base_cache)
    vecfield_from_vectorpotential!(v_tmp,ψtemp,base_cache)
    viscousflow_velocity_bc_op!(vb,v_tmp,sys)

    return vb
end

function viscousflow_velocity_bc_op!(vb,v,sys::ILMSystem)
    @unpack base_cache = sys
    interpolate!(vb,v,base_cache)
    return vb
end

#= Fields =#

vorticity!(masked_w::Nodes{Dual},w::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache) =
            masked_curlv_from_curlv_masked!(masked_w,w,dv,base_cache,wcache)

total_vorticity!(w::Nodes{Dual},masked_w::Nodes{Dual},dv::VectorData,base_cache::BasicILMCache,wcache::VectorPotentialCache) =
            curlv_masked_from_masked_curlv!(w,masked_w,dv,base_cache,wcache)

for f in [:vorticity,:total_vorticity]

    f! = Symbol(string(f)*"!")

    @eval function $f!(masked_w::Nodes{Dual},w::Nodes{Dual},τ,sys::ILMSystem,t)
        @unpack extra_cache, base_cache = sys
        @unpack dvb, velcache = extra_cache
        @unpack wcache = velcache

        prescribed_surface_jump!(dvb,t,sys)
        $f!(masked_w,w,dvb,base_cache,wcache)
    end

    @eval $f(w::Nodes{Dual},τ,sys::ILMSystem,t) = $f!(zeros_gridcurl(sys),w,τ,sys,t)

    @eval @snapshotoutput $f
end

function velocity!(v::Edges{Primal},curl_vmasked::Nodes{Dual},masked_divv::Nodes{Primal},dv::VectorData,vp,base_cache::BasicILMCache,velcache::VectorFieldCache,w_tmp::Nodes{Dual})
    @unpack wcache  = velcache

    masked_curlv = w_tmp

    masked_curlv_from_curlv_masked!(masked_curlv,curl_vmasked,dv,base_cache,wcache)
    vecfield_helmholtz!(v,masked_curlv,masked_divv,dv,vp,base_cache,velcache)

    return v
end

function velocity!(v::Edges{Primal},w::Nodes{Dual},sys::ILMSystem,t)
    @unpack forcing, phys_params, extra_cache, base_cache = sys
    @unpack dvb, velcache, divv_tmp, w_tmp = extra_cache

    prescribed_surface_jump!(dvb,t,sys)

    freestream_func = get_freestream_func(forcing)
    Vinf = freestream_func(t,phys_params)

    fill!(divv_tmp,0.0)
    velocity!(v,w,divv_tmp,dvb,Vinf,base_cache,velcache,w_tmp)
end

velocity(w::Nodes{Dual},τ,sys::ILMSystem,t) = velocity!(zeros_grid(sys),w,sys,t)

@snapshotoutput velocity

function streamfunction!(ψ::Nodes{Dual},w::Nodes{Dual},vp::Tuple,base_cache::BasicILMCache,wcache::VectorPotentialCache)
    @unpack stemp = wcache

    fill!(ψ,0.0)
    vectorpotential_from_curlv!(stemp,w,base_cache)
    ψ .+= stemp

    vectorpotential_uniformvecfield!(stemp,vp...,base_cache)
    ψ .+= stemp

    return ψ
end

function streamfunction!(ψ::Nodes{Dual},w::Nodes{Dual},sys::ILMSystem,t)
    @unpack phys_params, forcing, extra_cache, base_cache = sys
    @unpack velcache = extra_cache
    @unpack wcache = velcache

    freestream_func = get_freestream_func(forcing)
    Vinf = freestream_func(t,phys_params)

    streamfunction!(ψ,w,Vinf,base_cache,wcache)

end

streamfunction(w::Nodes{Dual},τ,sys::ILMSystem,t) = streamfunction!(zeros_gridcurl(sys),w,sys,t)

@snapshotoutput streamfunction

function pressure!(press::Nodes{Primal},w::Nodes{Dual},τ,sys::ILMSystem,t)
      @unpack extra_cache, base_cache = sys
      @unpack velcache, v_tmp, dv_tmp, dv, divv_tmp = extra_cache

      velocity!(v_tmp,w,sys,t)
      viscousflow_velocity_ode_rhs!(dv,v_tmp,sys,t)
      viscousflow_velocity_constraint_force!(dv_tmp,τ,sys)
      dv .-= dv_tmp

      divergence!(divv_tmp,dv,base_cache)
      inverse_laplacian!(press,divv_tmp,base_cache)
      return press
end

pressure(w::Nodes{Dual},τ,sys::ILMSystem,t) = pressure!(zeros_griddiv(sys),w,τ,sys,t)

@snapshotoutput pressure

#= Surface fields =#

function traction!(tract::VectorData{N},τ::VectorData{N},sys::ILMSystem,t) where {N}
    @unpack bc, extra_cache, base_cache, phys_params = sys
    @unpack vb_tmp, dvb = extra_cache
    @unpack sscalar_cache = base_cache

    prescribed_surface_average!(vb_tmp,t,sys)
    surface_velocity!(dvb,sys,t)
    vb_tmp .-= dvb

    nrm = normals(sys)
    pointwise_dot!(sscalar_cache,nrm,vb_tmp)
    prescribed_surface_jump!(dvb,t,sys)
    product!(tract,dvb,sscalar_cache)

    tract .+= τ
    return tract

end
traction!(out::VectorData{0},τ::VectorData{0},sys::ILMSystem,t) = out
traction(w::Nodes{Dual},τ::VectorData,sys::ILMSystem,t) = traction!(zeros_surface(sys),τ,sys,t)
@snapshotoutput traction

function pressurejump!(dpb::ScalarData{N},τ::VectorData{N},sys::ILMSystem,t) where {N}
    @unpack base_cache = sys
    @unpack sdata_cache = base_cache

    nrm = normals(sys)
    traction!(sdata_cache,τ,sys,t)
    pointwise_dot!(dpb,nrm,sdata_cache)
    dpb .*= -1.0
    return dpb
end
pressurejump!(out::VectorData{0},τ::VectorData{0},sys::ILMSystem,t) = out
pressurejump(w::Nodes{Dual},τ::VectorData,sys::ILMSystem,t) = pressurejump!(zeros_surfacescalar(sys),τ,sys,t)
@snapshotoutput pressurejump

#= Integrated metrics =#

force(w::Nodes{Dual},τ::VectorData{0},sys::ILMSystem{S,P,0},t,bodyi::Int) where {S,P} = nothing, nothing #Vector{Float64}(), Vector{Float64}()

function force(w::Nodes{Dual},τ::VectorData{N},sys::ILMSystem{S,P,N},t,bodyi::Int) where {S,P,N}
    @unpack base_cache = sys
    @unpack sdata_cache = base_cache
    traction!(sdata_cache,τ,sys,t)
    fx = integrate(sdata_cache.u,sys,bodyi)
    fy = integrate(sdata_cache.v,sys,bodyi)

    return fx, fy
end

@vectorsurfacemetric force

moment(w::Nodes{Dual},τ::VectorData{0},sys::ILMSystem{S,P,0},t,bodyi::Int;kwargs...) where {S,P} = nothing # Vector{Float64}()

function moment(w::Nodes{Dual},τ::VectorData{N},sys::ILMSystem{S,P,N},t,bodyi::Int;center=(0.0,0.0)) where {S,P,N}
    @unpack base_cache = sys
    @unpack sdata_cache, sscalar_cache = base_cache
    xc, yc = center
    pts = points(sys)
    pts.u .-= xc
    pts.v .-= yc

    traction!(sdata_cache,τ,sys,t)
    pointwise_cross!(sscalar_cache,pts,sdata_cache)
    mom = integrate(sscalar_cache,sys,bodyi)

    return mom
end

@scalarsurfacemetric moment

end


text:
"""
    Integrate1D(F::Function, Cube::HyperCube; tol::Real=1e-14, FullSol::Bool=false, meth=nothing)
Integrates `F` over a one-dimensional domain specified via a `HyperCube` by rephrasing the integral as an ODE and using `DifferentialEquations.jl`.
"""
function Integrate1D(F::Function, Cube::HyperCube; tol::Real=1e-14, FullSol::Bool=false, meth=nothing)
    length(Cube) != 1 && throw(ArgumentError("Cube dim = $(length(Cube)) instead of 1"))
    Integrate1D(F,(Cube.L[1],Cube.U[1]); tol=tol,FullSol=FullSol,meth=meth)
end
Integrate1D(F::Function, Interval::AbstractVector{<:Number}; tol::Real=1e-14, FullSol::Bool=false, meth=nothing) = Integrate1D(F, Tuple(Interval); tol=tol, FullSol=FullSol, meth=meth)
function Integrate1D(F::Function, Interval::Tuple{<:Number,<:Number}; tol::Real=1e-14, FullSol::Bool=false, meth=nothing)
    Interval = floatify(Interval)
    !(0. < tol < 1.) && throw("Integrate1D: tol unsuitable")
    Interval[1] > Interval[2] && throw(ArgumentError("Interval orientation wrong."))
    f(u,p,t) = F(t);    u0 = 0.
    if tol < 1e-15
        u0 = BigFloat(u0);        Interval = BigFloat.(Interval)
        meth = isnothing(meth) ? Vern9() : meth
    else
        meth = isnothing(meth) ? Tsit5() : meth
    end
    if FullSol
        return solve(ODEProblem(f,u0,Interval),meth; reltol=tol,abstol=tol)
    else
        return solve(ODEProblem(f,u0,Interval),meth; reltol=tol,abstol=tol,save_everystep=false,save_start=false,save_end=true).u[end]
    end
end

"""
    IntegrateND(F::Function,Cube::HyperCube; tol::Real=1e-12, WE::Bool=false, kwargs...)
Integrates the function `F` over `Cube` with the help of **HCubature.jl** to a tolerance of `tol`.
If `WE=true`, the result is returned as a `Measurement` which also contains the estimated error in the result.
"""
function IntegrateND(F::Function, Cube::HyperCube; tol::Real=1e-12, WE::Bool=false, kwargs...)
    if length(Cube) == 1
        val, uncert = hquadrature(F, Cube.L[1], Cube.U[1]; rtol=tol, atol=tol, kwargs...)
    else
        val, uncert = hcubature(F, Cube.L, Cube.U; rtol=tol, atol=tol, kwargs...)
    end
    if length(val) == 1
        return WE ? measurement(val[1],uncert[1]) : val[1]
    else
        return WE ? measurement.(val,uncert) : val
    end
end
IntegrateND(F::Function, L::AbstractVector{<:Number}, U::AbstractVector{<:Number}; tol::Real=1e-12, WE::Bool=false, kwargs...) = IntegrateND(F,HyperCube(L,U); tol=tol, WE=WE, kwargs...)
IntegrateND(F::Function, Interval::Union{AbstractVector{<:Number},Tuple{<:Number,<:Number}}; tol::Real=1e-12, WE::Bool=false, kwargs...) = IntegrateND(F,HyperCube(Interval); tol=tol, WE=WE, kwargs...)

"""
    IntegrateOverConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, Confnum::Real, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
Integrates a function `F` over the intersection of `Domain` and the confidence region of level `Confnum`.
"""
function IntegrateOverConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, Confnum::Real, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
    @assert length(Domain) == pdim(DM)
    # Multiply F with characteristic function for confidence region
    Threshold = LogLikeMLE(DM) - 0.5InvChisqCDF(pdim(DM), ConfVol(Confnum))
    InsideRegion(X::AbstractVector{<:Number}) = loglikelihood(DM, X; kwargs...) > Threshold
    Integrand(X::AbstractVector{<:Number}) = InsideRegion(X) ? F(X) : zero(suff(X))
    # Use HCubature instead of MonteCarlo
    MonteCarloArea(Integrand, Domain, N; WE=WE)
end

"""
    IntegrateOverApproxConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, sol::AbstractODESolution, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
    IntegrateOverApproxConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
Integrates a function `F` over the intersection of `Domain` and the polygon defined by `sol`.
"""
function IntegrateOverApproxConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, sol::AbstractODESolution, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
    @assert length(Domain) == pdim(DM) == length(sol.u[1]) == 2
    Integrand(X::AbstractVector{<:Number}) = ApproxInRegion(sol, X) ? F(X) : zero(suff(X))
    # Use HCubature instead of MonteCarlo
    MonteCarloArea(Integrand, Domain, N; WE=WE)
end

function IntegrateOverApproxConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, Tup::Tuple{<:AbstractVector{<:Plane},<:AbstractVector{<:AbstractODESolution}}, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
    IntegrateOverApproxConfidenceRegion(DM, Domain, Tup[1], Tup[2], F; N=N, WE=WE, kwargs...)
end
function IntegrateOverApproxConfidenceRegion(DM::AbstractDataModel, Domain::HyperCube, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, F::Function; N::Int=Int(1e5), WE::Bool=true, kwargs...)
    @assert length(Domain) == pdim(DM) == ConsistentElDims(Planes)
    @assert length(Planes) == length(sols)
    Integrand(X::AbstractVector{<:Number}) = ApproxInRegion(Planes, sols, X) ? F(X) : zero(suff(X))
    # Use HCubature instead of MonteCarlo
    MonteCarloArea(Integrand, Domain, N; WE=WE)
end

function MonteCarloArea(Test::Function,Cube::HyperCube,N::Int=Int(1e7); WE::Bool=false)
    if WE
        return CubeVol(Cube) * MonteCarloRatioWE(Test,Cube,N)
    else
        return CubeVol(Cube) * MonteCarloRatio(Test,Cube,N)
    end
end
function MonteCarloRatio(Test::Function,Cube::HyperCube,N::Int=Int(1e7))
    (1/N)* @distributed (+) for i in 1:N
        Test(rand.(Uniform.(Cube.L,Cube.U)))
    end
end

function MonteCarloRatioWE(Test::Function,LU::HyperCube,N::Int=Int(1e7); chunksize::Int=Int(N/20))
    chunksize > N && error("chunksize > N")
    if N%chunksize != 0
        N += Int(N%chunksize + 1)
    end
    chunks = Int(N/chunksize)
    # Output not normalized by chunksize
    function CarloLoop(Test::Function,LU::HyperCube,chunksize::Int)
        tot = [rand.(Uniform.(LU.L,LU.U)) for i in 1:chunksize] .|> Test
        res = sum(tot)
        [res, sum(abs2, (tot .- (res/chunksize)))]
    end
    Tot = @distributed (+) for i in 1:chunks
        CarloLoop(Test,LU,chunksize)
    end
    measurement(Tot[1]/N, sqrt(1/((N-1)*N) * Tot[2]))
end

# # From Cuba.jl docs
# function CompactDomainTransform(F::ModelOrFunction, Cube::HyperCube)
#     (!all(x->isfinite(x),Cube.L) || !all(x->isfinite(x),Cube.U)) && throw("Not applicable.")
#     if length(Cube) == 1
#         W = Cube.U[1] - Cube.L[1]
#         return x -> W * F(Cube.L[1] + W * x)
#         # Use mul! or something like that?
#     else
#         W = CubeWidths(Cube);   V = prod(CubeWidths)
#         return x -> V * F(Cube.L + W * x)
#     end
# end
#
# function HalfInfiniteTransform(F::ModelOrFunction, Cube::HyperCube)
#     # integral over [a,∞]
#     if Cube.U[1] == Inf && isfinite(Cube.L[1])
#         return x -> (1-x)^-2 * F(Cube.L[1] + x/(1-x))
#     end
# end
#
# function InfiniteDomainTransform(F::ModelOrFunction, Cube::HyperCube)
#     if Cube.L[1] == -Inf && Cube.L[1] == Inf
#         return x -> F((2x - 1.)/((1 - x)*x)) * (2x^2 - 2y + 1) / ((1-x^2)*x^2)
#     end
# end


text:
#!/usr/bin/julia

# Trizen
# Date: 18 August 2016
# https://github.com/trizen

# Count the number of partitions of n.

# See also: https://oeis.org/A000041
#           https://en.wikipedia.org/wiki/Partition_(number_theory)

function count_partitions(x::Int64)
    n = 2
    p = Int64[1]

    while (n <= x+1)
        i = 0
        q = 2
        push!(p, 0)

        while q <= n
            p[n] += (i % 4 > 1 ? -1 : 1) * p[n-q+1]
            i += 1
            j = div(i, 2) + 1
            isodd(i) && (j *= -1)
            q = div(j * (3j - 1), 2) + 1
        end

        n += 1
    end

    p[n-1]
end

println("P(200) = ", count_partitions(200))        # 3972999029388

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
@testset "node construction" begin
    dist = Normal()
    address = :z
    n = node(Float64, address, dist, false, NONSTANDARD)
    @test n.logprob == 0.0
end

@testset "trace construction" begin
    dist = Normal()
    address = :z
    n = node(Float64, address, dist, false, NONSTANDARD)
    t = trace()
    t[address] = n
    @test t[address] == n
end

@testset "sample construction" begin
    t = trace()
    s_d = sample(t, :count, Poisson(3.0))
    @info "Sampled $s_d"
    @info "Sample node is $(t[:count])"
    @test t[:count].value == s_d
end

function program!(t :: Trace, n :: Int)
    loc = sample(t, :loc, Normal())
    scale = sample(t, :scale, Gamma(2.0, 2.0))
    log_rate = sample(t, :log_rate, Normal(loc, scale))
    for i in 1:n
        sample(t, (:data, i), Poisson(exp(log_rate)))
    end
end

@testset "tracing program" begin
    t = trace()
    n_datapoints = 10
    program!(t, n_datapoints)
    logprob!(t)
    @info "Joint density of execution = $(t.logprob_sum)"
    @test :loc in keys(t)
    @test (:data, 3) in keys(t)
end