text:
using IICGP
using Test

cfg_filename = string(@__DIR__, "/dualcgpga_test.yaml")
game = "gravitar"
mcfg, ecfg, ccfg, reducer, bootstrap = IICGP.dualcgp_config(
    cfg_filename, game
)
logid = "2021-test-logid"
resdir = dirname(@__DIR__)
fit(e::SymInd, c::SymInd) = e.chromosome[1] * c.chromosome[1]

@testset "CGP GA Evolution" begin
    evo = IICGP.DualCGPGAEvo(mcfg, ecfg, ccfg, fit, logid, resdir)
    @test length(evo.encoder_sympop) == ecfg.n_population
    @test length(evo.controller_sympop) == ccfg.n_population
    sympops = [evo.encoder_sympop, evo.controller_sympop]
    cfgs = [ecfg, ccfg]
    types = ["encoder", "controller"]
    for l in eachindex(sympops)
        sympop = sympops[l]
        cfg = cfgs[l]
        type = types[l]
        n_chr = cfg.rows * cfg.columns * (3 + cfg.n_parameters) + cfg.n_out
        for i in eachindex(sympop)
            @test typeof(sympop[i].chromosome) <: Vector{Float64}
            @test length(sympop[i].chromosome) == n_chr
            @test sympop[i].index == i
            @test sympop[i].fitness == [-Inf]
            @test sympop[i].type == type
        end
    end
    mat_size = (ecfg.n_population, ccfg.n_population)
    @test size(evo.fitness_matrix) == mat_size
    @test size(evo.elites_matrix) == mat_size
    @test all(i -> i == -Inf, evo.fitness_matrix)
    @test all(i -> i == false, evo.elites_matrix)
end


text:
module BugReporting

export replay, make_interactive_report

using rr_jll
using Zstd_jll
using HTTP, JSON
using AWSCore, AWSS3
using Tar
using Pkg

const WSS_ENDPOINT = "wss://53ly7yebjg.execute-api.us-east-1.amazonaws.com/test"
const GITHUB_APP_ID = "Iv1.c29a629771fe63c4"
const TRACE_BUCKET = "julialang-dumps"

function check_rr_available()
    if isempty(rr_jll.rr_path)
        error("RR not available on this platform")
    end
end

function collect_inner_traces(trace_directory)
    # If this is already a trace directory, return just it
    if isfile(joinpath(trace_directory, "version"))
        return [trace_directory]
    end
    # If this is a directory containing traces, return those
    traces = String[]
    for f in readdir(trace_directory, join=true)
        if isfile(joinpath(f, "version"))
            push!(traces, f)
        end
    end
    return unique(realpath.(traces))
end

function find_latest_trace(trace_directory)
    # What _RR_TRACE_DIR calls a "trace directory" is not what `rr pack` calls a "trace directory"
    # This function allows us to normalize to the inner "latest-trace" directory, if necessary.
    latest_symlink = joinpath(trace_directory, "latest-trace")
    if !isfile(joinpath(trace_directory, "version")) && islink(latest_symlink)
        return realpath(latest_symlink)
    end
    return trace_directory
end

function rr_pack(trace_directory)
    check_rr_available()

    rr() do rr_path
        for dir in collect_inner_traces(trace_directory)
            @debug("rr pack'ing $(dir)")
            run(`$rr_path pack $(dir)`)
        end
    end
end

function rr_record(args...; trace_dir=nothing)
    check_rr_available()

    record_flags = split(get(ENV, "JULIA_RR_RECORD_ARGS", ""), ' ', keepempty=false)
    rr() do rr_path
        new_env = copy(ENV)
        if trace_dir !== nothing
            new_env["_RR_TRACE_DIR"] = trace_dir
        end
        # Intersperse all given arguments with spaces, then splat:
        rr_cmd = `$(rr_path) record $(record_flags)`
        for arg in args
            rr_cmd = `$(rr_cmd) $(arg)`
        end
        run(ignorestatus(setenv(rr_cmd, new_env)))
    end
end

function download_rr_trace(trace_url; verbose=true)
    Pkg.PlatformEngines.probe_platform_engines!()
    artifact_hash = Pkg.create_artifact() do dir
        mktempdir() do dl_dir
            # Download into temporary directory, unpack into artifact directory
            local_path = joinpath(dl_dir, "trace.tar.zst")
            Pkg.PlatformEngines.download(trace_url, local_path; verbose=verbose)
            proc = zstdmt() do zstdp
                open(`$zstdp --stdout -d $local_path`, "r+")
            end
            Tar.extract(proc, dir)
        end
    end
    return Pkg.artifact_path(artifact_hash)
end

function replay(trace_url)
    if startswith(trace_url, "s3://")
        trace_url = string("https://s3.amazonaws.com/julialang-dumps/", trace_url[6:end])
    end
    if startswith(trace_url, "https://")
        trace_url = download_rr_trace(trace_url)
    end

    if !isdir(trace_url)
        error("Invalid trace location: $(trace_url)")
    end

    rr() do rr_path
        run(`$(rr_path) replay $(find_latest_trace(trace_url))`)
    end
end

function make_interactive_report(report_type, ARGS=[])
    default_julia_args = `--history-file=no`
    if report_type == "justrr"
        rr_record(`$(Base.julia_cmd()) $default_julia_args`, ARGS)
        return
    elseif report_type == "rr"
        artifact_hash = Pkg.create_artifact() do trace_dir
            rr_record(`$(Base.julia_cmd()) $default_julia_args`, ARGS; trace_dir=trace_dir)
            @info "Preparing trace directory for upload (if your trace is large this may take a few minutes)"
            rr_pack(trace_dir)
        end
        upload_rr_trace(Pkg.artifact_path(artifact_hash))
        return
    end
    error("Unknown report type")
end

const S3_CHUNK_SIZE = 25 * 1024 * 1024 # 25 MB

include("sync_compat.jl")

function upload_rr_trace(trace_directory)
    # Auto-pack this trace directory if it hasn't already been:
    rr_pack(trace_directory)

    c = Channel()
    t = @async HTTP.WebSockets.open(WSS_ENDPOINT) do ws
        write(ws, "Hello Server, if it's not too much trouble, please send me S3 credentials")
        x = readavailable(ws)
        put!(c, JSON.parse(String(x))["connectionId"])
        # This will block until the user has completed the authentication flow
        x = readavailable(ws)
        push!(c, JSON.parse(String(x)))
    end
    bind(c, t)
    connectionId = take!(c)
    println()
    println("""
    ### IMPORTANT =============================================================
    You are about to upload a trace directory to a publicly accessible location.
    Such traces contain any information that was accessed by the traced
    executable during its execution. This includes any code loaded, any
    secrets entered, the contents of any configuration files, etc.
    DO NOT proceed, if you do not wish to make this information publicly available.
    By proceeding you explicitly agree to waive any privacy interest in the
    uploaded information.
    ### =======================================================================
    """)
    println("To upload a trace, please authenticate, by visiting:\n")
    println("\thttps://github.com/login/oauth/authorize?client_id=$GITHUB_APP_ID&state=$(HTTP.escapeuri(connectionId))")
    s3creds = take!(c)

    println()
    @info "Uploading Trace directory"

    creds = AWSCore.AWSCredentials(
        s3creds["AWS_ACCESS_KEY_ID"],
        s3creds["AWS_SECRET_ACCESS_KEY"],
        s3creds["AWS_SESSION_TOKEN"])
    aws = AWSCore.aws_config(creds = creds, region="us-east-1")

    # Tar it up
    proc = zstdmt() do zstdp
        open(`$zstdp -`, "r+")
    end

    t = @async begin try
        upload = s3_begin_multipart_upload(aws, TRACE_BUCKET, s3creds["UPLOAD_PATH"])
        tags = Vector{String}()
        i = 1
        @Base.Experimental.sync begin
            while isopen(proc)
                buf = Vector{UInt8}(undef, S3_CHUNK_SIZE)
                n = readbytes!(proc, buf)
                n < S3_CHUNK_SIZE && resize!(buf, n)
                let partno = i, buf=buf
                    @async begin
                        try
                            push!(tags, s3_upload_part(aws, upload, partno, buf))
                        catch e
                            close(proc)
                            rethrow(e)
                        end
                    end
                end
                i += 1
            end
        end
        s3_complete_multipart_upload(aws, upload, tags)
    catch e
        Base.showerror(stderr, e)
    end
    end

    # Start the Tar creation process, the file will be uploaded as it's created
    Tar.create(trace_directory, proc)
    close(proc.in)

    wait(t)
    println("Uploaded to https://s3.amazonaws.com/$TRACE_BUCKET/$(s3creds["UPLOAD_PATH"])")
end

end # module

ext: jl


text:
module QuasiTriangular
## QuasiUpperTriangular Matrix of real numbers

using LinearAlgebra
import Base.size
import Base.similar
import Base.getindex
import Base.require_one_based_indexing
import Base.setindex!
import Base.copy
import LinearAlgebra.checksquare
import LinearAlgebra.BlasInt
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra.BLAS.libblas

export QuasiUpperTriangular, I_plus_rA_ldiv_B!, I_plus_rA_plus_sB_ldiv_C!, A_ldiv_B!, A_rdiv_B!, A_rdiv_Bt!, A_mul_B!, A_mul_Bt!, At_mul_B!

abstract type AbstractQuasiTriangular{T, S <: AbstractMatrix} <: AbstractMatrix{T} end

struct QuasiUpperTriangular{T, S <: AbstractMatrix{T}} <: AbstractQuasiTriangular{T,S}
    data::S

    function QuasiUpperTriangular{T,S}(data) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(data)
        checksquare(data)
        new(data)
    end
end
QuasiUpperTriangular(A::QuasiUpperTriangular) = A
QuasiUpperTriangular{T}(A::QuasiUpperTriangular{T}) where {T} = A
function QuasiUpperTriangular(A::AbstractMatrix)
    return QuasiUpperTriangular{eltype(A), typeof(A)}(A)
end
function QuasiUpperTriangular{T}(A::AbstractMatrix) where T
    QuasiUpperTriangular(convert(AbstractMatrix{T}, A))
end

function QuasiUpperTriangular{T}(A::QuasiUpperTriangular) where T
    Anew = convert(AbstractMatrix{T}, A.data)
    QuasiUpperTriangular(Anew)
end
Matrix(A::QuasiUpperTriangular{T}) where {T} = Matrix{T}(A)

size(A::QuasiUpperTriangular, d) = size(A.data, d)
size(A::QuasiUpperTriangular) = size(A.data)

# For A<:AbstractTriangular, similar(A[, neweltype]) should yield a matrix with the same
# triangular type and underlying storage type as A. The following method covers these cases.
similar(A::QuasiUpperTriangular, ::Type{T}) where {T} = QuasiUpperTriangular(similar(parent(A), T))
# On the other hand, similar(A, [neweltype,] shape...) should yield a matrix of the underlying
# storage type of A (not wrapped in a triangular type). The following method covers these cases.
similar(A::QuasiUpperTriangular, ::Type{T}, dims::Dims{N}) where {T,N} = similar(parent(A), T, dims)

Array(A::AbstractQuasiTriangular) = Matrix(A)
parent(A:: AbstractQuasiTriangular) = A.data

copy(A::QuasiUpperTriangular) = QuasiUpperTriangular(copy(A.data))

real(A::QuasiUpperTriangular{<:Real}) = A
real(A::QuasiUpperTriangular{<:Complex}) = (B = real(A.data); QuasiUpperTriangular(B))

getindex(A::QuasiUpperTriangular{T,S}, i::Integer, j::Integer) where {T,S<:AbstractMatrix{T}} =
    i <= j + 1 ? A.data[i,j] : zero(A.data[j,i])

