text:
const PKG_DIR = abspath(@__DIR__, "..")
const GIT_DIR = joinpath(PKG_DIR, ".git")
const GIT_REV = try
    readchomp(`git --git-dir $GIT_DIR rev-parse --short HEAD`)
catch
    # Fallback to using the full SHA when git is not installed
    LibGit2.with(LibGit2.GitRepo(GIT_DIR)) do repo
        string(LibGit2.GitHash(LibGit2.GitObject(repo, "HEAD")))
    end
end

# Note: `kubectl apply` will only spawn a new job if there is a change to the job
# specification in the rendered manifest. If we simply used `GIT_REV` for as the image tag
# than any dirty changes may be ignored.
#
# Alternative solutions:
# - Use `imagePullPolicy: Always` (doesn't work with local images)
# - Introduce a unique label (modify the template)
# - Delete the old job first
# - Possibly switching to a pod for the manager instead of a job
const TAG = if !isempty(read(`git --git-dir $GIT_DIR status --short`))
    "$GIT_REV-dirty-$(getpid())"
else
    GIT_REV  # Re-runs here may just use the old job
end

const JOB_TEMPLATE = Mustache.load(joinpath(@__DIR__, "job.template.yaml"))
const TEST_IMAGE = get(ENV, "K8S_CLUSTER_MANAGERS_TEST_IMAGE", "k8s-cluster-managers:$TAG")

const POD_NAME_REGEX = r"Worker pod (?<worker_id>\d+): (?<pod_name>[a-z0-9.-]+)"

# Note: Regex should be generic enough to capture any stray output from the workers
const POD_OUTPUT_REGEX = r"From worker (?<worker_id>\d+):\s+(?<output>.*?)\r?\n"

# As a convenience we'll automatically build the Docker image when a user uses `Pkg.test()`.
# If the environmental variable is set we expect the Docker image has already been built.
if !haskey(ENV, "K8S_CLUSTER_MANAGERS_TEST_IMAGE")
    if success(`command -v minikube`) && !haskey(ENV, "MINIKUBE_ACTIVE_DOCKERD")
        @warn "minikube users should run `eval \$(minikube docker-env)` before executing " *
            "tests. Otherwise you may see pods fail with the reason \"ErrImageNeverPull\""
    end

    run(`docker build -t $TEST_IMAGE $PKG_DIR`)

    # Alternate build call which works on Apple Silicon
    # run(pipeline(`docker save $TEST_IMAGE`, `minikube ssh --native-ssh=false -- docker load`))
end

@testset "pod control" begin
    pod_control_manifest = YAML.load_file(joinpath(@__DIR__, "pod-control.yaml"))

    # Overwrite some parts of the specification
    # Note: We don't need to use the `TEST_IMAGE` here but it avoid having to download
    # another Docker image.
    pod_control_manifest["spec"]["containers"][1]["image"] = TEST_IMAGE
    pod_control_manifest["spec"]["containers"][1]["imagePullPolicy"] = "Never"

    @testset "named" begin
        manifest = deepcopy(pod_control_manifest)

        # Note: We do not want to use "generateName" here we want to use the same name when we
        # call `create_pod` multiple times. However, we do want to avoid conflicts with previous
        # `Pkg.test` executions.
        name = "test-pod-control-named-" * randsuffix()
        delete!(manifest["metadata"], "generateName")
        manifest["metadata"]["name"] = name

        @test_throws KubeError get_pod(name)
        @test_throws KubeError delete_pod(name)

        @info "Creating pod $name"
        create_pod(manifest)

        @test_throws KubeError create_pod(manifest)
        @test get_pod(name)["status"]["phase"] == "Pending"

        @info "Waiting for pod to start..."
        while get_pod(name)["status"]["phase"] == "Pending"
            sleep(1)
        end

        pod = get_pod(name)
        @test pod["status"]["phase"] == "Running"
        @test isempty(get(pod["metadata"], "labels", Dict()))

        label_pod(name, "testset" => "pod-control")

        pod = get_pod(name)
        @test get_pod(name)["metadata"]["labels"] == Dict("testset" => "pod-control")

        # Avoid deleting the pod if we've reached a unhandled phase. This allows for
        # investigation with `kubectl`.
        if pod["status"]["phase"] == "Running"
            @info "Deleting pod $name"
            # Avoid waiting for the pod to be deleted as this can take some time
            delete_pod(name, wait=false)

            # Note: Ideally we would be able to capture the "Terminating" status as reported by
            # `kubectl get pods -w`. Unforturnately I'm not sure how to retrieve this
            # information as it does not seem to be reported by `kubectl get pod` or
            # `kubectl get events`.
            reason = nothing
            while reason === nothing || reason == "Started"
                query = "involvedObject.name=$name"
                output = "jsonpath={.items[-1:].reason}"
                kubectl_cmd = `$(kubectl()) get events --field-selector $query -o=$output`
                reason = read(kubectl_cmd, String)
            end

            @test reason == "Killing"
        else
            @warn "Skipping deletion of pod $name"
        end

        # TODO: Would be nice to show details of the pod if any of these tests fail
    end

    @testset "generate name" begin
        manifest = deepcopy(pod_control_manifest)

        prefix = "test-pod-control-generate-name-"
        manifest["metadata"]["generateName"] = prefix

        name_a = create_pod(manifest)
        name_b = create_pod(manifest)

        @test name_a != name_b
        @test startswith(name_a, prefix)
        @test startswith(name_b, prefix)

        delete_pod(name_a, wait=false)
        delete_pod(name_b, wait=false)
    end

    @testset "exec" begin
        manifest = deepcopy(pod_control_manifest)
        manifest["metadata"]["generateName"] = "test-pod-control-exec-"

        name = create_pod(manifest)
        wait_for_running_pod(name; timeout=30)

        @test exec_pod(name, `bash -c "exit 0"`) === nothing

        try
            exec_pod(name, `kill -s INT 1`)
        catch e
            # Ensure that useful information from stdout/stderr is included
            if e isa KubeError
                @test occursin("\"kill\": executable file not found in \$PATH", e.msg)
                @test occursin("command terminated with exit code 126", e.msg)
            else
                rethrow()
            end
        end

        delete_pod(name; wait=false)
    end
end

let job_name = "test-success"
    @testset "$job_name" begin
        code = """
            using Distributed, K8sClusterManagers

            # Avoid trying to pull local-only image
            function configure(pod)
                pod["spec"]["containers"][1]["imagePullPolicy"] = "Never"
                return pod
            end
            pids = addprocs(K8sClusterManager(1; configure, pending_timeout=60, cpu="0.5", memory="300Mi"))

            println("Num Processes: ", nprocs())
            for i in workers()
                # Return the name of the pod via HOSTNAME
                println("Worker pod \$i: ", remotecall_fetch(() -> ENV["HOSTNAME"], i))
            end

            # Ensure that stdout/stderr on the worker is displayed on the manager
            @everywhere pids begin
                println(ENV["HOSTNAME"])
            end
            """

        command = ["julia"]
        args = ["-e", escape_yaml_string(code)]
        manifest = render(JOB_TEMPLATE; job_name, image=TEST_IMAGE, command, args)
        k8s_create(IOBuffer(manifest))

        # Wait for job to reach status: "Complete" or "Failed".
        #
        # There are a few scenarios where the job may become stuck:
        # - Insufficient cluster resources (pod stuck in the "Pending" status)
        # - Local Docker image does not exist (ErrImageNeverPull)
        @info "Waiting for $job_name job. This could take up to 4 minutes..."
        wait_job(job_name, condition=!isempty, timeout=4 * 60)

        manager_pod = first(pod_names("job-name" => job_name))
        worker_pod = first(pod_names("manager" => manager_pod))

        manager_log = pod_logs(manager_pod)
        call_matches = collect(eachmatch(POD_NAME_REGEX, manager_log))
        output_matches = collect(eachmatch(POD_OUTPUT_REGEX, manager_log))

        test_results = [
            @test get_job(job_name, jsonpath="{.status..type}") == "Complete"

            @test pod_exists(manager_pod)
            @test pod_exists(worker_pod)

            @test pod_phase(manager_pod) == "Succeeded"
            @test pod_phase(worker_pod) == "Succeeded"

            @test length(call_matches) == 1
            @test call_matches[1][:worker_id] == "2"
            @test call_matches[1][:pod_name] == worker_pod

            @test length(output_matches) == 1
            @test output_matches[1][:worker_id] == "2"
            @test output_matches[1][:output] == worker_pod

            # Ensure there are no unexpected error messages in the log
            @test !occursin(r"\bError\b"i, manager_log)
        ]

        # Display details to assist in debugging the failure
        if any(r -> !(r isa Test.Pass || r isa Test.Broken), test_results)
            report(job_name, "manager" => manager_pod, "worker" => worker_pod)
        end
    end
end

let job_name = "test-multi-addprocs"
    @testset "$job_name" begin
        code = """
            using Distributed, K8sClusterManagers

            # Avoid trying to pull local-only image
            function configure(pod)
                pod["spec"]["containers"][1]["imagePullPolicy"] = "Never"
                return pod
            end

            mgr = K8sClusterManager(1; configure, pending_timeout=60, cpu="0.5", memory="300Mi")
            addprocs(mgr)
            addprocs(mgr)

            println("Num Processes: ", nprocs())
            for i in workers()
                # Return the name of the pod via HOSTNAME
                println("Worker pod \$i: ", remotecall_fetch(() -> ENV["HOSTNAME"], i))
            end
            """

        command = ["julia"]
        args = ["-e", escape_yaml_string(code)]
        manifest = render(JOB_TEMPLATE; job_name, image=TEST_IMAGE, command, args)
        k8s_create(IOBuffer(manifest))

        # Wait for job to reach status: "Complete" or "Failed".
        @info "Waiting for $job_name job. This could take up to 4 minutes..."
        wait_job(job_name, condition=!isempty, timeout=4 * 60)

        manager_pod = first(pod_names("job-name" => job_name))
        worker_pods = pod_names("manager" => manager_pod; sort_by="{.metadata.labels.worker-id}")

        manager_log = pod_logs(manager_pod)
        matches = collect(eachmatch(POD_NAME_REGEX, manager_log))
        reported_worker_pods = map(m -> m[:pod_name], sort!(matches; by=m -> m[:worker_id]))

        test_results = [
            @test get_job(job_name, jsonpath="{.status..type}") == "Complete"

            @test pod_exists(manager_pod)
            @test length(worker_pods) == 2
            @test pod_exists(worker_pods[1])
            @test pod_exists(worker_pods[2])

            @test pod_phase(manager_pod) == "Succeeded"
            @test pod_phase(worker_pods[1]) == "Succeeded"
            @test pod_phase(worker_pods[2]) == "Succeeded"

            @test length(matches) == length(worker_pods)
            @test reported_worker_pods == worker_pods

            # Ensure there are no unexpected error messages in the log
            @test !occursin(r"\bError\b"i, manager_log)
        ]

        # Display details to assist in debugging the failure
        if any(r -> !(r isa Test.Pass || r isa Test.Broken), test_results)
            n = length(worker_pods)
            worker_pairs = map(enumerate(worker_pods)) do (i, w)
                "worker $i/$n" => w
            end

            report(job_name, "manager" => manager_pod, worker_pairs...)
        end
    end
end

