text:
import Base: strip
abstract type Runner; end

export default_host_platform

# Host platform _must_ match the C++ string ABI of the binaries we get from the
# repositories.  Note: when preferred_gcc_version=v"4" we can't really build for
# that C++ string ABI :-(
"""
    default_host_platform

The default host platform in the build environment.
"""
const default_host_platform = Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11")

function nbits(p::AbstractPlatform)
    if arch(p) in ("i686", "armv6l", "armv7l")
        return 32
    elseif arch(p) in ("x86_64", "aarch64", "powerpc64le")
        return 64
    else
        error("Unknown bitwidth for architecture $(arch(p))")
    end
end

function proc_family(p::AbstractPlatform)
    if arch(p) in ("x86_64", "i686")
        return "intel"
    elseif arch(p) in ("armv6l", "armv7l", "aarch64")
        return "arm"
    elseif arch(p) == "powerpc64le"
        return "power"
    else
        error("Unknown processor family for architecture $(arch(p))")
    end
end

# Convert platform to a triplet, but strip out the ABI parts.
# Also translate `armvXl` -> `arm` for now, since that's what most
# compiler toolchains call it.  :(
function aatriplet(p::AbstractPlatform)
    t = triplet(abi_agnostic(p))
    t = replace(t, "armv7l" => "arm")
    t = replace(t, "armv6l" => "arm")
    return t
end
# XXX: we want AnyPlatform to look like `x86_64-linux-musl` in the build environment.
aatriplet(p::AnyPlatform) = aatriplet(default_host_platform)