function setindex!(A::QuasiUpperTriangular, x, i::Integer, j::Integer)
    if i > j + 1
        x == 0 || throw(ArgumentError("cannot set index in the lower triangular part " *
            "($i, $j) of an UpperTriangular matrix to a nonzero value ($x)"))
    else
        A.data[i,j] = x
    end
    return A
end

## Generic quasi triangular right vector, multiplication
function A_mul_B!(a::QuasiUpperTriangular,b::AbstractVector,work::AbstractVector)
    if size(a,1) < 27
        # pure Julia is faster
        A_mul_B!(a,b)
        return
    end
    copy!(work,b)
    BLAS.trmv!('U','N','N',a.data,b)
    @inbounds @simd for i= 2:size(a,1)
        b[i] += a[i,i-1]*work[i-1]
    end
end

# pure Julia implementation right multiplication
function A_mul_B!(A::QuasiUpperTriangular, B::AbstractVector)
    m = length(B)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    @inbounds Bi2 = A.data[1,1]*B[1]
    @inbounds @simd for k = 2:m
        Bi2 += A.data[1,k]*B[k]
    end
    @inbounds for i = 2:m
        Bi1 = A.data[i,i-1]*B[i-1]
        @simd for k = i:m
            Bi1 += A.data[i,k]*B[k]
        end
        B[i-1] = Bi2
        Bi2 = Bi1
    end
    B[m] = Bi2
end

# pure Julia transpose implementation right vector multiplication
function At_mul_B!(A::QuasiUpperTriangular, B::AbstractVector)
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    @inbounds for j = 1:n
        Bij2 = A.data[m,m]'B[m,j]
        @simd for k = 1:m - 1
            Bij2 += A.data[k,m]'B[k,j]
        end
        for i = m-1:-1:1
            Bij1 = A.data[i+1,i]'B[i+1,j]
            @simd for k = 1:i
                Bij1 += A.data[k,i]'B[k,j]
            end
            B[i+1,j] = Bij2
            Bij2 = Bij1
        end
        B[1,j] = Bij2
    end
    B
end

# Pure Julia implementation left vector multiplication
function A_mul_B!(A::AbstractVector, B::QuasiUpperTriangular)
    m, n = size(A)
    if size(B, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    @inbounds for i = 1:m
        Aij2 = A[i,n]*B.data[n,n]
        @simd for k = 1:n - 1
            Aij2 += A[i,k]*B.data[k,n]
        end
        for j = n-1:-1:1
            Aij1 = A[i,j+1]*B.data[j+1,j]
            @simd for k = 1:j
                Aij1 += A[i,k]*B.data[k,j]
            end
            A[i,j+1] = Aij2
            Aij2 = Aij1
        end
        A[i,1] = Aij2
    end
    A
end

# Pure Julia implementation transpose left vector multiplication
function A_mul_Bt!(A::AbstractVector, B::QuasiUpperTriangular)
    m, n = size(A)
    if size(B, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    @inbounds for i = 1:m
        Aij2 = A[i,1]*B.data[1,1]'
        @simd for k = 2:n
            Aij2 += A[i,k]*B.data[1,k]'
        end
        for j = 2:n
            Aij1 = A[i,j-1]*B.data[j,j-1]'
            @simd for k = j:n
                Aij1 += A[i,k]*B.data[j,k]'
            end
            A[i,j-1] = Aij2
            Aij2 = Aij1
        end
        A[i,n] = Aij2
    end
    A
end

# right matrix multiplication c = a*b
function A_mul_B!(c::AbstractMatrix, a::QuasiUpperTriangular, b::AbstractMatrix)
    A_mul_B!(c, 1.0, a, b)
    c
end

# right matrix multiplication on view in VecOrMat
function A_mul_B!(c::AbstractVecOrMat, a::QuasiUpperTriangular, b::AbstractVecOrMat, nr::Int64, nc::Int64)
    A_mul_B!(c, 1.0, a, b, nr, nc)
    c
end

# c = alpha*a*b
function A_mul_B!(c::AbstractMatrix, alpha::Float64, a::QuasiUpperTriangular, b::AbstractMatrix)
    m, n = size(b)
    if size(a, 1) != m
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(b,1))"))
    end
    copy!(c,b)
    BLAS.trmm!('L','U','N','N',alpha,a.data,c)

    @inbounds for i= 2:m
        x = alpha*a[i,i-1]
        @simd for j=1:n
            c[i,j] += x*b[i-1,j]
        end
    end
    c
end

# c = alpha*a*b in view for VecOrMat
function A_mul_B!(c::AbstractVecOrMat, alpha::Float64, a::QuasiUpperTriangular, b::AbstractVecOrMat, nr::Int64, nc::Int64)
    m, n = size(a)
    copy!(c,b)
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('L'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(nr), Ref{BlasInt}(nc),
          Ref{Float64}(alpha), a.data, Ref{BlasInt}(nr), c, Ref{BlasInt}(nr))

    b1 = reshape(b,nr,nc)
    c1 = reshape(c,nr,nc)
    @inbounds for i= 2:m
        x = a[i,i-1]
        @simd for j=1:m
            c1[i,j] += x*b1[i-1,j]
        end
    end
    c
end

# right matrix multiplication
function A_mul_B!(A::QuasiUpperTriangular, B::AbstractMatrix)
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    @inbounds for j = 1:n
        Bij2 = A.data[1,1]*B[1,j]
        @simd for k = 2:m
            Bij2 += A.data[1,k]*B[k,j]
        end
        for i = 2:m
            Bij1 = A.data[i,i-1]*B[i-1,j]
            @simd for k = i:m
                Bij1 += A.data[i,k]*B[k,j]
            end
            B[i-1,j] = Bij2
            Bij2 = Bij1
        end
        B[m,j] = Bij2
    end
    B
end

function At_mul_B!(c::AbstractMatrix, a::QuasiUpperTriangular, b::AbstractMatrix)
    At_mul_B!(c,1.0,a,b)
    c
end

function At_mul_B!(c::AbstractMatrix, alpha::Float64, a::QuasiUpperTriangular, b::AbstractMatrix)
    m, n = size(b)
    if size(a, 1) != m
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(b,1))"))
    end
    copy!(c,b)
    BLAS.trmm!('L','U','T','N',alpha,a.data,c)

    @inbounds for i= 1:m-1
        x = alpha*a[i+1,i]
        @simd for j=1:n
            c[i,j] += x*b[i+1,j]
        end
    end
    c
end

function At_mul_B!(A::QuasiUpperTriangular, B::AbstractMatrix)
    m, n = size(B, 1), size(B, 2)
    if m != size(A, 1)
        throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m"))
    end
    @inbounds for j = 1:n
        Bij2 = A.data[m,m]'B[m,j]
        @simd for k = 1:m - 1
            Bij2 += A.data[k,m]'B[k,j]
        end
        for i = m-1:-1:1
            Bij1 = A.data[i+1,i]'B[i+1,j]
            @simd for k = 1:i
                Bij1 += A.data[k,i]'B[k,j]
            end
            B[i+1,j] = Bij2
            Bij2 = Bij1
        end
        B[1,j] = Bij2
    end
    B
end

function A_mul_B!(c::AbstractMatrix, a::AbstractMatrix, b::QuasiUpperTriangular)
    A_mul_B!(c,1.0,a,b)
    c
end

function A_mul_B!(c::AbstractVecOrMat, a::AbstractVecOrMat, b::QuasiUpperTriangular, nr::Int64, nc::Int64)
    A_mul_B!(c,1.0,a,b,nr,nc)
    c
end