let job_name = "test-interrupt"
    @testset "$job_name" begin
        code = """
            using Distributed, K8sClusterManagers

            # Avoid trying to pull local-only image
            function configure(pod)
                pod["spec"]["containers"][1]["imagePullPolicy"] = "Never"
                return pod
            end

            mgr = K8sClusterManager(1; configure, pending_timeout=60, cpu="0.5", memory="300Mi")
            pids = addprocs(mgr)
            interrupt(pids)
            """

        command = ["julia"]
        args = ["-e", escape_yaml_string(code)]
        manifest = render(JOB_TEMPLATE; job_name, image=TEST_IMAGE, command, args)
        k8s_create(IOBuffer(manifest))

        # Wait for job to reach status: "Complete" or "Failed".
        @info "Waiting for $job_name job. This could take up to 4 minutes..."
        wait_job(job_name, condition=!isempty, timeout=4 * 60)

        manager_pod = first(pod_names("job-name" => job_name))
        worker_pod = first(pod_names("manager" => manager_pod))

        manager_log = pod_logs(manager_pod)

        test_results = [
            @test pod_phase(manager_pod) == "Succeeded"
            @test pod_phase(worker_pod) == "Failed"

            # Ensure there are no unexpected error messages in the log
            @test length(collect(eachmatch(r"\bError\b"i, manager_log))) == 1
        ]

        # Display details to assist in debugging the failure
        if any(r -> !(r isa Test.Pass || r isa Test.Broken), test_results)
            report(job_name, "manager" => manager_pod, "worker" => worker_pod)
        end
    end
end

let job_name = "test-oom"
    @testset "$job_name" begin
        code = """
            using Distributed, K8sClusterManagers

            # Avoid trying to pull local-only image
            function configure(pod)
                pod["spec"]["containers"][1]["imagePullPolicy"] = "Never"
                return pod
            end

            mgr = K8sClusterManager(1; configure, pending_timeout=60, cpu="0.5", memory="300Mi")
            pids = addprocs(mgr)

            @everywhere begin
                function oom(T=Int64)
                    max_elements = Sys.total_memory() ÷ sizeof(T)
                    fill(zero(T), max_elements + 1)
                end
            end

            f = remotecall(oom, first(pids))

            # Use `wait` call to ensure that Distributed.jl notices the termination of the
            # worker. We can safely ignore the expected `ProcessExitedException`.
            try
                wait(f)
                error("`oom` call has not resulted in worker termination")
            catch e
                e isa ProcessExitedException || rethrow()
            end

            # Wait for the deregistration task which reports abnormal worker terminations.
            wait(K8sClusterManagers.DEREGISTER_ALERT)
            """

        command = ["julia"]
        args = ["-e", escape_yaml_string(code)]
        manifest = render(JOB_TEMPLATE; job_name, image=TEST_IMAGE, command, args)
        k8s_create(IOBuffer(manifest))

        # Wait for job to reach status: "Complete" or "Failed".
        @info "Waiting for $job_name job. This could take up to 4 minutes..."
        wait_job(job_name, condition=!isempty, timeout=4 * 60)

        manager_pod = first(pod_names("job-name" => job_name))
        worker_pod = first(pod_names("manager" => manager_pod))

        manager_log = pod_logs(manager_pod)
        worker_oom_msg = "Worker 2 on pod $worker_pod was terminated due to: OOMKilled"

        test_results = [
            @test pod_phase(manager_pod) == "Succeeded"
            @test pod_phase(worker_pod) == "Failed"

            @test pod_status(worker_pod) == ("terminated" => "OOMKilled")
            @test occursin(worker_oom_msg, manager_log)

            # Ensure there are no unexpected error messages in the log
            @test !occursin(r"\bError\b"i, manager_log)
        ]

        # Display details to assist in debugging the failure
        if any(r -> !(r isa Test.Pass || r isa Test.Broken), test_results)
            report(job_name, "manager" => manager_pod, "worker" => worker_pod)
        end
    end
end

let job_name = "test-pending-timeout"
    @testset "$job_name" begin
        code = """
            using Distributed, K8sClusterManagers

            # Avoid trying to pull local-only image
            function configure(pod)
                pod["spec"]["containers"][1]["imagePullPolicy"] = "Never"
                return pod
            end

            # Request 1 exbibyte of memory (should always fail)
            mgr = K8sClusterManager(1; configure, pending_timeout=1, memory="1Ei")
            pids = addprocs(mgr)

            println("Num Processes: ", nprocs())
            for i in pids
                # Return the name of the pod via HOSTNAME
                println("Worker pod \$i: ", remotecall_fetch(() -> ENV["HOSTNAME"], i))
            end
            """

        command = ["julia"]
        args = ["-e", escape_yaml_string(code)]
        manifest = render(JOB_TEMPLATE; job_name, image=TEST_IMAGE, command, args)
        k8s_create(IOBuffer(manifest))

        # Wait for job to reach status: "Complete" or "Failed".
        @info "Waiting for $job_name job. This could take up to 4 minutes..."
        wait_job(job_name, condition=!isempty, timeout=4 * 60)

        manager_pod = first(pod_names("job-name" => job_name))

        manager_log = pod_logs(manager_pod)
        call_matches = collect(eachmatch(POD_NAME_REGEX, manager_log))

        test_results = [
            @test get_job(job_name, jsonpath="{.status..type}") == "Complete"

            @test pod_exists(manager_pod)
            @test pod_phase(manager_pod) == "Succeeded"

            @test length(call_matches) == 0

            @test occursin(r"timed out after waiting for worker", manager_log)

            # Ensure there are no unexpected error messages in the log
            @test !occursin(r"\bError\b"i, manager_log)
        ]

        # Display details to assist in debugging the failure
        if any(r -> !(r isa Test.Pass || r isa Test.Broken), test_results)
            report(job_name, "manager" => manager_pod)
        end
    end
end


text:
###
# Reading an archive
#
# Basic outline for reading an archive in C:
#     1. Ask `archive_read_new` for an archive reader object.
#     2. Update any global properties as appropriate.
#        In particular, you'll certainly want to call appropriate
#        `archive_read_support_XXX` functions.
#     3. Call `archive_read_open_XXX` to open the archive
#     4. Repeatedly call `archive_read_next_header` to get information about
#        successive archive entries.  Call `archive_read_data` to extract
#        data for entries of interest.
#     5. Call `archive_read_free` to end processing.
#
# Basic outline for reading an archive in Julia:
#     1. Create a LibArchive.Reader
#     2. Update any global properties as appropriate.
#        In particular, you'll certainly want to call appropriate
#        LibArchive.support_XXX functions.
#     3. Read the entries (API TBD). LibArchive.jl will call libarchive
#        open functions if needed.
#     4. Call LibArchive.free to end processing.

abstract ReaderData

type Reader{T<:ReaderData} <: Archive
    data::T
    ptr::Ptr{Void}
    opened::Bool
    function Reader(data::T)
        ptr = ccall((:archive_read_new, libarchive), Ptr{Void}, ())
        ptr == C_NULL && throw(OutOfMemoryError())
        return Reader{T}(data, ptr, false)
    end
    # For pre-openned archive. (e.g. from LibALPM)
    function Reader(data::T, ptr::Ptr{Void}, opened::Bool)
        obj = new(data, ptr, opened)
        finalizer(obj, free)
        obj
    end
end

Reader{T<:ReaderData}(data::T) = Reader{T}(data)

function free(archive::Reader)
    ptr = archive.ptr
    ptr == C_NULL && return
    ccall((:archive_read_free, libarchive), Cint, (Ptr{Void},), ptr)
    archive.ptr = C_NULL
    nothing
end

Base.close(archive::Reader) =
    @_la_call(archive_read_close, (Ptr{Void},), archive)

function Base.cconvert(::Type{Ptr{Void}}, archive::Reader)
    archive.ptr == C_NULL && error("archive already freed")
    archive
end
Base.unsafe_convert(::Type{Ptr{Void}}, archive::Reader) = archive.ptr

function ensure_open(archive::Reader)
    archive.opened && return
    archive.opened = true
    try
        do_open(archive)
    catch
        archive.opened = false
        rethrow()
    end
    nothing
end

###
# Open filename

immutable ReadFileName{T} <: ReaderData
    name::T
    block_size::Int
end

Reader(fname::AbstractString; block_size=10240) =
    Reader(ReadFileName(fname, Int(block_size)::Int))
Reader(; block_size=10240) =
    Reader(ReadFileName(nothing, Int(block_size)::Int))

function do_open{T}(archive::Reader{ReadFileName{T}})
    data = archive.data
    @_la_call(archive_read_open_filename, (Ptr{Void}, Cstring, Csize_t),
              archive, data.name, data.block_size)
end

function do_open(archive::Reader{ReadFileName{Void}})
    data = archive.data
    @_la_call(archive_read_open_filename, (Ptr{Void}, Ptr{Void}, Csize_t),
              archive, C_NULL, data.block_size)
end

immutable ReadFD <: ReaderData
    fd::Cint
    block_size::Int
end

Reader(fd::Integer; block_size=10240) =
    Reader(ReadFD(Cint(fd), Int(block_size)::Int))

function do_open(archive::Reader{ReadFD})
    data = archive.data
    @_la_call(archive_read_open_fd, (Ptr{Void}, Cint, Csize_t),
              archive, data.fd, data.block_size)
end

# """
# Use this for reading multivolume files by filenames.
# NOTE: Must be NULL terminated. Sorting is NOT done.
# """
# int archive_read_open_filenames(struct archive*,
#                                 const char **_filenames, size_t block_size)

###
# Open memory

immutable ReadMemory{T,TO} <: ReaderData
    ptr::T
    obj::TO # Reference for GC
    size::Int
end

Reader{T<:Ptr,TO}(ptr::T, size, obj::TO=nothing) =
    Reader(ReadMemory{T,TO}(ptr, obj, size))
Reader(ary::Vector, size=sizeof(ary)) = Reader(pointer(ary), size, ary)

Base.unsafe_convert(::Type{Ptr{Void}}, mem_data::ReadMemory) = mem_data.ptr

function do_open{T<:ReadMemory}(archive::Reader{T})
    data = archive.data
    @_la_call(archive_read_open_memory, (Ptr{Void}, Ptr{Void}, Csize_t),
              archive, data, data.size)
end

###
# Generic reader

immutable GenericReadData{T} <: ReaderData
    data::T
    buff::Vector{UInt8}
end

Reader{T}(data::T, buff_size=10240) =
    Reader(GenericReadData{T}(data, Vector{UInt8}(buff_size)))

reader_open(archive::Reader, data) = nothing
function reader_readbytes end
function reader_skip end
function reader_seek end
reader_close(archive::Reader, data) = nothing

