text:
export build_tarballs, autobuild, print_buildjl, product_hashes_from_github_release, build
import GitHub: gh_get_json, DEFAULT_API
import SHA: sha256

"""
    build_tarballs(ARGS, src_name, src_version, sources, script, platforms,
                   products, dependencies; kwargs...)

This should be the top-level function called from a `build_tarballs.jl` file.
It takes in the information baked into a `build_tarballs.jl` file such as the
`sources` to download, the `products` to build, etc... and will automatically
download, build and package the tarballs, generating a `build.jl` file when
appropriate.  Note that `ARGS` should be the top-level Julia `ARGS` command-
line arguments object.  This function does some rudimentary parsing of the
`ARGS`, call it with `--help` in the `ARGS` to see what it can do.
"""
function build_tarballs(ARGS, src_name, src_version, sources, script,
                        platforms, products, dependencies; kwargs...)
    # See if someone has passed in `--help`, and if so, give them the
    # assistance they so clearly long for
    if "--help" in ARGS
        println(strip("""
        Usage: build_tarballs.jl [target1,target2,...] [--only-buildjl]
                                 [--verbose] [--help]

        Options:
            targets         By default `build_tarballs.jl` will build a tarball
                            for every target within the `platforms` variable.
                            To override this, pass in a list of comma-separated
                            target triplets for each target to be built.  Note
                            that this can be used to build for platforms that
                            are not listed in the 'default list' of platforms
                            in the build_tarballs.jl script.

            --part=n/m      For breaking up long builds into multiple jobs,
                            divides the `platforms` list into `m` roughly
                            equal parts and then only builds part `n`
                            (`1 ≤ n ≤ m`).  (Does not produce `build.jl`,
                            which must be generated by a later
                            `--only-buildjl` stage.)

            --verbose       This streams compiler output to stdout during the
                            build which can be very helpful for finding bugs.
                            Note that it is colorized if you pass the
                            --color=yes option to julia, see examples below.

            --debug         This causes a failed build to drop into an
                            interactive bash shell for debugging purposes.

            --only-buildjl  This disables building of any tarballs, and merely
                            reconstructs a `build.jl` file from a github
                            release.  This is mostly useful as a later stage in
                            a travis/github releases autodeployment setup.

            --help          Print out this message.

        Examples:
            julia --color=yes build_tarballs.jl --verbose
                This builds all tarballs, with colorized output.

            julia build_tarballs.jl x86_64-linux-gnu,i686-linux-gnu
                This builds two tarballs for the two platforms given, with a
                minimum of output messages.
        """))
        return nothing
    end

    function check_flag(flag)
        flag_present = flag in ARGS
        ARGS = filter!(x -> x != flag, ARGS)
        return flag_present
    end

    # This sets whether we should build verbosely or not
    verbose = check_flag("--verbose")

    # This sets whether we drop into a debug shell on failure or not
    debug = check_flag("--debug")

    # This flag skips actually building and instead attempts to reconstruct a
    # build.jl from a GitHub release page.  Use this to automatically deploy a
    # build.jl file even when sharding targets across multiple CI builds.
    only_buildjl = check_flag("--only-buildjl")

    # --part=n/m builds only part n out of m divisions
    # of the platforms list.
    i = findlast(x -> startswith(x, "--part="), ARGS)
    should_override_platforms = i !== nothing
    if should_override_platforms
        if i != findfirst(x -> startswith(x, "--part="), ARGS)
            error("multiple --part arguments are not allowed")
        end
        p = parse.(Int, split(ARGS[i][8:end], '/'))
        (length(p) == 2 && p[2] > 0 && 1 ≤ p[1] ≤ p[2]) || error("invalid argument ", ARGS[i])
        n = (length(platforms) + p[2]-1) ÷ p[2]
        platforms = platforms[n*(p[1]-1)+1:min(end,n*p[1])]
        deleteat!(ARGS, i)
    end

    # If the user passed in a platform (or a few, comma-separated) on the
    # command-line, use that instead of our default platforms
    if length(ARGS) > 0
        should_override_platforms = true
        platforms = platform_key_abi.(split(ARGS[1], ","))
    end

    # If we're running on CI (Travis, GitLab CI, etc...) and this is a
    # tagged release, automatically determine bin_path by building up a URL
    repo = get_repo_name()
    tag = get_tag_name()

    product_hashes = if !only_buildjl
        # If the user didn't just ask for a `build.jl`, go ahead and actually build
        @info("Building for $(join(triplet.(platforms), ", "))")

        # Build the given platforms using the given sources
        autobuild(pwd(), src_name, src_version, sources, script, platforms,
                         products, dependencies; verbose=verbose, debug=debug, kwargs...)
    else
        msg = strip("""
        Reconstructing product hashes from GitHub Release $(repo)/$(tag)
        """)
        @info(msg)

        # Reconstruct product_hashes from github
        product_hashes_from_github_release(repo, tag;
            product_filter="v$(src_version)",
            verbose=verbose
        )
    end

    # If we didn't override the default set of platforms OR we asked for only
    # a build.jl file, then write one out.  We don't write out when overriding
    # the default set of platforms because that is typically done either while
    # testing, or when we have sharded our tarball construction over multiple
    # invocations.
    if !should_override_platforms || only_buildjl
        # The location the binaries will be available from
        bin_path = "https://github.com/$(repo)/releases/download/$(tag)"

        # A dummy prefix to pass through products()
        dummy_products = products(Prefix(pwd()))
        print_buildjl(pwd(), src_name, src_version, dummy_products,
                      product_hashes, bin_path)

        if verbose
            @info("Writing out the following reconstructed build.jl:")
            print_buildjl(Base.stdout, dummy_products, product_hashes, bin_path)
        end
    end

    return product_hashes
end

# Helper function to get things from ENV, returning `nothing`
# if they either don't exist or are empty
function get_ENV(key)
    if !haskey(ENV, key)
        return nothing
    end

    if isempty(ENV[key])
        return nothing
    end

    return ENV[key]
end

function get_repo_name()
    # Helper function to synthesize repository slug from environment variables
    function get_gitlab_repo_name()
        owner = get_ENV("CI_REPO_OWNER")
        name = get_ENV("CI_REPO_NAME")
        if owner != nothing && name != nothing
            return "$(owner)/$(name)"
        end
        return nothing
    end

    # Helper function to guess repository slug from git remote URL
    function read_git_origin()
        try
            repo = LibGit2.GitRepo(".")
            url = LibGit2.url(LibGit2.get(LibGit2.GitRemote, repo, "origin"))
            owner = basename(dirname(url))
            if occursin(":", owner)
                owner = owner[findlast(isequal(':'), owner)+1:end]
            end
            name = basename(url)
            if endswith(name, ".git")
                name = name[1:end-4]
            end
            return "$(owner)/$(name)"
        catch
        end
        return nothing
    end

    return something(
        get_ENV("TRAVIS_REPO_SLUG"),
        get_gitlab_repo_name(),
        read_git_origin(),
        "<repo owner>/<repo name>",
    )
end

function get_tag_name()
    # Helper function to guess tag from current commit taggedness
    function read_git_tag()
        try
            repo = LibGit2.GitRepo(".")
            head_gitsha = LibGit2.GitHash(LibGit2.head(repo))
            for tag in LibGit2.tag_list(repo)
                tag_gitsha = LibGit2.GitHash(LibGit2.GitCommit(repo, tag))
                if head_gitsha == tag_gitsha
                    return tag
                end
            end
        catch
        end
        return nothing
    end

    return something(
        get_ENV("TRAVIS_TAG"),
        get_ENV("CI_COMMIT_TAG"),
        read_git_tag(),
        "<tag>",
    )
end

"""
    autobuild(dir::AbstractString, src_name::AbstractString,
              src_version::VersionNumber, sources::Vector,
              script::AbstractString, platforms::Vector,
              products::Function, dependencies::Vector;
              verbose::Bool = true, kwargs...)

Runs the boiler plate code to download, build, and package a source package
for a list of platforms.  `src_name` represents the name of the source package
being built (and will set the name of the built tarballs), `platforms` is a
list of platforms to build for, `sources` is a list of tuples giving
`(url, hash)` of all sources to download and unpack before building begins,
`script` is a string representing a `bash` script to run to build the desired
products, which are listed as `Product` objects within the vector returned by
the `products` function. `dependencies` gives a list of dependencies that
provide `build.jl` files that should be installed before building begins to
allow this build process to depend on the results of another build process.
Setting `debug` to `true` will cause a failed build to drop into an interactive
shell so that the build can be inspected easily.
"""
function autobuild(dir::AbstractString,
                   src_name::AbstractString,
                   src_version::VersionNumber,
                   sources::Vector,
                   script::AbstractString,
                   platforms::Vector,
                   products::Function,
                   dependencies::Vector;
                   verbose::Bool = true,
                   kwargs...)
    # If we're on CI and we're not verbose, schedule a task to output a "." every few seconds
    if (haskey(ENV, "TRAVIS") || haskey(ENV, "CI")) && !verbose
        run_travis_busytask = true
        travis_busytask = @async begin
            # Don't let Travis think we're asleep...
            @info("Brewing a pot of coffee for Travis...")
            while run_travis_busytask
                sleep(4)
                print(".")
            end
        end
    end

    # This is what we'll eventually return
    product_hashes = Dict()

    # If we end up packaging any local directories into tarballs, we'll store them here
    mktempdir() do tempdir
        # We must prepare our sources.  Download them, hash them, etc...
        sources = Any[s for s in sources]
        for idx in 1:length(sources)
            # If the given source is a local path that is a directory, package it up and insert it into our sources
            if typeof(sources[idx]) <: AbstractString
                if !isdir(sources[idx])
                    error("Sources must either be a pair (url => hash) or a local directory")
                end

                # Package up this directory and calculate its hash
                tarball_path = joinpath(tempdir, basename(sources[idx]) * ".tar.gz")
                package(sources[idx], tarball_path; verbose=verbose)
                tarball_hash = open(tarball_path, "r") do f
                    bytes2hex(sha256(f))
                end

                # Move it to a filename that has the hash as a part of it (to avoid name collisions)
                tarball_pathv = joinpath(tempdir, string(tarball_hash, "-", basename(sources[idx]), ".tar.gz"))
                mv(tarball_path, tarball_pathv)

                # Now that it's packaged, store this into sources[idx]
                sources[idx] = (tarball_pathv => tarball_hash)
            elseif typeof(sources[idx]) <: Pair
                src_url, src_hash = sources[idx]

                # If it's a .git url, clone it
                if endswith(src_url, ".git")
                    src_path = storage_dir("downloads", basename(src_url))

                    # If this git repository already exists, ensure that its origin remote actually matches
                    if isdir(src_path)
                        origin_url = LibGit2.with(LibGit2.GitRepo(src_path)) do repo
                            LibGit2.url(LibGit2.get(LibGit2.GitRemote, repo, "origin"))
                        end

                        # If the origin url doesn't match, wipe out this git repo.  We'd rather have a
                        # thrashed cache than an incorrect cache.
                        if origin_url != src_url
                            rm(src_path; recursive=true, force=true)
                        end
                    end

                    if isdir(src_path)
                        # If we didn't just mercilessly obliterate the cached git repo, use it!
                        LibGit2.with(LibGit2.GitRepo(src_path)) do repo
                            LibGit2.fetch(repo)
                        end
                    else
                        # If there is no src_path yet, clone it down.
                        repo = LibGit2.clone(src_url, src_path; isbare=true)
                    end
                else
                    if isfile(src_url)
                        # Immediately abspath() a src_url so we don't lose track of
                        # sources given to us with a relative path
                        src_path = abspath(src_url)

                        # And if this is a locally-sourced tarball, just verify
                        verify(src_path, src_hash; verbose=verbose)
                    else
                        # Otherwise, download and verify
                        src_path = storage_dir("downloads", basename(src_url))
                        download_verify(src_url, src_hash, src_path; verbose=verbose)
                    end
                end

                # Now that it's downloaded, store this into sources[idx]
                sources[idx] = (src_path => src_hash)
            else
                error("Sources must be either a `URL => hash` pair, or a path to a local directory")
            end
        end

        # Our build products will go into ./products
        out_path = joinpath(dir, "products")
        try mkpath(out_path) catch; end

        for platform in platforms
            target = triplet(platform)

            # We build in a platform-specific directory
            build_path = joinpath(pwd(), "build", target)
            try mkpath(build_path) catch; end

            src_paths, src_hashes = collect(zip(sources...))

            # Convert from tuples to arrays, if need be
            src_paths = collect(src_paths)
            src_hashes = collect(src_hashes)
            prefix, ur = setup_workspace(
                build_path,
                src_paths,
                src_hashes,
                dependencies,
                platform;
                verbose=verbose,
                downloads_dir=storage_dir("downloads"),
            )

            # Don't keep the downloads directory around
            rm(joinpath(prefix, "downloads"); force=true, recursive=true)

            # Collect dependency manifests so that our auditing doesn't touch these files that
            # were installed by dependencies
            manifest_dir = joinpath(prefix, "manifests")
            dep_manifests = if isdir(manifest_dir)
               [joinpath(prefix, "manifests", f) for f in readdir(manifest_dir)]
            else
                String[]
            end

            build(ur, src_name, products(prefix), script, platform, prefix;
                  verbose=verbose, ignore_manifests=dep_manifests, kwargs...)

            # Remove the files of any dependencies
            for dependency in dependencies
                dep_script = script_for_dep(dependency, prefix.path)[1]
                m = Module(:__anon__)
                Core.eval(m, quote
                    using BinaryProvider
                    # Override BinaryProvider functionality so that it doesn't actually install anything
                    platform_key() = $platform
                    platform_key_abi() = $platform
                    function write_deps_file(args...; kwargs...); end
                    function install(args...; kwargs...); end

                    # Include build.jl file to extract download_info
                    ARGS = [$(prefix.path)]
                end)
                include_string(m, dep_script)
                Core.eval(m, quote
                    # Grab the information we need in order to extract a manifest, then uninstall it
                    url, hash = choose_download(download_info, platform_key_abi())
                    manifest_path = BinaryProvider.manifest_from_url(url; prefix=prefix)
                    BinaryProvider.uninstall(manifest_path; verbose=$verbose)
                end)
            end

            # Once we're built up, go ahead and package this prefix out
            tarball_path, tarball_hash = package(
                prefix,
                joinpath(out_path, src_name),
                src_version;
                platform=platform,
                verbose=verbose,
                force=true,
            )
            product_hashes[target] = (basename(tarball_path), tarball_hash)

            # Destroy the workspace
            rm(dirname(prefix.path); recursive=true)

            # If the whole build_path is empty, then remove it too.  If it's not, it's probably
            # because some other build is doing something simultaneously with this target, and we
            # don't want to mess with their stuff.
            if isempty(readdir(build_path))
                rm(build_path; recursive=true)
            end
        end
    end

    if (haskey(ENV, "TRAVIS") || haskey(ENV, "CI")) && !verbose
        run_travis_busytask = false
        wait(travis_busytask)
        println()
    end

    # Return our product hashes
    return product_hashes