"""
    generate_compiler_wrappers!(platform::AbstractPlatform; bin_path::AbstractString,
                                host_platform::AbstractPlatform = $(repr(default_host_platform)),
                                compilers::Vector{Symbol} = [:c],
                                allow_unsafe_flags::Bool = false,
                                lock_microarchitecture::Bool = true)

We generate a set of compiler wrapper scripts within our build environment to force all
build systems to honor the necessary sets of compiler flags to build for our systems.
Note that while `platform_envs()` sets many environment variables, those values are
intended to be optional/overridable.  These values, while still overridable by directly
invoking a compiler binary directly (e.g. /opt/{target}/bin/{target}-gcc), are much more
difficult to override, as the flags embedded in these wrappers are absolutely necessary,
and even simple programs will not compile without them.
"""
function generate_compiler_wrappers!(platform::AbstractPlatform; bin_path::AbstractString,
                                     host_platform::AbstractPlatform = default_host_platform,
                                     compilers::Vector{Symbol} = [:c],
                                     allow_unsafe_flags::Bool = false,
                                     lock_microarchitecture::Bool = true,
                                     bootstrap::Bool = !isempty(bootstrap_list),
                                     )
    # Wipe that directory out, in case it already had compiler wrappers
    rm(bin_path; recursive=true, force=true)
    mkpath(bin_path)

    # Early-exit if we're bootstrapping
    if bootstrap
        return
    end

    target = aatriplet(platform)
    host_target = aatriplet(host_platform)

    function wrapper(io::IO,
                     prog::String;
                     # Flags that are always prepended
                     flags::Vector{String} = String[],
                     # Flags that are prepended if we think we're compiling (e.g. no `-x assembler`)
                     compile_only_flags::Vector = String[],
                     # Flags that are postpended if we think we're linking (e.g. no `-c`)
                     link_only_flags::Vector = String[],
                     allow_ccache::Bool = true,
                     no_soft_float::Bool = false,
                     hash_args::Bool = false,
                     extra_cmds::String = "",
                     env::Dict{String,String} = Dict{String,String}(),
                     unsafe_flags = String[])
        write(io, """
        #!/bin/bash
        # This compiler wrapper script brought into existence by `generate_compiler_wrappers!()`

        if [ "x\${SUPER_VERBOSE}" = "x" ]; then
            vrun() { "\$@"; }
        else
            vrun() { echo -e "\\e[96m\$@\\e[0m" >&2; "\$@"; }
        fi

        ARGS=( "\$@" )
        PRE_FLAGS=()
        POST_FLAGS=()
        """)

        # Sometimes we need to look at the hash of our arguments
        if hash_args
            write(io, """
            ARGS_HASH="\$(echo -n "\$*" | sha1sum | cut -c1-8)"
            """)
        end

        # If we're given always-prepend flags, include them
        if !isempty(flags)
            println(io)
            for cf in flags
                println(io, "PRE_FLAGS+=( $cf )")
            end
            println(io)
        end

        # If we're given compile-only flags, include them only if `-x assembler` is not provided
        if !isempty(compile_only_flags)
            println(io)
            println(io, "if [[ \" \${ARGS[@]} \" != *' -x assembler '* ]]; then")
            for cf in compile_only_flags
                println(io, "    PRE_FLAGS+=( $cf )")
            end
            println(io, "fi")
            println(io)
        end

        # If we're given link-only flags, include them only if `-c` or other link-disablers are not provided.
        if !isempty(link_only_flags)
            println(io)
            println(io, "if [[ \" \${ARGS[@]} \" != *' -c '* ]] && [[ \" \${ARGS[@]} \" != *' -E '* ]] && [[ \" \${ARGS[@]} \" != *' -M '* ]] && [[ \" \${ARGS[@]} \" != *' -fsyntax-only '* ]]; then")
            for lf in link_only_flags
                println(io, "    POST_FLAGS+=( $lf )")
            end
            println(io, "fi")
            println(io)
        end

        # Insert extra commands from the user (usually some kind of conditional setting
        # of PRE_FLAGS and POST_FLAGS)
        println(io)
        write(io, extra_cmds)
        println(io)

        for (name, val) in env
            write(io, "export $(name)=\"$(val)\"\n")
        end

        # TODO: improve this check
        if lock_microarchitecture
            write(io, raw"""
                      if [[ " ${ARGS[@]} " == *"-march="* ]]; then
                          echo "Cannot force an architecture" >&2
                          exit 1
                      fi
                      """)
            println(io)
        end

        if no_soft_float
            write(io, raw"""
                      if [[ " ${ARGS[@]} " == *"-mfloat-abi=soft"* ]]; then
                          echo "${target} platform does not support soft-float ABI" >&2
                          exit 1
                      fi
                      """)
            println(io)
        end

        if length(unsafe_flags) >= 1
            write(io, """
            if [[ "\${ARGS[@]}" =~ \"$(join(unsafe_flags, "\"|\""))\" ]]; then
                echo -e \"You used one or more of the unsafe flags: $(join(unsafe_flags, ", "))\\nPlease repent.\" >&2
                exit 1
            fi
            """)
            println(io)
        end

        if allow_ccache
            write(io, """
            if [ \${USE_CCACHE} == "true" ]; then
                CCACHE="ccache"
            fi
            """)
        end
        write(io, """
        vrun \${CCACHE} $(prog) "\${PRE_FLAGS[@]}" "\${ARGS[@]}" "\${POST_FLAGS[@]}"
        """)
    end

    # Helper invocations
    target_tool(io::IO, tool::String, args...; kwargs...) = wrapper(io, "/opt/$(target)/bin/$(target)-$(tool)", args...; kwargs...)
    llvm_tool(io::IO, tool::String, args...; kwargs...) = wrapper(io, "/opt/$(host_target)/bin/llvm-$(tool)", args...; kwargs...)

    # For now this is required for Clang, since apple spells aarch64 as "arm64".
    # Should probably be fixed upstream, but will do for now
    clang_target_triplet(p::AbstractPlatform) = replace(aatriplet(p), "aarch64" => "arm64")

    function clang_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        # Focus the clang targeting laser
        append!(flags, [
            # Set the `target` for `clang` so it generates the right kind of code
            "-target $(clang_target_triplet(p))",
            # Set our sysroot to the platform-specific location, dropping compiler ABI annotations
            "--sysroot=/opt/$(aatriplet(p))/$(aatriplet(p))/sys-root",
        ])
        # For MacOS and FreeBSD, we don't set `-rtlib`, and FreeBSD is special-cased within the LLVM source tree
        # to not allow for -gcc-toolchain, which means that we have to manually add the location of libgcc_s.  LE SIGH.
        # We do that within `clang_linker_flags()`, so that we don't get "unused argument" warnings all over the place.
        # https://github.com/llvm-mirror/clang/blob/f3b7928366f63b51ffc97e74f8afcff497c57e8d/lib/Driver/ToolChains/FreeBSD.cpp
        # For everything else, we provide `-rtlib=libgcc` because clang-builtins are broken (pending Valentin-based-magic),
        # and we also need to provide `-stdlib=libstdc++` to match Julia on these platforms.
        if !Sys.isbsd(p)
            append!(flags, [
                # Find GCC toolchain here (for things like libgcc_s)
                "--gcc-toolchain=/opt/$(aatriplet(p))"
                # Use libgcc as the C runtime library
                "-rtlib=libgcc"
                # Use libstdc++ as the C++ runtime library
                "-stdlib=libstdc++"
            ])
        end
        return flags
    end

    function min_macos_version_flag(p::AbstractPlatform)
        # If no `os_version` is specified in `p`, default to the oldest we support in the Julia world,
        # which is `10.8`, but if it is actually specified, then set that corresponding value.
        #min_macos_version = something(os_version(p), v"14.0.0")

        # Eventually, we'll take this in `os_version(p)`, but not just yet.  We need to fix the paths
        # to the compiler shards first, since right now they have `14` at the end
        min_macos_version = v"14.0.0"

        kernel_to_macos = Dict(
            12 => "10.8",
            13 => "10.9",
            14 => "10.10",
            15 => "10.11",
            16 => "10.12",
            17 => "10.13",
            18 => "10.14",
            19 => "10.15",
            20 => "11.0",
        )

        # Ask compilers to compile for a minimum macOS version
        return "-mmacosx-version-min=$(kernel_to_macos[min_macos_version.major])"
    end

    function clang_compile_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        if lock_microarchitecture
            append!(flags, get_march_flags(arch(p), march(p), "clang"))
        end
        if Sys.isapple(p)
            append!(flags, String[
                # On MacOS, we need to override the typical C++ include search paths, because it always includes
                # the toolchain C++ headers first.  Valentin tracked this down to:
                # https://github.com/llvm/llvm-project/blob/0378f3a90341d990236c44f297b923a32b35fab1/clang/lib/Driver/ToolChains/Darwin.cpp#L1944-L1978
                "-nostdinc++",
                "-isystem",
                "/opt/$(aatriplet(p))/$(aatriplet(p))/sys-root/usr/include/c++/v1",
                # We also add `-Wno-unused-command-line-argument` so that if someone does something like
                # `clang -Werror -o foo a.o b.o`, it doesn't complain due to the fact that that is using
                # `clang` as a linker (and we have no real way to detect that in the wrapper), which will
                # cause `clang` to complain about compiler flags being passed in.
                "-Wno-unused-command-line-argument",
                min_macos_version_flag(p),
            ])
        end
        return flags
    end

    function clang_link_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        # On macos and freebsd, we must pass in the `/lib` directory for some reason
        if Sys.isbsd(p)
            push!(flags, "-L/opt/$(aatriplet(p))/$(aatriplet(p))/lib")
        end

        # we want to use a particular linker with clang.  But we want to avoid warnings about unused
        # flags when just compiling, so we put it into "linker-only flags".
        push!(flags, "-fuse-ld=$(aatriplet(p))")

        # On macos, we need to pass `-headerpad_max_install_names` so that we have lots of space
        # for `install_name_tool` shenanigans during audit fixups.
        if Sys.isapple(p)
            push!(flags, "-headerpad_max_install_names")
        end
        return flags
    end

    function macos_gcc_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        # On macOS, if we're on an old GCC, the default -syslibroot that gets
        # passed to the linker isn't calculated correctly, so we have to manually set it.
        gcc_version, llvm_version = select_compiler_versions(p)
        if gcc_version.major in (4, 5)
            push!(flags, "-Wl,-syslibroot,/opt/$(aatriplet(p))/$(aatriplet(p))/sys-root")
        end
        push!(flags, min_macos_version_flag(p))
        return flags
    end

    function gcc_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        # Force proper cxx11 string ABI usage w00t w00t!
        if cxxstring_abi(p) == "cxx11"
            push!(flags, "-D_GLIBCXX_USE_CXX11_ABI=1")
        elseif cxxstring_abi(p) == "cxx03"
            push!(flags, "-D_GLIBCXX_USE_CXX11_ABI=0")
        end

        # Use hash of arguments to provide consistent, unique random seed
        push!(flags, "-frandom-seed=0x\${ARGS_HASH}")

        if Sys.isapple(p)
            macos_gcc_flags!(p, flags)
        end
        return flags
    end

    function gcc_compile_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        if lock_microarchitecture
            append!(flags, get_march_flags(arch(p), march(p), "gcc"))
        end
        return flags
    end

    function gcc_link_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        # Yes, it does seem that the inclusion of `/lib64` on `powerpc64le` was fixed
        # in GCC 6, broken again in GCC 7, and then fixed again for GCC 8 and 9
        gcc_version, llvm_version = select_compiler_versions(p)
        if arch(p) == "powerpc64le" && Sys.islinux(p) && gcc_version.major in (4, 5, 7)
            append!(flags, String[
                "-L/opt/$(aatriplet(p))/$(aatriplet(p))/sys-root/lib64",
                "-Wl,-rpath-link,/opt/$(aatriplet(p))/$(aatriplet(p))/sys-root/lib64",
            ])
        elseif Sys.isapple(p)
            push!(flags, "-headerpad_max_install_names")
        end
        return flags
    end

    function gcc_unsafe_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        if !allow_unsafe_flags
            return String["-Ofast", "-ffast-math", "-funsafe-math-optimizations"]
        end
        return String[]
    end

    function gcc_wrapper(io::IO, tool::String, p::AbstractPlatform, allow_ccache::Bool = true)
        return wrapper(io,
            "/opt/$(aatriplet(p))/bin/$(aatriplet(p))-$(tool)";
            flags=gcc_flags!(p),
            compile_only_flags=gcc_compile_flags!(p),
            link_only_flags=gcc_link_flags!(p),
            unsafe_flags=gcc_unsafe_flags!(p),
            hash_args = true,
            allow_ccache,
            no_soft_float=arch(p) in ("armv6l", "armv7l"),
        )
    end

    function clang_wrapper(io::IO, tool::String, p::AbstractPlatform, extra_flags::Vector{String} = String[])
        flags = clang_flags!(p)
        append!(flags, extra_flags)
        return wrapper(io,
            "/opt/$(host_target)/bin/$(tool)";
            flags=flags,
            compile_only_flags=clang_compile_flags!(p),
            link_only_flags=clang_link_flags!(p),
            no_soft_float=arch(p) in ("armv6l", "armv7l"),
        )
    end

    # C/C++/Fortran
    gcc(io::IO, p::AbstractPlatform)      = gcc_wrapper(io, "gcc", p)
    gxx(io::IO, p::AbstractPlatform)      = gcc_wrapper(io, "g++", p)
    gfortran(io::IO, p::AbstractPlatform) = gcc_wrapper(io, "gfortran", p, false)

    clang(io::IO, p::AbstractPlatform)    = clang_wrapper(io, "clang", p)
    clangxx(io::IO, p::AbstractPlatform)  = clang_wrapper(io, "clang++", p)
    objc(io::IO, p::AbstractPlatform)     = clang_wrapper(io, "clang", p, ["-x objective-c"])

    # Our general `cc`  points to `gcc` for most systems, but `clang` for MacOS and FreeBSD
    function cc(io::IO, p::AbstractPlatform)
        if Sys.isbsd(p)
            return clang(io, p)
        else
            return gcc(io, p)
        end
    end
    function cxx(io::IO, p::AbstractPlatform)
        if Sys.isbsd(p)
            return clangxx(io, p)
        else
            return gxx(io, p)
        end
    end
    fc(io::IO, p::AbstractPlatform) = gfortran(io, p)

    # Go stuff where we build an environment mapping each time we invoke `go-${target}`
    function GOOS(p::AbstractPlatform)
        if os(p) == "macos"
            return "darwin"
        end
        return os(p)
    end
    function GOARCH(p::AbstractPlatform)
        arch_mapping = Dict(
            "armv6l" => "arm",
            "armv7l" => "arm",
            "aarch64" => "arm64",
            "x86_64" => "amd64",
            "i686" => "386",
            "powerpc64le" => "ppc64le",
        )
        return arch_mapping[arch(p)]
    end
    function go(io::IO, p::AbstractPlatform)
        env = Dict(
            "GOOS" => GOOS(p),
            "GOROOT" => "/opt/$(host_target)/go",
            "GOARCH" => GOARCH(p),
        )
        return wrapper(io, "/opt/$(host_target)/go/bin/go"; env=env, allow_ccache=false)
    end

    # Rust stuff
    function rust_flags!(p::AbstractPlatform, flags::Vector{String} = String[])
        push!(flags, "--target=$(map_rust_target(p))")
        if Sys.islinux(p)
            push!(flags, "-Clinker=$(aatriplet(p))-gcc")

            # Add aarch64 workaround https://github.com/rust-lang/rust/issues/46651#issuecomment-402850885
            if arch(p) == "aarch64" && libc(p) == "musl"
                push!(flags, "-C link-arg=-lgcc")
            end
        elseif Sys.iswindows(p)
            # Rust on i686 mingw32 can't deal with exceptions
            # https://github.com/rust-lang/rust/issues/12859
            if arch(p) == "i686"
                push!(flags, "-C panic=abort")
            end
        end
        return flags
    end
    rustc(io::IO, p::AbstractPlatform) = wrapper(io, "/opt/$(host_target)/bin/rustc"; flags=rust_flags!(p), allow_ccache=false)
    rustup(io::IO, p::AbstractPlatform) = wrapper(io, "/opt/$(host_target)/bin/rustup"; allow_ccache=false)
    cargo(io::IO, p::AbstractPlatform) = wrapper(io, "/opt/$(host_target)/bin/cargo"; allow_ccache=false)

    # Meson REQUIRES that `CC`, `CXX`, etc.. are set to the host utils.  womp womp.
    function meson(io::IO, p::AbstractPlatform)
        meson_env = Dict(
            "AR"     => "$(host_target)-ar",
            "CC"     => "$(host_target)-cc",
            "CXX"    => "$(host_target)-c++",
            "FC"     => "$(host_target)-f77",
            "LD"     => "$(host_target)-ld",
            "NM"     => "$(host_target)-nm",
            "OBJC"   => "$(host_target)-objc",
            "RANLIB" => "$(host_target)-ranlib",
        )
        wrapper(io, "/usr/bin/meson"; allow_ccache=false, env=meson_env)
    end

    # Patchelf needs some page-alignment on aarch64 and ppc64le forced into its noggin
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/commit/cce4f8fdbb16425d245ab87a50f60d1a16d04948
    function patchelf(io::IO, p::AbstractPlatform)
        extra_cmds = ""
        if Sys.islinux(p) && arch(p) in ("aarch64", "powerpc64le")
            extra_cmds = raw"""
            if [[ " ${ARGS[@]} " != *'--page-size'* ]]; then
                PRE_FLAGS+=( '--page-size' '65536' )
            fi
            """
        end
        wrapper(io, "/usr/bin/patchelf"; allow_ccache=false, extra_cmds=extra_cmds)
    end

    # We pass `-D` to all `ar` invocations (unless `-U` is explicitly passed) for reproducibility
    function ar(io::IO, p::AbstractPlatform)
        ar_name = string(aatriplet(p), "-ar")
        if Sys.isapple(p)
            ar_name = "llvm-ar"
        end
        extra_cmds = raw"""
        if [[ " ${ARGS[0]} " =~ --* ]]; then
            # do nothing, it's probably --version or something
            true
        elif [[ " ${ARGS[0]} " != *'U'* ]]; then
            # Eliminate the `u` option, as it's incompatible with `D` and is just an optimization
            if [[ " ${ARGS[0]} " == *'u'* ]]; then
                ARGS[0]=$(echo "${ARGS[0]}" | tr -d u)
            fi

            # Add -D for "Deterministic mode"
            ARGS[0]="${ARGS[0]}D"
        else
            echo "Non-reproducibility alert: This 'ar' invocation uses the '-U' flag which embeds timestamps." >&2
            echo "ar flags: ${ARGS[@]}" >&2
            echo "Continuing build, but please repent." >&2
        fi
        """
        wrapper(io, string("/opt/", aatriplet(p), "/bin/", ar_name); allow_ccache=false, extra_cmds=extra_cmds)
    end

    function ranlib(io::IO, p::AbstractPlatform)
        if !Sys.isapple(p)
            ranlib_name = string(aatriplet(p), "-ranlib")
            extra_cmds = raw"""
            if [[ " ${ARGS[@]} " =~ "-[hHvVt]*U" ]]; then
                echo "Non-reproducibility alert: This `ranlib` invocation uses the `-U` flag which embeds timestamps." >&2
                echo "ranlib flags: ${ARGS[@]}" >&2
                echo "Continuing build, but please repent." >&2
            else
                PRE_FLAGS+=( '-D' )
            fi
            """
        else
            # llvm-ranlib is always reproducible
            ranlib_name = "llvm-ranlib"
            extra_cmds = ""
        end
        wrapper(io, string("/opt/", aatriplet(p), "/bin/", ranlib_name); allow_ccache=false, extra_cmds=extra_cmds)
    end

    function dlltool(io::IO, p::AbstractPlatform)
        extra_cmds = raw"""
        PRE_FLAGS+=( --temp-prefix /tmp/dlltool-${ARGS_HASH} )
        """
        wrapper(io, string("/opt/", aatriplet(p), "/bin/", string(aatriplet(p), "-dlltool")); allow_ccache=false, extra_cmds=extra_cmds, hash_args=true)
    end

    # Write out a bunch of common tools
    for tool in (:as, :cpp, :ld, :nm, :libtool, :objcopy, :objdump, :otool,
                 :strip, :install_name_tool, :dlltool, :windres, :winmc, :lipo)
        @eval $(tool)(io::IO, p::AbstractPlatform) = $(wrapper)(io, string("/opt/", aatriplet(p), "/bin/", aatriplet(p), "-", $(string(tool))); allow_ccache=false)
    end

    # c++filt is hard to write in symbols
    function cxxfilt(io::IO, p::AbstractPlatform)
        if Sys.isapple(p)
            # We must use `llvm-cxxfilt` on MacOS
            path = "/opt/$(aatriplet(p))/bin/llvm-cxxfilt"
        else
            path = "/opt/$(aatriplet(p))/bin/$(aatriplet(p))-c++filt"
        end
        return wrapper(io, path; allow_ccache=false)
    end

    function dsymutil(io::IO, p::AbstractPlatform)
        if !Sys.isapple(p)
            # Nobody except macOS has a `dsymutil`
            return (io, p) -> nothing
        end

        return wrapper(io, "/opt/$(aatriplet(p))/bin/dsymutil"; allow_ccache=false)
    end

    function readelf(io::IO, p::AbstractPlatform)
        if Sys.isapple(p)
            # macOS doesn't have a readelf
            return (io, p) -> nothing
        end
        return wrapper(io, "/opt/$(aatriplet(p))/bin/$(aatriplet(p))-readelf"; allow_ccache=false)
    end

    function write_wrapper(wrappergen, p, fname)
        file_path = joinpath(bin_path, triplet(p), fname)
        mkpath(dirname(file_path))
        open(io -> Base.invokelatest(wrappergen, io, p), file_path, "w")
        chmod(file_path, 0o775)
    end

    ## Generate compiler wrappers for both our host and our target.
    for p in unique((host_platform, platform))
        t = aatriplet(p)

        # Generate `:c` compilers
        if :c in compilers
            write_wrapper(cc, p, "$(t)-cc")
            write_wrapper(cxx, p, "$(t)-c++")

            # Generate `gcc`, `g++`, `clang` and `clang++`
            write_wrapper(gcc, p, "$(t)-gcc")
            write_wrapper(gxx, p, "$(t)-g++")
            write_wrapper(clang, p, "$(t)-clang")
            write_wrapper(clangxx, p, "$(t)-clang++")
            write_wrapper(objc, p, "$(t)-objc")

            # Someday, you will be split out
            write_wrapper(gfortran, p, "$(t)-f77")
            write_wrapper(gfortran, p, "$(t)-gfortran")
        end

        # Binutils (we always do these)
        write_wrapper(ar, p, "$(t)-ar")
        write_wrapper(as, p, "$(t)-as")
        write_wrapper(cpp, p, "$(t)-cpp")
        write_wrapper(cxxfilt, p, "$(t)-c++filt")
        write_wrapper(ld, p, "$(t)-ld")
        # ld wrappers for clang's `-fuse-ld=$(target)`
        if Sys.isapple(p)
            write_wrapper(ld, p, "ld64.$(t)")
        else
            write_wrapper(ld, p, "ld.$(t)")
        end
        write_wrapper(nm, p, "$(t)-nm")
        write_wrapper(libtool, p, "$(t)-libtool")
        write_wrapper(objcopy, p, "$(t)-objcopy")
        write_wrapper(objdump, p, "$(t)-objdump")
        write_wrapper(ranlib, p, "$(t)-ranlib")
        write_wrapper(readelf, p, "$(t)-readelf")
        write_wrapper(strip, p, "$(t)-strip")

        # Special mac stuff
        if Sys.isapple(p)
            write_wrapper(install_name_tool, p, "$(t)-install_name_tool")
            write_wrapper(lipo, p, "$(t)-lipo")
            write_wrapper(dsymutil, p, "$(t)-dsymutil")
            write_wrapper(otool, p, "$(t)-otool")
        end

        # Special Windows stuff
        if Sys.iswindows(p)
            write_wrapper(dlltool, p, "$(t)-dlltool")
            write_wrapper(windres, p, "$(t)-windres")
            write_wrapper(winmc, p, "$(t)-winmc")
        end

        # Generate go stuff
        if :go in compilers
            write_wrapper(go, p, "$(t)-go")
        end

        # Misc. utilities
        write_wrapper(patchelf, p, "$(t)-patchelf")
    end

    # Rust stuff doesn't use the normal "host" platform, it uses x86_64-linux-gnu, so we always have THREE around,
    # because clever build systems like `meson` ask Rust what its native system is, and it truthfully answers
    # `x86_64-linux-gnu`, while other build systems might say `x86_64-linux-musl` with no less accuracy.  So for
    # safety, we just ship all three all the time.
    if :rust in compilers
        for p in unique((platform, host_platform))
            t = aatriplet(p)
            write_wrapper(rustc, p, "$(t)-rustc")
            write_wrapper(rustup, p, "$(t)-rustup")
            write_wrapper(cargo, p, "$(t)-cargo")

            # For FreeBSD and macOS we need to create an unversioned link for
            # gcc because that's the linker our Rust toolchain expects:
            # https://github.com/JuliaPackaging/Yggdrasil/blob/fff0583bc2d8f32e450c427684f295524f38535d/0_RootFS/Rust/build_tarballs.jl#L115-L126.
            if Sys.isbsd(p) && os_version(p) !== nothing
                tmp_p = deepcopy(p)
                delete!(tags(tmp_p), "os_version")
                symlink("$(t)-gcc", joinpath(bin_path, "$(aatriplet(tmp_p))-gcc"))
            end
            # Currently our Rust toolchain expects the linker for armv7l and
            # armv6l with the platform "*l" suffix in the platform.  Until
            # https://github.com/JuliaPackaging/Yggdrasil/pull/2168 makes it to
            # the Rust toolchain, we create a symlink to work around this issue.
            if proc_family(p) == "arm" && nbits(p) == 32
                symlink("$(t)-gcc", joinpath(bin_path, "$(triplet(abi_agnostic(p)))-gcc"))
            end
        end
    end

    # Write a single wrapper for `meson`
    write_wrapper(meson, host_platform, "meson")

    default_tools = [
        # Binutils
        "ar", "as", "c++filt", "ld", "nm", "libtool", "objcopy", "ranlib", "readelf", "strip",

        # Misc. utilities
        "patchelf",
    ]

    if Sys.isapple(platform)
        append!(default_tools, ("dsymutil", "lipo", "otool", "install_name_tool"))
    elseif Sys.iswindows(platform)
        append!(default_tools, ("dlltool", "windres", "winmc"))
    end

    if :c in compilers
        append!(default_tools, ("cc", "c++", "cpp", "f77", "gfortran", "gcc", "clang", "g++", "clang++", "objc"))
    end
    if :rust in compilers
        append!(default_tools, ("rustc","rustup","cargo"))
    end
    if :go in compilers
        append!(default_tools, ("go",))
    end
    # Create symlinks for default compiler invocations, invoke target toolchain
    for tool in default_tools
        symlink("$(target)-$(tool)", joinpath(bin_path, triplet(platform), tool))
    end