function reader_open_callback(c_archive, archive)
    try
        clear_error(archive)
        reader_open(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

function reader_read_callback(c_archive, archive, buff::Ptr{Ptr{Void}})
    try
        clear_error(archive)
        bytes_read = reader_readbytes(archive, archive.data.data,
                                      archive.data.buff)
        unsafe_store!(buff, pointer(archive.data.buff))
        return Cssize_t(bytes_read)
    catch ex
        set_exception(archive, ex)
        return Cssize_t(0)
    end
end

function reader_skip_callback(c_archive, archive, request)
    try
        clear_error(archive)
        return Int64(reader_skip(archive, archive.data.data, request))
    catch ex
        set_exception(archive, ex)
        return Int64(0)
    end
end

function reader_seek_callback(c_archive, archive, request, whence)
    try
        clear_error(archive)
        return Int64(reader_seek(archive, archive.data.data, request, whence))
    catch ex
        set_exception(archive, ex)
        return Int64(0)
    end
end

function reader_close_callback(c_archive, archive)
    try
        clear_error(archive)
        reader_close(archive, archive.data.data)
        return Libc.errno(archive) == 0 ? Cint(0) : Status.WARN
    catch ex
        return set_exception(archive, ex)
    end
end

reader_readbytes(archive, io::IO, buff) = readbytes!(io, buff)
reader_skip(archive, io::IO, sz) =
    (skip(io, sz); sz)

function do_open{T<:GenericReadData}(archive::Reader{T})
    # Set various callbacks
    @_la_call(archive_read_set_callback_data,
              (Ptr{Void}, Any), archive, archive)
    @_la_call(archive_read_set_open_callback,
              (Ptr{Void}, Ptr{Void}), archive,
              to_open_callback(reader_open_callback, Reader{T}))
    @_la_call(archive_read_set_read_callback,
              (Ptr{Void}, Ptr{Void}), archive,
              to_read_callback(reader_read_callback, Reader{T}))
    @_la_call(archive_read_set_seek_callback,
              (Ptr{Void}, Ptr{Void}), archive,
              to_seek_callback(reader_seek_callback, Reader{T}))
    @_la_call(archive_read_set_skip_callback,
              (Ptr{Void}, Ptr{Void}), archive,
              to_skip_callback(reader_skip_callback, Reader{T}))
    @_la_call(archive_read_set_close_callback,
              (Ptr{Void}, Ptr{Void}), archive,
              to_close_callback(reader_close_callback, Reader{T}))
    @_la_call(archive_read_open1, (Ptr{Void},), archive)
end

function next_header(archive::Reader)
    ensure_open(archive)
    entry = Entry(archive)
    @_la_call(archive_read_next_header2, (Ptr{Void}, Ptr{Void}),
              archive, entry)
    entry
end

"""
Retrieve the byte offset in UNCOMPRESSED data where last-read
header started.
"""
header_position(archive::Reader) =
    ccall((:archive_read_header_position, libarchive),
          Int64, (Ptr{Void},), archive)

"Read data from the body of an entry.  Similar to read(2)."
@inline function unsafe_archive_read(archive::Reader, ptr::Ptr{UInt8}, sz::UInt)
    nb = ccall((:archive_read_data, libarchive), Cssize_t,
               (Ptr{Void}, Ptr{Void}, Csize_t), archive, ptr, sz)
    nb < 0 && _la_error(Cint(nb), archive)
    nb % Csize_t
end

"Read data from the body of an entry.  Similar to read(2)."
function unsafe_read(archive::Reader, ptr::Ptr{UInt8}, sz::UInt)
    unsafe_archive_read(archive, ptr, sz) != sz && throw(EOFError())
    nothing
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.read(archive::Reader, ::Type{UInt8})
    b = Ref{UInt8}()
    unsafe_read(archive, Base.unsafe_convert(Ptr{UInt8}, b), UInt(1))
    b[]
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.readbytes!(archive::Reader, b::Array{UInt8}, nb=length(b))
    nbread = unsafe_archive_read(archive, Base.unsafe_convert(Ptr{UInt8}, b),
                                 UInt(nb))
    nbread < nb && resize!(b, nbread)
    nbread
end

"Read data from the body of an entry.  Similar to read(2)."
function Base.read(archive::Reader)
    # Maybe we can optimize this with block read?
    block_size = UInt(1024)
    res = Vector{UInt8}(block_size)
    pos = 1
    while true
        ptr = Base.unsafe_convert(Ptr{UInt8}, Ref(res, pos))
        nbread = unsafe_archive_read(archive, ptr, block_size)
        if nbread < block_size
            resize!(res, nbread + pos - 1)
            return res
        end
        pos = Int(pos + block_size)
        resize!(res, pos + block_size - 1)
    end
end

"Read data from the body of an entry.  Similar to read(2)."
@inline Base.readavailable(archive::Reader) = read(archive)

"Seek within the body of an entry.  Similar to lseek(2)."
Base.seek(archive::Reader, offset, what) =
    ccall((:archive_seek_data, libarchive), Int64,
          (Ptr{Void}, Int64, Cint), archive, offset, what)

"Skips entire entry"
Base.skip(archive::Reader) =
    @_la_call(archive_read_data_skip, (Ptr{Void},), archive)

"Writes data to specified filedes"
read_into_fd(archive::Reader, fd::Integer) =
    @_la_call(archive_read_data_into_fd, (Ptr{Void}, Cint), archive, fd)

Reader(f::Function) = archive_guard(f, Reader())
Reader(f::Function, args...; kws...) =
    archive_guard(f, Reader(args...; kws...))

@deprecate file_reader(block_size=10240) Reader(block_size=block_size)
@deprecate file_reader(fname::AbstractString,
                       block_size=10240) Reader(fname, block_size=block_size)
@deprecate file_reader(fd::Integer,
                       block_size=10240) Reader(fd, block_size=block_size)
@deprecate mem_reader{T<:Vector}(data::T, size=sizeof(data)) Reader(data, size)
@deprecate mem_reader(obj, size=sizeof(obj)) Reader(pointer(obj), size, obj)
@deprecate gen_reader(data, buff_size=10240) Reader(data, buff_size)

# /*
#  * Set read options.
#  */
# /* Apply option to the format only. */
# int archive_read_set_format_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option to the filter only. */
# int archive_read_set_filter_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option to both the format and the filter. */
# int archive_read_set_option(struct archive *_a,
# 			    const char *m, const char *o,
# 			    const char *v);
# /* Apply option string to both the format and the filter. */
# int archive_read_set_options(struct archive *_a,
# 			    const char *opts);

# /*
#  * A zero-copy version of archive_read_data that also exposes the file offset
#  * of each returned block.  Note that the client has no way to specify
#  * the desired size of the block.  The API does guarantee that offsets will
#  * be strictly increasing and that returned blocks will not overlap.
#  */
# int archive_read_data_block(struct archive *a,
#                    const void **buff, size_t *size, int64_t *offset);

# /*-
#  * Convenience function to recreate the current entry (whose header
#  * has just been read) on disk.
#  *
#  * This does quite a bit more than just copy data to disk. It also:
#  *  - Creates intermediate directories as required.
#  *  - Manages directory permissions:  non-writable directories will
#  *    be initially created with write permission enabled; when the
#  *    archive is closed, dir permissions are edited to the values specified
#  *    in the archive.
#  *  - Checks hardlinks:  hardlinks will not be extracted unless the
#  *    linked-to file was also extracted within the same session. (TODO)
#  */

# /* The "flags" argument selects optional behavior, 'OR' the flags you want. */

# int archive_read_extract(struct archive *, struct archive_entry *,
# 		     int flags);
# int archive_read_extract2(struct archive *, struct archive_entry *,
# 		     struct archive * /* dest */);
# void	 archive_read_extract_set_progress_callback(struct archive *,
#      void (*_progress_func)(void *), void *_user_data);

# /* Record the dev/ino of a file that will not be written.  This is
#  * generally set to the dev/ino of the archive being read. */
# void		archive_read_extract_set_skip_file(struct archive *,
#      int64_t, int64_t);

ext: jl


text:
using Images, ImageFiltering, LinearAlgebra

"""
    gradient_magnitude(f)

Return the directional derivatives (components of the gradient)
and the gradient mag for both color and grayscale images. For
`spatialorder(f)=="xy"`, `grad[1]` is a derivative wrt x, the first
coordinate and `grad[2]` wrt to the second...

"""
function gradient_magnitude(f)
    @assert(length(size(f))==2) # only 2D supported for the moment
    s=ImageFiltering.Kernel.sobel()     # vytvorime Sobelovsky filtr ve smeru x a y
    grad=[Images.imfilter(f,s[i]) for i in 1:2]   # filtrujeme ve smeru x a y
    return grad, Images.magnitude(grad[1],grad[2]) # vracime gradient a jeho normu
end

"""
    find_edge_points(img::Image; spacing::Int=35, point_count::Int=35, sigma=1, mask::BitArray{2}=BitArray{2}(undef, 0,0))

Greadily locate `point_count` of points in `img` with the highest gradient magnitude. These points have will be atleast `spacing`
pixels appart, they are returned as an array of `CartesianIndex`. A `mask` can be used to exclude some parts of `img` from the search.

Note: `sigma` is the gaussian filter used to smooth `img` before looking for gradients.
"""
function find_edge_points(img::Image;
                          spacing::Int=10,
                          point_count::Int=500,
                          sigma=1,
                          mask::BitArray{2}=BitArray{2}(undef, 0,0),
                          debug::Bool=false)

    # calculate gradient mag
    grad, mag = gradient_magnitude(ImageFiltering.imfilter(img, ImageFiltering.Kernel.gaussian(sigma)))

    @assert typeof(mag) <: Array

    #const
    marker = zero(eltype(mag)) # marker for masked gradients we don't want keypoints at the points of 0 gradient
    # marker = 0.25
    if !isempty(mask)
        @assert( size(mag) == size(mask))
        for i = 1:length(mag)
            if mask[i] == 0
                mag[i] = marker
            end
        end
    end

    N = 2
    permutation_array = sortperm(reshape(mag,length(mag)), alg=PartialQuickSort(2*point_count*spacing^(N-1)),rev=true) # check here

    image_inds = CartesianIndices(size(mag))

    numkpts = 0 # point_count of keypoints found so far
    kpts = Array{CartesianIndex{2}, 1}(undef, point_count)

    # go through indices from the largest gradient
    for perm in permutation_array

        if numkpts >= point_count break; end
        if mag[perm] == marker continue; end

        # @assert(mag[perm] > 10 * eps(Float64))

        # good keypoint
        numkpts += 1

        ind = image_inds[perm]
        ppos = Tuple(ind) # position in pixel units

        kpts[numkpts] = ind

        # prevent the neighbors from being added
        fill_circle!(mag, ppos, spacing, marker)
    end # for

    if numkpts < point_count # shrink array if less keypoints were found
        resize!(kpts, numkpts)
    end
    if debug
        return kpts, mag
    end
    return kpts
end # function find_keypoints_from_gradients

"""
    fill_circle!(mag, pos, spacing, marker)

Fill a circle with a radius of `spacing` around the point `pos` into the matrix `mag` with `marker`.
"""
function fill_circle!(mag, pos, spacing, marker)
    mag_size = size(mag)
    for k = (pos[1]-spacing <= 1 ? 1 : pos[1]-spacing) : (pos[1]-spacing >= mag_size[1] ? pos[1]-spacing : mag_size[1]),
        l = (pos[2]-spacing <= 1 ? 1 : pos[2]-spacing) : (pos[2]-spacing >= mag_size[2] ? pos[2]-spacing : mag_size[2])
        # @assert (k >= 1) k
        if (sum(((k, l) .- pos).^2) <= spacing^2)
            mag[k, l] = marker
        end
    end
end

# function find_keypoints_from_gradients_p_field(f::Array{T, 2}; spacing=25, point_count::Int=100, sigma=ones(Float64,2), nlength=1., mask=[]) where {T}
#
#     # calculate gradient mag
#     grad, mag = gradient_magnitude(ImageFiltering.imfilter(f, ImageFiltering.Kernel.gaussian(sigma)))
#
#     mag_save = copy(mag)
#     mag = mag.^4
#
#     @assert typeof(mag) <: Array
#
#     #const
#     marker = zero(eltype(mag)) # marker for masked gradients we don't want keypoints at the points of 0 gradient
#     if mask != []
#         @assert( size(mag) == size(mask))
#         for i = 1:length(mag)
#             if mask[i] == 0
#                 mag[i] = marker
#             end
#         end
#     end
#     # sort indices by decreasing mag
#
#     kpts = Array{SimpleKeypoint}(undef, point_count)
#     # go through indices from the largest gradient and put each
#     for k in 1:point_count
#
#         ind = wsample(Array(1:length(mag)), mag[:])
#
#         # @assert(mag[ind] > 10 * eps(Float64))
#         # good keypoint, let"s add it
#
#         ppos = Tuple(CartesianIndices(size(mag))[ind]) # position in pixel units
#
#         N = 2
#         norm = LinearAlgebra.norm(Float64[ grad[l][ind] for l=1:N], 2) * nlength
#         @assert mag_save[ind] ≈ norm
#         kpts[k] = SimpleKeypoint(ppos, norm)
#         # prevent the neighbors from being added
#         fill_circle!(mag, ppos, spacing, marker)
#     end # for i
#
#     return kpts, mag_save, mag
# end # function find_keypoints_from_gradients


text:
module VerilogFIR

export generate_fir

# max and min values for n bit number in 2s complement
maxrep(n) = 2^(n-1)-1
minrep(n) = -2^(n-1)

fminabs(x...) = floor(Int64, minimum(abs.(x)))
clog2(x) = ceil(Int64, log2(x))
flog2(x) = floor(Int64, log2(x))

nbits_u(x) = clog2(x+1)
nbits_s(x) = x > 0 ? clog2(x+1)+1 : clog2(-x)+1

tab(n) = " "^(4n)

macro makedict(args...)
    return esc(:(Dict($([:($(QuoteNode(i)) => $i) for i in args]...))))
end

################################################################################
# Top level function
################################################################################
generate_fir(args...; kwargs...) = generate_fir(stdout, args...; kwargs...)

function generate_fir(file::String, args...; kwargs...)
    io = open(file, "w")
    generate_fir(io, args...; kwargs...)
    close(io)
end

function generate_fir(
        io::IO, coeffs;
        w_input   = 8,
        w_output  = 8,
        w_coeff   = 8,
    )

    # Scale coefficients
    scaled_coeffs, bitshift = scale_coefficients(coeffs, w_coeff)

    # Compute the number of taps
    taps  = length(scaled_coeffs)
    w_acc = acc_width(scaled_coeffs, w_input)

    d = @makedict w_input w_output w_coeff scaled_coeffs bitshift taps w_acc

    # Generate Verilog header
    header(io, d)
    declarations(io, d)
    initialize(io, d)
    register_input(io, d)
    gen_filter(io, d)
    output(io, d)
    endmodule(io, d)
end

@doc """
    generate_fir(io::Union{IO,String}, coeffs; w_input, w_output, w_coeff)

Generate Verilog for a multicycle FIR filter described by `coeffs`. If `io` is
an `IO`, output will be written directly to `io`. Otherwise, if `io` is a
`String`, a file will be created with the name of `io` and output will be
written to that file.

# Key word arguments
* `w_input`  - Bitwidth of the input. Default: `8`.
* `w_output` - Bitwidth of the output. Default: `8`.
* `w_coeff`  - Bitwidth of the coefficients. Default: `8`.
""" generate_fir

################################################################################
# Analysis Functions
################################################################################

function scale_coefficients(coeffs, w_coeff)
    minval, maxval = extrema(coeffs)
    # Determine scale factor
    pos_scale = maxrep(w_coeff)/maxval
    neg_scale = minrep(w_coeff)/minval
    # Get the maximum number we can multiply coefficients by to get while still
    # remaining within the bit width designated for the coefficients.
    maxfactor = fminabs(pos_scale, neg_scale)
    # Convert to the largest power of 2 less than maxfactor
    factor = 2^flog2(maxfactor)
    return round.(Int64, factor * coeffs), flog2(factor)
end

function acc_width(scaled_coeffs, w_input)
    minval = minrep(w_input)
    maxval = maxrep(w_input)
    # Compute the maximum and minimum values for the output of this filter
    maxoutput = 0
    minoutput = 0
    for c in scaled_coeffs
        if c > 0
            maxoutput += maxval*c
            minoutput += minval*c
        elseif c < 0
            maxoutput += minval*c
            minoutput += maxval*c
        end
    end
    bits_for_max = nbits_s(maxoutput)
    bits_for_min = nbits_s(minoutput)
    return max(bits_for_max, bits_for_min)
end

################################################################################
# Printing functions
################################################################################
function header(io, d)
    print(io,"""
        module filter(
                // Interface signals
                input clk,
                input start,
                output reg done,
                // Data Signals
                input  [$(d[:w_input]-1):0] data_in,
                output reg [$(d[:w_output]-1):0] data_out
            );
        """)
end

function declarations(io, d)
    print(io,"""
            // Coefficient Storage
            reg signed [$(d[:w_coeff]-1):0] coeff [$(d[:taps]-1):0];
            reg signed [$(d[:w_input]-1):0] data  [$(d[:taps]-1):0];
            // Counter for iterating through coefficients.
            reg [$(clog2(d[:taps])-1):0] count;
            // Accumulator
            reg signed [$(d[:w_acc]-1):0] acc;

            // State machine signals
            localparam IDLE = 0;
            localparam RUN  = 1;

            reg state;
        """)
end

function initialize(io, d)
    # Print initial statement
    print(io, "\n")
    print(io, tab(1), "initial begin\n")
    for (i,c) in enumerate(d[:scaled_coeffs])
        print(io, tab(2), "coeff[$(i-1)] = $c;\n")
    end
    # end initial block
    print(io, tab(1), "end\n")
end

function register_input(io, d)
    print(io, """
            always @(posedge clk) begin : capture
                integer i;
                if (start) begin
                for (i = 0; i < $(d[:taps]-1) ; i = i+1) begin
                        data[i+1] <= data[i];
                    end
                    data[0] <= data_in;
                end
            end
        """)
end

function gen_filter(io, d)
    print(io, """
            always @(posedge clk) begin
                case (state)
                    IDLE: begin
                        done <= 1'b0;
                        if (start) begin
                            count <= $(d[:taps]-1);
                            acc   <= 0;
                            state <= RUN;
                        end
                    end

                    RUN: begin
                        count <= count - 1'b1;
                        acc   <= acc + data[count] * coeff[count];
                        if (count == 0) begin
                            state <= IDLE;
                            done  <= 1'b1;
                        end
                    end
                endcase
            end
        """)
end

function output(io, d)
    output_left = d[:bitshift] + d[:w_output] - 1
    # Determine if saturation detection is needed for the output
    saturate = output_left < d[:w_acc] - 1

    if saturate
        print(io, """
                always @(posedge clk) begin
                    if (done) begin
                        // Saturate if necessary
                        if (acc >= 2 ** $(output_left)) begin
                            data_out <= $(maxrep(d[:w_output]));
                        end else if (acc < -(2 ** $(output_left))) begin
                            data_out <= $(minrep(d[:w_output]));
                        end else begin
                            data_out <= acc[$output_left:$(d[:bitshift])];
                        end
                    end
                end
            """)
    else
        print(io, """
                always @(posedge clk) begin
                    if (done) begin
                        data_out <= acc[$output_left:$(d[:bitshift])];
                    end
                end
            """)
    end
end

function endmodule(io, d)
    print(io, "endmodule")
end
end # module


text:
function AESCBC(plaintext, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache)
	len = length(plaintext)
	pad = 16 - (len % 16)
	ciphertext = similar(Array{UInt8, 1}, pad + len)
	AESCBC!(ciphertext, plaintext, iv, key, cache)
end

function AESCBC!(ciphertext, plaintext, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache)
	len = length(plaintext)
	pad = 16 - (len % 16)
	for i in 1:len
		ciphertext[i] = UInt8(plaintext[i])
	end
	for i in 1:pad
		ciphertext[len+i] = pad
	end
	iters = Int((len + pad) / 16)
	ciphertextblock = AESBlock(ciphertext)
	prevblock = AESBlock(ciphertext)
	# iteration 1 start
	ciphertextblock .⊻= iv
	AESEncryptBlock!(ciphertextblock, ciphertextblock, key, cache)
	increment!(ciphertextblock)
	# iteration 1 end
	for i in 2:iters
		@. ciphertextblock ⊻= prevblock
		AESEncryptBlock!(ciphertextblock, ciphertextblock, key, cache)
		increment!(ciphertextblock)
		increment!(prevblock)
	end
	ciphertext
end

function AESCBC_D(ciphertext::Array{UInt8, 1}, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache; remove_pad=true)
	len = length(ciphertext)
	iters = Int(len / 16)
	plaintext = similar(Array{UInt8, 1}, len)
	AESCBC_D!(plaintext, ciphertext, iv, key, cache; remove_pad=remove_pad)
end

function AESCBC_D!(plaintext, ciphertext::Array{UInt8, 1}, iv::Array{UInt8, 1}, key::AbstractAESKey, cache::AbstractAESCache; remove_pad=true)
	len = length(ciphertext)
	iters = Int(len / 16)
	plaintextblock = AESBlock(plaintext)
	ciphertextblock = AESBlock(ciphertext)
	prevblock = AESBlock(ciphertext)
	# iteration 1 start
	AESDecryptBlock!(plaintextblock, ciphertextblock, key, cache)
	plaintextblock .⊻= iv
	increment!(plaintextblock)
	increment!(ciphertextblock)
	# iteration 1 end
	for i in 2:iters
		AESDecryptBlock!(plaintextblock, ciphertextblock, key, cache)
		plaintextblock .⊻= prevblock
		increment!(plaintextblock)
		increment!(ciphertextblock)
		increment!(prevblock)
	end
	if remove_pad
		pad = plaintext[end]
		@view(plaintext[1:len-pad])
	else
		plaintext
	end
end


text:
# Function to check the eigen values

function checkEigenValues(EigenValues, elementNo, elemcount)

  counter = 0
  for i = 1:length(EigenValues)
    if isreal(EigenValues[i]) ==  true
      if EigenValues[i] < 0 && abs(EigenValues[i]) > 1e-12
        counter += 1
        elemcount += 1
      end
    else
      # Check if the real and imaginary parts are positive
      println
      if real(EigenValues[i]) < 0 && abs(real(EigenValues[i])) > 1e-12
        counter += 1
        elemcount += 1
      elseif imag(EigenValues[i]) < 0 && abs(imag(EigenValues[i])) > 1e-12
        counter += 1
        elemcount += 1
      end
    end  # End if isreal(EigenValues[i]) ==  true
  end # End for i = 1:length(EigenValues)
  if counter > 0
    println("Element ", elementNo, " has negative eigen values")
  end

  return nothing
end

function elementEigenValues(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh,Tsol, Tres, Tdim}

  res_0 = zeros(eqn.res_vec)
  res_0_norm = physicsRhs(mesh, sbp, eqn, opts, res_0, (evalResidual,))
  pert = 1e-6
  numDofPerElement = mesh.numDofPerNode*mesh.numNodesPerElement
  elem_eigen_vals = zeros(numDofPerElement)
  elem_count = 0 # Count the number of elements with negative eigen values
  #=
  for i = 1:mesh.numEl
    elem_jac = zeros(numDofPerElement,numDofPerElement)
    elem_res_vec = zeros(numDofPerElement) # Element wise residual
    elem_q_vec = zeros(elem_res_vec)  # 1D array containing element wise conservative variables
    elemRes0 = zeros(elem_res_vec) # Collects the original residual
    entry_orig = zero(eltype(eqn.q_vec)) # stores one element of original res_vec) at a time
    orig_res_vec = copy(eqn.res_vec)

    orig_q_vec = copy(eqn.res_vec)
    for j = 1:mesh.numNodesPerElement # Populating elemRes0 & elem_q_vec
      for k = 1:mesh.numDofPerNode
        elemRes0[(j-1)*mesh.numDofPerNode + k] = eqn.res[k,j,i]
        elem_q_vec[(j-1)*mesh.numDofPerNode + k] = eqn.q[k,j,i]
      end
    end

    # println("elemRes0 - eqn.res_vec = \n", elemRes0 - eqn.res_vec)
    # println("elem_q_vec - eqn.q_vec = \n", elem_q_vec - eqn.q_vec)

    for j = 1:numDofPerElement # Perturb all dofs in the element
      if j == 1
        entry_orig = elem_q_vec[j]
        elem_q_vec[j] += pert  # Perturb the jth element of q
      else
        elem_q_vec[j-1] = entry_orig
        entry_orig = elem_q_vec[j]
        elem_q_vec[j] += pert
      end
      #println("elem_q_vec - eqn.q_vec = \n", elem_q_vec - eqn.q_vec)
      # Disassembling q_vec to q at the element level
      for l = 1:mesh.numNodesPerElement
        for k = 1:mesh.numDofPerNode
          eqn.q[k,l,i] = elem_q_vec[(l-1)*mesh.numDofPerNode + k]
        end
      end
      evalResidual(mesh, sbp, eqn, opts) # Evaluate function with perturbed q

      # Populate elem_res_vec with perturbed residuals
      for l = 1:mesh.numNodesPerElement
        for k = 1:mesh.numDofPerNode
          elem_res_vec[(l-1)*mesh.numDofPerNode + k] = eqn.res[k,l,i]
        end
      end

      fill!(eqn.res_vec, 0.0)
      array3DTo1D(mesh, sbp, eqn, opts, eqn.res,  eqn.res_vec)
      # println("elem_res_vec - eqn.res_vec = \n", elem_res_vec - eqn.res_vec)
      # println("orig_res_vec - eqn.res_vec = \n", orig_res_vec - eqn.res_vec)

      nl_solvers.calcJacRow(sview(elem_jac, :, j), elemRes0, elem_res_vec, pert)
    end  # End for j = = 1:numDofPerElement

    # println("elem_res_vec = \n", elem_res_vec)
    eqn.q[mesh.numDofPerNode,mesh.numNodesPerElement,i] = entry_orig

    elem_eigen_vals = eigvals(elem_jac)
    checkEigenValues(elem_eigen_vals, i, elem_count)

    #= println("eigen values of element ", i)
    counter = 0
    for kappa = 1:length(elem_eigen_vals)
      if elem_eigen_vals[kappa] < 0 && abs(elem_eigen_vals[kappa]) > 1e-12
        counter += 1
      end
    end =#

    #println(elem_eigen_vals, '\n', '\n')
    #println("number of negative element eigen values = ", counter)

  end

  println("Number of elements with negative eigen values = ", elem_count)
  =#
  res_0_norm = physicsRhs(mesh, sbp, eqn, opts, res_0, (evalResidual,))
  jac = zeros(Float64, mesh.numDof, mesh.numDof)
  res_0 = copy(eqn.res_vec)
  nl_solvers.calcJacFD(mesh, sbp, eqn, opts, evalResidual, res_0, pert, jac)
  println("Jacobian successfully computed")

  eigen_values = eigvals(jac)
  counter = 0
  for i = 1:length(eigen_values)
    if eigen_values[i] < 0 && abs(eigen_values[i]) > 1e-12
      #println(eigen_values[i])
      counter += 1
    end
  end
  #println(eigen_values)
  println("number of negative eigen values for assembled residual = ", counter)
  println("Number of eigen values of assembled residual = ", length(eigen_values))
  printSolution("eigen_values.dat", eigen_values)

  return nothing
end


text:
function test_gausslegendre()
    n = 4
    d = 2n
    qs = Q_GaussLegendre(n)
    @test integral(qs, x -> x^(d-1), 1..2) ≈ 2^d/d-1^d/d
    @test integral(BestRule(n), x -> x^(d-2), LegendreMeasure()) ≈ 2/(d-1)
end

function test_gaussjacobi()
    n = 4
    α = 1/2
    β = 3/4
    qs = Q_GaussJacobi(n, α, β)
    I = integral(x -> cos(x)*(1+x)^α*(1-x)^β, -1..1)
    @test abs(integral(qs, cos, -1..1) - I) < 1e-5
    @test abs(integral(BestRule(n), cos, JacobiMeasure(α, β)) -I) < 1e-5
end

function test_gausslaguerre()
    n = 15
    α = 1.2
    qs = Q_GaussLaguerre(n, α)
    I = integral(x->cos(x)*x^(α)*exp(-x), 0..10)
    @test abs(integral(qs, cos, HalfLine()) - I) < 1e-4
    @test abs(integral(BestRule(n), cos, LaguerreMeasure(α)) - I) < 1e-4
end

function test_gausshermite()
    n = 5
    qs = Q_GaussHermite(n)
    I = integral(x->cos(x)*exp(-x^2), -5..5)
    @test abs(integral(qs, cos, DomainSets.FullSpace{Float64}()) - I) < 1e-5
    @test abs(integral(BestRule(n), cos, HermiteMeasure()) - I) < 1e-5
end

function test_gauss()
    test_gausslegendre()
    test_gaussjacobi()
    test_gausslaguerre()
    test_gausshermite()
end


text:
module ModU

using CustomUnitRanges: filename_for_urange
include(filename_for_urange)

end

using ModU: URange

@testset "URange" begin
    @testset "Bool" begin
        r = URange(false, true)
        @test !isempty(r)
        @test length(r) == 2
        @test r[1] == false
        @test r[2] == true
        r = URange(true, false)
        @test isempty(r)
        @test length(r) == 0
        r = URange(false, false)
        @test !isempty(r)
        @test length(r) == 1
        @test r[1] == false
        r = URange(true, true)
        @test !isempty(r)
        @test length(r) == 1
        @test r[1] == true
    end

    @testset "Int/Float" begin
        for T in (Int, Float64)
            r = URange(-T(3),-T(5))
            @test isempty(r)
            @test length(r) == 0
            @test size(r) == (0,)
            r = URange(-T(5),T(3))
            rref = -5:3
            @test !isempty(r)
            @test length(r) == length(rref)
            @test size(r) == size(rref)
            @test step(r) == 1
            @test first(r) == -5
            @test last(r) == 3
            @test minimum(r) == -5
            @test maximum(r) == 3
            @test r[1] == -5
            @test r[2] == -4
            @test r[8] == 2
            @test r[9] == 3
            @test_throws BoundsError r[10]
            @test_throws BoundsError r[0]
            @test r[2:8] === URange(T(-4), T(2))
            @test r+1 === -T(4):T(4)
            @test 2*r === -T(10):T(2):T(6)
            k = -6
            for i in r
                @test i == (k+=1)
            end
            if T<:Integer
                @test intersect(r, URange(-5,2)) === URange(-5,2)
                @test intersect(r, -1:5) === -1:3
                @test intersect(r, 2:5) === 2:3
                @test intersect(r, 2) === intersect(2, r) === 2:2
                @test intersect(r, 7) === intersect(7, r) === 7:6
            else
                @test intersect(r, URange(-5,2)) == collect(URange(-5,2))
                @test intersect(r, -1:5) == collect(-1:3)
                @test intersect(r, 2:5) == collect(2:3)
                @test intersect(r, 2) == intersect(2, r) == [2.0]
                @test intersect(r, 7) == intersect(7, r) == Float64[]
            end
            @test findin(r, 0:5) == (T<:Integer ? (6:9) : collect(6:9))
            @test findin(r, URange(0,5)) == (T<:Integer ? (6:9) : collect(6:9))
            x, y = promote(r, 1:3)
            @test convert(URange{Int16}, 0:3) === URange{Int16}(0,3)
            @test convert(URange{Int16}, Base.OneTo(5)) === URange{Int16}(1,5)
            @test convert(URange{Int16}, r) === URange{Int16}(-5, 3)
            @test x === r
            @test y === URange(T(1),T(3))
            io = IOBuffer()
            show(io, r)
            str = takebuf_string(io)
            @test str == "ModU.URange(-$(T(5)),$(T(3)))"
        end
    end

    @testset "smallint" begin
        r = URange{Int16}(-5, 3)
        @test eltype(r) == Int16
        k = -6
        for i in r
            @test i == (k+=1)
        end
    end

    @testset "promotion" begin
        r = URange{Int16}(-5, 3)
        ri = URange(2,20)
        x, y = promote(r, ri)
        @test isa(x, URange{Int})
        @test isa(y, URange{Int})
        x, y = promote(r, 2:20)
        @test isa(x, URange{Int})
        @test isa(y, URange{Int})
        x, y = promote(r, Base.OneTo(3))
        @test isa(x, URange{Int})
        @test isa(y, URange{Int})
    end
end


text:
"""
Brute-force (𝒪(n²)) kNN search algorithm.
Returns an KxN array of tuples (index, distance) of the k nearest neighbors
for each point in `data`.
"""
function brute_knn(data::Vector{V},
                   metric::M,
                   k::Int) where {V <: AbstractArray, M <: SemiMetric}

    np = length(data)
    Dtype = result_type(metric, data[1], data[1])
    distances = Matrix{NNTuple{Int, Dtype}}(undef, np, np)

    @inbounds @fastmath for j = 1:np, i = 1:np
        d = evaluate(metric, data[i], data[j])
        distances[i, j] = NNTuple(j, d)
    end

    # start at 2 as index 1 is the distance to itself (once sorted)
    knn_graph_tuples = sort(distances, dims=2)[1:end, 2:k+1]
    knn_graph = Matrix{Tuple{Int, Dtype}}(undef, k, np)
    for i = 1:k, j = 1:np
        knn_graph[i, j] = (knn_graph_tuples[j, i].idx,
                           knn_graph_tuples[j, i].dist)
    end
    return knn_graph
end

"""
    brute_search(data, queries, k, metric) -> idxs, dists

Search for the nearest `k` neighbors of `q` in `queries` in `data`.
"""
function brute_search(data::Vector{V},
                      queries::Vector{V},
                      k::Integer,
                      metric::M=Euclidean()) where {V <: AbstractVector,
                                                    M <: SemiMetric}
    np = length(data)
    nq = length(queries)
    Dtype = result_type(metric, data[1], data[1])
    distances = Matrix{NNTuple{Int, Dtype}}(undef, np, nq)

    @inbounds @fastmath for i = 1:nq, j = 1:np
        d = evaluate(metric, queries[i], data[j])
        distances[j, i] = NNTuple(j, d)
    end

    knn_graph_tuples = sort(distances, dims=1)[1:k, 1:end]
    knn_graph = Matrix{Tuple{Int, Dtype}}(undef, k, nq)
    for i = 1:k, j = 1:nq
        knn_graph[i, j] = (knn_graph_tuples[i, j].idx, knn_graph_tuples[i, j].dist)
    end
    return getindex.(knn_graph, 1), getindex.(knn_graph, 2)
end


text:
module KernelBenchmarks

export  KernelParam,
        # kernel types
        ReadOnly,
        WriteOnly,
        ReadWrite,
        # load/store types
        Standard,
        Nontemporal,
        # itertors
        Sequential,
        PseudoRandom,
        # exported function
        execute!,
        # interactive
        introspect

##### stdlib
using InteractiveUtils
import REPL
using REPL.TerminalMenus

##### "Internal" Packages
# These are packages I've developed for testing purposes and not part of the
# greater Julia ecosystem.
using MaxLFSR

##### External Packages
using DataFrames

# For vector load intrinsics
using SIMD

include("threaded.jl")
include("kernels.jl")
include("introspection.jl")
include("deprecated.jl")

function run(A)
    df = DataFrame()
    num_threads = [1, 4, 8, 10, 16, 20]
    kernels = [ReadOnly, WriteOnly]
    unroll = 8
    loops = [1, 2, 4, 8]

    iter = Iterators.product(num_threads, loops, kernels)
    for tup in iter
        nthreads = tup[1]
        nloops = tup[2]
        kernel = tup[3]
        num_iterations = kernel == ReadOnly ? 3 : 1

        params = KernelParam(;
            kernel = kernel,
            loadtype = Standard,
            storetype = Nontemporal,
            iterator = PseudoRandom,
            vectorsize = 16,
            unroll = unroll,
            loops = nloops,
        )

        # Warmup iteration.
        threaded(execute!, A, params)

        # Actual run.
        rt = @elapsed threaded(
            execute!,
            A,
            params;
            iterations = num_iterations,
            nthreads = nthreads,
        )

        bw = (num_iterations * sizeof(A)) / (rt)

        stats = Dict(
            :num_threads => nthreads,
            :num_bytes => 64 * unroll * nloops,
            :kernel => string(kernel),
            :bandwidth => bw,
        )
        display(stats),

        push!(df, stats; cols = :union)
    end
    return df
end

end # module


text:
type CRED
    user::String
    pass::String
    personal::String
    secret::String
    token::String
end

CRED(user, pass, personal, secret) = CRED(user, pass, personal, secret, "")
CRED(user, personal, secret) = CRED(user, chomp(readline()), personal, secret, "")

type Subreddit
	name::String
	sorting::String
	ids::Vector{String}
    responses::Vector{Response}
    count::Int64
end


text:
using TropicalTensors
using TropicalTensors: getbonds, getvertices, rand_Js, rand_hs
using Test
using Random

@testset "spinglass" begin
    lt = CubicLattice(4,4,4)
    sg = Spinglass(lt, rand_Js(Float32, Ferro(), lt), rand_hs(Float32, Zero(), lt))
    res = solve(sg; usecuda=false)
    @test res.n == 144
end

@testset "counting tropical" begin
    lt = CubicLattice(4,4,4)
    sg = Spinglass(lt, rand_Js(Float32, Ferro(), lt), rand_hs(Float32, Zero(), lt))
    res = solve(CountingTropicalF32, sg; usecuda=false)
    @test res.n == 144
    @test res.c == 2
end

using BitBasis
function count_degeneracy_exact(sg::Spinglass{LT,T}) where {LT, T}
    elist = T[]
    bonds = getbonds(sg.lattice)
    lt = sg.lattice
    LI = LinearIndices((lt.Nx, lt.Ny, lt.Nz))
    for b in basis(BitStr64{length(sg.lattice)})
        eng = zero(T)
        for ((i, j), J) in sg.Js
            eng += (b[LI[i...]] == b[LI[j...]]) ? J : -J
        end
        push!(elist, eng)
    end
    findall(==(minimum(elist)), elist) |> length
end

@testset "counting tropical 2" begin
    Random.seed!(3)
    lt = CubicLattice(2, 3, 3)
    sg = Spinglass(lt, rand_Js(Float32, Randpm(), lt), rand_hs(Float32, Zero(), lt))
    res = solve(CountingTropical{Int,Int}, sg; usecuda=false)
    @test res.c == count_degeneracy_exact(sg)
end


text:
using Test

m = test_page_model()
include("../src/climate_model.jl")

rt_g = m[:ClimateTemperature, :rt_g_globaltemperature]
rt_g_compare = readpagedata(m, "test/validationdata/rt_g_globaltemperature.csv")

@test rt_g ≈ rt_g_compare rtol=1e-4


text:
using Unitful, Dates
using PrettyTables

export CouplerState
export coupler_push!, coupler_pull!, coupler_put!, coupler_get
export coupler_add_field!

# TODO: Build constructor that uses a model component's grid.
struct CplGridInfo{GT, GP, GH, GE}
    gridtype::GT
    gridparams::GP
    gridhandle::GH
    gridencoding::GE
end

mutable struct CplFieldInfo{AT, GT <: CplGridInfo, UT <: Unitful.Units}
    data::AT
    gridinfo::GT
    units::UT
    datetime::DateTime
end

mutable struct CouplerState{DT}
    CplStateBlob::DT
end

"""
    CouplerState()

Type for holding coupler "state". This is the namespace through which coupled components
communicate. Its role is to provide a level of indirection so that components remain modular
and so that any data communication, interpolation, reindexing/unit conversions and filtering
etc... can be embeded in the intermdediate coupling layer.

A field is exported by one component and imported by one or more other components.
"""
function CouplerState()
    return CouplerState(Dict{Symbol, CplFieldInfo}())
end

"""
    coupler_add_field!(
            coupler::CouplerState,
            fieldname::Symbol,
            fieldvalue,
            grid,
            datetime::DateTime,
            units::Unitful.Units = Unitful.NoUnits,
        )

Add a field to the coupler that is accessible with key `fieldname`.

# Arguments
- `coupler`: coupler object the field is added to.
- `fieldname`: key to access the field in the coupler.
- `fieldvalue`: data array of field values.
- `grid`: grid the field is stored on.
- `datetime`: time associated with the field state.
- `units`: units associated with the field values. Dimensionless by default.
"""
function coupler_add_field!(
    coupler::CouplerState,
    fieldname::Symbol,
    fieldvalue,
    grid,
    datetime::DateTime,
    units::Unitful.Units = Unitful.NoUnits,
)
    gridinfo = CplGridInfo(nothing, nothing, nothing, nothing)
    push!(coupler.CplStateBlob, fieldname => CplFieldInfo(fieldvalue, gridinfo, units, datetime))
end

"""
    coupler_pull!(model, coupler::CouplerState)

Update model with fields retrieved from the coupler.

`coupler_pull!` is an adapter function to be implemented for each
model component using the coupler. It should get coupling fields via
`coupler_get` calls and perform any operations on these fields to prepare
them for use in the component model.
"""
function coupler_pull!(model, coupler::CouplerState) end

"""
    coupler_get(coupler::CouplerState, fieldname::Symbol, gridinfo, datetime::DateTime, units::Unitful.Units)

Retrieve data array corresponding to `fieldname`.

Returns data on the grid specified by `gridinfo` and in the units of `units`. Checks that
the coupler data field is the state at time `datetime`.
"""
function coupler_get(coupler::CouplerState, fieldname::Symbol, gridinfo, datetime::DateTime, units::Unitful.Units)
    cplfield = coupler.CplStateBlob[fieldname]

    # check that retrieving component and coupler are at same time
    datetime != cplfield.datetime &&
    throw(ErrorException("Retrieval time ($datetime) != coupler field time ($(cplfield.datetime))"))

    regriddata = regrid(cplfield.data, gridinfo, cplfield.gridinfo)

    return uconvert(regriddata, units, cplfield.units)
end

"""
    coupler_push!(coupler::CouplerState, model)

Update coupler with fields retrieved from the coupler.

`coupler_push!` is an adapter function to be implemented for each
model component using the coupler. It should send coupling fields via
`coupler_put!` calls and perform any operations on these fields to prepare
them for the coupler.
"""
function coupler_push!(coupler::CouplerState, model) end

"""
    coupler_put!(coupler::CouplerState, fieldname::Symbol, fieldvalue, gridinfo, datetime::DateTime, units::Unitful.Units)

Updates coupler field `fieldname` with `fieldvalue`, the field's value at time `datetime`.

`gridinfo` and `units` inform the coupler of the format of the inputted data
allowing conversion to match the grid and units of the coupler field.
"""
function coupler_put!(
    coupler::CouplerState,
    fieldname::Symbol,
    fieldvalue,
    gridinfo,
    datetime::DateTime,
    units::Unitful.Units,
)
    cplfield = coupler.CplStateBlob[fieldname]

    # map new data to grid of coupler field; new data -> coupler grid
    regriddata = regrid(fieldvalue, cplfield.gridinfo, gridinfo)

    # store new regridded data in correct units
    cplfield.data .= uconvert(regriddata, cplfield.units, units)

    # update timestamp of coupler field
    settime!(cplfield, datetime)

    return nothing
end

# TODO: maps cplfield data from `fromgrid` to `togrid`
function regrid(data, togrid, fromgrid)
    return data
end

# convert from `fromunits` to `tounits`
# ustrip throws error if the units aren't dimensionally compatible
function uconvert(data, tounits::Unitful.Units, fromunits::Unitful.Units)
    tounits == fromunits && return data
    return ustrip.(tounits, data .* fromunits)
end

# set the time of the coupled field
function settime!(cplfield::CplFieldInfo, newtime::DateTime)
    cplfield.datetime = newtime
end

# display table of registered fields & their info
function Base.show(io::IO, coupler::CouplerState)
    #=
    ---------------------------------------------------------------------------
    |  Field Name  |  Units  |  Field Time (Y-m-d H:M:S)   |     Grid Type    |
    ---------------------------------------------------------------------------
    |  :OceanSST   |   °C    |  2021-01-01T00:00:00        |  DG Cubed Sphere |
    |  :Field2     |    m    |  2021-01-01T00:00:00        |  CG Cubed Sphere |
    ---------------------------------------------------------------------------
    =#
    fields = coupler.CplStateBlob
    data = Array{Any}(undef, length(fields), 4)
    for (i, k) in enumerate(keys(coupler.CplStateBlob))
        entry = coupler.CplStateBlob[k]
        data[i, :] = [k entry.units entry.datetime entry.gridinfo.gridtype]
    end
    header = (["Field Name", "Units", "Field Time", "Grid Type"], ["", "", "Y-m-d H:M:S", ""])
    pretty_table(data, header = header)
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine

max_forks_repo_name: CliMA/CouplerMachine


text:
"""
Assign ranks to the values but keeps the rank the same if the values are within
tolerance of each other.
"""
function assign_ranks_within_tolerance(values; by = (x) -> x, tolerance = 1e-5, rev = false)
    perm = sortperm(values, by = by, rev = rev)
    ranked = Any[]
    rank = 1
    prev = by(values[perm[1]])
    num_of_this_rank = 0

    for i in eachindex(perm)
        r = values[perm[i]]
        v = by(r)
        if abs(prev - v) > tolerance
            rank += num_of_this_rank
            num_of_this_rank = 1
        else
            num_of_this_rank += 1
        end
        push!(ranked, (rank, r, v))
        prev = v
    end

    return ranked
end


text:
# linear systems
E = [0. 1; 1 0]
A = [1. 1; 1 -1]
B = Matrix([0.5 1.5]')
C = [1.; 1.]
D = [1. 2; 0 1]
X = Line([1., -1], 0.) # line x = y
U = Interval(0.9, 1.1)
W = BallInf(zeros(2), 2.0)
sd = 2
id = 1
nd = 2

# polynomial system
@polyvar x y
p = 2x^2 - 3x + y
sdp = 2

# blackbox system
add_one(x) = x .+ 1
add_one(x, u) = x .+ 1 .+ u
add_one(x, u, w) = x .+ 1 .+ u .+ w
state = [1.0; 2.0]
input = 1
noise = [3.0; 1.0]
statePlusOne = add_one(state)
stateInputPlusOne = add_one(state, input)
stateInputNoisePlusOne = add_one(state, input, noise)

# Scalar System
a = 1.; b = 2.; c = 0.1; d = 3.; Xs = 1; Us = 2; Ws = 3; e = 2.;
As = [a][:,:]; Bs = [b][:,:]; Cs = [c]; Ds = [d][:,:]; Es = [e][:,:]

@testset "Continuous identity system" begin
    s = ContinuousIdentitySystem(sd)
    @test state_matrix(s) == I(sd)
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous constrained identity system" begin
    s = ConstrainedContinuousIdentitySystem(sd, X)
    @test state_matrix(s) == I(sd)
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
end

@testset "Continuous linear system" begin
    s = LinearContinuousSystem(A)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = LinearContinuousSystem(a)
    @test scalar_sys == LinearContinuousSystem(As)
end

@testset "Continuous affine system" begin
    s = AffineContinuousSystem(A, C)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == C
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = AffineContinuousSystem(a, c)
    @test scalar_sys == AffineContinuousSystem(As, Cs)
end

@testset "Continuous linear control system" begin
    s = LinearControlContinuousSystem(A, B)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = LinearControlContinuousSystem(a, b)
    @test scalar_sys == LinearControlContinuousSystem(As, Bs)
end

@testset "Continuous constrained linear system" begin
    s = ConstrainedLinearContinuousSystem(A, X)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = ConstrainedLinearContinuousSystem(a, Xs)
    @test scalar_sys == ConstrainedLinearContinuousSystem(As, Xs)
end

@testset "Continuous constrained affine system" begin
    s = ConstrainedAffineContinuousSystem(A, C, X)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == C
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = ConstrainedAffineContinuousSystem(a, c, Xs)
    @test scalar_sys == ConstrainedAffineContinuousSystem(As, Cs, Xs)
end

@testset "Continuous affine control system with state constraints" begin
    s = ConstrainedAffineControlContinuousSystem(A, B, C, X, U)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == C
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = ConstrainedAffineControlContinuousSystem(a, b, c, Xs, Us)
    @test scalar_sys == ConstrainedAffineControlContinuousSystem(As, Bs, Cs, Xs, Us)
end

@testset "Continuous affine control system with state constraints" begin
    A = zeros(2, 2)
    B = Matrix([0.5 1.5]')
    c = [0.0, 1.0]
    s = AffineControlContinuousSystem(A, B, c)
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test noisedim(s) == 0
    for s = [s, typeof(s)]
        @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    a = 1.; b = 2.; c = 0.1
    A = [a][:,:]; B = [b][:,:]; C = [c]
    scalar_sys = AffineControlContinuousSystem(a, b, c)
    @test scalar_sys == AffineControlContinuousSystem(A, B, C)
end

@testset "Continuous constrained linear control system" begin
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end

    @testset "initial value problem composite type" begin
        x0 = Singleton([1.5, 2.0])
        ivp = InitialValueProblem(s, x0)
        @test initial_state(ivp) == x0
        @test system(ivp) == s
        @test inputset(ivp) == U
        @test statedim(ivp) == 2
        @test inputdim(ivp) == 1
    end

    # check that the matrices A and B need not be of the same type
    # with A dense and B a lazy adjoint
    B_adjoint = [0.5 1.5]'
    s = ConstrainedLinearControlContinuousSystem(A, B_adjoint, X, U)

    # check that the matrices A and B need not be of the same type
    # with A sparse and B a lazy adjoint
    A_sparse = sparse([1], [2], [1.0], 2, 2) # sparse matrix
    s = ConstrainedLinearControlContinuousSystem(A_sparse, B, X, U)
end

@testset "Continuous linear algebraic system" begin
    s = LinearAlgebraicContinuousSystem(A, E)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test s.E == E
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    # Scalar System
    scalar_sys = LinearAlgebraicContinuousSystem(a, e)
    @test scalar_sys == LinearAlgebraicContinuousSystem(As, Es)
end

@testset "Continuous constrained linear algebraic system" begin
    s = ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test s.E == E
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = ConstrainedLinearAlgebraicContinuousSystem(a, e, Xs)
    @test scalar_sys == ConstrainedLinearAlgebraicContinuousSystem(As, Es, Xs)

    @testset "Initial value problem" begin
        x0 = Singleton([1.5, 2.0])
        ivp = IVP(s, x0)
        # getter functions
        @test initial_state(ivp) == x0
    end
end

@testset "Polynomial system in continuous time" begin
    # default constructor for scalar p and
    s = PolynomialContinuousSystem(p)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    # @test s.p == p
    @test statedim(s) == sdp
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end

    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)

    # constructor for scalar p and given number of variables
    s = PolynomialContinuousSystem([p], 2)
end

@testset "Polynomial system in continuous time with state constraints" begin
    # default constructor for scalar p and
    s = ConstrainedPolynomialContinuousSystem(p, X)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    # @test s.p == p
    @test statedim(s) == sdp
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end

    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)
end

function vanderpol!(x, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1]
    return dx
end

@testset "Continuous system defined by a function" begin
    s = BlackBoxContinuousSystem(vanderpol!, 2)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]

    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    s = ConstrainedBlackBoxContinuousSystem(vanderpol!, 2, H)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == H
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]
    @test stateset(s) == H
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
end