end

function build(runner::Runner, name::AbstractString,
               products::Vector{P}, script::AbstractString,
               platform::Platform, prefix::Prefix;
               verbose::Bool = false, force::Bool = false,
               autofix::Bool = true, ignore_audit_errors::Bool = true,
               skip_audit::Bool = false,
               ignore_manifests::Vector = [], debug::Bool = false) where {P <: Product}
    # First, look to see whether our products are satisfied or not
    if isempty(products)
        # If we've been given no products, always build and always say it's satisfied
        s = product -> true
        should_build = true
    else
        s = p -> satisfied(p; verbose=verbose, platform=platform, isolate=true)
        should_build = !all(s(p) for p in products)
    end

    # If it is not satisfied, (or we're forcing the issue) build it
    if force || should_build
        # Verbose mode tells us what's going on
        if verbose
            if !should_build
                @info("Force-building $(name) despite its satisfaction")
            else
                @info("Building $(name) as it is unsatisfied")
            end
        end

        # We setup our bash session to do three things:
        #   - Save bash history on every command issued (so we have a fake ~/.bash_history)
        #   - Save environment on quit (so that we can regenerate all environment variables
        #     when we're debugging, etc...)
        #   - Copy the current srcdir over to disk if something breaks.  This is done so
        #     that if we want to debug a build halfway through, we can get at it.  We normally
        #     build with a tmpfs mounted at `$WORKSPACE/srcdir`.
        trapped_script = """
        alias ll='ls -la'

        vecho() {
            if [[ "$(verbose)" == "true" ]]; then
                echo "\$@"
            fi
        }
        vecho_red() {
            (vecho "\$@" >&2)
        }

        # Save bash history (and optionally echo it out as it happens)
        save_history() {
            vecho_red " ---> \$BASH_COMMAND"
            history -s "\$BASH_COMMAND"
            history -a
        }

        # Save our environment into `/meta/.env`, eliminating read-only variables
        # so that this file can be sourced upon entering a debug shell.
        save_env() {
            set +x
            set > /meta/.env
            # Ignore read-only variables
            for l in BASHOPTS BASH_VERSINFO UID EUID PPID SHELLOPTS; do
                grep -v "^\$l=" /meta/.env > /meta/.env2
                mv /meta/.env2 /meta/.env
            done
            echo "cd \$(pwd)" >> /meta/.env
        }

        # We do a little sleight-of-hand here; we want to build inside of a tmpfs
        # because `srcdir` might be mapped in through a networked filesystem, which
        # totally wrecks our I/O performance.  So what we do instead is bind-mount
        # `srcdir` to another location so that we can always get at it, copy its
        # contents to a new tmpfs we mount at the location of `srcdir`, then when
        # we exit on an error, we copy everything back over again
        tmpify_srcdir() {
            vecho "Copying srcdir to tmpfs..."
            mkdir -p \$WORKSPACE/.true_srcdir
            mount --bind \$WORKSPACE/srcdir \$WORKSPACE/.true_srcdir
            mount -t tmpfs tmpfs \$WORKSPACE/srcdir
            rsync -rlptD \$WORKSPACE/.true_srcdir/ \$WORKSPACE/srcdir

            # We may have changed what pwd() means out from underneath ourselves
            cd \$(pwd)
        }

        # Copy our tmpfs version of `srcdir` back onto disk.
        save_srcdir() {
            vecho_red "Saving srcdir due to previous error..."
            rsync -rlptD \$WORKSPACE/srcdir/ \$WORKSPACE/.true_srcdir --delete
        }

        # If /meta is mounted, then we want to save history and environment
        if [[ -d /meta ]]; then
            trap save_history DEBUG
            trap "save_env" EXIT
            trap "save_env; save_srcdir" INT TERM ERR
        fi

        # Stop if we hit any errors.
        set -e

        # Swap out srcdir from underneath our feet
        tmpify_srcdir

        $(script)
        """

        logpath = joinpath(logdir(prefix), "$(name).log")
        did_succeed = run(runner, `/bin/bash -c $(trapped_script)`, logpath; verbose=verbose)
        if !did_succeed
            if debug
                @warn("Build failed, launching debug shell")
                run_interactive(runner, `/bin/bash --init-file /meta/.env`)
            end
            msg = "Build for $(name) on $(triplet(platform)) did not complete successfully\n"
            error(msg)
        end

        # Run an audit of the prefix to ensure it is properly relocatable
        if !skip_audit
            audit_result = audit(prefix; platform=platform,
                                 verbose=verbose, autofix=autofix,
                                 ignore_manifests=ignore_manifests)
            if !audit_result && !ignore_audit_errors
                msg = replace("""
                Audit failed for $(prefix.path).
                Address the errors above to ensure relocatability.
                To override this check, set `ignore_audit_errors = true`.
                """, '\n' => ' ')
                error(strip(msg))
            end
        end

        # Finally, check to see if we are now satisfied
        if !all(s(p) for p in products)
            if verbose
                @warn("Built $(name) but still unsatisfied!")
            end
        end
    elseif !should_build && verbose
        @info("Not building as $(name) is already satisfied")
    end
    return true
end

function print_buildjl(io::IO, products::Vector, product_hashes::Dict,
                       bin_path::AbstractString)
    print(io, """
    using BinaryProvider # requires BinaryProvider 0.3.0 or later

    # Parse some basic command-line arguments
    const verbose = "--verbose" in ARGS
    const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
    """)

    # Print out products
    print(io, "products = [\n")
    for prod in products
        print(io, "    $(repr(prod)),\n")
    end
    print(io, "]\n\n")

    # Print binary locations/tarball hashes
    print(io, """
    # Download binaries from hosted location
    bin_prefix = "$bin_path"

    # Listing of files generated by BinaryBuilder:
    """)

    println(io, "download_info = Dict(")
    for platform in sort(collect(keys(product_hashes)))
        fname, hash = product_hashes[platform]
        pkey = platform_key_abi(platform)
        println(io, "    $(pkey) => (\"\$bin_prefix/$(fname)\", \"$(hash)\"),")
    end
    println(io, ")\n")

    print(io, """
    # Install unsatisfied or updated dependencies:
    unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
    dl_info = choose_download(download_info, platform_key_abi())
    if dl_info === nothing && unsatisfied
        # If we don't have a compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform (\\\"\$(Sys.MACHINE)\\\", parsed as \\\"\$(triplet(platform_key_abi()))\\\") is not supported by this package!")
    end

    # If we have a download, and we are unsatisfied (or the version we're
    # trying to install is not itself installed) then load it up!
    if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
        # Download and install binaries
        install(dl_info...; prefix=prefix, force=true, verbose=verbose)
    end

    # Write out a deps.jl file that will contain mappings for our products
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
    """)
end

function print_buildjl(build_dir::AbstractString, src_name::AbstractString,
                       src_version::VersionNumber, products::Vector,
                       product_hashes::Dict, bin_path::AbstractString)
    mkpath(joinpath(build_dir, "products"))
    open(joinpath(build_dir, "products", "build_$(src_name).v$(src_version).jl"), "w") do io
        print_buildjl(io, products, product_hashes, bin_path)
    end
end

"""
If you have a sharded build on Github, it would be nice if we could get an auto-generated
`build.jl` just like if we build serially.  This function eases the pain by reconstructing
it from a releases page.
"""
function product_hashes_from_github_release(repo_name::AbstractString, tag_name::AbstractString;
                                            product_filter::AbstractString = "",
                                            verbose::Bool = false)
    # Get list of files within this release
    release = gh_get_json(DEFAULT_API, "/repos/$(repo_name)/releases/tags/$(tag_name)", auth=github_auth())

    # Try to extract the platform key from each, use that to find all tarballs
    function can_extract_platform(filename)
        # Short-circuit build.jl because that's quite often there.  :P
        if startswith(filename, "build") && endswith(filename, ".jl")
            return false
        end

        unknown_platform = typeof(extract_platform_key(filename)) <: UnknownPlatform
        if unknown_platform && verbose
            @info("Ignoring file $(filename); can't extract its platform key")
        end
        return !unknown_platform
    end
    assets = [a for a in release["assets"] if can_extract_platform(a["name"])]
    assets = [a for a in assets if occursin(product_filter, a["name"])]

    # Download each tarball, hash it, and reconstruct product_hashes.
    product_hashes = Dict()
    mktempdir() do d
        for asset in assets
            # For each asset (tarball), download it
            filepath = joinpath(d, asset["name"])
            url = asset["browser_download_url"]
            BinaryProvider.download(url, filepath; verbose=verbose)

            # Hash it
            hash = open(filepath) do file
                return bytes2hex(sha256(file))
            end

            # Then fit it into our product_hashes
            file_triplet = triplet(extract_platform_key(asset["name"]))
            product_hashes[file_triplet] = (asset["name"], hash)

            if verbose
                @info("Calculated $hash for $(asset["name"])")
            end
        end
    end

    return product_hashes
end


text:
function cpp_to_julia(var::CxxScalar{T}) where {T}
    var = Vec{T, 4}(var)
    return (var[1], var[2], var[3], var[4])
end

function cpp_to_julia(var::CxxVec{T, N}) where {T, N}
    return Vec{T, N}(var)
end

function julia_to_cpp(sc::Scalar)
    if size(sc,1)==0
        return CxxScalar{Float64}(0,0,0,0)
    elseif size(sc, 1) == 1
        return CxxScalar{Float64}(Float64(sc[1]), 0, 0, 0)
    elseif size(sc,1) == 2
        return CxxScalar{Float64}(Float64(sc[1]), Float64(sc[2]), 0, 0)
    elseif size(sc,1) == 3
        return CxxScalar{Float64}(Float64(sc[1]), Float64(sc[2]), Float64(sc[3]), 0)
    end
    return CxxScalar{Float64}(Float64(sc[1]), Float64(sc[2]), Float64(sc[3]), Float64(sc[4]))
end

function julia_to_cpp(vec::Vec{T, N}) where {T, N}
    return CxxVec{T, N}(Base.pointer(vec))
end

function julia_to_cpp(var::Array{T, 1}) where {T <: Scalar}
    ret = CxxWrap.StdVector{CxxScalar}()
    for x in var
        push!(ret, julia_to_cpp(x))
    end
    return ret
end

function julia_to_cpp(var::Array{Vec{T, N}, 1}) where {T, N}
    ret = CxxWrap.StdVector{CxxVec{T, N}}()
    for x in var
        push!(ret, julia_to_cpp(x))
    end
    return ret
end

function cpp_to_julia(var::CxxWrap.StdVector{T}) where {T <: CxxScalar}
    ret = Array{Scalar, 1}()
    for x in var
        push!(ret, cpp_to_julia(x))
    end
    return ret
end

function cpp_to_julia(var::CxxWrap.StdVector{CxxVec{T, N}}) where {T, N}
    ret = Array{Vec{T, N}, 1}()
    for x in var
        push!(ret, cpp_to_julia(x))
    end
    return ret
end

ext: jl


