text:
import GeometryTypes: ⊖, support_vector_max
import GeometryTypes: gjk0
import GeometryTypes: type_immutable, with_immutable, make_immutable

context("gjk") do

    context("gjk examples") do
        context("two simplices") do
            c1 = Simplex(Vec(-1.))
            c2 = Simplex(Vec(4.))
            @fact gjk(c1,c2) --> roughly(5)
            @fact min_euclidean(c1,c2) --> roughly(5)

            c1 = Simplex(Vec(-1.,0,0))
            c2 = Simplex(Vec(4.,0,0))
            @fact gjk(c1,c2) --> roughly(5)
            @fact min_euclidean(c1,c2) --> roughly(5)

            c1 = FlexibleConvexHull([Vec(0.,0), Vec(0.,1), Vec(1.,0),Vec(1.,1)])
            c2 = Simplex(Vec(4.,0.5))
            @fact gjk(c1, c2) --> roughly(3)
            @fact min_euclidean(c1,c2) --> roughly(3)

            pt1 = Vec(1,2,3.)
            pt2 = Vec(3,4,5.)
            @fact gjk(pt1, pt2) --> roughly(norm(pt1-pt2))
            @fact min_euclidean(pt1, pt2) --> roughly(norm(pt1-pt2))
        end

        context("gjk intersecting lines") do
            c1 = Simplex(Vec(1,1.), Vec(1, 2.))
            @fact gjk(c1, c1) --> 0.
            @fact min_euclidean(c1,c1) --> 0.

            c2 = Simplex(Vec(1,1.), Vec(10, 2.))
            @fact gjk(c1, c2) --> 0.
            @fact min_euclidean(c1,c2) --> 0.

            c3 = Simplex(Vec(0, 1.), Vec(2,2.))
            md = vertices(c1 ⊖ c3)
            @fact md --> [Vec(1.0,0.0),Vec(-1.0,-1.0),Vec(1.0,1.0),Vec(-1.0,0.0)]
            @fact gjk0(FlexibleConvexHull(md)) --> (Vec(0, 0.), 0.)
            @fact gjk(c1, c3) --> 0.
            @fact min_euclidean(c1,c3) --> 0.

        end

        context("Cube") do
            c = HyperCube(Vec(0.5,0.5,0.5), 1.)
            @fact min_euclidean(Vec(2,2,2.), c) ≈ gjk(Vec(2,2,2.), c) ≈ √(3/4) --> true

            s = Simplex(Vec(1, 0.5, 0.5), Vec(1,2,3.))
            @fact 0 <= min_euclidean(s, c) <= 1e-14 --> true
            @fact 0 <= gjk(s, c) <= 1e-14 --> true

            s = Simplex(Vec(2,2,2.), Vec(2,3,2.), Vec(2,2,7.), Vec(3,4,5.))
            @fact min_euclidean(s,c) ≈ gjk(s, c) ≈ √(3/4) --> true
        end

    end

    context("support_vector_max") do
        r = HyperRectangle(Vec(-0.5, -1.), Vec(1., 2.))
        @fact support_vector_max(r, Vec(1,0.)) --> (Vec(0.5,-1.), 0.5)
        @fact support_vector_max(r, Vec(2,0.)) --> (Vec(0.5,-1.), 1.)
        @fact support_vector_max(r, Vec(-1,0.)) --> (Vec(-0.5,-1.), 0.5)
        @fact support_vector_max(r, Vec(0, 1.)) --> (Vec(-0.5,1.), 1.)
        @fact support_vector_max(r, Vec(1, 1.)) --> (Vec(0.5,1.), 1.5)
        @fact support_vector_max(FlexibleConvexHull(r), Vec(1, 1.)) --> (Vec(0.5,1.), 1.5)

        c1 = Simplex(Vec(1,1.), Vec(1, 2.))
        c3 = Simplex(Vec(0, 1.), Vec(2,2.))
        md = c1 ⊖ c3
        fh = FlexibleConvexHull(md)
        for v in [Vec(1,0.), Vec(12,-10.), Vec(0,-1.), Vec(1,1.)]
            v_md, s_md = support_vector_max(md, v)
            v_fh, s_fh = support_vector_max(fh, v)
            @fact s_md ≈ s_fh --> true
            @fact v_md ≈ v_fh --> true
        end
    end

    context("make immutable") do
        T = FixedSizeArrays.Vec{2,Float64}
        n = 3
        S = Simplex{n, T}
        FS = FlexibleSimplex{T}
        fs = FS([Vec(0.0,0.0), Vec(0.0,1.0), Vec(1.0,0.0)])
        s = S((Vec(0.0,0.0), Vec(0.0,1.0), Vec(1.0,0.0)))
        @fact (@inferred type_immutable(FS, Val{3})) --> S
        @fact type_immutable(fs) --> S
        @fact (@inferred type_immutable(fs, Val{3})) --> S
        @fact make_immutable(fs) --> s
        @fact (@inferred make_immutable(fs, Val{3})) --> s
        @fact with_immutable(identity, fs) --> s
    end

end


text:
I = (1,1)
G = (1,2)
αs = linspace(1/6,4,100)

N = 400
h = 1e-6
σoff = -0.2

NEV = 4
L = 20

@everywhere include("exciton-hamiltonian.jl")

######################### Just a Test #####################

@everywhere function energies(α,I,G,p0,L,N,NEV,σoff)
    BLAS.set_num_threads(1)

    σ = σoff - 2*α^2 - maximum(I)^2/2 - maximum(I)*maximum(p0)

    println("σ = $σ")
    println("p0 = $p0")
    run(`hostname`)

    @time Hm = Hamiltonian(α,I,G,p0,L,N)
    @time solution = eigs(Hm,nev=NEV,tol=1e-6,maxiter=10000,which=:LM,sigma=σ)

    return real(solution[1])
end

p0 = [((0,0,0),α) for α in αs]
px = [((h,0,0),α) for α in αs]
py = [((0,h,0),α) for α in αs]

@sync begin
    @async global E0 = pmap(x->energies(x[2],I,G,x[1],L,N,NEV,σoff),p0)
    @async global Ex = pmap(x->energies(x[2],I,G,x[1],L,N,NEV,σoff),px)
    @async global Ey = pmap(x->energies(x[2],I,G,x[1],L,N,NEV,σoff),py)
end

DEx = (Ex - E0)/h
DEy = (Ey - E0)/h

using JLD
save("outdir/derivatives-$N-$NEV.jld","αs",αs,"DEx",DEx,"DEy",DEy)

ext: jl

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
using Random
using MPI, MPIAsyncPools
using CodedComputing
using HDF5, LinearAlgebra
using ArgParse
using Dates
using Distributions
using DataStructures

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0
const rank = MPI.Comm_rank(comm)
const isroot = MPI.Comm_rank(comm) == root
const data_tag = 0
const control_tag = 1
const COMMON_BYTES = 8 # number of bytes reserved for communication from worker_main to coordinator_main

"""

Parse command line arguments. One can add implementation-specific arguments or parsing by
defining the functions `update_argsettings!` and `update_parsed_args!`, respectively. See
the function body for details. Returns a dictionary containing the parsed arguments.
"""
function parse_commandline(isroot::Bool)

    # setup error handling so that only the root node prints error messages
    s = ArgParseSettings("Principal component analysis MPI kernel", autofix_names=true)
    function root_handler(settings::ArgParseSettings, err, err_code::Integer=1)
        println(Base.stderr, err.text)
        println(Base.stderr, usage_string(settings))
        MPI.Finalize()
        exit(err_code)
    end
    function worker_handler(settings::ArgParseSettings, err)
        MPI.Finalize()
        exit(0)
    end
    if isroot
        s.exc_handler = root_handler
    else
        s.exc_handler = worker_handler
        s.add_help = false
    end

    # command-line arguments common to all pca implementations
    @add_arg_table s begin
        "inputfile"
            help = "HDF5 file containing the input data set"
            required = true
            arg_type = String
            range_tester = HDF5.ishdf5
        "outputfile"
            help = "HFD5 file to write the output to"
            required = true
            arg_type = String
            range_tester = (x) -> !isfile(x) || HDF5.ishdf5(x)
        "--niterations"
            help = "Number of iterations to run the algorithm for"
            default = 10
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--inputdataset"
            help = "Input dataset name"
            default = "X"
            arg_type = String
        "--outputdataset"
            help = "Output dataset name"
            default = "V"
            arg_type = String
        "--saveiterates"
            help = "Save all intermediate iterates to the output file"
            action = :store_true
        "--nsavediterates"
            help = "Maximum number of iterates to save to disk"
            default = 20
            arg_type = Int
            range_tester = (x) -> 0 < x
        "--enablegc"
            help = "Enable garbage collection while running the computation (defaults to disabled)"
            action = :store_true
        "--loadbalance"
            help = "Enable the load-balancer"
            action = :store_true
        "--lbtimelimit"
            help = "Time limit in seconds of the load-balancing optimization"
            default = 1.0
            arg_type = Float64
            range_tester = (x) -> 0 < x
        "--lbminimprovement"
            help = "Minimum improvement required for the load-balancer to accept a new solution"
            default = 0.98
            arg_type = Float64
            range_tester = (x) -> x <= 1
        "--lbaggressive"
            help = "[DEPRECATED] Set the load-balancer in aggressive mode, in which it doesn't make workers slower"
            action = :store_true
        "--profilerwindowsize"
            help = "Number of seconds that latency statistics are computed over"
            default = 10
            arg_type = Int
            range_tester = (x) -> 0 < x
        "--profilerminsamples"
            help = "Minimum number of samples required for the profiler to compute mean and variance"
            default = 10
            arg_type = Int
            range_tester = (x) -> 1 < x
        "--randomseed"
            help = "Used to compute a deterministic random seed that is unique to the coordinator and each worker"
            arg_type = UInt
            default = rand(UInt)
        "--nslow"
            help = "Number of workers that are slowed down artificially in each iteration"
            arg_type = Int
            default = 0
            range_tester = (x) -> 0 <= x
        "--slowprob"
            help = "Each worker is artificially slowed down with this probability (works in addition to --nslow)"
            arg_type = Float64
            default = 0.0
            range_tester = (x) -> 0 <= x <= 1
    end

    # optionally add implementation-specific arguments
    if @isdefined update_argsettings!
        update_argsettings!(s)
    end

    # common parsing
    parsed_args = parse_args(s, as_symbols=true)

    # store the number of workers
    nworkers = MPI.Comm_size(comm) - 1
    parsed_args[:nworkers] = nworkers

    # optional implementation-specific parsing
    if @isdefined update_parsed_args!
        update_parsed_args!(s, parsed_args)
    end

    return parsed_args
end

"""

Spin on `time_ns()` for the specified number of seconds before returning.
"""
function highres_sleep(seconds)
    t0 = time_ns()
    elapsed = 0.0
    while elapsed < seconds
        elapsed = (time_ns() - t0) / 1e9
    end
    return
end