function vanderpol_controlled!(x, u, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1] + u[1]
    return dx
end

@testset "Continuous control black-box system" begin
    add_one(x) = x + 1
    s = BlackBoxControlContinuousSystem(vanderpol_controlled!, 2, 1)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test noisedim(s) == 0
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous control system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    U = Interval(-0.1, 0.1)
    s = ConstrainedBlackBoxControlContinuousSystem(vanderpol_controlled!, 2, 1, H, U)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test noisedim(s) == 0
    @test stateset(s) == H
    @test inputset(s) == U
    @test noiseset(s) == nothing
    f! = (x, u, dx) -> s.f(x, u, dx)
    x = [1.0, 0.0]
    u = an_element(U)
    dx = similar(x)
    f!(x, u, dx)
    @test dx ≈ [0.0, -1.0 + u[1]]
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end
end

# ==============
# Noisy systems
# ==============

@testset "Noisy continuous linear system" begin
    s = NoisyLinearContinuousSystem(A, D)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == nd
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyLinearContinuousSystem(a, d)
    @test scalar_sys == NoisyLinearContinuousSystem(As, Ds)
end

@testset "Noisy Continuous constrained linear system" begin
    s = NoisyConstrainedLinearContinuousSystem(A, D, X, W)
    @test state_matrix(s) == A
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == 0
    @test noisedim(s) == nd
    @test stateset(s) == X
    @test inputset(s) == nothing
    @test noiseset(s) == W
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyConstrainedLinearContinuousSystem(a, d, Xs, Ws)
    @test scalar_sys == NoisyConstrainedLinearContinuousSystem(As, Ds, Xs, Ws)
