text:
# Ideally we would register a single job definition as part of the CloudFormation template
# and use overrides to change the image used. Unfortunately, this is not a supported
# override so we're left with a dilemma:
#
# 1. Define the job definition in CFN and change the image referenced by uploading a new
#    Docker image with the same name.
# 2. Create a new job definition for each Docker image

# The first option is problematic when executing tests in parallel using the same stack. If
# two CI pipelines are running concurrently then the last Docker image to be pushed will be
# the one used by both pipelines (assumes the push completes before the batch job starts).
#
# We went with the second option as it is safer for parallel pipelines and allows us to
# still use overrides to modify other parts of the job definition.
function batch_node_job_definition(;
    job_definition_name::AbstractString="$(STACK_NAME)-node",
    image::AbstractString=TEST_IMAGE,
    job_role_arn::AbstractString=STACK["TestBatchNodeJobRoleArn"],
)
    manager_code = """
        using AWSClusterManagers, Distributed, Memento
        setlevel!(getlogger("root"), "debug", recursive=true)

        addprocs(AWSBatchNodeManager())

        println("NumProcs: ", nprocs())
        for i in workers()
            println("Worker job \$i: ", remotecall_fetch(() -> ENV["AWS_BATCH_JOB_NODE_INDEX"], i))
        end
        """

    bind_to = "--bind-to \$(ip -o -4 addr list eth0 | awk '{print \$4}' | cut -d/ -f1)"
    worker_code = """
        using AWSClusterManagers
        start_batch_node_worker()
        """

    return Dict(
        "jobDefinitionName" => job_definition_name,
        "type" => "multinode",
        "nodeProperties" => Dict(
            "numNodes" => 3,
            "mainNode" => 0,
            "nodeRangeProperties" => [
                Dict(
                    "targetNodes" => "0",
                    "container" => Dict(
                        "image" => image,
                        "jobRoleArn" => job_role_arn,
                        "vcpus" => 1,
                        "memory" => 1024,  # MiB
                        "command" => ["julia", "-e", manager_code],
                    ),
                ),
                Dict(
                    "targetNodes" => "1:",
                    "container" => Dict(
                        "image" => image,
                        "jobRoleArn" => job_role_arn,
                        "vcpus" => 1,
                        "memory" => 1024,  # MiB
                        "command" => [
                            "bash",
                            "-c",
                            "julia $bind_to -e $(bash_quote(worker_code))",
                        ],
                    ),
                ),
            ],
        ),
        # Retrying to handle EC2 instance failures and internal AWS issues:
        # https://docs.aws.amazon.com/batch/latest/userguide/job_retries.html
        "retryStrategy" => Dict("attempts" => 3),
    )
end

# Use single quotes so that no shell interpolation occurs.
bash_quote(str::AbstractString) = string("'", replace(str, "'" => "'\\''"), "'")

# AWS Batch parallel multi-node jobs will only run on on-demand clusters. When running
# on spot the jobs will remain stuck in the RUNNABLE state
let ce = describe_compute_environment(STACK["ComputeEnvironmentArn"])
    if ce["computeResources"]["type"] != "EC2"  # on-demand
        error(
            "Aborting as compute environment $(STACK["ComputeEnvironmentArn"]) is not " *
            "using on-demand instances which are required for AWS Batch multi-node " *
            "parallel jobs.",
        )
    end
end

const BATCH_NODE_INDEX_REGEX = r"Worker job (?<worker_id>\d+): (?<node_index>\d+)"
const BATCH_NODE_JOB_DEF = register_job_definition(batch_node_job_definition())  # ARN

# Spawn all of the AWS Batch jobs at once in order to make online tests run faster. Each
# job spawned below has a corresponding testset. Ideally, the job spawning would be
# contained within the testset bun unfortunately that doesn't seem possible as `@sync` and
# `@testset` currently do not work together.
const BATCH_NODE_JOBS = Dict{String,BatchJob}()

let job_name = "test-worker-spawn-success"
    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF
    )
end

let job_name = "test-worker-spawn-failure"
    overrides = Dict(
        "numNodes" => 2,
        "nodePropertyOverrides" => [
            Dict(
                "targetNodes" => "1:",
                "containerOverrides" => Dict("command" => ["bash", "-c", "exit 0"]),
            ),
        ],
    )

    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF, node_overrides=overrides
    )
end

let job_name = "test-worker-link-local"
    overrides = Dict(
        "numNodes" => 2,
        "nodePropertyOverrides" => [
            Dict(
                "targetNodes" => "1:",
                "containerOverrides" => Dict(
                    "command" => [
                        "julia",
                        "-e",
                        """
                        using AWSClusterManagers, Memento
                        setlevel!(getlogger(), "debug", recursive=true)
                        try
                            start_batch_node_worker()
                        catch e   # Prevents the job from failing so we can retry AWS errors
                            showerror(stderr, e, catch_backtrace())
                        end
                        """,
                    ],
                ),
            ),
        ],
    )

    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF, node_overrides=overrides
    )
end

let job_name = "test-worker-link-local-bind-to"
    bind_to = "--bind-to \$(ip -o -4 addr list ecs-eth0 | awk '{print \$4}' | cut -d/ -f1)"
    worker_code = """
        using AWSClusterManagers, Memento
        setlevel!(getlogger("root"), "debug", recursive=true)
        try
            start_batch_node_worker()
        catch e   # Prevents the job from failing so we can retry AWS errors
            showerror(stderr, e, catch_backtrace())
        end
        """

    overrides = Dict(
        "numNodes" => 2,
        "nodePropertyOverrides" => [
            Dict(
                "targetNodes" => "1:",
                "containerOverrides" => Dict(
                    "command" => [
                        "bash",
                        "-c",
                        "julia $bind_to -e $(bash_quote(worker_code))",
                    ],
                ),
            ),
        ],
    )

    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF, node_overrides=overrides
    )
end

let job_name = "test-slow-manager"
    # Should match code in `batch_node_job_definition` but with an added delay
    manager_code = """
        using AWSClusterManagers, Distributed, Memento
        setlevel!(getlogger("root"), "debug", recursive=true)

        sleep(120)
        addprocs(AWSBatchNodeManager())

        println("NumProcs: ", nprocs())
        for i in workers()
            println("Worker job \$i: ", remotecall_fetch(() -> ENV["AWS_BATCH_JOB_NODE_INDEX"], i))
        end
        """

    overrides = Dict(
        "numNodes" => 2,
        "nodePropertyOverrides" => [
            Dict(
                "targetNodes" => "0",
                "containerOverrides" =>
                    Dict("command" => ["julia", "-e", manager_code]),
            ),
        ],
    )

    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF, node_overrides=overrides
    )
end

let job_name = "test-worker-timeout"
    # The default duration a worker process will wait for the manager to connect. For this
    # test this is also the amount of time between when the early worker checks in with the
    # manager and the late worker checks in.
    #
    # Note: Make sure to ignore any modification made on the local system that will not be
    # present for the batch job.
    worker_timeout = withenv("JULIA_WORKER_TIMEOUT" => nothing) do
        Distributed.worker_timeout()  # In seconds
    end

    # Amount of time it from the job start to executing `start_worker` (when the worker
    # timeout timer starts).
    start_delay = 10  # In seconds

    # Modify the manager to extend the worker check-in time. Ensures that the worker doesn't
    # timeout before the late worker checks in.
    manager_code = """
        using AWSClusterManagers, Dates, Distributed, Memento
        using AWSClusterManagers: AWS_BATCH_NODE_TIMEOUT
        setlevel!(getlogger(), "debug", recursive=true)

        check_in_timeout = Second(AWS_BATCH_NODE_TIMEOUT) + Second($(worker_timeout + start_delay))
        addprocs(AWSBatchNodeManager(timeout=check_in_timeout))

        println("NumProcs: ", nprocs())
        for i in workers()
            println("Worker job \$i: ", remotecall_fetch(() -> ENV["AWS_BATCH_JOB_NODE_INDEX"], i))
        end

        # Failure to launch all workers should trigger a retry via a non-zero exit code
        nworkers() == 2 || exit(2)
        """

    bind_to = "--bind-to \$(ip -o -4 addr list eth0 | awk '{print \$4}' | cut -d/ -f1)"

    # Note: The worker code logic tries to ensure that the execution of
    # `start_batch_node_worker` occurs around the same time for all workers.
    #
    # Requires that worker jobs have external network access and have permissions for IAM
    # access `batch:DescribeJobs`.
    worker_code = """
        using AWSBatch, AWSClusterManagers, Memento
        setlevel!(getlogger(), "debug")
        setlevel!(getlogger("AWSClusterManager"), "debug")
        node_index = parse(Int, ENV["AWS_BATCH_JOB_NODE_INDEX"])
        sibling_index = node_index % 2 + 1  # Assumes only 2 workers

        function wait_job_start(job::BatchJob)
            started_at = nothing
            while started_at === nothing
                started_at = get(describe(job), "startedAt", nothing)
                sleep(10)
            end
        end

        sibling_job_id = replace(
            ENV["AWS_BATCH_JOB_ID"],
            "#\$node_index" => "#\$sibling_index",
        )
        @info "Waiting for sibling job: \$sibling_job_id"
        wait_job_start(BatchJob(sibling_job_id))

        # Delaying a worker from reporting to the manager will delay the manager and cause
        # the workers that did report in to encounter the worker timeout. The delay here
        # should be less than the manager timeout to allow the delayed worker to still
        # report in.
        if node_index == 2  # The late worker will
            @info "Sleeping"
            sleep($(worker_timeout + start_delay))
        end

        start_batch_node_worker()
        """

    overrides = Dict(
        "numNodes" => 3,
        "nodePropertyOverrides" => [
            Dict(
                "targetNodes" => "0",
                "containerOverrides" =>
                    Dict("command" => ["julia", "-e", manager_code]),
            ),
            Dict(
                "targetNodes" => "1:",
                "containerOverrides" => Dict(
                    "command" => [
                        "bash",
                        "-c",
                        "julia $bind_to -e $(bash_quote(worker_code))",
                    ],
                ),
            ),
        ],
    )

    BATCH_NODE_JOBS[job_name] = submit_job(;
        job_name=job_name, job_definition=BATCH_NODE_JOB_DEF, node_overrides=overrides
    )
end