end

# Translation mappers for our target names to cargo-compatible ones
map_rust_arch(p::AbstractPlatform) = replace(arch(p), "armv7l" => "armv7")
function map_rust_target(p::AbstractPlatform)
    if Sys.isapple(p)
        return "$(map_rust_arch(p))-apple-darwin"
    elseif Sys.isfreebsd(p)
        return "$(map_rust_arch(p))-unknown-freebsd"
    elseif Sys.iswindows(p)
        return "$(map_rust_arch(p))-pc-windows-gnu"
    else
        libc_str = libc(p) == "glibc" ? "gnu" : libc(p)
        call_abi_str = something(call_abi(p), "")
        return "$(map_rust_arch(p))-unknown-linux-$(libc_str)$(call_abi_str)"
    end
end

"""
    platform_envs(platform::AbstractPlatform)

Given a `platform`, generate a `Dict` mapping representing all the environment
variables to be set within the build environment to force compiles toward the
defined target architecture.  Examples of things set are `PATH`, `CC`,
`RANLIB`, as well as nonstandard things like `target`.
"""
function platform_envs(platform::AbstractPlatform, src_name::AbstractString;
                       host_platform = default_host_platform,
                       bootstrap::Bool=!isempty(bootstrap_list),
                       verbose::Bool = false)
    global use_ccache

    # Convert platform to a triplet, but strip out the ABI parts
    target = aatriplet(platform)
    host_target = aatriplet(host_platform)

    # Prefix, libdir, etc...
    prefix = "/workspace/destdir"
    host_prefix = "/workspace/$(triplet(host_platform))/destdir"
    if Sys.iswindows(platform)
        libdir = "$(prefix)/bin"
    else
        libdir = "$(prefix)/lib"
    end
    if Sys.iswindows(host_platform)
        host_libdir = "$(host_prefix)/bin"
    else
        host_libdir = "$(host_prefix)/lib"
    end

    if get(stdout, :color, false)
        PS1 = string(
            raw"\[",
            Base.text_colors[:light_blue],
            raw"\]",
            "sandbox",
            raw"\[",
            Base.text_colors[:normal],
            raw"\]",
            ":",
            raw"\[",
            Base.text_colors[:yellow],
            raw"\]",
            raw"${PWD//$WORKSPACE/$\{WORKSPACE\}}",
            raw"\[",
            Base.text_colors[:normal],
            raw"\]",
            raw" \$ ",
        )
    else
        PS1 = raw"sandbox:${PWD//$WORKSPACE/$\{WORKSPACE\}} $ "
    end

    # Base mappings
    mapping = Dict(
        # Platform information (we save a `bb_target` because sometimes `target` gets
        # overwritten in `./configure`, and we want tools like `uname` to still see it)
        "bb_target" => target,
        "target" => target,
        # Also give people access to the full `-cxxYY` and `-libgfortranX` triplet if they really want it
        "bb_full_target" => triplet(platform),
        "rust_target" => map_rust_target(platform),
        "rust_host" => map_rust_target(host_platform),
        "nproc" => "$(get(ENV, "BINARYBUILDER_NPROC", Sys.CPU_THREADS))",
        "nbits" => string(nbits(platform)),
        "proc_family" => string(proc_family(platform)),
        "dlext" => platform_dlext(platform),
        "exeext" => platform_exeext(platform),
        "PATH" => "/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin",
        "MACHTYPE" => "x86_64-linux-musl",

        # Set location parameters
        "WORKSPACE" => "/workspace",
        "prefix" => prefix,
        "bindir" => "$(prefix)/bin",
        "libdir" => libdir,
        "includedir" => "$(prefix)/include",

        # Set variables for the host prefix
        "host_prefix" => host_prefix,
        "host_bindir" => "$(host_prefix)/bin",
        "host_libdir" => host_libdir,
        "host_includedir" => "$(host_prefix)/include",

        # Fancyness!
        "USER" => get(ENV, "USER", "julia"),
        # Docker filters out `PS1` so we route around it
        "HIDDEN_PS1" => PS1,
        "VERBOSE" => "$(verbose)",
        "V" => "$(verbose)",
        "HISTFILE"=>"/meta/.bash_history",
        "TERM" => "screen",
        "SRC_NAME" => src_name,
    )

    # If we're bootstrapping, that's it, quit out.
    if bootstrap
        return mapping
    end

    # Helper for generating the library include path for a target.  MacOS, as usual,
    # puts things in slightly different place.
    function target_lib_dir(p::AbstractPlatform)
        t = aatriplet(p)
        if Sys.isapple(p)
            return "/opt/$(t)/$(t)/lib:/opt/$(t)/lib"
        else
            return "/opt/$(t)/$(t)/lib64:/opt/$(t)/$(t)/lib"
        end
    end

    function GOARM(p::AbstractPlatform)
        # See https://github.com/golang/go/wiki/GoArm#supported-architectures
        if arch(p) == "armv6l"
            return "6"
        elseif arch(p) == "armv7l"
            return "7"
        else
            return ""
        end
    end

    function csl_paths(p::AbstractPlatform)
        libcs = if Sys.islinux(p) && proc_family(p) == "intel" && libc(p) == "musl"
            # We need to push musl directories before glibc ones
            ("musl", "glibc")
        else
            ("glibc", "musl")
        end

        archs = if Sys.islinux(p) && proc_family(p) == "intel" && arch(p) == "i686"
            # We need to push i686 directories before x86_64 ones
            ("i686", "x86_64")
        else
            ("x86", "i686_64")
        end

        return join(["/usr/lib/csl-$(libc)-$(arch)" for libc in libcs, arch in archs], ":")
    end

    merge!(mapping, Dict(
        "PATH" => join((
            # First things first, our compiler wrappers trump all
            "/opt/bin/$(triplet(platform))",
            # Allow users to use things like x86_64-linux-gnu here
            "/opt/$(target)/bin",
            # Also wrappers for the host
            "/opt/bin/$(triplet(host_platform))",
            "/opt/$(host_target)/bin",
            # Default alpine PATH
            mapping["PATH"],
            # Host tools, installed as `HostBuildDependency`
            mapping["host_bindir"],
            # Finally, dependency tools
            mapping["bindir"],
        ), ":"),

        "LD_LIBRARY_PATH" => join((
           # Start with a default path
           "/usr/local/lib64:/usr/local/lib:/usr/lib64:/usr/lib",
            # Add our loader directories
            "/lib64:/lib",
            # Add our CSL libraries for all architectures that can natively run within this environment
            csl_paths(host_platform),
            # Add our target/host-specific library directories for compiler support libraries
            target_lib_dir(host_platform),
            target_lib_dir(platform),
            # Finally, dependencies
            "$(prefix)/lib64:$(prefix)/lib",
        ), ":"),

        # Default mappings for some tools
        "CC" => "cc",
        "CXX" => "c++",
        "OBJC" => "objc",
        "FC" => "gfortran",
        "GO" => "go",
        "RUSTC" => "rustc",
        "CARGO" => "cargo",

        # Go stuff
        "GOCACHE" => "/workspace/.gocache",
        "GOPATH" => "/workspace/.gopath",
        "GOARM" => GOARM(platform),

        # Rust stuff
        "CARGO_BUILD_TARGET" => map_rust_target(platform),
        "CARGO_HOME" => "/opt/$(host_target)",
        "RUSTUP_HOME" => "/opt/$(host_target)",
        # TODO: we'll need a way to parameterize this toolchain number
        "RUSTUP_TOOLCHAIN" => "1.43.0-$(map_rust_target(host_platform))",

        # We conditionally add on some compiler flags; we'll cull empty ones at the end
        "USE_CCACHE" => "$(use_ccache)",
        "LLVM_TARGET" => target,
        "LLVM_HOST_TARGET" => host_target,

        # Let the user parameterize their scripts for toolchain locations
        "CMAKE_HOST_TOOLCHAIN" => "/opt/toolchains/$(triplet(host_platform))/host_$(host_target).cmake",
        "CMAKE_TARGET_TOOLCHAIN" => "/opt/toolchains/$(triplet(platform))/target_$(target).cmake",
        "MESON_HOST_TOOLCHAIN" => "/opt/toolchains/$(triplet(host_platform))/$(host_target).meson",
        "MESON_TARGET_TOOLCHAIN" => "/opt/toolchains/$(triplet(platform))/$(target).meson",

        # We should always be looking for packages already in the prefix
        "PKG_CONFIG_PATH" => "$(prefix)/lib/pkgconfig:$(prefix)/lib64/pkgconfig:$(prefix)/share/pkgconfig",
        "PKG_CONFIG_SYSROOT_DIR" => prefix,

        # ccache options
        "CCACHE_COMPILERCHECK" => "content",

        # Things to help us step closer to reproducible builds; eliminate timestamp
        # variability within our binaries.
        "SOURCE_DATE_EPOCH" => "0",
        "ZERO_AR_DATE" => "1",
    ))

    # If we're on macOS, we give a hint to things like `configure` that they should use this as the linker
    if Sys.isapple(platform)
        mapping["LD"] = "/opt/$(target)/bin/ld64.macos"
        mapping["MACOSX_DEPLOYMENT_TARGET"] = "10.8"
    end

    # There is no broad agreement on what host compilers should be called,
    # so we set all the environment variables that we've seen them called
    # and hope for the best.
    for host_map in (tool -> "HOST$(tool)", tool -> "$(tool)_FOR_BUILD", tool -> "BUILD_$(tool)", tool -> "$(tool)_BUILD")
        # Use full path to avoid collisions when the target is similar to the
        # host (e.g., `x86_64-linux-musl-cxx03` and `x86_64-linux-musl-cxx11`)
        host_bin_dir = "/opt/bin/$(triplet(host_platform))"

        # First, do the simple tools where it's just X => $(host_target)-x:
        for tool in ("AR", "AS", "LD", "LIPO", "NM", "RANLIB", "READELF", "OBJCOPY", "OBJDUMP", "STRIP")
            mapping[host_map(tool)] = "$(host_bin_dir)/$(host_target)-$(lowercase(tool))"
        end

        # Next, the more custom tool mappings
        for (env_name, tool) in (
            "CC" => "$(host_bin_dir)/$(host_target)-gcc",
            "CXX" => "$(host_bin_dir)/$(host_target)-g++",
            "DSYMUTIL" => "dsymutil",
            "FC" => "$(host_bin_dir)/$(host_target)-gfortran"
           )
            mapping[host_map(env_name)] = tool
        end
    end

    # There is no easy way to automatically determine the version of glibc for
    # all platforms, but some build systems want to know it.  Let's emulate with
    # the `GNU_LIBC_VERSION` environment variable what `getconf
    # GNU_LIBC_VERSION` would return, if it worked.
    if libc(platform) === "glibc"
        # This should be kept in sync with the version of glibc used in
        # https://github.com/JuliaPackaging/Yggdrasil/blob/master/0_RootFS/gcc_common.jl
        if arch(platform) in ("x86_64", "i686")
            mapping["GNU_LIBC_VERSION"] = "glibc 2.12.2"
        elseif arch(platform) in ("armv7l", "aarch64")
            mapping["GNU_LIBC_VERSION"] = "glibc 2.19"
        elseif arch(platform) === "powerpc64le"
            mapping["GNU_LIBC_VERSION"] = "glibc 2.17"
        end
    end

    return mapping