end

@testset "Noisy continuous control linear system" begin
    s = NoisyLinearControlContinuousSystem(A, B, D)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == nothing
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyLinearControlContinuousSystem(a, b, d)
    @test scalar_sys == NoisyLinearControlContinuousSystem(As, Bs, Ds)
end

@testset "Noisy Continuous constrained control linear system" begin
    s = NoisyConstrainedLinearControlContinuousSystem(A, B, D, X, U, W)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == nothing
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    for s = [s, typeof(s)]
        @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyConstrainedLinearControlContinuousSystem(a, b, d, Xs, Us, Ws)
    @test scalar_sys == NoisyConstrainedLinearControlContinuousSystem(As, Bs, Ds, Xs, Us, Ws)
end

@testset "Noisy continuous control affine system" begin
    s = NoisyAffineControlContinuousSystem(A, B, C, D)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == C
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyAffineControlContinuousSystem(a, b, c, d)
    @test scalar_sys == NoisyAffineControlContinuousSystem(As, Bs, Cs, Ds)
end

@testset "Noisy Continuous constrained control affine system" begin
    s = NoisyConstrainedAffineControlContinuousSystem(A, B, C, D, X, U, W)
    @test state_matrix(s) == A
    @test input_matrix(s) == B
    @test affine_term(s) == C
    @test noise_matrix(s) == D
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    for s = [s, typeof(s)]
        @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end
    # Scalar System
    scalar_sys = NoisyConstrainedAffineControlContinuousSystem(a, b, c, d, Xs, Us, Ws)
    @test scalar_sys == NoisyConstrainedAffineControlContinuousSystem(As, Bs, Cs, Ds, Xs, Us, Ws)