@testset "AWSBatchNodeManager (online)" begin
    # Note: Alternatively we could test report via Mocking but since the function is only
    # used for online testing and this particular test doesn't require an additional AWS
    # Batch job we'll test it here instead
    @testset "Report" begin
        job = BATCH_NODE_JOBS["test-worker-spawn-success"]

        manager_job = BatchJob(job.id * "#0")

        wait_finish(job)

        # Validate the report contains important information
        report_log = report(manager_job)
        test_results = [
            @test occursin(manager_job.id, report_log)
            @test occursin(string(status(manager_job)), report_log)
            @test occursin(status_reason(manager_job), report_log)
        ]

        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$report_log"
        end
    end

    @testset "Success" begin
        job = BATCH_NODE_JOBS["test-worker-spawn-success"]

        manager_job = BatchJob(job.id * "#0")
        worker_jobs = BatchJob.(job.id .* ("#1", "#2"))

        wait_finish(job)

        @test status(manager_job) == AWSBatch.SUCCEEDED
        @test all(status(w) == AWSBatch.SUCCEEDED for w in worker_jobs)

        # Expect 2 workers to check in and the worker ID order to match the node index order
        manager_log = log_messages(manager_job)
        matches = collect(eachmatch(BATCH_NODE_INDEX_REGEX, manager_log))
        test_results = [
            @test length(matches) == 2
            @test matches[1][:worker_id] == "2"
            @test matches[1][:node_index] == "1"
            @test matches[2][:worker_id] == "3"
            @test matches[2][:node_index] == "2"
        ]

        # Display the logs for all the jobs if any of the log tests fail
        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for worker 1:\n$(report(worker_jobs[1]))"
            @info "Details for worker 2:\n$(report(worker_jobs[2]))"
        end
    end

    @testset "Worker spawn failure" begin
        # Simulate a batch job which failed to start
        job = BATCH_NODE_JOBS["test-worker-spawn-failure"]

        manager_job = BatchJob(job.id * "#0")
        worker_job = BatchJob(job.id * "#1")

        wait_finish(job)

        # Even though the worker failed to spawn the cluster manager continues with the
        # subset of workers that reported in.
        @test status(manager_job) == AWSBatch.SUCCEEDED
        @test status(worker_job) == AWSBatch.SUCCEEDED

        manager_log = log_messages(manager_job)
        worker_log = log_messages(worker_job; retries=0)
        test_results = [
            @test occursin("Only 0 of the 1 workers job have reported in", manager_log)
            @test isempty(worker_log)
        ]

        # Display the logs for all the jobs if any of the log tests fail
        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for worker:\n$(report(worker_job))"
        end
    end

    @testset "Worker using link-local address" begin
        # Failing to specify a `--bind-to` address results in the link-local address being
        # reported from the workers which cannot be used by the manager to connect.
        job = BATCH_NODE_JOBS["test-worker-link-local"]

        manager_job = BatchJob(job.id * "#0")
        worker_job = BatchJob(job.id * "#1")

        wait_finish(job)

        @test status(manager_job) == AWSBatch.SUCCEEDED
        # Note: In practice worker jobs would actually fail but we catch the failure so that
        # we can retry the jobs for other AWS failure cases
        @test status(worker_job) == AWSBatch.SUCCEEDED

        manager_log = log_messages(manager_job)
        worker_log = log_messages(worker_job)
        test_results = [
            @test occursin("Only 0 of the 1 workers job have reported in", manager_log)
            @test occursin("Aborting due to use of link-local address", worker_log)
        ]

        # Display the logs for all the jobs if any of the log tests fail
        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for worker:\n$(report(worker_job))"
        end
    end

    @testset "Worker using link-local bind-to address" begin
        # Accidentially specifying the link-local address in `--bind-to`.
        job = BATCH_NODE_JOBS["test-worker-link-local-bind-to"]

        manager_job = BatchJob(job.id * "#0")
        worker_job = BatchJob(job.id * "#1")

        wait_finish(job)

        @test status(manager_job) == AWSBatch.SUCCEEDED
        # Note: In practice worker jobs would actually fail but we catch the failure so that
        # we can retry the jobs for other AWS failure cases
        @test status(worker_job) == AWSBatch.SUCCEEDED

        manager_log = log_messages(manager_job)
        worker_log = log_messages(worker_job)
        test_results = [
            @test occursin("Only 0 of the 1 workers job have reported in", manager_log)
            @test occursin("Aborting due to use of link-local address", worker_log)
        ]

        # Display the logs for all the jobs if any of the log tests fail
        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for worker:\n$(report(worker_job))"
        end
    end

    @testset "Worker connects before manager is ready" begin
        # If the workers manage to start and attempt to connect to the manager before the
        # manager is listening for connections the worker should attempt to reconnect.
        job = BATCH_NODE_JOBS["test-slow-manager"]

        manager_job = BatchJob(job.id * "#0")
        worker_job = BatchJob(job.id * "#1")

        wait_finish(job)

        @test status(manager_job) == AWSBatch.SUCCEEDED
        @test status(worker_job) == AWSBatch.SUCCEEDED

        manager_log = log_messages(manager_job)
        worker_log = log_messages(worker_job)

        test_results = [
            @test occursin("All workers have successfully reported in", manager_log)
        ]

        # Display the logs for all the jobs if any of the log tests fail
        # When the worker fails to connect the following error will occur:
        # `ERROR: IOError: connect: connection refused (ECONNREFUSED)`
        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for worker:\n$(report(worker_job))"
        end
    end

    @testset "Worker timeout" begin
        # In order to modify the Julia worker ordering the manager needs to wait for all
        # worker to check-in before the manager connects to any worker. Due to this
        # restriction it is possible that worker may wait longer than the
        # `Distributed.worker_timeout()` at which time the worker will self-terminate. If
        # this occurs during the cluster setup the manager node will hang and not proceed
        # past the initial cluster setup.
        #
        # If enough workers exist it is possible that the worker will self-terminate due to
        # hitting a timeout (The issue has been seen with as low as 30 nodes with a timeout
        # of 60 seconds). We typically shouldn't encounter this due to increasing the worker
        # timeout.
        #
        # Log example:
        # ```
        # [debug | AWSClusterManagers]: Worker connected from node B
        # [debug | AWSClusterManagers]: Worker connected from node A
        # ...
        # [debug | AWSClusterManagers]: All workers have successfully reported in
        # Worker X terminated.
        # Worker Y terminated.
        # ...
        # <manager stalled>
        # ```

        job = BATCH_NODE_JOBS["test-worker-timeout"]

        manager_job = BatchJob(job.id * "#0")
        early_worker_job = BatchJob(job.id * "#1")
        late_worker_job = BatchJob(job.id * "#2")

        wait_finish(job)

        test_results = [
            @test status(manager_job) == AWSBatch.SUCCEEDED
            @test status(early_worker_job) == AWSBatch.SUCCEEDED
            @test status(late_worker_job) == AWSBatch.SUCCEEDED
        ]

        if any(r -> !(r isa Test.Pass), test_results)
            @info "Details for manager:\n$(report(manager_job))"
            @info "Details for early worker:\n$(report(early_worker_job))"
            @info "Details for late worker:\n$(report(late_worker_job))"
        end
    end
end


text:
#!/usr/bin/env julia

const GLOBALS = Dict(:empCount=>0)

struct Employee
    name::String
    salary::Float64

    function Employee(name::AbstractString, sal::Real)
        GLOBALS[:empCount] += 1
        new(name, sal)
    end
end

display_employee_count() = println("Total Employees: $(GLOBALS[:empCount])")
# could also  do
# display_count(::Type{Employee})

function Base.show(io::IO, emp::Employee)
    println(io, "Employee (Name: $(emp.name), Salary: $(emp.salary))")
end

emp1 = Employee("Zara", 2000)
emp2 = Employee("Manni", 5000)

show(emp1)
show(emp2)

display_employee_count()

ext: jl

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
### using DArrays ###

function hpl_par2(A::Matrix, b::Vector)
    n = size(A,1)
    A = [A b]

    C = distribute(A, 2)
    nB = length(C.pmap)

    ## case if only one processor
    if nB <= 1
        x = A[1:n, 1:n] \ A[:,n+1]
        return x
    end

    depend = Array(RemoteRef, nB, nB)

    #pmap[i] is where block i's stuff is
    #block i is dist[i] to dist[i+1]-1
    for i = 1:nB
        #println("C=$(convert(Array, C))") #####
        ##panel factorization
        panel_p = remotecall_fetch(C.pmap[i], panel_factor_par2, C, i, n)

        ## Apply permutation from pivoting
        for j = (i+1):nB
           depend[i,j] = remotecall(C.pmap[j], permute, C, i, j, panel_p, n, false)
        end
        ## Special case for last column
        if i == nB
           depend[nB,nB] = remotecall(C.pmap[nB], permute, C, i, nB+1, panel_p, n, true)
        end

        ##Trailing updates
        (i == nB) ? (I = (C.dist[i]):n) :
                    (I = (C.dist[i]):(C.dist[i+1]-1))
        C_II = C[I,I]
        L_II = tril(C_II, -1) + eye(length(I))
        K = (I[length(I)]+1):n
        if length(K) > 0
            C_KI = C[K,I]
        else
            C_KI = zeros(0)
        end

        for j=(i+1):nB
            dep = depend[i,j]
            depend[j,i] = remotecall(C.pmap[j], trailing_update_par2, C, L_II, C_KI, i, j, n, false, dep)
        end

        ## Special case for last column
        if i == nB
            dep = depend[nB,nB]
            remotecall_fetch(C.pmap[nB], trailing_update_par2, C, L_II, C_KI, i, nB+1, n, true, dep)
        else
            #enforce dependencies for nonspecial case
            for j=(i+1):nB
                wait(depend[j,i])
            end
        end
    end

    A = convert(Array, C)
    x = triu(A[1:n,1:n]) \ A[:,n+1]
end ## hpl_par2()