end

runner_override = ""
function preferred_runner()
    global runner_override
    if runner_override != ""
        if runner_override in ["userns", "privileged"]
            return UserNSRunner
        elseif runner_override in ["docker"]
            return DockerRunner
        end
    end

    @static if Sys.islinux()
        return UserNSRunner
    else
        return DockerRunner
    end
end

"""
    runshell(platform::AbstractPlatform = HostPlatform())

Launch an interactive shell session within the user namespace, with environment
setup to target the given `platform`.
"""
function runshell(platform::AbstractPlatform = HostPlatform(); kwargs...)
    runshell(preferred_runner(), platform; kwargs...)
end

function runshell(r::Runner, args...; kwargs...)
    run_interactive(r, `/bin/bash -l`, args...; kwargs...)
end

function runshell(::Type{R}, platform::AbstractPlatform = HostPlatform(); verbose::Bool=false,kwargs...) where {R <: Runner}
    return runshell(R(pwd(); cwd="/workspace/", platform=platform, verbose=verbose, kwargs...); verbose=verbose)
end


text:
using LinearAlgebra: issymmetric

pbar = Progress(12, 0.2, "interface")

@testset "Interface" begin

@testset "From SimpleGraph" begin
    for ((gs, info), E_VAL) in product(make_testgraphs(SimpleGraph), test_edgeval_types)
        gv = SimpleValueGraph(gs, E_VAL)
        @test eltype(gs) == eltype(gv)
        @test edgeval_type(gv) == E_VAL
        @test all(edges(gv)) do e
            val(e) == default_edgeval(E_VAL)
        end
        @test nv(gs) == nv(gv)
        @test ne(gs) == ne(gv)
        @test all(edges(gs)) do e
            has_edge(gv, src(e), dst(e))
        end
    end
end

next!(pbar)

@testset "From SimpleDiGraph" begin
    for ((gs, info), E_VAL) in product(make_testgraphs(SimpleDiGraph), test_edgeval_types)
        gv = SimpleValueDiGraph(gs, E_VAL)
        @test eltype(gs) == eltype(gv)
        @test edgeval_type(gv) == E_VAL
        @test all(edges(gv)) do e
            val(e) == default_edgeval(E_VAL)
        end
        @test nv(gs) == nv(gv)
        @test ne(gs) == ne(gv)
        @test all(edges(gs)) do e
            has_edge(gv, src(e), dst(e))
        end
    end
end

next!(pbar)

@testset "add_edge! undirected" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        m = 25
        gs = SimpleGraph{V}(n)
        gv = SimpleValueGraph{V, E_VAL}(n)
        for i = 1:m
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            add_edge!(gs, u, v)
            add_edge!(gv, u, v, w)
            @test ne(gs) == ne(gv)
            @test has_edge(gv, u, v)
            @test has_edge(gv, u, v, w)
            @test get_edgeval(gv, u, v) == w
        end
    end
end

next!(pbar)

@testset "add_edge! directed" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        m = 25
        gs = SimpleDiGraph{V}(n)
        gv = SimpleValueDiGraph{V, E_VAL}(n)
        for i = 1:m
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            add_edge!(gs, u, v)
            add_edge!(gv, u, v, w)
            @test ne(gs) == ne(gv)
            @test has_edge(gv, u, v)
            @test has_edge(gv, u, v, w)
            @test get_edgeval(gv, u, v) == w
        end
    end
end

next!(pbar)

@testset "rem_edge! undirected" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        k = 25
        gs = CompleteGraph(V(n))
        gv = SimpleValueGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            rem_edge!(gs, u, v)
            rem_edge!(gv, u, v)
            @test ne(gs) == ne(gv)
            @test !has_edge(gv, u, v)
            @test get_edgeval(gv, u, v) == nothing
        end
    end
end

next!(pbar)

@testset "rem_edge! directed" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        k = 25
        gs = CompleteDiGraph(V(n))
        gv = SimpleValueDiGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            rem_edge!(gs, u, v)
            rem_edge!(gv, u, v)
            @test ne(gs) == ne(gv)
            @test !has_edge(gv, u, v)
            @test get_edgeval(gv, u, v) == nothing
        end
    end
end

next!(pbar)

@testset "get_edgeval & set_edgeval! undirected" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        m = 6
        k = 10
        gs = erdos_renyi(V(5), 6)
        gv = SimpleValueGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            set_edgeval!(gv, u, v, w)
            @test get_edgeval(gv, u, v) == (has_edge(gs, u, v) ? w : nothing)
        end
    end
end

next!(pbar)

@testset "get_edgeval & set_edgeval! with key undirected" begin
    for (V, E_VAL) in product(test_vertex_types, filter(T -> T <: TupleOrNamedTuple, test_edgeval_types))
        n = 5
        m = 6
        k = 10
        gs = erdos_renyi(V(5), 6)
        gv = SimpleValueGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            key = rand(1:length(E_VAL.parameters))
            set_edgeval!(gv, u, v, key, w[key])
            @test get_edgeval(gv, u, v, key) == (has_edge(gs, u, v) ? w[key] : nothing)
        end
    end
end

next!(pbar)

@testset "get_edgeval & set_edgeval! directed" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        n = 5
        m = 6
        k = 10
        gs = erdos_renyi(V(5), 6, is_directed=true)
        gv = SimpleValueDiGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            set_edgeval!(gv, u, v, w)
            @test get_edgeval(gv, u, v) == (has_edge(gs, u, v) ? w : nothing)
        end
    end
end

next!(pbar)

@testset "get_edgeval & set_edgeval! with key directed" begin
    for (V, E_VAL) in product(test_vertex_types, filter(T -> T <: TupleOrNamedTuple, test_edgeval_types))
        n = 5
        m = 6
        k = 10
        gs = erdos_renyi(V(5), 6, is_directed=true)
        gv = SimpleValueDiGraph(gs, E_VAL)
        for i = 1:k
            u = rand(1:n)
            v = rand(1:n)
            w = rand_sample(E_VAL)
            key = rand(1:length(E_VAL.parameters))
            set_edgeval!(gv, u, v, key, w[key])
            @test get_edgeval(gv, u, v, key) == (has_edge(gs, u, v) ? w[key] : nothing)
        end
    end
end

next!(pbar)

@testset "add_vertex! undirected" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        g = SimpleValueGraph(CompleteGraph(V(4)), E_VAL)
        n = nv(g)
        m = ne(g)
        for i in 1:10
            add_vertex!(g)
            @test nv(g) == n + i
            @test ne(g) == m
        end
    end
end

next!(pbar)

@testset "add_vertex! directed" begin
    for (V, E_VAL) in product(test_vertex_types, test_edgeval_types)
        g = SimpleValueDiGraph(CompleteDiGraph(V(4)), E_VAL)
        n = nv(g)
        m = ne(g)
        for i in 1:10
            add_vertex!(g)
            @test nv(g) == n + i
            @test ne(g) == m
        end
    end