function A_mul_B!(c::AbstractMatrix, alpha::Float64, a::AbstractMatrix, b::QuasiUpperTriangular)
    m, n = size(a)
    if size(b, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    copy!(c,a)
    BLAS.trmm!('R','U','N','N',alpha,b.data,c)

    @inbounds for i= 2:n
        x = alpha*b[i,i-1]
        @simd for j=1:m
            c[j,i-1] += x*a[j,i]
        end
    end
    c
end

function A_mul_B!(c::AbstractVecOrMat, alpha::Float64, a::AbstractVecOrMat, b::QuasiUpperTriangular, nr::Int64, nc::Int64)
    m, n = size(b)
    copy!(c,a)
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('R'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(nr), Ref{BlasInt}(nc),
          Ref{Float64}(alpha), b.data, Ref{BlasInt}(nc), c, Ref{BlasInt}(nr))

    a1 = reshape(a,nr,nc)
    c1 = reshape(c,nr,nc)
    @inbounds for i= 2:m
        x = b[i,i-1]
        @simd for j=1:m
            c1[j,i-1] += x*a1[j,i]
        end
    end
    c
end

function A_mul_B!(A::AbstractMatrix, B::QuasiUpperTriangular)
    m, n = size(A)
    if size(B, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    @inbounds for i = 1:m
        Aij2 = A[i,n]*B.data[n,n]
        @simd for k = 1:n - 1
            Aij2 += A[i,k]*B.data[k,n]
        end
        for j = n-1:-1:1
            Aij1 = A[i,j+1]*B.data[j+1,j]
            @simd for k = 1:j
                Aij1 += A[i,k]*B.data[k,j]
            end
            A[i,j+1] = Aij2
            Aij2 = Aij1
        end
        A[i,1] = Aij2
    end
    A
end

function A_mul_Bt!(c::AbstractMatrix, a::AbstractMatrix, b::QuasiUpperTriangular)
    A_mul_Bt!(c,1.0,a,b)
    c
end

function A_mul_Bt!(c::AbstractMatrix, alpha::Float64, a::AbstractMatrix, b::QuasiUpperTriangular)
    m, n = size(a)
    if size(b, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    copy!(c,a)
    BLAS.trmm!('R','U','T','N',alpha,b.data,c)
    @inbounds for j= 2:n
        x = alpha*b[j,j-1]
        @simd for i=1:m
            c[i,j] += x*a[i,j-1]
        end
    end
    c
end

function A_mul_Bt!(A::AbstractMatrix, B::QuasiUpperTriangular)
    m, n = size(A)
    if size(B, 1) != n
        throw(DimensionMismatch("right hand side B needs first dimension of size $n, has size $(size(B,1))"))
    end
    @inbounds for i = 1:m
        Aij2 = A[i,1]*B.data[1,1]'
        @simd for k = 2:n
            Aij2 += A[i,k]*B.data[1,k]'
        end
        for j = 2:n
            Aij1 = A[i,j-1]*B.data[j,j-1]'
            @simd for k = j:n
                Aij1 += A[i,k]*B.data[j,k]'
            end
            A[i,j-1] = Aij2
            Aij2 = Aij1
        end
        A[i,n] = Aij2
    end
    A
end

A_mul_B!(c::QuasiUpperTriangular,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = A_mul_B!(c.data,a,b.data)
At_mul_B!(c::QuasiUpperTriangular,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = At_mul_B!(c.data,a,b.data)
A_mul_Bt!(c::QuasiUpperTriangular,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = A_mul_Bt!(c.data,a,b.data)
A_mul_B!(c::AbstractMatrix,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = A_mul_B!(c,a,b.data)
At_mul_B!(c::AbstractMatrix,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = At_mul_B!(c,a,b.data)
A_mul_Bt!(c::AbstractMatrix,a::QuasiUpperTriangular,b::QuasiUpperTriangular) = A_mul_Bt!(c,a,b.data)

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64}, offset_a::Int64,
                  ma::Int64, na::Int64, b::QuasiUpperTriangular, offset_b::Int64, nb::Int64)
    m, n = size(b)
    copyto!(c, offset_c, a, offset_a, ma*na)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('R'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(ma), Ref{BlasInt}(na),
          Ref{Float64}(alpha), b.data, Ref{BlasInt}(na), Ref(c, offset_c), Ref{BlasInt}(ma))

    inda = offset_a + ma
    indc = offset_c
    @inbounds for i = 2:na
        x = b[i,i-1]
        @simd for j = 1:ma
            c[indc] += x*a[inda]
            inda += 1
            indc += 1
        end
    end
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::AbstractVecOrMat{Float64}, offset_a::Int64,
                  ma::Int64, na::Int64, b::QuasiUpperTriangular)
    nb = size(b, 2)
    A_mul_B!(c, offset_c, a, offset_a, ma, na, b, 1, nb)
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::QuasiUpperTriangular, offset_a::Int64,
                  ma::Int64, na::Int64, b::AbstractVecOrMat{Float64}, offset_b::Int64, nb::Int64)
    copyto!(c, offset_c, b, offset_b, na*nb)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('L'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(ma), Ref{BlasInt}(nb),
          Ref{Float64}(alpha), a.data, Ref{BlasInt}(ma), Ref(c, offset_c), Ref{BlasInt}(na))

    @inbounds for i = 2:ma
        x = a[i,i-1]
        indb = offset_b
        indc = offset_c + 1
        @simd for j = 1:nb
            c[indc] += x*b[indb]
            indb += ma
            indc += ma
        end
    end
end

function A_mul_B!(c::AbstractVecOrMat{Float64}, offset_c::Int64, a::QuasiUpperTriangular,
                  b::AbstractVecOrMat{Float64}, offset_b::Int64, nb::Int64)
    ma, na = size(a)
    A_mul_B!(c, offset_c, a, 1, ma, na, b, offset_b, nb)
end

#=
function A_mul_B!(c::Array{Float64,1}, offset_c::Int64,
                  a::QuasiUpperTriangular, offset_a::Int64, ma::Int64, na::Int64,
                  b::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, offset_b::Int64, nb::Int64)
    copyto!(c, offset_c, b, offset_b, na*nb)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('L'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(na), Ref{BlasInt}(nb),
          Ref{Float64}(alpha), a.data, Ref{BlasInt}(ma), Ref(c, offset_c), Ref{BlasInt}(na))

    @inbounds for i= 2:ma
        x = a[i,i-1]
        indb = offset_b
        indc = offset_c + 1
        @simd for j=1:nb
            c[indc] += x*b[indb]
            indb += ma
            indc += ma
        end
        offset_b +=  1
        offset_c += 1
    end
    c
end

function A_mul_B!(c::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, offset_c::Int64,
                  a::QuasiUpperTriangular, offset_a::Int64, ma::Int64, na::Int64,
                  b::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, offset_b::Int64, nb::Int64)
    copyto!(c, offset_c, b, offset_b, na*nb)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('L'), Ref{UInt8}('U'), Ref{UInt8}('N'), Ref{UInt8}('N'), Ref{BlasInt}(na), Ref{BlasInt}(nb),
          Ref{Float64}(alpha), a.data, Ref{BlasInt}(ma), Ref(c, offset_c), Ref{BlasInt}(na))

    offsetb_orig = offset_b
    offsetc_orig = offset_c
    @inbounds for i= 2:ma
        x = a[i,i-1]
        indb = offset_b
        indc = offset_c +1
        @simd for j=1:nb
            c[indc] += x*b[indb]
            indb += ma
            indc += ma
        end
        offset_b += 1
        offset_c += 1
    end
    c
end
=#

function At_mul_B!(c::VecOrMat{Float64}, offset_c::Int64, a::QuasiUpperTriangular{Float64}, offset_a::Int64,
                  ma::Int64, na::Int64, b::VecOrMat{Float64}, offset_b::Int64, nb::Int64)
    copyto!(c, offset_c, b, offset_b, ma*nb)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          Ref{UInt8}('L'), Ref{UInt8}('U'), Ref{UInt8}('T'), Ref{UInt8}('N'), Ref{BlasInt}(ma), Ref{BlasInt}(nb),
          Ref{Float64}(alpha), a.data, Ref{BlasInt}(ma), Ref(c, offset_c), Ref{BlasInt}(ma))

    @inbounds for i = 2:ma
        x = a[i,i-1]
        indb = offset_b + i - 1
        indc = offset_c + i - 2
        @simd for j = 1:nb
            c[indc] += x*b[indb]
            indb += ma
            indc += ma
        end
    end
end

function A_mul_Bt!(c::AbstractVector{Float64}, offset_c::Int64, a::AbstractVector{Float64}, offset_a::Int64, ma::Int64, na::Int64, b::QuasiUpperTriangular{Float64}, offset_b::Int64, nb::Int64)
    copyto!(c, offset_c, a, offset_a, ma*na)
    alpha = 1.0
    ccall((@blasfunc(dtrmm_), libblas), Cvoid,
          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
           Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}),
          'R', 'U', 'T', 'N', ma, nb,
          alpha, b.data, nb, Ref(c, offset_c), ma)

    inda = offset_a
    indc = offset_c + ma
    @inbounds for j= 2:na
        x = alpha*b[j,j-1]
        @simd for i=1:ma
            c[indc] += x*a[inda]
            inda += 1
            indc += 1
        end
    end
    c
end

function A_mul_Bt!(c::AbstractVector{Float64}, offset_c::Int64, a::AbstractVector{Float64}, offset_a::Int64, ma::Int64, na::Int64, b::QuasiUpperTriangular{Float64})
    nb = size(b, 1)
    A_mul_Bt!(c, offset_c, a, offset_a, ma, na, b, 1, nb)
end

# solver by substitution
function A_ldiv_B!(a::QuasiUpperTriangular, b::AbstractMatrix)
    m, n = size(a)
    nb, p = size(b)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    j = n
    @inbounds while j > 0
        if j == 1 || a.data[j,j-1] == 0
            a.data[j,j] == zero(a.data[j,j]) && throw(SingularException(j))
            for k = 1:p
                xj = b[j, k] = a.data[j,j] \ b[j, k]
                @simd for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                    b[i, k] -= a.data[i,j] * xj
                end
            end
            j -= 1
        else
            a11, a21, a12, a22 = a.data[j-1:j,j-1:j]
            det = a11*a22 - a12*a21
            det == zero(a.data[j,j]) && throw(SingularException(j))
            m1 = -a21/det
            m2 = a11/a21
            for k = 1:p
                x2 = m1 * (b[j-1,k] - m2 * b[j,k])
                x1 = b[j-1,k] = a21 \ (b[j,k] - a22 * x2)
                b[j,k] = x2
                @simd for i in j-2:-1:1
                    b[i,k] -= a.data[i,j-1] * x1 + a.data[i,j] * x2
                end
            end
            j -= 2
        end
    end
end

function A_rdiv_B!(a::AbstractMatrix, b::QuasiUpperTriangular)
    m, n = size(a)
    nb, p = size(b)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    @inbounds for i = 1:m
        j = 1
        while  j <= n
            if j == n || b.data[j+1,j] == 0
                b.data[j,j] == zero(b.data[j,j]) && throw(SingularException(j))
                aij = a[i,j]
                @simd for k in 1:j-1
                   aij -= a[i,k] * b.data[k,j]
                end
                a[i,j] = aij/b.data[j,j]
                j += 1
            else
                b11, b21, b12, b22 = b.data[j:j+1,j:j+1]
                det = b11*b22 - b12*b21
                det == zero(b.data[j,j]) && throw(SingularException(j))
                a1 = a[i,j]
                a2 = a[i,j+1]
                @simd for k in 1:j-1
                    a1 -= a[i,k]*b.data[k,j]
                    a2 -= a[i,k]*b.data[k,j+1]
                end
                m1 = -b21/det
                m2 = b22/b21
                a[i,j] = m1 * (a2 - m2 * a1)
                a[i,j+1] = b21 \ (a1 - b11 * a[i,j])
                j += 2
            end
        end
    end
end

function A_rdiv_Bt!(a::AbstractMatrix, b::QuasiUpperTriangular)
    m, n = size(a)
    nb, p = size(b)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    @inbounds for i = 1:m
        j = n
        while  j > 0
            if j == 1 || b.data[j,j-1] == 0
                b.data[j,j] == zero(b.data[j,j]) && throw(SingularException(j))
                aij = a[i,j]
                @simd for k = j + 1:n
                    aij -= a[i,k] * b.data[j,k]
                end
                a[i,j] = aij/b.data[j,j]
                j -= 1
            else
                b11, b21, b12, b22 = b.data[j-1:j,j-1:j]
                det = b11*b22 - b12*b21
                det == zero(b.data[j,j]) && throw(SingularException(j))
                a1 = a[i,j-1]
                a2 = a[i,j]
                @simd for k in j + 1:n
                    a1 -= a[i,k]*b.data[j-1,k]
                    a2 -= a[i,k]*b.data[j,k]
                end
                m1 = -b21/det
                m2 = b11/b21
                a[i,j] = m1 * (a1 - m2 * a2)
                a[i,j-1] = b21 \ (a2 - b22 * a[i,j])
                j -= 2
            end
        end
    end
end

function I_plus_rA_ldiv_B!(r::Float64,a::QuasiUpperTriangular, b::AbstractVector)
    m, n = size(a)
    nb, = length(b)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    j = n
    @inbounds while j > 0
        if j == 1 || r*a.data[j,j-1] == 0
            pivot = 1.0 + r*a.data[j,j]
            pivot == zero(a.data[j,j]) && throw(SingularException(j))
            b[j] = pivot \ b[j]
            xj = r*b[j]
            @simd for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                b[i] -= a.data[i,j] * xj
            end
            j -= 1
        else
            a11 = 1.0 + r*a.data[j-1,j-1]
            a21 = r*a.data[j,j-1]
            a12 = r*a.data[j-1,j]
            a22 = 1.0 + r*a.data[j,j]
            det = a11*a22 - a12*a21
            det == zero(a.data[j,j]) && throw(SingularException(j))
            m1 = -a21/det
            m2 = a11/a21
            x2 = m1 * (b[j-1] - m2 * b[j])
            x1 = b[j-1] = a21 \ (b[j] - a22 * x2)
            b[j] = x2
            x1 *= r
            x2 *= r
            @simd for i in j-2:-1:1
                b[i] -= a.data[i,j-1] * x1 + a.data[i,j] * x2
            end
            j -= 2
        end
    end
end

function I_plus_rA_plus_sB_ldiv_C!(r::Float64, s::Float64,a::QuasiUpperTriangular, b::QuasiUpperTriangular, c::AbstractVector)
    m, n = size(a)
    nb = length(c)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    j = n
    @inbounds while j > 0
        if j == 1 || r*a.data[j,j-1] + s*b.data[j,j-1] == 0
            pivot = 1.0 + r*a.data[j,j] + s*b.data[j,j]
            pivot == zero(a.data[j,j]) && throw(SingularException(j))
            c[j] = pivot \ c[j]
            xj1 = r*c[j]
            xj2 = s*c[j]
            @simd for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                c[i] -= a.data[i,j]*xj1 + b.data[i,j]*xj2
            end
            j -= 1
        else
            a11 = 1.0 + r*a.data[j-1,j-1] + s*b.data[j-1,j-1]
            a21 = r*a.data[j,j-1] + s*b.data[j,j-1]
            a12 = r*a.data[j-1,j] + s*b.data[j-1,j]
            a22 = 1.0 + r*a.data[j,j] + s*b.data[j,j]
            det = a11*a22 - a12*a21
            det == zero(a.data[j,j]) && throw(SingularException(j))
            m1 = -a21/det
            m2 = a11/a21
            x2 = m1 * (c[j-1] - m2 * c[j])
            x1 = c[j-1] = a21 \ (c[j] - a22 * x2)
            c[j] = x2
            x11 = r*x1
            x12 = s*x1
            x21 = r*x2
            x22 = s*x2
            @simd for i in j-2:-1:1
                c[i] -= a.data[i,j-1]*x11 + b.data[i,j-1]*x12 + a.data[i,j]*x21 + b.data[i,j]*x22
            end
            j -= 2
        end
    end