function panel_factor_par2(C, i, n)
    (C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
                           (I = (C.dist[i]):(C.dist[i+1]-1))
    K = I[1]:n
    C_KI = C[K,I]
    #(C_KI, panel_p) = lu!(C_KI) #economy mode
    panel_p = lu!(C_KI)[2]
    C[K,I] = C_KI

    return panel_p
end ##panel_factor_par2()

function permute(C, i, j, panel_p, n, flag)
    if flag
        K = (C.dist[i]):n
        J = (n+1):(n+1)
        C_KJ = C[K,J]

        C_KJ = C_KJ[panel_p,:]
        C[K,J] = C_KJ
    else
        K = (C.dist[i]):n
        J = (C.dist[j]):(C.dist[j+1]-1)
        C_KJ = C[K,J]

        C_KJ = C_KJ[panel_p,:]
        C[K,J] = C_KJ
    end
end ##permute()

function trailing_update_par2(C, L_II, C_KI, i, j, n, flag, dep)
    if isa(dep, RemoteRef); wait(dep); end
    if flag
        #(C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
        #                       (I = (C.dist[i]):(C.dist[i+1]-1))
        I = C.dist[i]:n
        J = (n+1):(n+1)
        K = (I[length(I)]+1):n
        C_IJ = C[I,J]
        if length(K) > 0
            C_KJ = C[K,J]
        else
            C_KJ = zeros(0)
        end
        ## Compute blocks of U
        C_IJ = L_II \ C_IJ
        C[I,J] = C_IJ
    else
        #(C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
        #                       (I = (C.dist[i]):(C.dist[i+1]-1))

        I = (C.dist[i]):(C.dist[i+1]-1)
        J = (C.dist[j]):(C.dist[j+1]-1)
        K = (I[length(I)]+1):n
        C_IJ = C[I,J]
        if length(K) > 0
            C_KJ = C[K,J]
        else
            C_KJ = zeros(0)
        end

        ## Compute blocks of U
        C_IJ = L_II \ C_IJ
        C[I,J] = C_IJ
        ## Trailing submatrix update - All flops are here
        if !isempty(C_KJ)
            cm, ck = size(C_KI)
            cn = size(C_IJ,2)
            blas_gemm('N','N',cm,cn,ck,-1.0,C_KI,cm,C_IJ,ck,1.0,C_KJ,cm)
            #C_KJ = C_KJ - C_KI*C_IJ
            C[K,J] = C_KJ
        end
    end
end ## trailing_update_par2()

## Test n*n matrix on np processors
## Prints 5 numbers that should be close to zero
function test(n, np)
    A = rand(n,n); b = rand(n);
    X = (@elapsed x = A \ b);
    Y = (@elapsed y = hpl_par(A,b, max(1,div(n,np))));
    Z = (@elapsed z = hpl_par2(A,b));
    for i=1:(min(5,n))
        print(z[i]-y[i], " ")
    end
    println()
    return (X,Y,Z)
end

## test k times and collect average
function test(n,np,k)
    sum1 = 0; sum2 = 0; sum3 = 0;
    for i = 1:k
        (X,Y,Z) = test(n,np)
        sum1 += X
        sum2 += Y
        sum3 += Z
    end
    return (sum1/k, sum2/k, sum3/k)
end


text:
"""
    Code for a super-Gaussian mask shape for use with CCF
Author: Eric Ford
Created: September 2020
"""

using QuadGK

"""   SuperGaussianCCFMask
A truncated Gaussian mask with two parameters, it's stdandard deviation and where to truncate it, both as a velocity in m/s.
Mask weights are stored separately in a line list.

TODO: Repalce Gaussian with super-Gaussian
Warning:  Not implemented/tested yet.
"""
struct SuperGaussianCCFMask <: AbstractCCFMaskShape
    σ_sqrt2::Float64
    power::Float64
    half_width_truncation::Float64
    normalization::Float64

    """ `SuperGaussianCCFMask( σ ; half_truncation_width_in_σ=2 )` """
    function SuperGaussianCCFMask(σ::Real, p::Real, w::Real=2 )
        @assert 0 < σ <= 300000   # 300km/s is arbitrary choice for an upper limit
        @assert 1 <= p <= 2
        @assert 0 < w <= 4
        function integrand(x)
            exp(-(0.5*(x/σ)^2)^p)
        end
        integral = quadgk(integrand, -w*σ, w*σ)[1]
        norm = 1.0/integral
        new(σ*sqrt(2.0),p,σ*w,norm)
    end
end

""" `SuperGaussianCCFMask( inst ; scale_factor )` """
function SuperGaussianCCFMask(inst::InstT; power::Real = default_supergaussian_ccf_exponent, fwhm::Real = default_supergaussian_ccf_fwhm,
        σ_scale_factor::Real = 1, truncation_scale_factor::Real = default_supergaussian_ccf_truncation_scale_factor ) where { InstT<:AbstractInstrument }
    # TODO: Update default values
    σ = σ_scale_factor * fwhm/sqrt(8 * log(2)^(1/power))  # From Ryan Petersburg email 9/11/2020
    w = truncation_scale_factor
    SuperGaussianCCFMask(σ,power,w/2)
end

λ_min(m::SuperGaussianCCFMask,λ::Real) = λ/calc_doppler_factor(m.half_width_truncation)
λ_max(m::SuperGaussianCCFMask,λ::Real) = λ*calc_doppler_factor(m.half_width_truncation)

function integrate(m::SuperGaussianCCFMask, v_lo::Real,v_hi::Real)
    quadgk(m, v_lo, v_hi, atol=1e-4)[1]
end

""" Functor for returning PSF for Δv <= half_width.  """
function (m::SuperGaussianCCFMask)(Δv::Real)
    if abs2(Δv) > abs2(m.half_width_truncation)
        return zero(Δv)
    else
        return m.normalization*exp(-((Δv/m.σ_sqrt2)^2)^m.power)
    end
end

function mask_with_increased_fwhm(m::SuperGaussianCCFMask, Δfwhm::Real )
    # WARNING: Haven't checked how good of an approximation this is
    σ = m.σ_sqrt2/sqrt(2)
    fwhm_orig = σ * sqrt(8 * log(2)^(1/m.power))
    fwhm_new = sqrt(fwhm_orig^2 + Δfwhm^2)
    σ_new = fwhm_new/sqrt(8 * log(2)^(1/m.power))
    return SuperGaussianCCFMask(σ_new, m.power, m.half_width_truncation * (σ_new * sqrt(2) / σ_sqrt2) )
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
__precompile__()
module COSMO

using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg

export assemble!, warmStart!, empty_model!

const DefaultFloat = Float64
const DefaultInt   = LinearAlgebra.BlasInt

include("./kktsolver.jl")
# optional dependencies
if in("Pardiso",keys(Pkg.installed()))
    include("./kktsolver_pardiso.jl")
end

include("./algebra.jl")
include("./projections.jl")
include("./settings.jl")            # TODO: unmodified - revisit
include("./trees.jl")
include("./types.jl")               # some types still need tidying
include("./constraint.jl")          # TODO: unmodified - revisit
include("./parameters.jl")          # TODO: unmodified - revisit
include("./residuals.jl")           # TODO: unmodified - revisit
include("./scaling.jl")             # TODO: set scaling / E scaling is broken
include("./infeasibility.jl")       # TODO: stylistic fixes needed
include("./chordal_decomposition.jl")
include("./printing.jl")            # TODO: unmodified - revisit
include("./setup.jl")               # TODO: unmodified - revisit (short - consolidate?)
include("./solver.jl")              # TODO: unmodified - revisit
include("./interface.jl")           # TODO: unmodified - revisit
include("./MOIWrapper.jl")

end #end module


text:
include("../testpolymodule.jl")

using DynamicPolynomials
@polyvar x

using SemialgebraicSets
using PolyJuMP
const NonNeg = TestPolyModule.NonNeg{MonomialBasis, typeof(@set x^2 ≤ 0),
                                     monomialtype(x), monovectype(x)}

MOIU.@model(
    PolyNonNegModel,
    (), (MOI.LessThan,), (NonNeg,), (),
    (), (), (MOI.VectorOfVariables,), (MOI.VectorAffineFunction,))

config = MOI.Test.TestConfig()
optimize!(mock) = MOIU.mock_optimize!(mock, [1.0, 1.0],
    (MOI.VectorOfVariables, NonNeg) => [[0.0]],
    (MOI.VectorAffineFunction{Float64}, NonNeg) => [[0.0, -0.5],
                                                    [0.0, 0.5],
                                                    [0.0]])
mock = bridged_mock(optimize!; model = PolyNonNegModel{Float64}())
Tests.plus_minus_test(mock, config; polymodule = TestPolyModule)


text:
# import VideoIO
import FFMPEG

export record

function record(f_canvas, filename, timestamps;
        fps = 60, px_per_pt = 1, ffmpeg_options = [])

    nframes = length(timestamps)
    bname = basename(filename)

    mktempdir() do dir
        for (i, t) in enumerate(timestamps)
            canv = f_canvas(t)
            png(canv, joinpath(dir, "$(lpad(i, 5, '0')).png"); px_per_pt = 1)
        end

        FFMPEG.exe(` -i $(joinpath(dir, "%5d.png")) -framerate $fps $ffmpeg_options $(joinpath(dir, bname))`)
        cp(joinpath(dir, bname), filename, force = true)
    end
end


text:
using PEG

function do_op(x)
  local v
  local op_v_pairs
  v, op_v_pairs = x
  for op_v ∈ op_v_pairs
    op, v2 = op_v
    v = eval(Symbol(op))(v, v2)
  end
  v
end

@rule num = r"\d+"w |> x->parse(Int, x) , r"\("p & add & r"\)"p |> x->x[2]
@rule mul = num & (r"[/*]"p & num)[*] |> do_op
@rule add = mul & (r"[+-]"p & mul)[*] |> do_op

using Test
@test add("2 + 3 * 5") == (17, "")
@test add("2 * 3 + 5") == (11, "")
@test add("(2 + 3) * 5") == (25, "")

@test try
  parse_whole(add, "1 + 2 ^ 3")
catch x
  split(x.msg, "\n")[1]
end == "On line 1, at column 7 (byte 7):"


text:
module JSONSerializer

using JSON

function serialize(data)::String
  JSON.json(data)
end

function deserialize(json)
  JSON.parse(json)
end

end


text:
@with_kw struct FastTree
    tolerance::Float64 = 1.0
end

# This is the main function

function solve(solver::FastTree, problem::Problem)
    center = problem.input.center
    radius = problem.input.radius[1]
    (W, b) = (problem.network.layers[1].weights, problem.network.layers[1].bias)

    input = forward_affine_map(solver, W, b, problem.input)
    lower, upper = low(input), high(input)
    local_lower, local_upper = similar(lower), similar(lower)

    np = pyimport("numpy")
    q, e, p = np.linalg.svd(W)

    k_1 = size(W, 1)
    k_0 = size(W, 2)

    C = vcat(W, -W)
    d = zeros(2k_1)

    hps = Array{Hyperplane, 1}(undef, k_0-k_1)

    for i in k_1+1:k_0
        hps[i-k_1] = Hyperplane(p[i,:], dot(p[i,:], center))
    end

    S = hps[1]
    if k_0-k_1 > 1
        for i in 2:k_0-k_1
            S = intersection(S, hps[i])
        end
    end

    stack = Vector{Hyperrectangle}(undef, 0)
    push!(stack, input)
    #count = 1
    while !isempty(stack)
        interval = popfirst!(stack)
        reach = forward_network(solver, problem.network, interval)
        if issubset(reach, problem.output)
            continue
        else
            if get_largest_width(interval) > solver.tolerance
                sections = bisect(interval)
                for i in 1:2
                    local_lower = low(sections[i])
                    local_upper = high(sections[i])

                    d = vcat(local_upper - b, b - local_lower)
                    P_i = HPolyhedron(C, d)

                    inter = intersection(problem.input, P_i)

                    if isempty(inter) == false
                        O_i = box_approximation(intersection(P_i, S))
                        inner = true
                        for j in 1:k_0
                            if low(O_i)[j] ≤ low(problem.input)[j]
                                inner = false
                                break
                            end

                            if high(problem.input)[j] ≤ high(O_i)[j]
                                inner = false
                                break
                            end
                        end
                        if inner == false
                            push!(stack, sections[i])
                            #count += 1
                        end
                    end
                end
            else
                return BasicResult(:violated)
            end
        end
    end
    #print("\n$(count)\n")
    return BasicResult(:holds)
end

function forward_network(solver::FastTree, nnet::Network, input::Hyperrectangle)
    layers = nnet.layers
    act = layers[1].activation
    reach = Hyperrectangle(low = act.(low(input)), high = act.(high(input)))

    for i in 2:length(layers)
        reach = forward_layer(solver, layers[i], reach)
    end
    return reach
end

function forward_layer(solver::FastTree, L::Layer, input::Hyperrectangle)
    (W, b, act) = (L.weights, L.bias, L.activation)
    center = zeros(size(W, 1))
    gamma  = zeros(size(W, 1))
    for j in 1:size(W, 1)
        node = Node(W[j,:], b[j], act)
        center[j], gamma[j] = forward_node(solver, node, input)
    end
    return Hyperrectangle(center, gamma)
end

function forward_node(solver::FastTree, node::Node, input::Hyperrectangle)
    output    = node.w' * input.center + node.b
    deviation = sum(abs.(node.w) .* input.radius)
    βmax = node.act(output + deviation)
    βmin = node.act(output - deviation)
    return ((βmax + βmin)/2, (βmax - βmin)/2)
end

function forward_affine_map(solver::FastTree, W::Matrix, b::Vector, input::Hyperrectangle)
    center = W * input.center + b
    radius = abs.(W) * input.radius
    return Hyperrectangle(center, radius)
end

#to determine whether
function ishull(x::Hyperrectangle, y::Hyperrectangle, W::Matrix, b::Vector)
    center = y.center
    radius = y.radius[1]
    k_1 = size(W, 1) #length(lower)
    #I = zeros(k_1, k_1) #Identity matrix
    #IN = zeros(k_1, k_1) #negative Identity matrix
    #for k = 1:k_1
        #I[k, k] = 1.0
        #IN[k, k] = -1.0
    #end
    #Hi
    #C_1 = vcat(I, IN)
    C_1 = vcat(Diagonal(ones(k_1)), -Diagonal(ones(k_1)))
    d_1 = vcat(high(x), -low(x))
    #Pi
    C = C_1 * W #size(2k_1, k_0)
    d = d_1 - C_1 * b #size(2k_1)

    #print("C: ")
    #println(C)
    #print("d: ")
    #println(d)

    inter = true #if the intersection of Pi and the input set is empty
    for l = 1:length(d)
        sum = 0.0
        for m = 1:size(C, 2)
            if C[l,m] > 0
                sum += low(y, m) * C[l,m]
            elseif C[l,m] < 0.0
                sum += high(y, m) * C[l,m]
            end
        end
        if sum > d[l]
            inter = false
            break
        end
    end

    if inter
        hull = false #if Hi is a hull of Z, or if Pi and center of the input set are close
        for l = 1:length(d)
            #print(l)
            #print(" distance: ")
            #println(distance(center, C[l,:], d[l]))
            #print("radius: ")
            #println(radius)
            if distance(center, C[l,:], d[l]) > radius
                return true
            end
        end
    end
    return false
end


text:
"""
    Recorder

Struct used to record tracking results to a video.

# Keyword Arguments:
- `filename = "trackingresult.mp4"`
- `framerate = 30`
"""
Base.@kwdef mutable struct Recorder
    filename = "trackingresult.mp4"
    framerate = 30
    encoder = nothing
    index::Int = 1
end

function Recorder(img; kwargs...)
    encoder_options = (crf=22, preset="medium")
    filename = "trackingresult.mp4"
    Recorder(;encoder=open_video_out(filename, img; framerate=framerate, encoder_options=encoder_options), kwargs...)
end

"""
    showblobs(img::AbstractMatrix{T}, result::TrackingResult, m::Measurement; rad=8, recorder=nothing, display=true) where T

Overlay found blobs on `img`

#Arguments:
- `img`: an image
- `result`: a `TrackingResult`
- `m`: a `Measurement`
- `rad`: radius of blobs to draw
- `recorder`: an optional `Recorder`
- `display = Base.display`: function to display image. Use `display=nothing` to not display.
"""
function showblobs(img::AbstractMatrix{T},result,m;rad=8, recorder=nothing, display=Base.display, ignoreempty=false) where T
    # img = copy(img)
    display === nothing && recorder === nothing && return
    if ignoreempty
        isempty(result.blobs) && isempty(m.coordinates) && return
    end
    blobs = result.blobs
    foreach(blobs) do blob
        blobcoord = location(blob)
        c = length(blob.trace) == 1 ? RGB{eltype(T)}(0.,1.,0.) : RGB{eltype(T)}(0.,0.,1.)
        draw!(img, ImageDraw.CirclePointRadius(blobcoord, rad),c)
        mcoord = blob.trace[end]
        mcoord == OOB && return
        draw!(img, ImageDraw.LineSegment(blobcoord, mcoord))
    end
    foreach(enumerate(m.coordinates)) do (ci,coord)
        c = RGB{eltype(T)}(1.,0.,0.) # ci ∈ newcoordinds ? RGB{T}(0,1,0) : RGB{T}(1,0,0)
        draw!(img, ImageDraw.CirclePointRadius(coord, rad/2), c)
    end
    record(img, recorder)
    if display !== nothing
        display(img)
    end
    img
end

record(_, ::Nothing) = nothing

function record(img, recorder)
    r = recorder
    if r.encoder === nothing
        encoder_options = (crf=22, preset="medium")
        r.encoder = open_video_out(recorder.filename, img, framerate=r.framerate, encoder_options=encoder_options)
    end
    VideoIO.write(r.encoder, img, r.index)
    r.index += 1
end

Base.finalize(::Nothing) = nothing
function Base.finalize(r::Recorder)
    close_video_out!(r.encoder)
end

ImageDraw.draw!(img,b::Images.BlobLoG;kwargs...) = draw!(img,location(b);kwargs...)
ImageDraw.draw!(img,coord::CartesianIndex;a=3,c=RGB(1.,0.,0.)) = draw!(img, CirclePointRadius(coord, float(a)), eltype(img)(c))

function ImageDraw.draw!(img, trace::Trace, args...; c=RGB(1,1,1), kwargs...)
    trace = skipmissing(trace)
    draw!(img, trace[1], c=RGB(0,1,0))
    draw!(img, trace[end], c=RGB(1,0,0))
    draw!(img, ImageDraw.Path(trace), eltype(img)(c), args...; kwargs...)
end

function ImageDraw.draw!(img, traces::Vector{Trace}, args...; kwargs...)
    for trace in traces
        ImageDraw.draw!(img, trace, args...; kwargs...)
    end
    img
end


text:
abstract ScheduledEvent

immutable StopEvent <: ScheduledEvent
	time::TimeOfDay
end

# this encompasses everything we need to know about when/how to do the scheduled event
immutable NormalEvent{T} <: ScheduledEvent
	time::TimeOfDay   				# when this should run
	pub::Publisher            # the publisher to be scheduled (which contains both the function call and receivers/listeners)
	args::T 											# the arguments to the function
end

# Base.string(event::ScheduledEvent) = (isdefined(event,:pub) ? "Event{time=$(event.time), $(event.pub), $(event.args)}" : "Event{time=$(event.time)}")
# @createIoMethods ScheduledEvent
# Base.print(io::IO, event::ScheduledEvent) = print(io, string(event))
# Base.show(io::IO, event::ScheduledEvent) = print(io, string(event))

compareEventTimes(a1::ScheduledEvent, a2::ScheduledEvent) = return (a1.time < a2.time ? -1 : (a1.time == a2.time ? 0 : 1))

# --------------------------------------------------------------

# this type represents something that can handle scheduling publishing sometime in the future
# it will be the "main loop".  In simulation, it will keep track of "simulation time" and
# properly order events while updating the time.
# In production/live scheduling, it will use clock time and sleep/wake at appropriate times, while
# ideally polling various sockets for market data, orders, etc
abstract Scheduler

# some helper functions
isQueueEmpty(scheduler::Scheduler) = isempty(scheduler.eventList)
firstTimeInQueue(scheduler::Scheduler) = peek(scheduler.eventList).time

# --------------------------------------------------------------

type LiveScheduler
	timeOfLastEvent::TimeOfDay
	eventList::SinglyLinkedList{ScheduledEvent}
	LiveScheduler() = (s = new(); empty!(s); s)
end

NOW(scheduler::LiveScheduler) = currentTimeOfDay()

function Base.empty!(scheduler::LiveScheduler)
	scheduler.timeOfLastEvent = 0
	scheduler.eventList = SinglyLinkedList(ScheduledEvent)
end

# TODO:
function publish(scheduler::LiveScheduler, event::ScheduledEvent)
	return
end
function publish(scheduler::LiveScheduler, time::TimeOfDay, pub::Publisher, args...)
	return
end

# --------------------------------------------------------------

type SimulationScheduler <: Scheduler
	timeOfLastEvent::TimeOfDay
	eventList::SinglyLinkedList{ScheduledEvent}
	SimulationScheduler() = (s = new(); empty!(s); s)
end

NOW(scheduler::SimulationScheduler) = scheduler.timeOfLastEvent

function Base.empty!(scheduler::SimulationScheduler)
	scheduler.timeOfLastEvent = 0
	scheduler.eventList = SinglyLinkedList(ScheduledEvent)
end

function publish(scheduler::SimulationScheduler, event::ScheduledEvent)
	scheduler.timeOfLastEvent = event.time
	publish(event.pub, event.args...)
	return
end

function publish(scheduler::SimulationScheduler, time::TimeOfDay, pub::Publisher, args...)
	scheduler.timeOfLastEvent = time
	publish(pub, args...)
	return
end

# --------------------------------------------------------------

# this is a global constant (also a singleton), so that schedule can easily be called from anywhere without passing
# around the SCHEDULER object
const SCHEDULER = SimulationScheduler()  # TODO: allow choice between sim/live somehow

NOW() = NOW(SCHEDULER)
initScheduler() = empty!(SCHEDULER)
stopScheduler(time::TimeOfDay = zero(TimeOfDay)) = (insertSorted!(SCHEDULER.eventList, StopEvent(time), compareEventTimes); return)

# --------------------------------------------------------------

# returns true when we should stop
processEvent(event::NormalEvent) = (publish(SCHEDULER, event); false)
processEvent(event::StopEvent) = (SCHEDULER.timeOfLastEvent = event.time; true)

function process_until(time::TimeOfDay)
	# clear out the queue until the current time
	while !isQueueEmpty(SCHEDULER) && time >= firstTimeInQueue(SCHEDULER)
		nextevent = shift!(SCHEDULER.eventList)  # gets the first in the list
		processEvent(nextevent) && return true # stopped... return true
	end
	return false
end

# these functions should be called from non-main-loop areas... schedule the event, and don't worry about how it gets run
function Base.schedule(time::TimeOfDay, pub::Publisher, args...)
	event = NormalEvent(time, pub, args)
	insertSorted!(SCHEDULER.eventList, event, compareEventTimes)
	event
end

Base.schedule(time::TimeOfDay, cb::Callback, args...) = schedule(time, cb.pub, args...)

# call this method from the "main loop"... it either:
#		1) runs the function immediately, or
#		2) schedules the event and clears the queue up to and including that event
# note: returns true if we hit a StopEvent
# note: this function is really intended for the situation that a stream of events is "driving time forward",
#				such as historical market data.  if there are no other events waiting, just keep processing the
#				main event stream as quickly as possible
function schedule_do(time::TimeOfDay, pub::Publisher, args...)

	# # clear out the queue until the current time
	# while !isQueueEmpty(SCHEDULER) && time >= firstTimeInQueue(SCHEDULER)
	# 	nextevent = shift!(SCHEDULER.eventList)  # gets the first in the list
	# 	processEvent(nextevent) && return true # stopped... return true
	# end
	process_until(time) && return true

	# now publish it
	publish(SCHEDULER, time, pub, args...)
	false
end

schedule_do(time::TimeOfDay, cb::Callback, args...) = schedule_do(time, cb.pub, args...)

# --------------------------------------------------------------

# handle scheduling the same thing multiple times

function Base.schedule(times::AbstractVector{TimeOfDay}, args...)
	for t in times
		schedule(t, args...)
	end
end

function schedule_do(times::AbstractVector{TimeOfDay}, args...)
	for t in times
		schedule_do(t, args...)
	end
end

# --------------------------------------------------------------

# process events in order until either the queue is empty (returns false) or
# we hit a StopEvent (returns true)
function processAllEvents()
	while !isQueueEmpty(SCHEDULER)
		nextevent = shift!(SCHEDULER.eventList)  # gets the first in the list
		processEvent(nextevent) && return true # stopped... return true
	end
	return false
end


text:
const ymeyer = [3.478e4, 2.861e4, 2.365e4, 1.963e4, 1.637e4, 1.372e4,
                1.154e4, 9.744e3, 8.261e3, 7.030e3, 6.005e3, 5.147e3,
                4.427e3, 3.820e3, 3.307e3, 2.872e3]

const meyer = let res_init=zeros(16), jac_init=zeros(16,3), x_init=[0.02, 4.e3, 2.5e2]

    function res(x, r)
        for i = 1:16
            ti = Float64(5*i + 45)
            r[i] = x[1] * exp(x[2] / (ti + x[3])) - ymeyer[i]
        end
        return r
    end

    function jac(x, j)
        for i = 1:16
            ti = Float64(5*i + 45)
            u = ti + x[3]
            t = exp(x[2] / u)
            j[i, 1] = t
            j[i, 2] = x[1] * t/u
            j[i, 3] = -x[1] * x[2] * t / (u * u)
        end
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("meyer", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
# new_LP_test.jl
# New LP formulation for travel time estimation - test for metropolis
# Authored by Arthur J Delarue on 7/9/15

# VARMAP = {
# 	:Basic => 0,
# 	:Superbasic => -3,
# 	:NonbasicAtUpper => -2,
# 	:NonbasicAtLower => -1
# }

# PROB = 0.7
# MODEL = "metropolis_$(PROB)_relaxed"
# MAX_ROUNDS = 54   #DON'T CHOOSE 30 or 55
# MIN_RIDES = 1

# TURN_COST = 2.0

# LAMBDA = [1e3]#[1e9,1e3,2e3,5e3,1e4,2e4,5e4,1e5,2e5,5e5,1e6,2e6,5e6,1e25]
# DELTA_BOUND = [0.6,0.5,0.4]

function new_LP(
	graph::SimpleGraph,
	travelTimes::Array{Float64,2},
	numRides::Array{Int,2},
	roadTimes::AbstractArray{Float64, 2},
	distances::AbstractArray{Float64, 2},
	positions::Array{Coordinates};
	model_type::AbstractString=MODEL,
	max_rounds::Int=MAX_ROUNDS,
	min_rides::Int=MIN_RIDES,
	turnCost::Float64=TURN_COST,
	delta_bound::Array{Float64}=DELTA_BOUND)
	"""
	Runs our iterative LP formulation on a Metropolis object
	"""

	roads = edges(graph)
	nodes = vertices(graph)
	out = [copy(out_neighbors(graph,i)) for i in nodes]
	inn = [copy(in_neighbors(graph,i)) for i in nodes]

	TMAX = 1000
	TESTDIR = "minr$(min_rides)_i$(max_rounds)_$(model_type)"

	# Create directory if necessary:
	if !isdir("Outputs/$TESTDIR")
		mkdir("Outputs/$TESTDIR")
	end

	# Create output csv file to save algorithm data
	println("-- Saving outputs to directory Outputs/$(TESTDIR)/")

	# Create JuMP model for LP and QP
	println("**** Creating LP instance ****")
	m = Model(solver=GurobiSolver(TimeLimit=10000, Method=2, Crossover=1, OutputFlag=0))
	# m2 = Model(solver=GurobiSolver(TimeLimit=10000, Method=1, OutputFlag=1, InfUnbdInfo=1))

	# Add one variable for each road
	@defVar(m, t[i=nodes,j=out[i]] >= roadTimes[i,j])

	# Add decision variables for first LP
	pairs = [find(travel_times[i,:]) for i=nodes]
	@defVar(m, epsilon[i=nodes,j=pairs[i]] >= 0)
	@defVar(m, T[i=nodes, j=pairs[i]] >= 0)
	@addConstraint(m, TLowerBound[i=nodes,j=pairs[i]], T[i,j] - travelTimes[i,j] >= - epsilon[i,j])
	@addConstraint(m, TUpperBound[i=nodes,j=pairs[i]], T[i,j] - travelTimes[i,j] <= epsilon[i,j])

	# Define objective variables and constraints for m2
	@defVar(m, delta2[i=nodes,j=out[i]] >= 0)
	@addConstraint(m, objConstrLower[i=nodes,j=out[i]], -1 * t[i,j]/distances[i,j] + 1/(length(inn[i]) + length(out[j])) * (sum{1/distances[j,k] * t[j,k], k = out[j]} + sum{1/distances[h,i] * t[h,i], h=inn[i]}) <= delta2[i,j])
	@addConstraint(m, objConstrUpper[i=nodes,j=out[i]], t[i,j]/distances[i,j] - 1/(length(inn[i]) + length(out[j])) * (sum{1/distances[j,k] * t[j,k], k = out[j]} + sum{1/distances[h,i] * t[h,i], h=inn[i]}) <= delta2[i,j])

	# Create handles for all constraints
	numDataPoints = sum([length(pairs[i]) for i=nodes])
	@defConstrRef path[1:numDataPoints, 1:max_rounds]

	# Create arrays to keep track of constraints that are already in the set of constraints
	# Map from hashed paths to constraint indices
	hashedPathIndices = fill((Uint64 => Int64)[], numDataPoints)
	# Set of hashed paths for quick lookup
	hashedPaths = fill(Set(Uint64[]), numDataPoints)
	# Keep track of last path
	previousPaths = [0 for i = 1:numDataPoints]

	status = 0
	newTimes = roadTimes
	old_objective = 0
	totalPathConstraints = 0

	# Compute shortest paths (with turn cost)
	println("**** Computing shortest paths ****")
	@time new_graph, new_edge_dists, new_nodes = modifyGraphForDijkstra(graph, newTimes, positions, turn_cost=turnCost)
	old_nodes = getInverseMapping(new_nodes, nv(new_graph))
	@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)

	# Run over all pairs of nodes that have data
	srcs = Int[]
	dsts = Int[]
	sizehint(srcs, 12000)
	sizehint(dsts, 12000)
	for i in nodes, j in nodes
		if travelTimes[i,j] > 0
			# Load info to reconstruct shortest paths in parallel
			push!(srcs, i)
			push!(dsts, j)
		end
	end
	totalNumExpensiveTurns = Array(Int, (numDataPoints, max_rounds))

	# Create variable bound arrays for MathProgBase manual updating
	variableLowerBounds = zeros(2 * numDataPoints + 2 * length(roads))
	variableUpperBounds = zeros(2 * numDataPoints + 2 * length(roads))

	l = 1
	while l <= max_rounds
		actual_l = l
		if split(model_type, "_")[end] == "relaxed"
			delta = (1/2)^(div(actual_l-1,length(delta_bound)))*delta_bound[((actual_l-1)%length(delta_bound)) + 1]
			# Avoid numerical issues
			if delta < 1e-4
				delta = 0
			end
		end
		setSolver(m, GurobiSolver(TimeLimit=10000, Method=2, Crossover=1, OutputFlag = 0))
		@setObjective(m, Min, sum{ sqrt(numRides[i,j]/travel_times[i,j]) * epsilon[i,j], i=nodes, j=pairs[i]})

		println("###### ROUND $l ######")

		# Create upper and lower bound arrays for MathProgBase manual constraint management
		pathLowerBounds = zeros(totalPathConstraints + 2 * numDataPoints + 2 * length(roads))
		pathUpperBounds = zeros(totalPathConstraints + 2 * numDataPoints + 2 * length(roads))

		# Fill in upper and lower bounds for non-path constraints
		for i = nodes, j = pairs[i]
			pathLowerBounds[TLowerBound[i,j].idx] = travelTimes[i,j]
			pathLowerBounds[TUpperBound[i,j].idx] = -Inf
			pathUpperBounds[TLowerBound[i,j].idx] = Inf
			pathUpperBounds[TUpperBound[i,j].idx] = travelTimes[i,j]
		end
		for i = nodes, j = out[i]
			pathLowerBounds[objConstrLower[i,j].idx] = -Inf
			pathUpperBounds[objConstrLower[i,j].idx] = 0.0
			pathLowerBounds[objConstrUpper[i,j].idx] = -Inf
			pathUpperBounds[objConstrUpper[i,j].idx] = 0.0
		end

		# Fill in variable bounds for first LP
		for i = nodes, j = pairs[i]
			variableLowerBounds[T[i,j].col] = 0.0
			variableUpperBounds[T[i,j].col] = Inf
			variableLowerBounds[epsilon[i,j].col] = 0.0
			variableUpperBounds[epsilon[i,j].col] = Inf
		end
		for i = nodes, j = out[i]
			variableLowerBounds[t[i,j].col] = roadTimes[i,j]
			variableUpperBounds[t[i,j].col] = Inf
			variableLowerBounds[delta2[i,j].col] = 0.0
			variableUpperBounds[delta2[i,j].col] = Inf
		end

		# Add path constraints
		println("**** Adding constraints ****")
		tic()
		paths, numExpensiveTurns = reconstructMultiplePathsWithExpensiveTurnsParallel(new_sp.previous, srcs, dsts, old_nodes, new_sp.real_destinations, newTimes, new_edge_dists)

		for i=1:numDataPoints
			# Update old paths
			if previousPaths[i] != 0
				for hashedPath in hashedPaths[i]
					index = hashedPathIndices[i][hashedPath]
					if split(model_type, "_")[end] == "relaxed"
						pathLowerBounds[path[i,index].idx] = - turnCost * totalNumExpensiveTurns[i,index] - delta * travelTimes[srcs[i],dsts[i]]
					elseif split(model_type, "_")[end] == "strict"
						pathLowerBounds[path[i,index].idx] = - turnCost * totalNumExpensiveTurns[i,index]
					end
					pathUpperBounds[path[i,index].idx] = Inf
				end
			end
			# If path is already in model
			if hash(paths[i]) in hashedPaths[i]
				# Choose this path to be the good one
				index = hashedPathIndices[i][hash(paths[i])]
				pathLowerBounds[path[i,index].idx] = - turnCost * numExpensiveTurns[i]
				pathUpperBounds[path[i,index].idx] = - turnCost * numExpensiveTurns[i]
			# If this is the first path for this pair of nodes
			elseif length(hashedPaths[i]) == 0
				hashedPaths[i] = Set([hash(paths[i])])
				hashedPathIndices[i] = Dict{Uint64, Int64}(hash(paths[i]) => 1)
				path[i,1] = @addConstraint(m, sum{t[paths[i][a],paths[i][a+1]], a=1:(length(paths[i])-1)} - T[srcs[i],dsts[i]] == - turnCost * numExpensiveTurns[i])
				push!(pathLowerBounds, - turnCost * numExpensiveTurns[i])
				push!(pathUpperBounds, - turnCost * numExpensiveTurns[i])
				totalNumExpensiveTurns[i,1] = numExpensiveTurns[i]
				totalPathConstraints += 1
			# If this is a new path but not the first one
			else
				push!(hashedPaths[i], hash(paths[i]))
				len = length(hashedPaths[i])
				hashedPathIndices[i][hash(paths[i])] = len
				path[i,len] = @addConstraint(m, sum{t[paths[i][a],paths[i][a+1]], a=1:(length(paths[i])-1)} - T[srcs[i],dsts[i]] == - turnCost * numExpensiveTurns[i])
				push!(pathLowerBounds, - turnCost * numExpensiveTurns[i])
				push!(pathUpperBounds, - turnCost * numExpensiveTurns[i])
				totalNumExpensiveTurns[i,len] = numExpensiveTurns[i]
				totalPathConstraints += 1
			end
		end
		toc()

		# Solve LP
		println("**** Solving LP ****")
		buildInternalModel(m)
		im = getInternalModel(m)

		MathProgBase.setvarLB!(im, variableLowerBounds)
		MathProgBase.setvarUB!(im, variableUpperBounds)
		MathProgBase.setconstrLB!(im, pathLowerBounds)
		MathProgBase.setconstrUB!(im, pathUpperBounds)
		MathProgBase.updatemodel!(im)

		MathProgBase.optimize!(im)
		status = MathProgBase.status(im)
		objective = MathProgBase.getobjval(im)
		if status == :Infeasible
			print_iis_gurobi(m, im)
			break
		end

		println("**** Setting up second LP ****")

		# Get delta values
		result = MathProgBase.getsolution(im)
		for i = nodes, j = pairs[i]
			variableLowerBounds[T[i,j].col] = result[T[i,j].col]
			variableUpperBounds[T[i,j].col] = result[T[i,j].col]
		end
		if split(model_type, "_")[end] == "relaxed"
			# Find delta values from how much constraints were actually violated
			deltaValues = zeros(length(nodes), length(nodes))
			deltaPercentage = zeros(length(nodes), length(nodes))
			# Get values of vector Ax
			constraintValues = MathProgBase.getconstrsolution(im)
			for i = 1:numDataPoints
				difference = 0.0
				# Look at all paths (except equality)
				for hashedPath in hashedPaths[i]
					if hashedPath != hash(paths[i])
						# Find value of violation. deltaValue is maximum violation
						index = hashedPathIndices[i][hashedPath]
						newDiff = - constraintValues[path[i,index].idx] - totalNumExpensiveTurns[i,index] * turnCost
						if newDiff > difference
							difference = newDiff
						end
					end
				end
				deltaValues[srcs[i], dsts[i]] = difference
				deltaPercentage[srcs[i], dsts[i]] = difference/travelTimes[srcs[i],dsts[i]]
			end
			println("Max delta value: ", maximum(deltaPercentage))
		end

		setSolver(m, GurobiSolver(TimeLimit=10000, Method=1, OutputFlag=0))
		@setObjective(m, Min, sum{delta2[i,j], i=nodes, j=out[i]})

		println("**** Adding constraints ****")
		# Set up second LP, add constraints
		@time for i=1:numDataPoints
			if previousPaths[i] != 0 && split(model_type, "_")[end] == "relaxed"
				for hashedPath in hashedPaths[i]
					index = hashedPathIndices[i][hashedPath]
					if hashedPath != hash(paths[i])
						if split(model_type, "_")[end] == "relaxed"
							value = - turnCost * totalNumExpensiveTurns[i,index] - deltaValues[srcs[i],dsts[i]]
						else #if split(model_type, "_")[end] == "strict"
							value = - turnCost * totalNumExpensiveTurns[i,index]
						end
						pathLowerBounds[path[i,index].idx] = value
					else
						previousPaths[i] = hashedPathIndices[i][hash(paths[i])]
					end
				end
			else
				previousPaths[i] = hashedPathIndices[i][hash(paths[i])]
			end
		end

		# Solve second LP
		println("**** Solving second LP ****")
		buildInternalModel(m)
		im = getInternalModel(m)

		MathProgBase.setvarLB!(im, variableLowerBounds)
		MathProgBase.setvarUB!(im, variableUpperBounds)
		MathProgBase.setconstrLB!(im, pathLowerBounds)
		MathProgBase.setconstrUB!(im, pathUpperBounds)
		MathProgBase.updatemodel!(im)

		MathProgBase.writeproblem(im, "one.lp")

		MathProgBase.optimize!(im)
		status = MathProgBase.status(im)

		# Debug if infeasible
		if status == :Infeasible
			println("!!!!!!!!!!!!!!!!!!!!!!!")
			println("!!!! Computing IIS !!!!")
			println("!!!!!!!!!!!!!!!!!!!!!!!")
			print_iis_gurobi(m, im)
			break
		# Prepare output
		elseif status == :Optimal || status == :Suboptimal
			result2 = MathProgBase.getsolution(im)
			newTimes = spzeros(length(nodes), length(nodes))
			for i = nodes, j=out[i]
				newTimes[i, j] = result2[t[i,j].col]
			end
			# Save updated Manhattan road times to file
			saveRoadTimes(newTimes, "$TESTDIR/metropolis-times-$l")
			if abs(objective - old_objective)/old_objective < 1e-10
				save("Outputs/$TESTDIR/end.jld", "num_iter", l)
				l = max_rounds
			else
				old_objective = objective
			end
		elseif status == :UserLimit
			println("!!!! User time limit exceeded !!!!")
			break
		end
		if l <= max_rounds
			println("**** Computing shortest paths ****")
			@time new_graph, new_edge_dists, new_nodes = modifyGraphForDijkstra(graph, newTimes, positions, turn_cost=turnCost)
			@time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
			println("**** Error calculation ****")
			printErrorStatsToFile("Outputs/$TESTDIR/errorstats-$(actual_l).txt", real_times, travel_times, new_sp.traveltime, numRides, meanRoadTime, newTimes, num_nodes = nv(graph))
		end
		l += 1
	end
	return status, newTimes
end

function compute_error(realData::AbstractArray{Float64, 2}, inputData::AbstractArray{Float64, 2}, algorithmData::AbstractArray{Float64}, numRides::AbstractArray{Int, 2}, realLinkTimes::AbstractArray{Float64,2}, newLinkTimes::AbstractArray{Float64,2}; num_nodes::Int = 192, out=Array{Int}[])
	"""
	Given a matrix of TABs from real Data, input Data and our algorithm, returns some statistical measures of error
	"""
	num_links = 0
	num_results = 0
	num_nonzero_results = 0
	average_squared_error = 0
	average_relative_error = 0
	average_bias = 0
	average_nonzero_input_error = 0
	average_nonzero_result_error = 0
	average_link_relative_error = 0
	average_link_bias = 0
	for i = 1:num_nodes, j=1:num_nodes
		if realData[i,j] > 0 && (length(out) == 0 || j in out[i])
			average_squared_error = (num_results * average_squared_error + (realData[i,j] - algorithmData[i,j]) ^ 2) / (num_results + 1)
			average_relative_error = (num_results * average_relative_error + abs(realData[i,j] - algorithmData[i,j])/realData[i,j])/(num_results + 1)
			average_bias = (num_results * average_bias + (algorithmData[i,j] - realData[i,j]))/(num_results + 1)
			num_results += 1
			if numRides[i,j] > 0
				average_nonzero_input_error = (num_nonzero_results * average_nonzero_input_error + abs(realData[i,j] - inputData[i,j])/realData[i,j])/(num_nonzero_results + 1)
				average_nonzero_result_error = (num_nonzero_results * average_nonzero_result_error + abs(realData[i,j] - algorithmData[i,j])/realData[i,j])/(num_nonzero_results + 1)
				num_nonzero_results += 1
			end
		end
		if realLinkTimes[i,j] > 0
			average_link_relative_error = (num_links * average_link_relative_error + abs(newLinkTimes[i,j]-realLinkTimes[i,j])/realLinkTimes[i,j])/(num_links + 1)
			average_link_bias = (num_links * average_link_bias + newLinkTimes[i,j] - realLinkTimes[i,j])/(num_links + 1)
			num_links += 1
		end
	end
	println("------------------------------------------------")
	println("Average squared error: \t\t\t\t\t", average_squared_error)
	println("Average relative error: \t\t\t\t", average_relative_error)
	println("Average bias: \t\t\t\t\t\t", average_bias)
	println("Average original relative error on nonzero rides: \t", average_nonzero_input_error)
	println("Average result relative error on nonzero rides: \t", average_nonzero_result_error)
	println("Average relative error on link times: \t\t\t", average_link_relative_error)
	println("Average bias on link times:\t\t\t\t", average_link_bias)
	println("------------------------------------------------")
	return average_squared_error, average_relative_error, average_bias, average_nonzero_input_error, average_nonzero_result_error, average_link_relative_error, average_link_bias
end

function printErrorStatsToFile(fileName::AbstractString, realData::AbstractArray{Float64, 2}, inputData::AbstractArray{Float64, 2}, algorithmData::AbstractArray{Float64}, numRides::AbstractArray{Int, 2}, realLinkTimes::AbstractArray{Float64,2}, newLinkTimes::AbstractArray{Float64,2}; num_nodes::Int = 192, out=Array{Int}[])
	average_squared_error, average_relative_error, average_bias, average_nonzero_input_error, average_nonzero_result_error, average_link_relative_error, average_link_bias = compute_error(realData, inputData, algorithmData, numRides, realLinkTimes, newLinkTimes, num_nodes=num_nodes, out=out)
	"""
	Same as compute_error, except prints to file as well
	"""
	file = open(fileName, "w")
	write(file, string("Average squared error: \t\t\t\t\t", average_squared_error, "\n"))
	write(file, string("Average relative error: \t\t\t\t", average_relative_error, "\n"))
	write(file, string("Average bias: \t\t\t\t\t\t", average_bias, "\n"))
	write(file, string("Average original relative error on nonzero rides: \t", average_nonzero_input_error, "\n"))
	write(file, string("Average result relative error on nonzero rides: \t", average_nonzero_result_error, "\n"))
	write(file, string("Average relative error on link times: \t\t\t", average_link_relative_error, "\n"))
	write(file, string("Average bias on link times:\t\t\t\t", average_link_bias, "\n"))
	close(file)
end

# # This code loads the data and calls the LP
# adjacencyList = load("Inputs/input-graph.jld", "adjList")
# graph = DiGraph(length(adjacencyList))
# for i = 1:nv(graph), j in adjacencyList[i]
# 	add_edge!(graph, i, j)
# end
# real_times = load("Inputs/input-realTimes.jld", "realTimes")
# minRoadTime = load("Inputs/input-speedLimits.jld", "speedLimits")
# meanRoadTime = load("Inputs/input-meanTimes.jld", "meanTimes")
# distances = load("Inputs/input-distances.jld", "distances")
# coordinates = load("Inputs/input-positions.jld", "coordinates")
# positions = Coordinates[]
# for i = 1:nv(graph)
# 	push!(positions, Coordinates(coordinates[i,1], coordinates[i,2]))
# end

# travel_times = load("Inputs/input-travelTimes-$(PROB).jld", "travelTimes")
# num_rides = load("Inputs/input-numRides-$(PROB).jld", "numRides")

function compute_cost_function(realData::AbstractArray{Float64,2}, inputData::AbstractArray{Float64,2}, numRides::AbstractArray{Int,2}; num_nodes = 192)
	cost = 0
	for i = 1:num_nodes, j=1:num_nodes
		if numRides[i,j] > 0
			cost += abs(realData[i,j] - inputData[i,j]) * sqrt(numRides[i,j]/inputData[i,j])
		end
	end
	return cost
end

# @time status, new_times = new_LP(graph, travel_times, num_rides, minRoadTime, distances, positions)
# @time new_graph, new_edge_dists, new_nodes = modifyGraphForDijkstra(graph, new_times, positions, turn_cost=TURN_COST)
# @time new_sp = parallelShortestPathsWithTurnsAuto(graph, new_graph, new_edge_dists, new_nodes)
# print("")
# compute_cost_function(real_times, travel_times, num_rides, num_nodes = nv(graph))

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
const YaoCompilerJob = CompilerJob{<:YaoCompileTarget}

function ci_cache_populate(
    cache::CodeCache,
    target::YaoCompileTarget,
    mi::MethodInstance,
    min_world,
    max_world,
    options::HardwareFreeOptions,
)
    interp = YaoInterpreter(; target, options, cache)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)
    wvc = WorldView(cache, min_world, max_world)
    @assert Core.Compiler.haskey(wvc, mi)

    # if src is rettyp_const, the codeinfo won't cache ci.inferred
    # (because it is normally not supposed to be used ever again).
    # to avoid the need to re-infer, set that field here.
    ci = Core.Compiler.getindex(wvc, mi)
    if ci !== nothing && ci.inferred === nothing
        ci.inferred = src
    end
    return ci