end

next!(pbar)

end # testset Interface

ext: jl


text:
"""
    multipathfinder(
        logp,
        ∇logp,
        θ₀s::AbstractVector{AbstractVector{<:Real}},
        ndraws::Int;
        kwargs...
    )

Filter samples from a mixture of multivariate normal distributions fit using `pathfinder`.

For `nruns=length(θ₀s)`, `nruns` parallel runs of pathfinder produce `nruns` multivariate
normal approximations ``q_k = q(\\phi | \\mu_k, \\Sigma_k)`` of the posterior. These are
combined to a mixture model ``q`` with uniform weights.

``q`` is augmented with the component index to generate random samples, that is, elements
``(k, \\phi)`` are drawn from the augmented mixture model
```math
\\tilde{q}(\\phi, k | \\mu, \\Sigma) = K^{-1} q(\\phi | \\mu_k, \\Sigma_k),
```
where ``k`` is a component index, and ``K=`` `nruns`. These draws are then resampled with
replacement. Discarding ``k`` from the samples would reproduce draws from ``q``.

If `importance=true`, then Pareto smoothed importance resampling is used, so that the
resulting draws better approximate draws from the target distribution ``p`` instead of
``q``.

# Arguments
- `logp`: a callable that computes the log-density of the target distribution.
- `∇logp`: a callable that computes the gradient of `logp`.
- `θ₀s`: vector of length `nruns` of initial points of length `dim` from which each
    single-path Pathfinder run will begin
- `ndraws`: number of approximate draws to return

# Keywords
- `ndraws_per_run::Int=5`: The number of draws to take for each component before resampling.
- `importance::Bool=true`: Perform Pareto smoothed importance resampling of draws.

# Returns
- `q::Distributions.MixtureModel`: Uniformly weighted mixture of ELBO-maximizing
    multivariate normal distributions
- `ϕ::AbstractMatrix{<:Real}`: approximate draws from target distribution with size
    `(dim, ndraws)`
- `component_inds::Vector{Int}`: Indices ``k`` of components in ``q`` from which each column
in `ϕ` was drawn.
"""
function multipathfinder(
    logp,
    ∇logp,
    θ₀s,
    ndraws;
    ndraws_per_run::Int=5,
    rng::Random.AbstractRNG=Random.default_rng(),
    importance::Bool=true,
    kwargs...,
)
    if ndraws > ndraws_per_run * length(θ₀s)
        @warn "More draws requested than total number of draws across replicas. Draws will not be unique."
    end

    # run pathfinder independently from each starting point
    # TODO: allow to be parallelized
    prev_solutions = []
    res = filter(!isnothing, map(θ₀s) do θ₀
        try
            return pathfinder(logp, ∇logp, θ₀, ndraws_per_run; prev_solutions=prev_solutions, rng=rng, kwargs...)
        catch err
            if err.msg == "Insufficient history in trajectory."
                return nothing
            else
                rethrow(err)
            end
        end
    end)
    qs = reduce(vcat, first.(res))
    ϕs = reduce(hcat, getindex.(res, 2))

    # draw samples from augmented mixture model
    inds = axes(ϕs, 2)
    sample_inds = if importance
        logqϕs = reduce(vcat, last.(res))
        log_ratios = map(((ϕ, logqϕ),) -> logp(ϕ) - logqϕ, zip(eachcol(ϕs), logqϕs))
        resample(rng, inds, log_ratios, ndraws)
    else
        resample(rng, inds, ndraws)
    end

    q = Distributions.MixtureModel(qs)
    ϕ = ϕs[:, sample_inds]

    # get component ids (k) of draws in ϕ
    component_ids = cld.(sample_inds, ndraws_per_run)

    return q, ϕ, component_ids
end


text:
function create(c::Client, ::Type{Message}, ch::DiscordChannel; kwargs...)
    return create_message(c, ch.id; compkwfix(; kwargs...)...)
end

function retrieve(c::Client, ::Type{Message}, ch::DiscordChannel, message::Integer)
    return get_channel_message(c, ch.id, message)
end
function retrieve(
    c::Client,
    ::Type{Message},
    ch::DiscordChannel;
    pinned::Bool=false,
    kwargs...,
)
    return if pinned
        get_pinned_messages(c, ch.id)
    else
        get_channel_messages(c, ch.id; kwargs...)
    end
end

update(c::Client, m::Message; kwargs...) = edit_message(c, m.channel_id, m.id; kwargs...)

function delete(c::Client, m::Message; pinned::Bool=false)
    return if pinned
        delete_pinned_channel_message(c, m.channel_id, m.id)
    else
        delete_message(c, m.channel_id, m.id)
    end
end
function delete(c::Client, ms::Vector{Message})
    return bulk_delete_messages(c, ms[1].channel_id; messages=map(m -> m.id, ms))
end


text:
#= type Time_Anomaly <: Time_Operator =#
#=   DefaultDim :: (Int64..,) =#
#= end =#

#= Anomally ( sst :: Array, Dim :: Tuple{Int64,Vararg{Int64}}=(3,4)) =
#broadcast( - , sst, mean(sst, Dim)) =#
#= export Anomally =#


text:
"""
Provides functions to access NFL data in a tabular format.
"""
module NFLTables

# include submodules
include("enumerations.jl")
include("schedules.jl")
include("NFLScrapRData.jl")
include("NFLFastR.jl")

# export the API
export POST, PRE, REG, SeasonPart, SuperBowl
# export submodules
export NFLFastR, NFLScrapR, Schedules

end # module


text:
# Join works on collections of collections (e.g. a table is a collection of
# rows).

innerjoin(left, right) = innerjoin(identity, identity, left, right)
innerjoin(lkey::Callable, rkey::Callable, left, right) = innerjoin(lkey, rkey, merge, left, right)
innerjoin(lkey::Callable, rkey::Callable, f::Callable, left, right) = innerjoin(lkey, rkey, f, isequal, left, right)

const ⨝ = innerjoin

"""
    innerjoin(lkey, rkey, f, comparison, left, right)

Performs a relational-style join operation between iterables `left` and `right`, returning
a collection of elements `f(l, r)` for which `comparison(lkey(l), rkey(r))` is `true` where
`l ∈ left`, `r ∈ right.`

# Example

```jldoctest
julia> innerjoin(iseven, iseven, tuple, ==, [1,2,3,4], [0,1,2])
6-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (2, 0)
 (2, 2)
 (3, 1)
 (4, 0)
 (4, 2)
```
"""
function innerjoin(lkey::Callable, rkey::Callable, f::Callable, comparison::Callable, left, right)
    # TODO Do this inference-free, like comprehensions...
    T = promote_op(f, eltype(left), eltype(right))
    out = empty(left, T)

    innerjoin!(out, lkey, rkey, f, comparison, left, right)
    return out
end

function innerjoin!(out, lkey::Callable, rkey::Callable, f::Callable, comparison::Callable, left, right)
    # The O(length(left)*length(right)) generic method when nothing about `comparison` is known
    for a ∈ left
        for b ∈ right
            if comparison(lkey(a), rkey(b))
                push!(out, f(a, b))
            end
        end
    end
    return out
end

function innerjoin!(out, lkey::Callable, rkey::Callable, f::Callable, ::typeof(isequal), left, right)
    # isequal heralds a hash-based approach, roughly O(length(left) * log(length(right)))

    if Base.IteratorSize(left) === Base.SizeUnknown() || Base.IteratorSize(right) === Base.SizeUnknown() || length(right) <= length(left)
        K = promote_op(rkey, eltype(right))
        V = eltype(right)
        dict = Dict{K,Vector{V}}() # maybe a different stategy if right is unique?
        for b ∈ right
            key = rkey(b)
            push!(get!(Vector{V}, dict, key), b)
        end

        @inbounds for a ∈ left
            key = lkey(a)
            dict_index = Base.ht_keyindex(dict, key)
            if dict_index > 0 # -1 if key not found
                for b ∈ dict.vals[dict_index]
                    push!(out, f(a, b))
                end
            end
        end
    else
        K = promote_op(lkey, eltype(left))
        V = eltype(left)
        dict = Dict{K,Vector{V}}() # maybe a different stategy if left is unique?
        for a ∈ left
            key = lkey(a)
            push!(get!(Vector{V}, dict, key), a)
        end

        @inbounds for b ∈ right
            key = rkey(b)
            dict_index = Base.ht_keyindex(dict, key)
            if dict_index > 0 # -1 if key not found
                for a ∈ dict.vals[dict_index]
                    push!(out, f(a, b))
                end
            end
        end
    end

    return out
end

# For arrays, the assumptions around the below are
#  * Accessing arrays via indices and mapviews should be similar in speed to the above
#  * We can specialize these methods on particular arrays - works well for TypedTables.Table
#    of AcceleratedArrays.AcceleratedArray

function innerjoin!(out, lkey::Callable, rkey::Callable, f::Callable, comparison::Callable, left::AbstractArray, right::AbstractArray)
    _innerjoin!(out, mapview(lkey, left), mapview(rkey, right), productview(f, left, right), comparison)
end

function innerjoin!(out, lkey::Callable, rkey::Callable, f::Callable, ::typeof(isequal), left::AbstractArray, right::AbstractArray)
    _innerjoin!(out, mapview(lkey, left), mapview(rkey, right), productview(f, left, right), isequal)
end

function _innerjoin!(out, l::AbstractArray, r::AbstractArray, v::AbstractArray, comparison::Callable)
    @boundscheck if (axes(l)..., axes(r)...) != axes(v)
        throw(DimensionMismatch("innerjoin arrays do not have matching dimensions"))
    end

    @inbounds for i_l in keys(l)
        for i_r in keys(r)
            if comparison(l[i_l], r[i_r])
                push!(out, v[Tuple(i_l)..., Tuple(i_r)...])
            end
        end
    end

    return out
end

function _innerjoin!(out, l::AbstractArray, r::AbstractArray, v::AbstractArray, ::typeof(isequal))
    @boundscheck if (axes(l)..., axes(r)...) != axes(v)
        throw(DimensionMismatch("innerjoin arrays do not have matching dimensions"))
    end

    if length(r) <= length(l)
        rkeys = keys(r)
        V = eltype(rkeys)
        dict = Dict{eltype(r), Vector{V}}()
        @inbounds for i_r ∈ rkeys
            push!(get!(Vector{V}, dict, r[i_r]), i_r)
        end

        @inbounds for i_l in keys(l)
            l_value = l[i_l]
            dict_index = Base.ht_keyindex(dict, l_value)
            if dict_index > 0 # -1 if key not found
                for i_r ∈ dict.vals[dict_index]
                    push!(out, v[_tuple(i_l)..., _tuple(i_r)...])
                end
            end
        end
    else
        lkeys = keys(l)
        V = eltype(lkeys)
        dict = Dict{eltype(l), Vector{V}}()
        @inbounds for i_l ∈ lkeys
            push!(get!(Vector{V}, dict, l[i_l]), i_l)
        end

        @inbounds for i_r in keys(r)
            r_value = r[i_r]
            dict_index = Base.ht_keyindex(dict, r_value)
            if dict_index > 0 # -1 if key not found
                for i_l ∈ dict.vals[dict_index]
                    push!(out, v[_tuple(i_l)..., _tuple(i_r)...])
                end
            end
        end
    end

    return out
end

_tuple(i::Integer) = (i,)
_tuple(i::Any) = Tuple(i)

# TODO more specialized methods for comparisons: ==, <, isless, etc - via sorting strategies

# TODO perhaps a better version would be:
# function innerjoin(left, right; lkey = identity, rkey = identity, f = tuple, comparison = isequal)
#     ...
# end


text:
#
using SpectralElements.Spectral

nr = 8
ns = 12