end

function I_plus_rA_ldiv_B!(r::Float64,a::QuasiUpperTriangular, b::AbstractMatrix)
    m, n = size(a)
    nb, p = size(b)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    j = n
    @inbounds while j > 0
        if j == 1 || r*a.data[j,j-1] == 0
            1.0 + r*a.data[j,j] == zero(a.data[j,j]) && throw(SingularException(j))
            for k = 1:p
                xj = b[j, k] = (1.0 + r*a.data[j,j]) \ b[j, k]
                @simd for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                    b[i, k] -= r*a.data[i,j] * xj
                end
            end
            j -= 1
        else
            a11 = 1.0 + r*a.data[j-1,j-1]
            a21 = r*a.data[j,j-1]
            a12 = r*a.data[j-1,j]
            a22 = 1.0 + r*a.data[j,j]
            det = a11*a22 - a12*a21
            det == zero(a.data[j,j]) && throw(SingularException(j))
            m1 = -a21/det
            m2 = a11/a21
            for k = 1:p
                x2 = m1 * (b[j-1,k] - m2 * b[j,k])
                x1 = b[j-1,k] = a21 \ (b[j,k] - a22 * x2)
                b[j,k] = x2
                @simd for i in j-2:-1:1
                    b[i,k] -= r*(a.data[i,j-1] * x1 + a.data[i,j] * x2)
                end
            end
            j -= 2
        end
    end
end

function I_plus_rA_plus_sB_ldiv_C!(r::Float64, s::Float64,a::QuasiUpperTriangular, b::QuasiUpperTriangular, c::AbstractMatrix)
    m, n = size(a)
    nb, p = size(c)
    if nb != n
        throw(DimensionMismatch("right hand side b needs first dimension of size $n, has size $(size(b,1))"))
    end
    j = n
    @inbounds while j > 0
        if j == 1 || r*a.data[j,j-1] + s*b.data[j,j-1] == 0
            1.0 + r*a.data[j,j] + s*b.data[j,j] == zero(a.data[j,j]) && throw(SingularException(j))
            for k = 1:p
                xj = c[j, k] = (1.0 + r*a.data[j,j] + s*b.data[j,j]) \ c[j, k]
                @simd for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                    c[i, k] -= (r*a.data[i,j] + s*b.data[i,j])* xj
                end
            end
            j -= 1
        else
            a11 = 1.0 + r*a.data[j-1,j-1] + s*b.data[j-1,j-1]
            a21 = r*a.data[j,j-1] + s*b.data[j,j-1]
            a12 = r*a.data[j-1,j] + s*b.data[j-1,j]
            a22 = 1.0 + r*a.data[j,j] + s*b.data[j,j]
            det = a11*a22 - a12*a21
            det == zero(a.data[j,j]) && throw(SingularException(j))
            m1 = -a21/det
            m2 = a11/a21
            for k = 1:p
                x2 = m1 * (c[j-1,k] - m2 * c[j,k])
                x1 = c[j-1,k] = a21 \ (c[j,k] - a22 * x2)
                c[j,k] = x2
                @simd for i in j-2:-1:1
                    c[i,k] -= (r*a.data[i,j-1] + s*b.data[i,j-1]) * x1 + (r*a.data[i,j] + s*b.data[i,j])* x2
                end
            end
            j -= 2
        end
    end
end

end


text:
export
    AutoLevelSetEvaluator,
    evaluate_parameters!

const HARD_BRAKE_INDEX = 4

function build_generators_evaluator(flags)
    feature_dim = 0
    flags["num_scenarios"] = 1
    flags["num_monte_carlo_runs"] = 100
    flags["prime_time"] = 1.
    flags["sampling_period"] = .1
    flags["sampling_time"] = 5.
    flags["min_num_vehicles"] = 20
    flags["max_num_vehicles"] = 20
    flags["num_lanes"] = 3
    flags["roadway_length"] = 300.
    flags["roadway_radius"] = 50.
    flags["response_time"] = .4
    flags["min_base_speed"] = 15.
    flags["max_base_speed"] = 25.

    # roadway
    roadway = gen_stadium_roadway(flags["num_lanes"],
        length = flags["roadway_length"],
        radius = flags["roadway_radius"])

    # scene
    scene = Scene(flags["max_num_vehicles"])
    scene_gen = HeuristicSceneGenerator(
        flags["min_num_vehicles"],
        flags["max_num_vehicles"],
        flags["min_base_speed"],
        flags["max_base_speed"],
        flags["min_vehicle_length"],
        flags["max_vehicle_length"],
        flags["min_vehicle_width"],
        flags["max_vehicle_width"],
        flags["min_init_dist"])

    # behavior
    context = IntegratedContinuous(flags["sampling_period"], 1)
    params = [get_aggressive_behavior_params(
                            lon_σ = flags["lon_accel_std_dev"],
                            lat_σ = flags["lat_accel_std_dev"],
                            response_time = flags["response_time"]),
                        get_passive_behavior_params(
                            lon_σ = flags["lon_accel_std_dev"],
                            lat_σ = flags["lat_accel_std_dev"],
                            response_time = flags["response_time"]),
                        get_normal_behavior_params(
                            lon_σ = flags["lon_accel_std_dev"],
                            lat_σ = flags["lat_accel_std_dev"],
                            response_time = flags["response_time"])]
                weights = WeightVec([.3,.3,.4])
    behavior_gen = PredefinedBehaviorGenerator(context, params, weights)

    # evaluator
    max_num_scenes = Int(ceil(
        (flags["prime_time"] + flags["sampling_time"]) / flags["sampling_period"]))
    rec = SceneRecord(max_num_scenes, flags["sampling_period"],
        flags["max_num_vehicles"])
    features = Array{Float64}(feature_dim, 1, flags["max_num_vehicles"])
    targets = Array{Float64}(flags["target_dim"], flags["max_num_vehicles"])
    agg_targets = Array{Float64}(flags["target_dim"], flags["max_num_vehicles"])
    ext = EmptyExtractor()
    eval = MonteCarloEvaluator(
        ext,
        flags["num_monte_carlo_runs"],
        context,
        flags["prime_time"],
        flags["sampling_time"],
        flags["veh_idx_can_change"],
        rec,
        features,
        targets,
        agg_targets)

    return (roadway, scene_gen, behavior_gen, eval)
end

type AutoLevelSetEvaluator <: LevelSetEvaluator
    roadway::Roadway
    scene_gen::SceneGenerator
    behavior_gen::BehaviorGenerator
    eval::Evaluator
    scene::Scene
    models::Dict{Int, DriverModel}
    function AutoLevelSetEvaluator(
            roadway::Roadway,
            scene_gen::SceneGenerator,
            behavior_gen::BehaviorGenerator,
            eval::Evaluator,
            scene = Scene(),
            models = Dict{Int, DriverModel}())
        new(roadway, scene_gen, behavior_gen, eval, scene, models)
    end
end

"""
Description:
    - alternate constructor directly from flags
"""
function AutoLevelSetEvaluator()
    parse_flags!(FLAGS, ARGS)
    roadway, scene_gen, behavior_gen, eval = build_generators_evaluator(FLAGS)
    return AutoLevelSetEvaluator(roadway, scene_gen, behavior_gen, eval)
end

function update_parameters!(lseval::AutoLevelSetEvaluator, lon_σ::Float64,
        lat_σ::Float64)
    for params in lseval.behavior_gen.params
        params.idm.σ = lon_σ
        params.lat.σ = lat_σ
    end
end

"""
Description:
    - Run simulations to evaluate the probability of a hard brake given the
        parameter values

Args:
    - lse: automotive level set estimator
    - params: array of parameter values at which to evaluate the target value
        params[1] = standard deviation of longitudinal acceleration (lon_σ)
            - min: 0. m/s^2
            - max: 2. m/s^2
        params[1] = standard deviation of lateral acceleration (lat_σ)
            - min: 0. m/s^2
            - max: .5 m/s^2
    - scenario_seed: random seed from which to generate the scenario
    - simulation_seed: random seed to use for simulating the scenario

Returns:
    - value of the target (hard brake probability)
        - min: 0.
        - max: 1.
"""
function evaluate_parameters!(lse::AutoLevelSetEvaluator,
        params::Array{Float64}, scenario_seed::Int = 1, simulation_seed::Int = 1)
    update_parameters!(lse, params[1], params[2])
    reset!(lse.scene_gen, lse.scene, lse.roadway, scenario_seed)
    reset!(lse.behavior_gen, lse.models, lse.scene, scenario_seed)
    evaluate!(lse.eval, lse.scene, lse.models, lse.roadway, simulation_seed)
    return mean(lse.eval.agg_targets, 2)[HARD_BRAKE_INDEX]
end


text:
### Concrete Types

struct RODESolution{T,N,uType,uType2,DType,tType,randType,P,A,IType,DE} <: AbstractRODESolution{T,N}
  u::uType
  u_analytic::uType2
  errors::DType
  t::tType
  W::randType
  prob::P
  alg::A
  interp::IType
  dense::Bool
  tslocation::Int
  destats::DE
  retcode::Symbol
  seed::UInt64
end
(sol::RODESolution)(t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(t,idxs,deriv,sol.prob.p,continuity)
(sol::RODESolution)(v,t,deriv::Type=Val{0};idxs=nothing,continuity=:left) = sol.interp(v,t,idxs,deriv,sol.prob.p,continuity)

function build_solution(
        prob::Union{AbstractRODEProblem,AbstractSDDEProblem},
        alg,t,u;W=nothing,timeseries_errors = length(u) > 2,
        dense = false,dense_errors=dense,calculate_error=true,
        interp = LinearInterpolation(t,u),
        retcode = :Default,
        seed = UInt64(0), destats=nothing, kwargs...)

  T = eltype(eltype(u))
  N = length((size(prob.u0)..., length(u)))

  if typeof(prob.f) <: Tuple
    f = prob.f[1]
  else
    f = prob.f
  end

  if has_analytic(f)
    u_analytic = Vector{typeof(prob.u0)}()
    errors = Dict{Symbol,real(eltype(prob.u0))}()
    sol = RODESolution{T,N,typeof(u),typeof(u_analytic),typeof(errors),typeof(t),typeof(W),
                       typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(
                       u,u_analytic,errors,t,W,prob,alg,interp,dense,0,destats,retcode,seed)

    if calculate_error
      calculate_solution_errors!(sol;timeseries_errors=timeseries_errors,dense_errors=dense_errors)
    end

    return sol
  else
    return RODESolution{T,N,typeof(u),Nothing,Nothing,typeof(t),
                        typeof(W),typeof(prob),typeof(alg),typeof(interp),typeof(destats)}(
                        u,nothing,nothing,t,W,prob,alg,interp,dense,0,destats,retcode,seed)
  end
end

function calculate_solution_errors!(sol::AbstractRODESolution;fill_uanalytic=true,
                                    timeseries_errors=true,dense_errors=true)

  if typeof(sol.prob.f) <: Tuple
    f = sol.prob.f[1]
  else
    f = sol.prob.f
  end

  if fill_uanalytic
    for i in 1:length(sol)
      push!(sol.u_analytic,f.analytic(sol.prob.u0,sol.prob.p,sol.t[i],sol.W[i]))
    end
  end

  if !isempty(sol.u_analytic)
    sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end]-sol.u_analytic[end])))
    if timeseries_errors
      sol.errors[:l∞] = norm(maximum(vecvecapply((x)->abs.(x),sol.u-sol.u_analytic)))
      sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply((x)->float.(x).^2,sol.u-sol.u_analytic))))
    end
    if dense_errors
      densetimes = collect(range(sol.t[1], stop=sol.t[end], length=100))
      interp_u = sol(densetimes)
      interp_analytic = [f.analytic(sol.u[1],sol.prob.p,t,sol.W(t)[1]) for t in densetimes]
      sol.errors[:L∞] = norm(maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)))
      sol.errors[:L2] = norm(sqrt(recursive_mean(vecvecapply((x)->float.(x).^2,interp_u-interp_analytic))))
    end
  end