end

function compile_method_instance(
    cache,
    target::YaoCompileTarget,
    method_instance::MethodInstance,
    options::HardwareFreeOptions,
    world,
)
    if GPUCompiler.ci_cache_lookup(cache, method_instance, world, typemax(Cint)) === nothing
        ci_cache_populate(cache, target, method_instance, world, typemax(Cint), options)
    end

    # set-up the compiler interface
    jl_options = Base.JLOptions()
    debug_info_kind = if jl_options.debug_level == 0
        LLVM.API.LLVMDebugEmissionKindNoDebug
    elseif jl_options.debug_level == 1
        LLVM.API.LLVMDebugEmissionKindLineTablesOnly
    elseif jl_options.debug_level >= 2
        LLVM.API.LLVMDebugEmissionKindFullDebug
    end

    lookup_fun =
        (mi, min_world, max_world) -> GPUCompiler.ci_cache_lookup(cache, mi, min_world, max_world)
    lookup_cb = @cfunction($lookup_fun, Any, (Any, UInt, UInt))
    params = Base.CodegenParams(;
        track_allocations = false,
        code_coverage = false,
        prefer_specsig = true,
        gnu_pubnames = false,
        debug_info_kind = Cint(debug_info_kind),
        lookup = Base.unsafe_convert(Ptr{Nothing}, lookup_cb),
    )

    # generate IR
    GC.@preserve lookup_cb begin
        native_code = ccall(
            :jl_create_native,
            Ptr{Cvoid},
            (Vector{MethodInstance}, Base.CodegenParams, Cint),
            [method_instance],
            params,
            1,
        ) #=extern policy=#
        @assert native_code != C_NULL
        llvm_mod_ref = ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef, (Ptr{Cvoid},), native_code)
        @assert llvm_mod_ref != C_NULL
        llvm_mod = LLVM.Module(llvm_mod_ref)
    end

    # get the top-level code
    code = GPUCompiler.ci_cache_lookup(cache, method_instance, world, typemax(Cint))

    # get the top-level function index
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(
        :jl_get_function_id,
        Nothing,
        (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
        native_code,
        code,
        llvm_func_idx,
        llvm_specfunc_idx,
    )
    @assert llvm_func_idx[] != -1
    @assert llvm_specfunc_idx[] != -1

    # get the top-level function)
    llvm_func_ref = ccall(
        :jl_get_llvm_function,
        LLVM.API.LLVMValueRef,
        (Ptr{Cvoid}, UInt32),
        native_code,
        llvm_func_idx[] - 1,
    )
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(
        :jl_get_llvm_function,
        LLVM.API.LLVMValueRef,
        (Ptr{Cvoid}, UInt32),
        native_code,
        llvm_specfunc_idx[] - 1,
    )
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)

    # # configure the module
    # triple!(llvm_mod, llvm_triple(job.target))
    # if julia_datalayout(job.target) !== nothing
    #     datalayout!(llvm_mod, julia_datalayout(job.target))
    # end
    return llvm_specfunc, llvm_func, llvm_mod