"""

Main loop run by each worker.
"""
function worker_loop(localdata, recvbuf, sendbuf; nslow::Integer, slowprob::Real, nworkers::Integer, kwargs...)

    @info "Worker $rank is running on $(gethostname())"

    # Markov model artificial latency
    ## timestamp of the most recent Markov model state transition
    transition_time = time_ns()

    ## model
    burst_latency_increase = 1.103
    transition_interval = 10 # seconds
    P = [0.989 0.011; 0.286 0.714]

    ## current state
    ## (initialized according to the steady state distribution)
    burst_state = 1
    if rand() < 0.037037
        burst_state = 2
    end

    # control channel, to tell the workers when to exit
    crreq = MPI.Irecv!(zeros(1), root, control_tag, comm)

    # first iteration (initializes state)
    rreq = MPI.Irecv!(recvbuf, root, data_tag, comm)
    index, _ = MPI.Waitany!([crreq, rreq])
    if index == 1 # exit message on control channel
        return
    end
    t = @elapsed state = worker_task!(recvbuf, sendbuf, localdata; nworkers, kwargs...)
    reinterpret(Float64, view(sendbuf, 1:COMMON_BYTES))[1] = t # send the recorded compute latency to the coordinator
    MPI.Isend(sendbuf, root, data_tag, comm)

    # remaining iterations
    i = 1
    start_time = time_ns()
    while true

        # # make a state transition every transition_interval seconds
        # current_time = time_ns()
        # if (current_time - transition_time) / 1e9 >= transition_interval
        #     transition_time = current_time
        #     if rand() < P[burst_state, 1]
        #         if burst_state == 2
        #             @info "worker $rank transitioning 2 => 1"
        #         end
        #         burst_state = 1
        #     else
        #         if burst_state == 1
        #             @info "worker $rank transitioning 1 => 2"
        #         end
        #         burst_state = 2
        #     end
        # end

        # worker task
        rreq = MPI.Irecv!(recvbuf, root, data_tag, comm)
        index, _ = MPI.Waitany!([crreq, rreq])
        if index == 1 # exit message on control channel
            MPI.Cancel!(rreq) # mark the rreq for cancellation
            MPI.Test!(rreq) # cleanup the data receive request
            break
        end
        # t = @elapsed state = worker_task!(recvbuf, sendbuf, localdata; state, nworkers, kwargs...)
        t = @elapsed begin
            t0 = @elapsed state = worker_task!(recvbuf, sendbuf, localdata; state, nworkers, kwargs...)

            # simulate using a slower computer
            latency_increase_factor = 1
            highres_sleep(t0 * (latency_increase_factor - 1))

            # introduce artificial latency variation between workers
            seconds = 0.0
            if (time_ns() - start_time) / 1e9 < 1 || rank < 40
                seconds = t0 * latency_increase_factor * (0.4 * rank / nworkers)
            end
            highres_sleep(seconds)

            # # during the first 25 iterations (corresponding to one pass over the data), slow down
            # # the first 9 workers by 25% more to simulate a latency burst
            # if i < 25 && rank < 10
            #     highres_sleep((t0 + t0 * 0.4 * rank / nworkers) * 0.25)
            # end

            # # if in the latency burst state, sleep to make up the latency difference
            # if burst_state == 2
            #     highres_sleep(t0 * (burst_latency_increase - 1))
            # end

            # # the nslow first workers are artificially slowed down
            # # (counted as comp. delay)
            # if rank <= nslow
            #     highres_sleep(t0)
            # end

            # # workers 1-3 are slow for the first 100 iterations
            # # the remaining workers are always slow
            # if i < 100 || rank > 3
            #     highres_sleep(t0)
            # end

            # # workers 4-6 become even slower after 50 iterations
            # if i >= 50 && 3 < rank <= 6
            #     highres_sleep(t0)
            # end
        end

        # # workers are artificially slowed down with prob. slowprob.
        # # (counted as comm. delay)
        # if !iszero(slowprob) && rand() < slowprob
        #     highres_sleep(t0)
        # end

        # send response to coordinator
        reinterpret(Float64, view(sendbuf, 1:COMMON_BYTES))[1] = t # send the recorded compute latency to the coordinator
        MPI.Send(sendbuf, root, data_tag, comm)
        i += 1
    end
    return
end

"""

Main function run by each worker.
"""
function worker_main()
    nworkers = MPI.Comm_size(comm) - 1
    parsed_args = parse_commandline(isroot)
    Random.seed!(parsed_args[:randomseed] + rank)
    try
        localdata, recvbuf, sendbuf = worker_setup(rank, nworkers; parsed_args...)
        GC.gc()
        GC.enable(parsed_args[:enablegc])
        @info "(rank $rank) ready"
        MPI.Barrier(comm)
        worker_loop(localdata, recvbuf, sendbuf; parsed_args...)
    catch e
        print(Base.stderr, "rank $rank exiting due to $e")
        exit(0) # only the root exits with non-zero status in case of error
    end
end

"""

Send a stop signal to each worker.
"""
function shutdown(pool::MPIAsyncPool)
    for i in pool.ranks
        MPI.Send(zeros(1), i, control_tag, comm)
    end
end

"""

Select the next partition to process for the `i`-th worker.
"""
function select_partitions!(partition_indices::AbstractVector{<:Integer}, nsubpartitions_all::AbstractVector{<:Integer}, i::Integer)
    0 < i <= length(partition_indices) || throw(DimensionMismatch("i is $i, but partition_indices has dimension $(length(partition_indices))"))
    0 < i <= length(nsubpartitions_all) || throw(DimensionMismatch("i is $i, but nsubpartitions_all has dimension $(length(nsubpartitions_all))"))
    0 < nsubpartitions_all[i] || throw(ArgumentError("Invalid number of partitions for worker $i: $nsubpartitions_all"))
    partition_indices[i] = mod(partition_indices[i], nsubpartitions_all[i]) + 1
end

"""

For each worker, choose the sub-partition it should process.
"""
function select_partitions!(partition_indices::AbstractVector{<:Integer}, nsubpartitions_all::AbstractVector{<:Integer})
    length(partition_indices) == length(nsubpartitions_all) || throw(DimensionMismatch("partition_indices has dimension $(length(partition_indices)), but nsubpartitions_all has dimension $(size(nsubpartitions_all))"))
    for i in 1:length(partition_indices)
        select_partitions!(partition_indices, nsubpartitions_all, i)
    end
    partition_indices
end

"""

For each worker, write the number of sub-partitions and the index of the sub-partition to process to `sendbuf`.
"""
function write_partitions!(sendbuf::AbstractVector{UInt8}; partition_indices::AbstractVector{<:Integer}, nsubpartitions_all::AbstractVector{<:Integer})
    length(partition_indices) == length(nsubpartitions_all) || throw(DimensionMismatch("partition_indices has dimension $(length(partition_indices)), but nsubpartitions_all has dimension $(size(nsubpartitions_all))"))
    nworkers = length(partition_indices)
    s = sizeof(UInt16) * 2 * nworkers
    s <= length(sendbuf) || throw(DimensionMismatch("sendbuf is of length $(length(sendbuf)), but nworkers is $nworkers"))
    vs = reinterpret(Tuple{UInt16,UInt16}, view(sendbuf, 1:s))
    for i in 1:nworkers
        vs[i] = (nsubpartitions_all[i], 0) # sub-partition_index is chosen by the workers
    end
    sendbuf
end

"""

Return the metadata recorded in a buffer received from a worker.
"""
function metadata_from_recvbuf(recvbuf::AbstractVector{UInt8})
    metadata = reinterpret(UInt16, metadata_view(recvbuf))
    length(metadata) == 4 || throw(ArgumentError("metadata is $metadata"))
    canary, worker_rank, nsubpartitions, subpartition_index = metadata
    canary == CANARY_VALUE || throw(ArgumentError("metadata is $metadata"))
    worker_rank, nsubpartitions, subpartition_index
end

"""

Return the compute delay recorded in a buffer received from a worker.
"""
function comp_delay_from_recvbuf(recvbuf::AbstractVector{UInt8})
    reinterpret(Float64, view(recvbuf, 1:COMMON_BYTES))[1]
end