end

@testset "Noisy continuous control black-box system" begin
    s = NoisyBlackBoxControlContinuousSystem(add_one, sd, id, nd)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test s.f(state, input, noise) ≈ stateInputNoisePlusOne
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == nothing
    @test inputset(s) == nothing
    @test noiseset(s) == nothing
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Noisy Continuous constrained control blackbox system" begin
    s = NoisyConstrainedBlackBoxControlContinuousSystem(add_one, sd, id, nd, X, U, W)
    @test state_matrix(s) == nothing
    @test input_matrix(s) == nothing
    @test affine_term(s) == nothing
    @test noise_matrix(s) == nothing
    @test s.f(state, input, noise) ≈ stateInputNoisePlusOne
    @test statedim(s) == sd
    @test inputdim(s) == id
    @test noisedim(s) == nd
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    for s = [s, typeof(s)]
        @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
    end
end

@testset "Second order linear systems" begin
    M = [1. 0; 0 2]
    C = [0.1 0; 0 0.2]
    K = [2. 1; 0 1]
    b = [1., 0]
    B = hcat([1., 0])
    d = [1., 0]
    X = BallInf(zeros(2), 1.0)
    U = Singleton(ones(2))

    sys = SecondOrderLinearContinuousSystem(M, C, K)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C && stiffness_matrix(sys) == K

    sys = SecondOrderAffineContinuousSystem(M, C, K, b)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C && stiffness_matrix(sys) == K
    @test affine_term(sys) == b

    sys = SecondOrderConstrainedLinearControlContinuousSystem(M, C, K, B, X, U)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C && stiffness_matrix(sys) == K
    @test input_matrix(sys) == B && stateset(sys) == X && inputset(sys) == U

    sys = SecondOrderConstrainedAffineControlContinuousSystem(M, C, K, B, d, X, U)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C && stiffness_matrix(sys) == K
    @test affine_term(sys) == d && input_matrix(sys) == B && stateset(sys) == X && inputset(sys) == U