end

function build_solution(sol::AbstractRODESolution{T,N},u_analytic,errors) where {T,N}
  RODESolution{T,N,typeof(sol.u),typeof(u_analytic),typeof(errors),typeof(sol.t),
               typeof(sol.W),typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
               sol.u,u_analytic,errors,sol.t,sol.W,sol.prob,sol.alg,sol.interp,
               sol.dense,sol.tslocation,sol.destats,sol.retcode,sol.seed)
end

function solution_new_retcode(sol::AbstractRODESolution{T,N},retcode) where {T,N}
  RODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),typeof(sol.t),
               typeof(sol.W),typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
               sol.u,sol.u_analytic,sol.errors,sol.t,sol.W,sol.prob,sol.alg,sol.interp,
               sol.dense,sol.tslocation,sol.destats,retcode,sol.seed)
end

function solution_new_tslocation(sol::AbstractRODESolution{T,N},tslocation) where {T,N}
  RODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),typeof(sol.t),
               typeof(sol.W),typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
               sol.u,sol.u_analytic,sol.errors,sol.t,sol.W,sol.prob,sol.alg,sol.interp,
               sol.dense,tslocation,sol.destats,sol.retcode,sol.seed)
end

function solution_slice(sol::AbstractRODESolution{T,N},I) where {T,N}
  RODESolution{T,N,typeof(sol.u),typeof(sol.u_analytic),typeof(sol.errors),typeof(sol.t),
               typeof(sol.W),typeof(sol.prob),typeof(sol.alg),typeof(sol.interp),typeof(sol.destats)}(
               sol.u[I],
               sol.u_analytic === nothing ? nothing : sol.u_analytic,
               sol.errors,sol.t[I],
               sol.W,sol.prob,
               sol.alg,sol.interp,
               false,sol.tslocation,sol.destats,
               sol.retcode,sol.seed)
end


text:
# convert dead reckoning back to odo

using Caesar, RoME
using Interpolations
using Distributions

using RoMEPlotting, ApproxManifoldProducts
using Gadfly, Fontconfig, Cairo
Gadfly.set_default_plot_size(35cm,25cm)

using DataFrames
using ProgressMeter
using DelimitedFiles

# mute LCM broadcasting beyond this computer
run(`$(ENV["HOME"])/mutelcm.sh`)

## load the original ascii data
include(joinpath(@__DIR__,"Plotting.jl"))
include(joinpath(@__DIR__,"SandsharkUtils.jl"))

odonoise = Matrix(Diagonal((20*[0.1;0.1;0.005]).^2))

epochs, odoDict, ppbrDict, ppbDict, pprDict, NAV = doEpochs(timestamps, rangedata, azidata, interp_x, interp_y, interp_yaw, odonoise, TSTART=50, TEND=200, SNRfloor=0.001, STRIDE=1)

# Build interpolators for x, y from LBL data
itpl_lblx = LinearInterpolation(lblkeys, lblX)
itpl_lbly = LinearInterpolation(lblkeys, lblY)

# # quick look at initial location
# itpl_lblx[epochs[1]]
# itpl_lbly[epochs[1]]
# interp_yaw[epochs[1]]
#
# DeltaX = itpl_lblx[epochs[1]] - interp_x[epochs[1]]
# DeltaY = itpl_lbly[epochs[1]] - interp_y[epochs[1]]

## differentiate odoDict back to deltas only

newSampleRate = 50
totT = (-(epochs[1]-epochs[end])*1e-9)
odoCount = round(Int,totT*newSampleRate)
odoT = range(epochs[1], epochs[end], length=odoCount) |> collect

rsmpX = map(t->interp_x(t), odoT)
rsmpY = map(t->interp_y(t), odoT)
rsmpT = map(t->interp_yaw(t), odoT)

rsmpLX = map(t->itpl_lblx(t), odoT)
rsmpLY = map(t->itpl_lbly(t), odoT)

# Gadfly.plot(x=rsmpX, y=rsmpY, Geom.path())
# Gadfly.plot(x=odoT, y=rsmpT, Geom.path())

# Gadfly.plot(x=rsmpLX, y=rsmpLY, Geom.path())

DX = extractDeltaOdo(rsmpX, rsmpY, rsmpT)

# filter put spikes
DX[1, 0.6 .< abs.(DX[1,:])] .= 0.0
DX[2, 0.6 .< abs.(DX[1,:])] .= 0.0
DX[3, 0.6 .< abs.(DX[1,:])] .= 0.0
DX[1, 0.6 .< abs.(DX[2,:])] .= 0.0
DX[2, 0.6 .< abs.(DX[2,:])] .= 0.0
DX[3, 0.6 .< abs.(DX[2,:])] .= 0.0

# Gadfly.plot(y=DX[1,:], Geom.path()) # |> PDF("/tmp/caesar/random/test.pdf");
# Gadfly.plot(y=DX[2,:], Geom.path())
# Gadfly.plot(y=DX[3,:], Geom.path())

## Back to jump free dead reckoning