"""

Main loop run by the coordinator.
"""
function coordinator_main()

    # setup
    latency_tol = 0.02
    parsed_args = parse_commandline(isroot)
    niterations::Int = parsed_args[:niterations]
    saveiterates::Bool = parsed_args[:saveiterates]
    nworkers::Int = parsed_args[:nworkers]
    nwait::Int = parsed_args[:nwait]
    loadbalance::Bool = parsed_args[:loadbalance]
    nsubpartitions_all = fill(parsed_args[:nsubpartitions], nworkers)
    partition_indices = zeros(UInt16, nworkers)
    to_worker_common_bytes = sizeof(UInt16) * 2 * nworkers
    prev_repochs = zeros(Int, nworkers)
    Random.seed!(parsed_args[:randomseed])
    @info "Coordinator is running with $(Threads.nthreads()) threads on $(gethostname())"
    @info "Main task is running on thread $(Threads.threadid())"

    # create the output directory if it doesn't exist, and make sure we can write to the output file
    mkpath(dirname(parsed_args[:outputfile]))
    close(h5open(parsed_args[:outputfile], "w"))

    # worker pool and communication buffers
    pool = MPIAsyncPool(nworkers)
    V, recvbuf, sendbuf, state, gradient_recvf = coordinator_setup(nworkers; parsed_args...)
    mod(length(recvbuf), nworkers) == 0 || error("the length of recvbuf must be divisible by the number of workers")
    ∇ = similar(V)
    ∇ .= 0
    isendbuf = similar(sendbuf, nworkers*length(sendbuf))
    irecvbuf = similar(recvbuf)

    # views into recvbuf corresponding to each worker
    recvbufs = [view(recvbuf, partition(length(recvbuf), nworkers, i)) for i in 1:nworkers]

    # number of iterates to save
    nsavediterates::Int = parsed_args[:nsavediterates]
    savediterateindices = unique(round.(Int, exp.(range(log(1), log(niterations), length=nsavediterates))))
    nsavediterates = length(savediterateindices)
    savediterateindex = 1

    # optionally store intermediate iterates
    if saveiterates
        iterates = zeros(Float64, size(V)..., nsavediterates)
    else
        iterates = zeros(Float64, size(V)..., 0)
    end

    # benchmark/analysis data
    responded = zeros(Int, nworkers, niterations)       # which workers respond in each iteration
    latency = zeros(nworkers, niterations)              # per-iteration latency of each worker
    compute_latency = zeros(nworkers, niterations)      # per-iteration computation latency of each worker
    ts_compute = zeros(niterations)                     # overall per-iteration latency
    ts_update = zeros(niterations)                      # per-iteration latency of the coordinator update
    loadbalanced = zeros(niterations)                   # indicates which iterations the load-balancer made changes in

    # 2-argument fwait needed for asyncmap!
    f = (epoch, repochs) -> fwait(epoch, repochs; parsed_args...)

    # 4-argument recvf needed for asyncmap!
    recvf = isnothing(gradient_recvf) ? nothing : (args...) -> gradient_recvf(args...; state, parsed_args...)

    # manually call the GC now, and optionally turn off GC, to avoid pauses later during execution
    GC.gc()
    GC.enable(parsed_args[:enablegc])

    # setup the latency profiler
    profiler_chin, profiler_chout = CodedComputing.setup_profiler_channels()
    profiler_task = Threads.@spawn CodedComputing.latency_profiler(
        profiler_chin, profiler_chout;
        nworkers,
        windowsize=Second(parsed_args[:profilerwindowsize]),
        minsamples=parsed_args[:profilerminsamples],
        )

    # vector of communication latencies for each worker
    # (to use instead of the measured value for stale results)
    communication_latency = fill(NaN, nworkers)
    communication_latency_samples = CircularBuffer{Float64}(nworkers)

    # setup the load-balancer
    _, loadbalancer_chout = CodedComputing.setup_loadbalancer_channels()
    loadbalancer_task = Threads.@spawn CodedComputing.load_balancer(
        profiler_chout, loadbalancer_chout;
        nsubpartitions=parsed_args[:nsubpartitions],
        nwait, nworkers,
        min_improvement=parsed_args[:lbminimprovement],
        )

    # ensure all workers have finished compiling before starting the computation
    # (this is only necessary when benchmarking)
    @info "Waiting for workers to get ready..."
    MPI.Barrier(comm)
    @info "Optimization starting"

    # remaining iterations
    for epoch in 1:niterations

        # set to NaN to remove old latency results
        pool.latency .= NaN

        # gradient_state, iterate_state = state
        # ∇i, tg = gradient_state
        # @info "epoch $epoch, tg.ninit / tg.n: $(tg.ninit / tg.n)"
        # @info "epoch $epoch"

        ## force an error if the profiler or load-balancer task has failed
        if istaskfailed(profiler_task)
            wait(profiler_task)
        end
        if istaskfailed(loadbalancer_task)
            wait(loadbalancer_task)
        end

        ## get updated partitioning from the load-balancer
        while isready(loadbalancer_chout)
            vout = take!(loadbalancer_chout)
            0 < vout.worker <= nworkers || throw(ArgumentError("incorrect $vout"))
            0 < vout.p || throw(ArgumentError("incorrect $vout"))
            if nsubpartitions_all[vout.worker] != vout.p
                # @info "epoch $epoch, worker $(vout.worker) nsubpartitions $(nsubpartitions_all[vout.worker]) => $(vout.p)"
                nsubpartitions_all[vout.worker] = vout.p
            end
            loadbalanced[epoch] = true
        end

        ## update partitioning
        select_partitions!(partition_indices, nsubpartitions_all)
        write_partitions!(sendbuf; partition_indices, nsubpartitions_all)

        ## worker task
        ts_compute[epoch] = @elapsed begin
            repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm; nwait, epoch=epoch, tag=data_tag, recvf, latency_tol)
        end
        responded[:, epoch] .= repochs

        ## record latency
        latency[:, epoch] .= pool.latency
        for i in 1:nworkers
            t = repochs[i] == epoch ? reinterpret(Float64, view(recvbufs[i], 1:COMMON_BYTES))[1] : NaN
            compute_latency[i, epoch] = t
        end

        # record new communication latency values
        for i in 1:nworkers
            if repochs[i] == epoch # fresh result
                communication_latency[i] = latency[i, epoch] - compute_latency[i, epoch]
                push!(communication_latency_samples, communication_latency[i])
            end
        end

        ## send any new latency information to the profiling sub-system
        if loadbalance
            communication_latency_mean = mean(communication_latency_samples)
            timestamp = Time(now())
            for i in 1:nworkers
                if repochs[i] != prev_repochs[i]
                    _, nsubpartitions, subpartition_index = metadata_from_recvbuf(recvbufs[i])
                    comp_latency = comp_delay_from_recvbuf(recvbufs[i])
                    comm_latency = communication_latency[i]
                    if repochs[i] != epoch # don't trust comm. latency for stale results
                        comm_latency = communication_latency_mean
                    end
                    v = CodedComputing.ProfilerInput(i, 1/nworkers, 1/nsubpartitions, timestamp, comp_latency, comm_latency)
                    push!(profiler_chin, v)
                end
                prev_repochs[i] = repochs[i]
            end
        end

        ## coordinator task
        ts_update[epoch] = @elapsed begin
            state = coordinator_task!(V, ∇, recvbufs, sendbuf, epoch, repochs; state, parsed_args...)
        end

        ## optinally record the iterate
        if saveiterates && epoch == savediterateindices[savediterateindex]
            ndims = length(size(V))
            selectdim(iterates, ndims+1, savediterateindex) .= V
            savediterateindex += 1
        end
    end

    @info "Optimization finished; writing output to disk"

    # write output
    h5open(parsed_args[:outputfile], "w") do fid

        # write program parameters
        for (key, val) in parsed_args
            if isnothing(val)
                continue
            end
            fid["parameters/$key"] = val
        end

        # write the computation output
        fid[parsed_args[:outputdataset]] = V

        # optionally save all iterates
        if saveiterates
            fid["iterates"] = iterates
            fid["savediterateindices"] = savediterateindices
        end

        # write benchmark data
        fid["benchmark/t_compute"] = ts_compute
        fid["benchmark/t_update"] = ts_update
        fid["benchmark/responded"] = responded
        fid["benchmark/latency"] = latency
        fid["benchmark/compute_latency"] = compute_latency
        fid["benchmark/loadbalanced"] = loadbalanced
    end
    @info "Output written to disk; stopping tasks"

    # signal all workers to stop
    waitall!(pool, recvbuf, irecvbuf)
    shutdown(pool)

    # stop the profiler and load-balancer
    close(profiler_chin)
    close(profiler_chout)
    close(loadbalancer_chout)

    @info "Tasks stopped; exiting"
    return
end

if isroot

    # check that all necessary functions are defined
    (@isdefined coordinator_setup) || error("function coordinator_setup must be defined")
    (@isdefined coordinator_task!) || error("function coordinator_task! must be defined")
    (@isdefined worker_setup) || error("function worker_setup must be defined")
    (@isdefined worker_task!) || error("function worker_task! must be defined")
    (@isdefined fwait) || error("function fwait must be defined")

    # start the coordinator loop
    coordinator_main()
else
    worker_main()
end
MPI.Barrier(comm)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to TEME transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/transformations/teme/teme.jl
# ========================================