end

@testset "Second order systems" begin
    M = [1. 0; 0 2]
    C = [0.1 0; 0 0.2]
    fi(x) = x + x.^2 + ones(2)
    fe = zeros(2)
    sys = SecondOrderContinuousSystem(M, C, fi, fe)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C

    fe = zeros(2)
    X = Universe(2)
    U = Universe(2)
    sys = SecondOrderConstrainedContinuousSystem(M, C, fi, fe, X, U)
    @test mass_matrix(sys) == M && viscosity_matrix(sys) == C
    @test stateset(sys) === X && inputset(sys) === U
end


text:
export femidx, fvmidx, get_edge_normal,
plot_u,bcnode,bcedge, layer_model, gauss_nodes, fem_nodes, fvm_nodes, subdomain, interior_node,
cholesky_outproduct, cholesky_factorize, cholesky_logdet, fem_randidx, fem_to_gauss_points, fem_to_fvm

"""
    femidx(d::Int64, m::Int64)

Returns the FEM index of the dof `d`. Basically, `femidx` is the inverse of

```
(i,j) → d = (j-1)*(m+1) + i
```
"""
function femidx(i::Int64, m::Int64)
    ii = mod(i, m+1)
    ii = ii == 0 ? m+1 : ii
    jj = div(i-1, m+1) + 1
    return ii, jj