nXYT = devAccumulateOdoPose2(DX' |> collect, [rsmpX[1];rsmpY[1];rsmpT[1]] )' |> collect

# Gadfly.plot(x=nXYT[1,:], y=nXYT[2,:], Geom.path()) # |> PDF("/tmp/caesar/random/test.pdf");
# Gadfly.plot(x=odoT, y=nXYT[3,:], Geom.path())
# # test new plotting function
# plotTrajectoryArrayPose2(nXYT)

# Gadfly.plot(
#   Gadfly.layer(x=nXYT[1,:], y=nXYT[2,:], Geom.path()),
#   Gadfly.layer(x=rsmpLX, y=rsmpLY, Geom.path(), Theme(default_color=colorant"red"))
# ) |> PDF("/tmp/caesar/random/test.pdf");

# |> PDF("/tmp/caesar/random/test.pdf");
# @async run(`evince /tmp/caesar/random/test.pdf`)

## add gyro bias

gyro_noise = 5e-8*randn(size(DX,2))
gyrobias = cumsum(gyro_noise)
# Gadfly.plot(y=gyrobias, Geom.line)

DXgyro = deepcopy(DX)
DXgyro[3,:] += gyrobias

mpp = MutablePose2Pose2Gaussian(MvNormal([rsmpX[1];rsmpY[1];rsmpT[1]], 1e-3*Matrix(LinearAlgebra.I, 3,3)))
dt = 1.0
nXYTgyrobias = zeros(3,size(DX,2))
Qc = 1e-6*Matrix(LinearAlgebra.I, 3,3)
for i in 1:size(DX,2)
  RoME.accumulateDiscreteLocalFrame!(mpp,DXgyro[:,i],Qc,dt)
  nXYTgyrobias[:,i] .= mpp.Zij.μ
end

# Gadfly.plot(
#   Gadfly.layer(x=nXYTgyrobias[1,:], y=nXYTgyrobias[2,:], Geom.path()),
#   Gadfly.layer(x=rsmpLX, y=rsmpLY, Geom.path(), Theme(default_color=colorant"red"))
# ) # |> PDF("/tmp/caesar/random/example.pdf");

# gfile = joinpath(ENV["HOME"], "Downloads", "gyrobias.csv")
# writedlm(gfile, gyrobias, ',')

## look at bearing and range measurents

ppbrKeys = keys(ppbrDict) |> collect |> sort

## write the new odom values to log file for real time processing.

#
using BotCoreLCMTypes, LCMCore
using Dates
#
#
lcm=LCM()
#

msg = pose_t()
rmsg = raw_t()

ppbrIdx = 1
startT = round(odoT[1]*1e-6)*1e-3 |> unix2datetime
offsetT = now() - startT
@showprogress "lcm publish loop" for idx in 1:length(odoT)
  global ppbrIdx

  # reset msg values for next round of pubs
  msg.rotation_rate = zeros(3)
  msg.vel = zeros(3)
  msg.accel = zeros(3)
  msg.orientation = [1,0,0,0.0]

  # get next timestamp to combine msgs
  msgstamp = round(odoT[idx]*1e-6)*1e-3 |> unix2datetime
  msgstamp += offsetT
  while now() < msgstamp
    sleep(0.001)
  end
  # delta odometry
  msg.utime = round(Int,odoT[idx]*1e-3)
  msg.pos = DX[1:3,idx]
  bytes = encode(msg)
  publish(lcm, "AUV_ODOMETRY",bytes)
  # odo with gyro bias
  msg.pos = DXgyro[1:3,idx]
  bytes = encode(msg)
  publish(lcm, "AUV_ODOMETRY_GYROBIAS",bytes)
  # publish the gyro bias separately
  msg.pos = zeros(3)
  msg.rotation_rate = [0;0;gyrobias[idx]]
  bytes = encode(msg)
  publish(lcm, "AUV__GYROBIAS",bytes)
  # publish ground truth
  msg.pos = [rsmpLX[idx]; rsmpLY[idx]; 0.0]
  bytes = encode(msg)
  publish(lcm, "AUV_LBL_INTERPOLATED",bytes)
  # publish magnetometer (was previously baked in)
  msg.pos = zeros(3)
  qMag = convert(Quaternion, Euler([0,0,rsmpT[idx]]))
  msg.orientation = [qMag.s; qMag.v]
  bytes = encode(msg)
  publish(lcm, "AUV_MAGNETOMETER",bytes)

  ## send possible ranges and bearings too
  epT = ppbrKeys[ppbrIdx]
  rmsgstamp = round(epT*1e-6)*1e-3 |> unix2datetime
  # @show msgstamp - (rmsgstamp + offsetT)
  if (rmsgstamp + offsetT) <= msgstamp
    rmsg.utime = round(Int64, epT*1e-3)
    # @show msg.utime, rmsg.utime
    # range data
    data = zeros(length(ppbrDict[epT].range.domain),2)
    data[:,1] = ppbrDict[epT].range.domain
    data[:,2] = ppbrDict[epT].range.weights.values
    rmsg.data = reinterpret(UInt8, vec(reshape(data,:,1))) |> collect
    rmsg.length = length(rmsg.data)
    bytes = encode(rmsg)
    publish(lcm, "AUV_RANGE_CORRL",bytes)
    # bearing data
    rmsg.data = reinterpret(UInt8, getPoints(ppbrDict[epT].bearing)[:]) |> collect
    rmsg.length = length(rmsg.data)
    bytes = encode(rmsg)
    publish(lcm, "AUV_BEARING_CORRL",bytes)

    # advance the counter
    ppbrIdx += 1
  end
end

close(lcm)

#


text:
module NSGAII

import YAML
using Cambrian
import JSON

include("evolution.jl")
include("crossover.jl")

end


text:
using IntervalArithmetic
using IntervalConstraintProgramming
using ModelingToolkit

using Test

@testset "Utilities" begin
    @test IntervalConstraintProgramming.unify_variables([:a, :c], [:c, :b]) == ([:a,:b,:c], [1,3], [3,2], [1,0,2], [0,2,1])
end
@testset "BasicContractor" begin
    @variables x y
    C = BasicContractor(x^2 + y^2)

    @test C(-∞..1, IntervalBox(-∞..∞,2)) == IntervalBox(-1..1, -1..1)

    X =IntervalBox(-1..1,2)
    @test C(X) == 0..2

    @test C((1,2)) == 5
end
@testset "Contractor" begin
    x = y = -∞..∞
    X = IntervalBox(x, y)

    C = @contractor x^2 + y^2

    @test C(-∞..1, X) == IntervalBox(-1..1, -1..1)

    X =IntervalBox(-1..1,2)
    @test C(X) == 0..2

    @test C((1,2)) == 5

end

@testset "Contractor without using macro" begin
    @variables x y
    C = Contractor(x^2 + y^2)
    @test  C(-∞..1, IntervalBox(-∞..∞, 2)) == IntervalBox(-1..1, -1..1)

end

@testset "Contractor (with macros and without macros) specifying variables explicitly" begin
    X =IntervalBox(0.5..1.5,3)
    A=-Inf..1

    C1 = @contractor(x+y, [x,y,z])
    @test C1(A,X) == IntervalBox(0.5..0.5, 0.5..0.5, 0.5..1.5)

    C2 = @contractor(y+z, [x,y,z])
    @test C2(A,X) == IntervalBox(0.5..1.5, 0.5..0.5, 0.5..0.5)

    vars = @variables x y z

    C1 = Contractor(vars, x+y)
    @test C1(A,X) == IntervalBox(0.5..0.5, 0.5..0.5, 0.5..1.5)

    C2 = Contractor(vars, y+z)
    @test C2(A,X) == IntervalBox(0.5..1.5, 0.5..0.5, 0.5..0.5)

    C1 = Contractor([x, y, z], x+y)
    @test C1(A,X) == IntervalBox(0.5..0.5, 0.5..0.5, 0.5..1.5)

    C2 = Contractor([x, y, z], y+z)
    @test C2(A,X) == IntervalBox(0.5..1.5, 0.5..0.5, 0.5..0.5)

    vars = @variables x1 x3 x2

    C = Contractor(vars, x1+x2)
    @test C(A, X) == IntervalBox(0.5..0.5, 0.5..1.5, 0.5..0.5)

end

@testset "Contractor is created by function name " begin
  vars = @variables x y
  g(x, y) = x + y
  C = Contractor(vars, g)

  @test C(-Inf..1, IntervalBox(0.5..1.5, 2)) == IntervalBox(0.5..0.5, 2)

end

@testset "Separators" begin
    II = -100..100
    X = IntervalBox(II, II)
    S = @constraint x^2 + y^2 <= 1

    @test typeof(S) <: IntervalConstraintProgramming.ConstraintSeparator

    inner, outer = S(X)
    @test inner == IntervalBox(-1..1, -1..1)
    @test outer == IntervalBox(II, II)

    X = IntervalBox(-∞..∞, -∞..∞)
    inner, outer = S(X)
    @test inner == IntervalBox(-1..1, -1..1)
    @test outer == IntervalBox(-∞..∞, -∞..∞)

end

@testset "Separators without using macros" begin
    II = -100..100
    X = IntervalBox(II, II)
    vars = @variables x y

    S = Separator(vars, x^2 + y^2 < 1)

    inner, outer = S(X)
    @test inner == IntervalBox(-1..1, -1..1)
    @test outer == IntervalBox(II, II)

    X = IntervalBox(-∞..∞, -∞..∞)
    inner, outer = S(X)
    @test inner == IntervalBox(-1..1, -1..1)
    @test outer == IntervalBox(-∞..∞, -∞..∞)

end

@testset "Separator is created by function name " begin
  vars = @variables x y
  g(x, y) = x + y < 1
  S = Separator(vars, g)

  @test S(IntervalBox(0.5..1.5, 2)) == (IntervalBox(0.5..0.5, 2), IntervalBox(0.5..1.5, 2))

end

@testset "pave" begin
    S1a = @constraint(x > 0 , [x, y])
    S1b = @constraint(y > 0 , [x, y])

    S1 = S1a ∩ S1b
    paving = pave(S1, IntervalBox(-3..3, -3..3), 2.0, 0.5)

    @test Set(paving.inner) == Set([IntervalBox(1.5..3, 0..3), IntervalBox(0..1.5, 1.5..3)])
    @test length(paving.inner) == 2
    @test isempty(paving.boundary) == false

    S2 = S1a ∪ S1b
    paving = pave(S2, IntervalBox(-3..3, -3..3), 0.1)
    @test Set(paving.inner) == Set([IntervalBox(-3..0, 0..3), IntervalBox(0..3, -3..3)])
    @test length(paving.inner) == 2
    @test isempty(paving.boundary) == false

    S3 = @constraint x^2 + y^2 <= 1
    X = IntervalBox(-∞..∞, -∞..∞)
    paving = pave(S3, X, 1.0, 0.5)

    @test Set(paving.inner) == Set([IntervalBox(Interval(0.0, 0.5), Interval(0.0, 0.8660254037844386)),
                            IntervalBox(Interval(0.0, 0.5), Interval(-0.8660254037844386, 0.0)),
                            IntervalBox(Interval(-0.5, 0.0), Interval(0.0, 0.8660254037844386)),
                            IntervalBox(Interval(-0.5, 0.0), Interval(-0.8660254037844386, 0.0))])

    @test Set(paving.boundary) == Set([ IntervalBox(Interval(0.5, 1.0), Interval(0.0, 0.8660254037844387)),
                                IntervalBox(Interval(0.0, 0.5), Interval(0.8660254037844386, 1.0)),
                                IntervalBox(Interval(0.5, 1.0), Interval(-0.8660254037844387, 0.0)),
                                IntervalBox(Interval(0.0, 0.5), Interval(-1.0, -0.8660254037844386)),
                                IntervalBox(Interval(-0.5, 0.0), Interval(0.8660254037844386, 1.0)),
                                IntervalBox(Interval(-1.0, -0.5), Interval(0.0, 0.8660254037844387)),
                                IntervalBox(Interval(-0.5, 0.0), Interval(-1.0, -0.8660254037844386)),
                                IntervalBox(Interval(-1.0, -0.5), Interval(-0.8660254037844387, 0.0))])

    @test length(paving.inner) == 4
    @test length(paving.boundary) == 8
end

# @testset "Constants" begin
#     x = y = -∞..∞
#     X = IntervalBox(x, y)

#     a = 3
#     S4 = @constraint x^2 + y^2 - $a <= 0
#     paving = pave(S4, X)

#     @test paving.ϵ == 0.01
#     @test length(paving.inner) == 1532
#     length(paving.boundary) == 1536

# end

@testset "Paving a 3D torus" begin
    S5 = @constraint (3 - sqrt(x^2 + y^2))^2 + z^2 <= 1
    x = y = z = -∞..∞
    X = IntervalBox(x, y, z)

    paving = pave(S5, X, 1.)

    @test typeof(paving) == IntervalConstraintProgramming.Paving{3, Float64}

end

@testset "Volume" begin
    x = 3..5
    @test Vol(x).bounds == 2

    V = Vol(IntervalBox(-1..1.5, 2..3.5))
    @test typeof(V) == IntervalConstraintProgramming.Vol{2, Float64}
    @test V.bounds == 3.75

end

@testset "Functions" begin
    @function f(x) = 4x;
    C1 = @contractor f(x);
    A = IntervalBox(0.5..1);
    x = IntervalBox(0..1);

    @test C1(A, x) == IntervalBox(0.125..0.25)   # x such that 4x ∈ A=[0.5, 1]

    C2 = @constraint f(x) ∈ [0.5, 0.6]
    X = IntervalBox(0..1)

    paving = pave(C2, X)
    @test length(paving.inner) == 2
    @test length(paving.boundary) == 2

    C3 = @constraint f(f(x)) ∈ [0.4, 0.8]
    @test length(paving.inner) == 2
    @test length(paving.boundary) == 2

end

@testset "Nested functions" begin
    @function f(x) = 2x
    @function g(x) = ( a = f(x); a^2 )
    @function g2(x) = ( a = f(f(x)); a^2 )

    C = @contractor g(x)
    C2 = @contractor g2(x)

    A = IntervalBox(0.5..1)
    x = IntervalBox(0..1)

    @test C(A, x) == IntervalBox(sqrt(A[1] / 4))
    @test C2(A, x) == IntervalBox(sqrt(A[1] / 16))

end

@testset "Iterated functions" begin

    @function f(x) = 2x

    A = IntervalBox(0.5..1)
    x = IntervalBox(0..1)

    @function g3(x) = ( a = (f↑2)(x); a^2 )
    C3 = @contractor g3(x)

    @test C3(A, x) == IntervalBox(sqrt(A[1] / 16))
end


text:
module Wrap

using CxxWrap

const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
const DEPSDIR = joinpath(PKGDIR, "deps")

@wrapmodule(joinpath(DEPSDIR, "libthreadpool.so"))

function __init__()
    @initcxx
end

end


text:
type Ant
  score::Int
  moves::Int
  orientation::Int # North, East, South, West
  pos_x::Int
  pos_y::Int
  lim_x::Int
  lim_y::Int
  max_moves::Int
  max_score::Int
  map::Array{Bool, 2} #(x,y)
  sources::Vector{Tuple{Int, Int}}

  Ant(max_moves::Int, map::Array{Bool, 2}) =
    build(new(0, 0, 1, 1, 1, size(map, 1), size(map, 2), max_moves, 0, map))
end

function left(a::Ant)
  a.moves += 1
  a.orientation = max(min(a.orientation - 1, 3), 0)
end

function right(a::Ant)
  a.moves += 1
  a.orientation = max(min(a.orientation + 1, 3), 0)
end

food_ahead(a::Ant) =
  a.map[next_cell_x(a), next_cell_y(a)]

next_cell_x(a::Ant) = if a.orientation == 1
  min(a.pos_x + 1, a.lim_x)
elseif a.orientation == 3
  max(a.pos_x - 1, 1)
else
  a.pos_x
end

next_cell_y(a::Ant) = if a.orientation == 0
  max(a.pos_y - 1, 1)
elseif a.orientation == 2
  min(a.pos_y + 1, a.lim_y)
else
  a.pos_y
end

function move(a::Ant)
  a.moves += 1
  a.pos_x = next_cell_x(a)
  a.pos_y = next_cell_y(a)
  if a.map[a.pos_x, a.pos_y] && a.moves <= a.max_moves && a.score <= a.max_score
    a.map[a.pos_x, a.pos_y] = false
    a.score += 1
  end
end

function reset!(a::Ant)
  a.pos_x = 1
  a.pos_y = 1
  a.orientation = 1
  a.moves = 0
  a.score = 0
  for (x,y) in a.sources
    a.map[x,y] = true
  end
end

function build(a::Ant)
  a.sources = []
  for x in 1:size(a.map, 1)
    for y in 1:size(a.map, 2)
      if a.map[x, y]
        a.max_score += 1
        push!(a.sources, (x,y))
      end
    end
  end
  return a
end

draw(a::Ant, x::Int, y::Int, b::Bool) = if x == a.pos_x && y == a.pos_y
  'X'
elseif b
  '#'
else
  ' '
end

describe(a::Ant) =
  join([join([draw(a, x, y, a.map[x,y]) for x in 1:a.lim_x],"") for y in 1:a.lim_y], "\n")

function load_trail(path::AbstractString)

  # Load the contents of the trail file.
  f = open(path)
  lines = readlines(f)
  close(f)

  # Find the dimensions of the trail.
  dims = [parse(Int, i) for i in Base.split(Base.strip(shift!(lines)), " ")]

  # Construct the trail, row by row.
  trail = Array(Bool, dims[1], dims[2])
  for y in 1:dims[2]
    for x in 1:dims[1]
      trail[x,y] =
        !(y > length(lines)) &&
        !(x > length(lines[y])) &&
        lines[y][x] == '#'
    end
  end
  return trail

end

"""
An evaluator, specialised to the artifical ant problem.
"""
type AntEvaluator <: SimpleEvaluator
  ant::Ant
  AntEvaluator(ant::Ant) = new(ant)
end

#function evaluate!(e::AntEvaluator, s::State, sc::FitnessScheme, c::Individual)
#  while e.ant.moves < e.ant.max_moves && e.ant.score < e.ant.max_score
#    execute(get(c.genome), e.ant)
#  end
#  score = e.ant.score
#  reset!(e.ant)
#  fitness(sc, score)
#end

"""
An evaluator, specialised for the artifical ant problem.

**Properties:**

* `moves::Int`, the maximum number of moves the ant is permitted to make.
* `trail::AbstractString`, the path to the file describing the trail for this problem.
"""
ant(s::Dict{Any,Any}) = AntEvaluator(Ant(s["moves"], load_trail(s["trail"])))


text:
#
# brent.jl --
#
# Find a local root or a local minimum of an univariate function by Brent's
# methods described in:
#
# [1] Richard Brent, "Algorithms for minimization without derivatives,"
#     Prentice-Hall, inc. (1973).
#
#-------------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 1973, Richard Brent.
# Copyright (C) 2015-2020, Éric Thiébaut.
#

module Brent

export
    fzero,
    fmin

"""
    fzero_atol(T)
    fzero_rtol(T)

yields default absolute and relative tolerances for Brent's `fzero` method.
Argument `T` is the floating-point type used for the computations.

"""
fzero_atol(::Type{T}) where {T<:AbstractFloat} = floatmin(T)
fzero_rtol(::Type{T}) where {T<:AbstractFloat} = eps(T)

"""
    fmin_atol(T)
    fmin_rtol(T)

yields default absolute and relative tolerances for Brent's `fmin` method.
Argument `T` is the floating-point type used for the computations.

"""
fmin_atol(::Type{T}) where {T<:AbstractFloat} = floatmin(T)
fmin_rtol(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))