arr = rand(nr,ns)
u   = Field(arr)

@test eachindex(vec(arr)) == eachindex(u)

@inferred similar(u)
@inferred +u
@inferred -u
@inferred 2u
@inferred 2 .+ u
@inferred 2 .- u
@inferred u .+ 2
@inferred u .- 2


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

include("model_AccountSasParameters.jl")
include("model_ActiveDirectoryProperties.jl")
include("model_AzureFilesIdentityBasedAuthentication.jl")
include("model_BlobRestoreParameters.jl")
include("model_BlobRestoreRange.jl")
include("model_BlobRestoreStatus.jl")
include("model_CheckNameAvailabilityResult.jl")
include("model_CustomDomain.jl")
include("model_DateAfterCreation.jl")
include("model_DateAfterModification.jl")
include("model_Dimension.jl")
include("model_Encryption.jl")
include("model_EncryptionScope.jl")
include("model_EncryptionScopeKeyVaultProperties.jl")
include("model_EncryptionScopeListResult.jl")
include("model_EncryptionScopeProperties.jl")
include("model_EncryptionService.jl")
include("model_EncryptionServices.jl")
include("model_Endpoints.jl")
include("model_ErrorResponse.jl")
include("model_GeoReplicationStats.jl")
include("model_IPRule.jl")
include("model_Identity.jl")
include("model_KeyVaultProperties.jl")
include("model_ListAccountSasResponse.jl")
include("model_ListServiceSasResponse.jl")
include("model_ManagementPolicy.jl")
include("model_ManagementPolicyAction.jl")
include("model_ManagementPolicyBaseBlob.jl")
include("model_ManagementPolicyDefinition.jl")
include("model_ManagementPolicyFilter.jl")
include("model_ManagementPolicyProperties.jl")
include("model_ManagementPolicyRule.jl")
include("model_ManagementPolicySchema.jl")
include("model_ManagementPolicySnapShot.jl")
include("model_MetricSpecification.jl")
include("model_NetworkRuleSet.jl")
include("model_ObjectReplicationPolicies.jl")
include("model_ObjectReplicationPolicy.jl")
include("model_ObjectReplicationPolicyFilter.jl")
include("model_ObjectReplicationPolicyProperties.jl")
include("model_ObjectReplicationPolicyRule.jl")
include("model_Operation.jl")
include("model_OperationDisplay.jl")
include("model_OperationListResult.jl")
include("model_OperationProperties.jl")
include("model_PrivateEndpoint.jl")
include("model_PrivateEndpointConnection.jl")
include("model_PrivateEndpointConnectionListResult.jl")
include("model_PrivateEndpointConnectionProperties.jl")
include("model_PrivateEndpointConnectionProvisioningState.jl")
include("model_PrivateEndpointServiceConnectionStatus.jl")
include("model_PrivateLinkResource.jl")
include("model_PrivateLinkResourceListResult.jl")
include("model_PrivateLinkResourceProperties.jl")
include("model_PrivateLinkServiceConnectionState.jl")
include("model_Resource.jl")
include("model_Restriction.jl")
include("model_RoutingPreference.jl")
include("model_SKUCapability.jl")
include("model_ServiceSasParameters.jl")
include("model_ServiceSpecification.jl")
include("model_Sku.jl")
include("model_SkuInformation.jl")
include("model_SkuName.jl")
include("model_StorageAccount.jl")
include("model_StorageAccountCheckNameAvailabilityParameters.jl")
include("model_StorageAccountCreateParameters.jl")
include("model_StorageAccountInternetEndpoints.jl")
include("model_StorageAccountKey.jl")
include("model_StorageAccountListKeysResult.jl")
include("model_StorageAccountListResult.jl")
include("model_StorageAccountMicrosoftEndpoints.jl")
include("model_StorageAccountProperties.jl")
include("model_StorageAccountPropertiesCreateParameters.jl")
include("model_StorageAccountPropertiesUpdateParameters.jl")
include("model_StorageAccountRegenerateKeyParameters.jl")
include("model_StorageAccountUpdateParameters.jl")
include("model_StorageSkuListResult.jl")
include("model_TagFilter.jl")
include("model_Tier.jl")
include("model_TrackedResource.jl")
include("model_Usage.jl")
include("model_UsageListResult.jl")
include("model_UsageName.jl")
include("model_VirtualNetworkRule.jl")


text:
function get_eop_history(obj::UT1Scale)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_name(obj::UT1Scale)
    return jcall(obj, "getName", JString, ())
end

function get_utc_scale(obj::UT1Scale)
    return jcall(obj, "getUTCScale", UTCScale, ())
end

function offset_from_tai(obj::UT1Scale, arg0::AbsoluteDate)
    return jcall(obj, "offsetFromTAI", jdouble, (AbsoluteDate,), arg0)
end

function offset_from_tai(obj::UT1Scale, arg0::FieldAbsoluteDate)
    return jcall(obj, "offsetFromTAI", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function to_string(obj::UT1Scale)
    return jcall(obj, "toString", JString, ())
end


text:
type RTRefSpace{T} <: RefSpace{T,3} end

valuetype{T}(ref::RTRefSpace{T}, charttype) = SVector{3,Tuple{Pt{universedimension(charttype),T},T}}

function (ϕ::RTRefSpace)(mp)

    u, v = parametric(mp)
    j = jacobian(mp)

    tu = tangents(mp,1)
    tv = tangents(mp,2)

    d = 2/j

    return SVector((
        ((tu*(u-1) + tv*v    ) / j, d),
        ((tu*u     + tv*(v-1)) / j, d),
        ((tu*u     + tv*v    ) / j, d)
    ))
end

divergence(ref::RTRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el))

"""
    ntrace(refspace, element, localindex, face)

Compute the normal trace of all local shape functions on `elements` belonging to
`refspace` on `face`. This function returns a matrix expressing the traces of local
shape functions in `refspace` as linear combinations of functions in the local
trace space. Cf. `restrict`. `localindex` is the index of `face` in the enumeration
of faces of `elements`. In many special cases knowing this index allows for highly
optimised implementations.
"""
function ntrace(x::RTRefSpace, el, q, fc)
    t = zeros(scalartype(x),1,3)
    t[q] = 1 / volume(fc)
    return t
end

function restrict{T}(ϕ::RTRefSpace{T}, dom1, dom2)

    K = numfunctions(ϕ)
    D = dimension(dom1)

    @assert K == 3
    @assert D == 2
    @assert D == dimension(dom2)

    Q = zeros(T,K,K)
    for i in 1:K

        # find the center of edge i of dom2
        a = dom2.vertices[mod1(i+1,D+1)]
        b = dom2.vertices[mod1(i+2,D+1)]
        c = (a + b) / 2

        # find the outer binormal there
        t = b - a
        l = norm(t)
        n = dom2.normals[1]
        m = cross(t, n) / l

        u = carttobary(dom1, c)
        x = neighborhood(dom1, u)

        y = ϕ(x)

        for j in 1:K
            Q[j,i] = dot(y[j][1], m) * l
        end
    end

    return Q
end


text:
using Enzyme

function main()
    function speelpenning(y::AbstractVector{VT}, x::AbstractVector{VT}) where {VT}
        y[1] = reduce(*, x)
        return nothing
    end

    y = [0.0]
    n = 10
    x = [i/(1.0+i) for i in 1:n]
    speelpenning(y,x)

    dx = zeros(n)
    dy = [1.0]
    autodiff(speelpenning, Duplicated(y,dy), Duplicated(x,dx))
    y = [0.0]
    speelpenning(y,x)

    errg = 0.0
    for (i, v) in enumerate(x)
        errg += abs(dx[i]-y[1]/v)
    end

    return (y[1]-1/(1.0+n)), errg
end


text:
function condition_io(u, t, integrator)
    UnPack.@unpack TS, Stats = integrator.p
    TS.dt_io += TS.dt
    io_flag = false
    if TS.dt_io > Stats.frequency
        TS.dt_io = 0
        io_flag = true
    end
    return io_flag || t ≈ 0 || t ≈ TS.t_max
end

condition_every_iter(u, t, integrator) = true

function affect_io!(integrator)
    UnPack.@unpack edmf, aux, grid, io_nt, diagnostics, case, gm, TS, Stats, skip_io = integrator.p
    skip_io && return nothing

    state = TC.State(integrator.u, aux, integrator.du)

    param_set = TC.parameter_set(gm)
    # TODO: is this the best location to call diagnostics?
    compute_diagnostics!(edmf, gm, grid, state, diagnostics, case, TS)

    # TODO: remove `vars` hack that avoids
    # https://github.com/Alexander-Barth/NCDatasets.jl/issues/135
    # opening/closing files every step should be okay. #removeVarsHack
    # TurbulenceConvection.io(sim) # #removeVarsHack
    TC.write_simulation_time(Stats, TS.t) # #removeVarsHack

    TC.io(io_nt.aux, Stats, state)
    TC.io(io_nt.diagnostics, Stats, diagnostics)

    TC.io(gm, grid, state, Stats) # #removeVarsHack
    TC.io(case, grid, state, Stats) # #removeVarsHack
    TC.io(edmf, grid, state, Stats, TS, param_set) # #removeVarsHack

    ODE.u_modified!(integrator, false) # We're legitamately not mutating `u` (the state vector)
end

function affect_filter!(integrator)
    UnPack.@unpack edmf, grid, gm, aux, case, TS = integrator.p
    state = TC.State(integrator.u, aux, integrator.du)
    TC.affect_filter!(edmf, grid, state, gm, case, TS)

    # We're lying to OrdinaryDiffEq.jl, in order to avoid
    # paying for an additional `∑tendencies!` call, which is required
    # to support supplying a continuous representation of the
    # solution.
    ODE.u_modified!(integrator, false)
end

function adaptive_dt!(integrator)
    UnPack.@unpack edmf, TS, dt_min = integrator.p
    TS.dt = min(TS.dt_max, max(edmf.dt_max, dt_min))
    SciMLBase.set_proposed_dt!(integrator, TS.dt)
    ODE.u_modified!(integrator, false)
end

function dt_max!(integrator)
    UnPack.@unpack gm, grid, edmf, aux, TS = integrator.p
    state = TC.State(integrator.u, aux, integrator.du)
    Δz = grid.Δz
    CFL_limit = TS.cfl_limit
    N_up = TC.n_updrafts(edmf)

    dt_max = TS.dt_max # initialize dt_max

    aux_tc = TC.center_aux_turbconv(state)
    aux_up_f = TC.face_aux_updrafts(state)
    aux_en_f = TC.face_aux_environment(state)
    KM = aux_tc.KM
    KH = aux_tc.KH

    # helper to calculate the rain velocity
    # TODO: assuming gm.W = 0
    # TODO: verify translation
    term_vel_rain = aux_tc.term_vel_rain
    term_vel_snow = aux_tc.term_vel_snow

    @inbounds for k in TC.real_face_indices(grid)
        TC.is_surface_face(grid, k) && continue
        @inbounds for i in 1:N_up
            dt_max = min(dt_max, CFL_limit * Δz / (abs(aux_up_f[i].w[k]) + eps(Float32)))
        end
        dt_max = min(dt_max, CFL_limit * Δz / (abs(aux_en_f.w[k]) + eps(Float32)))
    end
    @inbounds for k in TC.real_center_indices(grid)
        vel_max = max(term_vel_rain[k], term_vel_snow[k])
        # Check terminal rain/snow velocity CFL
        dt_max = min(dt_max, CFL_limit * Δz / (vel_max + eps(Float32)))
        # Check diffusion CFL (i.e., Fourier number)
        dt_max = min(dt_max, CFL_limit * Δz^2 / (max(KH[k], KM[k]) + eps(Float32)))
    end
    edmf.dt_max = dt_max

    ODE.u_modified!(integrator, false)
end