end

function jit_compile(job::YaoCompilerJob)
    world = Base.get_world_counter()
    mi = method_instance(job.source.f, job.source.tt, world)
    cache = get!(GLOBAL_CI_CACHE, job.target, GPUCompiler.CodeCache())
    llvm_specfunc, llvm_func, llvm_mod =
        compile_method_instance(cache, job.target, mi, job.params, world)

    specfunc_name = LLVM.name(llvm_specfunc)
    func_name = LLVM.name(llvm_func)
    linkage!(llvm_func, LLVM.API.LLVMExternalLinkage)
    linkage!(llvm_specfunc, LLVM.API.LLVMExternalLinkage)
    run_pipeline!(llvm_mod)
    return llvm_mod, func_name, specfunc_name
end

struct TargetHostKernel{F,TT,Target<:YaoCompileTarget}
    f::F
    target::Target
    specfunc::Ptr{Cvoid}
    func::Ptr{Cvoid}
end

function jit_link(job::CompilerJob, (llvm_mod, func_name, specfunc_name))
    fspec = job.source
    jitted_mod = compile!(orc[], llvm_mod)
    specfunc_addr = addressin(orc[], jitted_mod, specfunc_name)
    specfunc_ptr = pointer(specfunc_addr)
    func_addr = addressin(orc[], jitted_mod, func_name)
    func_ptr = pointer(func_addr)
    if specfunc_ptr === C_NULL || func_ptr === C_NULL
        @error "Compilation error" fspec specfunc_ptr func_ptr
    end
    TargetHostKernel{typeof(fspec.f),fspec.tt,typeof(job.target)}(
        fspec.f,
        job.target,
        specfunc_ptr,
        func_ptr,
    )