# goldstep = 1/φ^2 = 2 - φ ≈ 0.3812
import Base.MathConstants: φ
goldstep(::Type{T}) where {T<:AbstractFloat} = T(2 - φ)

"""
# Van Wijngaarden–Dekker–Brent method for finding a zero of a function

    fzero([T=Float64,] f, a, b; atol=floatmin(T), rtol=eps(T)) -> (x, fx)

seeks a local root of the function `f(x)` in the interval `[a,b]`.

It is assumed that `f(a)` and `f(b)` have opposite signs (an error is raised if
this does not hold).  `fzero` returns a zero `x` in the given interval `[a,b]`
to within a tolerance: `rtol*abs(x) + atol`.

If the function value at the endpoints `a` and `b` of the search interval are
known, the method can be called as:

    fzero([T=Float64,] f, a, fa, b, fb; atol=floatmin(T), rtol=eps(T))

to save some computations.

This function has been derived from Richard Brent's F77 code ZEROIN which
itself is a slightly modified translation of the Algol 60 procedure ZERO
given in:

> Richard Brent, "Algorithms for minimization without derivatives",
> Prentice-Hall, inc. (1973).

## Arguments

* `T` is the floating-point data type to use for computations.  Defaulting to
  `Float64`.  The function `f(x)` must return a result of type `T` or which can
  be automatically promoted to type `T`.

* `f` - The user-supplied function whose local root is being sought.  Called as
  `f(x)` to evaluate the function at any `x` in the interval `[a,b]`.

* `a`, `b` - The endpoints of the initial search interval.

* `fa`, `fb` - The function values at the endpoints of the initial search
  interval.

## Keywords

* `atol` is the absolute tolerance for the solution.  The default value
  for `atol` is `floatmin(T)`.

* `rtol` is the relative tolerance for the solution.  The recommended (and
  default) value for `rtol` is `eps(T)` where `eps(T)` is the relative
  machine precision defined as the smallest representable number such that `1 +
  eps > 1'.

## Tolerance

In original Brent's method, the tolerance (denoted `tol` in the code) is
defined as:

    δ = 2⋅ϵ⋅|b| + t

where `ϵ = eps(T)/2` is the relative machine precision halved for rounded
arithmetic (see Eq. (2.9) p. 51 in Brent's book), `b` (initially one of the
bounds) is the estimated solution and `t > 0` is chosen by the caller.

The value of `ϵ` should not be decreased below `eps(T)/2`, the relative machine
precision halved for rounded arithmetic, for then rounding errors might prevent
convergence.  If `ϵ = eps(T)/2`, the error is approximately bounded by:

    abs(b - x) ≤ 6⋅ϵ⋅abs(x) + 2⋅t

with `x` and `b` the exact and estimated solutions (see Eq. (2.18) p. 52 in
Brent's book).

In this version of the code, the tolerance is computed as:

    tol = δ = rtol*abs(b) + atol

where the absolute and relative tolerances are related to the parameters `ϵ`
and `t` of Brent's method by:

    atol ≡ t
    rtol ≡ 2⋅ϵ ≥ eps(T)

## Result

The returned value is the 2-tuple `(x, fx)` with `x` the estimated value of an
abscissa for which `f` is approximately zero in `[a,b]`; `fx = f(x)` is the
function value at `x`.

"""
fzero(f, a::Real, b::Real; kwds...) = fzero(Float64, f, a, b; kwds...)

fzero(f, a::Real, fa::Real, b::Real, fb::Real; kwds...) =
    fzero(Float64, f, a, fa, b, fb; kwds...)

function fzero(::Type{T}, f, a::Real, b::Real;
               atol::Real = fzero_atol(T),
               rtol::Real = fzero_rtol(T)) where {T<:AbstractFloat}
    _fzero1(f, T(a), T(b), T(atol), T(rtol))
end

function fzero(::Type{T}, f, a::Real, fa::Real, b::Real, fb::Real;
               atol::Real = fzero_atol(T),
               rtol::Real = fzero_rtol(T)) where {T<:AbstractFloat}
    _fzero1(f, T(a), T(fa), T(b), T(fb), T(atol), T(rtol))
end

# 2 helpers are used by fzero: _fzero1 for early return and _fzero2 to refine
# the interval by Brent's method.  In these 2 functions, the type of the
# floating-point variables should be stable and equal to `T`.  Type stability
# is one of the reasons for this splitting of the code.

function _fzero1(f, a::T, b::T,
                 atol::T, rtol::T) where {T<:AbstractFloat}
    # Compute the function value at the endpoints.  Return as early as
    # possible.
    fa = T(f(a))
    fa == 0 && return (a, fa)
    fb = T(f(b))
    fb == 0 && return (b, fb)

    # Iterate to reduce the interval.
    return _fzero2(f, a, fa, b, fb, atol, rtol)
end

function _fzero1(f, a::T, fa::T, b::T, fb::T,
                 atol::T, rtol::T) where {T<:AbstractFloat}
    # Return immediately if possible.
    fa == 0 && return (a, fa)
    fb == 0 && return (b, fb)

    # Iterate to reduce the interval.
    return _fzero2(f, a, fa, b, fb, atol, rtol)
end

function _fzero2(f, a::T, fa::T, b::T, fb::T,
                 atol::T, rtol::T) where {T<:AbstractFloat}

    # Check the assumptions and the tolerance parameters.
    (fa > 0) == (fb > 0) && error("f(a) and f(b) must have different signs")
    @assert atol > 0
    @assert 0 < rtol < 1

    # Initialize.
    c, fc = a, fa
    e = d = b - a

    # Loop to improve the interval bracketing the root.
    while true
        # Make sure B is the best point so far.
        if abs(fc) < abs(fb)
            a, fa = b, fb
            b, fb = c, fc
            c, fc = a, fa
        end

        # Compute tolerance.
        tol = rtol*abs(b) + atol

        # Check for convergence.
        m = (c - b)/2
        (abs(m) ≤ tol || fb == 0) && break

        # See if a bisection is forced.
        if abs(e) < tol || abs(fa) ≤ abs(fb)
            # Bounds decreasing too slowly, use bisection.
            d = e = m
        else
            s = fb/fa
            if a == c
                # Linear interpolation.
                p = 2*m*s
                q = 1 - s
            else
                # Inverse quadratic interpolation.
                q = fa/fc
                r = fb/fc
                p = (2*m*q*(q - r) - (b - a)*(r - 1))*s
                q = (q - 1)*(r - 1)*(s - 1)
            end
            if p > 0
                q = -q
            else
                p = -p
            end
            if 2*p < min(3*m*q - tol*abs(q), abs(e*q))
                # Take the interpolation point.
                e = d
                d = p/q
            else
                # Force a bisection.
                d = e = m
            end
        end
        a, fa = b, fb
        if abs(d) > tol
            b += d
        elseif m > 0
            b += tol
        else
            b -= tol
        end
        fb = T(f(b))
        fb == 0 && break
        if (fb > 0) == (fc > 0)
            # Drop point C (make it coincident with point A) and adjust bounds
            # of interval.
            c, fc = a, fa
            e = d = b - a
        end
    end
    return (b, fb)
end

"""
# Brent's method for finding a minimum of a function

    fmin([T=Float64,] f, a, b; atol=..., rtol=...) -> (x, fx, xlo, xhi)

seeks a local minimum of a function `f(x)` in the interval `[a,b]`.

The method used is a combination of golden section search and successive
parabolic interpolation.  Convergence is never much slower than that for a
Fibonacci search.  If `f` has a continuous second derivative which is positive
at the minimum (which is not at `a` or `b`), then convergence is superlinear,
and usually of the order of about 1.3247.

The keywords `rtol` and `atol` specify a tolerance `tol = rtol*abs(x) + atol`.
The function `f` is never evaluated at two points closer than `tol`.

If `f` is a unimodal function and the computed values of `f` are always
unimodal when separated by at least `sqrt(eps)*abs(x) + (atol/3)`, then `fmin`
approximates the abscissa of the global minimum of `f` on the interval `[a,b]`
with an error less than `3*sqrt(eps)*abs(fmin) + atol`.

If `f` is not unimodal, then `fmin` may approximate a local, but perhaps
non-global, minimum to the same accuracy.

This function has been derived from Richard Brent's FORTRAN 77 code FMIN which
itself is a slightly modified translation of the Algol 60 procedure LOCALMIN
given in:

> Richard Brent, "Algorithms for minimization without derivatives,"
> Prentice-Hall, inc. (1973).

## Arguments

* `f` - The user-supplied function whose local minimum is being sought.  Called
  as `f(x)` to evaluates the function at `x`.

* `a`, `b` - The endpoints of the interval.

## Keywords

* `atol` - A positive absolute error tolerance.

* `rtol` - A positive relative error tolerance.  `rtol` should be no smaller
  than twice the relative machine precision, and preferably not much less than
  the square root of the relative machine precision.

* `T` - The floating-point data type to use for computations.  Defaulting to
  `Float64`.  The function `f(x)` must return a result of type `T` or which can
  be automatically promoted to type `T`.

## Result

The returned value is the 4-tuple `(x, fx, xlo, xhi)` where `x` is the
estimated value of an abscissa for which `f` attains a local minimum value in
`[a,b]`, `fx` is the function value at `x`, `xlo` and `xhi` are the bounds for
the position of the local minimum.

"""
fmin(f, a::Real, b::Real; kwds...) = fmin(Float64, f, a, b; kwds...)

function fmin(::Type{T}, f, a::Real, b::Real;
              atol::Real = fmin_atol(T),
              rtol::Real = fmin_rtol(T)) where {T<:AbstractFloat}
    fmin0(f, T(a), T(b), T(atol), T(rtol))
end