end

"""
    fvmidx(d::Int64, m::Int64)

Returns the FVM index of the dof `d`. Basically, `femidx` is the inverse of

```
(i,j) → d = (j-1)*m + i
```
"""
function fvmidx(i::Int64, m::Int64)
    ii = mod(i, m)
    ii = ii == 0 ? m : ii
    jj = div(i-1, m) + 1
    return ii, jj
end

function plot_u(u::Array{Float64}, m::Int64, n::Int64, h::Float64)
    x = LinRange(0,m*h,m+1)
    y = LinRange(0,n*h,n+1)
    X, Y = np.meshgrid(x, y)
    U = reshape(u[1:(m+1)*(n+1)], m+1, n+1)'|>Array
    V = reshape(u[(m+1)*(n+1)+1:end], m+1, n+1)'|>Array
    X, Y, U, V
end

@doc raw"""
    layer_model(u::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)

Convert the vertical profile of a quantity to a layer model.
The input `u` is a length $n$ vector, the output is a length $4mn$ vector, representing the $4mn$ Gauss points.
"""
function layer_model(u::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)
    repeat(reshape(u, :, 1), 1, 4m)'[:]
end

@doc raw"""
    layer_model(u::PyObject, m::Int64, n::Int64, h::Float64)

A differential kernel for [`layer_model`](@ref).
"""
function layer_model(u::PyObject, m::Int64, n::Int64, h::Float64)
    reshape(repeat(u, 1, 4m), (-1,))
end

"""
    get_edge_normal(edge::Array{Int64,1}, m::Int64, n::Int64, h::Float64)

Returns the normal vector given edge `edge`.
"""
function get_edge_normal(edge::Array{Int64,1}, m::Int64, n::Int64, h::Float64)
    i1, j1 = femidx(edge[1], m)
    i2, j2 = femidx(edge[2], m)
    if i1==i2==1
        return [-1.0;0.0]
    elseif i1==i2==m+1
        return [1.0;0.0]
    elseif j1==j2==1
        return [0.0;-1.0]
    elseif j1==j2==n+1
        return [0.0;1.0]
    else
        error("Invalid edge index $edge")
    end
end

function get_edge_normal(edges::Array{Int64,2}, m::Int64, n::Int64, h::Float64)
    out = zeros(size(edges,1), 2)
    for i = 1:size(edges,1)
        out[i,:] = get_edge_normal(edges[i,:], m, n, h)
    end
end

"""
    bcedge(desc::String, m::Int64, n::Int64, h::Float64)

Returns the edge indices for description. See [`bcnode`](@ref)
"""
function bcedge(desc::String, m::Int64, n::Int64, h::Float64)
    desc=="all" && (desc="left|right|upper|lower")
    descs = strip.(split(desc, '|'))
    edges = bcedge_.(descs, m, n, h)
    unique(vcat(edges...), dims=1)
end

function bcedge_(desc::AbstractString, m::Int64, n::Int64, h::Float64)
    nodes = bcnode_(desc, m, n, h)
    [nodes[1:end-1] nodes[2:end]]
end

"""
    bcnode(desc::String, m::Int64, n::Int64, h::Float64)

Returns the node indices for the description. Multiple descriptions can be concatented via `|`

```
                upper
        |------------------|
left    |                  | right
        |                  |
        |__________________|

                lower
```

# Example
```julia
bcnode("left|upper", m, n, h)
```
"""
function bcnode(desc::String, m::Int64, n::Int64, h::Float64)
    desc=="all" && (desc="left|right|upper|lower")
    descs = strip.(split(desc, '|'))
    nodes = bcnode_.(descs, m, n, h)
    unique(vcat(nodes...))
end

function bcnode_(desc::AbstractString, m::Int64, n::Int64, h::Float64)
    nodes = Int64[]
    if desc=="upper" || desc=="top"
        for i = 1:m+1
            push!(nodes, i)
        end
    elseif desc=="lower" || desc == "bottom"
        for i = 1:m+1
            push!(nodes, i+n*(m+1))
        end
    elseif desc=="left"
        for j = 1:n+1
            push!(nodes, 1+(j-1)*(m+1))
        end
    elseif desc=="right"
        for j = 1:n+1
            push!(nodes, m+1+(j-1)*(m+1))
        end
    else
        error("$desc is not a valid specification. Only `upper` (top), `lower` (bottom), `left`, `right`, `all` or their combination via `|`, is accepted.")
    end
    return nodes
end

"""
    interior_node(desc::String, m::Int64, n::Int64, h::Float64)

In contrast to [`bcnode`](@ref), `interior_node` returns the nodes that are not specified by `desc`, including those on the boundary.
"""
function interior_node(desc::String, m::Int64, n::Int64, h::Float64)
    bc = bcnode(desc, m, n, h)
    collect(setdiff(Set(1:(m+1)*(n+1)), bc))
end

"""
    fem_to_fvm(u::Union{PyObject, Array{Float64}}, m::Int64, n::Int64, h::Float64)

Interpolates the nodal values of `u` to cell values.
"""
function fem_to_fvm(u::Union{PyObject, Array{Float64}}, m::Int64, n::Int64, h::Float64)
    idx1 = zeros(Int64, m*n)
    idx2 = zeros(Int64, m*n)
    idx3 = zeros(Int64, m*n)
    idx4 = zeros(Int64, m*n)
    for i = 1:m
        for j = 1:n
            idx1[(j-1)*m+i] = (j-1)*(m+1)+i
            idx2[(j-1)*m+i] = (j-1)*(m+1)+i+1
            idx3[(j-1)*m+i] = j*(m+1)+i
            idx4[(j-1)*m+i] = j*(m+1)+i+1
        end
    end
    (u[idx1] + u[idx2] + u[idx3] + u[idx4])/4
end

@doc raw"""
    gauss_nodes(m::Int64, n::Int64, h::Float64)

Returns the node matrix of Gauss points for all elements. The matrix has a size $4mn\times 2$
"""
function gauss_nodes(m::Int64, n::Int64, h::Float64)
    xs = zeros(4*m*n, 2)
    k = 1
    for i = 1:m
        for j = 1:n
            idx = (j-1)*m + i
            x1 = (i-1)*h
            y1 = (j-1)*h
            for p = 1:2
                for q = 1:2
                    k = (idx-1)*4 + 2*(q-1) + p
                    ξ = pts[p]; η = pts[q]
                    x = x1 + ξ*h; y = y1 + η*h
                    xs[k,:] = [x;y]
                    k += 1
                end
            end
        end
    end
    xs
end

@doc raw"""
    fem_nodes(m::Int64, n::Int64, h::Float64)

Returns the FEM node matrix of size $(m+1)(n+1)\times 2$
"""
function fem_nodes(m::Int64, n::Int64, h::Float64)
    xs = zeros((m+1)*(n+1), 2)
    k = 1
    for i = 1:m+1
        for j = 1:n+1
            idx = (j-1)*(m+1) + i
            x1 = (i-1)*h
            y1 = (j-1)*h
            xs[idx,:] = [x1;y1]
            k += 1
        end
    end
    xs
end

@doc raw"""
    fvm_nodes(m::Int64, n::Int64, h::Float64)

Returns the FVM node matrix of size $(m+1)(n+1)\times 2$
"""
function fvm_nodes(m::Int64, n::Int64, h::Float64)
    xs = zeros(m*n, 2)
    k = 1
    for i = 1:m
        for j = 1:n
            idx = (j-1)*m + i
            x1 = (i-1/2)*h
            y1 = (j-1/2)*h
            xs[idx,:] = [x1;y1]
            k += 1
        end
    end
    xs
end

@doc raw"""
    fem_randidx(N::Int64, m::Int64, n::Int64, h::Float64)

Returns $N$ random index
"""
function fem_randidx(N::Int64, m::Int64, n::Int64, h::Float64)
    idx = Set([])
    for k = 1:N
        while true
            i = rand(2:m)
            j = rand(2:n)
            if j*(m+1)+i in idx
                continue
            else
                push!(idx, j*(m+1)+i)
                break
            end
        end
    end
    sort([k for k in idx])
end

"""
    subdomain(f::Function, m::Int64, n::Int64, h::Float64)

Returns the subdomain defined by `f(x, y)==true`.
"""
function subdomain(f::Function, m::Int64, n::Int64, h::Float64)
    nodes = fem_nodes(m, n, h)
    idx = f.(nodes[:,1], nodes[:,2])
    findall(idx)
end

@doc raw"""
    cholesky_outproduct(L::Union{Array{<:Real,2}, PyObject})

Returns
$$A = LL'$$
where `L` (length=6) is a vectorized form of $L$
$$L = \begin{matrix}
l_1 & 0 & 0\\
l_4 & l_2 & 0 \\
l_5 & l_6 & l_3
\end{matrix}$$
and `A` (length=9) is also a vectorized form of $A$
"""
function cholesky_outproduct(A::Union{Array{<:Real,2}, PyObject})
    @assert size(A,2)==6
    op_ = load_op_and_grad(libadfem,"cholesky_backward_op")
    A = convert_to_tensor([A], [Float64]); A = A[1]
    L = op_(A)
end

@doc raw"""
    cholesky_factorize(A::Union{Array{<:Real,2}, PyObject})

Returns the cholesky factor of `A`. See [`cholesky_outproduct`](@ref) for details.
"""
function cholesky_factorize(A::Union{Array{<:Real,2}, PyObject})
    @assert size(A,2)==9
    op_ = load_op_and_grad(libadfem,"cholesky_forward_op")
    A = convert_to_tensor([A], [Float64]); A = A[1]
    L = op_(A)
end

@doc raw"""
    cholesky_logdet(A::Union{Array{<:Real,2}, PyObject})

Returns the cholesky factor of `A` as well as the log determinant. See [`cholesky_outproduct`](@ref) for details.
"""
function cholesky_logdet(A::Union{Array{<:Real,2}, PyObject})
    op_ = load_op_and_grad("$(@__DIR__)/../depsCholeskyOp/build/libCholeskyOp/build","cholesky_logdet")
    A = convert_to_tensor(A, dtype=Float64)
    L, J = op_(A)
end

@doc raw"""
    fem_to_gauss_points(u::PyOject, m::Int64, n::Int64, h::Float64)
    fem_to_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64)

Given a vector of length $(m+1)(n+1)$, `u`, returns the function values at each Gauss point.

Returns a vector of length $4mn$.
"""
function fem_to_gauss_points(u::PyObject, m::Int64, n::Int64, h::Float64)
    fem_to_gauss_points_ = load_op_and_grad(libadfem,"fem_to_gauss_points")
    u,m_,n_,h = convert_to_tensor(Any[u,m,n,h], [Float64,Int64,Int64,Float64])
    out = fem_to_gauss_points_(u,m_,n_,h)
    out = set_shape(out, (4*m*n,))
end

function fem_to_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64)
    rhs = zeros(4*m*n)
    for j = 1:n
        for i = 1:m
            cell_idx = (j-1)*m + i
            uA = u[[(j-1)*(m+1)+i; (j-1)*(m+1)+i+1; j*(m+1)+i; j*(m+1)+i+1]]
            for q = 1:2
                for p = 1:2
                    ξ = pts[p]; η = pts[q]
                    k = (cell_idx-1)*4 + 2*(q-1) + p
                    rhs[k] = uA[1] * (1-ξ)*(1-η) +
                             uA[2] * ξ*(1-η) +
                             uA[3] * (1-ξ)*η +
                             uA[4] *  ξ  * η
                end
            end
        end
    end
    return rhs
end