text:
"""
    svmtrain(svm::LSSVC, x::AbstractMatrix, y::AbstractVector) -> Tuple

Solves a Least Squares Support Vector **classification** problem using the Conjugate
Gradient method. In particular, it uses the Lanczos process due to the fact that the
matrices are symmetric.

# Arguments
- `svm::LSSVC`: The Support Vector Machine that contains the hyperparameters, as well as the kernel to be used.
- `x::AbstractMatrix`: The data matrix with the features. It is expected that this array is already standardized, i.e. the mean for each feature is zero and its standard deviation is one.
- `y::AbstractVector`: A vector that contains the classes. It is expected that there are only two classes, -1 and 1.

# Returns
- `Tuple`: A tuple containing `x`, `y` and the following two elements:
- `b`: Contains the bias for the decision function.
- `α`: Contains the weights for the decision function.
"""
function svmtrain(svm::LSSVC, x::AbstractMatrix, y::AbstractVector)
    n = size(y, 1)
    # Specify the keyword arguments
    kwargs = _kwargs2dict(svm)
    # We build the kernel matrix and the omega matrix
    kern_mat = _build_kernel_matrix(x; kwargs...)
    Ω = (y .* y') .* kern_mat
    H = Ω + I / svm.γ

    # * Start solving the subproblems
    # First, solve for eta
    (η, stats) = cg_lanczos(H, y)
    # Then, solve for nu
    (ν, stats) = cg_lanczos(H, ones(n))

    # We then compute s
    s = dot(y, η)

    # Finally, we solve the problem for alpha and b
    b = dot(η, ones(n)) / s
    α = ν .- (η * b)

    return (x, y, α, b)
end

"""
    svmpredict(svm::LSSVC, fits, xnew::AbstractMatrix) -> AbstractArray

Uses the information obtained from `svmtrain` such as the bias and weights to construct a decision function and predict new class values. For the _classification_ problem only.

# Arguments
- `svm::LSSVC`: The Support Vector Machine that contains the hyperparameters, as well as the kernel to be used.
- `fits`: It can be any container data structure but it must have four elements: `x`, the data matrix; `y`, the labels vector; `α`, the weights; and `b`, the bias.
- `xnew::AbstractMatrix`: The data matrix that contains the new instances to be predicted.

# Returns
- `Array`: The labels corresponding to the prediction to each of the instances in `xnew`.
"""
function svmpredict(svm::LSSVC, fits, xnew::AbstractMatrix)
    x, y, α, b = fits
    kwargs = _kwargs2dict(svm)
    @assert size(x, 1) == size(xnew, 1)
    kern_mat = _build_kernel_matrix(x, xnew; kwargs...)
    result = sum(@. kern_mat * y * α; dims=1) .+ b
    # We need to remove the trailing dimension
    result = reshape(result, size(result, 2))

    return sign.(result)
end

function svmtrain_mc(svm::LSSVC, x, y, nclass)
    num_classifiers = nclass * (nclass - 1) / 2
    # Create a collection to store the learned parameters
    class_parameters = Vector{Tuple}(undef, Int(num_classifiers))
    # Create a collection to store the class pairs
    class_pairs = similar(class_parameters)
    c_idx = 1 # For keeping track of the classifiers

    for idx = 1:nclass - 1
        # Get the elements and indices for the first class
        a_class, a_idxs = _find_and_copy(idx, y)
        a_class .= 1.0 # The first class is encoded as 1.0
        for jdx = (idx + 1):nclass
            # Get the elements and indices for the second class
            b_class, b_idxs = _find_and_copy(jdx, y)
            b_class .= -1.0 # The second class is encoded as -1.0

            # Train a binary classification problem with the first and second classes
            all_indxs = vcat(a_idxs, b_idxs) # Join all indices
            all_classes = vcat(a_class, b_class) # Join both classes

            # Always copy the model and use views for the samples
            samples = @view(x[:, all_indxs])
            # Solve the binary classification problem between classes idx and jdx
            new_svm = deepcopy(svm)
            fits = svmtrain(new_svm, samples, all_classes)

            # We save the parameters of that classifier
            class_parameters[c_idx] = fits
            # We save the class pairs for decoding them later
            class_pairs[c_idx] = (idx, jdx)
            c_idx += 1
        end
    end

    return (deepcopy(svm), class_parameters, class_pairs)
end

function svmpredict_mc(fits, Xnew::AbstractMatrix)
    # Extract the model, parameters and class codes
    svm, params, pairs = fits
    nclass = length(pairs) # The third element are the class pairs
    # This will hold all the predictions for the classifiers
    pooled_predictions = Matrix{eltype(Xnew)}(undef, nclass, size(Xnew, 2))

    for (idx, p, s) in zip(1:nclass, params, pairs)
        # For prediction, we just need the fitted parameters
        prediction = svmpredict(svm, p, Xnew)

        # Now that we have predictions, decode them using the class pairs
        broadcast!(x -> x == 1.0 ? s[1] : s[2], prediction, prediction)

        # Save the decoded predictions to our pooling collection
        pooled_predictions[idx, :] = prediction
    end

    # Apply the voting scheme
    results = _predictions_by_votes(pooled_predictions)

    return results
end

"""
    svmtrain(svm::LSSVR, x::AbstractMatrix, y::AbstractVector) -> Tuple

Solves a Least Squares Support Vector **regression** problem using the Conjugate Gradient
method. In particular, it uses the Lanczos process due to the fact that the matrices are
symmetric.

# Arguments
- `svm::LSSVR`: The Support Vector Machine that contains the hyperparameters, as well as the kernel to be used.
- `x::AbstractMatrix`: The data matrix with the features. It is expected that this array is already standardized, i.e. the mean for each feature is zero and its standard deviation is one.
- `y::AbstractVector`: A vector that contains the continuous value of the function estimation. The elements can be any subtype of `Real`.

# Returns
- `Tuple`: A tuple containing `x` and the following two elements:
- `b`: Contains the bias for the decision function.
- `α`: Contains the weights for the decision function.
"""
function svmtrain(svm::LSSVR, x::AbstractMatrix, y::AbstractVector)
    n = size(y, 1)
    # Specify the keyword arguments
    kwargs = _kwargs2dict(svm)
    # We build the kernel matrix and the omega matrix
    kern_mat = _build_kernel_matrix(x; kwargs...)
    H = kern_mat + I / svm.γ

    # * Start solving the subproblems
    # First, solve for eta
    (η, stats) = cg_lanczos(H, ones(n))
    # Then, solve for nu
    (ν, stats) = cg_lanczos(H, y)

    # We then compute s
    s = dot(ones(n), η)

    # Finally, we solve the problem for alpha and b
    b = dot(η, y) / s
    α = ν .- (η * b)

    return (x, α, b)
end

"""
    svmpredict(svm::LSSVR, fits, xnew::AbstractMatrix) -> AbstractArray

Uses the information obtained from `svmtrain` such as the bias and weights to construct a
decision function and predict the new values of the function. For the _regression_
problem only.

# Arguments
- `svm::LSSVR`: The Support Vector Machine that contains the hyperparameters, as well as the kernel to be used.
- `fits`: It can be any container data structure but it must have four elements: `x`, the data matrix; `y`, the labels vector; `α`, the weights; and `b`, the bias.
- `xnew::AbstractMatrix`: The data matrix that contains the new instances to be predicted.

# Returns
- `Array`: The labels corresponding to the prediction to each of the instances in `xnew`.
"""
function svmpredict(svm::LSSVR, fits, xnew::AbstractMatrix)
    x, α, b = fits
    @assert size(x, 1) == size(xnew, 1)

    kwargs = _kwargs2dict(svm)
    kern_mat = _build_kernel_matrix(x, xnew; kwargs...)
    result = sum(kern_mat .* α; dims=1) .+ b

    # We need to remove the trailing dimension
    result = reshape(result, size(result, 2))

    return result
end

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
using Documenter, AMDGPU

makedocs(
    sitename="AMDGPU.jl",
    pages = [
        "Home" => "index.md",
        "Quick Start" => "quickstart.md",
        "Kernel Dependencies" => "kernel_deps.md",
        "Global Variables" => "globals.md",
        "Exceptions" => "exceptions.md",
        "Memory" => "memory.md",
        "Intrinsics" => [
            "Execution Control" => "execution_control.md",
        ],
        "API Reference" => "api.md"
    ]
)

ext: jl


text:
@testset "s00" begin
    x = 0.5791308486706011000e-3
    y = 0.4020579816732961219e-4
    s = s00(2400000.5, 53736.0, x, y)
    @test s ≈ -0.1220036263270905693e-7 atol=1e-18
end

@testset "s00a" begin
    s = s00a(2400000.5, 52541.0)
    @test s ≈ -0.1340684448919163584e-7 atol=1e-18
end

@testset "s00b" begin
    s = s00b(2400000.5, 52541.0)
    @test s ≈ -0.1340695782951026584e-7 atol=1e-18
end

@testset "s06" begin
    x = 0.5791308486706011000e-3
    y = 0.4020579816732961219e-4
    s = s06(2400000.5, 53736.0, x, y)
    @test s ≈ -0.1220032213076463117e-7 atol=1e-18
end

@testset "s06a" begin
    s = s06a(2400000.5, 52541.0)
    @test s ≈ -0.1340680437291812383e-7 atol=1e-18
end

@testset "s2c" begin
    c = s2c(3.0123, -0.999)
    @test c[1] ≈ -0.5366267667260523906 atol=1e-12
    @test c[2] ≈ 0.0697711109765145365 atol=1e-12
    @test c[3] ≈ -0.8409302618566214041 atol=1e-12
end

@testset "s2p" begin
    p = s2p(-3.21, 0.123, 0.456)
    @test p[1] ≈ -0.4514964673880165228 atol=1e-12
    @test p[2] ≈ 0.0309339427734258688 atol=1e-12
    @test p[3] ≈ 0.0559466810510877933 atol=1e-12
end

@testset "s2pv" begin
    pv = s2pv(-3.21, 0.123, 0.456, -7.8e-6, 9.01e-6, -1.23e-5)
    @test pv[1][1] ≈ -0.4514964673880165228 atol=1e-12
    @test pv[1][2] ≈ 0.0309339427734258688 atol=1e-12
    @test pv[1][3] ≈ 0.0559466810510877933 atol=1e-12
    @test pv[2][1] ≈ 0.1292270850663260170e-4 atol=1e-16
    @test pv[2][2] ≈ 0.2652814182060691422e-5 atol=1e-16
    @test pv[2][3] ≈ 0.2568431853930292259e-5 atol=1e-16
end

@testset "s2xpv" begin
    s1 = 2.0
    s2 = 3.0
    pv = [[0.3,1.2,-2.5],
          [0.5,2.3,-0.4]]
    spv = ERFA._s2xpv(s1, s2, pv)
    @test spv[1][1] ≈ 0.6 atol=1e-12
    @test spv[1][2] ≈ 2.4 atol=1e-12
    @test spv[1][3] ≈ -5.0 atol=1e-12
    @test spv[2][1] ≈ 1.5 atol=1e-12
    @test spv[2][2] ≈ 6.9 atol=1e-12
    @test spv[2][3] ≈ -1.2 atol=1e-12
    pve = [[1.2,-2.5],
          [0.5,2.3,-0.4]]
    @test_throws ArgumentError ERFA._s2xpv(s1, s2, pve)
    @test spv == [s1 .* pv[1], s2 .* pv[2]]
end

@testset "sepp" begin
    a = [1.,0.1,0.2]
    b = [-3.,1e-3,0.2]
    s = sepp(a, b)
    @test s ≈ 2.860391919024660768 atol=1e-12
    @test_throws ArgumentError sepp(a[1:2], b)
    @test_throws ArgumentError sepp(a, b[1:2])
end

@testset "seps" begin
    s = seps(1., .1, .2, -3.)
    @test s ≈ 2.346722016996998842 atol=1e-14
end

@testset "sp00" begin
    s = sp00(2400000.5, 52541.0)
    @test s ≈ -0.6216698469981019309e-11 atol=1e-12
end

@testset "starpm" begin
    ra1 =   0.01686756
    dec1 = -1.093989828
    pmr1 = -1.78323516e-5
    pmd1 =  2.336024047e-6
    px1 =   0.74723
    rv1 = -21.6
    ra2, dec2, pmr2, pmd2, px2, rv2 = starpm(ra1, dec1, pmr1, pmd1, px1, rv1,
                                                  2400000.5, 50083.0, 2400000.5, 53736.0)
    @test ra2 ≈ 0.01668919069414256149 atol=1e-13
    @test dec2 ≈ -1.093966454217127897 atol=1e-13
    @test pmr2 ≈ -0.1783662682153176524e-4 atol=1e-17
    @test pmd2 ≈ 0.2338092915983989595e-5 atol=1e-17
    @test px2 ≈ 0.7473533835317719243 atol=1e-13
    @test rv2 ≈ -21.59905170476417175 atol=1e-11
end

@testset "starpv" begin
    ra =   0.01686756
    dec = -1.093989828
    pmr = -1.78323516e-5
    pmd =  2.336024047e-6
    px =   0.74723
    rv = -21.6
    pv = starpv(ra, dec, pmr, pmd, px, rv)
    @test pv[1][1] ≈ 126668.5912743160601 atol=1e-10
    @test pv[1][2] ≈ 2136.792716839935195 atol=1e-12
    @test pv[1][3] ≈ -245251.2339876830091 atol=1e-10
    @test pv[2][1] ≈ -0.4051854008955659551e-2 atol=1e-13
    @test pv[2][2] ≈ -0.6253919754414777970e-2 atol=1e-15
    @test pv[2][3] ≈ 0.1189353714588109341e-1 atol=1e-13
end

@testset "sxp" begin
    s = 2.0
    p = [0.3,1.2,-2.5]
    sp = ERFA._sxp(s, p)
    @test sp[1] ≈ 0.6 atol=0.0
    @test sp[2] ≈ 2.4 atol=0.0
    @test sp[3] ≈ -5.0 atol=0.0
    @test_throws ArgumentError ERFA._sxp(s, p[1:2])
    @test sp == s .* p
end

@testset "sxpv" begin
    s = 2.0
    pv = [[0.3,1.2,-2.5],[0.5,3.2,-0.7]]
    spv = ERFA._sxpv(s, pv)
    @test spv[1][1] ≈ 0.6 atol=0.0
    @test spv[1][2] ≈ 2.4 atol=0.0
    @test spv[1][3] ≈ -5.0 atol=0.0
    @test spv[2][1] ≈ 1.0 atol=0.0
    @test spv[2][2] ≈ 6.4 atol=0.0
    @test spv[2][3] ≈ -1.4 atol=0.0
    pve = [[1.2,-2.5],[0.5,3.2,-0.7]]
    @test_throws ArgumentError ERFA._sxpv(s, pve)
    @test spv == s .* pv
end


text:
# Reader Helper
# =============
#
# Utilities to generate file readers in BioJulia packages.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

module ReaderHelper

import Automa
import BufferedStreams

@inline function anchor!(stream::BufferedStreams.BufferedInputStream, p, immobilize = true)
    stream.anchor = p
    stream.immobilized = immobilize
    return stream
end

@inline function upanchor!(stream::BufferedStreams.BufferedInputStream)
    @assert stream.anchor != 0 "upanchor! called with no anchor set"
    anchor = stream.anchor
    stream.anchor = 0
    stream.immobilized = false
    return anchor
end

function ensure_margin!(stream::BufferedStreams.BufferedInputStream)
    if stream.position * 20 > length(stream.buffer) * 19
        BufferedStreams.shiftdata!(stream)
    end
    return nothing
end

@inline function resize_and_copy!(dst::Vector{UInt8}, src::Vector{UInt8}, r::UnitRange{Int})
    return resize_and_copy!(dst, 1, src, r)
end

@inline function resize_and_copy!(dst::Vector{UInt8}, dstart::Int, src::Vector{UInt8}, r::UnitRange{Int})
    rlen = length(r)
    if length(dst) != dstart + rlen - 1
        resize!(dst, dstart + rlen - 1)
    end
    copyto!(dst, dstart, src, first(r), rlen)
    return dst
end

@inline function append_from_anchor!(dst::Vector{UInt8}, dstart::Int, stream::BufferedStreams.BufferedInputStream, p::Int)
    return resize_and_copy!(dst, dstart, stream.buffer, upanchor!(stream):p)
end

function generate_index_function(record_type, machine, init_code, actions; kwargs...)
    kwargs = Dict(kwargs)
    context = Automa.CodeGenContext(
        generator = get(kwargs, :generator, :goto),
        checkbounds = get(kwargs, :checkbounds, false),
        loopunroll = get(kwargs, :loopunroll, 0)
    )
    quote
        function index!(record::$(record_type))
            data = record.data
            p = 1
            p_end = p_eof = sizeof(data)
            initialize!(record)
            $(init_code)
            cs = $(machine.start_state)
            $(Automa.generate_exec_code(context, machine, actions))
            if cs != 0
                throw(ArgumentError(string("failed to index ", $(record_type), " ~>", repr(String(data[p:min(p+7,p_end)])))))
            end
            @assert isfilled(record)
            return record
        end
    end
end

function generate_readheader_function(reader_type, metainfo_type, machine, init_code, actions, finish_code=:())
    quote
        function readheader!(reader::$(reader_type))
            _readheader!(reader, reader.state)
        end

        function _readheader!(reader::$(reader_type), state::BioCore.Ragel.State)
            stream = state.stream
            BioCore.ReaderHelper.ensure_margin!(stream)
            cs = state.cs
            linenum = state.linenum
            data = stream.buffer
            p = stream.position
            p_end = stream.available
            p_eof = -1
            finish_header = false
            record = $(metainfo_type)()

            $(init_code)

            while true
                $(Automa.generate_exec_code(Automa.CodeGenContext(generator=:table), machine, actions))

                state.cs = cs
                state.finished = cs == 0
                state.linenum = linenum
                stream.position = p

                if cs < 0
                    error("$($(reader_type)) file format error on line ", linenum)
                elseif finish_header
                    $(finish_code)
                    break
                elseif p > p_eof ≥ 0
                    error("incomplete $($(reader_type)) input on line ", linenum)
                else
                    hits_eof = BufferedStreams.fillbuffer!(stream) == 0
                    p = stream.position
                    p_end = stream.available
                    if hits_eof
                        p_eof = p_end
                    end
                end
            end
        end
    end
end

function generate_read_function(reader_type, machine, init_code, actions; kwargs...)
    kwargs = Dict(kwargs)
    context = Automa.CodeGenContext(
        generator=get(kwargs, :generator, :goto),
        checkbounds=get(kwargs, :checkbounds, false),
        loopunroll=get(kwargs, :loopunroll, 0)
    )
    quote
        function Base.read!(reader::$(reader_type), record::eltype($(reader_type)))::eltype($(reader_type))
            return _read!(reader, reader.state, record)
        end

        function _read!(reader::$(reader_type), state::BioCore.Ragel.State, record::eltype($(reader_type)))
            stream = state.stream
            BioCore.ReaderHelper.ensure_margin!(stream)
            cs = state.cs
            linenum = state.linenum
            data = stream.buffer
            p = stream.position
            p_end = stream.available
            p_eof = -1
            found_record = false
            initialize!(record)

            $(init_code)

            if state.finished
                throw(EOFError())
            end

            while true
                $(Automa.generate_exec_code(context, machine, actions))

                state.cs = cs
                state.finished |= cs == 0
                state.linenum = linenum
                stream.position = p

                if cs < 0
                    error($(reader_type), " file format error on line ", linenum, " ~>", repr(String(data[p:min(p+7,p_end)])))
                elseif found_record
                    break
                elseif cs == 0
                    throw(EOFError())
                elseif p > p_eof ≥ 0
                    error("incomplete $($(reader_type)) input on line ", linenum)
                elseif BufferedStreams.available_bytes(stream) < 64
                    hits_eof = BufferedStreams.fillbuffer!(stream) == 0
                    p = stream.position
                    p_end = stream.available
                    if hits_eof
                        p_eof = p_end
                    end
                end
            end

            @assert isfilled(record)
            return record
        end
    end
end

end


text:
module SparseDeconvolution

using BenchmarkTools
using StructuredOptimization
using JuMP, MathProgBase, SCS, ECOS
using AbstractOperators
using RIM
using PyPlot

function set_up(;opt="")

	srand(123)
	Fs = 4000        # sampling frequency
	SNR = 15
	env = AcEnv(Fs)

	Lx,Ly,Lz = 4.,5.,3. # room dimensions
	T60 = 0.3           # reverberation time
	geo = CuboidRoom(Lx,Ly,Lz,T60,env)

	xs = [0.5 0.5 0.5]'            # src pos (in meters)
	xr = [Lx-0.1 Ly-0.3 Lz-0.2]'   # mic pos
	Nh = div(Fs,4)                 # IR samples
	Nx = div(Fs,2)                 # x  samples

	t = linspace(0,1/Fs*Nx,Nx)

	h = rim(xs,xr,Nh,geo,env)[:] # Impuse Response

	x = full(sprandn(Nx, 0.05))  # sparse input

	y0 = conv(x,h)                                      # output signal
	y = y0+10^(-SNR/10)*sqrt(var(y0))*randn(length(y0)) # add noise

	H = Conv(Float64,size(x),h)                              # Abstract Operator
	T = hcat([[zeros(i);h;zeros(Nx-1-i)] for i = 0:Nx-1]...) # Full Matrix
	lambda = 1e-2*vecnorm(H'*y,Inf)

    x_test = randn(Nx)
    @assert norm(T*x_test-H*x_test) < 1e-8

	if opt == "MatrixFree"
		setup = y, H, lambda, Nx
	elseif opt == "JuMP"
		M = Model()
		@variables M begin
			x0[1:Nx]
			tt[1:Nx]
			w
		end
		@objective(M,Min,[0.5;lambda*ones(Nx)]'*[w;tt])
		@constraint(M, soc, norm( [1-w;2*(T*x0-y)] ) <= 1+w)
		@constraint(M,  x0 .<= tt)
		@constraint(M, -tt .<= x0)
		setup = y, T, M, x0
	else
		setup = y, T, lambda, Nx
	end
	return setup, t, x

end

function run_demo()

	setup, t, x = set_up()
    tol = 1e-6
    verb = 1
    solver = PANOC

	println("Solving Unregularized problem")
	xu = setup[2]\setup[1]

	println("Solving Regularized problem with Full Matrix")
	slv = solver(tol = tol, verbose = verb)
	@time x0, = solve_problem(slv, setup...)
	xm = copy(~x0)

	setup, t, x = set_up(opt = "MatrixFree")

	println("Solving Regularized problem with Abstract Operator")
	slv = solver(tol = tol, verbose = verb)
	@time x0, = solve_problem(slv, setup...)
	x1 = copy(~x0)

	println("  regularized MSE: $( 20*log10(norm(x1-x)/norm(x)) )")
	println("unregularized MSE: $( 20*log10(norm(xu-x)/norm(x)) )")

	return t, x, x1, xu
end

function run_demo_JuMP()

	setup, t, x = set_up(opt = "JuMP")

	println("Solving Unregularized problem")
	xu = setup[2]\setup[1]

	println("Solving Regularized problem with conic solver")
	slv = SCSSolver()
	@time x1, = solve_problem(slv, setup...)

	println("  regularized MSE: $( 20*log10(norm(x1-x)/norm(x)) )")
	println("unregularized MSE: $( 20*log10(norm(xu-x)/norm(x)) )")

	return t, x, x1, xu
end

#StructuredOptimization Matrix Free
function solve_problem(slv::S, y, H::A, lambda, Nx) where {S <: StructuredOptimization.ForwardBackwardSolver, A <: AbstractOperator}
	x0 = StructuredOptimization.Variable(Nx)
	it, = @minimize ls(H*x0-y)+lambda*norm(x0,1) with slv
	return x0, it
end

#StructuredOptimization non-Matrix Free
function solve_problem(slv::S, y, T::A, lambda, Nx) where {S <: StructuredOptimization.ForwardBackwardSolver, A <: AbstractMatrix }
	x0 = StructuredOptimization.Variable(Nx)
	it, = @minimize ls(T*x0-y)+lambda*norm(x0,1) with slv
	return x0, it
end

#JuMP non-Matrix Free
function solve_problem(slv::S, y, T, M, x0) where {S <: MathProgBase.SolverInterface.AbstractMathProgSolver}
	JuMP.setsolver(M, slv)
	JuMP.solve(M)
	return getvalue(x0), 0
end

function benchmark(;verb = 0, samples = 5, seconds = 100, tol = 1e-6)

	suite = BenchmarkGroup()

	opt     = [
           #"JuMP",
		   #"JuMP",
		   "",
		   "",
		   "",
		   "",
          ]
	solvers = [
		   #"ECOSSolver",
		   #"SCSSolver",
		   "PANOC",
		   "ZeroFPR",
		   "FPG",
		   "PG"]
	slv_opt = [
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)",
		   "(verbose = $verb, tol = $tol)"]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup, t, x = set_up(opt = opt[i])
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable((x0,it) = solve_problem(solver, setup...),
			       setup = (
					it = 0;
					setup = deepcopy($setup);
					solver = deepcopy($solver) ),
			       teardown = (
					  $its[$solvers[$i]] = it;
					  ),
			       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("SparseDeconvolution its")
	println(its)
	return results
end

function benchmarkMatrixFree(;verb = 0, samples = 5, seconds = 100, tol = 1e-6)

	suite = BenchmarkGroup()

	solvers = ["PANOC", "ZeroFPR", "FPG", "PG"]
	slv_opt = ["(verbose = $verb, tol = $tol)",
               "(verbose = $verb, tol = $tol)",
               "(verbose = $verb, tol = $tol)",
               "(verbose = $verb, tol = $tol)"]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup, t, x = set_up(opt = "MatrixFree")
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable((x0,it) = solve_problem(solver, setup...),
			       setup = (
					it = 0;
					setup = deepcopy($setup);
					solver = deepcopy($solver) ),
			       teardown = (
					  $its[$solvers[$i]] = it;
					  ),
			       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("SparseDeconvolution Matrix Free its")
	println(its)
	return results
end

function show_results(t, x, x1, xu)

	figure()
	plot(t,x )
	plot(t,x1)
	plot(t,xu)

end

end


text:
# the largest license in `https://github.com/spdx/license-list-data/blob/v3.10/text`
# is the `APL-1.0` license which is 45 KB.
# We take a factor of 10 larger, to allow for
# compound licenses.
const MAX_LICENSE_SIZE_IN_BYTES = 45958 * 10

# Based on
# https://github.com/ericphanson/LicenseCheck.jl/issues/2#issue-805995984
const LICENSE_NAMES = let
    names = ("LICENSE", "LICENCE", "COPYING", "NOTICE", "COPYRIGHT") ∪ OSI_LICENSES
    name_cases = (uppercase, uppercasefirst, lowercase)
    ext_cases = (uppercase, lowercase)
    extensions = (".md", ".txt", "", ".rst")
    Set(string(case(name), extcase(ext))
        for name in names, case in name_cases, ext in extensions, extcase in ext_cases)
end

# only ~500 items this way
const LOWERCASE_LICENSE_NAMES = Set(lowercase(lic) for lic in LICENSE_NAMES)

const LICENSE_TABLE_TYPE = Vector{NamedTuple{(:license_filename, :licenses_found,
                                              :license_file_percent_covered),
                                             Tuple{String,Vector{String},Float64}}}
const LICENSE_TABLE_TYPE_STRING = "Vector{@NamedTuple{license_filename::String, licenses_found::Vector{String}, license_file_percent_covered::Float64}}"

# like `readdir`, but returns only files
readfiles(dir) = filter!(f -> isfile(joinpath(dir, f)), readdir(dir))

# constructs a table of `licensecheck` results
function license_table(dir, names; validate_strings=true, validate_paths=true)
    table = LICENSE_TABLE_TYPE()
    for lic in names
        path = joinpath(dir, lic)
        validate_paths && (isfile(path) || continue)
        text = read(path, String)
        # In licensecheck we're going to convert the string to a `Cstring` with
        # `unsafe_convert`, so the string can't have any NUL character
        validate_strings && ((isvalid(String, text) && !Base.containsnul(text)) || continue)
        lc = licensecheck(text)
        if lc.license_file_percent_covered > 0
            push!(table, (; license_filename=lic, lc...))
        end
    end
    sort!(table; by=x -> x.license_file_percent_covered, rev=true)
    return table
end

"""
    find_licenses_by_list_intersection(dir) -> $(LICENSE_TABLE_TYPE_STRING)

Checks to see if any license name in `LicenseCheck.LICENSE_NAMES` exists
in `dir`, and if so, calls [`licensecheck`](@ref) on it. Returns the results
of all existing licenses and their `licensecheck` values, sorted from highest
`license_file_percent_covered` to lowest.

Operates by filtering the results of `readdir`, which should be efficient
for small and moderately sized directories. See [`find_licenses_by_list`](@ref)
for an alternate approach for very large directories.
"""
function find_licenses_by_list_intersection(dir; files=readfiles(dir))
    names = filter!(lic -> lowercase(lic) ∈ LOWERCASE_LICENSE_NAMES, files)
    return license_table(dir, names; validate_paths=false)
end

"""
    find_licenses_by_list(dir) -> $(LICENSE_TABLE_TYPE_STRING)

Checks to see if any license name in `LicenseCheck.LICENSE_NAMES` exists
in `dir`, and if so, calls [`licensecheck`](@ref) on it. Returns the results
of all existing licenses and their `licensecheck` values, sorted from highest
`license_file_percent_covered` to lowest.

This function does not ever call `readdir(dir)` and instead just checks if each
of the $(length(LICENSE_NAMES)) names in `LicenseCheck.LICENSE_NAMES` exists,
so it can work on extremely large directories.

On case-insensitive filesystems, this can return multiple results for the same file,
with different cases.
"""
find_licenses_by_list(dir) = license_table(dir, LICENSE_NAMES)

"""
    find_licenses_by_bruteforce(dir; max_bytes = LicenseCheck.MAX_LICENSE_SIZE_IN_BYTES) -> $(LICENSE_TABLE_TYPE_STRING)

Calls [`licensecheck`](@ref) on every plaintext file in `dir` whose size is less than `max_bytes`,
returning the results as a table. The parameter `max_bytes` defaults to $(MAX_LICENSE_SIZE_IN_BYTES ÷ 1000) KiB.
"""
function find_licenses_by_bruteforce(dir; max_bytes=MAX_LICENSE_SIZE_IN_BYTES,
                                     files=readfiles(dir))
    names = filter!(file -> stat(joinpath(dir, file)).size < max_bytes, files)
    return license_table(dir, names; validate_paths=false)
end

const CUTOFF = 100

"""
    find_licenses(dir; allow_brute=true, max_bytes = MAX_LICENSE_SIZE_IN_BYTES) -> $(LICENSE_TABLE_TYPE_STRING)

Compiles a table of possible licenses at the top-level of a directory `dir` with their path and the results of [`licensecheck`](@ref), sorted by `license_file_percent_covered`. Uses [`find_licenses_by_bruteforce`](@ref) for directories
with size less than $CUTOFF and [`find_licenses_by_list_intersection`](@ref) for larger directories.

Simply acts as an alias for [`find_licenses_by_list_intersection`](@ref) if `allow_brute=false`.

## Example

```julia
julia> find_licenses(".")
1-element Vector{NamedTuple{(:license_filename, :licenses_found, :license_file_percent_covered), Tuple{String, Vector{String}, Float64}}}:
 (license_filename = "LICENSE", licenses_found = ["MIT"], license_file_percent_covered = 98.82352941176471)

```
"""
function find_licenses(dir; allow_brute=true, max_bytes=MAX_LICENSE_SIZE_IN_BYTES)
    files = readfiles(dir)
    if allow_brute && (length(files) < CUTOFF)
        return find_licenses_by_bruteforce(dir; files=files, max_bytes=max_bytes)
    else
        return find_licenses_by_list_intersection(dir; files=files)
    end
end

"""
    find_license(dir; max_bytes = MAX_LICENSE_SIZE_IN_BYTES) -> Union{Nothing, @NamedTuple{license_filename::String, licenses_found::Vector{String}, license_file_percent_covered::Float64}}

Returns the license with the highest `license_file_percent_covered` from [`find_licenses`](@ref). If file
is found with any license content, returns `nothing`.

## Example
```julia
julia> find_license(".")
(license_filename = "LICENSE", licenses_found = ["MIT"], license_file_percent_covered = 98.82352941176471)

```
"""
function find_license(dir; kwargs...)
    table = find_licenses(dir; kwargs...)
    return isempty(table) ? nothing : table[1]
end


text:
@testset "Modules" begin
    # Default constructors
    fam = FAM()
    dam = DAM()
    sfam = SFAM()
    dvfa = DVFA()
    ddvfa = DDVFA()

    # Specify constructors
    fam_2 = FAM(opts_FAM())
    dam_2 = DAM(opts_DAM())
    sfam_2 = SFAM(opts_SFAM())
    dvfa_2 = DVFA(opts_DVFA())
    ddvfa_2 = DDVFA(opts_DDVFA())
end # @testset "Modules"


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module PrimeSwingFactorial

using PrimesIterator, Products, Nemo

export ModulePrimeSwingFactorial
export PSfactorial, Swing

"""
 Cf. P. Luschny, [Swing, divide and conquer the factorial](https://oeis.org/A000142/a000142.pdf), excerpt.

* PSfactorial, Swing
"""
const ModulePrimeSwingFactorial = ""

const SwingOddpart = [1,1,1,3,3,15,5,35,35, 315, 63, 693, 231,
   3003, 429, 6435, 6435, 109395,12155,230945,46189,969969,
   88179,2028117, 676039,16900975,1300075,35102025,5014575,
   145422675,9694845,300540195,300540195]

"""
Computes the odd part of the swinging factorial ``n≀`` (cf. A163590).
"""
function swing_oddpart(n::Int)
    n < 33 && return ZZ(SwingOddpart[n+1])

    sqrtn = isqrt(n)
    factors = Primes(div(n,2) + 1, n)
    P = Primes(sqrtn + 1, div(n, 3))
    s = [p for p in P if isodd(div(n, p))]

    for prime in Primes(3, sqrtn)
        p, q = 1, n
        while true
            q = div(q, prime)
            q == 0 && break
            isodd(q) && (p *= prime)
        end
        p > 1 && push!(s, p)
    end

    return ∏(factors)*∏(s)
end

"""
Computes the swinging factorial (a.k.a. Swing numbers n≀) (cf. A056040).
"""
function Swing(n::Int)
    sh = count_ones(div(n, 2))
    swing_oddpart(n) << sh
end

const FactorialOddPart = [1, 1, 1, 3, 3, 15, 45, 315, 315, 2835, 14175, 155925,
    467775, 6081075, 42567525, 638512875, 638512875, 10854718875, 97692469875,
    1856156927625, 9280784638125, 194896477400625, 2143861251406875,
    49308808782358125, 147926426347074375, 3698160658676859375]

"""
Return the largest odd divisor of ``n!``. Cf. A049606.
"""
function factorial_oddpart(n::Int)
    n < length(FactorialOddPart) && return ZZ(FactorialOddPart[n+1])
    swing_oddpart(n)*(factorial_oddpart(div(n,2))^2)
end

"""
Return the factorial ``n! = 1×2× ... ×n``, which is the order of the symmetric group S_n or the number of permutations of n letters (cf. A000142).
"""
function PSfactorial(n::Int)
    n < 0 && ArgumentError("Argument must be >= 0")
    sh = n - count_ones(n)
    factorial_oddpart(n) << sh
end

#START-TEST-########################################################

using Test, BenchmarkTools

function test()
    @testset "PrimeSwingF" begin
        for n in 0:999
            S = PSfactorial(n)
            B = Base.factorial(BigInt(n))
            #@test S == B
            @test all(S .== B)
        end
    end
end

function demo()
end

"""
10^1  128.016 ns (3        allocations:     48 bytes)
10^2    8.148 μs (78       allocations:   3.31 KiB)
10^3   68.933 μs (469      allocations:  16.53 KiB)
10^4  832.578 μs (3780     allocations: 106.00 KiB)
10^5   17.717 ms (31657    allocations: 816.70 KiB)
10^6  288.445 ms (263237   allocations:   6.53 MiB)
10^7    5.484  s (2209289  allocations:  53.95 MiB)
10^8   86.723  s (18926351 allocations: 446.32 MiB)
10^9 1218.261  s (167.46 M allocations: 3.828 GiB)
"""
function superperf()
    GC.gc()
    println(10)
    @btime PSfactorial(10)
    GC.gc()
    println(100)
    @btime PSfactorial(100)
    GC.gc()
    println(1000)
    @btime PSfactorial(1000)
    GC.gc()
    println(10000)
    @btime PSfactorial(10000)
    GC.gc()
    println(100000)
    @btime PSfactorial(100000)
    GC.gc()
    println(1000000)
    @btime PSfactorial(1000000)
    GC.gc()
    println(10000000)
    @btime PSfactorial(10000000)
    GC.gc()
    println(100000000)
    @btime PSfactorial(100000000)
    GC.gc()
    println(1000000000)
    @time PSfactorial(1000000000)
end

function perf()
    GC.gc()
    @time Base.factorial(BigInt(1000000))
    @time PSfactorial(1000000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
mutable struct DeepQPolicy <: QPolicy
    primaryNetwork      # Neural Network for Deep Q function
end

get_QValues(policy::DeepQPolicy, inputs) = policy.primaryNetwork(inputs)

get_target(policy::DeepQPolicy, γ, r, done, s′) = dropdims(r .+ γ.*(1.0 .- done).*maximum(policy.primaryNetwork(s′); dims=1); dims=1)

get_params(policy::DeepQPolicy) = Flux.params(policy.primaryNetwork)

action(policy::MooreDDQN.DeepQPolicy, r, s::S, A) where {S} = argmax(policy.primaryNetwork(s()))[1]


text:
struct DijkstraHeapEntry{T,U<:Integer}
    vertex::U
    dist::T
end

isless(e1::DijkstraHeapEntry, e2::DijkstraHeapEntry) = e1.dist < e2.dist

"""
    struct DijkstraState{T, U}

An [`AbstractPathState`](@ref) designed for Dijkstra shortest-paths calculations.
"""
struct DijkstraState{T,U<:Integer} <: AbstractPathState
    parents::Vector{U}
    dists::Vector{T}
    predecessors::Vector{Vector{U}}
    pathcounts::Vector{U}
    closest_vertices::Vector{U}
end

"""
    dijkstra_shortest_paths(g, srcs, distmx=weights(g));

Perform [Dijkstra's algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)
on a graph, computing shortest distances between `srcs` and all other vertices.
Return a [`DijkstraState`](@ref) that contains various traversal information.

### Optional Arguments
- `allpaths=false`: If true, returns a [`DijkstraState`](@ref) that keeps track of all
predecessors of a given vertex.
"""
function dijkstra_shortest_paths(
    g::AbstractGraph,
    srcs::Vector{U},
    distmx::AbstractMatrix{T}=weights(g);
    allpaths=false,
    trackvertices=false
    ) where T where U<:Integer

    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    preds = fill(Vector{U}(), nvg)
    visited = zeros(Bool, nvg)
    pathcounts = zeros(Int, nvg)
    H = Vector{DijkstraHeapEntry{T,U}}()  # this should be Vector{T}() in 0.4, I think.
    dists[srcs] = zero(T)
    pathcounts[srcs] = 1
    closest_vertices = Vector{U}()  # Maintains vertices in order of distances from source

    sizehint!(H, nvg)
    sizehint!(closest_vertices, nvg)

    for v in srcs
        heappush!(H, DijkstraHeapEntry{T,U}(v, dists[v]))
        visited[v] = true
    end

    while !isempty(H)
        hentry = heappop!(H)
            # info("Popped H - got $(hentry.vertex)")
        u = hentry.vertex

        if trackvertices
          push!(closest_vertices, u)
        end

        for v in out_neighbors(g, u)
            alt = (dists[u] == typemax(T)) ? typemax(T) : dists[u] + distmx[u, v]

            if !visited[v]
                dists[v] = alt
                parents[v] = u
                pathcounts[v] += pathcounts[u]
                visited[v] = true
                if allpaths
                    preds[v] = [u;]
                end
                heappush!(H, DijkstraHeapEntry{T,U}(v, alt))
                # info("Pushed $v")
            else
                if alt < dists[v]
                    dists[v] = alt
                    parents[v] = u
                    #615
                    pathcounts[v] = 0
                    preds[v] = []
                    heappush!(H, DijkstraHeapEntry{T,U}(v, alt))
                end
                if alt == dists[v]
                    pathcounts[v] += pathcounts[u]
                    if allpaths
                        push!(preds[v], u)
                    end
                end
            end
        end
    end

    if trackvertices
      for s in vertices(g)
        if !visited[s]
          push!(closest_vertices, s)
        end
      end
    end

    pathcounts[srcs] = 1
    parents[srcs] = 0
    for src in srcs
        preds[src] = []
    end

    return DijkstraState{T,U}(parents, dists, preds, pathcounts, closest_vertices)
end

dijkstra_shortest_paths(g::AbstractGraph, src::Integer, distmx::AbstractMatrix = weights(g); allpaths=false, trackvertices=false) =
dijkstra_shortest_paths(g, [src;], distmx; allpaths=allpaths, trackvertices=trackvertices)

"""
    struct MultipleDijkstraState{T, U}

An [`AbstractPathState`](@ref) designed for multisource_dijkstra_shortest_paths calculation.
"""
struct MultipleDijkstraState{T,U<:Integer} <: AbstractPathState
    dists::Matrix{T}
    parents::Matrix{U}
end

@doc_str """
    parallel_multisource_dijkstra_shortest_paths(g, sources=vertices(g), distmx=weights(g))

Compute the shortest paths between all pairs of vertices in graph `g` by running
[`dijkstra_shortest_paths`] for every vertex and using an optional list of source vertex `sources` and
an optional distance matrix `distmx`. Return a [`MultipleDijkstraState`](@ref) with relevant
traversal information.
"""

function parallel_multisource_dijkstra_shortest_paths(
    g::AbstractGraph,
    sources::AbstractVector = vertices(g),
    distmx::AbstractMatrix{T} = weights(g)
    ) where T

    U = eltype(g)
    n_v = nv(g)
    r_v = length(sources)

    dists   = SharedArray(zeros(T, r_v, n_v))
    parents = SharedArray(zeros(U, r_v, n_v))

    @sync @parallel for i in 1:r_v
      state = dijkstra_shortest_paths(g, sources[i], distmx)
      dists[i, :] = state.dists
      parents[i, :] = state.parents
    end

    result = MultipleDijkstraState(Matrix(dists), Matrix(parents))
    return result
end


text:
const tile_fg_names = Dict{Any, Any}(
    '0' => "Air",
    '1' => "Dirt",
    '3' => "Snow",
    '4' => "Girder",
    '5' => "Tower",
    '6' => "Stone",
    '7' => "Cement",
    '8' => "Rock",
    '9' => "Wood",
    'a' => "Wood Stone",
    'b' => "Cliffside",
    'c' => "Pool Edges",
    'd' => "Temple A",
    'e' => "Temple B",
    'f' => "Cliffside Alt",
    'g' => "Reflection",
    'G' => "Reflection Alt",
    'h' => "Grass",
    'i' => "Summit",
    'j' => "Summit No Snow",
    'k' => "Core",
    'l' => "Deadgrass",
    'm' => "Lost Levels",
    'n' => "Scifi",
    'z' => "Template",

    "Air" => '0',
    "Dirt" => '1',
    "Snow" => '3',
    "Girder" => '4',
    "Tower" => '5',
    "Stone" => '6',
    "Cement" => '7',
    "Rock" => '8',
    "Wood" => '9',
    "Wood Stone" => 'a',
    "Cliffside" => 'b',
    "Pool Edges" => 'c',
    "Temple A" => 'd',
    "Temple B" => 'e',
    "Cliffside Alt" => 'f',
    "Reflection" => 'g',
    "Reflection Alt" => 'G',
    "Grass" => 'h',
    "Summit" => 'i',
    "Summit No Snow" => 'j',
    "Core" => 'k',
    "Deadgrass" => 'l',
    "Lost Levels" => 'm',
    "Scifi" => 'n',
    "Template" => 'z'
)

const tile_bg_names = Dict{Any, Any}(
    '0' => "Air",
    '1' => "Dirt",
    '2' => "Brick",
    '3' => "Brick Ruined",
    '4' => "Wood",
    '5' => "Resort Stone",
    '6' => "Cliffside",
    '7' => "Pool",
    '8' => "Temple A",
    '9' => "Temple B",
    'a' => "Reflection",
    'b' => "Snow",
    'c' => "Summit",
    'd' => "Core",
    'e' => "Lost Levels",
    'z' => "Template",

    "Air" => '0',
    "Dirt" => '1',
    "Brick" => '2',
    "Brick Ruined" => '3',
    "Wood" => '4',
    "Resort Stone" => '5',
    "Cliffside" => '6',
    "Pool" => '7',
    "Temple A" => '8',
    "Temple B" => '9',
    "Reflection" => 'a',
    "Snow" => 'b',
    "Summit" => 'c',
    "Core" => 'd',
    "Lost Levels" => 'e',
    "Template" => 'z'
)

const valid_fg_tiles = Char[
    '0', '1', '3', '4', '5', '6', '7',
    '8', '9', 'a', 'b', 'c', 'd', 'e',
    'f', 'g', 'G', 'h', 'i', 'j', 'k',
    'l'
]

const valid_bg_tiles = Char[
    '0', '1', '2', '3', '4', '5', '6',
    '7', '8', '9', 'a', 'b', 'c', 'd'
]

@valueequals mutable struct Tiles
    data::Array{Char, 2}

    Tiles() = new(fill('0', (40, 23)))

    Tiles(d::Array{Char, 2}) = new(d)

    function Tiles(s::String)
        s = chomp(replace(s, "\r\n" => "\n"))
        lines = split(s, "\n")

        cols = maximum(length.(lines))
        rows = length(lines)

        res = fill('0', rows, cols)

        for (i, line) in enumerate(lines)
            res[i, 1:length(line)] = collect(line)
        end

        return new(res)
    end
end

# Removes illegal characters from the tileset
# Makes it possible to put simple entity map and fg in same string
function FgTiles(tiles::Tiles, valid::Array{Char, 1}=valid_fg_tiles)
    tiles.data = [c in valid ? c : '0' for c in tiles.data]

    return tiles
end

FgTiles(tiles::Array{Char, 2}, valid::Array{Char, 1}=valid_fg_tiles) = FgTiles(Tiles(tiles), valid)
FgTiles(s::String, valid::Array{Char, 1}=valid_fg_tiles) = FgTiles(Tiles(s), valid)

BgTiles(tiles::Tiles, valid::Array{Char, 1}=valid_bg_tiles) = FgTiles(tiles, valid)
BgTiles(tiles::Array{Char, 2}, valid::Array{Char, 1}=valid_bg_tiles) = FgTiles(Tiles(tiles), valid)
BgTiles(s::String, valid::Array{Char, 1}=valid_bg_tiles) = FgTiles(Tiles(s), valid)

function getRelevantColumnCount(data::Matrix{T}, emptyValue::T) where T
    rows, cols = size(data)
    columns = fill(0, rows)

    for y in 1:rows
        columns[y] = something(findlast(n -> n != emptyValue, view(data, y, 1:cols)), 0)
    end

    return columns
end

function getRelevantRowCount(columnLengths::Array{Int, 1})
    return something(findlast(n -> n > 0, columnLengths), 0)
end

# All trailing '0's don't need to be written to be loadable
# We can also stop when all remaining rows are '0's
function minimizeTilesString(data::Matrix{T}, emptyValue::T, separator::Union{Char, String}=',') where T
    columnLengths = getRelevantColumnCount(data, emptyValue)
    relevantRows = getRelevantRowCount(columnLengths)

    lines = Array{String, 1}(undef, relevantRows)

    for row in 1:relevantRows
        # String on Char[] is faster than joining with empty string
        if separator == ""
            lines[row] = String(view(data, row, 1:columnLengths[row]))

        else
            lines[row] = join(view(data, row, 1:columnLengths[row]), separator)
        end
    end

    return join(lines, '\n')
end

function Base.string(t::Tiles)
    return join(String.([t.data[i, :] for i in 1:size(t.data, 1)]), "\n")
end

Base.size(t::Tiles) = reverse(size(t.data)) .* 8

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
"""
$(TYPEDEF)
"""
abstract type AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractScreeningDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractFactorialDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractOptimalDesign <: AbstractDesign end

"""
$(TYPEDEF)
"""
abstract type AbstractRandomDesign <: AbstractDesign end

"""
$(TYPEDEF)

Encapsulates  a  Plackett-Burman  screening  design  constructed  using  Paley's
method.  Factor  levels are  encoded as  `:high` and  `:low` symbols,  and extra
dummy variables, possibly  aliasing interactions between factors,  will be added
to pad a design.

$(TYPEDFIELDS)
"""
struct PlackettBurman <: AbstractScreeningDesign
    matrix::DataFrame
    factors::Tuple
    dummy_factors::Tuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(@formula(y ~ x1 + x2 + x3 + x4))
PlackettBurman(8×7 DataFrames.DataFrame
│ Row │ x1     │ x2     │ x3     │ x4     │ dummy1 │ dummy2 │ dummy3 │
│     │ Symbol │ Symbol │ Symbol │ Symbol │ Symbol │ Symbol │ Symbol │
├─────┼────────┼────────┼────────┼────────┼────────┼────────┼────────┤
│ 1   │ high   │ high   │ high   │ high   │ high   │ high   │ high   │
│ 2   │ low    │ high   │ low    │ high   │ high   │ low    │ low    │
│ 3   │ high   │ low    │ high   │ high   │ low    │ low    │ low    │
│ 4   │ low    │ high   │ high   │ low    │ low    │ low    │ high   │
│ 5   │ high   │ high   │ low    │ low    │ low    │ high   │ low    │
│ 6   │ high   │ low    │ low    │ low    │ high   │ low    │ high   │
│ 7   │ low    │ low    │ low    │ high   │ low    │ high   │ high   │
│ 8   │ low    │ low    │ high   │ low    │ high   │ high   │ low    │, (:x1, :x2, :x3, :x4), (:dummy1, :dummy2, :dummy3), y ~ x1 + x2 + x3 + x4 + dummy1 + dummy2 + dummy3)

```
"""
function PlackettBurman(formula::FormulaTerm)
    symbol_factors = Tuple(r.sym for r in formula.rhs)

    initial_design = plackettburman(length(symbol_factors))
    categorical_design = similar(initial_design, Symbol)

    for i = 1:size(initial_design, 1), j = 1:size(initial_design, 2)
        categorical_design[i, j] = initial_design[i, j] == 1.0 ? :high : :low
    end

    design = DataFrame(categorical_design)

    dummy_factors = Tuple(Symbol("dummy" * string(i)) for i = 1:(length(names(design)) -
                                                                 length(symbol_factors)))

    design_names = (symbol_factors..., dummy_factors...)
    names!(design, collect(design_names))

    PlackettBurman(design,
                   symbol_factors,
                   dummy_factors,
                   term(formula.lhs) ~ sum(term.(design_names)))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> PlackettBurman(4)
PlackettBurman(8×7 DataFrames.DataFrame
│ Row │ factor1 │ factor2 │ factor3 │ factor4 │ dummy1 │ dummy2 │ dummy3 │
│     │ Symbol  │ Symbol  │ Symbol  │ Symbol  │ Symbol │ Symbol │ Symbol │
├─────┼─────────┼─────────┼─────────┼─────────┼────────┼────────┼────────┤
│ 1   │ high    │ high    │ high    │ high    │ high   │ high   │ high   │
│ 2   │ low     │ high    │ low     │ high    │ high   │ low    │ low    │
│ 3   │ high    │ low     │ high    │ high    │ low    │ low    │ low    │
│ 4   │ low     │ high    │ high    │ low     │ low    │ low    │ high   │
│ 5   │ high    │ high    │ low     │ low     │ low    │ high   │ low    │
│ 6   │ high    │ low     │ low     │ low     │ high   │ low    │ high   │
│ 7   │ low     │ low     │ low     │ high    │ low    │ high   │ high   │
│ 8   │ low     │ low     │ high    │ low     │ high   │ high   │ low    │, (:factor1, :factor2, :factor3, :factor4), (:dummy1, :dummy2, :dummy3), response ~ factor1 + factor2 + factor3 + factor4 + dummy1 + dummy2 + dummy3)

```
"""
function PlackettBurman(factors::Int)
    PlackettBurman(term(:response) ~ sum(term.(Symbol("factor" * string(i)) for i = 1:factors)))
end

"""
$(TYPEDEF)

Encapsulates a full  factorial design. Can contain an explicit  design where all
experiments are completely computed and stored, or an iterator of a design, from
which experiments may be extracted as needed.

$(TYPEDFIELDS)
"""
struct FullFactorial <: AbstractFactorialDesign
    matrix::Union{DataFrame, Missing}
    iterator::Base.Iterators.ProductIterator
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]), @formula(y ~ A + B +C), explicit = true)
FullFactorial(12×3 DataFrames.DataFrame
│ Row │ A   │ B   │ C    │
│     │ Any │ Any │ Any  │
├─────┼─────┼─────┼──────┤
│ 1   │ 1   │ a   │ 1.0  │
│ 2   │ 2   │ a   │ 1.0  │
│ 3   │ 4   │ a   │ 1.0  │
│ 4   │ 1   │ b   │ 1.0  │
│ 5   │ 2   │ b   │ 1.0  │
│ 6   │ 4   │ b   │ 1.0  │
│ 7   │ 1   │ a   │ -1.0 │
│ 8   │ 2   │ a   │ -1.0 │
│ 9   │ 4   │ a   │ -1.0 │
│ 10  │ 1   │ b   │ -1.0 │
│ 11  │ 2   │ b   │ -1.0 │
│ 12  │ 4   │ b   │ -1.0 │, Base.Iterators.ProductIterator{Tuple{Array{Int64,1},Array{Symbol,1},Array{Float64,1}}}(([1, 2, 4], Symbol[:a, :b], [1.0, -1.0])), (A = [1, 2, 4], B = Symbol[:a, :b], C = [1.0, -1.0]), y ~ A + B + C)

```
"""
function FullFactorial(factors::NamedTuple, formula::FormulaTerm; explicit::Bool = false)
    iterator = fullfactorial(values(factors))
    matrix = missing

    if explicit
        matrix = DataFrame(explicit_fullfactorial(iterator))
        names!(matrix, [r.sym for r in formula.rhs])
    end

    FullFactorial(matrix, iterator, factors, formula)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial((A = [1, 2, 4], B = [:a, :b], C = [1.0, -1.0]), explicit = true)
FullFactorial(12×3 DataFrames.DataFrame
│ Row │ A   │ B   │ C    │
│     │ Any │ Any │ Any  │
├─────┼─────┼─────┼──────┤
│ 1   │ 1   │ a   │ 1.0  │
│ 2   │ 2   │ a   │ 1.0  │
│ 3   │ 4   │ a   │ 1.0  │
│ 4   │ 1   │ b   │ 1.0  │
│ 5   │ 2   │ b   │ 1.0  │
│ 6   │ 4   │ b   │ 1.0  │
│ 7   │ 1   │ a   │ -1.0 │
│ 8   │ 2   │ a   │ -1.0 │
│ 9   │ 4   │ a   │ -1.0 │
│ 10  │ 1   │ b   │ -1.0 │
│ 11  │ 2   │ b   │ -1.0 │
│ 12  │ 4   │ b   │ -1.0 │, Base.Iterators.ProductIterator{Tuple{Array{Int64,1},Array{Symbol,1},Array{Float64,1}}}(([1, 2, 4], Symbol[:a, :b], [1.0, -1.0])), (A = [1, 2, 4], B = Symbol[:a, :b], C = [1.0, -1.0]), response ~ A + B + C)

```
"""
function FullFactorial(factors::NamedTuple; explicit::Bool = false)
    FullFactorial(factors, term(:response) ~ sum(term.(keys(factors))), explicit = explicit)
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> FullFactorial(([1, 2, 4], [:a, :b], [1.0, -1.0]), explicit = true)
FullFactorial(12×3 DataFrames.DataFrame
│ Row │ factor1 │ factor2 │ factor3 │
│     │ Any     │ Any     │ Any     │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ 1       │ a       │ 1.0     │
│ 2   │ 2       │ a       │ 1.0     │
│ 3   │ 4       │ a       │ 1.0     │
│ 4   │ 1       │ b       │ 1.0     │
│ 5   │ 2       │ b       │ 1.0     │
│ 6   │ 4       │ b       │ 1.0     │
│ 7   │ 1       │ a       │ -1.0    │
│ 8   │ 2       │ a       │ -1.0    │
│ 9   │ 4       │ a       │ -1.0    │
│ 10  │ 1       │ b       │ -1.0    │
│ 11  │ 2       │ b       │ -1.0    │
│ 12  │ 4       │ b       │ -1.0    │, Base.Iterators.ProductIterator{Tuple{Array{Int64,1},Array{Symbol,1},Array{Float64,1}}}(([1, 2, 4], Symbol[:a, :b], [1.0, -1.0])), (factor1 = [1, 2, 4], factor2 = Symbol[:a, :b], factor3 = [1.0, -1.0]), response ~ factor1 + factor2 + factor3)

```
"""
function FullFactorial(factors::Tuple; explicit::Bool = false)
    FullFactorial(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors), explicit = explicit)
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct FractionalFactorial <: AbstractFactorialDesign
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)

Encapsulates a random design generator.  Receives a `NamedTuple` of factor names
associated  with `Distribution`s  from the  the *Distributions*  package.  After
instantiating  a  `RandomDesign`,  you  must   request  samples  from  it  using
[`rand`](@ref).

"""
struct RandomDesign <: AbstractRandomDesign
    factors::NamedTuple
    formula::FormulaTerm
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> RandomDesign((f1 = Distributions.Uniform(2, 3), f2 = Distributions.DiscreteUniform(-1, 5), f3 = Distributions.Uniform(5, 10)))
RandomDesign((f1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), f2 = Distributions.DiscreteUniform(a=-1, b=5), f3 = Distributions.Uniform{Float64}(a=5.0, b=10.0)), response ~ f1 + f2 + f3)

```
"""
function RandomDesign(factors::NamedTuple)
    RandomDesign(factors, term(:response) ~ sum(term.(keys(factors))))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> RandomDesign((Distributions.Uniform(2, 3), Distributions.DiscreteUniform(-1, 5), Distributions.Uniform(5, 10)))
RandomDesign((factor1 = Distributions.Uniform{Float64}(a=2.0, b=3.0), factor2 = Distributions.DiscreteUniform(a=-1, b=5), factor3 = Distributions.Uniform{Float64}(a=5.0, b=10.0)), response ~ factor1 + factor2 + factor3)

```
"""
function RandomDesign(factors::Tuple)
    RandomDesign(NamedTuple{Tuple(Symbol("factor" * string(i)) for i = 1:length(factors))}(factors))
end

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> rand(RandomDesign((f1 = Distributions.Uniform(2, 3), f2 = Distributions.DiscreteUniform(-1, 5), f3 = Distributions.Uniform(5, 10))), 12)
12×3 DataFrames.DataFrame
│ Row │ f1      │ f2   │ f3      │
│     │ Real    │ Real │ Real    │
├─────┼─────────┼──────┼─────────┤
│ 1   │ 2.04922 │ 5    │ 8.85741 │
│ 2   │ 2.25659 │ -1   │ 6.57617 │
│ 3   │ 2.86526 │ 4    │ 5.41422 │
│ 4   │ 2.81201 │ -1   │ 5.40944 │
│ 5   │ 2.92412 │ 1    │ 6.22902 │
│ 6   │ 2.20051 │ 5    │ 8.80749 │
│ 7   │ 2.69459 │ 3    │ 6.34626 │
│ 8   │ 2.28902 │ 2    │ 8.72759 │
│ 9   │ 2.95173 │ 0    │ 5.42517 │
│ 10  │ 2.35163 │ 4    │ 5.89413 │
│ 11  │ 2.97614 │ 5    │ 9.99911 │
│ 12  │ 2.4577  │ 0    │ 8.303   │

```
"""
function rand(design::RandomDesign, n::Int = 1)
    DataFrame(random_design(values(design.factors), n), collect(keys(design.factors)))
end

"""
$(TYPEDEF)

$(TYPEDFIELDS)
"""
struct OptimalDesign <: AbstractOptimalDesign
    matrix::DataFrame
    candidates::AbstractDesign
    criteria::Dict{Symbol, Float64}
    formula::FormulaTerm
end


text:
module RocketTakeUntilOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: take_until()" begin

    println("Testing: operator take_until()")

    run_proxyshowcheck("TakeUntil", take_until(never()), args = (check_subscription = true, ))

    run_testset([
        (
            source = from(1:5) |> take_until(of(1)),
            values = @ts([ 1:5, c ])
        ),
        (
            source = from(1:5) |> async(0) |> take_until(of(1)),
            values = @ts(c)
        ),
        (
            source = timer(10, 1000) |> take_until(timer(100)),
            values = @ts([ 0 ] ~ c)
        ),
        (
            source = completed() |> take_until(of(1)),
            values = @ts(c)
        ),
        (
            source = completed() |> async(0) |> take_until(of(1)),
            values = @ts(c)
        ),
        (
            source = completed() |> take_until(timer(100)),
            values = @ts(c)
        ),
        (
            source = completed() |> async(0) |> take_until(timer(100)),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> take_until(of(1)),
            values = @ts(e(1))
        ),
        (
            source = faulted(1) |> async(0) |> take_until(of(1)),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> take_until(timer(100)),
            values = @ts(e(1))
        ),
        (
            source = faulted(1) |> async(0) |> take_until(timer(100)),
            values = @ts(e(1))
        ),
        (
            source = never() |> take_until(of(1)),
            values = @ts(c)
        ),
        (
            source = never() |> take_until(timer(100)),
            values = @ts(100 ~ c)
        )
    ])

end

end


text:
using MultiJuMP, JuMP
using Test
#using AmplNLWriter
using Ipopt
using Clp

@testset "Utopia and Nadir points" begin
    mdata = MultiJuMP.MultiData()
    push!(mdata.objectives, SingleObjective(;sense = MOI.MIN_SENSE))
    push!(mdata.objectives, SingleObjective(;sense = MOI.MAX_SENSE))
    mdata.Phi .= [1 2; 3 -1]
    nadir = [2.0, -3.0]
    utopia = [1.0, 1.0]

    @test getnadir(mdata) == nadir
    @test getutopia(mdata) == utopia
end

@testset "NBI optimisation" begin
    m = multi_model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    @variable(m, x[i=1:5])
    @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = get_multidata(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    optimize!(m, method = MultiJuMP.NBI(false))

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0    , 7.16982 , 4.48665 , 2.08008]
    f2true = [2.13057 , -4.01115, -2.78066, -1.45458, 0.0513689]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "WS optimisation" begin
    m = multi_model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    @variable(m, x[i=1:5])
    @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = get_multidata(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    optimize!(m, method = MultiJuMP.WeightedSum())

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0, 10.0, 2.88267, 0.726877]
    f2true = [2.13057, -4.01115, -4.01115, -0.510963, 1.48613]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "EPS optimisation" begin
    m = multi_model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
    @variable(m, x[i=1:5])
        @NLexpressions m begin
        f1, sum(x[i]^2 for i=1:5)
        f2, 3x[1]+2x[2]-x[3]/3+0.01*(x[4]-x[5])^3
    end
    @NLconstraints m begin
        x[1]+2x[2]-x[3]-0.5x[4]+x[5]==2
        4x[1]-2x[2]+0.8x[3]+0.6x[4]+0.5x[5]^2 == 0
        sum(x[i]^2 for i=1:5) <= 10
    end

    obj1 = SingleObjective(f1)
    obj2 = SingleObjective(f2)

    multim = get_multidata(m)
    multim.objectives = [obj1, obj2]
    multim.pointsperdim = 5
    optimize!(m, method = MultiJuMP.EpsilonCons())

    f1arr = convert(Array{Float64},
                    [multim.paretofront[i][1] for i in 1:multim.pointsperdim])
    f2arr = convert(Array{Float64},
                    [multim.paretofront[i][2] for i in 1:multim.pointsperdim])
    f1true = [0.555081, 10.0, 7.63877, 5.27754, 2.91631]
    f2true = [2.13057, -4.01115, -2.99356, -1.86918, -0.532781]

    @test isapprox(f1arr, f1true, atol=1e-2)
    @test isapprox(f2arr, f2true, atol=1e-2)
end

@testset "WS linear" begin
    mmodel = multi_model(Clp.Optimizer, linear = true)
    y = @variable(mmodel, 0 <= y <= 10.0)
    z = @variable(mmodel, 0 <= z <= 10.0)
    @constraint(mmodel, y + z <= 15.0)

    # objectives
    exp_obj1 = @expression(mmodel, -y +0.05 * z)
    exp_obj2 = @expression(mmodel, 0.05 * y - z)
    obj1 = SingleObjective(exp_obj1)
    obj2 = SingleObjective(exp_obj2)

    # setting objectives in the MultiData
    multim = get_multidata(mmodel)
    multim.objectives = [obj1, obj2]

    optimize!(mmodel, method = MultiJuMP.WeightedSum())
    @test termination_status(mmodel) == MOI.OPTIMAL

    true_par_vals = [-10.0, -9.75, -4.5, 0.5]
    for idx in eachindex(true_par_vals)
        tf1 = true_par_vals[idx]
        tf2 = true_par_vals[4 - idx + 1]
        @test any(multim.paretofront) do v
            v[1] ≈ tf1 && v[2] ≈ tf2
        end
    end
end

@testset "EPS linear" begin
    mmodel = multi_model(Clp.Optimizer, linear = true)
    y = @variable(mmodel, 0 <= y <= 10.0)
    z = @variable(mmodel, 0 <= z <= 10.0)
    @constraint(mmodel, y + z <= 15.0)

    # objectives
    exp_obj1 = @expression(mmodel, -y +0.05 * z)
    exp_obj2 = @expression(mmodel, 0.05 * y - z)
    obj1 = SingleObjective(exp_obj1)
    obj2 = SingleObjective(exp_obj2)

    # setting objectives in the MultiData
    multim = get_multidata(mmodel)
    multim.objectives = [obj1, obj2]

    optimize!(mmodel, method = MultiJuMP.EpsilonCons())
    @test termination_status(mmodel) == MOI.OPTIMAL
    true_par_pos = begin
        v = [10.0, 10.0, 5.0, 0.0]
        zip(v, reverse(v)) |> collect
    end

    for (yv, zv) in true_par_pos
        @test any(multim.paretovarvalues) do v
            ≈(v[1], yv, atol = 0.01) && ≈(v[2], zv, atol = 0.01)
        end
    end
end

@testset "WS linear - alternative model construction" begin
    m = Model(Clp.Optimizer)
    mmodel = multi_model(m, linear = true)
    y = @variable(mmodel, 0 <= y <= 10.0)
    z = @variable(mmodel, 0 <= z <= 10.0)
    @constraint(mmodel, y + z <= 15.0)

    # objectives
    exp_obj1 = @expression(mmodel, -y +0.05 * z)
    exp_obj2 = @expression(mmodel, 0.05 * y - z)
    obj1 = SingleObjective(exp_obj1)
    obj2 = SingleObjective(exp_obj2)

    # setting objectives in the MultiData
    multim = get_multidata(mmodel)
    multim.objectives = [obj1, obj2]

    optimize!(mmodel, method = MultiJuMP.WeightedSum())
    @test termination_status(mmodel) == MOI.OPTIMAL

    true_par_vals = [-10.0, -9.75, -4.5, 0.5]
    for idx in eachindex(true_par_vals)
        tf1 = true_par_vals[idx]
        tf2 = true_par_vals[4 - idx + 1]
        @test any(multim.paretofront) do v
            v[1] ≈ tf1 && v[2] ≈ tf2
        end
    end
end

@testset "TODO: Initial value test" begin
    # TODO: Let f1 have two minima and
    # check that initial value gives correct minima
    # TODO: Add anonymous variables
    @test 1 == 1
end

@testset "TODO: NBI Inequality test" begin
    # TODO: Test that the inequality-extension of NBI works
    # We can use the nc_paper_2.jl example
    @test 1 == 1
end


text:
type SparseMatrix{T} <: AbstractMatrix{T}
	obj::Ptr{Void}
end


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

## Observed Data
three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")
three_path_par = Feather.read("test/comparisons/three_path_par.feather")
three_path_start = Feather.read("test/comparisons/three_path_start.feather")

semobserved = SemObsCommon(data = Matrix(three_path_dat))

diff_fin = SemFiniteDiff(BFGS(), Optim.Options())

## Model definition
@variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

loss = Loss([SemML(semobserved, [0.0], similar(start_val))])

imply = ImplySymbolic(A, S, F, x, start_val)

model_fin = Sem(semobserved, imply, loss, diff_fin)

A = Vector{Any}(undef, 6)
for i in 1:6
    global A[i] = deepcopy(imply)
end

##

rowind = [1:25, collect(26:50), collect(51:75)]

differ_group = [
    fill(1, 31),
    vcat(fill(1, 14), fill(2, 17)),
    vcat(fill(1, 7), fill(2, 7), fill(3, 10), fill(2, 7))
    ]

function parsubset(differ_group, start_val)
    differ_group = hcat(differ_group...)
    npar, ngroup = size(differ_group)

    parunique = Vector{Int64}(undef, npar)
    for i in 1:npar
        parunique[i] = length(unique(differ_group[i, :]))
    end

    npar_effective = sum(parunique)

    locations = zeros(Int64, ngroup, npar)
    for i in 1:ngroup
        locations[i, :] .=
            cumsum(parunique) .- parunique .+ differ_group[:, i]
    end

    start_val_long = vcat(fill.(start_val, parunique)...)

    parsubsets = zeros(Bool, ngroup, npar_effective)
    for i in 1:ngroup
        parsubsets[i, locations[i, :]] .= true
    end

    return parsubsets, start_val_long
end

parsubset(differ_group, start_val)

function make_onelement_array(A)
    isa(A, Array) ? nothing : (A = [A])
    return A
end

function semvec(observed, imply, loss, diff)

    observed = make_onelement_array(observed)
    imply = make_onelement_array(imply)
    loss = make_onelement_array(loss)
    diff = make_onelement_array(diff)

    sem_vec = Array{AbstractSem}(undef, maximum(length.([observed, imply, loss, diff])))
    sem_vec .= Sem.(observed, imply, loss, diff)

    return sem_vec
end

function get_observed(rowind, data, semobserved;
            args = (),
            kwargs = NamedTuple())
    observed_vec = Vector{semobserved}(undef, length(rowind))
    for i in 1:length(rowind)
        observed_vec[i] = semobserved(
                            args...;
                            data = Matrix(data[rowind[i], :]),
                            kwargs...)
    end
    return observed_vec
end

function start_val(multigroup_sem, start_val)
end

function SemMG(
    data, start_val, differ_group, rowind,
    observed, imply, loss, diff;
    obs_args = (),
    imply_agrs = (),
    loss_args = (),
    diff_args = ())

    #diff = diff(diff_args...)
    #imply = imply(imply_args...)
    #loss = loss(loss_args...)
    #observed = observed(obs_args...)

    sem_vec = semvec(observed, imply, loss, diff)

    par_subsets, start_val_long = parsubset(differ_group, start_val)

    return SemMG2(sem_vec, par_subsets), start_val_long
end

obs_list = get_observed(rowind, three_path_dat, SemObsCommon)

sem_mg, start_val_mg = SemMG(three_path_dat, start_val, differ_group, rowind, obs_list,
    imply, loss, diff_fin)

sem_mg

function (semmg::SemMG)(par)
    F = zero(eltype(par))
    for i in 1:length(semmg.sem_vec)
        F += semmg.sem_vec[i](par[sem.par_subsets[i,:]])
    end
    return F
end

start_val_mg[sem_mg.par_subsets[3,:]] == start_val

## Testing Zone

function MultigroupSem(
    data, algorithm, options, A, S, F, x, start_val,
    differ_group, rowind,
    semobserved, imply, loss, diff;
    obs_args = NamedTuple(),
    imply_agrs = NamedTuple(),
    loss_args = NamedTuple(),
    diff_args = NamedTuple())

    diff = diff(algorithm, options, )
end

function fitfun(par, models, parsubsets)
    F = zero(eltype(par))
    for i in 1:length(models)
        F += models[i](par[parsubsets[i,:]])
    end
    return F
end

function test(diff, algo; args1 = NamedTuple())
    return diff(;algo, args1...)
end

function test2(;algo = nothing, data = nothing)
    return data
end

test(test2, "a")

Dict([("data", 4)])

x = (data = 10,)

x

function test3(;data)
    return data
end

test3(;x...)

maximum(length.([1 2 3]))

length(semobserved)

function f(constructor, value, length)
    vec = Vector{constructor}(undef, length)
    for i = 1:length
        vec[i] = constructor(value)
    end
end

struct SemMG2{V <: Tuple{AbstractSem}, D <: Array{Bool, 2}} <: AbstractSem
    sem_vec::V
    par_subsets::D
end


text:
let
    deg = 45.0
    rad = deg*pi/180
    r = Rx(deg, use_degrees=true)

    tol = 1e-8
    @test isapprox(r[1, 1], 1.0,       atol=tol)
    @test isapprox(r[1, 2], 0.0,       atol=tol)
    @test isapprox(r[1, 3], 0.0,       atol=tol)

    @test isapprox(r[2, 1], 0.0,       atol=tol)
    @test isapprox(r[2, 2], +cos(rad), atol=tol)
    @test isapprox(r[2, 3], +sin(rad), atol=tol)

    @test isapprox(r[3, 1], 0.0,       atol=tol)
    @test isapprox(r[3, 2], -sin(rad), atol=tol)
    @test isapprox(r[3, 3], +cos(rad), atol=tol)
end

let
    deg = 45.0
    rad = deg*pi/180
    r = Ry(deg, use_degrees=true)

    tol = 1e-8
    @test isapprox(r[1, 1], +cos(rad), atol=tol)
    @test isapprox(r[1, 2], 0.0,       atol=tol)
    @test isapprox(r[1, 3], -sin(rad), atol=tol)

    @test isapprox(r[2, 1], 0.0,       atol=tol)
    @test isapprox(r[2, 2], 1.0,       atol=tol)
    @test isapprox(r[2, 3], 0.0,       atol=tol)

    @test isapprox(r[3, 1], +sin(rad), atol=tol)
    @test isapprox(r[3, 2], 0.0,       atol=tol)
    @test isapprox(r[3, 3], +cos(rad), atol=tol)
end

let
    deg = 45.0
    rad = deg*pi/180
    r = Rz(deg, use_degrees=true)

    tol = 1e-8
    @test isapprox(r[1, 1], +cos(rad), atol=tol)
    @test isapprox(r[1, 2], +sin(rad), atol=tol)
    @test isapprox(r[1, 3], 0.0,       atol=tol)

    @test isapprox(r[2, 1], -sin(rad), atol=tol)
    @test isapprox(r[2, 2], +cos(rad), atol=tol)
    @test isapprox(r[2, 3], 0.0,       atol=tol)

    @test isapprox(r[3, 1], 0.0,       atol=tol)
    @test isapprox(r[3, 2], 0.0,       atol=tol)
    @test isapprox(r[3, 3], 1.0,       atol=tol)
end

###########################
# Quaternion Constructors #
###########################

let
    # Vector initialization
    q = Quaternion([1.1, 0.0, 0.0, 0.0])

    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # Test vector initialization errors
    @test_throws ArgumentError Quaternion([1 2 3 4 5])
end

let
    # Matrix initialization Quaternion Fomulation 1
    mat = [1.0 0.0 0.0;
           0.0 1.0 0.0;
           0.0 0.0 1.0]

    q = Quaternion(mat)

    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # Matrix initialization Quaternion Fomulation 2
    mat = [1.0 0.0 0.0;
           0.0 -1.0 0.0;
           0.0 0.0 -1.0]

    q = Quaternion(mat)

    tol = 1e-12
    @test isapprox(q.q0, 0.0, atol=tol)
    @test isapprox(q.q1, 1.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # Matrix initialization Quaternion Fomulation 3
    mat = [-1.0 0.0 0.0;
           0.0 1.0 0.0;
           0.0 0.0 -1.0]

    q = Quaternion(mat)

    tol = 1e-12
    @test isapprox(q.q0, 0.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 1.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # Matrix initialization Quaternion Fomulation 4
    mat = [-1.0 0.0 0.0;
           0.0 -1.0 0.0;
           0.0 0.0 1.0]

    q = Quaternion(mat)

    tol = 1e-12
    @test isapprox(q.q0, 0.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 1.0, atol=tol)
end

let
    # EulerAngle initialization - 121
    e = EulerAngle(121, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 123
    e = EulerAngle(123, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 131
    e = EulerAngle(131, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 132
    e = EulerAngle(132, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 212
    e = EulerAngle(212, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 213
    e = EulerAngle(213, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 231
    e = EulerAngle(231, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 232
    e = EulerAngle(232, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 312
    e = EulerAngle(312, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 313
    e = EulerAngle(313, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 321
    e = EulerAngle(321, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)

    # EulerAngle initialization - 323
    e = EulerAngle(323, 0.0, 0.0, 0.0)

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)
end

let
    # EulerAngle initialization - 323
    e = EulerAngle(323, 0.0, 0.0, 0.0)

    # Change to 999 to confirm error
    e.seq = 999

    @test typeof(e) == EulerAngle
    @test_throws ArgumentError Quaternion(e)
end

let
    # EulerAxis initialization
    e = EulerAxis(0.0, [0.0, 0.0, 0.0])

    q = Quaternion(e)
    tol = 1e-12
    @test isapprox(q.q0, 1.0, atol=tol)
    @test isapprox(q.q1, 0.0, atol=tol)
    @test isapprox(q.q2, 0.0, atol=tol)
    @test isapprox(q.q3, 0.0, atol=tol)
end

#########################
# Quaternion Operations #
#########################

let
    q = Quaternion(1.0, 1.0, 1.0, 1.0)

    v = q[:]

    tol = 1e-12
    array_isapprox(v, 0.5, atol=tol)

    # Test Access
    @test isapprox(q[1], 0.5, atol=tol)
    @test isapprox(q[2], 0.5, atol=tol)
    @test isapprox(q[3], 0.5, atol=tol)
    @test isapprox(q[4], 0.5, atol=tol)

    @test length(q[2:3]) == 2

    @test_throws BoundsError q5 = q[5]
end

let
    q = Quaternion(1.0, 1.0, 1.0, 1.0)

    v = as_vector(q)

    tol = 1e-12
    array_isapprox(v, 0.5, atol=tol)
end

let
    q = Quaternion(1.0, 0.0, 0.0, 0.0)

    mat = as_matrix(q)

    tol = 1.0e-12
    array_isapprox(mat, one(mat), atol=tol)
end

let
    q1 = Quaternion(randn(4))
    q2 = copy(q1)

    @test pointer_from_objref(q1) != pointer_from_objref(q2)
end

let
    q1 = Quaternion(randn(4))
    q2 = deepcopy(q1)

    @test pointer_from_objref(q1) != pointer_from_objref(q2)
end

let
    q = Quaternion([1.0 0.0 0.0 0.0])

    # Circumvent normalization during Quaternion initialization
    q.q0 = 1.0
    q.q1 = 1.0
    q.q2 = 1.0
    q.q3 = 1.0

    # Norm
    @test norm(q) == 2.0

    # Normalize Quaternion
    normalize(q)

    tol = 1e-12
    array_isapprox(q[:], 0.5, atol=tol)

    @test normalize(q) == nothing

    @test_throws ArgumentError Quaternion([1 2 3 4 5])
end

let
    q = Quaternion(randn(4))

    qc = conj(q)

    tol = 1e-12
    @test isapprox(qc.q0,  q.q0, atol=tol)
    @test isapprox(qc.q1, -q.q1, atol=tol)
    @test isapprox(qc.q2, -q.q2, atol=tol)
    @test isapprox(qc.q3, -q.q3, atol=tol)
end

let
    q = Quaternion(randn(4))

    qi = inv(q)

    tol = 1e-12
    @test isapprox(qi.q0,  q.q0, atol=tol)
    @test isapprox(qi.q1, -q.q1, atol=tol)
    @test isapprox(qi.q2, -q.q2, atol=tol)
    @test isapprox(qi.q3, -q.q3, atol=tol)
end

let
    q  = Quaternion(randn(4))
    q2 = -q

    tol = 1e-12
    @test isapprox(q2.q0, -q.q0, atol=tol)
    @test isapprox(q2.q1, -q.q1, atol=tol)
    @test isapprox(q2.q2, -q.q2, atol=tol)
    @test isapprox(q2.q3, -q.q3, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])
    q2 = Quaternion([1.0 1.0 1.0 1.0])

    q = q2 - q1

    tol = 1e-12
    array_isapprox(q, 0.0, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])
    q2 = Quaternion([1.0 1.0 1.0 1.0])

    q = q2 + q1

    tol = 1e-12
    array_isapprox(q, 1.0, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])

    q = q1 + 1.5

    tol = 1e-12
    array_isapprox(q, 2.0, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])

    q = 1.5 + q1

    tol = 1e-12
    array_isapprox(q, 2.0, atol=tol)
end

let
    q1 = Quaternion([1.0 0.0 0.0 0.0])
    q2 = Quaternion([1.0 0.0 0.0 0.0])

    qp = q2*q1

    tol = 1e-12
    isapprox(qp.q0, 1.0, atol=tol)
    isapprox(qp.q1, 0.0, atol=tol)
    isapprox(qp.q2, 0.0, atol=tol)
    isapprox(qp.q3, 0.0, atol=tol)

    qp = q1*q2

    tol = 1e-12
    isapprox(qp.q0, 1.0, atol=tol)
    isapprox(qp.q1, 0.0, atol=tol)
    isapprox(qp.q2, 0.0, atol=tol)
    isapprox(qp.q3, 0.0, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])

    q = q1*4

    tol = 1e-12
    array_isapprox(q, 2.0, atol=tol)
end

let
    q1 = Quaternion([1.0 1.0 1.0 1.0])

    q = 4*q1

    tol = 1e-12
    array_isapprox(q, 2.0, atol=tol)
end

let
    q1  = Quaternion([1.0 0.0 0.0 0.0])
    q2  = Quaternion([0.0 0.0 0.0 1.0])
    tol = 1e-12

    q = slerp(q1, q2, 0.0)
    array_isapprox(q[:], [1.0 0.0 0.0 0.0], atol=tol)

    q = slerp(q1, q2, 1.0)
    array_isapprox(q[:], [0.0 0.0 0.0 1.0], atol=tol)

    q = slerp(q1, q2, 0.5)
    array_isapprox(q[:], [sqrt(2)/2.0 0.0 0.0 sqrt(2)/2.0], atol=tol)

    # Edge Cases
    q1  = Quaternion([1.0 0.0 0.0 0.0])
    q2  = Quaternion([-1.0 0.0 0.0 0.0])

    q = slerp(q1, q2, 0.0)
    array_isapprox(q[:], [1.0 0.0 0.0 0.0], atol=tol)
end

###########################
# EulerAngle Constructors #
###########################

let
    # Test invalid constructors
    @test_throws ArgumentError EulerAngle(999, 0.1, 0.2, 0.3)
    @test_throws ArgumentError EulerAngle(123, [0.1, 0.2, 0.3, 0.4])
    @test_throws ArgumentError EulerAngle(123, zeros(Float64, 4, 4))
end

let
    # Test valid constructors
    e = EulerAngle(123, 0.1, 0.2, 0.3)

    tol = 1.0e-12
    @test isapprox(e.phi, 0.1, atol=tol)
    @test isapprox(e.theta, 0.2, atol=tol)
    @test isapprox(e.psi, 0.3, atol=tol)

    # Test valid constructors
    e = EulerAngle(123, [0.4, 0.5, 0.6])

    tol = 1.0e-12
    @test isapprox(e.phi, 0.4, atol=tol)
    @test isapprox(e.theta, 0.5, atol=tol)
    @test isapprox(e.psi, 0.6, atol=tol)
end

# Test alternate constructors
let
    e = EulerAngle(123, Quaternion(-1.0, 0.0, 0.0, 0.0))
    tol = 1.0e-12
    array_isapprox((e[:]), 0.0, atol=tol)

    e = EulerAngle(123, EulerAxis(0.0, 0.1, 0.2, 0.3))
    tol = 1.0e-12
    array_isapprox((e[:]), 0.0, atol=tol)
end

let
    # Test error cases
    mat = zeros(Float64, 4, 4)
    @test_throws ArgumentError EulerAngle(123, mat)

    mat = [1.0 0.0 0.0;
           0.0 1.0 0.0;
           0.0 0.0 1.0]
    tol = 1.0e-12

    @test_throws ArgumentError EulerAngle(999, mat)

    # Test angle construction
    e = EulerAngle(121, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, pi, atol=tol)

    e = EulerAngle(123, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(131, mat)
    @test isapprox(e.phi, pi, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(132, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(212, mat)
    @test isapprox(e.phi, pi, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(213, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(231, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(232, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, pi, atol=tol)

    e = EulerAngle(312, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(313, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, pi, atol=tol)

    e = EulerAngle(321, mat)
    @test isapprox(e.phi, 0.0, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)

    e = EulerAngle(323, mat)
    @test isapprox(e.phi, pi, atol=tol)
    @test isapprox(e.theta, 0.0, atol=tol)
    @test isapprox(e.psi, 0.0, atol=tol)
end

# Test Operators
let
    e = EulerAngle(123, 0.1, 0.2, 0.3)

    v = e[:]

    tol = 1e-12
    array_isapprox(e[:], [0.1, 0.2, 0.3], atol=tol)

    # Test Access
    @test isapprox(e[1], 0.1, atol=tol)
    @test isapprox(e[2], 0.2, atol=tol)
    @test isapprox(e[3], 0.3, atol=tol)

    @test length(e[2:3]) == 2

    @test_throws BoundsError e[4]
end

let
    q = EulerAngle(123, 0.5, 0.5, 0.5)

    v = as_vector(q)

    tol = 1.0e-12
    array_isapprox(v, 0.5, atol=tol)
end

let
    mat = as_matrix(EulerAngle(123, 0.0, 0.0, 0.0))

    tol = 1.0e-12
    array_isapprox(mat, one(mat), atol=tol)
end

let
    e1 = EulerAngle(123, randn(3))
    e2 = copy(e1)

    @test pointer_from_objref(e1) != pointer_from_objref(e2)
end

let
    e1 = EulerAngle(123, randn(3))
    e2 = deepcopy(e1)

    @test pointer_from_objref(e1) != pointer_from_objref(e2)
end

##########################
# EulerAxis Constructors #
##########################

# Constructor Exceptions
let
    @test_throws ArgumentError EulerAxis([0.1, 0.2, 0.3])
    @test_throws ArgumentError EulerAxis([0.1, 0.2, 0.3])
    @test_throws ArgumentError EulerAxis(zeros(Float64, 4, 4))
end

# Vector Constructor
let
    e = EulerAxis([0.1, 0.2, 0.3, 0.4])

    tol = 1.0e-12
    array_isapprox(e[:], [0.1, 0.2, 0.3, 0.4], atol=tol)

    e = EulerAxis(0.1, [0.2, 0.3, 0.4])

    tol = 1.0e-12
    array_isapprox(e[:], [0.1, 0.2, 0.3, 0.4], atol=tol)
end

# Euler Angle constructor
let
    e = EulerAxis(EulerAngle(123, 0.0, 0.0, 0.0))

    tol = 1.0e-12
    array_isapprox(e[:], 0.0, atol=tol)
end

# Quaternion constructor
let
    mat = [1.0 0.0 0.0;
           0.0 1.0 0.0;
           0.0 0.0 1.0]
    e = EulerAxis(mat)

    tol = 1.0e-12
    array_isapprox(e[:], 0.0, atol=tol)
end

# Matrix constructor
let
    q = Quaternion(1.0, 0.0, 0.0, 0.0)

    e = EulerAxis(q)

    tol = 1.0e-12
    array_isapprox(e[:], 0.0, atol=tol)
end

# Test Operators
let
    e = EulerAxis(0.4, 0.1, 0.2, 0.3)

    v = e[:]

    tol = 1e-12
    array_isapprox((e[:]), [0.4, 0.1, 0.2, 0.3], atol=tol)

    # Test Access
    @test isapprox(e[1], 0.4, atol=tol)
    @test isapprox(e[2], 0.1, atol=tol)
    @test isapprox(e[3], 0.2, atol=tol)
    @test isapprox(e[4], 0.3, atol=tol)

    @test length(e[2:3]) == 2

    @test_throws BoundsError e[5]
end

let
    e = EulerAxis(0.5, 0.5, 0.5, 0.5)

    v = as_vector(e)

    tol = 1.0e-12
    array_isapprox(v, 0.5, atol=tol)
end

let
    e = EulerAxis(0.0, 0.0, 0.0, 0.0)

    mat = as_matrix(e)

    tol = 1.0e-12
    array_isapprox(mat, one(mat), atol=tol)
end

let
    e1 = EulerAxis(randn(4))
    e2 = copy(e1)

    @test pointer_from_objref(e1) != pointer_from_objref(e2)
end

let
    e1 = EulerAxis(randn(4))
    e2 = deepcopy(e1)

    @test pointer_from_objref(e1) != pointer_from_objref(e2)
end


text:
@testset "Test CLI" begin
    for (args, expected_output) in [
        (["help"], m.HELP_INFO * "\n"),
        (["ask"], m.HELP_INFO * "\n"),
        (["run"], "Usage: monkey run <file> [--vm | --jl]\n"),
        (["run", "fixtures/hello_world.m"], "ERROR: Only .mo files are supported!\n"),
        (["run", "fixtures/not_exist.mo"], "ERROR: File not found!\n"),
        (["run", "fixtures/hello_world.mo"], "Hello, world!\n"),
        (["run", "fixtures/hello_world.mo", "--vm"], "Hello, world!\n"),
        (["run", "fixtures/hello_world.mo", "--jl"], "Hello, world!\n"),
        (
            ["run", "fixtures/hello_world.mo", "--invalid"],
            "Usage: monkey run <file> [--vm | --jl]\n",
        ),
        (["repl"], m.REPL_PRELUDE * "\n>> " * m.REPL_FAREWELL * "\n"),
        (["repl", "--vm"], m.REPL_PRELUDE * "\n>> " * m.REPL_FAREWELL * "\n"),
    ]
        input = IOBuffer()
        output = IOBuffer(UInt8[], read = true, write = true)

        for arg in args
            push!(ARGS, arg)
        end

        m.julia_main(; input = input, output = output)
        @test String(output.data) == expected_output

        for _ = 1:length(args)
            pop!(ARGS)
        end
    end
end