end

const jit_compiled_cache = Dict{UInt,Any}()

compile(target::YaoCompileTarget, f, tt) = compile(target, f, tt, HardwareFreeOptions())

# NOTE: by default we just compile target to generic Julia code as fallback
# this makes it compatible with the old simulators, symbolic engines etc.
function compile(
    target::YaoCompileTarget,
    f::F,
    tt::TT = Tuple{},
    options::HardwareFreeOptions = HardwareFreeOptions(),
) where {F,TT<:Type}
    fspec = FunctionSpec(f, tt, false, nothing)
    job = CompilerJob(target, fspec, options)
    return GPUCompiler.cached_compilation(jit_compiled_cache, job, jit_compile, jit_link)
end

@generated function (kernel::TargetHostKernel{F,TT})(args...) where {F,TT}
    expr = quote
        args = Any[args...]
        ccall(kernel.func, Any, (Any, Ptr{Any}, Int32), kernel.f, args, length(args))
    end
    return expr
end


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

mutable struct DedicatedHostInstanceViewWithName <: SwaggerModel
    assetId::Any # spec type: Union{ Nothing, String } # spec name: assetId
    availableCapacity::Any # spec type: Union{ Nothing, DedicatedHostAvailableCapacity } # spec name: availableCapacity
    statuses::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: statuses

    function DedicatedHostInstanceViewWithName(;assetId=nothing, availableCapacity=nothing, statuses=nothing)
        o = new()
        validate_property(DedicatedHostInstanceViewWithName, Symbol("assetId"), assetId)
        setfield!(o, Symbol("assetId"), assetId)
        validate_property(DedicatedHostInstanceViewWithName, Symbol("availableCapacity"), availableCapacity)
        setfield!(o, Symbol("availableCapacity"), availableCapacity)
        validate_property(DedicatedHostInstanceViewWithName, Symbol("statuses"), statuses)
        setfield!(o, Symbol("statuses"), statuses)
        o
    end