# Functions rTEMEtoTOD and rTODtoTEME
# -----------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#   v_tod =   -4.7460885670 i +    0.7860772220 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Functions rTEMEtoTOD and rTODtoTEME" begin
    JD_TT   = 2453101.828154745

    ## rTEMEtoTOD
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = rTEMEtoTOD(JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_tod = D_TOD_TEME*r_teme
    v_tod = D_TOD_TEME*v_teme

    @test r_tod[1] ≈ +5094.51620300 atol=1e-7
    @test r_tod[2] ≈ +6127.36527840 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460883850  atol=1e-9
    @test v_tod[2] ≈ +0.7860783240  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_TEME = rTEMEtoTOD(Quaternion,
                            JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_tod = vect(conj(q_TOD_TEME)*r_teme*q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME)*v_teme*q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-7
    @test r_tod[2] ≈ +6127.36527840 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460883850  atol=1e-9
    @test v_tod[2] ≈ +0.7860783240  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = rTEMEtoTOD(JD_TT)

    r_tod = D_TOD_TEME*r_teme
    v_tod = D_TOD_TEME*v_teme

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_TEME = rTEMEtoTOD(Quaternion, JD_TT)

    r_tod = vect(conj(q_TOD_TEME)*r_teme*q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME)*v_teme*q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## rTODtoTEME
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = rTODtoTEME(JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_teme = D_TEME_TOD*r_tod
    v_teme = D_TEME_TOD*v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_TOD = rTODtoTEME(Quaternion,
                            JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_teme = vect(conj(q_TEME_TOD)*r_tod*q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD)*v_tod*q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]
    v_tod = [-4.7460885670; 0.7860772220; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = rTODtoTEME(JD_TT)

    r_teme = D_TEME_TOD*r_tod
    v_teme = D_TEME_TOD*v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_TOD = rTODtoTEME(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_TOD)*r_tod*q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD)*v_tod*q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions rTEMEtoMOD and rMODtoTEME
# -----------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#   v_mod =   -4.7462624950 i +    0.7860141490 j +    5.5317910250 k [km/s]
#
################################################################################

@testset "Functions rTEMEtoMOD and rMODtoTEME" begin
    JD_TT   = 2453101.828154745

    ## rTEMEtoMOD
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = rTEMEtoMOD(JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_mod = D_MOD_TEME*r_teme
    v_mod = D_MOD_TEME*v_teme

    @test r_mod[1] ≈ +5094.02837450 atol=1e-7
    @test r_mod[2] ≈ +6127.87081640 atol=1e-7
    @test r_mod[3] ≈ +6380.24851640 atol=1e-7

    @test v_mod[1] ≈ -4.7462630520  atol=1e-9
    @test v_mod[2] ≈ +0.7860140450  atol=1e-9
    @test v_mod[3] ≈ +5.5317905620  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TEME = rTEMEtoMOD(Quaternion,
                            JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_mod = vect(conj(q_MOD_TEME)*r_teme*q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME)*v_teme*q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-7
    @test r_mod[2] ≈ +6127.87081640 atol=1e-7
    @test r_mod[3] ≈ +6380.24851640 atol=1e-7

    @test v_mod[1] ≈ -4.7462630520  atol=1e-9
    @test v_mod[2] ≈ +0.7860140450  atol=1e-9
    @test v_mod[3] ≈ +5.5317905620  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = rTEMEtoMOD(JD_TT)

    r_mod = D_MOD_TEME*r_teme
    v_mod = D_MOD_TEME*v_teme

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TEME = rTEMEtoMOD(Quaternion, JD_TT)

    r_mod = vect(conj(q_MOD_TEME)*r_teme*q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME)*v_teme*q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## rMODtoTEME
    ## ==========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_TEME_MOD = rMODtoTEME(JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_teme = D_TEME_MOD*r_mod
    v_teme = D_TEME_MOD*v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_MOD = rMODtoTEME(Quaternion,
                            JD_TT,
                            -0.003875*pi/(180*3600),
                            -0.052195*pi/(180*3600))

    r_teme = vect(conj(q_TEME_MOD)*r_mod*q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD)*v_mod*q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]
    v_mod = [-4.7462624950; 0.7860141490; 5.5317910250]

    ## DCM
    ## ---

    D_TEME_MOD = rMODtoTEME(JD_TT)

    r_teme = D_TEME_MOD*r_mod
    v_teme = D_TEME_MOD*v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_MOD = rMODtoTEME(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_MOD)*r_mod*q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD)*v_mod*q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions rTEMEtoGCRF and rGCRFtoTEME
# -------------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = -0.003875"
#   δΔψ   = -0.052195"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_gcrf = 5102.50895790  i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =  -4.7432201570 i + 0.7905364970    j + 5.5337557270    k [km/s]
#
# Furthermore, using:
#
#   JD_TT = 2453101.828154745
#   δΔϵ   = 0"
#   δΔψ   = 0"
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
################################################################################

@testset "Functions rTEMEtoGCRF and rGCRFtoTEME" begin
    JD_TT   = 2453101.828154745

    ## rTEMEtoGCRF
    ## ===========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_GCRF_TEME = rTEMEtoGCRF(JD_TT,
                              -0.003875*pi/(180*3600),
                              -0.052195*pi/(180*3600))

    r_gcrf = D_GCRF_TEME*r_teme
    v_gcrf = D_GCRF_TEME*v_teme

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-7
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-7
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-7

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-9

    ## Quaternion
    ## ----------

    q_GCRF_TEME = rTEMEtoGCRF(Quaternion,
                              JD_TT,
                              -0.003875*pi/(180*3600),
                              -0.052195*pi/(180*3600))

    r_gcrf = vect(conj(q_GCRF_TEME)*r_teme*q_GCRF_TEME)
    v_gcrf = vect(conj(q_GCRF_TEME)*v_teme*q_GCRF_TEME)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-7
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-7
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-7

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-9
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-9
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_J2000_TEME = rTEMEtoGCRF(JD_TT)

    r_j2000 = D_J2000_TEME*r_teme
    v_j2000 = D_J2000_TEME*v_teme

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## Quaternion
    ## ----------

    q_J2000_TEME = rTEMEtoGCRF(Quaternion, JD_TT)

    r_j2000 = vect(conj(q_J2000_TEME)*r_teme*q_J2000_TEME)
    v_j2000 = vect(conj(q_J2000_TEME)*v_teme*q_J2000_TEME)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## rGCRFtoTEME
    ## ===========

    ## -------------------------------------------------------------------------
    ##                             First Test
    ## -------------------------------------------------------------------------

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_TEME_GCRF = rGCRFtoTEME(JD_TT,
                              -0.003875*pi/(180*3600),
                              -0.052195*pi/(180*3600))

    r_teme = D_TEME_GCRF*r_gcrf
    v_teme = D_TEME_GCRF*v_gcrf

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_GCRF = rGCRFtoTEME(Quaternion,
                              JD_TT,
                              -0.003875*pi/(180*3600),
                              -0.052195*pi/(180*3600))

    r_teme = vect(conj(q_TEME_GCRF)*r_gcrf*q_TEME_GCRF)
    v_teme = vect(conj(q_TEME_GCRF)*v_gcrf*q_TEME_GCRF)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## -------------------------------------------------------------------------
    ##                            Second Test
    ## -------------------------------------------------------------------------

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_TEME_J2000 = rGCRFtoTEME(JD_TT)

    r_teme = D_TEME_J2000*r_j2000
    v_teme = D_TEME_J2000*v_j2000

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_J2000 = rGCRFtoTEME(Quaternion, JD_TT)

    r_teme = vect(conj(q_TEME_J2000)*r_j2000*q_TEME_J2000)
    v_teme = vect(conj(q_TEME_J2000)*v_j2000*q_TEME_J2000)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end

# Functions rTEMEtoPEF and rPEFtoTEME
# -----------------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400
#   LOD    = 0.0015563 s
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# one gets the following data:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_pef  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Functions rTEMEtoPEF and rPEFtoTEME" begin
    JD_UT1 = DatetoJD(2004,4,6,7,51,28.386009) - 0.4399619/86400
    LOD    = 0.0015563
    w      = 7.292115146706979e-5*(1-LOD/86400)

    ## rTEMEtoPEF
    ## ==========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_PEF = rTEMEtoPEF(JD_UT1)

    r_pef = D_TEME_PEF*r_teme
    v_pef = D_TEME_PEF*v_teme - [0;0;w] × r_pef

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7

    @test v_pef[1] ≈ -3.2256327470  atol=1e-9
    @test v_pef[2] ≈ -2.8724425110  atol=1e-9
    @test v_pef[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_PEF = rTEMEtoPEF(Quaternion, JD_UT1)

    r_pef = vect(conj(q_TEME_PEF)*r_teme*q_TEME_PEF)
    v_pef = vect(conj(q_TEME_PEF)*v_teme*q_TEME_PEF) - [0;0;w] × r_pef

    @test r_pef[1] ≈ -1033.47503130 atol=1e-7
    @test r_pef[2] ≈ +7901.30558560 atol=1e-7
    @test r_pef[3] ≈ +6380.34453270 atol=1e-7

    @test v_pef[1] ≈ -3.2256327470  atol=1e-9
    @test v_pef[2] ≈ -2.8724425110  atol=1e-9
    @test v_pef[3] ≈ +5.5319312880  atol=1e-9

    ## rPEFtoTEME
    ## ==========

    r_pef  = [-1033.47503130; 7901.30558560; 6380.34453270]
    v_pef  = [-3.2256327470; -2.8724425110; +5.5319312880]

    ## DCM
    ## ---

    D_PEF_TEME = rPEFtoTEME(JD_UT1)

    r_teme = D_PEF_TEME*r_pef
    v_teme = D_PEF_TEME*(v_pef + [0;0;w] × r_pef)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_PEF_TEME = rPEFtoTEME(Quaternion, JD_UT1)

    r_teme = vect(conj(q_PEF_TEME)*r_pef*q_PEF_TEME)
    v_teme = vect(conj(q_PEF_TEME)*(v_pef + [0;0;w] × r_pef)*q_PEF_TEME)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Some{T}

A wrapper type used in `Union{Some{T}, Nothing}` to distinguish between the absence
of a value ([`nothing`](@ref)) and the presence of a `nothing` value (i.e. `Some(nothing)`).

Use [`something`](@ref) to access the value wrapped by a `Some` object.
"""
struct Some{T}
    value::T
end

promote_rule(::Type{Some{T}}, ::Type{Some{S}}) where {T, S<:T} = Some{T}
promote_rule(::Type{Some{T}}, ::Type{Nothing}) where {T} = Union{Some{T}, Nothing}

convert(::Type{Some{T}}, x::Some) where {T} = Some{T}(convert(T, x.value))
convert(::Type{Union{Some{T}, Nothing}}, x::Some) where {T} = convert(Some{T}, x)

convert(::Type{Union{T, Nothing}}, x::Any) where {T} = convert(T, x)
convert(::Type{Nothing}, x::Any) = throw(MethodError(convert, (Nothing, x)))
convert(::Type{Nothing}, x::Nothing) = nothing

function show(io::IO, x::Some)
    if get(io, :typeinfo, Any) == typeof(x)
        show(io, x.value)
    else
        print(io, "Some(")
        show(io, x.value)
        print(io, ')')
    end
end

"""
    notnothing(x)

Throw an error if `x === nothing`, and return `x` if not.
"""
notnothing(x::Any) = x
notnothing(::Nothing) = throw(ArgumentError("nothing passed to notnothing"))

"""
    something(x, y...)

Return the first value in the arguments which is not equal to [`nothing`](@ref),
if any. Otherwise throw an error.
Arguments of type [`Some`](@ref) are unwrapped.

# Examples
```jldoctest
julia> something(nothing, 1)
1

julia> something(Some(1), nothing)
1

julia> something(missing, nothing)
missing

julia> something(nothing, nothing)
ERROR: ArgumentError: No value arguments present
```
"""
function something end

something() = throw(ArgumentError("No value arguments present"))
something(x::Nothing, y...) = something(y...)
something(x::Some, y...) = x.value
something(x::Any, y...) = x

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using BasicFWI
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using DelimitedFiles

# if nworkers()==1
	# addprocs(2)
# end

println("IM here")
#@everywhere begin
	using BasicFWI
	using jInv.Mesh
	using jInv.Utils
	using jInv.ForwardShare
#end

# setup model and attenuation function
nx = 164; nz = 80
domain = [0.0,10.0,0.0,4.0]
Mr = getRegularMesh(domain,[nx,nz])

# the slowness model
m = 3.0*ones(nx,nz);
m[20:70,10:32] .= 5.0; # velocity...
m = 1.0./(m.^2); # convert to slowness squared

# m = readdlm("SEGmodel2Dsalt.dat");
# m = m*1e-3;
# m = m';
# m = 1.0./m^2

# attenuation for BC
padx = 20; padz = 20
a    = 2.0;
xc = getCellCenteredGrid(Mr)
# gammaxL = a*(xc[:,1] .- xc[padx,1]).^2;
# gammaxL[padx+1:end,:] .= 0
# gammaxR = a*(xc[:,1] .- xc[end-padx+1,1]).^2
# gammaxR[1:end-padx,:] .= 0

# gammax = gammaxL + gammaxR

# gammaz = a*(xc[:,2] .- xc[end-padz+1,2]).^2
# gammaz[:,1:end-padz] .= 0

# gamma = gammax + gammaz

gamma = getABL(Mr,true,[padx;padz],a);

# parameters for the Helmholtz (units in km)
h = Mr.h;
n = [nx; nz;]
omega = 2*pi*[1.0;1.2;1.4;1.5;2.0;3.0]
nfreq = length(omega)

# generate sources
q = zeros(tuple(n.+1...)); q[padx+1:2:end-padx-1,1] .= 1e4
Q = sdiag(vec(q))
Q = Q[:,(LinearIndices(sum(Q,dims=2) .!= 0))[findall(sum(Q,dims=2) .!= 0)]]
nsrc = size(Q,2)

# receivers
p = zeros(tuple(n.+1...)); p[padx+1:end-padx-1,1] .= 1
P = sdiag(vec(p))
P = P[:,(LinearIndices(sum(P,dims=2) .!= 0))[findall(sum(P,dims=2) .!= 0)]]
nrec = size(P,2)

pFor  = getFWIparam(omega,gamma,Q,P,Mr)
pForp = getFWIparam(omega,gamma,Q,P,Mr,true)

# inversion mesh and forward mesh are the same here
M2Mp = ones(length(pForp))


text:
using Test
using TightBindingLattice

@testset "symmetryoperation" begin
    iden = IdentityOperation(Int, 2)

    @testset "IdentityOperation" begin
        @test iden == IdentityOperation{Int}(2)

        @test iden != IdentityOperation(Int, 3)
        @test iden * iden == iden
        n=0;  @test iden^n == iden
        n=4;  @test iden^n == iden
        n=-3; @test iden^n == iden
        @test inv(iden) == iden
        @test isidentity(iden)
        @test istranslation(iden)
        @test ispoint(iden)
        @test domaintype(iden) == Int
        @test dimension(iden) == 2
        @test hash(iden) == hash(IdentityOperation(Int, 2))
        @test hash(iden) != hash(IdentityOperation(Int, 3))

        @test apply_operation(iden, [3,2]) == [3,2]
        @test iden([3,2]) == [3,2]

        @test_throws DimensionMismatch iden * IdentityOperation(Int, 3)
        @test_throws DimensionMismatch iden([3,2,1])
        @test_throws DimensionMismatch apply_operation(iden, [3,2,1])
    end

    @testset "TranslationOperation" begin
        t0  = TranslationOperation([0,0])
        t1  = TranslationOperation([1, -1])
        t1p = TranslationOperation{Int}([1, -1])
        t2  = TranslationOperation([2,4])

        @testset "properties" begin
            @test domaintype(t1) == Int
            @test dimension(t1) == 2
            @test dimension(t2) == 2
            @test  isidentity(t0)
            @test !isidentity(t1)
            @test  istranslation(t0)
            @test  istranslation(t1)
            @test  ispoint(t0)
            @test !ispoint(t1)
        end

        @testset "equality" begin
            @test t1 == t1p
            @test t1 !== t1p
            @test hash(t1) == hash(t1p)
            @test hash(t1) != hash(t2)
        end

        @testset "equality-heterotype" begin
            @test t0 == IdentityOperation(Int, 2)
            @test IdentityOperation(Int, 2) == t0
        end

        @testset "operator" begin
            @test t1 * iden == t1
            @test iden * t1 == t1

            @test_throws DimensionMismatch t1 * IdentityOperation(Int, 5)
            @test_throws DimensionMismatch IdentityOperation(Int, 5) * t1

            @test t2^2 == TranslationOperation([4,8])

            t3 = t1*t2
            @test t3.displacement == [3, 3]
            @test t3 == TranslationOperation([3, 3])
            @test inv(t3) == TranslationOperation([-3, -3])

            t4 = TranslationOperation([1,2,3,4])
            @test_throws DimensionMismatch t1 * t4
            @test_throws DimensionMismatch apply_operation(t1, [1,2,3])
            @test_throws DimensionMismatch t1([1,2,3])
        end

        @testset "apply" begin
            @test apply_operation(t1, [5,6]) == [6, 5]
            @test t1([5,6]) == [6, 5]
        end

        @testset "promotion" begin
            @test promote_type(typeof(t1), typeof(iden)) == TranslationOperation{Int}
            arr = TranslationOperation{Int}[]
            push!(arr, IdentityOperation(Int, 2))
            push!(arr, [3, 4])
            @test isa(arr[1], TranslationOperation{Int})
            @test isa(arr[2], TranslationOperation{Int})
            @test arr[1].displacement == [0,0]
            @test arr[2].displacement == [3,4]
        end
    end

    @testset "PointOperation" begin
        p0 = PointOperation([1 0; 0 1])
        p1 = PointOperation([0 -1; 1 -1])
        p1p = PointOperation{Int}([0 -1; 1 -1])
        p2 = PointOperation([0 1; 1 0])

        @testset "properties" begin
            @test domaintype(p0) == Int
            @test dimension(p0) == 2
            @test dimension(p1) == 2
            @test  isidentity(p0)
            @test !isidentity(p1)
            @test  istranslation(p0)
            @test !istranslation(p1)
            @test  ispoint(p0)
            @test  ispoint(p1)
        end

        @testset "equality" begin
            @test p1 == p1p
            @test p1 !== p1p
            @test hash(p1) == hash(p1p)
            @test hash(p1) != hash(p2)
        end

        @testset "equality-heterotype" begin
            @test p0 == iden
            @test iden == p0
            @test p1 != iden
            @test iden != p1
            @test p1 != TranslationOperation([1,0])
            @test p0 == TranslationOperation([0,0])
            @test TranslationOperation([1,0]) != p1
            @test TranslationOperation([0,0]) == p0
        end

        @testset "operator" begin
            @test p1*p2 == PointOperation([0 -1; 1 -1] * [0 1; 1 0])
            @test inv(p1) * p1 == PointOperation([1 0; 0 1])
            @test p1 * inv(p1) == PointOperation([1 0; 0 1])

            @test p1^3 == p1 * p1 * p1
            @test p1^(-1) == inv(p1)
            @test p1^(-2) == inv(p1)*inv(p1)
            @test p2^2 == PointOperation([1 0; 0 1])

            n =  3;  @test p1^n == p1 * p1 * p1
            n = -1;  @test p1^n == inv(p1)
            n = -2;  @test p1^n == inv(p1)*inv(p1)
            n =  2;  @test p2^n == PointOperation([1 0; 0 1])
        end

        @testset "apply" begin
            @test apply_operation(p2, [5, 6]) == [6, 5]
            @test p2([5, 6]) == [6, 5]

            p3 = PointOperation([0 1 0; 1 0 0; 0 0 1])
            @test_throws DimensionMismatch p1*p3
            @test_throws DimensionMismatch apply_operation(p1, [1,2,3])
        end

        @testset "promotion" begin
            @test promote_type(typeof(p1), typeof(iden)) == PointOperation{Int}
            arr = PointOperation{Int}[]
            push!(arr, IdentityOperation(Int, 2))
            push!(arr, [0 1; 1 0])
            @test isa(arr[1], PointOperation{Int})
            @test arr[1].matrix == [1 0; 0 1]
            @test isa(arr[2], PointOperation{Int})
            @test arr[2].matrix == [0 1; 1 0]
        end
    end

    @testset "SpaceOperation" begin
        @test_throws DimensionMismatch SpaceOperation([1 0 0; 0 1 0], [1,2,3])
        @test_throws DimensionMismatch SpaceOperation([1 0; 0 1], [1,2,3])

        let sop = SpaceOperation{Int, Int}(2)
            @test sop.matrix == [1 0; 0 1]
            @test sop.displacement == [0, 0]
        end

        c4p = PointOperation([ 0 -1;  1  0])
        c4m = PointOperation([ 0  1; -1  0])
        m10 = PointOperation([-1  0;  0  1])
        t10 = TranslationOperation([ 1,  0])
        t01 = TranslationOperation([ 0,  1])

        @testset "properties" begin
            s0 = SpaceOperation(Int, Int, 2)
            st = SpaceOperation([1 0; 0 1], [1, 0])
            sp = SpaceOperation([0 1; 1 0], [0, 0])
            s2 = SpaceOperation([0 1; 1 0], [1, 0])

            @test domaintype(s0) == Int
            @test domaintype(st) == Int
            @test domaintype(sp) == Int
            @test domaintype(s2) == Int

            @test  isidentity(s0)
            @test !isidentity(st)
            @test !isidentity(sp)
            @test !isidentity(s2)

            @test  istranslation(s0)
            @test  istranslation(st)
            @test !istranslation(sp)
            @test !istranslation(s2)

            @test  ispoint(s0)
            @test !ispoint(st)
            @test  ispoint(sp)
            @test !ispoint(s2)
        end

        @testset "equality" begin
            sop = SpaceOperation(c4p)
            @test sop.matrix == [0 -1; 1 0] && sop.displacement == [0,0]

            sop = SpaceOperation(t10)
            @test sop.matrix == [1 0; 0 1] && sop.displacement == [1,0]

            sop = SpaceOperation(c4p, t10)
            @test sop.matrix == [0 -1; 1 0] && sop.displacement == [1, 0]
            @test sop == SpaceOperation([0 -1; 1 0], [1, 0])
            @test sop != SpaceOperation([0 -1; 1 0], [0, 1])
            @test sop != SpaceOperation([0  1; 1 0], [0, 1])

            @test hash(sop) == hash(SpaceOperation([0 -1; 1 0], [1,0]))
            @test hash(sop) != hash(SpaceOperation([0  1; 1 0], [1,0]))
            @test hash(sop) != hash(SpaceOperation([0 -1; 1 0], [0,0]))
        end

        @testset "space" begin
            c4p_t10 = c4p*t10
            t10_c4p = t10*c4p

            @test c4p_t10 == SpaceOperation([0 -1; 1 0], [1,  0])
            @test t10_c4p == SpaceOperation([0 -1; 1 0], [0, -1])

            @test c4p * c4p_t10 == SpaceOperation([-1 0; 0 -1], [ 1,  0])
            @test c4p * t10_c4p == SpaceOperation([-1 0; 0 -1], [ 0, -1])
            @test c4p_t10 * c4p == SpaceOperation([-1 0; 0 -1], [ 0, -1])
            @test t10_c4p * c4p == SpaceOperation([-1 0; 0 -1], [-1,  0])

            # associativity
            @test all((A * B) * C == A * (B * C)
                        for A in AbstractSpaceSymmetryOperation[t10, t01, m10, c4p, c4m],
                            B in AbstractSpaceSymmetryOperation[t10, t01, m10, c4p, c4m],
                            C in AbstractSpaceSymmetryOperation[t10, t01, m10, c4p, c4m])

            # inverse
            @test all(inv(A * B) == inv(B) * inv(A)
                          for A in AbstractSpaceSymmetryOperation[t10, t01, m10, c4p, c4m],
                              B in AbstractSpaceSymmetryOperation[t10, t01, m10, c4p, c4m])
        end

        @testset "comparison with other types" begin
            @test SpaceOperation(c4p) == c4p
            @test c4p == SpaceOperation(c4p)
            @test isidentity(SpaceOperation(c4p * c4m))
            @test SpaceOperation(c4p * c4m) == IdentityOperation(Int, 2)
            @test IdentityOperation(Int, 2) == SpaceOperation(c4p * c4m)

            @test m10 * m10 * t10 == t10
            @test c4p * t10 * c4m == TranslationOperation([0, 1])
            @test TranslationOperation([0, 1]) == c4p * t10 * c4m
            @test t10 * c4p * c4p * t10 == PointOperation([-1 0; 0 -1])
            @test PointOperation([-1 0; 0 -1]) == t10 * c4p * c4p * t10
        end

        @testset "promotion" begin
            @test promote_type(IdentityOperation{Int}, SpaceOperation{Int, Int}) == SpaceOperation{Int, Int}
            @test promote_type(TranslationOperation{Int}, SpaceOperation{Int, Int}) == SpaceOperation{Int, Int}
            @test promote_type(PointOperation{Int}, SpaceOperation{Int, Int}) == SpaceOperation{Int, Int}

            @test promote_type(IdentityOperation{Int}, SpaceOperation{Int, Rational{Int}}) == SpaceOperation{Int, Rational{Int}}
            @test promote_type(TranslationOperation{Int}, SpaceOperation{Int, Rational{Int}}) == SpaceOperation{Int, Rational{Int}}
            @test promote_type(PointOperation{Int}, SpaceOperation{Int, Rational{Int}}) == SpaceOperation{Int, Rational{Int}}

            @test promote_type(SpaceOperation{Int, Int}, SpaceOperation{Int, Rational{Int}}) == SpaceOperation{Int, Rational{Int}}

            arr = SpaceOperation{Int, Rational{Int}}[]
            push!(arr, PointOperation([0 1; 1 0]))
            push!(arr, TranslationOperation([1, 0]))
            push!(arr, IdentityOperation(Int, 2))
            push!(arr, [1,0])
            push!(arr, [0 1; 1 0])

            @test arr[1].matrix == [0 1; 1 0] && arr[1].displacement == [0,0]
            @test arr[2].matrix == [1 0; 0 1] && arr[2].displacement == [1,0]
            @test arr[3].matrix == [1 0; 0 1] && arr[3].displacement == [0,0]
            @test arr[4].matrix == [1 0; 0 1] && arr[4].displacement == [1,0]
            @test arr[5].matrix == [0 1; 1 0] && arr[5].displacement == [0,0]

        end

        @testset "apply" begin
            t = TranslationOperation([1, 0])
            p = PointOperation([0 1; 1 0])
            tp = t * p
            pt = p * t

            @test tp^3 == t * p * t * p * t * p
            @test apply_operation(tp, [5,0]) == [1, 5]
            @test apply_operation(pt, [5,0]) == [0, 6]
            @test tp([5,0]) == [1, 5]
            @test pt([5,0]) == [0, 6]
            # @test tp([5,0]) == tpc([5,0])
            # @test pt([5,0]) == ptc([5,0])
            @test tp^0 == IdentityOperation(Int, 2)
            @test tp^1 == tp
            @test tp^2 == t * p * t * p
            @test tp^(-2) == inv(p) * inv(t) * inv(p) * inv(t)

            n = 0;   @test tp^n == IdentityOperation(Int, 2)
            n = 1;   @test tp^n == tp
            n = 2;   @test tp^n == t * p * t * p
            n = -2;  @test tp^n == inv(p) * inv(t) * inv(p) * inv(t)
        end
    end

    # @testset "product" begin
    #     r0 = ProductOperation()

    #     t = TranslationOperation([2, 4])
    #     p = PointOperation([0 -1; 1 -1])

    #     rt = ProductOperation(t)
    #     rp = ProductOperation{Int}(p)

    #     @test r0 * t == rt
    #     @test t * r0 == rt

    #     tp = t * p
    #     pt = p * t

    #     tpc = canonize(tp)
    #     ptc = canonize(pt)

    #     @test tp^3 == t * p * t * p * t * p
    #     @test tpc.factors[1] == ptc.factors[1]
    #     @test isa(tpc.factors[1], PointOperation) && isa(tpc.factors[2], TranslationOperation)
    #     @test isa(ptc.factors[1], PointOperation) && isa(ptc.factors[2], TranslationOperation)
    #     @test apply_operation(tp, [5,0]) == apply_operation(tpc, [5,0])
    #     @test apply_operation(pt, [5,0]) == apply_operation(ptc, [5,0])
    #     @test tp([5,0]) == tpc([5,0])
    #     @test pt([5,0]) == ptc([5,0])

    #     @test tp^0 == ProductOperation()
    #     @test tp^1 == tp
    #     @test tp^2 == t * p * t * p
    #     @test tp^(-2) == inv(p) * inv(t) * inv(p) * inv(t)

    #     n = 0
    #     @test tp^n == ProductOperation()
    #     n = 1
    #     @test tp^n == tp
    #     n = 2
    #     @test tp^n == t * p * t * p
    #     n = -2
    #     @test tp^n == inv(p) * inv(t) * inv(p) * inv(t)

    #     @test canonize(tp^3 * inv(tp^3)) == IdentityOperation()
    #     @test iscanonical(pt)
    #     @test !iscanonical(tp)
    #     @test iscanonical(tpc)
    #     @test iscanonical(ptc)

    #     pp = ProductOperation(p, p)
    #     @test pp != p * p
    #     @test canonize(pp) == p*p

    #     @test domaintype(pp) == Int
    #     @test dimension(pp) == 2
    # end
end


text:
using Distributions
using Random
using LinearAlgebra
using Base.Threads
using Base.Iterators

"""
    MeasurementModel

`MeasurementModel` is an abstract type representing a vector function of a vector input
with correlated uncertainties.  A `MeasurementModel` is responsible for calculating
vector function outputs and the associated Jacobian matrix both evaluated at the
input values.

To emphasize that a MeasurementModel represents a vector function that takes an
UncertainValues object representing the input values into an UncertainValues object
representing the output values, the function operator has been overloaded:

    (mm::MeasurementModel)(uvs::UncertainValues)::UncertainValues

or

    (mm::MeasurementModel)(uvs::Dict{<:Label, UncertainValue})::UncertainValues

This helps you to think of a MeasurementModel as the equivalent of a function for
UncertainValues.

Furthermore, you can compose / pipe / chain `MeasurementModel`s together to implement a
multi-step calculation using the ∘ operator.  If `mm1` and `mm2` are `MeasurementModel`s
such that the inputs for `mm2` come from `mm1` (and the original inputs) then:

    (mm2 ∘ mm1)(inputs) = mm2(mm1(inputs))

Equivalently, if `mm3` and `mm4` both take the same inputs they may be computed in
parallel using:

    (mm3 | mm4)(inputs) = cat(mm3(inputs), mm4(inputs))

You can combine `∘` and `|` together like this:

    (mm2 ∘ (mm3 | mm4) ∘ mm1)(inputs) = mm2(cat(mm3(mm1(inputs)),mm4(mm1(inputs))))

which first computes `o1=mm1(inputs)` then computes `o3=mm3(inputs+o1)` and
`o4=mm4(inputs+o1)` finally  computes `mm2(inputs+o1+o3+o4)`.

Composing `MeasurementModel`s is where the magic occurs. For a complex measurement model,
it may not be feasible/possible to calculate an analytical form of the partial derivative.
However, if the model can be broken into a series of smaller composable steps for which the
analytical partial derivative can be calculated, the steps can be combined to build the
full model.  Of course this isn't magic, it results from this property of Jacobians:

``\\mathbf{J}_{g(f(x))} |_x = \\mathbf{J}_g(y) |_{f(x)} \\mathbf{J}_{f(x)} |_x``

To use this framework on your measurement model, you must implement the function
`compute(...)` for `YourMeasurementModel`.

    compute(mm::YourMeasurementModel, inputs::LabeledValues, withJac::Boolean)::MMResult

where

    MMResult = Tuple{LabeledValues, Union{Missing,AbstractMatrix{Float64}}}

The function `compute(...)` is responsible for computing the models output values and the
Jacobian of the output values with respect to the input values.  The Jacobian is a matrix
with one row per output value and one column per input value. The contents of the r,c-th
element is the partial derivative of the r-th output value with respect to the c-th input
value.

`MMResult` represents the function output values (`LabeledValues`) and the Jacobian as
an `AbstractMatrix{Float64}`.  To optimize the `compute(...)` function when used to
simply compute the vector valued function, the Jacobian must calculated only when
`withJac=true`.  If `withJac=false` then returning 'missing' for the Jacobian is
encouraged.  This facilates efficiency when using 'compute(...)' in situations like
Monte Carlo propagation (using `mcpropagate(...)`) where it is wasteful to compute
the Jacobian but you'd otherwise like to use identical code to compute the function values.
"""
abstract type MeasurementModel end

"""
    MMResult = Tuple{LabeledValues,Union{Missing,AbstractMatrix{Float64}}}

THe type of return values from the `compute(...)` function.
"""
const MMResult = Tuple{LabeledValues,Union{Missing,AbstractMatrix{Float64}}}

# Allows measurement models to be used like functions
function (mm::MeasurementModel)(inputs::UncertainValues)::UncertainValues
    return propagate(mm, inputs)
end

# Allows measurement models to be used like functions
(mm::MeasurementModel)(inputs::Dict{<:Label,UncertainValue})::UncertainValues =
    propagate(mm, uvs(inputs))
(mm::MeasurementModel)(inputs::LabeledValues)::LabeledValues = compute(mm, inputs, false)[1]

"""
    propagate(mm::MeasurementModel, uvs::UncertainValues)::UncertainValues

Propagate the input (measured values as UncertainValues) through the MeasurementModel to
produce the output values (as UncertainValues).
"""
function propagate(mm::MeasurementModel, inputs::UncertainValues)::UncertainValues
    (outvals, jac) = compute(mm, LabeledValues(labels(inputs), values(inputs)), true)
    return uvs(labels(outvals), values(outvals), jac * inputs.covariance * transpose(jac))
end

struct MCSampler
    inputs::UncertainValues
    zeros::LabeledValues
    nonzeros::Dict{Label,Int}
    mvnorm::MvNormal
    rng::AbstractRNG

    function MCSampler(inp::UncertainValues, rng::AbstractRNG, eps = 1.0e-10)
        f(lbl) = abs(fractional(inp[lbl])) < eps
        z = filter(f, labels(inp))
        zeros = LabeledValues(z, Float64[value(inp, lbl) for lbl in z])
        nz = filter(lbl -> !(lbl in z), labels(inp))
        nzuvs = extract(inp, nz)
        # CSV.write(joinpath(homedir(),"Desktop\\znuvs.csv"), asa(DataFrame,nzuvs))
        nonzeros = Dict(lbl => indexin(nzuvs, lbl) for lbl in nz)
        mvnorm = MvNormal(nzuvs.values, nzuvs.covariance)
        return new(inp, zeros, nonzeros, mvnorm, rng)
    end
end

function sample(samp::MCSampler)::LabeledValues
    reslbls = copy(labels(samp.inputs))
    resvals = zeros(length(samp.inputs))
    rvals = rand(samp.rng, samp.mvnorm)
    for (i, lbl) in enumerate(reslbls)
        if lbl in keys(samp.zeros)
            resvals[i] = samp.zeros[lbl]
        else
            resvals[i] = rvals[samp.nonzeros[lbl]]
        end
    end
    return LabeledValues(reslbls, resvals)
end

"""
    mcpropagate(mm::MeasurementModel, inputs::UncertainValues, n::Int, parallel::Bool=true, rng::AbstractRNG = Random.GLOBAL_RNG)::UncertainValues

Propagate the `inputs` through the `MeasurementModel` using a MonteCarlo algorithm in which
the inputs are assumed to be represented by a `MvNormal` distribution with covariance
from `inputs`.  Performs 'n' iterations and multi-thread if `parallel=true`.
"""
function mcpropagate(
    mm::MeasurementModel,
    inputs::UncertainValues,
    n::Int;
    parallel = false,
    rng::AbstractRNG = Random.GLOBAL_RNG,
)
    perform(inps) = values(compute(mm, inps, false)[1])
    mcs = MCSampler(inputs, rng, 1.0e-8)
    # mvn = MvNormal(inputs.values, inputs.covariance)
    (outvals, _) = compute(mm, LabeledValues(labels(inputs), values(inputs)), false)
    samples, inlabels = Array{Float64}(undef, (length(outvals), n)), labels(inputs)
    if parallel && (nthreads() > 1)
        @threads for i = 1:n
            samples[:, i] = perform(sample(mcs)) # rand(mvn))
        end
    else
        for i = 1:n
            samples[:, i] = perform(sample(mcs))
        end
    end
    return estimated(labels(outvals), samples)
end

"""
    compute(mm::MeasurementModel, uvs::LabeledValues)::LabeledValues

Calculate the output values for the specified set of input LabeledValues.
"""
compute(mm::MeasurementModel, inputs::LabeledValues)::LabeledValues =
    compute(mm, inputs, false)[1]

"""
    MaintainInputs <: MeasurementModel

Carry over a subset of the input variables to the next step in a calculation.  Typically
used in consort with a `ParallelMeasurementModel` when inputs to this step will also
be required in subsequent steps.
"""
struct MaintainInputs <: MeasurementModel
    labels::Vector{Label}
    all::Bool

    MaintainInputs(labels::AbstractVector{<:Label}, all::Bool = false) = new(labels, all)
    MaintainInputs(ty, labels, all::Bool = false) = new(labelsByType(ty, labels), all)
    MaintainInputs(label::Label, all::Bool = false) = new([label], all)
end

function compute(mm::MaintainInputs, inputs::LabeledValues, withJac::Bool)::MMResult
    if mm.all
        jac = withJac ? Array{Float64}(I, length(inputs), length(inputs)) : missing
        return (inputs, jac)
    else
        vals = [inputs[lbl] for lbl in mm.labels]
        jac = withJac ? zeros(length(mm.labels), length(inputs)) : missing
        if withJac
            for (out, lbl) in enumerate(mm.labels)
                jac[out, indexin(inputs, lbl)] = 1.0
            end
        end
        return (LabeledValues(mm.labels, vals), jac)
    end
end

"""
    AllInputs <: MeasurementModel

Carry over all of the input variables to the next step in a calculation.  Typically
used in consort with a `ParallelMeasurementModel` when inputs to this step will also
be required in subsequent steps.
"""
struct AllInputs <: MeasurementModel end

function compute(ai::AllInputs, inputs::LabeledValues, withJac::Bool)::MMResult
    jac = withJac ? Array{Float64}(I, length(inputs), length(inputs)) : missing
    return (inputs, jac)
end

"""
    ParallelMeasurementModel <: MeasurementModel

A `ParallelMeasurementModel` is a collection of `MeasurementModel`s that can be executed
simultaneously since they share the same inputs (or a subset of the same inputs).
`ParallelMeasurementModel`s should be favored over `ComposedMeasurementModel`s since
they can be threaded and combining Jacobians is more computationally efficient.

    ParallelMeasurementModel(models::AbstractVector{MeasurementModel}, multithread=false)

While `ParallelMeasurementModel` supports multi-threading, multi-threading should be used
only when the cost of the calculation is going to exceed the overhead necessary to create
and manage the thread.  Usually this means, only use `multithread=true` on one of the
outer-most steps of a large calculation where splitting the calculation can keep the
Jacobians as small as possible until the last possible moment.

The result of the ParallelMeasurementModel is the union of the outputs from each model.
"""
struct ParallelMeasurementModel <: MeasurementModel
    models::Vector{MeasurementModel}
    multi::Bool

    ParallelMeasurementModel(
        models::AbstractVector{<:MeasurementModel},
        multithread = false,
    ) = new(models, multithread)
end

function compute(
    mm::ParallelMeasurementModel,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    results = Vector{MMResult}(undef, length(mm.models))
    if mm.multi && (nthreads() > 1)
        @threads for i in eachindex(mm.models)
            results[i] = compute(mm.models[i], inputs, withJac)
        end
    else
        for i in eachindex(mm.models)
            results[i] = compute(mm.models[i], inputs, withJac)
        end
    end
    outvals = mapreduce(mmr -> mmr[1], merge, results)
    jac = withJac ? vcat(map(mmr -> mmr[2], results)...) : missing
    return (outvals, jac)
end

"""
    ComposedMeasurementModel <: MeasurementModel

A `ComposedMeasurementModel` is used when the output from step-i will be used in a
subsequent step.  Favor `ParallelMeasurementModel` whenever a collection of
`MeasurementModel`s share the same inputs (or subsets of the same inputs) since
the `ParallelMeasurementModel`s can be run on multiple threads and the Jacobians
can be concatenated rather than multiplied. `ParallelMeasurementModel`s and
`ComposedMeasurementModel`s can be combined to produce efficient calculation models.

The final result of the ComposedMeasurementModel is the output of the final step.
"""
struct ComposedMeasurementModel <: MeasurementModel
    models::Vector{MeasurementModel}
end

function compute(
    smm::ComposedMeasurementModel,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # N inputs
    nextinp, nextjac = inputs, withJac ? I : missing
    for model in smm.models
        (nextinp, jac) = compute(model, nextinp, withJac)
        nextjac = withJac ? jac * nextjac : missing
    end
    return (nextinp, nextjac)
end

"""
    (∘)(mm1::MeasurementModel, mm2::MeasurementModel)

Implements composition of `MeasurementModel`s.

Examples:

    (g ∘ f)(x) == propagate(ComposedMeasurementModel([f, g]), x)
    (g ∘ f)(x) == g(f(x))
    (h ∘ g ∘ f)(x) == propagate(ComposedMeasurementModel([f, g, h]), x)
    (h ∘ g ∘ f)(x) == h(g(f(x)))

Note:

    (h ∘ (g ∘ f))(x) == ((h ∘ g) ∘ f)(x) == (h ∘ g ∘ f)(x)
"""
Base.:∘(mm1::ComposedMeasurementModel, mm2::MeasurementModel) =
    ComposedMeasurementModel([mm2, mm1.models...])

Base.:∘(mm1::MeasurementModel, mm2::ComposedMeasurementModel) =
    ComposedMeasurementModel([mm2.models..., mm1])

Base.:∘(mm1::ComposedMeasurementModel, mm2::ComposedMeasurementModel) =
    ComposedMeasurementModel([mm2.models..., mm1.models...])

Base.:∘(mm1::MeasurementModel, mm2::MeasurementModel) = ComposedMeasurementModel([mm2, mm1])

# So missing measurement models compile down to a NoOP
Base.:∘(mm1::MeasurementModel, mm2::Missing) = mm1
Base.:∘(mm1::Missing, mm2::MeasurementModel) = mm2

"""
    (|)(mm1::MeasurementModel, mm2::MeasurementModel)
    (^)(mm1::MeasurementModel, mm2::MeasurementModel)

Implements a mechanism to combine `MeasurementModel`s that work on the same input to
produce output that is the combination of the outputs of all the measurement models.  This
is useful when a calculation forks into 2 or more distinct calculations which are later
composed as is shown in the examples.

Examples:

    j = f | g # Creates a ParallelMeasurementModel([f,g], false)
    jp = f ^ g # Creates a ParallelMeasurementModel([f,g], true)
    y = j(x) # where y combines the outputs of f(x) and h(x)
    z = (f | g | h)(x) # Conceptually like combine(f(x), g(x), h(x)) into a single output
    zp = (f ^ g ^ h)(x) # ParallelMeasurementModel([f,g,h], true)  Conceptually like combine(f(x), g(x), h(x)) into a single output
    (k ∘ (f | g | h))(x) == k(z) # Conceptually like k(f(x),g(x),h(x))
"""
Base.:|(mm1::ParallelMeasurementModel, mm2::MeasurementModel) =
    ParallelMeasurementModel([mm1.models..., mm2])

Base.:|(mm1::MeasurementModel, mm2::ParallelMeasurementModel) =
    ParallelMeasurementModel([mm1, mm2.models...])

Base.:|(mm1::ParallelMeasurementModel, mm2::ParallelMeasurementModel) =
    ParallelMeasurementModel([mm1.models..., mm2.models...])

Base.:|(mm1::MeasurementModel, mm2::MeasurementModel) = ParallelMeasurementModel([mm1, mm2])

# So missing measurement models compile down to a NoOP
Base.:|(mm1::MeasurementModel, ::Missing) = mm1
Base.:|(::Missing, mm2::MeasurementModel) = mm2

Base.:^(mm1::ParallelMeasurementModel, mm2::MeasurementModel) =
    ParallelMeasurementModel([mm1.models..., mm2], true)

Base.:^(mm1::MeasurementModel, mm2::ParallelMeasurementModel) =
    ParallelMeasurementModel([mm1, mm2.models...], true)

Base.:^(mm1::ParallelMeasurementModel, mm2::ParallelMeasurementModel) =
    ParallelMeasurementModel([mm1.models..., mm2.models...], true)

Base.:^(mm1::MeasurementModel, mm2::MeasurementModel) =
    ParallelMeasurementModel([mm1, mm2], true)

# So missing measurement models compile down to a NoOP
Base.:^(mm1::MeasurementModel, ::Missing) = mm1
Base.:^(::Missing, mm2::MeasurementModel) = mm2

"""
    parallel(mm1::MeasurementModel, models::MeasurementModel...)

Implements a mechanism to combine `MeasurementModel`s that work on the same input to
produce output that is the combination of the outputs of all the measurement models.  This
is useful when a calculation forks into 2 or more distinct calculations which are later
composed as is shown in the examples.

Examples:

    j = parallel(f,g) # Creates a ParallelMeasurementModel([f,g], true)
    y = j(x) # where y combines the outputs of f(x) and h(x)
    z = parallel(f, g, h)(x) # Conceptually like combine(f(x), g(x), h(x)) into a single output
    (k ∘ parallel(f, g, h))(x) == k(z) # Conceptually like k(f(x),g(x),h(x))
"""
parallel(mm1::MeasurementModel, models::MeasurementModel...) =
    ParallelMeasurementModel([mm1, models...])

"""
    filter(labels::AbstractVector{<:Label}, mmr::MMResult)::MMResult

Trims an MMResult down to only those labels in `labels`.  This can optimize calculations
which are carrying a lot of intermediary results that are no longer needed.  All Jacobian
input columns are maintained but only those output rows associated with a label in `labels`
will be retained.  This method reorders the output values to match the order in `labels`.

If the input 'mmr' is N functions of M variables and `length(labels)=P` then the result
will be P functions of M variables.
"""
function Base.filter(labels::AbstractVector{<:Label}, mmr::MMResult)::MMResult
    indexes = map(lbl -> indexin(mmr[1], lbl), labels)
    return (
        LabeledValues(labels(mmr[1])[indexes], values(mmr[1])[indexes]),
        mmr[2][indexes, :],
    )
end


text:
export SolverLevMarq

struct SolverLevMarq{T<:Real} <: AbstractNextBasisSolver
  maxiter::Integer
  xtol::T
  gtol::T
  ftol::T
  λinit::T
  λinc::T
  λdec::T
  min_step_quality::T
  good_step_quality::T
  λmax::T
  λmin::T
  diagmin::T
end

function SolverLevMarq(;
  maxiter::Integer=100,
  xtol::T=1e-8,
  gtol::T=1e-8,
  ftol::T=1e-5,
  λinit::T=1.0,
  λinc::T=10.,
  λdec::T=10.,
  min_step_quality::T=1e-3,
  good_step_quality::T=.75,
  λmax::T = 1e12,
  λmin::T = 1e-12,
  diagmin::T = 1e-6
  ) where T
  @assert 0 ≤ min_step_quality < good_step_quality ≤ 1
  SolverLevMarq{T}(maxiter, xtol, gtol, ftol,
                λinit, λinc, λdec,
                min_step_quality, good_step_quality,
                λmax, λmin, diagmin)
end

function (self::SolverLevMarq{T})(prob::Problem{L, B, CS},
                                  r::Vector{T},
                                  uinit::Vector{T};
                                  useprojections::Bool=true,
                                  savetrace::Bool=false) where {T, L, B, CS}
  u = copy(uinit)
  λ = self.λinit
  xconverged, gconverged, λconverged, fconverged = false, false, false, false
  Δu = zeros(u)
  tr = DataFrame(i=Int[], sse=T[], prsse=T[], trsse=T[], ρ=T[],
                 λ=Int[], lo=String[], ri=String[],
                 normΔu=T[], normg=T[], expdecr=T[]
                 )  # TODO: Time consuming
  local F::Vector{T}, J::Matrix{T}

  for iter in 1:self.maxiter
    F, J = valueandjacobian(prob, u)          # TODO: Time consuming
    F -= r
    sse = sum(abs2, F)
    JJ = J'*J
    D = Diagonal(JJ)                # TODO: Time consuming
    D[0 .< D .< self.diagmin] = self.diagmin  # TODO: -do-
    Jf = -J'*F                      # TODO: transp is Time consuming
                                    # This is also the neg gradient
    normg = norm(Jf, Inf)
    gconverged = normg < self.gtol
    gconverged && break

    accept_step = false
    while !(accept_step || λconverged || fconverged)
      JJD = JJ + λ*D                # TODO: Time consuming
      Δu = JJD\Jf
      if useprojections
        Δu = projectΩ(prob, u+Δu)-u
      end

      # Try this step and compare with linear prediction
      pred_sse = sum(abs2, J*Δu+F)
      trialf = predictor(prob, u+Δu) - r    # TODO: Time consuming   #X
      trial_sse = sum(abs2, trialf)
      expd_decr = 1 - pred_sse/sse          # Positive
      fconverged = expd_decr < self.ftol

      ρ = (sse-trial_sse)/(sse-pred_sse)
      lower_trust = ρ ≤  self.min_step_quality
      raise_trust = ρ > self.good_step_quality

      normΔu = norm(Δu)
      logrelλ = round(Int, log(self.λinc, λ))
      trentry = [iter, sse, pred_sse, trial_sse, ρ,
                 logrelλ,
                 lower_trust ? "↑" : " ",
                 raise_trust ? "↓" : " ",
                 normΔu, normg, expd_decr]
      push!(tr, trentry)

      if lower_trust
        if λ==self.λmax
          λconverged = true
        else
          λ = min(λ*self.λinc, self.λmax)
        end
      else        #Accept this step
        accept_step = true
        u += Δu
        if raise_trust
          λ = max(λ/self.λdec, self.λmin)
        end
      end
    end   #Inner Loop over λ

    xconverged = norm(Δu) < self.xtol
    (xconverged || gconverged || λconverged || fconverged) && break
  end    #Outer Loop over Δu

  savetrace && print(tr)

  u, sum(abs2, predictor(prob, u)-r)  #X
end


text:
struct GlobalsContext
    invokes::Set{Any}
end
GlobalsContext() = GlobalsContext(Set())

"""
    fix_globals!(mod::LLVM.Module)

Replace function addresses in `mod` with references to global data structures.
For each global variable, two LLVM global objects are created:

* `jl.global.data` -- An LLVM 'i8' vector holding a serialized version of the Julia object.
* `jl.global` -- A pointer to the unserialized Julia object.

The `inttopt` with the function address is replaced by `jl.global`.

A function `jl_init_globals` is added to `mod`. This function deserializes the data in
`jl.global.data` and updates `jl.global`.
"""

_opcode(x::LLVM.ConstantExpr) = LLVM.API.LLVMGetConstOpcode(LLVM.ref(x))

function fix_globals!(mod::LLVM.Module)
    # Create a `jl_init_globals` function.
    jl_init_globals_func = LLVM.Function(mod, "jl_init_globals",
                                         LLVM.FunctionType(julia_to_llvm(Cvoid), LLVMType[]))
    jl_init_global_entry = BasicBlock(jl_init_globals_func, "entry", context(mod))

    # Definitions for utility functions
    func_type = LLVM.FunctionType(julia_to_llvm(Any), LLVMType[LLVM.PointerType(julia_to_llvm(Int8))])
    deserialize_funs = Dict()

    uint8_t = julia_to_llvm(UInt8)

    ctx = SerializeContext()
    es = []
    objs = Set()
    gptridx = Dict()
    instrs = []
    gptrs = []
    j = 1   # counter for position in gptridx
    Builder(context(mod)) do builder
        toinstr!(x) = x
        function toinstr!(x::LLVM.ConstantExpr)
            if _opcode(x) == LLVM.API.LLVMAddrSpaceCast
                val = toinstr!(first(operands(x)))
                ret = addrspacecast!(builder, val, llvmtype(x))
                return ret
            elseif _opcode(x) == LLVM.API.LLVMGetElementPtr
                ops = operands(x)
                val = toinstr!(first(ops))
                ret = gep!(builder, val, [ops[i] for i in 2:length(ops)])
                return ret
            elseif _opcode(x) == LLVM.API.LLVMBitCast
                ops = operands(x)
                val = toinstr!(first(ops))
                ret = pointercast!(builder, val, llvmtype(x))
                return ret
            elseif _opcode(x) == LLVM.API.LLVMIntToPtr
                ptr = Ptr{Any}(convert(Int, first(operands(x))))
                obj = unsafe_pointer_to_objref(ptr)
                if !in(obj, objs)
                    push!(es, serialize(ctx, obj))
                    push!(objs, obj)
                    # Create pointers to the data.
                    gptr = GlobalVariable(mod, julia_to_llvm(Any), "jl.global")
                    linkage!(gptr, LLVM.API.LLVMInternalLinkage)
                    LLVM.API.LLVMSetInitializer(LLVM.ref(gptr), LLVM.ref(null(julia_to_llvm(Any))))
                    push!(gptrs, gptr)
                    gptridx[obj] = j
                    j += 1
                end
                gptr = gptrs[gptridx[obj]]
                gptr2 = load!(builder, gptr)
                ret = pointercast!(builder, gptr2, llvmtype(x))
                return ret
            end
            return x
        end
        for fun in functions(mod)
            if startswith(LLVM.name(fun), "jfptr")
                unsafe_delete!(mod, fun)
                continue
            end

            for blk in blocks(fun), instr in instructions(blk)
                # Set up functions to walk the operands of the instruction
                # and convert appropriate ConstantExpr's to instructions.
                # Look for `LLVMIntToPtr` expressions.
                position!(builder, instr)
                ops = operands(instr)
                N = opcode(instr) == LLVM.API.LLVMCall ? length(ops) - 1 : length(ops)
                if opcode(instr) == LLVM.API.LLVMCall && name(last(operands(instr))) == "jl_type_error"
                    continue
                end
                for i in 1:N
                    try
                        if opcode(instr) == LLVM.API.LLVMPHI
                            position!(builder, last(instructions(LLVM.incoming(instr)[i][2])))
                        end
                        ops[i] = toinstr!(ops[i])
                    catch x
                    end
                end
            end
        end
    end
    nglobals = length(es)
    #@show mod
    #verify(mod)
    for i in 1:nglobals
        # Assign the appropriate function argument to the appropriate global.
        es[i] = :(unsafe_store!($((Symbol("global", i))), $(es[i])))
    end
    # Define the deserializing function.
    fune = quote
        function _deserialize_globals(Vptr, $((Symbol("global", i) for i in 1:nglobals)...))
            $(ctx.init...)
            $(es...)
            return
        end
    end
    # @show fune
    # Execute the deserializing function.
    deser_fun = eval(fune)
    v = take!(ctx.io)
    gv_typ = LLVM.ArrayType(uint8_t, length(v))
    data = LLVM.GlobalVariable(mod, gv_typ, "jl.global.data")
    linkage!(data, LLVM.API.LLVMExternalLinkage)
    constant!(data, true)
    LLVM.API.LLVMSetInitializer(LLVM.ref(data),
                                LLVM.API.LLVMConstArray(LLVM.ref(uint8_t),
                                                        [LLVM.ref(ConstantInt(uint8_t, x)) for x in v],
                                                        UInt32(length(v))))
    Builder(context(mod)) do builder
        dataptr = gep!(builder, data, [ConstantInt(0, context(mod)), ConstantInt(0, context(mod))])

        # Create the Julia object from `data` and include that in `init_fun`.
        position!(builder, jl_init_global_entry)
        gfunc_type = LLVM.FunctionType(julia_to_llvm(Cvoid),
                                       LLVMType[LLVM.PointerType(julia_to_llvm(Int8)),
                                                Iterators.repeated(LLVM.FunctionType(julia_to_llvm(Any)), nglobals)...])
        deserialize_globals_func = LLVM.Function(mod, "_deserialize_globals", gfunc_type)
        LLVM.linkage!(deserialize_globals_func, LLVM.API.LLVMExternalLinkage)
        for i in 1:nglobals
            # The following fix is to match the argument types which are an integer, not a %jl_value_t**.
            gptrs[i] = LLVM.ptrtoint!(builder, gptrs[i], julia_to_llvm(Csize_t))
        end
        LLVM.call!(builder, deserialize_globals_func, LLVM.Value[dataptr, gptrs...])
        ret!(builder)
    end
    tt = Tuple{Ptr{UInt8}, Iterators.repeated(Ptr{Any}, nglobals)...}
    deser_mod = irgen(deser_fun, tt, overdub = false, fix_globals = false, optimize_llvm = false)
    d = find_ccalls(deser_fun, tt)
    fix_ccalls!(deser_mod, d)
    # rename deserialization function to "_deserialize_globals"
    fun = first(TypedCodeUtils.filter(x -> LLVM.name(x) == "_deserialize_globals", functions(deser_mod)))[2]
    # LLVM.name!(fun, "_deserialize_globals")
    linkage!(fun, LLVM.API.LLVMExternalLinkage)
    # link into the main module
    LLVM.link!(mod, deser_mod)
    return
end


text:
# Autogenerated wrapper script for NativeFileDialog_jll for x86_64-linux-gnu
export libnfd

using GTK3_jll
JLLWrappers.@generate_wrapper_header("NativeFileDialog")
JLLWrappers.@declare_library_product(libnfd, "libnfd.so")
function __init__()
    JLLWrappers.@generate_init_header(GTK3_jll)
    JLLWrappers.@init_library_product(
        libnfd,
        "lib/libnfd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Plots
using CSV
using DataFrames
using ArgParse

function parse_commandline()

	s = ArgParseSettings()

	@add_arg_table s begin
		"--function"
		required = true
		"--file"
		required = true
		"--dpi"
		arg_type = Int
		default = 600
	end

	return parse_args(s)
end

function get_fitness_function(args)

	if cmp(args["function"], "rastrigin") == 0
		fun = x -> 10length(x) + sum(x.^2 .- 10 .* cos.(2π .* x))
	elseif cmp(args["function"], "spheric") == 0
		fun = x -> sum(x.^2)
	elseif cmp(args["function"], "rosenbrock") == 0
		fun = x -> (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
	end

	return fun
end

function main()
	args = parse_commandline()
	x = y = range(-0.5, stop=0.5, length=200)

	fun = get_fitness_function(args)

	l = Vector{Float64}()
	for i in x
	   for j in y
	       v = [i j]
	       append!(l, fun(v))
	   end
	end

	csv_file = CSV.read(args["file"], DataFrame)
	filename = split(args["file"], "/")[2]

	contour(x,y,l, dpi=args["dpi"])
	for (i, idx) in enumerate(csv_file.hive)
		if idx == -1
			label="WoC"
		else
			label="Hive $idx"
		end
		scatter!([csv_file.x[i]], [csv_file.y[i]], label=label)
	end
	png("./"*args["function"]*"/"*filename*".png")

end

main()