function monitor_cfl!(integrator)
    UnPack.@unpack gm, grid, edmf, aux, TS = integrator.p
    state = TC.State(integrator.u, aux, integrator.du)
    Δz = grid.Δz
    Δt = TS.dt
    CFL_limit = TS.cfl_limit

    aux_tc = TC.center_aux_turbconv(state)

    # helper to calculate the rain velocity
    # TODO: assuming gm.W = 0
    # TODO: verify translation
    term_vel_rain = aux_tc.term_vel_rain
    term_vel_snow = aux_tc.term_vel_snow

    @inbounds for k in TC.real_center_indices(grid)
        # check stability criterion
        CFL_out_rain = Δt / Δz * term_vel_rain[k]
        CFL_out_snow = Δt / Δz * term_vel_snow[k]
        if TC.is_toa_center(grid, k)
            CFL_in_rain = 0.0
            CFL_in_snow = 0.0
        else
            CFL_in_rain = Δt / Δz * term_vel_rain[k + 1]
            CFL_in_snow = Δt / Δz * term_vel_snow[k + 1]
        end
        if max(CFL_in_rain, CFL_in_snow, CFL_out_rain, CFL_out_snow) > CFL_limit
            error("Time step is too large for rain fall velocity!")
        end
    end

    ODE.u_modified!(integrator, false)
end


text:
using GPR
using ConstrainedDynamics
using ConstrainedDynamicsVis
using Plots
using Random
using GaussianProcesses
using LineSearches
using Optim

include("utils/utils.jl")
include("utils/data/simulations.jl")
include("utils/data/datasets.jl")

g = 9.81
m = 1.
l = 1.
Is = 1/12*m*l^2 + m*(.5l)^2

function energy(θ, ω) where T
    Ekin = .5Is*ω^2
    Epot = m*g*.5l*(1-cos(θ))
    return Epot + Ekin
end

function explicitEuler(steps, Δt, θ, ω)
    dω(θ, ω) = return -1.5sin(θ)*g/l  # 0.5lmgsin(θ)/I
    dθ(θ, ω) = return ω
    θs, ωs = zeros(steps), zeros(steps)
    for i in 1:steps
        Δθ = dθ(θ, ω)*Δt
        Δω = dω(θ, ω)*Δt
        θ += Δθ
        ω += Δω
        θs[i], ωs[i] = θ, ω
    end
    return θs, ωs
end

function gpVariationalIntegration(steps, θ, nsamples)
    Δtsim = 0.0001
    offset = Int(0.01/Δtsim)
    traindf = DataFrame(sold = Vector{Vector{State}}(), scurr = Vector{Vector{State}}())
    s = simplependulum2D(200*offset, Δt=Δtsim, θstart=θ)[1]
    for _ in 1:nsamples
        i = rand(1:length(s.x[1])-offset)
        push!(traindf, [getStates(s, i), getStates(s, i+offset)])
    end

    mechanism = simplependulum2D(1, Δt=0.01)[2]  # Reset Δt to 0.01 in mechanism
    # Create train and testsets
    xtrain_old = reduce(hcat, [CState(x) for x in traindf.sold])

    xtrain_curr = [CState(x) for x in traindf.scurr]
    vωindices = [9, 10, 11]
    ytrain = [[cs[i] for cs in xtrain_curr] for i in vωindices]

    predictedstates = Vector{CState{Float64,1}}()
    params = Vector{Float64}()  # declare in outer scope
    open(joinpath(dirname(@__FILE__), "config", "config.json"),"r") do f
        params = JSON.parse(f)["P1_MAX16"]
    end
    gps = Vector{GPE}()
    cache = MDCache()
    for (id, yi) in enumerate(ytrain)
        kernel = SEArd(log.(params[2:end]), log(params[1]))
        mean = MeanZero()
        gp = GP(xtrain_old, yi, mean, kernel)
        GaussianProcesses.optimize!(gp, LBFGS(linesearch = BackTracking(order=2)), Optim.Options(time_limit=10.))
        push!(gps, gp)
    end

    oldstates = CState([0, .5sin(θ)l, -.5cos(θ)l, q2vec(UnitQuaternion(RotX(θ)))..., zeros(6)...])
    setstates!(mechanism, oldstates)
    getvω(μ) = return [SVector(0,μ[1:2]...)], [SVector(μ[3],0,0)]
    for _ in 1:steps
        obs = reshape(oldstates.state, :, 1)
        μ = [predict_y(gp, obs)[1][1] for gp in gps]
        vcurr, ωcurr = getvω(μ)
        projectv!(vcurr, ωcurr, mechanism)
        foreachactive(updatestate!, mechanism.bodies, mechanism.Δt)  # Now at xcurr, vcurr
        oldstates = CState(mechanism)
        push!(predictedstates, oldstates)
    end
    θs, ωs = zeros(steps), zeros(steps)
    for (i, cstate) in enumerate(predictedstates)
        q = UnitQuaternion(cstate[4], cstate[5:7])
        θs[i] = Rotations.rotation_angle(q)*sign(q.x)*sign(q.w)  # Signum for axis direction
        ωs[i] = cstate[11]
    end
    return θs, ωs
end

function visualizePendulum(θs)
    mechanism = simplependulum2D(1)[2]
    s = Storage{Float64}(length(θs), 1)
    for i in 1:length(θs)
        x = [0, .5sin(θs[i]), -.5cos(θs[i])]
        q = UnitQuaternion(RotX(θs[i]))
        overwritestorage(s, CState([x..., q2vec(q)..., zeros(6)...]), i)
    end
    ConstrainedDynamicsVis.visualize(mechanism, s; showframes = true, env = "editor")
end

trueenergy = m*g*.5l*(1-cos(π/2))

θeuler, ωeuler = explicitEuler(100*60*100, 0.01, π/2, 0)
θgp, ωgp = gpVariationalIntegration(100*60*100, π/2, 20)
plt = plot([1:length(θeuler)]./100, [[abs(energy(θ, ω)/trueenergy-1)*100 for (θ, ω) in zip(θeuler, ωeuler)], [abs(energy(θ, ω)/trueenergy-1)*100 for (θ, ω) in zip(θgp, ωgp)]],
           title="Energy change comparison", label = ["Explicit Euler" "Implicit GP 20 samples"], xlabel="Time in seconds", ylabel = "Relative % error of total energy")
savefig(plt, "energyComparison100m")

θeuler, ωeuler = explicitEuler(100*10, 0.01, π/2, 0)
θgp, ωgp = gpVariationalIntegration(100*10, π/2, 20)
θgp2, ωgp2 = gpVariationalIntegration(100*10, π/2, 10)
plt = plot([1:length(θeuler)]./100, [[abs(energy(θ, ω)/trueenergy-1)*100 for (θ, ω) in zip(θeuler, ωeuler)], [abs(energy(θ, ω)/trueenergy-1)*100 for (θ, ω) in zip(θgp, ωgp)], [abs(energy(θ, ω)/trueenergy-1)*100 for (θ, ω) in zip(θgp2, ωgp2)]],
            title="Energy change comparison", label = ["Explicit Euler" "Implicit GP 20 samples" "Implicit GP 10 samples"], xlabel="Time in seconds", ylabel = "Relative % error of total energy")
savefig(plt, "energyComparison10s")


text:
using CSV, GLM, Plots, Random; pyplot()
Random.seed!(0)

df = CSV.read("../data/weightHeight.csv", copycols = true)
n = size(df)[1]
df[shuffle(1:n),:] = df
df[[10,40,60,130,140,175,190,200],:Sex] .= "O1"
df[[9,44,63,132,138,172,192,199],:Sex] .= "O2"

model = lm(@formula(Height ~ Weight + Sex), df,
        contrasts=Dict(:Sex=>DummyCoding(base="F",levels=["M","O1","O2","F"])))
b0, b1, b2, b3, b4  = coef(model)
pred(weight,sex) = b0+b1*weight+b2*(sex=="M")+b3*(sex=="O1")+b3*(sex=="O2")
println(model)

males = df[df.Sex .== "M",:]
females = df[df.Sex .== "F",:]
other = df[(df.Sex .!= "M") .& (df.Sex .!= "F"),:]

xlim = [minimum(df.Weight), maximum(df.Weight)]
scatter(males.Weight, males.Height, c=:blue, msw=0, label="Males")
plot!(xlim, pred.(xlim,"M"), c=:blue, label="Male model")

scatter!(females.Weight, females.Height, c=:red, msw=0, label="Females")
plot!(xlim, pred.(xlim,"F"),
	c=:red, label="Female model", xlims=(xlim),
	xlabel="Weight (kg)", ylabel="Height (cm)", legend=:topleft)

scatter!(other.Weight, other.Height, c=:green, msw=0, label="Other")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
## Cartpole / Inverted Pendulum
#TODO test
import TrajectoryOptimization
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics/urdf")
urdf_cartpole = joinpath(urdf_folder, "cartpole.urdf")

model_urdf = Model(urdf_cartpole,[1.;0.]) # underactuated, only control of slider