text:
@testset "error functions" begin
    @testset "real argument" begin
        @test erf(Float16(1))   ≈ 0.84270079294971486934 rtol=2*eps(Float16)
        @test erf(Float32(1))   ≈ 0.84270079294971486934 rtol=2*eps(Float32)
        @test erf(Float64(1))   ≈ 0.84270079294971486934 rtol=2*eps(Float64)

        @test erfc(Float16(1))  ≈ 0.15729920705028513066 rtol=2*eps(Float16)
        @test erfc(Float32(1))  ≈ 0.15729920705028513066 rtol=2*eps(Float32)
        @test erfc(Float64(1))  ≈ 0.15729920705028513066 rtol=2*eps(Float64)

        @test_throws MethodError erfcx(Float16(1))
        @test erfcx(Float32(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float32)
        @test erfcx(Float64(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float64)

        @test_throws MethodError logerfc(Float16(1))
        @test_throws MethodError logerfc(Float16(-1))
        @test logerfc(Float32(-100)) ≈ 0.6931471805599453 rtol=2*eps(Float32)
        @test logerfc(Float64(-100)) ≈ 0.6931471805599453 rtol=2*eps(Float64)
        @test logerfc(Float32(1000)) ≈ -1.0000074801207219e6 rtol=2*eps(Float32)
        @test logerfc(Float64(1000)) ≈ -1.0000074801207219e6 rtol=2*eps(Float64)
        @test logerfc(1000) ≈ -1.0000074801207219e6 rtol=2*eps(Float32)
        @test logerfc(Float32(10000)) ≈ log(erfc(BigFloat(10000, 100))) rtol=2*eps(Float32)
        @test logerfc(Float64(10000)) ≈ log(erfc(BigFloat(10000, 100))) rtol=2*eps(Float64)

        @test_throws MethodError logerfcx(Float16(1))
        @test_throws MethodError logerfcx(Float16(-1))
        @test iszero(logerfcx(0))
        @test logerfcx(Float32(1)) ≈ -0.849605509933248248576017509499 rtol=2eps(Float32)
        @test logerfcx(Float64(1)) ≈ -0.849605509933248248576017509499 rtol=2eps(Float32)
        @test logerfcx(Float32(-1)) ≈ 1.61123231767807049464268192445 rtol=2eps(Float32)
        @test logerfcx(Float64(-1)) ≈ 1.61123231767807049464268192445 rtol=2eps(Float32)
        @test logerfcx(Float32(-100)) ≈ 10000.6931471805599453094172321 rtol=2eps(Float32)
        @test logerfcx(Float64(-100)) ≈ 10000.6931471805599453094172321 rtol=2eps(Float64)
        @test logerfcx(Float32(100)) ≈ -5.17758512266433257046678208395 rtol=2eps(Float32)
        @test logerfcx(Float64(100)) ≈ -5.17758512266433257046678208395 rtol=2eps(Float64)
        @test logerfcx(Float32(-1000)) ≈ 1.00000069314718055994530941723e6 rtol=2eps(Float32)
        @test logerfcx(Float64(-1000)) ≈ 1.00000069314718055994530941723e6 rtol=2eps(Float64)
        @test logerfcx(Float32(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float32)
        @test logerfcx(Float64(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float64)

        @test_throws MethodError erfi(Float16(1))
        @test erfi(Float32(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float32)
        @test erfi(Float64(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float64)

        @test erfinv(Integer(0)) == 0
        @test_throws MethodError erfinv(Float16(1))
        @test erfinv(Float32(0.84270079294971486934)) ≈ 1 rtol=2*eps(Float32)
        @test erfinv(Float64(0.84270079294971486934)) ≈ 1 rtol=2*eps(Float64)

        @test erfcinv(Integer(1)) == 0
        @test_throws MethodError erfcinv(Float16(1))
        @test erfcinv(Float32(0.15729920705028513066)) ≈ 1 rtol=2*eps(Float32)
        @test erfcinv(Float64(0.15729920705028513066)) ≈ 1 rtol=2*eps(Float64)

        @test_throws MethodError dawson(Float16(1))
        @test dawson(Float32(1)) ≈ 0.53807950691276841914 rtol=2*eps(Float32)
        @test dawson(Float64(1)) ≈ 0.53807950691276841914 rtol=2*eps(Float64)
    end

    @testset "complex arguments" begin
        @test erf(ComplexF16(1+2im)) ≈ -0.53664356577856503399-5.0491437034470346695im
        @test erf(ComplexF32(1+2im)) ≈ -0.53664356577856503399-5.0491437034470346695im
        @test erf(ComplexF64(1+2im)) ≈ -0.53664356577856503399-5.0491437034470346695im

        @test erfc(ComplexF16(1+2im)) ≈ 1.5366435657785650340+5.0491437034470346695im
        @test erfc(ComplexF32(1+2im)) ≈ 1.5366435657785650340+5.0491437034470346695im
        @test erfc(ComplexF64(1+2im)) ≈ 1.5366435657785650340+5.0491437034470346695im

        @test_throws MethodError erfcx(ComplexF16(1))
        @test erfcx(ComplexF32(1+2im)) ≈ 0.14023958136627794370-0.22221344017989910261im
        @test erfcx(ComplexF64(1+2im)) ≈ 0.14023958136627794370-0.22221344017989910261im

        @test_throws MethodError erfi(ComplexF16(1))
        @test erfi(ComplexF32(1+2im)) ≈ -0.011259006028815025076+1.0036063427256517509im
        @test erfi(ComplexF64(1+2im)) ≈ -0.011259006028815025076+1.0036063427256517509im

        @test_throws MethodError erfinv(Complex(1))

        @test_throws MethodError erfcinv(Complex(1))

        @test_throws MethodError dawson(ComplexF16(1))
        @test dawson(ComplexF32(1+2im)) ≈ -13.388927316482919244-11.828715103889593303im
        @test dawson(ComplexF64(1+2im)) ≈ -13.388927316482919244-11.828715103889593303im
    end

    @testset "BigFloat arguments" begin
        @test erf(BigFloat(1))  ≈ 0.84270079294971486934 rtol=2*eps()

        @test erfc(BigFloat(1)) ≈ 0.15729920705028513066 rtol=2*eps()

        @test erfcx(BigFloat(1))        ≈ 0.42758357615580700442    rtol=2*eps()
        @test erfcx(BigFloat(2e9))      ≈ 2.820947917738782e-10     rtol=2*eps()
        @test erfcx(BigFloat(1.8))      ≈ erfcx(1.8)                rtol=4*eps()
        @test erfcx(BigFloat(1.8e8))    ≈ erfcx(1.8e8)              rtol=4*eps()
        @test erfcx(BigFloat(1.8e88))   ≈ erfcx(1.8e88)             rtol=4*eps()
        @test isnan(erfcx(BigFloat(NaN)))

        @test logerfc(BigFloat(1000, 100)) ≈ -1.0000074801207219e6 rtol=2*eps(Float64)
        @test isnan(logerfc(BigFloat(NaN)))

        @test_throws MethodError erfi(big(1.0))

        @test_throws MethodError erfinv(BigFloat(1))

        @test_throws MethodError erfcinv(BigFloat(1))

        @test_throws MethodError dawson(BigFloat(1))
    end

    @testset "inverse" begin
        for elty in [Float32,Float64]
            for x in exp10.(range(-200, stop=-0.01, length=50))
                @test isapprox(erf(erfinv(x)), x, atol=1e-12*x)
                @test isapprox(erf(erfinv(-x)), -x, atol=1e-12*x)
                @test isapprox(erfc(erfcinv(2*x)), 2*x, atol=1e-12*x)
                if x > 1e-20
                    xf = Float32(x)
                    @test isapprox(erf(erfinv(xf)), xf, atol=1e-5*xf)
                    @test isapprox(erf(erfinv(-xf)), -xf, atol=1e-5*xf)
                    @test isapprox(erfc(erfcinv(2xf)), 2xf, atol=1e-5*xf)
                end
            end
            @test erfinv(one(elty))  == Inf
            @test erfinv(-one(elty)) == -Inf
            @test_throws DomainError erfinv(convert(elty,2.0))

            @test erfcinv(zero(elty)) == Inf
            @test_throws DomainError erfcinv(-one(elty))
        end

        @test erfinv(one(Int))  == erfinv(1.0)
        @test erfcinv(one(Int)) == erfcinv(1.0)
    end
end


text:
include("roadway/runtests.jl")
include("vehicles/runtests.jl")
include("test_actions.jl")
include("utils/runtests.jl")


text:
# DAY 4 PROBLEM 1 ADVENT OF CODE
# Get the winning board and calculate the output.

data = """1,76,38,96,62,41,27,33,4,2,94,15,89,25,66,14,30,0,71,21,48,44,87,73,60,50,77,45,29,18,5,99,65,16,93,95,37,3,52,32,46,80,98,63,92,24,35,55,12,81,51,17,70,78,61,91,54,8,72,40,74,68,75,67,39,64,10,53,9,31,6,7,47,42,90,20,19,36,22,43,58,28,79,86,57,49,83,84,97,11,85,26,69,23,59,82,88,34,56,13

85 23 65 78 93
27 53 10 12 26
 5 34 83 25  6
56 40 73 29 54
33 68 41 32 82

 8 31 14 70 91
53 49 86 13 21
66 28 76 78 93
39 63 80 43 23
56 25 60 67 72

67 78 36 64 14
46 16 80 23 94
22 47 51 65 57
33 76 21 92 97
31 95 54 27 20

 1 77 86 43 30
28 88  7  5 60
66 24  3 57 33
38 23 59 84 44
74 47 17 29 85

21 50 86  2 70
85 19 22 93 25
99 38 74 30 65
81  0 47 78 63
34 11 51 88 64

45 15 29 81 30
75 21 88 91 49
39 20  4 17 78
10 12 38 11  7
98  6 65 69 86

36 20 31 44 69
30 65 55 88 64
74 85 82 61  5
57 17 90 43 54
58 83 52 23  7

42 16 82 86 76
60 26 27 59 55
 7 53 22 78  5
18 61 10 15 17
28 46 14 87 77

21 43 15 47 61
24 76 28  3 27
19 62 69 82 93
49 29 97 74 41
92 36 37 99 40

31  4  3 62 51
24 57 78 67 53
13  5 76 38 55
79  9 75 98 71
65  1 39 18 47

59  4 38 95 99
85 68 69 93 43
83 57 48 42 15
47 50 80 79 90
56 87 78 64 25

21 37 14 67 95
88 39 26 38 49
89 83 54 77 96
48 86 94 19 20
43 41  8 74 58

 1 36 12 90 91
63 21 98 82 66
39 86  7 52 77
80 81 44 33 58
78 30 11 51 28

81 74  7 33 96
75 60 87 47 91
39 73 30 50 13
 4 41  9 43 77
34 82 72 48 12

93 63 74 25 57
29 76  9 45 70
98 77 71 16 41
47 54 18 14 55
31 89 67 87 83

 8 72 45 93 68
74 26 69 94 65
28  9 20 47 41
46 54 21 56 22
84 62 18 15 48

20 51 81 40 69
71 10 13 93 75
44 86  0 95 37
99 39 76 80 66
14 64 49 62 27

75  7 51 86 79
43 30 61 39 16
85 63 90 28 96
88 78 72 31 73
98 87 23 19 58

20 95 47 97 12
92 25 68 87 91
37 10 78 23 63
74 93 58 39  5
76 51 48 72 16

37 18 32 34 85
22 31 98 42 19
29 72 48 76 25
47  1 21  7 53
79 82 86 52 78

20 16 47 78 92
88 15 71 67  2
 5 52 90 70  9
22 49 28 82 27
 6 19 61 73 48

71 26  7 11 79
52 30 47  1 31
17 75 94 91 28
81 98 23 55 21
77 15 39 24 16

 5 75 44 88 65
89 45 23 69 19
41 61 67 52 54
47 38 57 12 98
62 70 26 87 53

50  4 65 77 25
 6 21  5 27 92
39 63 97 75 79
60 34 87 26 74
99 24 44 85  2

13 64 38 78 21
74 17 83 57 94
25 39 69 53  4
54 33 81 50 76
42 75 19 77 26

63 31 70 19 39
38 87 15 90 75
61 98  6 29 86
78 62 32 11 60
55 97 13 73 82

51 63 68 84 36
12 33 37 31  8
18 41 34 74 23
72 39 85 48 60
24 19 29 88  0

46 51 17 23 13
20 93 97 99 81
57 47 33 84 44
28 96  2 43 56
68 36 62 15  5

81 99  5 30 10
38 62 57  8 37
 7 86 98  3 54
46 82 96 15 72
83  1 75 25 50

47 57 11 61 27
53 10 31 91 98
76 85 55 38 23
 6 81 67 71 70
35 29 17 50 56

24 65 15  1 89
45 60 97 23 14
84 56 58  5 54
 3 72 51 46 79
67 70 78 34 77

38 11 54 23  2
33 14 10 96 63
43  5 36 20 30
70 53 66 71  9
91 90 21  7 88

94 44  4 86 26
39 70 54 50 30
55 40 12 72 71
68  7 66 47 91
31 24 13  1 96

79 14 40 87 68
16 32 53 46 98
38 95 21 89 69
62 60 19 81 33
70 52 28 83  0

62 42 38 48 64
61 79 78 97 98
89  7  3 29 68
92 76 14 67  1
41 99 72 47 60

 5 75 18 42 33
72 61 36 31 29
19 58  1 34 94
54 84 92 99 38
76 68 79 53 37

14 91 37  5 98
68 29 34 76 43
75  0 67 33 69
81 47 58 30 93
88 92 42 77 54

64 24 28 54 53
72 68  3 73  4
83  6 59 66 94
87 80 55 20 16
13 82 74 31 70

63 92 71  0 83
98 40 50 55  2
88  5 85 30 23
10 75 81 58 68
51 31 14 89  1

67 93 94 54 53
38 71 34 40 24
31 63 30 99 75
 4 57 86 19 70
60 49 87 68 74

56 94 79 53  7
24 12 19  6 99
82 51 41 46 43
17 49 52 78 55
75 48 61 70 87

14 55 32 21 31
88 83 23 44  4
 1 77 45 90 85
46 81 51 27 62
60 24 29 18  0

95 92 91 27 26
22 43 45 64 62
83 23 25 85 94
84 53 72 28 20
75 60 52 18 73

95 41  7 21 32
58 65 16 56 97
68 25 91 83 24
66 89 15 55  6
 2 30 84 10 90

58 86 44 19 74
57 89 17  6 83
77 35 60 32 13
97 63 62 28 76
55 31 11  0 52

33 39 59 42 45
61 50 92  9 79
15  0 28  5 72
91 24 21 29 87
86 76 43 31 93

63 11 86 45 85
96 74 66 93 32
95 30 99 23 18
69 97 48 15  1
42 87 47 83 80

93  5 40 64  2
44 51 15 54 83
69 77 90 58 11
 0 48 43 30 55
25 72 38 73 52

89 58 71 68 15
23 65  9 36 74
21 29 42 79 98
55 47 33 39 28
16 75 91 69 57

13 79 12 71  2
60 94 99 43 82
84 89 29 91 87
74 80 25 32 21
70 14 68 92 11

78  1 16 51 87
58 94 59 15 43
79 41 50 47 39
53 37  9 28 72
34 63 89 35 18

31 67 70 42 43
60  2 89 49 22
56 17 81 24 74
20 65  1 96 51
68  7  0 38 25

59 14 29 53 19
 9  2 11 33 44
81  6 10 47 58
20 34 62 55 40
71 38 69 45 78

59 36 70 42 21
 3 16 49 79 98
74 25  8 84 19
61 80 47 65 64
91 62 52  9 40

 1 85 63  7  2
 0 20 61 26 77
99 37 74 42 76
25 94 19 78 60
79 72 95 22 11

51 21 79 76 32
55 23 69 19 61
71 54 94 47 92
 5 64  6 68 16
91 81  9 99 30

61 69 82 86 68
66 81 28 38 36
26 29 31 11  8
72 51 12 95 63
18 30 88 17 32

34  8 14 42 67
66 79 65 20 52
37 87 74 24  3
59 54 21 32 89
31  4 62 76 30

11 93  8 92 55
38 72 99  3 83
12 75  0 41 46
17 25  5 39 48
14 18 86 29 84

 6 20 41 51 48
 5 67 30 24 47
 3  8 92 22 39
 4 56 36 31 75
 2 45 85 81 96

47 43 72 22  3
19 87 53 12 60
29 40 56 68 18
66 97 70 33 39
85 37  0 90 98

61 35 81 84 94
11  1 58 45 77
 6 99 67 36 43
 5  7  0 87 80
44 78 39 70 20

58 34 49 29 75
17 15 28 23 84
59 25 92 48  0
20 81 47  3 71
68 60  5 22 87

90 32 41 39  6
36 78 67 24 50
55 72 52 75 44
87 15 92 31 58
83 89 68 19 43

99 44 53 68 25
71 67 16 19 36
35 58 14 86 48
88 18 61 24 23
87  9 91 37 15

37  5 63 68 28
41 50 76 99 64
34 92 78 94 71
11 96 97 42 58
33 45  0 93 48

33 68  9 12 81
60 98 28  8 99
14 17  6 82 15
57 69 43 38 29
47 84 76 22 18

79 70 92 38 47
12 82 98 46  0
76 15 53 59 97
18 52 49 29 96
44 64 68 89 24

95 14 17 27 42
55 43 57 29 25
34 73 86 50 16
69 37 75 63 39
78 79  3  4 30

27 31 15 92 46
36 23 72 40 50
51 99 55 89 21
12 70 84 63 85
78 88 77 75  0

15 67 40 39 28
 9 79 22 52 75
96 65 86 98 14
97 87 44 84 68
36 26 89 43 27

79 59 48 27 36
85 92 93 76 24
 2 25  7 42 90
23 29 74 35 86
58 60 31 75 57

10 43 83 75  8
88 12 38 30  9
60 67 59 76  6
55 45 74 34 25
97 49 65 96 69

59 86 15  3 19
89  4 74 61 23
52 98  8 79 39
95 17 22 14 51
50 18 94 30 84

19 63 58 72 67
35 93 29 91  0
39 26 43 84 21
70 42  2 53 12
59 99  8  1 86

23 86 34 22 65
71 10 16 50 91
66 89 49 81 43
40  7 26 75 61
62 59  2 46 95

24 21  0 49 25
92 42 48 12  7
81 93 59 68  3
14 23 63 39 29
35 43  6 44 89

67 74 95 34 10
39 90 59 44 51
17 16 97 24 62
20 54 76 63 88
87 66 14 78 82

96 86 67 59 79
66  3 30 77 71
 2 91 99 82 31
48 65 75 98 53
63 54 64 76  1

85 96 40 98 24
16 20 10 23 17
79 59 53 42 65
67  2  5 80 75
62 38 19 74 73

43 10 79 92  8
52 36  4  5 67
56 29 33 24 97
85 17 53 75 65
62 64  1 21 83

93 92 79 17 12
40 88  6 82 34
90 96 53 25 43
14 62 54 10 39
49 68 41 16 44

67 99 24 58 76
43 53 59 54 51
47  6 61  8  2
80 68 90 14  4
29 46 94 89 50

14 45 19 33 43
 6 55  4 31 80
51  2 69 68 61
71 70 79 91 93
66 18 54 13 87

 8 45 61 54 30
85 16 19 82 37
56 39 11 47  4
74 70 10 60 91
21 63 95 53 72

71 21 63 86 27
53 52 40 23 81
 2 47 92 68 15
46 45 31  8  1
34 80 37 11 69

96  0 15 90 66
65 43 92 83 18
 3 47 19  8 32
71 26 42 34 28
62 99 55  5 12

37 99 30 21  3
63 18 68 47 27
57  0 65 85 20
 7 58 40 92 43
15 19  5  4 53

46 16 45 95 68
 6 44 31 47 73
84 82 71 75 94
26 25 17 32 49
18 96 13 58  9

71 36 13 68 10
84  7 60 79 41
 1 83 43 81 97
90 53 80 19 38
48 25 32 42 29

37 68 86 44 78
87 67 77 70 60
45 34 27 15 47
12 21 13 55 26
81 41 63 40 74

24 50 93 94 57
99  4 56  5 28
42 31 22  6 76
90 89 16 49 59
 9  7 43 71 54

69 75 94 38 46
52 64 50 72 42
76 63 13 60 10
99 80 43 33 17
25 31  4 89 22

88 57 22 66 34
85 16 87 95 59
73  2 46  5 29
25 69 53  6 14
96 77 19 91 43

46 99 52 47 76
89 53 24 13 59
45  5  1 30 19
68 25 22 10 73
42 27 31  0 94

42 44 98 89 87
65 10 80 56 41
 3 35 95 48 43
85 97 83 12 94
50 38 93 47 17

16 73 18 81 89
 6 48 54 93 19
35 52 88 49 31
43 79 83 14 28
50 62 98 26 22

38 47  7 20 35
45 76 63 96 24
98 53  2 87 80
83 86 92 48  1
73 60 26 94  6

80 50 29 53 92
66 90 79 98 46
40 21 58 38 60
35 13 72 28  6
48 76 51 96 12

79 80 24 37 51
86 70  1 22 71
52 69 10 83 13
12 40  3  0 30
46 50 48 76  5"""

"""
"""
function parse_board(board::AbstractString) ::Vector{Vector{Integer}}
    return [[parse(Int16, n) for n in split(row, ' ') if !isempty(n)]
            for row in split(board, '\n')]
end

"""
"""
function parse_boards(boards::AbstractString) ::Tuple{Vector{Integer}, Vector{Vector{Vector{Integer}}}}
    # Get individual boards and number draws.
    boards = split(boards, "\n\n")
    # Get number draws.
    numbers = [parse(Int16, n) for n in split(first(boards), ',')]
    # return number draws, the parsed boards.
    return numbers, [parse_board(b) for b in boards[2:end]]
end

"""
"""
function get_board_dimensions(board::Vector{Vector{Integer}}) ::Tuple{Integer, Integer}
    return length(board), length(first(board))
end

"""
"""
function get_valid_board_start_indexes(board::Vector{Vector{Integer}}; n::Integer=5) ::Vector{Tuple{Integer, Integer}}
    # Get the dimensions of table and maximum x and y to get n adjacent items.
    y_dim, x_dim = get_board_dimensions(board)
    x_max = x_dim - (n-1)
    y_max = y_dim - (n-1)
    # Get indexes to start from to get n successive items.
    indexes::Vector{Tuple{Integer, Integer}} = []
    # Get x indexes.
    for (x,y) in Iterators.product(1:x_max, 1:y_dim)
        push!(indexes, (x,y))
    end
    # Get y indexes.
    for (x,y) in Iterators.product(2:x_dim, 1:y_max)
        push!(indexes, (x,y))
    end
    return indexes
end

"""
"""
function get_sequences(board::Vector{Vector{Integer}}, x::Integer, y::Integer; n::Integer=5) ::Vector{Vector{Integer}}
    # Establish return.
    returnval::Vector{Vector{Integer}} = []
    # Get the dimensions of table and maximum x and y to get n adjacent items.
    y_dim::Integer, x_dim::Integer = get_board_dimensions(board)
    x_max::Integer = x + n - 1
    y_max::Integer = y + n - 1
    # Check if we have enough items left to get the column and row to get a sequence of n.
    row_valid::Bool = x_max <= x_dim
    col_valid::Bool = y_max <= y_dim
    # Add items.
    if row_valid
        push!(returnval, board[y][x:x_max])
    end
    if col_valid
        push!(returnval, [board[i][x] for i=y:y_max])
    end
    if (row_valid && col_valid)
        push!(returnval, [board[y+i][x+i] for i=0:n-1])
    end
    # handle up and right diagonal case.
    if (row_valid && (y - n + 1) >= 1)
        push!(returnval, [board[y-i][x+i] for i=0:n-1])
    end
    return returnval
end

"""
TODO solve this function.
"""
function get_board_axes(board::Vector{Vector{Integer}}; n::Integer=5) # ::Vector{Vector{Integer}}
    """
    1 1 1 1 1
    1 0 0 0 0
    1 0 0 0 0
    1 0 0 0 0
    1 0 0 0 0

    1 1 1 1 1 1
    1 1 0 0 0 0
    1 1 0 0 0 0
    1 1 0 0 0 0
    1 1 0 0 0 0

    If we come up with an iteration method that hits only the cells marked 1
    for each array, we can use them to get items on the diagonal, horizontal,
    and vertical with length 5.
    """
    return [i
            for (x, y) in get_valid_board_start_indexes(board, n=n)
            for i in get_sequences(board, x, y, n=n)]
end

"""
"""
function is_winner(board::Vector{Vector{Integer}}, draws::Vector{<:Integer}; n::Integer=5) ::Bool
    return any(all(x in draws for x in axis)
               for axis in get_board_axes(board, n=n))
end

"""
"""
function get_winning_board(boards::Vector{Vector{Vector{Integer}}}, draws::Vector{Integer}; n::Integer=5) ::Tuple{Vector{Vector{Integer}}, Vector{Integer}}
    # Iterate from n to all numbers drawn.
    for i=n:length(draws)
        drawn = draws[begin:i]
        # See if any boards win given this set of draws.
        for board in boards
            if is_winner(board, drawn, n=n)
                return board, drawn
            end
        end
    end
    return nothing
end

"""
"""
function calculate_board_score(board::Vector{Vector{Integer}}, drawn::Vector{Integer}) ::Integer
    score = sum(x for x in Iterators.flatten(board) if !(x in drawn))
    return score * drawn[end]
end

draws, boards = parse_boards(data)
winning_board, drawn = get_winning_board(boards, draws)
calculate_board_score(winning_board, drawn)


text:
# Miscellaneous utilities

"""
```
crossover(x::Array{T}, y::Array{T}) where {T<:Real}
```

Find where `x` crosses over `y` (returns boolean vector where crossover occurs)
"""
function crossover(x::Array{T}, y::Array{T}) where {T<:Real}
    @assert size(x,1) == size(y,1)
    out = falses(size(x))
    @inbounds for i in 2:size(x,1)
        out[i] = ((x[i] > y[i]) && (x[i-1] < y[i-1]))
    end
    return out
end

"""
```
crossunder(x::Array{T}, y::Array{T}) where {T<:Real}
```

Find where `x` crosses under `y` (returns boolean vector where crossunder occurs)
"""
function crossunder(x::Array{T}, y::Array{T}) where {T<:Real}
    @assert size(x,1) == size(y,1)
    out = falses(size(x))
    @inbounds for i in 2:size(x,1)
        out[i] = ((x[i] < y[i]) && (x[i-1] > y[i-1]))
    end
    return out
end


text:
export derivatives, der

"Outer product"
⊗(a::Vector{Float64},b::Vector{Float64}) = a*b'

"Dictionary associating functions with their derivatives"
derivatives = Dict{Function, Function}()

"Derivative of a function"
der(f::Function) = derivatives[f]


text:
"""
    panas(f,c1,c2,a,b,c,m,n)
berisi program untuk mencari solusi
persamaan differensial parsial yaitu persamaan panas. Program ini secara default berisi
8 masukan, yaitu fungsi nilai awal `f(x)`, koefisien nilai batas `c1` dan `c2`, nilai batas solusi
`a` dan `b`, koefisien persamaan gelombang `c`, dan jumlah grid solusi x dan t yaitu `m` dan
`n`.
# Example
```jl
julia> f(x) = 4*x - 4*x.^2;

julia> U,R = panas(f,0,0,1,0.2,1,6,11);

julia> R
0.4999999999999999

julia> U
11×6 Adjoint{Float64,Array{Float64,2}}:
 0.0  0.64       0.96      0.96      0.64       0.0
 0.0  0.48       0.8       0.8       0.48       0.0
 0.0  0.4        0.64      0.64      0.4        0.0
 0.0  0.32       0.52      0.52      0.32       0.0
 ⋮                                              ⋮
 0.0  0.1375     0.2225    0.2225    0.1375     0.0
 0.0  0.11125    0.18      0.18      0.11125    0.0
 0.0  0.09       0.145625  0.145625  0.09       0.0
 0.0  0.0728125  0.117813  0.117813  0.0728125  0.0
```
return solusi persamaan panas `U` dan rasio `R`.
"""
function panas(f,c1,c2,a,b,c,m,n)
  h = a/(m-1);
  k = b/(n-1);
  r = c^2*k/h^2;
  U = zeros(m,n);
  U[1,:] .= c1;
  U[m,:] .= c2;
  U[2:m-1,1] = f.(h:h:(m-2)*h)';
  for j = 2:n
    for i = 2:m-1
      U[i,j]=(1-2*r)*U[i,j-1]+ r*(U[i-1,j-1]+U[i+1,j-1]);
    end
  end
  U=U';
  return U,r
end


text:
function getswiphase(data, options)
    mask = robustmask(view(data.mag,:,:,:,1))
    savenii(mask, "maskforphase", options.writesteps, data.header)
    # TODO output not readable
    combined = getcombinedphase(data, options, mask)
    swiphase = createphasemask!(combined, mask, options.phase_scaling_type, options.phase_scaling_strength)
    savenii(swiphase, "swiphase", options.writesteps, data.header)
    return swiphase
end

function createphasemask!(swiphase, mask, phase_scaling_type, phase_scaling_strength)
    swiphase[.!mask] .= 0
    pos = swiphase .> 0

    if phase_scaling_type == :negativetanh
        swiphase = -swiphase
        phase_scaling_type = :tanh
    end

    if phase_scaling_type == :tanh
        m = median(swiphase[mask .& (swiphase .> 0)])
        m *= 10 / phase_scaling_strength
        f(x) = (1 + tanh(1 - x/m)) / 2
        swiphase .= f.(swiphase)

    elseif phase_scaling_type == :positive
        swiphase[.!pos] .= 1
        swiphase[pos] .= robustrescale(swiphase[pos], 1, 0) .^ phase_scaling_strength

    elseif phase_scaling_type == :negative
        swiphase[pos] .= 1
        swiphase[.!pos] .= robustrescale(swiphase[.!pos], 0, 1) .^ phase_scaling_strength

    elseif phase_scaling_type == :triangular
        swiphase[pos] .= robustrescale(swiphase[pos], 1, 0) .^ phase_scaling_strength
        swiphase[.!pos] .= robustrescale(swiphase[.!pos], 0, 1) .^ phase_scaling_strength

    else
        error("$phase_scaling_type not defined for SWI")
    end
    swiphase[swiphase .< 0] .= 0
    return swiphase
end

function getcombinedphase(data, options, mask)
    phase = data.phase
    mag = data.mag
    TEs = to_dim(data.TEs, 4)
    σ = options.phase_hp_σ

    unwrapped = similar(phase)
    if options.phase_unwrap == :laplacian
        for iEco in 1:size(phase, 4)
            unwrapped[:,:,:,iEco] .= laplacianunwrap(view(phase,:,:,:,iEco))
        end
        savenii(unwrapped, "unwrappedphase", options.writesteps, data.header)

        for iEco in 1:size(phase, 4)
            smoothed = gaussiansmooth3d(unwrapped[:,:,:,iEco], σ; mask, dims=1:3)
            unwrapped[:,:,:,iEco] .-= smoothed
        end
        combined = combine_phase(unwrapped, mag, TEs)
        savenii(unwrapped, "filteredphase", options.writesteps, data.header)
        savenii(combined, "combinedphase", options.writesteps, data.header)

    elseif options.phase_unwrap == :laplacianslice
        for iEco in 1:size(phase, 4), iSlc in 1:size(phase, 3)
            unwrapped[:,:,iSlc,iEco] .= laplacianunwrap(view(phase,:,:,iSlc,iEco))
            smoothed = gaussiansmooth3d(unwrapped[:,:,iSlc,iEco], σ; mask=mask[:,:,iSlc], dims=1:2)
            unwrapped[:,:,iSlc,iEco] .-= smoothed
        end
        combined = combine_phase(unwrapped, mag, TEs)
        savenii(unwrapped, "unwrappedphase", options.writesteps, data.header)
        savenii(combined, "combinedphase", options.writesteps, data.header)

    elseif options.phase_unwrap == :romeo
        unwrapped = romeo(phase; mag, TEs)#, mask = mask)
        savenii(unwrapped, "unwrappedphase", options.writesteps, data.header)

        combined = combine_phase(unwrapped, mag, TEs)
        savenii(combined, "combinedphase", options.writesteps, data.header)

        combined .-= gaussiansmooth3d(combined, σ; mask, dims=1:2)
        savenii(combined, "filteredphase", options.writesteps, data.header)

    else
        error("Unwrapping $options.phase_unwrap ($(typeof(options.phase_unwrap))) not defined!")
    end

    return combined
end

combine_phase(unwrapped::AbstractArray{T,3}, mag, TEs) where T = copy(unwrapped) # one echo
function combine_phase(unwrapped::AbstractArray{T,4}, mag, TEs) where T
    dims = 4
    TEs = reshape(TEs, ones(Int, dims-1)..., length(TEs)) # size = (1,1,1,nEco)

    combined = sum(unwrapped .* TEs .* mag .* mag; dims)
    combined ./= sum(mag .* mag .* float.(TEs).^2; dims)
    return dropdims(combined; dims)
end


text:
@testset "Rectilinear Grid" begin
    @testset "Rectilinear Grid - 2D" begin
        for T in Floats
            x = Vec{5,T}(0, 1//4, 1//2, 3//4, 1)
            y = Vec{4,T}(0, 1//3, 2//3, 1)
            rg = RectilinearGrid(x, y)
            @test rg.x == x
            @test rg.y == y
            @test xmin(rg) == 0
            @test ymin(rg) == 0
            @test xmax(rg) == 1
            @test ymax(rg) == 1
            @test RectilinearGrid(Vec{2,T}(3//4, 1), Vec{2,T}(0, 1//3)) ⊆ rg
        end
    end

    @testset "Rectilinear Grid - 3D" begin
        for T in Floats
            x = Vec{5,T}(0, 1//4, 1//2, 3//4, 1)
            y = Vec{4,T}(0, 1//3, 2//3, 1)
            z = Vec{4,T}(0, 1//3, 2//3, 1)
            rg = RectilinearGrid(x, y, z)
            @test rg.x == x
            @test rg.y == y
            @test rg.z == z
            @test xmin(rg) == 0
            @test ymin(rg) == 0
            @test zmin(rg) == 0
            @test xmax(rg) == 1
            @test ymax(rg) == 1
            @test zmax(rg) == 1
            @test RectilinearGrid(Vec{2,T}(3//4, 1),
                                  Vec{2,T}(0, 1//3),
                                  Vec{2,T}(0,    1)) ⊆ rg
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    Q, R, C, prow, pcol, info = spqr_wrapper(A, B, tol, get_details)

SPQR_WRAPPER wrapper around spqr to get additional statistics.
The ouput `R` has reduced rown count, if final diagonal elements would become below `tol`.
Note: In addition to the original, a row pivoting array is returned
Maybe that should be absorbed into Q in the sparse case of `qr()`.
"""
function spqr_wrapper(A::AbstractArray{T}, B::AbstractArray{T}, tol, get_details) where T
    # Copyright 2012, Leslie Foster and Timothy A Davis.

    if get_details != 0
        # get detailed statistics for time and memory usage
        t = time_ns()
    else
        t = UInt(0)
    end

    m, n = size(A)
    info = NamedTuple()

    gt(tol::AbstractFloat) = x::AbstractFloat -> abs(x) > tol

    function _spqr_wrapper(A::AbstractSparseArray, tol)
        m, n = size(A)
        qrf = qr(A, tol=tol)
        Q, R, prow, pcol = qrf.Q, qrf.R, qrf.prow, qrf.pcol
        rk = maximum(R.rowval)
        norm_E_fro = sqrt(sum((sum(abs2,A,dims=1)[pcol] - sum(abs2,R,dims=1)[:])[rk+1:end]))
        Q, R, prow, pcol, rk, norm_E_fro
    end

    _spqr_wrapper(A::Adjoint{<:Any,<:AbstractSparseMatrix}, t) = _spqr_wrapper(sparse(A), t)

    function _spqr_wrapper(A::AbstractArray, tol)
        qrf = qr(A, Val(true)) # enable column pivoting
        Q, R, prow, pcol = qrf.Q, qrf.R, [1:size(A, 1)...], qrf.p
        rm, rn = size(R)
        k = findlast(gt(tol), diag(R))
        rk = k == nothing ? 0 : k
        norm_E_fro = norm(view(R, k+1:rm, k+1:rn))
        # norm_E_fro1 = norm(A[prow,pcol] - Q * R, 2) much more effort!
        R[k+1:end,:] .= 0
        Q, R, prow, pcol, rk, norm_E_fro
    end

    Q, R, prow, pcol, rank, norm_E_fro = _spqr_wrapper(A, tol)
    R = striprows(R, rank)

    if isempty(B)
        # C is empty
        C = zeros(T, m, 0)
    else
        # also compute C = Q'*B if B is present
        C = Q' * B
    end

    t = get_details != 0 ? time_ns() - t : t
    info = (rank_A_estimate=rank, norm_E_fro=norm_E_fro, time=t)

    Q, R, C, prow, pcol, info
end

function striprows(A::SparseMatrixCSC, k::Int)
    m, n = size(A)
    if k < m
        SparseMatrixCSC(k, n, A.colptr, A.rowval, A.nzval)
    else
        A
    end
end

function striprows(A::AbstractMatrix, k::Int)
    m, n = size(A)
    if k < m
        A[1:k,:]
    else
        A
    end
end

"""
    Calculate `norm((Q' * A[prow,pcol])[k+1:end,k+1:end])`
"""
function errornorm(Q::AbstractMatrix, A::AbstractMatrix{T}, k::Int, prow::Vector{Int}, pcol::Vector{Int}) where T<:Number

    m, n = size(A)
    B = Vector{T}(undef, m)
    sum = real(T)(0)
    if k < m
        for i = k+1:n
            copyto!(B, view(A, prow, pcol[i]))
            if norm(B, Inf) != 0
                lmul!(Q', B)
                j = k
                while j < m
                    j += 1
                    sum += abs2(B[j])
                end
            end
        end
    end
    sqrt(sum)
end


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

import Test

# Tests type inference on Julia >= v1.6 only:
macro tinferred(expr)
    if VERSION >= v"1.6"
        esc(:(Test.@inferred($expr)))
    else
        esc(:($expr))
    end
end

using LinearAlgebra, StaticArrays

sum_pow2(x) = sum(map(x -> x^2, x))
sum_pow3(x) = sum(map(x -> x^3, x))
f = (xs...) -> (sum(map(sum_pow2, xs)), sum(map(sum_pow3, xs)), 42)
xs = (2f0, (3, 4f0), SVector(5f0, 6f0, 7f0))
ΔΩ = (10, 20, 30)

fy = let f = f; (xs...) -> SVector(f(xs...)); end
ΔΩy = SVector(ΔΩ)


text:
using Arblib

"""
    parse_arbdoc(filename)
Parse a .rst file from the Arb documentation. Returns the title of the
document together with a list of sections with their titles and the
functions documented in them.
"""
function parse_arbdoc(filename)
    return open(filename) do file
        lines = readlines(file)

        title = lines[3]
        @assert startswith(lines[4], "======")

        # Quick and dirty way to find all the sections
        splits = findall(str -> startswith(str, "------"), lines)

        sections = Vector{Tuple{String,Vector{String}}}()
        for n = 1:length(splits)
            section_title = lines[splits[n]-1]
            push!(sections, (section_title, []))

            start = splits[n]
            stop = n == length(splits) ? length(lines) : splits[n+1] - 1

            for line in lines[start:stop]
                if startswith(line, ".. function::")
                    push!(sections[n][2], strip(line[length(".. function::")+1:end]))
                end
            end
        end

        return title, sections
    end
end

"""
    generate_file(title, sections; verbose, manual_overrides)
    generate_file(filename, title, sections; verbose, manual_overrides)
Given a title and a list of sections, as returned by
[`parse_arbdoc`](@ref), return a string with the Julia code to load
all of those functions. If a file name is given then also write the
string to that file.

For each function it checks that it is correctly parsed. Those that do
not parse due to types which are not supported and likely never will
be are commented out in the code with `ns` (not supported). Those that
do not parse due to types that are currently not implemented, though
might be so in the future, are commented out with `ni` (not
implemented). If they do not parse for some other reason an error is
thrown.

If the argument `manual_overrides` is given then comment out functions
which exists as keys in the dictionary. It adds the value of the key
as an extra comment.
"""
function generate_file(
    filename,
    title,
    sections;
    verbose = false,
    manual_overrides = Dict{String,String}(),
)
    str = generate_file(
        title,
        sections,
        verbose = verbose,
        manual_overrides = manual_overrides,
    )
    open(filename, "w") do file
        write(file, str)
    end
    return str
end

function generate_file(
    title,
    sections;
    verbose = false,
    manual_overrides = Dict{String,String}(),
)
    str = "###\n"
    str *= "### " * title * "\n"
    str *= "###\n"

    num_functions = 0
    num_unsupparg = 0
    num_keyerror = 0

    unsuppargs = Set{String}()
    keyerrors = Set{String}()

    for (section_title, functions) in sections
        str *= "\n### " * section_title * "\n"

        for s in functions
            num_functions += 1
            if s in keys(manual_overrides)
                str *= "#mo arbcall\"" * s * "\" # $(manual_overrides[s])\n"
            else
                try

                    f = Arblib.Arbfunction(s)

                    s == Arblib.arbsignature(f) ||
                        @warn( "Expected signature: $s\n Obtained signature: $(Arblib.arbsignature(f))")

                    str *= "arbcall\"" * s * "\"\n"
                catch e
                    if e isa Arblib.UnsupportedArgumentType
                        push!(unsuppargs, e.key)
                        num_unsupparg += 1
                        str *= "#ns arbcall\"" * s * "\"\n"
                    elseif e isa KeyError
                        push!(keyerrors, e.key)
                        num_keyerror += 1
                        str *= "#ni arbcall\"" * s * "\"\n"
                    else
                        rethrow(e)
                    end
                end
            end
        end
    end

    if verbose
        correctly_parsed = num_functions - num_unsupparg - num_keyerror
        @info "Correctly parsed functions: $correctly_parsed/$num_functions"
        @info "Unsupported types: $(collect(unsuppargs))"
        @info "Not implemented types: $(collect(keyerrors))"
    end

    return str
end

"""
    parse_and_generate_arbdoc(arb_doc_dir, out_dir = "src/arbcalls/")
Parses the Arb documentation and generates corresponding Julia files.
The value of `arb_doc_dir` should be a path to the directory
`doc/source/` in the Arb directory.
"""
function parse_and_generate_arbdoc(arb_doc_dir, out_dir = "src/arbcalls/")
    filenames = (
        "mag",
        "arf",
        "arb",
        "acb",
        "arb_poly",
        "acb_poly",
        "arb_mat",
        "acb_mat",
        "acb_hypgeom",
        "arb_hypgeom",
        "acb_elliptic",
        "acb_dirichlet",
    )

    unused_filenames = (
        "arb_fmpz_poly",
        "acb_dft",
        "acb_modular",
        "dirichlet",
        "bernoulli",
        "hypgeom",
        "partitions",
        "arb_calc",
        "acb_calc",
        "fmpz_extras",
        "bool_mat",
        "dlog",
        "fmpr",
    )

    manual_overrides = Dict{String,String}(
        "void mag_print(const mag_t x)" => "clashes with Base.print",
        "double arf_get_d(const arf_t x, arf_rnd_t rnd)" => "clashes with arf_get_si",
        "slong arf_get_si(const arf_t x, arf_rnd_t rnd)" => "clashes with arf_get_d",
        "void arf_print(const arf_t x)" => "clashes with Base.print",
        "arb_ptr _arb_vec_init(slong n)" => "clashes with similar method for acb",
        "double _arb_vec_estimate_allocated_bytes(slong len, slong prec)" =>
            "clashes with similar method for acb",
        "void arb_print(const arb_t x)" => "clashes with Base.print",
        "int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd)" =>
            "clashes with arb_can_round_mpfr",
        "int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd)" =>
            "clashes with arb_can_round_arf",
        "void arb_root(arb_t z, const arb_t x, ulong k, slong prec)" =>
            "alias to arb_root_ui",
        "acb_ptr _acb_vec_init(slong n)" => "clashes with similar method for arb",
        "double _acb_vec_estimate_allocated_bytes(slong len, slong prec)" =>
            "clashes with similar method for arb",
        "void acb_print(const acb_t x)" => "clashes with Base.print",
        "slong arb_poly_length(const arb_poly_t poly)" => "clashes with Base.length",
        "void _arb_poly_add(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)" =>
            "clashes with _arb_vec_add",
        "void _arb_poly_sub(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)" =>
            "clashes with _arb_vec_sub",
        "slong acb_poly_length(const acb_poly_t poly)" => "clashes with Base.length",
        "void _acb_poly_add(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)" =>
            "clashes with _acb_vec_add",
        "void _acb_poly_sub(acb_ptr C, acb_srcptr A, slong lenA, acb_srcptr B, slong lenB, slong prec)" =>
            "clashes with _acb_vec_sub",
        "void _acb_poly_elliptic_k_series(acb_ptr res, acb_srcptr z, slong zlen, slong len, slong prec)" =>
            "alias to _acb_elliptic_k_series",
        "void acb_poly_elliptic_k_series(acb_poly_t res, const acb_poly_t z, slong n, slong prec)" =>
            "alias to acb_elliptic_k_series",
        "void _acb_poly_elliptic_p_series(acb_ptr res, acb_srcptr z, slong zlen, const acb_t tau, slong len, slong prec)" =>
            "alias to _acb_elliptic_p_series",
        "void acb_poly_elliptic_p_series(acb_poly_t res, const acb_poly_t z, const acb_t tau, slong n, slong prec)" =>
            "alias to acb_elliptic_p_series",
    )

    for filename in filenames
        @info "Generating $filename.jl"
        title, sections = parse_arbdoc(joinpath(arb_doc_dir, "$filename.rst"))
        generate_file(
            joinpath(out_dir, "$filename.jl"),
            title,
            sections,
            verbose = true,
            manual_overrides = manual_overrides,
        )
        println("")
    end
end


text:
"""
@knet Op nce(q,s) computes the activation function for Noise
Contrastive Estimation.  Given q[K] and s[K,B] the output is p[K,B]
where p[i,j]=exp(s[i,j])/(exp(s[i,j])+q[i]).
"""
type NCE <: Op; NCE(;o...)=new(); end

ninputs(::NCE)=2
canoverwrite(::NCE)=true
back_reads_x(::NCE)=false
back_reads_y(::NCE)=true

# kq[K,1], s[K,B], p[i,j]=exp(s[i,j])/(exp(s[i,j])+kq[i])
nceforw!(kq::Array,s::Array,p::Array)=(for i=1:size2(s,1), j=1:size2(s,2); p[i,j]=exp(s[i,j])/(exp(s[i,j])+kq[i]); end; p)

@gpu nceforw!{T}(kq::CudaArray{T},s::CudaArray{T},p::CudaArray{T})=
    (T <: Float32 ? ccall((:nceforw32,libknet),Void,(Cint,Cint,Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}), size2(s,1),size2(s,2),kq,s,p) :
     T <: Float64 ? ccall((:nceforw64,libknet),Void,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),size2(s,1),size2(s,2),kq,s,p) :
     error("$T not supported"); gpusync(); p)

forw(::NCE, kq, s, p; o...)=(size2(s) != size2(p) ? throw(DimensionMismatch()) :
                             size2(s,1) != length(kq) ? throw(DimensionMismatch()) :
                             nceforw!(kq,s,p))

# ds[i,j] = dp[i,j]*p[i,j]*(1-p[i,j]) same as sigmback!
function back(::NCE, dp, dq, ds; y=nothing, o...) # y is p
    dq != nothing && (Base.warn_once("Taking gradient of constant"); fillsync!(dq,0))
    ds != nothing && (length(y)==length(ds)==length(dp)||throw(DimensionMismatch()); sigmback(y,dp,ds))
end

function infersize(a::NCE, q, s, p)
    q = (q == nothing ? [0,1] : length(q) == 2 ? [q...] : throw(DimensionMismatch("q=$q")))
    s = (s == nothing ? [0,0] : length(s) == 2 ? [s...] : throw(DimensionMismatch("s=$s")))
    p = (p == nothing ? [0,0] : length(p) == 2 ? [p...] : throw(DimensionMismatch("p=$p")))
    q[2] = (q[2]==1 ? 1 : q[2]==0 ? 1 : throw(DimensionMismatch("q=$q")))
    p1 = (q[1] > 0 ? q[1] : s[1] > 0 ? s[1] : p[1])
    p2 = (s[2] > 0 ? s[2] : p[2])
    if p1 > 0
        q[1]==p1 ? nothing : q[1]==0 ? q[1]=p1 : throw(DimensionMismatch())
        s[1]==p1 ? nothing : s[1]==0 ? s[1]=p1 : throw(DimensionMismatch())
        p[1]==p1 ? nothing : p[1]==0 ? p[1]=p1 : throw(DimensionMismatch())
    end
    if p2 > 0
        s[2]==p2 ? nothing : s[2]==0 ? s[2]=p2 : throw(DimensionMismatch())
        p[2]==p2 ? nothing : p[2]==0 ? p[2]=p2 : throw(DimensionMismatch())
    end
    return (tuple(q...), tuple(s...), tuple(p...))
end

### DEAD CODE:
# nce(q,s,p)=(NCE(),q,s,p)


text:
# poweriter.jl

using MIRT: poweriter

using Random: seed!
using LinearAlgebra: opnorm, norm
using Test: @test, @inferred

seed!(0)
M,N = 30,20

A = randn(M,N) # real
s0 = opnorm(A)
chat = false
_,s1 = @inferred poweriter(A ; tol=1e-9, chat=chat)
@test s0 ≈ s1

A = randn(ComplexF32, M, N) # complex
s0 = opnorm(A)
_,s1 = @inferred poweriter(A ; tol=1e-9, chat=chat)
@test s0 ≈ s1


text:
# cache such that creating svgs from LaTeX don't need to be created every time
# this is also used for test cases such that `tex2svg` doesn't need to be installed on Github Objects
include("latexsvgfile.jl")
latex(text::LaTeXString) = latex(text, O)
latex(text::LaTeXString, pos::Point) = latex(text, pos, :stroke)
latex(text::LaTeXString, x, y) = latex(text, Point(x, y), :stroke)

"""
    latex(text::LaTeXString, pos::Point, object::Symbol)

Add the latex string `text` to the top left corner of the LaTeX path.
Can be added to `Luxor.jl` graphics via [`Video`](@ref).

**NOTES:**
- **This only works if `tex2svg` is installed.**
    It can be installed using the following command (you may have to prefix this command with `sudo` depending on your installation):

        npm install -g mathjax-node-cli

- **The `latex` method must be called from within an [`Object`](@ref).**

# Arguments
- `text::LaTeXString`: a LaTeX string to render.
- `pos::Point`: position of the upper left corner of the latex text. Default: `O`
    - can be written as `x, y` instead of `Point(x, y)`
- `object::Symbol`: graphics objects defined by `Luxor.jl`. Default `:stroke`.
Available objects:
  - `:stroke` - Draws the latex string on the canvas. For more info check `Luxor.strokepath`
  - `:path` - Creates the path of the latex string but does not render it to the canvas.

# Throws
- `IOError`: mathjax-node-cli is not installed

# Example

```
using Javis
using LaTeXStrings

function ground(args...)
    background("white")
    sethue("black")
end

function draw_latex(video, object, frame)
    fontsize(50)
    x = 100
    y = 120
    latex(L"\\sqrt{5}", x, y)
end

demo = Video(500, 500)
Background(1:2, ground)
Object(draw_latex)
render(demo; pathname = "latex.gif")
```

"""
function latex(text::LaTeXString, pos::Point, draw_object::Symbol)
    object = CURRENT_OBJECT[1]
    opts = object.opts
    t = get(opts, :draw_text_t, 1.0)
    return animate_latex(text, pos, t, draw_object)
end

function animate_latex(text, pos::Point, t, object)
    svg = get_latex_svg(text)
    object == :stroke && (object = :fill)
    if t >= 1
        translate(pos)
        pathsvg(svg)
        do_action(object)
        translate(-pos)
        return
    end

    pathsvg(svg)
    polygon = pathtopoly()
    w, h = polywh(polygon)

    translate(pos)
    pathsvg(svg)
    do_action(:clip)
    r = t * sqrt(w^2 + h^2)
    circle(O, r, :fill)
    translate(-pos)
end

function get_latex_svg(text::LaTeXString)
    # check if it's cached
    if haskey(LaTeXSVG, text)
        svg = LaTeXSVG[text]
    else
        # remove the $
        ts = text.s[2:(end - 1)]
        command = `tex2svg $ts`
        try
            svg = read(command, String)
        catch e
            @warn "Using LaTeX needs the program `tex2svg` which might not be installed"
            @info "It can be installed using `npm install -g mathjax-node-cli`"
            throw(e)
        end
        LaTeXSVG[text] = svg
    end
    return svg
end


text:
"""
    corrplot

This plot type shows the correlation among input variables.
A correlation plot may be produced by a matrix.

A correlation matrix can also be created from the columns of a `DataFrame`
using the [`@df`](@ref) macro like so:

```julia
@df iris corrplot([:SepalLength :SepalWidth :PetalLength :PetalWidth])
```

The marker color in scatter plots reveals the degree of correlation.
Pass the desired colorgradient to `markercolor`.

With the default gradient positive correlations are blue, neutral are yellow
and negative are red. In the 2d-histograms, the color gradient shows the frequency
of points in that bin (as usual, controlled by `seriescolor`).
"""
@userplot CorrPlot

recipetype(::Val{:corrplot}, args...) = CorrPlot(args)

"""
    to_corrplot_matrix(mat)

Transforms the input into a correlation plot matrix.
Meant to be overloaded by other types!
"""
to_corrplot_matrix(x) = x

function update_ticks_guides(d::KW, labs, i, j, n)
    # d[:title]  = (i==1 ? _cycle(labs,j) : "")
    # d[:xticks] = (i==n)
    d[:xguide] = (i==n ? _cycle(labs,j) : "")
    # d[:yticks] = (j==1)
    d[:yguide] = (j==1 ? _cycle(labs,i) : "")
end

@recipe function f(cp::CorrPlot)
    mat = to_corrplot_matrix(cp.args[1])
    n = size(mat,2)
    C = cor(mat)
    labs = pop!(plotattributes, :label, [""])

    link := :x  # need custom linking for y
    layout := (n,n)
    legend := false
    foreground_color_border := nothing
    margin := 1mm
    titlefont := font(11)
    fillcolor --> Plots.fg_color(plotattributes)
    linecolor --> Plots.fg_color(plotattributes)
    markeralpha := 0.4
    grad = cgrad(get(plotattributes, :markercolor, :RdYlBu))
    indices = reshape(1:n^2, n, n)'
    title = get(plotattributes,:title,"")
    title_location = get(plotattributes, :title_location, :center)
    title := ""

    # histograms on the diagonal
    for i=1:n
        @series begin
            if title != "" && title_location == :left && i == 1
                title := title
            end
            seriestype := :histogram
            subplot := indices[i,i]
            grid := false
            xformatter --> ((i == n) ? :auto : (x -> ""))
            yformatter --> ((i == 1) ? :auto : (y -> ""))
            update_ticks_guides(plotattributes, labs, i, i, n)
            view(mat,:,i)
        end
    end

    # scatters
    for i=1:n
        ylink := setdiff(vec(indices[i,:]), indices[i,i])
        vi = view(mat,:,i)
        for j = 1:n
            j==i && continue
            vj = view(mat,:,j)
            subplot := indices[i,j]
            update_ticks_guides(plotattributes, labs, i, j, n)
            if i > j
                #below diag... scatter
                @series begin
                    seriestype := :scatter
                    markercolor := grad[0.5 + 0.5C[i,j]]
                    smooth := true
                    markerstrokewidth --> 0
                    xformatter --> ((i == n) ? :auto : (x -> ""))
                    yformatter --> ((j == 1) ? :auto : (y -> ""))
                    vj, vi
                end
            else
                #above diag... hist2d
                @series begin
                    seriestype := get(plotattributes, :seriestype, :histogram2d)
                    if title != "" && i == 1 && ((title_location == :center && j == div(n,2)+1) || (title_location == :right && j == n))
                        if iseven(n)
                            title_location := :left
                        end
                        title := title
                    end
                    xformatter --> ((i == n) ? :auto : (x -> ""))
                    yformatter --> ((j == 1) ? :auto : (y -> ""))
                    vj, vi
                end
            end
        end
    end
end


text:
function main()
    n = 600851475143

    i = 2
    while (n != 1)
        if (n % i == 0)
            n = n / i
        else
            i = i + 1
        end
    end

    println(i)
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
#####
##### other functions for computing Jacobians
#####

"""
    banded_function_evals(mat_fun::Fxn, U, Fargs...)

Computes block-banded matrix whose bands are entries of matrix-valued
function evals (e.g., a Jacobian). Returns SMatrix whose blocks correspond to function
components evaluated at values of U.

## Example:
```julia
julia> mat_fun(U) = [U[1] U[2]; U[2] U[1]]
julia> U = (randn(10),randn(10))
julia> banded_function_evals(mat_fun,U)
```
"""
function banded_function_evals(mat_fun::Fxn,U,Fargs...) where Fxn
    n = length(first(U))
    Nfields = length(U)
    A = SMatrix{Nfields,Nfields}([spzeros(n,n) for i = 1:Nfields, j = 1:Nfields])
    banded_function_evals!(A, mat_fun, U, Fargs...)
    return A
end

"""
    banded_function_evals!(A,mat_fun::Fxn, U, Fargs ...) where Fxn

Mutating version of [`banded_function_evals`](@ref).
"""
function banded_function_evals!(A,mat_fun::Fxn,U,Fargs...) where {Fxn}
    Nfields = size(A,2)
    num_pts = length(first(U))

    for i = 1:num_pts
        mat_i = mat_fun(getindex.(U,i),getindex.(Fargs,i)...)
        for n = 1:Nfields, m = 1:Nfields
            A[m,n][i,i] = mat_i[m,n] # TODO: replace with fast sparse constructor
        end
    end
end


text:
import Base: @pure

flatten(xs) = reduce(vcat, vec.(xs))

function indicesinvolved(expr)
    @match expr begin
        A_[idx__] => [(A => idx)]
        f_(ex__)  => reduce(vcat, [indicesinvolved(x) for x in ex])
        _ => []
    end
end

# Duplicate?
# """
# Type Domain `promote_op`
#
# The first argument is the type of the function
# """
# @generated function promote_op_t(::Type{F}, ::Type{T}) where {F,T}
#     :(_promote_op_t(F, $(T.parameters...)))
# end

# Dummy type to imitate a splatted function application
# This is passed into Inference.return_type (--> Compiler.return_type in v0.7)
struct _F{G}
    g::G
end
(f::_F{G})(x) where {G} = f.g(x...)

f(g) = x->g(x...)

# the equivalent to Base._promote_op
function _promote_op_t(F, @nospecialize T)
    G = Tuple{Base.Generator{Tuple{T},F}}
    Core.Compiler.return_type(first, G)
end

Zip2{I1,I2} = Base.Iterators.Zip{Tuple{I1,I2}}

function _promote_op_t(F, @nospecialize(R), @nospecialize(S))
    G = Tuple{Base.Generator{Zip2{Tuple{R},Tuple{S}},_F{F}}}
    Core.Compiler.return_type(first, G)
end

function _promote_op_t(F, @nospecialize(R), @nospecialize(S), @nospecialize T...)
    _promote_op_t(F, _promote_op_t(F, R, S), T...)
end

@pure function eltypes(::Type{T}) where T<:Tuple
    Tuple{map(eltype, T.parameters)...}
end

"""
Type Domain `promote_op`

The first argument is the type of the function
"""
@pure @generated function promote_op_t(::Type{F}, ::Type{T}) where {F,T}
    :(_promote_op_t(F, $(T.parameters...)))
end

"""
`promote_arraytype(F::Type, Ts::Type...)`

Returns an output array type for the result of applying a function of type `F`
on arrays of type `Ts`.
"""
@pure function promote_arraytype(::Type{F}, T...) where F
    length(T) == 1 && return T[1]
    promote_arraytype(F, promote_arraytype(F, T[1]), T[2:end]...)
end

@pure function promote_arraytype(::Type{F}, ::Type{T}, ::Type{S}) where {F, T<:Array, S<:Array}
    Array{_promote_op_t(F, eltype(T), eltype(S)), max(ndims(T), ndims(S))}
end

using SparseArrays

idxtype(::Type{SparseMatrixCSC{X,I}}) where {X,I} = I
@pure function promote_arraytype(::Type{F}, ::Type{T}, ::Type{S}) where {F, T<:SparseMatrixCSC, S<:SparseMatrixCSC}
    SparseMatrixCSC{_promote_op_t(F, eltype(T), eltype(S)),
                    promote_type(idxtype(T), idxtype(S))}
end

@pure function promote_arraytype(::Type{typeof(+)}, ::Type{T}, ::Type{S}) where {T<:Array, S<:SparseMatrixCSC}
    Array{_promote_op_t(F, eltype(T), eltype(S))}
end

@pure function promote_arraytype(::Type{typeof(*)}, ::Type{T}, ::Type{S}) where {T<:Array, S<:SparseMatrixCSC}
    SparseMatrixCSC{_promote_op_t(F, eltype(T), eltype(S)),
                    promote_type(idxtype(T), idxtype(S))}
end

@pure function promote_arraytype(::Type{typeof(+)}, ::Type{T}, ::Type{S}) where {T<:SparseMatrixCSC, S<:Array}
    Array{_promote_op_t(F, eltype(T), eltype(S))}
end

"""
`reduction_identity(f, T::Type)`

Identity value for reducing a collection of `T` with function `f`
"""
reduction_identity(f::Union{typeof(+), typeof(-)}, ::Type{T}) where {T} = zero(T)
reduction_identity(f::typeof(min), ::Type{T}) where {T} = typemax(T)
reduction_identity(f::typeof(max), ::Type{T}) where {T} = typemin(T)
reduction_identity(f::typeof(*), ::Type{T}) where {T} = one(T)
reduction_identity(f::typeof(push!), ::Type{T}) where {T} = T[]

function merge_dictofvecs(dicts...)
    merged_dict = Dict()
    for dict in dicts
        for (k, v) in dict
            if k in keys(merged_dict)
                merged_dict[k] = vcat(merged_dict[k], v)
            else
                merged_dict[k] = v
            end
        end
    end
    merged_dict
end


text:
#!/usr/bin/env julia

import Celeste: ParallelRun
import Celeste.SDSSIO: RunCamcolField

const datadir = joinpath(Pkg.dir("Celeste"), "test", "data")
wd = pwd()
cd(datadir)
run(`make RUN=3900 CAMCOL=2 FIELD=453`)
run(`make RUN=4382 CAMCOL=1 FIELD=70`)
cd(wd)

function benchmark_infer()
    # Warm up---this compiles the code
    box_compile = ParallelRun.BoundingBox(164.39, 164.40, 39.11, 39.13)
    field_triplets_compile = [RunCamcolField(3900, 6, 269),]
    ParallelRun.one_node_infer(field_triplets_compile, datadir; box=box_compile)

    # This box should only contain objid = 1237662224072638610, a faint
    # galaxy that overlaps with a really bright galaxy.
    # See http://skyserver.sdss.org/dr10/en/tools/chart/navi.aspx?ra=200.00958464398&dec=38.1880678852474
    # Update: Looks like the run time is all PSF fitting.
    # That's not worth fixing since PSF fitting is going away.
    # Let's update this script with another objid that's slow once we come across one.
    box = ParallelRun.BoundingBox(200.00957, 200.00959, 38.18806, 38.18807)
    field_triplets = [RunCamcolField(3900, 2, 453), RunCamcolField(4382, 1, 70)]
    @time ParallelRun.one_node_infer(field_triplets, datadir; box=box)

    Profile.init(10^6, 0.01)
    Profile.clear_malloc_data()
    # about half the run time is psf fitting, the other half is elbo evaluation
    @profile ParallelRun.one_node_infer(field_triplets, datadir; box=box)
    Profile.print(format=:flat, sortedby=:count)
end

benchmark_infer()


text:
using Lynx
using Luxor
using Colors
using Observables

include("utils.jl")

canvas  = Canvas(hexpand = true, vexpand = true)
name    = Observable("image.png")
image   = Image(asset"image.png")

# image = slice(image, 0, 0, 300, 300)
# image = image[0, 0, 300, 300]

mouse  = Mouse(canvas)
pos    = Point(200, 100)
zoom   = 1.0

zoom_btn = SpinButton(0:10:300, start = 100)

on(zoom_btn) do x
    global zoom = x / 100
end

function setimage!(path::String)
    global image;
    if !isempty(path)
        name[] = basename(path)
        image = Image(path)
    end
end

function updateimg(args...)
    path = Lynx.file!("Choose an image", filters = ["*.png"])
    setimage!(path)
end

rotr(args...) = rotl!(image)
rotl(args...) = rotr!(image)

onupdate(canvas) do dt, c
    background("#fff")
    origin(pos.x, pos.y)
    scale(zoom)
    drawimage(image)
end

onscroll(canvas) do e
    z = 10 * e.direction
    zoom_btn[] = clamp(zoom_btn[] + z, 10, 300)
end

onmousedrag(canvas) do event
    global pos += Point(mouse.dx, mouse.dy)
end

app = Window("Images", 900, 800,
    Box(:v,
        Paned(:v, position = 750,
            Paned(:v, position = 60,
                Box(:h, margin = 10, spacing = 5,
                    Button(updateimg, "insert-image", :Button, tooltip_text = "Select image"),
                    Button(rotl, "object-rotate-left", :Button, tooltip_text = "Rotate left"),
                    Button(rotr, "object-rotate-right", :Button, tooltip_text = "Rotate right"),
                ), # Box
                canvas,
            ), # Paned,
            Box(:h, spacing = 10, margin = 10,
                Label(name[], label = name, halign = Start),
                Box(:h, halign = End, hexpand = true, spacing = 10,
                    Label("Zoom "),
                    zoom_btn
                ) # Box
            ) # Box
        ) # Paned
    ) # Box
) # Window

@showall app; nothing


text:
# currently monster of a spaghetti code mess, multiple types interacting at the same
# time. Safest at this point in development is to just get this code running and
# will refactor once the dust has settled.
# current code is the result of several unit tests between IIF and RoME.jl
# a unified test to follow -- after which refactoring can start
function evalPotentialSpecific(
      fnc::T,
      Xi::Vector{Graphs.ExVertex},
      gwp::GenericWrapParam{T},
      solvefor::Int;
      N::Int=100  ) where {T <: FunctorPairwise}
  #
  # TODO -- enable partial constraints

  # TODO -- this part can be collapsed into common generic solver component, could be constructed and maintained at addFactor! time
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, N, solvefor)
  gwp.params = ARR
  gwp.varidx = sfidx
  gwp.measurement = gwp.samplerfnc(gwp.usrfnc!, maxlen)
  zDim = size(gwp.measurement[1],1) # TODO -- zDim aspect desperately needs to be redone
  if gwp.specialzDim
    zDim = gwp.usrfnc!.zDim[sfidx]
  end
  # gwp.zDim[sfidx] ??
  fr = FastRootGenericWrapParam{T}(gwp.params[sfidx], zDim, gwp)
  # and return complete fr/gwp
  # TODO -- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
  for n in 1:maxlen
    gwp.particleidx = n
    # gwp(x, res)
    numericRootGenericRandomizedFnc!( fr )
        # r = nlsolve( gwp, ARR[sfidx][:,gwp.particleidx] )
        # remember this is a deepcopy of original sfidx, since we are generating a proposal distribution
        # and not directly replacing the existing variable belief estimate
        # gwp.params[gwp.varidx][:,gwp.particleidx] = r.zero[:]
  end

  return gwp.params[gwp.varidx]
  # return evalPotential(typ, Xi, solvefor)
end

function evalPotentialSpecific(
      fnc::T,
      Xi::Vector{Graphs.ExVertex},
      gwp::GenericWrapParam{T},
      solvefor::Int;
      N::Int=100,
      spreadfactor::Float64=10.0  ) where {T <: FunctorPairwiseNH}
  #
  # TODO -- enable partial constraints

  # TODO -- this part can be collapsed into common generic solver component, could be constructed and maintained at addFactor! time
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, N, solvefor)
  gwp.params = ARR
  gwp.varidx = sfidx
  gwp.measurement = gwp.samplerfnc(gwp.usrfnc!, maxlen)
  zDim = size(gwp.measurement[1],1) # TODO -- zDim aspect desperately needs to be redone
  if gwp.specialzDim
    zDim = gwp.usrfnc!.zDim[sfidx]
  end
  fr = FastRootGenericWrapParam{T}(gwp.params[sfidx], zDim, gwp)
  # and return complete fr/gwp

  # nullhypothesis
  nhc = rand(gwp.usrfnc!.nullhypothesis, maxlen) - 1
  val = gwp.params[gwp.varidx]
  d = size(val,1)
  var = Base.var(val,2) + 1e-3
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*diagm(var[:]))

  # TODO --  Threads.@threads see area4 branch
  for n in 1:maxlen
    # gwp(x, res)
    if nhc[n] != 0
      gwp.particleidx = n
      numericRootGenericRandomizedFnc!( fr )
    else
      gwp.params[gwp.varidx][:,n] += rand(ENT)
    end
  end

  return gwp.params[gwp.varidx]
  # return evalPotential(typ, Xi, solvefor)
end

function evalPotentialSpecific(
      fnc::T,
      Xi::Vector{Graphs.ExVertex},
      gwp::GenericWrapParam{T},
      solvefor::Int;
      N::Int=100  ) where {T <: FunctorPairwiseMinimize}
  #
  # TODO -- this part can be collapsed into common generic solver component, could be constructed and maintained at addFactor! time
  ARR = Array{Array{Float64,2},1}()
  maxlen, sfidx = prepareparamsarray!(ARR, Xi, N, solvefor)
  gwp.params = ARR
  gwp.varidx = sfidx
  gwp.measurement = gwp.samplerfnc(gwp.usrfnc!, maxlen)
  zDim = size(gwp.measurement[1],1) # TODO -- zDim aspect desperately needs to be redone
  if gwp.specialzDim
    zDim = gwp.usrfnc!.zDim[sfidx]
  end

  fr = FastRootGenericWrapParam{T}(gwp.params[sfidx], zDim, gwp)

  # TODO -- once Threads.@threads have been optmized JuliaLang/julia#19967, also see area4 branch
  for n in 1:maxlen
    gwp.particleidx = n
    # gwp(x, res)
    # implement minimization here
    res = zeros(fr.xDim)
    gg = (x) -> fr.gwp(res, x)
    r = optimize(  gg, fr.X[1:fr.xDim,fr.gwp.particleidx] )
    # TODO -- clearly lots of optmization to be done here
    fr.Y[1:fr.xDim] = r.minimizer
    fr.X[:,fr.gwp.particleidx] = fr.Y

    # error("not implemented yet")
  end
  return gwp.params[gwp.varidx]
end

function evalPotentialSpecific(
      fnc::T,
      Xi::Vector{Graphs.ExVertex},
      generalwrapper::GenericWrapParam{T},
      solvefor::Int;
      N::Int=100  ) where {T <: FunctorSingleton}
  #
  generalwrapper.measurement = generalwrapper.samplerfnc(generalwrapper.usrfnc!, N)
  if !generalwrapper.partial
    return generalwrapper.measurement[1]
  else
    val = deepcopy(getVal(Xi[1]))
    i = 0
    for dimnum in fnc.partial
      i += 1
      val[dimnum,:] = generalwrapper.measurement[1][i,:]
    end
    return val
  end
end

function evalPotentialSpecific(
      fnc::T,
      Xi::Vector{Graphs.ExVertex},
      generalwrapper::GenericWrapParam{T},
      solvefor::Int;
      N::Int=100,
      spreadfactor::Float64=10.0  ) where {T <: FunctorSingletonNH}
  #

  # determine amount share of null hypothesis particles
  generalwrapper.measurement = generalwrapper.samplerfnc(generalwrapper.usrfnc!, N)
  # values of 0 imply null hypothesis
  # generalwrapper.usrfnc!.nullhypothesis::Distributions.Categorical
  nhc = rand(generalwrapper.usrfnc!.nullhypothesis, N) - 1

  val = getVal(Xi[1])
  d = size(val,1)
  var = Base.var(val,2) + 1e-3
  ENT = Distributions.MvNormal(zeros(d), spreadfactor*diagm(var[:]))

  for i in 1:N
    if nhc[i] == 0
      generalwrapper.measurement[1][:,i] = val[:,i] + rand(ENT)
    end
  end
  # TODO -- returning to memory location inside
  return generalwrapper.measurement[1]
end

# function evalPotentialSpecific{T <: FunctorPartialSingleton}(
#       fnc::T,
#       Xi::Vector{Graphs.ExVertex},
#       generalwrapper::GenericWrapParam{T},
#       solvefor::Int;
#       N::Int=100  )
#   #
#   generalwrapper.measurement = generalwrapper.samplerfnc(generalwrapper.usrfnc!, N)
#   return generalwrapper.measurement[1]
# end

# Multiple dispatch occurs internally, resulting in factor graph potential evaluations
function evalFactor2(fgl::FactorGraph, fct::Graphs.ExVertex, solvefor::Int; N::Int=100)
  # return evalPotential(fct.attributes["data"].fnc, solvefor) #evalPotential(fct.attributes["fnc"], solvefor)

  # TODO -- this build up of Xi is excessive and should be reduced
  # could happen at addFactor time
  Xi = Graphs.ExVertex[]
  for id in getData(fct).fncargvID
  # for id in fct.attributes["data"].fncargvID
    push!(Xi, getVert(fgl,id) ) # TODO localapi
    # push!(Xi, dlapi.getvertex(fgl,id))
  end
  fnctype = getData(fct).fnc
  # fnctype = fct.attributes["data"].fnc
  return evalPotentialSpecific(fnctype.usrfnc!, Xi, fnctype, solvefor, N=N)
  # return evalPotentialSpecific(modulefnc, Xi, fnctype, solvefor, N=N)
end

"""
    findRelatedFromPotential(fg, fctvert, varnodeid, N)

Compute proposal belief on varnodeid through fctvert representing some constraint in factor graph.
Always full dimension of variable node, where partial constraints will only influence directed
subset of variable dimensions. Remaining dimensions will keep existing variable values.
"""
function findRelatedFromPotential(fg::FactorGraph, idfct::Graphs.ExVertex, vertid::Int, N::Int) # vert
  # assuming it is properly initialized TODO
  ptsbw = evalFactor2(fg, idfct, vertid, N=N);
  # sum(abs(ptsbw)) < 1e-14 ? error("findRelatedFromPotential -- an input is zero") : nothing  # NOTE -- disable this validation test

  # TODO -- better to upsample before the projection
  Ndim = size(ptsbw,1)
  Npoints = size(ptsbw,2)
  # Assume we only have large particle population sizes, thanks to addNode!
  p = kde!(ptsbw, "lcv")
  if Npoints != N # this is where we control the overall particle set size
      p = resample(p,N)
  end
  return p
end


text:
using RobustAndOptimalControl, ControlSystems

# example 19.1
A = [
    -1 0 4
    0 2 0
    0 0 -3
]
B1 = [
    0 0
    1 0
    0 0
]
C1 = B1'
B2 = [
    1
    1
    1
]
C2 = B2'
D11 = 0I(2)
D12 = [
    0
    1
]
D21 = D12'
D22 = 0I(1)

P = ExtendedStateSpace(A, B1, B2, C1, C2, D11, D12, D21, D22)
for γ ∈ (nothing, 1000)
    K, Cl = h2synthesize(P, γ) # the example uses h2syn
    @test norm(Cl) ≈ 8.9582 rtol = 1e-3

    K_ML = let
        A = [-1.0 -4.236148004700622 4.0; -4.23607198996497 -6.472215758593602 -4.23607198996497; 0.0 -4.236148004700622 -3.0]
        B = [0.0; 4.23607198996497; 0.0]
        C = [0.0 -4.236148004700622 0.0]
        D = [0.0]
        sys = ss(A,B,C,D)
    end
    @test K ≈ K_ML rtol=1e-3

    Cl_ML = let
        A = [-1.0 0.0 4.0 0.0 -4.236148004700622 0.0; 0.0 2.0 0.0 0.0 -4.236148004700622 0.0; 0.0 0.0 -3.0 0.0 -4.236148004700622 0.0; 0.0 0.0 0.0 -1.0 -4.236148004700622 4.0; 4.23607198996497 4.23607198996497 4.23607198996497 -4.23607198996497 -6.472215758593602 -4.23607198996497; 0.0 0.0 0.0 0.0 -4.236148004700622 -3.0]
        B = [0.0 0.0; 1.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 4.23607198996497; 0.0 0.0]
        C = [0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 -4.236148004700622 0.0]
        D = [0.0 0.0; 0.0 0.0]
        ss(A,B,C,D)
    end

    @test Cl_ML ≈ Cl rtol=1e-3

    # s = tf("s")
    # K⁺ = ss(-148.79(s+1)*(s+3) / ((s+31.74)*(s+3.85)))
    # RobustAndOptimalControl.controller_reduction_weight(P, K⁺)
    # K̂⁺ = controller_reduction(P, K⁺, 1, true)
    # @test K̂⁺ ≈ tf(-117.085, [1, 34.526])
end

# Test conversion between H2 problem and LQG
K, Cl = h2synthesize(P)
#=
    The controller above is given by -L(sI-Ae)\K where L is feedback gain and K is Kalman gain.
    We have the following equivalences

    C1'C1    = Q
    D12'D12  = R

    B1*B1'   = R1
    D21*D21' = R2
=#
Q = C1'C1
R = D12'D12
R1 = B1*B1'
R2 = D21*D21'

L = lqr(system_mapping(P), Q, R)
@test L ≈ -K.C

Kal = kalman(system_mapping(P), R1, R2)
@test Kal ≈ K.B

l = LQGProblem(P)
@test lqr(l) == L
@test kalman(l) == Kal

# Implement normalized coprime fact from sec 13.8


text:
using Test
using SpikingNeuralNetworks
using ClearStacktrace
using Plots
unicodeplots()
SNN = SpikingNeuralNetworks
SNN.@load_units

for tests in [
    "spike_current_search_test.jl"
    ]
    include(tests)
end