end # type DedicatedHostInstanceViewWithName

const _property_map_DedicatedHostInstanceViewWithName = Dict{Symbol,Symbol}(Symbol("assetId")=>Symbol("assetId"), Symbol("availableCapacity")=>Symbol("availableCapacity"), Symbol("statuses")=>Symbol("statuses"))
const _property_types_DedicatedHostInstanceViewWithName = Dict{Symbol,String}(Symbol("assetId")=>"String", Symbol("availableCapacity")=>"DedicatedHostAvailableCapacity", Symbol("statuses")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ DedicatedHostInstanceViewWithName }) = collect(keys(_property_map_DedicatedHostInstanceViewWithName))
Swagger.property_type(::Type{ DedicatedHostInstanceViewWithName }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DedicatedHostInstanceViewWithName[name]))}
Swagger.field_name(::Type{ DedicatedHostInstanceViewWithName }, property_name::Symbol) =  _property_map_DedicatedHostInstanceViewWithName[property_name]

function check_required(o::DedicatedHostInstanceViewWithName)
    true
end

function validate_property(::Type{ DedicatedHostInstanceViewWithName }, name::Symbol, val)
end


text:
"""
LazyJSON.jl

See README.md
"""

module LazyJSON

const JSON = LazyJSON

abstract type Value end

const Key = Union{Int, Base.CodeUnits{UInt8,Base.String}}

parse(s::AbstractString; path=[]) = Base.parse(JSON.Value, s;
                                               path=mkpath(path))

mkpath(p) = Key[x isa AbstractString ? codeunits(x) : x for x in p]

const Bytes = Base.CodeUnits{UInt8,Base.String}

const enable_assertions = false

const disable_coroutine = false

const wait_after_each_value = false

# Parser Data Structures

"""
The JSON parser encountered invalid input.
"""

struct ParseError <: Exception
    bytes::Base.String
    index::Int
    message::Base.String
end

# JSON Value Data Structures

"""
`JSON.String` stores the start and end indexes of a string within a JSON text.
"""

struct String # FIXME <: AbstractString ?
    bytes::Bytes
    first::Int
    last::Int
    has_escape::Bool
end

"""
`JSON.Number` stores the start and end indexes of a number within a JSON text.
"""

struct Number <: Base.Number
    bytes::Bytes
    first::Int
    last::Int
end

"""
Fields:
 - `task`, runs `parse_value` function as a co-routine.
 - `result`, the top-level JSON value.
 - `error`, exception thrown by `parse_value`.
"""

mutable struct Parser
    task::Task
    result
    error
    wait_count::Int
    Parser() = (p = new(); p.error = nothing; p.wait_count = 0; p)
end

"""
`getjson(v, i)` is used by the `AbstractArray` and `AbstractDict` API methods
to retrieve values from `JSON.Object` and `JSON.Array`.
"""

function getjson end

"""
If `use_promotejson` is `false`: The `getjson(v, i)` function is simply `v[i]`.
In this case numbers and strings are returned as `JSON.Number` and `JSON.String`.

If `use_promotejson` is `true` the `getjson(v, i)` function calls `promotejson`
on the value it retrieves from a collection and replaces the original value
with the promoted value.  In this case numbers and strings are returned as
`Base.Int64`, `Base.Float64` and `Base.SubString`.
"""

const use_promotejson = false

@static if use_promotejson

    function getjson(v, i)
        x = v[i]
        j = promotejson(x)
        if !(x === j)
            v[i] = j
        end
        return j
    end

else

    getjson(v, i) = @inbounds(v[i])

end

promotejson(v) = v
promotejson(v::JSON.String) = convert(SubString, v)
promotejson(v::JSON.Number) = convert(Base.Number, v)

"""
`JSON.Array` implements the `AbstractArray` interface.

The `AbstractArray` methods call `parse_more!` as needed to ensure that the
requested items are available in the array.

The parser task calls `push!` to add parsed values to the array.
`push!` calls `wait` to return control to the main task.
"""

mutable struct Array <: AbstractArray{Any, 1}
    parser::Parser
    v::Vector{Any}
    iscomplete::Bool
    Array(p) = new(p, [], false)
end

Base.push!(a::JSON.Array, x) = (push!(a.v, x); wait_for_parse_more(a.parser))

"""
`JSON.Object` implements the `AbstractDict` interface.

The keys an values are stored in a single `Vector`:
`[key1, value2, key2, value2, ...]`

The `AbstractDict` methods call `parse_more!` as needed to ensure that the
requested items are available in the array.

The parser task calls `push!` to add parsed values to the array.
`push!` calls `wait` to return control to the main task.
"""

mutable struct Object <: AbstractDict{AbstractString, Any}
    parser::Parser
    v::Vector{Any}
    iscomplete::Bool
    Object(p) = new(p, [], false)
end

function Base.push!(o::Object, x)
    push!(o.v, x)
    wait_for_parse_more(o.parser)
end

# JSON Parser