function cartpole_dynamics!(Xdot, X, U)
    mc = 2;   # mass of the cart in kg (10)
    mp = 0.5;    # mass of the pole (point mass at the end) in kg
    l = 0.5;   # length of the pole in m
    g = 9.81;  # gravity m/s^2

    q = X[1:2];
    qd = X[3:4];

    if isfinite(q[2])
        s = sin(q[2])
        c = cos(q[2])
    else
        s = Inf
        c = Inf
    end

    H = [mc+mp mp*l*c; mp*l*c mp*l^2];
    C = [0 -mp*qd[2]*l*s; 0 0];
    G = [0; mp*g*l*s];
    B = [1; 0];

    qdd = -H\(C*qd + G - B*U');

    Xdot[1:2] = qd
    Xdot[3:4] = qdd
    return nothing
end

function cartpole_dynamics_udp!(Xdot, X, U)
    mc = 10.0;   # mass of the cart in kg (10)
    mp = 1.0;    # mass of the pole (point mass at the end) in kg
    l = 0.5;   # length of the pole in m
    g = 9.81;  # gravity m/s^2

    q = X[1:2];
    qd = X[3:4];

    if isfinite(q[2])
        s = sin(q[2])
        c = cos(q[2])
    else
        s = Inf
        c = Inf
    end

    H = [mc+mp mp*l*c; mp*l*c mp*l^2];
    C = [0 -mp*qd[2]*l*s; 0 0];
    G = [0; mp*g*l*s];
    B = [1; 0];

    qdd = -H\(C*qd + G - B*U');

    Xdot[1:2] = qd
    Xdot[3:4] = qdd
    return nothing
end

n,m = 4,1
model_analytical = Model(cartpole_dynamics!,n,m)
model_udp = Model(cartpole_dynamics_udp!,n,m)

# initial and goal states
x0 = [0.;0.;0.;0.]
xf = [0.;pi;0.;0.]

x0_flipped = [0.;pi;0.;0.]
xf_flipped = [0.;0.;0.;0.]

# costs
Q = 0.01*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = 0.01*Diagonal(I,m)

# simulation
tf = 5.0
dt = 0.1

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_uncon_flipped = LQRObjective(Q, R, Qf, tf, x0_flipped, xf_flipped)

cartpole = [model_urdf, obj_uncon_flipped]
cartpole_analytical = [model_analytical, obj_uncon]
cartpole_udp = [model_udp, obj_uncon]
cartpole_mechanism = MechanismState(parse_urdf(Float64,urdf_cartpole))


text:
module LazyHelp

using Lazy

export @take
macro take(amt, l)
	:(take($(esc(amt)), $(esc(l))))
end

export @takewhile
macro takewhile(cond, l)
	:(takewhile($(esc(cond)), $(esc(l))))
end

export @map
macro map(f, l)
	:(map($(esc(f)), $(esc(l))))
end

export @length
macro length(l)
	:(length($(esc(l))))
end

export @<|
macro <|(f, args...)
	:(($(esc(f)))(($(esc(args)))...))
end

export @filter
macro filter(cond, l)
	:(filter($(esc(cond)), $(esc(l))))
end

export @head
macro head(l)
	:(first($(esc(l))))
end

export @tail
macro tail(l)
	:(tail($(esc(l))))
end

function fib()
	function rec(n::Int)
		(n == 1 || n == 2) && return 1
		return memo[n - 1] + memo[n - 2]
	end
	memo = @map rec Lazy.range(1)
	memo
end

include("lazyarray.jl")

end


text:
@testset "edgevals_type($G{$V, $E_VALS})" for
    G in (EdgeValGraph, EdgeValOutDiGraph, EdgeValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL

    g = G{V, E_VALS}(0)
    @assert g isa G{V, E_VALS}

    @test edgevals_type(g) == E_VALS
    @test edgevals_type(typeof(g)) == E_VALS
    @test edgevals_type(G{V, E_VALS}) == E_VALS

end

@testset "edgetype($G{$G, $E_VALS})" for
    G in (EdgeValGraph, EdgeValOutDiGraph, EdgeValDiGraph),
    V in TEST_VERTEX_TYPES_SMALL,
    E_VALS in TEST_EDGEVAL_TYPES_SMALL

    g = G{V, E_VALS}(0)
    @assert g isa G{V, E_VALS}

    ET_should_be = (is_directed(g) ? ValDiEdge : ValEdge){V, E_VALS}

    @test edgetype(g) == ET_should_be
    @test edgetype(typeof(g)) == ET_should_be
    @test edgetype(G{V, E_VALS}) == ET_should_be

end


text:
@cell arc(Point2f(0), 10f0, 0f0, pi, linewidth=20)

@cell mesh(Recti(0, 0, 200, 200))

@cell poly(Recti(0, 0, 200, 200), strokewidth=20, strokecolor=:red, color=(:black, 0.4))

@cell begin
    scene = poly([Rect(0, 0, 20, 20)])
    scatter!(Rect(0, 0, 20, 20), color=:red, markersize=20)
    current_figure()
end

@cell begin
    lines(Rect(0, 0, 1, 1), linewidth=4)
    scatter!([Point2f(0.5, 0.5)], markersize=1, markerspace=SceneSpace, marker='I')
    current_figure()
end

@cell lines(RNG.rand(10), RNG.rand(10), color=RNG.rand(10), linewidth=10)
@cell lines(RNG.rand(10), RNG.rand(10), color=RNG.rand(RGBAf, 10), linewidth=10)
@cell scatter(0..1, RNG.rand(10), markersize=RNG.rand(10) .* 20)
@cell scatter(LinRange(0, 1, 10), RNG.rand(10))
@cell scatter(RNG.rand(10), LinRange(0, 1, 10))

@cell begin
    angles = range(0, stop=2pi, length=20)
    pos = Point2f.(sin.(angles), cos.(angles))
    scatter(pos, markersize=0.2, markerspace=SceneSpace, rotations=-angles, marker='▲', axis=(;aspect = DataAspect()))
    scatter!(pos, markersize=10, color=:red)
    current_figure()
end

@cell heatmap(RNG.rand(50, 50), colormap=:RdBu, alpha=0.2)

@cell contour(RNG.rand(10, 100))
@cell contour(RNG.rand(100, 10))
@cell contour(RNG.randn(100, 90), levels=3)

@cell contour(RNG.randn(100, 90), levels=[0.1, 0.5, 0.8])
@cell contour(RNG.randn(100, 90), levels=[0.1, 0.5, 0.8], color=:black)
@cell contour(RNG.randn(33, 30), levels=[0.1, 0.5, 0.9], color=[:black, :green, (:blue, 0.4)], linewidth=2)
@cell contour(RNG.randn(33, 30), levels=[0.1, 0.5, 0.9], colormap=:Spectral)
@cell contour(
    RNG.rand(33, 30) .* 6 .- 3, levels=[-2.5, 0.4, 0.5, 0.6, 2.5],
    colormap=[(:black, 0.2), :red, :blue, :green, (:black, 0.2)],
    colorrange=(0.2, 0.8)
)

@cell lines(Circle(Point2f(0), Float32(1)))

@cell begin
    v(x::Point2{T}) where T = Point2{T}(x[2], 4 * x[1])
    streamplot(v, -2..2, -2..2)
end

@cell lines(-1..1, x -> x^2)
@cell scatter(-1..1, x -> x^2)

@cell begin
    r = range(-3pi, stop=3pi, length=100)
    fig, ax, vplot = Makie.volume(r, r, r, (x, y, z) -> cos(x) + sin(y) + cos(z), algorithm=:iso, isorange=0.1f0, show_axis=false)
    v2 = volume!(ax, r, r, r, (x, y, z) -> cos(x) + sin(y) + cos(z), algorithm=:mip,
                 transformation=(translation=Vec3f(6pi, 0, 0),))
    fig
end

@cell meshscatter(RNG.rand(10), RNG.rand(10), RNG.rand(10), color=RNG.rand(10))
@cell meshscatter(RNG.rand(10), RNG.rand(10), RNG.rand(10), color=RNG.rand(RGBAf, 10))

@cell begin
    s1 = uv_mesh(Sphere(Point3f(0), 1f0))
    mesh(uv_mesh(Sphere(Point3f(0), 1f0)), color=RNG.rand(50, 50))
    # ugh, bug In GeometryTypes for UVs of non unit spheres.
    s2 = uv_mesh(Sphere(Point3f(0), 1f0))
    s2.position .= s2.position .+ (Point3f(0, 2, 0),)
    mesh!(s2, color=RNG.rand(RGBAf, 50, 50))
    current_figure()
end

@cell "Unequal x and y sizes in surface" begin
    #
    NL = 15
    NR = 31
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    lspace = range(-10, stop=10, length=NL)
    rspace = range(-10, stop=10, length=NR)

    z = Float32[xy_data(x, y) for x in lspace, y in rspace]
    l = range(0, stop=3, length=NL)
    r = range(0, stop=3, length=NR)
    surface(
        l, r, z,
        colormap=:Spectral
    )
end

@cell "Matrices of data in surfaces" begin
    NL = 30
    NR = 31
    function xy_data(x, y)
        r = sqrt(x^2 + y^2)
        r == 0.0 ? 1f0 : (sin(r) / r)
    end
    lspace = range(-10, stop=10, length=NL)
    rspace = range(-10, stop=10, length=NR)

    z = Float32[xy_data(x, y) for x in lspace, y in rspace]
    l = range(0, stop=3, length=NL)
    r = range(0, stop=3, length=NR)
    surface(
        [l for l in l, r in r], [r for l in l, r in r], z,
        colormap=:Spectral
    )
end

@cell "heatmaps & surface" begin
    data =
        hcat(LinRange(2, 3, 4), LinRange(2, 2.5, 4), LinRange(2.5, 3, 4), [1, NaN, NaN, 5])

    fig = Figure()
    heatmap(
        fig[1, 1],
        data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
    )
    surface(
        fig[1, 2],
        zeros(size(data)),
        color = data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
        shading = false,
    )
    surface!(
        Axis(fig[2, 2]),
        data,
        colorrange = (2, 3),
        highclip = :red,
        lowclip = :black,
        nan_color = (:green, 0.5),
        shading = false,
    )
    fig
end

@cell "log10 heatmap" begin
    heatmap(RNG.rand(10, 5), axis = (yscale = log10, xscale=log10))
end

@cell "reverse range heatmap" begin
    x = [1 0
         2 3]
    heatmap(1:2, 1:-1:0, x)
end


text:
struct Dmnorm
    mu
    Uinv
    det
end

"""
    dmnorm(X = nothing; mu = nothing, S = nothing)
Compute the normal probability density of
multivariate observations.
* `X` : X-data used to estimate the mean and
    the covariance matrix of the population.
    If `nothing`, `mu` and `S` must be provided.
* `mu` : Mean vector of the normal distribution.
    If `nothing`, `mu` is computed by the column-means of `X`.
* `S` : Covariance matrix of the normal distribution.
    If `nothing`, `S` is computed by cov(`X`).
"""
function dmnorm(X = nothing; mu = nothing, S = nothing)
    isnothing(S) ? zS = nothing : zS = copy(S)
    dmnorm!(X; mu = mu, S = zS)
end

function dmnorm!(X = nothing; mu = nothing, S = nothing)
    !isnothing(X) ? X = ensure_mat(X) : nothing
    if isnothing(mu)
        mu = vec(mean(X, dims = 1))
    end
    if isnothing(S)
        S = cov(X)
    end
    U = cholesky!(Hermitian(S)).U # This modifies S only if S is provided
    zd = det(U)^2
    zd == 0 ? zd = 1e-20 : nothing
    Uinv = LinearAlgebra.inv!(U)
    #cholesky!(S)
    #U = sqrt(diag(diag(S), nrow = p))
    #Uinv = solve(diag(diag(S), nrow = p))
    Dmnorm(mu, Uinv, zd)
end

function predict(object::Dmnorm, X)
    X = ensure_mat(X)
    p = size(X, 2)
    mu = reshape(object.mu, 1, length(object.mu))
    d = mahsqchol(X, mu, object.Uinv)
    ds = (2 * pi)^(-p / 2) * (1 / sqrt(object.det)) * exp.(-.5 * d)
    (pred = ds,)
end


text:
# # Chebyshev center
# Boyd & Vandenberghe, "Convex Optimization"
# Joëlle Skaf - 08/16/05
#
# Adapted for Convex.jl by Karanveer Mohan and David Zeng - 26/05/14
#
# The goal is to find the largest Euclidean ball (i.e. its center and
# radius) that lies in a polyhedron described by affine inequalites in this
# fashion: $P = \{x : a_i'*x \leq b_i, i=1,\ldots,m \}$ where $x \in \mathbb{R}^2$.

using Convex, LinearAlgebra, SCS

# Generate the input data
a1 = [ 2;  1];
a2 = [ 2; -1];
a3 = [-1;  2];
a4 = [-1; -2];
b = ones(4, 1);

# Create and solve the model
r = Variable(1)
x_c = Variable(2)
p = maximize(r)
p.constraints += a1' * x_c + r * norm(a1, 2) <= b[1];
p.constraints += a2' * x_c + r * norm(a2, 2) <= b[2];
p.constraints += a3' * x_c + r * norm(a3, 2) <= b[3];
p.constraints += a4' * x_c + r * norm(a4, 2) <= b[4];
solve!(p, () -> SCS.Optimizer(verbose=0))
p.optval

# Generate the figure
x = range(-1.5, stop=1.5, length=100);
theta = 0:pi/100:2*pi;
using Plots
plot(x, x -> -x * a1[1] / a1[2] + b[1] / a1[2])
plot!(x, x -> -x * a2[1]/ a2[2] + b[2] / a2[2])
plot!(x, x -> -x * a3[1]/ a3[2] + b[3] / a3[2])
plot!(x, x -> -x * a4[1]/ a4[2] + b[4] / a4[2])
plot!(evaluate(x_c)[1] .+ evaluate(r) * cos.(theta), evaluate(x_c)[2] .+ evaluate(r) * sin.(theta), linewidth = 2)
plot!(title ="Largest Euclidean ball lying in a 2D polyhedron", legend = nothing)


text:
module IDS

using Spec
import ..Util: Counter, increment!, Box, reset!
using DataStructures: LinkedList, cons, nil, list, head, tail
export append, base, combine, increment!, increment,
       Paired, pair
export combine, append, base, increment, firstelem
export withscope, scope, scope!
export ID, @uid, uid, toid

include("idinterface.jl")
include("idtypes.jl")
include("id.jl")

end


text: println("Hello wolrd!")