"""
    fmin0(f, a, b, atol, rtol)

runs Brent's algorithm given a search interval `[a,b]` and tolerances.
Floating-point type is guessed from the type of the arguments.

"""
function fmin0(f, a::T, b::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered.  Initialize the search and call
    # the real worker.
    if a > b
        a, b = b, a
    end
    x = a + goldstep(T)*(b - a)
    fx = T(f(x))
    _fmin(f, a, b, x, fx, x, fx, x, fx, atol, rtol)
end

"""
    fmin1(f, a, b, x, fx, atol, rtol)

runs Brent's algorithm with a single given initial point, `x`, inside the
search interval `[a,b]` and with one known function value: `fx = f(x)`.
Floating-point type is guessed from the type of the arguments.

"""
function fmin1(f, a::T, b::T,
               x::T, fx::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered and check that given point is in
    # the interval.  Then call main loop of Brent's algorithm.
    if a > b
        a, b = b, a
    end
    a ≤ x ≤ b || error("given point outside search interval")
    _fmin(f, a, b, x, fx, x, fx, x, fx, atol, rtol)
end

"""
    fmin2(f, a, b, x, fx, w, fw, atol, rtol)

runs Brent's algorithm with two given initial points `x` and `w` (not
necessarily distinct) inside the search interval `[a,b]` and with known
function values `fx = f(x)` and `fw = f(w)`.  Floating-point type is guessed
from the type of the arguments.

"""
function fmin2(f, a::T, b::T,
               x::T, fx::T,
               w::T, fw::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered and check that given points are in
    # the interval.
    if a > b
        a, b = b, a
    end
    (a ≤ x ≤ b && a ≤ w ≤ b) ||
        error("given point(s) outside search interval")

    # Reorder the points as assumed by Brent's algorithm.
    if fw < fx
        x, fx, w, fw = w, fw, x, fx
    end
    _fmin(f, a, b, x, fx, w, fw, w, fw, atol, rtol)
end

"""
    fmin3(f, a, b, x, fx, w, fw, v, fv, atol, rtol)

runs Brent's algorithm with three given initial points `x`, `w` and `v` (not
necessarily distinct) inside the search interval `[a,b]` and with known
function values `fx = f(x)`, `fw = f(w)` and `fv = f(v)`.  Floating-point type
is guessed from the type of the arguments.

"""
function fmin3(f, a::T, b::T,
               x::T, fx::T,
               w::T, fw::T,
               v::T, fv::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered and check that given points are in
    # the interval.
    if a > b
        a, b = b, a
    end
    (a ≤ x ≤ b && a ≤ w ≤ b && a ≤ v ≤ b) ||
        error("given point(s) outside search interval")

    # Reorder the points as assumed by Brent's algorithm.
    if fw < fx
        x, fx, w, fw = w, fw, x, fx
    end
    if fv < fx
        x, fx, v, fv = v, fv, x, fx
    end
    if fv < fw
        w, fw, v, fv = v, fv, w, fw
    end
    _fmin(f, a, b, x, fx, w, fw, v, fv, atol, rtol)
end

"""
    _fmin(f, a, b, x, fx, w, fw, v, fv, atol, rtol)

performs the main loop of Brent's algorithm assuming that all parameters are
properly set (as explained below).

Original Brent's algorithm assumes that the minimum is in the open interval
`(a,b)` with `a ≤ b` and keeps track of the following variables:

- `x`, `fx = f(x)`: position and least function value found so far;
- `w`, `fw = f(w)`: previous values of `x` and `fx`;
- `v`, `fv = f(v)`: previous values of `w` and fw`;
- `d`: computed step (new try is: u = x + d, unless d too small);
- `e`: the previous value of `d`, if a parabolic step is taken; the
  difference between the most distant current endpoint and `x`, if a
  golden step is taken.

Other variables need not be saved, notably:

- `u`, `fu = f(u)` the next point to try and its function value.

Thus the main loop of Brent's algorithm can be entered with any `x`, `w`, `v`
(not necessarily distinct) which are in `[a,b]` and such that:

    fx = f(x)  ≤  fw = f(w)  ≤  fv = f(v)

other internal variables are:

    d = x - w
    e = w - v

"""
function _fmin(f, a::T, b::T,
               x::T, fx::T,
               w::T, fw::T,
               v::T, fv::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Constant for golden step.
    c = goldstep(T)

    # Check tolerances. (Other arguments are assumed to be checked by the
    # caller.)
    @assert atol ≥ 0
    @assert 0 < rtol < 1

    # Initialize.
    d = x - w
    e = w - v

    while true

        # Compute mid-point and check the stopping criterion.
        m = (a + b)/2
        tol = rtol*abs(x) + atol
        tol2 = 2*tol
        if abs(x - m) ≤ tol2 - (b - a)/2
            return (x, fx, a, b)
        end

        # Determine next step to take.
        take_golden_step = true
        if abs(e) > tol
            # Fit a parabola (make sure final Q ≥ 0).
            r = (x - w)*(fx - fv)
            q = (x - v)*(fx - fw)
            if q > r
                p = (x - w)*r - (x - v)*q
                q = 2*(q - r)
            else
                p = (x - v)*q - (x - w)*r
                q = 2*(r - q)
            end
            if 2*abs(p) < q*abs(e) && q*(a - x) < p < q*(b - x)
                # Take the parabolic interpolation step.
                take_golden_step = false
                e = d
                d = p/q
                s = x + d
                # F must not be evaluated too close to A or B.
                if s - a < tol2 || b - s < tol2
                    d = (x < m ? tol : -tol)
                end
            end
        end
        if take_golden_step
            # Take a golden-section step.
            e = (x < m ? b : a) - x
            d = c*e
        end

        # F must not be evaluated too close to X.
        if abs(d) ≥ tol
            u = x + d
        elseif d > 0
            u = x + tol
        else
            u = x - tol
        end
        fu = T(f(u))

        # Update A, B, V, W, and X.
        if fu ≤ fx
            if u < x
                b = x
            else
                a = x
            end
            v, fv = w, fw
            w, fw = x, fx
            x, fx = u, fu
        else
            if u < x
                a = u
            else
                b = u
            end
            if fu ≤ fw || w == x
                v, fv = w, fw
	        w, fw = u, fu
            elseif fu ≤ fv || v == x || v == w
                v, fv = u, fu
            end
        end
    end
end

"""
    fminbrkt(f, x, fx, w, fw, v, fv, atol, rtol)

runs Brent's algorithm to minimize function `f` with a bracket of the minimum
defined by 3 points (`x`, `w` and `v`) with known function values (`fx`, `fw`
and `fv`) and such that the least function value is at `x` which is inside the
interval `[v,w]`.

"""
function fminbrkt(f, x::T, fx::T, w::T, fw::T,
                  v::T, fv::T, atol::T, rtol::T) where {T<:AbstractFloat}
    # Reorder the points as assumed by Brent's algorithm.
    if fv < fw
        v, fv, w, fw = w, fw, v, fv
    end
    if v < w
        a, b = v, w
    else
        a, b = w, v
    end
    (a ≤ x ≤ b && fx ≤ fw) || error("illegal bracket")
    _fmin(f, a, b, x, fx, w, fw, v, fv, atol, rtol)
end

end # module


text: □(x,y,z) = x*x + y*y + z*z # convenient function for square sum


text:
using Compose
using FileIO
using ImageIO

struct LabeledView
    label::String
    image::Array
end

"""
    labeledimage(label::String, image::Array)

Displays a string label below an image.  Useful in Weave, Juno or Jupyter.
"""
function labeledimage(label::String, image::Array)
    io = IOBuffer(maxsize = 10 * 1024 * 1024)
    save(Stream(format"PNG", io), image)
    pix = max(size(image, 1), size(image, 2))
    scaleX, scaleY = size(image, 1) / pix, size(image, 2) / pix
    return compose(
        context(),
        (
            context(
                0.5 * (1.0 - 0.8 * scaleX),
                0.5 * (1.0 - 0.8 * scaleY),
                0.8 * scaleX,
                0.8 * scaleY,
            ),
            bitmap("image/png", take!(io), 0.0, 0.0, scaleX, scaleY),
        ),
        (
            context(0.1, 0.8, 0.8, 0.2),
            text(0.5, 0.9, label, hcenter, vbottom),
            stroke("transparent"),
            fill("black"),
            fontsize(0.06h),
        ),
    )
end

"""
    labeledimages(labels::AbstractVector{<:AbstractString}, images::AbstractVector{<:AbstractArray}; ncols=3, halign = hleft)

Create a matrix of labeled images.  Useful in Weave, Juno or Jupyter.
"""
function labeledimages(
    labels::AbstractVector{<:AbstractString},
    images::AbstractVector{<:AbstractArray};
    ncols = 3,
    halign = hleft,
)
    @assert length(labels) == length(images)
    nrows = (length(images) + ncols - 1) ÷ ncols
    sc = 1.0 / max(ncols, nrows)
    tmp = []
    for i in eachindex(labels)
        r, c = (i - 1) ÷ ncols, (i - 1) % ncols
        push!(tmp, (context(c * sc, r * sc, sc, sc), labeledimage(labels[i], images[i])))
    end
    return compose(context(), tmp...)
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
module Helloworld

# Write your package code here.
"""
  greet() -> nothing

docstring of greet()
ほげほげ
"""
greet() = println("Hello World");

end


text:
a = 10

if a > 10
a = a - 5
else
a = a * 2 + 4
end

println(a)


text:
using MultipleViewGeometry, Test
using StaticArrays, LinearAlgebra

# Tests for a set of two-dimensional Cartesian points represented by homogeneous
# coordinates.
ℳ = map(Point2D, [(-10.0, -10.0),
                   (-10.0,  10.0),
                   ( 10.0, -10.0),
                   ( 10.0,  10.0)])

ℳʹ, 𝐓 = hartley_normalization(ℳ)

@test ℳʹ == map(Point2D,  [(-1.0,-1.0),
                            (-1.0, 1.0),
                            (1.0, -1.0),
                            (1.0,  1.0,)])

 @test 𝐓 == [0.1 0.0 -0.0;
             0.0 0.1 -0.0;
             0.0 0.0 1.0]
 @test hartley_transformation(ℳ) == [0.1 0.0 -0.0;
                                      0.0 0.1 -0.0;
                                      0.0 0.0 1.0]

# Tests for a set of three-dimensional Cartesian points represented by homogeneous
# coordinates.
ℳ = map(Point3D,
           [(-10.0, -10.0, -10.0),
            (-10.0, -10.0,  10.0),
            (-10.0,  10.0, -10.0),
            (-10.0,  10.0,  10.0),
            ( 10.0, -10.0, -10.0),
            ( 10.0, -10.0,  10.0),
            ( 10.0,  10.0, -10.0),
            ( 10.0,  10.0,  10.0)])

ℳʹ, 𝐓 = hartley_normalization(ℳ)
@test ℳʹ == map(Point3D,  [(-1.0,-1.0, -1.0),
                            (-1.0,-1.0,  1.0),
                            (-1.0, 1.0, -1.0),
                            (-1.0, 1.0,  1.0),
                            (1.0, -1.0, -1.0),
                            (1.0, -1.0,  1.0),
                            (1.0,  1.0, -1.0),
                            (1.0,  1.0,  1.0)])

@test 𝐓 == [0.1 0.0 0.0 -0.0;
            0.0 0.1 0.0 -0.0;
            0.0 0.0 0.1 -0.0;
            0.0 0.0 0.0 1.0]
@test hartley_transformation(ℳ) ==  [0.1 0.0 0.0 -0.0;
                                      0.0 0.1 0.0 -0.0;
                                      0.0 0.0 0.1 -0.0;
                                      0.0 0.0 0.0 1.0]