"""
Parse JSON text.

Returns a `Bool`, `Number`, `AbstractString`, `AbstractArray`, `AbstractDict` or
`Nothing`.

The parser runs in a lazy co-routine. If the returned value is a collection
type it will be incomplete. Each call to `parse_more!(::Parser)` causes another
an additional value to be parsed and added to the result collection.
`parse_more!` is called automatically by the `AbstractArray` and `AbstractDict`
interface methods.
"""

function Base.parse(::Type{JSON.Value},
                    as::AbstractString;
                    path::Vector{Key}=[])
    s = codeunits(as)                   # Convert string to byte-vector
    l = sizeof(s)
    p = Parser()                        # Run parse_value() as a co-routine.
    v = isempty(path) ? p : path
    @static if disable_coroutine
    parse_value(p, s, l, v, 1)
    else
    p.task = @task try                  # Pass the parser as the parent value
        i = parse_value(p, s, l, v, 1)  # and start at index 1.
        if enable_assertions
            check_end(p, s, l, i)
        end
    catch e                             # Capturue parse errors to be rethrown
        p.error = e                     # in the main task by parse_more()
    end
    end
    parse_more!(p)                      # Ensure that the top level JSON value
    return p.result                     # is push!-ed into p.result.
end

iscomplete(x) = true
iscomplete(x::JSON.Object) = x.iscomplete
iscomplete(x::JSON.Array) = x.iscomplete

function Base.push!(p::Parser, value)
    p.result = value                    # Store the result value and wait for
    disable_coroutine || wait()         # the next call to parse_more!().
end

@static if disable_coroutine
parse_more!(p::Parser) = nothing
else
function parse_more!(p::Parser)
    if p.error != nothing               # Don't try to parse more if the last
        throw(p.error)                  # call resulted in an error.
    end

    schedule(current_task())            # Ensure that the main task will restart
    yieldto(p.task)                     # then run the parse_value() co-routine.

    if p.error != nothing               # Rethrow error captured from the
        throw(p.error)                  # parse_value() co-routine.
    end
end
end

function parse_more!(collection, index)
    while !collection.iscomplete &&     # Keep parsing until a specified index
           index > length(collection.v) # is available, or the end is reached.
        parse_more!(collection.parser)
    end
end

parse_all!(c) = while !c.iscomplete parse_more!(c.parser) end

function wait_for_parse_more(p)
    @static if disable_coroutine
        return
    end
    @static if wait_after_each_value
        wait()
    else
        p.wait_count += 1
        if p.wait_count > 100           # Keep parsing for a while before
            p.wait_count = 0            # yeilding control back to the main
            wait()                      # task to avoid task switching overhead.
        end
    end
end

function check_end(p::Parser, s::Bytes, l, i)
    if !iscomplete(p.result)
        throw(JSON.ParseError(s.s, i, "unexpected end of input"))
    end
    if i <= l
        i, x = skip_ws(s, l, i)
        if i <= l && !isws(x)
            throw(JSON.ParseError(s.s, i, "'$(Char(x))' after end of input"))
        end
    end
end

"""
Throw a JSON.ParseError if the parse function's current character is not `x`.
"""

macro assume(x)
    if !enable_assertions
        return :()
    end
    s = typeof(x) == Char ? "'$x'" : string(x)
    esc(quote
        if i > sizeof(s) || !(s[i] in UInt8[($x)...])
            throw(JSON.ParseError(s.s, i, "expected $($s)"))
        end
    end)
end

"""
Whitespace

https://tools.ietf.org/html/rfc7159#section-2

      ws = *(
              %x20 /              ; Space
              %x09 /              ; Horizontal tab
              %x0A /              ; Line feed or New line
              %x0D )              ; Carriage return
"""

isws(x) = x == 0x20 ||
          x == 0x09 ||
          x == 0x0A ||
          x == 0x0D

get_x(s, l, i) = i > l ? 0x00 : @inbounds(s[i])

next_x(s, l, i) = (i += 1; (i, get_x(s, l, i)))

function skip_ws(s, l, i)
    x = get_x(s, l, i)
    while isws(x)
        i, x = next_x(s, l, i)
    end
    return i, x
end

"""
Values

https://tools.ietf.org/html/rfc7159#section-3

      value = false / null / true / object / array / number / string

      false = %x66.61.6c.73.65   ; false

      null  = %x6e.75.6c.6c      ; null

      true  = %x74.72.75.65      ; true
"""

parse_false(p, s, l, v, i) = (v != nothing && push!(v, false);   i + 5)
parse_null(p, s, l, v, i)  = (v != nothing && push!(v, nothing); i + 4)
parse_true(p, s, l, v, i)  = (v != nothing && push!(v, true);    i + 4)

"""
Objects

https://tools.ietf.org/html/rfc7159#section-4

      object = begin-object [ member *( value-separator member ) ]
               end-object

      member = string name-separator value
"""

function parse_object(p, s, l, v, i)

    if v isa AbstractVector{Key}
        key = v[1]
        if key isa Base.CodeUnits
            vl = length(v)
            path = vl > 1 ? view(v, 2:vl) : p
        else
            key = nothing
        end
    else
        key = nothing
    end
    if v == nothing || v isa AbstractVector{Key}
        o = nothing
    else
        o = Object(p)
        push!(v, o)
    end

    @assume '{'
    i += 1
    i, x = skip_ws(s, l, i)
    while x != UInt8('}')

        i, x = skip_ws(s, l, i)
        k = i
        i = parse_string(p, s, l, o, i)
        oo = key != nothing && view(s, k+1:i-2) == key ? path : o
        i, x = skip_ws(s, l, i)
        @assume ':'
        i += 1

        i = parse_value(p, s, l, oo, i)
        x = get_x(s, l, i)
        if x == UInt8(',')
            i += 1
        else
            @assume '}'
        end
    end

    if o isa Object
        o.iscomplete = true
    end

    return i + 1
end

"""
Arrays

https://tools.ietf.org/html/rfc7159#section-5

   array = begin-array [ value *( value-separator value ) ] end-array
"""

function parse_array(p, s, l, v, i)

    if v isa AbstractVector{Key}
        key = v[1]
        if key isa Int
            vl = length(v)
            path = vl > 1 ? view(v, 2:vl) : p
        else
            key = 0
        end
    else
        key = 0
    end
    if v == nothing || v isa AbstractVector{Key}
        a = nothing
    else
        a = JSON.Array(p)
        push!(v, a)
    end

    @assume '['
    i += 1
    i, x = skip_ws(s, l, i)
    count = 1
    while x != UInt8(']')
        aa = key != 0 && count == key ? path : a
        i = parse_value(p, s, l, aa, i)
        x = get_x(s, l, i)
        if x == UInt8(',')
            i += 1
        else
            @assume ']'
        end
        count += 1
    end

    if a isa JSON.Array
        a.iscomplete = true
    end

    return i + 1
end

"""
Numbers

https://tools.ietf.org/html/rfc7159#section-6

      number = [ minus ] int [ frac ] [ exp ]

      decimal-point = %x2E       ; .

      digit1-9 = %x31-39         ; 1-9

      e = %x65 / %x45            ; e E

      exp = e [ minus / plus ] 1*DIGIT

      frac = decimal-point 1*DIGIT

      int = zero / ( digit1-9 *DIGIT )

      minus = %x2D               ; -

      plus = %x2B                ; +

      zero = %x30                ; 0
"""

function parse_number(p, s, l, v, i)

    start = i

    x = @inbounds(s[i])
    if x == UInt8('-')
        i, x = next_x(s, l, i)
    end
    @assume '0':'9'
    prev = x
    i, x = next_x(s, l, i)
    if prev == UInt8('0')
        @assume '.', 'e', 'E', ']', '}', ','
    end

    while !isws(x) &&
          x != UInt8(']') &&
          x != UInt8('}') &&
          x != UInt8(',') &&
          x != 0x00
        prev = x
        i, x = next_x(s, l, i)
        if prev == UInt8('.')
            @assume '0':'9'
        end
    end

    i -= 1
    @assume '0':'9'
    if v != nothing
        push!(v, JSON.Number(s, start, i))
    end

    return i + 1
end

"""
Strings

https://tools.ietf.org/html/rfc7159#section-7

      string = quotation-mark *char quotation-mark

      char = unescaped /
          escape (
              %x22 /          ; "    quotation mark  U+0022
              %x5C /          ; \\    reverse solidus U+005C
              %x2F /          ; /    solidus         U+002F
              %x62 /          ; b    backspace       U+0008
              %x66 /          ; f    form feed       U+000C
              %x6E /          ; n    line feed       U+000A
              %x72 /          ; r    carriage return U+000D
              %x74 /          ; t    tab             U+0009
              %x75 4HEXDIG )  ; uXXXX                U+XXXX

      escape = %x5C              ; \\

      quotation-mark = %x22      ; "

      unescaped = %x20-21 / %x23-5B / %x5D-10FFFF
"""

function parse_string(p, s, l, v, i)

    start = i
    has_escape = false

    @assume '"'
    i, x = next_x(s, l, i)

    while x != UInt8('"') && x != 0x00
        if x == UInt8('\\')
            has_escape = true
            i, x = next_x(s, l, i)
            @assume '"', '\\', '/', 'b', 'f', 'n', 'r', 't', 'u'
        end
        i, x = next_x(s, l, i)
    end
    @assume '"'

    if v != nothing
        push!(v, JSON.String(s, start+1, i-1, has_escape))
    end

    return i + 1
end

"""
    parse_value(p::Parser, string::Bytes, parent_value, index)

Parse a [JSON value](https://tools.ietf.org/html/rfc7159#section-3) from a
`string` of bytes starting at `index`. The parsed value is appended to the
parent value (e.g. `push!(parent_value, parsed_value)`).
"""

function parse_value(p::Parser, s::Bytes, l, v, i)

    i, x = skip_ws(s, l, i)

    i = if x == UInt8('{')                      parse_object(p, s, l, v, i)
    elseif x == UInt8('[')                       parse_array(p, s, l, v, i)
    elseif x == UInt8('"')                      parse_string(p, s, l, v, i)
    elseif x >= UInt8('0')  &&
           x <= UInt8('9')  ||
           x == UInt8('-')                      parse_number(p, s, l, v, i)
    elseif x == UInt8('f')                       parse_false(p, s, l, v, i)
    elseif x == UInt8('n')                        parse_null(p, s, l, v, i)
    elseif x == UInt8('t')                        parse_true(p, s, l, v, i)
    else
        throw(JSON.ParseError(s.s, i, "invalid input"))
    end

    i, x = skip_ws(s, l, i)
    return i
end

# Message Display

Base.show(io::IO, s::JSON.String) =
    print(io, SubString(s.bytes.s, s.first-1, s.last+1))

Base.show(io::IO, n::JSON.Number) = print(io, string(n))

function Base.show(io::IO, e::JSON.ParseError)

    s = e.bytes
    l = findprev(isequal('\n'), s, e.index)
    l = l != nothing ? l + 1 : 1
    r = findnext(isequal('\n'), s, max(1, e.index-1))
    r = r != nothing ? r - 1 : length(s)
    l = min(length(s), l)
    line_number = length(split(SubString(s, 1, l), '\n'))
    col_number = e.index - l + 1
    print(io, "JSON.ParseError: ", e.message,
              " at line ", line_number, ", col ", col_number, "\n",
              SubString(s, l, r), "\n",
              lpad("", col_number - 1, " "), "^")
end

include("AbstractString.jl")
include("Number.jl")
include("AbstractArray.jl")
include("AbstractDict.jl")
include("LazyerJSON.jl")

using .LazyerJSON

end # module JSON