<|start_filename|>src/main.jl<|end_filename|>
####
#### CLI settings
####

const ALLOWED_FILE_SUFFIXES = (".mat", ".nii", ".nii.gz", ".par", ".xml", ".rec")
const ALLOWED_FILE_SUFFIXES_STRING = join(ALLOWED_FILE_SUFFIXES, ", ", ", and ")

const CLI_SETTINGS = ArgParseSettings(
    prog = "",
    fromfile_prefix_chars = "@",
    error_on_conflict = false,
    exit_after_help = false,
)

@add_arg_table! CLI_SETTINGS begin
    "input"
        nargs = '+' # At least one input is required
        arg_type = String
        required = true
        help = "one or more input filenames. Valid file types are limited to: $ALLOWED_FILE_SUFFIXES_STRING"
    "--mask", "-m"
        nargs = '+' # If --mask is passed, at least one input is required
        arg_type = String
        help = "one or more mask filenames. Masks are loaded and subsequently applied to the corresponding input files via elementwise multiplication. The number of mask files must equal the number of input files. Valid file types are the same as for input files, and are limited to: $ALLOWED_FILE_SUFFIXES_STRING"
    "--output", "-o"
        arg_type = String
        help = "output directory. If not specified, output file(s) will be stored in the same location as the corresponding input file(s). Outputs are stored with the same basename as the input files with additional suffixes; see --T2map and --T2part"
    "--T2map"
        action = :store_true
        help = "call T2mapSEcorr to compute T2 distributions from 4D multi spin-echo input images. T2 distributions and T2 maps produced by T2mapSEcorr are saved as MAT files with extensions .t2dist.mat and .t2maps.mat"
    "--T2part"
        action = :store_true
        help = "call T2partSEcorr to analyze 4D T2 distributions to produce parameter maps. If --T2map is also passed, input 4D arrays are interpreted as multi spin-echo images and T2 distributions are first computed by T2mapSEcorr. If only --T2part is passed, input 4D arrays are interpreted as T2 distributions and only T2partSEcorr is called. Output T2 parts are saved as a MAT file with extension .t2parts.mat"
    "--quiet", "-q"
        action = :store_true
        help = "suppress printing to the terminal. Note: all terminal outputs, including errors and warnings, are still printed to the log file"
    "--dry"
        action = :store_true
        help = "execute dry run of processing without saving any results"
    "--legacy"
        action = :store_true
        help = "use legacy settings and algorithms from the original MATLAB pipeline. This ensures that the same T2-distributions and T2-parts will be produced as those from MATLAB. Note that execution time will be much slower, and less robust algorithms will be used."
end

add_arg_group!(CLI_SETTINGS,
    "T2mapSEcorr/T2partSEcorr arguments",
    :t2_map_part,
)

@add_arg_table! CLI_SETTINGS begin
    "--MatrixSize"
        nargs = 3
        arg_type = Int
        help = "first three dimensions of input 4D image. Inferred automatically"
        group = :t2_map_part
    "--nTE"
        arg_type = Int
        help = "number of echoes in input signal. Inferred automatically when --T2map is passed"
        group = :t2_map_part
    "--TE"
        arg_type = Float64
        help = "inter-echo spacing (Units: seconds). Required when --T2map is passed"
        group = :t2_map_part
    "--nT2"
        arg_type = Int
        help = "number of T2 components used in the multi-exponential analysis. Required when --T2map is passed. Inferred from fourth dimension of input image if only --T2part is passed"
        group = :t2_map_part
    "--T2Range"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values (Units: seconds). T2 components are logarithmically spaced between these bounds. Required parameter."
        group = :t2_map_part
    "--SPWin"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values of the short peak window (Units: seconds). Required parameter when --T2part is passed"
        group = :t2_map_part
    "--MPWin"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values of the middle peak window (Units: seconds). Required parameter when --T2part is passed"
        group = :t2_map_part
    "--T1"
        arg_type = Float64
        help = "assumed value of longitudinal T1 relaxation (Units: seconds)."
        group = :t2_map_part
    "--Reg"
        arg_type = String
        help = "routine used for choosing regularization parameter. One of \"none\", \"chi2\", \"gcv\", or \"lcurve\", representing no regularization, --Chi2Factor based Tikhonov regularization, generalized cross-validation based regularization, and L-curve based regularization, respectively."
        group = :t2_map_part
    "--Chi2Factor"
        arg_type = Float64
        help = "if --Reg=\"chi2\", the T2 distribution is regularized such that the chi^2 goodness of fit is increased by a multiplicative factor --Chi2Factor relative to the unregularized solution"
        group = :t2_map_part
    "--Sigmoid"
        arg_type = Float64
        help = "replace the hard upper limit cutoff time of the short peak window, --SPWin[2], with a smoothed sigmoidal cutoff function, f, scaled and shifted such that f(--SPWin[2] +/- --Sigmoid) = 0.5 -/+ 0.4. --Sigmoid is the time scale of the smoothing (Units: seconds)"
        group = :t2_map_part
    "--Threshold"
        arg_type = Float64
        help = "first echo intensity cutoff for empty voxels. Processing is skipped for voxels with intensity <= --Threshold"
        group = :t2_map_part
end

add_arg_group!(CLI_SETTINGS,
    "B1 correction and stimulated echo correction",
    :B1_SE_corr,
)

@add_arg_table! CLI_SETTINGS begin
    "--nRefAngles"
        arg_type = Int
        help = "in estimating the local refocusing flip angle to correct for B1 inhomogeneities, up to --nRefAngles angles in the range [--MinRefAngle, 180] are explicitly checked. The optimal angle is then determined through interpolation from these --nRefAngles observations."
        group = :B1_SE_corr
    "--nRefAnglesMin"
        arg_type = Int
        help = "initial number of angles to check during flip angle estimation before refinement near likely optima. Setting --nRefAnglesMin equal to --nRefAngles forces all angles to be checked."
        group = :B1_SE_corr
    "--MinRefAngle"
        arg_type = Float64
        help = "minimum refocusing angle for flip angle estimation (Units: degrees)."
        group = :B1_SE_corr
    "--SetFlipAngle"
        arg_type = Float64
        help = "to skip B1 inhomogeneity correction, use --SetFlipAngle to assume a fixed refocusing flip angle for all voxels (Units: degrees)."
        group = :B1_SE_corr
    "--SetRefConAngle"
        arg_type = Float64
        help = "refocusing pulse control angle for stimulated echo correction. Unlike B1 inhomogeneity correction, stimulated echo correction must be performed manually. By default, --SetRefConAngle is set to 180 degrees, equivalent to no stimulated echo correction (Units: degrees)."
        group = :B1_SE_corr
    "--RefConAngle"
        arg_type = Float64
        help = "Deprecated flag; see --SetRefConAngle."
        group = :B1_SE_corr
end

add_arg_group!(CLI_SETTINGS,
    "Additional save options",
    :save_opts,
)

@add_arg_table! CLI_SETTINGS begin
    "--SaveDecayCurve"
        action = :store_true
        help = "include a 4D array of the time domain decay curves resulting from the NNLS fits in the output maps dictionary"
        group = :save_opts
    "--SaveNNLSBasis"
        action = :store_true
        help = "include a 5D (or 2D if --SetFlipAngle is used) array of NNLS basis matrices in the output maps dictionary. Note: this 5D array is extremely large for typical image sizes; in most cases, this flag should only be set when debugging small datasets"
        group = :save_opts
    "--SaveRegParam"
        action = :store_true
        help = "include 3D arrays of the regularization parameters and resulting chi^2-factors in the output maps dictionary"
        group = :save_opts
    "--SaveResidualNorm"
        action = :store_true
        help = "include a 3D array of the l2-norms of the residuals from the NNLS fits in the output maps dictionary"
        group = :save_opts
end

add_arg_group!(CLI_SETTINGS,
    "BET arguments",
    :bet_args,
)

@add_arg_table! CLI_SETTINGS begin
    "--bet"
        action = :store_true
        help = "use the BET brain extraction tool from the FSL library of analyis tools to automatically create a binary brain mask. Only voxels within the binary mask will be analyzed. Note that if a mask is passed explicitly with the --mask flag, this mask will be used and the --bet flag will be ignored"
        group = :bet_args
    "--betargs"
        arg_type = String
        default = "-m -n -f 0.25 -R"
        help = "BET command line interface arguments. Must be passed as a single string with arguments separated by spaces, e.g. '-m -n'. The flag '-m' indicates that a binary mask should be computed, and therefore will be added to the list of arguments if not provided"
        group = :bet_args
    "--betpath"
        arg_type = String
        default = "bet"
        help = "path to BET executable."
        group = :bet_args
end

"""
    main(command_line_args = ARGS)

Entry point function for command line interface, parsing the command line arguments `ARGS` and subsequently calling one or both of `T2mapSEcorr` and `T2partSEcorr` with the parsed settings.
See the [Arguments](@ref) section for available options.

See also:
* [`T2mapSEcorr`](@ref)
* [`T2partSEcorr`](@ref)
"""
function main(command_line_args::Vector{String} = ARGS)

    # Parse command line arguments
    opts = parse_args(command_line_args, CLI_SETTINGS; as_symbols = true)
    opts === nothing && return # Help message was triggered. Return nothing instead of exit(0)
    opts = handle_cli_deprecations!(opts)

    # Get input file list and output folder list
    for file_info in get_file_infos(opts)
        # Make output path
        if !opts[:dry]
            mkpath(file_info[:outputfolder])
        end

        # Save settings files
        if !opts[:dry]
            map(filter(s -> startswith(s, "@"), command_line_args)) do settingsfile
                src = settingsfile[2:end] # drop "@" character
                dst = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * "." * basename(src))
                cp(src, dst; force = true)
            end
        end

        # Main processing
        tee_capture(
            logfile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".log"),
            suppress_terminal = opts[:quiet],
            suppress_logfile = opts[:dry],
        ) do io
            try
                main_(io, file_info, opts)
            catch e
                @warn "Error during processing of file: $(file_info[:inputfile])"
                @warn sprint(showerror, e, catch_backtrace())
            end
        end
    end

    return nothing
end

function main_(io::IO, file_info::Dict{Symbol,Any}, opts::Dict{Symbol,Any})

    # Starting message/starting time
    t_start = tic()
    printheader(io, "Starting DECAES with $(Threads.nthreads()) threads")

    # Load image(s)
    image = @showtime(io,
        "Loading input file: $(file_info[:inputfile])",
        load_image(file_info[:inputfile]),
    )

    # Apply mask
    if file_info[:maskfile] !== nothing
        @showtime(io,
            "Applying mask from file: '$(file_info[:maskfile])'",
            try_apply_maskfile!(image, file_info[:maskfile]),
        )
    elseif opts[:bet]
        @showtime(io,
            "Making and applying BET mask with args: '$(opts[:betargs])'",
            try_apply_bet!(image, opts[:betpath], opts[:betargs]),
        )
    end

    # Compute T2 distribution from input 4D multi-echo image
    if opts[:T2map]
        maps, dist = @showtime(io,
            "Running T2mapSEcorr on file: $(file_info[:inputfile])",
            T2mapSEcorr(io, image, t2map_options(image, opts)),
        )

        # Save T2-distribution to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2dist.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 distribution to file: $savefile",
                MAT.matwrite(savefile, Dict{String,Any}("dist" => dist)),
            )
        end

        # Save T2-maps to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2maps.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 parameter maps to file: $savefile",
                MAT.matwrite(savefile, maps),
            )
        end
    else
        # Input image is the T2 distribution
        dist = image
    end

    # Analyze T2 distribution to produce parameter maps
    if opts[:T2part]
        parts = @showtime(io,
            "Running T2partSEcorr",
            T2partSEcorr(io, dist, t2part_options(dist, opts)),
        )

        # Save T2-parts to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2parts.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 parts maps to file: $savefile",
                MAT.matwrite(savefile, parts),
            )
        end
    end

    # Done message
    printheader(io, "Finished ($(round(toc(t_start); digits = 2)) seconds)")

    return nothing
end

"""
Entrypoint function for compiling DECAES into an executable [app](https://julialang.github.io/PackageCompiler.jl/dev/apps/).
"""
function julia_main()::Cint
    try
        main(ARGS)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

####
#### Helper functions
####

function handle_cli_deprecations!(opts)
    handle_renamed_cli_flag!(opts, :RefConAngle => :SetRefConAngle)
end

function handle_renamed_cli_flag!(opts, oldnew::Pair{Symbol, Symbol})
    oldflag, newflag = oldnew
    if opts[oldflag] !== nothing
        if opts[newflag] !== nothing
            error("The flag --$newflag and the deprecated flag --$oldflag were both passed; use --$newflag only.")
        else
            @warn "The flag --$oldflag is deprecated and will be removed in future releases; use --$newflag instead."
        end
        opts[newflag] = opts[oldflag]
        delete!(opts, oldflag)
    end
    return opts
end

function t2map_options(image::Array, opts::Dict{Symbol,Any})
    fields = fieldsof(T2mapOptions, Set)
    kwargs = Dict{Symbol,Any}()
    for (k,v) in opts
        (v === nothing) && continue # filter unset cli args
        (v isa AbstractVector && isempty(v)) && continue # filter unset cli args (empty vectors are unset cli varargs)
        (k ∉ fields) && continue # filter T2mapOptions fields
        kwargs[k] = v isa AbstractVector ? tuple(v...) : v
    end
    T2mapOptions(image; kwargs...)
end

function t2part_options(dist::Array, opts::Dict{Symbol,Any})
    fields = fieldsof(T2partOptions, Set)
    kwargs = Dict{Symbol,Any}()
    for (k,v) in opts
        (v === nothing) && continue # filter unset cli args
        (v isa AbstractVector && isempty(v)) && continue # filter unset cli args (empty vectors are unset cli varargs)
        (k === :nT2 && !opts[:T2map]) && continue # nT2 must be explicitly passed, unless not performing T2-mapping, in which case it is inferred from `dist`
        (k ∉ fields) && continue # filter T2mapOptions fields
        kwargs[k] = v isa AbstractVector ? tuple(v...) : v
    end
    T2partOptions(dist; kwargs...)
end

function get_file_infos(opts::Dict{Symbol,Any})
    @unpack input, output, mask = opts

    # Read in input files
    inputfiles = String[path for path in input if is_allowed_suffix(path)]

    if isempty(inputfiles)
        msg = if !isempty(input) && isfile(input[1])
            "No valid file types were found for processing, but a file name was passed.\n" *
            "Perhaps you meant to prepend an '@' character to a settings file, e.g. @$(input[1])?\n" *
            "If not, note that only $ALLOWED_FILE_SUFFIXES_STRING file types are supported"
        else
            "No valid files were found for processing. Note that currently only $ALLOWED_FILE_SUFFIXES_STRING file types are supported"
        end
        error(msg)
    end

    # Get output folders
    outputfolders = if output === nothing
        dirname.(inputfiles)
    else
        [output for _ in 1:length(inputfiles)]
    end

    # Get mask files
    maskfiles = if isempty(mask)
        fill(nothing, length(inputfiles))
    elseif length(mask) == length(inputfiles)
        String.(mask)
    else
        error("Number of mask files passed does not equal the number of input image files passed")
    end

    # Create file_info dictionaries
    file_info = Dict{Symbol,Any}[]
    for (inputfile, outputfolder, maskfile) in zip(inputfiles, outputfolders, maskfiles)
        d = eltype(file_info)(
            :inputfile => inputfile,
            :outputfolder => outputfolder,
            :maskfile => maskfile,
            :choppedinputfile => chop_allowed_suffix(basename(inputfile)),
        )
        push!(file_info, d)
    end

    return file_info
end

function load_image(filename, ::Val{N}) where {N}
    image = if maybe_get_suffix(filename) == ".mat"
        # Load first `N`-dimensional array which is found, or throw an error if none are found
        data = MAT.matread(filename)
        array_keys = findall(x -> x isa AbstractArray{T,N} where {T}, data)
        if isempty(array_keys)
            error("No 4D array was found in the input file: $filename")
        end
        if length(array_keys) > 1
            array_keys = sort(array_keys)
            @warn "Multiple possible images found in file: $(filename)\nChoosing field $(repr(array_keys[1])) out of the following options: $(join(repr.(array_keys), ", "))"
        end
        data[array_keys[1]]

    elseif maybe_get_suffix(filename) ∈ (".nii", ".nii.gz")
        # Check slope field; if scl_slope == 0, data is not scaled and raw data should be returned
        #   See e.g. https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/scl_slopeinter.html
        # Loaded data is coerced to a `N`-dimensional array
        data = NIfTI.niread(filename)
        if data.header.scl_slope == 0
            data.raw[ntuple(_ -> Colon(), N)...] # Return raw data
        else
            data[ntuple(_ -> Colon(), N)...] # Return scaled data (getindex from the NIfTI package handles scaling)
        end

    elseif maybe_get_suffix(filename) ∈ (".par", ".xml", ".rec")
        # Load PAR/REC or XML/REC file, coercing resulting data into a `N`-dimensional array
        _, data = ParXRec.parxrec(filename)
        data[ntuple(_ -> Colon(), N)...]

    else
        error("Currently, only $ALLOWED_FILE_SUFFIXES_STRING files are supported")
    end

    # Currently, the pipeline is ~twice as fast on Float64 arrays than Float32 arrays (unclear why).
    # However, the MATLAB toolbox converts images to double as well, so here we simply do the same
    image = convert(Array{Float64,N}, image)

    return image
end
load_image(filename; ndims::Int = 4) = load_image(filename, Val(ndims))

function try_apply_maskfile!(image, maskfile)
    try
        image .*= load_image(maskfile, Val(3))
    catch e
        @warn "Error while loading mask file: $maskfile"
        @warn sprint(showerror, e, catch_backtrace())
    end
    return image
end

function try_apply_bet!(image, betpath, betargs)
    try
        image .*= make_bet_mask(image, betpath, betargs)
    catch e
        @warn "Error while making mask using BET"
        @warn sprint(showerror, e, catch_backtrace())
    end
    return image
end

function make_bet_mask(image::Array{T,3}, betpath, betargs) where {T}
    # Split betargs, and ensure that "-m" (make binary mask) is among args
    args = convert(Vector{String}, filter!(!isempty, split(betargs, " ")))
    if "-m" ∉ args
        push!(args, "-m")
    end

    # Create mask using BET and return mask
    mask = mktempdir() do temppath
        tempbase = basename(tempname())
        nifti_imagefile = joinpath(temppath, tempbase * ".nii")
        nifti_maskfile = joinpath(temppath, tempbase * ".bet")
        NIfTI.niwrite(nifti_imagefile, NIfTI.NIVolume(image)) # create nifti file for bet
        run(Cmd([
            betpath;
            nifti_imagefile;
            nifti_maskfile;
            args
        ]))
        # BET appends "_mask" and ".nii.gz" to output file name.
        # Find this file, ensure it is unique, then load and return it
        bet_maskfiles = filter!(file -> startswith(file, tempbase * ".bet_mask"), readdir(temppath))
        @assert length(bet_maskfiles) == 1 # ensure unique; this should never be false using a temp filename
        load_image(joinpath(temppath, bet_maskfiles[1]), Val(3))
    end

    return mask
end
make_bet_mask(image::Array{T,4}, args...; kwargs...) where {T} = make_bet_mask(image[:,:,:,1], args...; kwargs...) # use first echo

maybe_get_first(f, xs) = findfirst(f, xs) |> I -> I === nothing ? nothing : xs[I]
maybe_get_suffix(filename) = maybe_get_first(ext -> endswith(lowercase(filename), ext), ALLOWED_FILE_SUFFIXES) # case-insensitive
is_allowed_suffix(filename) = maybe_get_suffix(filename) !== nothing

function chop_allowed_suffix(filename::AbstractString)
    ext = maybe_get_suffix(filename)
    if ext !== nothing
        return filename[1:end-length(ext)]
    else
        error("Currently only $ALLOWED_FILE_SUFFIXES_STRING file types are supported")
    end
end


<|start_filename|>src/lsqnonneg_reg.jl<|end_filename|>
####
#### Lazy wrappers for LHS matrix and RHS vector for augmented Tikhonov-regularized NNLS problems
####

struct PaddedVector{T, Vd <: AbstractVector{T}} <: AbstractVector{T}
    d::Vd
    pad::Int
end
Base.size(x::PaddedVector) = (length(x.d) + x.pad,)

function Base.copyto!(y::AbstractVector{T}, x::PaddedVector{T}) where {T}
    @assert size(x) == size(y)
    @unpack d, pad = x
    m = length(x.d)
    @inbounds @simd ivdep for i in 1:m
        y[i] = d[i]
    end
    @inbounds @simd ivdep for i in m+1:m+pad
        y[i] = zero(T)
    end
    return y
end

struct TikhonovPaddedMatrix{T, MC <: AbstractMatrix{T}} <: AbstractMatrix{T}
    C::MC
    μ::Base.RefValue{T}
end
TikhonovPaddedMatrix(C::AbstractMatrix, μ) = TikhonovPaddedMatrix(C, Ref(μ))
Base.size(A::TikhonovPaddedMatrix) = ((m,n) = size(A.C); return (m+n,n))

function Base.copyto!(B::AbstractMatrix{T}, A::TikhonovPaddedMatrix{T}) where {T}
    @assert size(A) == size(B)
    @unpack C, μ = A
    m, n = size(C)
    @inbounds @simd ivdep for j in 1:n
        @simd ivdep for i in 1:m
            B[i,j] = C[i,j]
        end
        @simd ivdep for i in m+1:m+n
            B[i,j] = ifelse(i == m+j, μ[], zero(T))
        end
    end
    return B
end

cache!(cache::AbstractArray, x, f) = cache !== nothing && push!(cache, (; x, f))
cache!(cache::Dict, x, f) = cache !== nothing && (cache[x] = f)

lin_interp(x, x₁, x₂, y₁, y₂) = y₁ + (y₂ - y₁) * (x - x₁)  / (x₂ - x₁)
exp_interp(x, x₁, x₂, y₁, y₂) = y₁ + log1p(expm1(y₂ - y₁) * (x - x₁)  / (x₂ - x₁))

####
#### Tikhonov regularized NNLS problem
####

struct NNLSTikhonovRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W
end
function NNLSTikhonovRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(TikhonovPaddedMatrix(C, T(NaN)), PaddedVector(d, n))
    NNLSTikhonovRegProblem(C, d, m, n, nnls_work)
end

mu(work::NNLSTikhonovRegProblem) = work.nnls_work.C.μ[]
mu!(work::NNLSTikhonovRegProblem, μ) = work.nnls_work.C.μ[] = μ

function solve!(work::NNLSTikhonovRegProblem, μ)
    # Set regularization parameter and solve NNLS problem
    mu!(work, μ)
    solve!(work.nnls_work)
    return solution(work)
end

solution(work::NNLSTikhonovRegProblem) = solution(work.nnls_work)

resnorm_sq(work::NNLSTikhonovRegProblem) = chi2(work)
resnorm(work::NNLSTikhonovRegProblem) = sqrt(resnorm_sq(work))

seminorm_sq(work::NNLSTikhonovRegProblem) = solution(work) ⋅ solution(work)
seminorm(work::NNLSTikhonovRegProblem) = sqrt(seminorm_sq(work))

loss(work::NNLSTikhonovRegProblem) = chi2(work.nnls_work)

reg(work::NNLSTikhonovRegProblem) = mu(work)^2 * seminorm_sq(work)

chi2(work::NNLSTikhonovRegProblem) = max(loss(work) - reg(work), 0)

# Extract columns spanning solution space of regularized NNLS problem
function subproblem(work::NNLSTikhonovRegProblem{T}) where {T}
    @unpack C, d, m, n, nnls_work = work
    μ = mu(work)
    x = solution(work)
    Jnz = findall(!≈(0), x)
    nnz = length(Jnz)
    return TikhonovPaddedMatrix(C[:,Jnz], μ), PaddedVector(d, nnz), x[Jnz]
end

function ∇reg(work::NNLSTikhonovRegProblem)
    @unpack C, d = work
    μ = mu(work)
    x = solution(work)
    C′, d′, x′ = subproblem(work)
    b′ = PaddedVector(C′.C' \ x′, d′.pad)
    ∇nnls_work = NNLSProblem(C′, b′)
    solve!(∇nnls_work)
    ∇x′ = -solution(∇nnls_work)
    ∇μ² = 2 * ((C′.C * x′ - d′.d) ⋅ (C′.C * ∇x′)) #+ (x′ ⋅ x′) + 2 * μ^2 * (x′ ⋅ ∇x′)
    ∇μ = 2 * μ * ∇μ²
    return ∇μ
end

function chi2factor_relerr!(work::NNLSTikhonovRegProblem, logμ, ∇logμ = nothing; χ²target)
    μ = exp(logμ)
    solve!(work, μ)
    χ² = chi2(work)
    relerr = log(χ² / χ²target) # better behaved than χ² / χ²target - 1 for large χ²
    if ∇logμ !== nothing && length(∇logμ) > 0
        ∂χ²_∂logμ = ∇reg(work)
        ∂relerr_∂logμ = μ * ∂χ²_∂logμ / χ²
        @inbounds ∇logμ[1] = ∂relerr_∂logμ
    end
    return relerr
end
chi2factor_relerr⁻¹(relerr; χ²target) = χ²target * exp(relerr)

function chi2factor_loss!(work::NNLSTikhonovRegProblem, logμ, ∇logμ = nothing; χ²target)
    relerr = chi2factor_relerr!(work, logμ, ∇logμ; χ²target)
    loss = abs(relerr)
    if ∇logμ !== nothing && length(∇logμ) > 0
        @inbounds ∂relerr_∂logμ = ∇logμ[1]
        ∂loss_∂logμ = sign(relerr) * ∂relerr_∂logμ
        @inbounds ∇logμ[1] = ∂loss_∂logμ
    end
    return loss
end

"""
Helper struct which wraps `N` caches of `NNLSTikhonovRegProblem` workspaces.
Useful for optimization problems where the last function call may not be
the optimium, but perhaps it was one or two calls previous and is still in the
`NNLSTikhonovRegProblemCache` and a recomputation can be avoided.
"""
struct NNLSTikhonovRegProblemCache{N,W}
    cache::NTuple{N,W}
    idx::Base.RefValue{Int}
end
function NNLSTikhonovRegProblemCache(C::AbstractMatrix{T}, d::AbstractVector{T}, ::Val{N} = Val(5)) where {T,N}
    cache = ntuple(_ -> NNLSTikhonovRegProblem(C, d), N)
    idx = Ref(1)
    return NNLSTikhonovRegProblemCache(cache, idx)
end
increment_cache_index!(work::NNLSTikhonovRegProblemCache{N}) where {N} = (work.idx[] = mod1(work.idx[] + 1, N))
set_cache_index!(work::NNLSTikhonovRegProblemCache{N}, i) where {N} = (work.idx[] = mod1(i, N))
reset_cache!(work::NNLSTikhonovRegProblemCache{N}) where {N} = foreach(w -> mu!(w, NaN), work.cache)
get_cache(work::NNLSTikhonovRegProblemCache) = work.cache[work.idx[]]

function solve!(work::NNLSTikhonovRegProblemCache, μ)
    i = findfirst(w -> μ == mu(w), work.cache)
    if i === nothing
        increment_cache_index!(work)
        solve!(get_cache(work), μ)
    else
        set_cache_index!(work, i)
    end
    return solution(get_cache(work))
end


####
#### Chi2 method for choosing Tikhonov regularization parameter
####

struct NNLSChi2RegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSChi2RegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSChi2RegProblem(C, d, m, n, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSChi2RegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_chi2(C::AbstractMatrix, d::AbstractVector, Chi2Factor::Real)

Returns the regularized NNLS solution, X, that incurrs an increase in ``\\chi^2`` approximately by a factor of `Chi2Factor`.
The regularized NNLS problem solved internally is:

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||x||_2^2
```

where ``\\mu`` is determined by approximating a solution to the nonlinear equation

```math
\\frac{\\chi^2(\\mu)}{\\chi^2_{min}} = \\mathrm{Chi2Factor}
\\quad
\\text{where}
\\quad
\\chi^2_{min} = \\chi^2(\\mu = 0)
```

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data
- `Chi2Factor::Real`: Desired ``\\chi^2`` increase due to regularization

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Actual increase ``\\chi^2(\\mu)/\\chi^2_{min}``, which will be approximately equal to the input `Chi2Factor`
"""
function lsqnonneg_chi2(C, d, Chi2Factor; kwargs...)
    work = lsqnonneg_chi2_work(C, d)
    lsqnonneg_chi2!(work, Chi2Factor; kwargs...)
end
lsqnonneg_chi2_work(C, d) = NNLSChi2RegProblem(C, d)

function lsqnonneg_chi2!(work::NNLSChi2RegProblem{T}, Chi2Factor::T; bisection = true, legacy = false) where {T}
    # Non-regularized solution
    solve!(work.nnls_work)
    chi2_min = chi2(work.nnls_work)

    # Prepare to solve
    χ²target = Chi2Factor * chi2_min
    reset_cache!(work.nnls_work_smooth_cache)

    if legacy
        # Use the legacy algorithm: double μ starting from an initial guess, then interpolate the root using a cubic spline fit 
        mu_final, chi2_final = chi2factor_search_from_minimum(chi2_min, Chi2Factor; legacy = legacy) do μ
            μ == 0 && return chi2_min
            solve!(work.nnls_work_smooth_cache, μ)
            return chi2(get_cache(work.nnls_work_smooth_cache))
        end
        if mu_final == 0
            x_final = solution(work.nnls_work)
        else
            x_final = solve!(work.nnls_work_smooth_cache, mu_final)
        end

    elseif bisection
        # Find bracketing interval containing root, then perform bisection search
        f = function (logμ)
            increment_cache_index!(work.nnls_work_smooth_cache)
            return chi2factor_relerr!(get_cache(work.nnls_work_smooth_cache), logμ; χ²target)
        end
        cache = NamedTuple{(:x, :f), NTuple{2,T}}[]
        a, b, fa, fb = bracketing_interval(f, T(-4.0), T(1.0), T(1.5); maxiters = 6, cache)
        bisect(f, a, b, fa, fb; xtol = T(0.05), ftol = (Chi2Factor-1)/100, cache)

        # Spline rootfinding on evaluated points to improve accuracy
        sort!(cache; by = d -> d.x)
        logmu_root = spline_root([d.x for d in cache], [d.f for d in cache]; deg_spline = 1)
        logmu_root !== nothing && !any(d -> d.x ≈ logmu_root, cache) && cache!(cache, logmu_root, f(logmu_root))

        # Return regularization which minimizes relerr
        _, i = findmin([abs(d.f) for d in cache])
        logmu_final, relerr_final = cache[i]
        mu_final, chi2_final = exp(logmu_final), chi2factor_relerr⁻¹(relerr_final; χ²target)
        x_final = solve!(work.nnls_work_smooth_cache, mu_final)

    else
        # Instead of rootfinding, reformulate as a minimization problem. Solve using NLopt.
        alg = :LN_COBYLA # local, gradient-free, linear approximation of objective
        # alg = :LN_BOBYQA # local, gradient-free, quadratic approximation of objective
        # alg = :GN_AGS # global, gradient-free, hilbert curve based dimension reduction
        # alg = :LN_NELDERMEAD # local, gradient-free, simplex method
        # alg = :LN_SBPLX # local, gradient-free, subspace searching simplex method
        # alg = :LD_CCSAQ # local, first-order (rough ranking: [:LD_MMA, :LD_SLSQP, :LD_LBFGS, :LD_CCSAQ, :LD_AUGLAG])

        opt = NLopt.Opt(alg, 1)
        opt.lower_bounds  = -8.0
        opt.upper_bounds  = 2.0
        opt.xtol_rel      = 0.05
        opt.min_objective = function (logμ, ∇logμ)
            @inbounds _logμ = logμ[1]
            increment_cache_index!(work.nnls_work_smooth_cache)
            loss = chi2factor_loss!(get_cache(work.nnls_work_smooth_cache), _logμ, ∇logμ; χ²target)
            return Float64(loss)
        end
        minf, minx, ret   = NLopt.optimize(opt, [-4.0])

        mu_final = exp(T(minx[1]))
        x_final = solve!(work.nnls_work_smooth_cache, mu_final)
        chi2_final = chi2(get_cache(work.nnls_work_smooth_cache))
    end

    return (x = x_final, mu = mu_final, chi2factor = chi2_final/chi2_min)
end

function chi2factor_search_from_minimum(f, χ²min::T, χ²fact::T, μmin::T = T(1e-3), μfact = T(2.0); legacy = false) where {T}
    # Minimize energy of spectrum; loop to find largest μ that keeps chi-squared in desired range
    μ_cache = T[zero(T)]
    χ²_cache = T[χ²min]
    μnew = μmin
    while true
        # Cache function value at μ = μnew
        χ²new = f(μnew)
        push!(μ_cache, μnew)
        push!(χ²_cache, χ²new)

        # Break when χ²fact reached, else increase regularization
        (χ²new >= χ²fact * χ²min) && break
        μnew *= μfact
    end

    # Solve χ²(μ) = χ²fact * χ²min using a spline fitting root finding method
    if legacy
        # Legacy algorithm fits spline to all (μ, χ²) values observed, including for μ=0.
        # This poses several problems:
        #   1) while unlikely, it is possible for the spline to return a negative regularization parameter
        #   2) the μ values are exponentially spaced, leading to poorly conditioned splines
        μ = spline_root_legacy_slow(μ_cache, χ²_cache, χ²fact * χ²min)
        μ = μ === nothing ? μmin : μ
    else
        if length(μ_cache) == 2
            # Solution is contained in [0,μmin]; `spline_root` with two points performs root finding via simple linear interpolation
            μ = spline_root(μ_cache, χ²_cache, χ²fact * χ²min; deg_spline = 1)
            μ = μ === nothing ? μmin : μ
        else
            # Perform spline fit on log-log scale on data with μ > 0. This solves the above problems with the legacy algorithm:
            #   1) Root is found in terms of logμ, guaranteeing μ > 0
            #   2) logμ is linearly spaced, leading to well-conditioned splines
            logμ = spline_root(log.(μ_cache[2:end]), log.(χ²_cache[2:end]), log(χ²fact * χ²min); deg_spline = 1)
            μ = logμ === nothing ? μmin : exp(logμ)
        end
    end

    # Compute the final regularized solution
    χ² = f(μ)

    return μ, χ²
end

function chi2factor_search_from_guess(f, χ²min::T, χ²fact::T, μ₀::T = T(1e-2), μfact = T(1.5), μmin::T = T(1e-4)) where {T}
    # Find interval containing χ²target = χ²fact * χ²min
    χ²target = χ²fact * χ²min
    μnew = μ₀
    χ²new = f(μnew)
    logμ_cache = T[log(μnew)]
    logχ²_cache = T[log(χ²new)]

    if χ²new > χ²target
        while μnew > μmin && !(μnew ≈ μmin)
            μnew = max(μnew/μfact, μmin)
            χ²new = f(μnew)
            pushfirst!(logμ_cache, log(μnew))
            pushfirst!(logχ²_cache, log(χ²new))
            (χ²new < χ²target && length(logμ_cache) >= 3) && break
        end
    else
        while true
            μnew *= μfact
            χ²new = f(μnew)
            push!(logμ_cache, log(μnew))
            push!(logχ²_cache, log(χ²new))
            (χ²new > χ²target && length(logμ_cache) >= 3) && break
        end
    end

    if logμ_cache[1] ≈ log(μmin) && logχ²_cache[1] > log(χ²target)
        # μ decreased to μmin but χ²target was not reached; linearly interpolate between μ=0 and μ=μmin points
        μ = spline_root([zero(T), μmin], [log(χ²min), logχ²_cache[1]], log(χ²target))
        μ = μ === nothing ? μmin : μ
        χ² = f(μ)
    else
        # Find optimal μ and evaluate at the interpolated solution
        logμ = spline_root(logμ_cache, logχ²_cache, log(χ²target))
        μ = logμ === nothing ? μmin : exp(logμ)
        χ² = f(μ)
    end

    return μ, χ²
end

####
#### Rootfinding methods. `secant_method` and `bisection_method` are modified codes from Roots.jl:
####
####    https://github.com/JuliaMath/Roots.jl/blob/8a5ff76e8e8305d4ad5719fe1dd665d8a7bd7ec3/src/simple.jl
####
#### The MIT License (MIT) Copyright (c) 2013 <NAME>
#### Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#### The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#### THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

function secant_method(f, xs; xtol = zero(float(real(first(xs)))), xrtol = 8*eps(one(float(real(first(xs))))), maxiters = 1000, cache = nothing)
    if length(xs) == 1 # secant needs a, b; only a given
        a  = float(xs[1])
        h  = eps(one(real(a)))^(1/3)
        da = h * oneunit(a) + abs(a) * h^2 # adjust for if eps(a) > h
        b  = a + da
    else
        a, b = promote(float(xs[1]), float(xs[2]))
    end
    fa, fb = f(a), f(b)
    cache!(cache, a, fa)
    cache!(cache, b, fb)
    secant(f, a, b, fa, fb; xtol, xrtol, maxiters, cache)
end

function secant(f, a::T, b::T, fa::T, fb::T; xtol = zero(T), xrtol = 8eps(T), maxiters = 1000, cache = nothing) where {T}
    # No function change; return arbitrary endpoint
    if fb == fa
        return (a, fa)
    end

    cnt = 0
    mbest = abs(fa) < abs(fb) ? a : b
    fbest = min(abs(fa), abs(fb))
    uatol = xtol / oneunit(xtol) * oneunit(real(a))
    adjustunit = oneunit(real(fb)) / oneunit(real(b))

    while cnt < maxiters
        m = b - (b - a) * fb / (fb - fa)
        fm = f(m)
        cache!(cache, m, fm)

        abs(fm) < abs(fbest) && ((mbest, fbest) = (m, fm))
        iszero(fm) && return (m, fm)
        isnan(fm) || isinf(fm) && return (mbest, fbest) # function failed; bail out
        abs(fm) <= adjustunit * max(uatol, abs(m) * xrtol) && return (m, fm)
        fm == fb && return (m, fm)

        a, b, fa, fb = b, m, fb, fm
        cnt += 1
    end

    return (mbest, fbest) # maxiters reached
end

function bisection_method(f, a::Number, b::Number; xtol = nothing, xrtol = nothing, ftol = nothing, cache = nothing, maxiters = 1000)
    x₁, x₂ = float.((a,b))
    y₁, y₂ = f(x₁), f(x₂)
    cache!(cache, x₁, y₁)
    cache!(cache, x₂, y₂)

    T = eltype(x₁)
    xtol = xtol === nothing ? zero(T) : abs(xtol)
    xrtol = xrtol === nothing ? zero(one(T)) : abs(xrtol)
    ftol = ftol === nothing ? zero(T) : abs(ftol)

    bisect(f, x₁, x₂, y₁, y₂; xtol, xrtol, ftol, cache, maxiters)
end

function bisect(f, x₁::T, x₂::T, y₁::T, y₂::T; xtol = zero(T), xrtol = zero(one(T)), ftol = zero(T), cache = nothing, maxiters = 1000) where {T}
    # No sign change; return arbitrary endpoint
    if y₁ * y₂ >= 0
        return (x₁, y₁)
    end

    if y₂ < 0
        x₁, x₂, y₁, y₂ = x₂, x₁, y₂, y₁
    end

    xₘ = (x₁ + x₂)/2
    yₘ = f(xₘ)
    cache!(cache, xₘ, yₘ)

    cnt = 1
    while cnt < maxiters
        if iszero(yₘ) || isnan(yₘ) || abs(x₁ - x₂) <= xtol + max(abs(x₁), abs(x₂)) * xrtol || abs(yₘ) <= ftol
            return (xₘ, yₘ)
        end

        if yₘ < 0
            x₁, y₁ = xₘ, yₘ
        else
            x₂, y₂ = xₘ, yₘ
        end

        xₘ = (x₁ + x₂)/2
        yₘ = f(xₘ)
        cache!(cache, xₘ, yₘ)

        cnt += 1
    end

    return (xₘ, yₘ)
end

function bracketing_interval(f, a, δ, dilate = 1; maxiters = 1000, cache = nothing)
    # Initialize cache
    fa = f(a)
    cache!(cache, a, fa)
    cnt = 0
    if fa > 0
        b = a - δ
        fb = f(b)
        cache!(cache, b, fb)
        δ *= dilate
        while fb > 0 && cnt < maxiters
            a, fa = b, fb
            b = a - δ
            fb = f(b)
            cache!(cache, b, fb)
            δ *= dilate
            cnt += 1
        end
    else
        b = a + δ
        fb = f(b)
        cache!(cache, b, fb)
        δ *= dilate
        while fb < 0 && cnt < maxiters
            a, fa = b, fb
            b = a + δ
            fb = f(b)
            cache!(cache, b, fb)
            δ *= dilate
            cnt += 1
        end
    end
    return a, b, fa, fb
end

####
#### L-curve method for choosing Tikhonov regularization parameter
####

struct NNLSLCurveRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSLCurveRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSLCurveRegProblem(C, d, m, n, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSLCurveRegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_lcurve(C::AbstractMatrix, d::AbstractVector)

Returns the regularized NNLS solution, X, of the equation

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||L x||_2^2
```

where ``L`` is the identity matrix and ``\\mu`` is chosen by locating the corner of the "L-curve".

Details of L-curve theory and the Generalized Cross-Validation (GCV) method can be found in:
[<NAME>., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580](https://doi.org/10.1137/1034115)

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Resulting increase in ``\\chi^2`` relative to unregularized (``\\mu = 0``) solution
"""
function lsqnonneg_lcurve(C, d)
    work = lsqnonneg_lcurve_work(C, d)
    lsqnonneg_lcurve!(work)
end
lsqnonneg_lcurve_work(C, d) = NNLSLCurveRegProblem(C, d)

function lsqnonneg_lcurve!(work::NNLSLCurveRegProblem{T,N}) where {T,N}
    # Compute the regularization using the L-curve method
    reset_cache!(work.nnls_work_smooth_cache)
    logmu_bounds = (T(-8), T(2))
    logmu_final = lcurve_corner(logmu_bounds...) do μ
        solve!(work.nnls_work_smooth_cache, μ)
        ξ = log(resnorm(get_cache(work.nnls_work_smooth_cache)))
        η = log(seminorm(get_cache(work.nnls_work_smooth_cache)))
        return SA{T}[ξ, η]
    end

    # Return the final regularized solution
    mu_final = exp(logmu_final)
    x_final = solve!(work.nnls_work_smooth_cache, mu_final)
    x_unreg = solve!(work.nnls_work)
    chi2factor_final = chi2(get_cache(work.nnls_work_smooth_cache)) / chi2(work.nnls_work)

    return (x = x_final, mu = mu_final, chi2factor = chi2factor_final)
end

"""
    lcurve_corner(f, xlow, xhigh)

Find the corner of the L-curve via curvature maximization using Algorithm 1 from:

<NAME> and <NAME>, “A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems”.
IOPSciNotes, vol. 1, no. 2, p. 025004, Aug. 2020, doi: 10.1088/2633-1357/abad0d
"""
function lcurve_corner(f, xlow::T = -8.0, xhigh::T = 2.0; xtol::T = 0.05, Ptol::T = 0.05, Ctol::T = 0.01, cache = nothing, refine = false, backtracking = true, verbose = false, kwargs...) where {T}
    # Initialize state
    msg(s, state) = verbose && (@info "$s: [x⃗, P⃗, C⃗] = "; display(hcat(state.x⃗, state.P⃗, [cache[x].C for x in state.x⃗])))
    cache === nothing && (cache = Dict{T, NamedTuple{(:P, :C), Tuple{SVector{2,T}, T}}}())
    state = LCurveCornerState(f, T(xlow), T(xhigh), cache; kwargs...)
    state_cache = backtracking ? [state] : nothing

    # Tolerances are relative to initial curve size
    Ptopleft, Pbottomright = state.P⃗[1], state.P⃗[4]
    Ptol = norm(Ptopleft - Pbottomright) * Ptol
    Ctol = norm(Ptopleft - Pbottomright) * Ctol

    # For very small regularization, points on the L-curve may be extremely close,
    # and curvature calculations will be highly ill-posed. Similarly, while it should
    # not typically occur, additionally check for points on the L-curve becoming
    # extremely close for large regularization. Points which don't satisfy `Pfilter`
    # are assigned curvature -Inf.
    Pfilter = P -> min(norm(P - Ptopleft), norm(P - Pbottomright)) > Ctol
    update_curvature!(state, cache; Pfilter)
    msg("Starting", state)

    while true
        if backtracking
            # Find state with minimum diameter which contains the current best estimate maximum curvature point
            (x, (P, C)), _, _ = mapfindmax(((x, (P, C)),) -> C, collect(cache))
            for s in state_cache
                if (s.x⃗[2] == x || s.x⃗[3] == x) && abs(s.x⃗[4] - s.x⃗[1]) <= abs(state.x⃗[4] - state.x⃗[1])
                    state = s
                end
            end
        end

        # Move state toward region of lower curvature
        if cache[state.x⃗[2]].C > cache[state.x⃗[3]].C
            state = move_left(f, state, cache)
            update_curvature!(state, cache; Pfilter)
            msg("C₂ > C₃; moved left", state)
        else
            state = move_right(f, state, cache)
            update_curvature!(state, cache; Pfilter)
            msg("C₃ ≥ C₂; moved right", state)
        end
        backtracking && push!(state_cache, state)
        is_converged(state; xtol, Ptol) && break
    end

    if refine
        x = refine!(f, state, cache; Pfilter)
        msg("Converged; refined solution", state)
    else
        x = cache[state.x⃗[2]].C > cache[state.x⃗[3]].C ? state.x⃗[2] : state.x⃗[3]
        msg("Converged", state)
    end

    return x
end

struct LCurveCornerState{T, V <: SVector{2,T}}
    x⃗::SVector{4,T} # grid of regularization parameters
    P⃗::SVector{4,V} # points (residual norm, solution seminorm) evaluated at x⃗
end
function LCurveCornerState(f, x₁::T, x₄::T, cache = nothing) where {T}
    x₂   = (T(φ) * x₁ + x₄) / (T(φ) + 1)
    x₃   = x₁ + (x₄ - x₂)
    x⃗    = SA{T}[x₁, x₂, x₃, x₄]
    P⃗    = f.(exp.(x⃗))
    foreach(i -> cache!(cache, x⃗[i], (P = P⃗[i], C = T(-Inf))), 1:4)
    return LCurveCornerState(x⃗, P⃗)
end

is_converged(state::LCurveCornerState; xtol, Ptol) = abs(state.x⃗[4] - state.x⃗[1]) < xtol || norm(state.P⃗[1] - state.P⃗[4]) < Ptol

function maybecall!(f, x::T, state::LCurveCornerState{T}, cache) where {T}
    for (_x, (_P, _C)) in cache
        x ≈ _x && return _P
    end
    P = f(exp(x)) # only recalculate if not cached
    cache!(cache, x, (P = P, C = T(-Inf)))
    return P
end

function move_left(f, state::LCurveCornerState{T}, cache) where {T}
    @unpack x⃗, P⃗ = state
    x⃗ = SA[x⃗[1], (T(φ) * x⃗[1] + x⃗[3]) / (T(φ) + 1), x⃗[2], x⃗[3]]
    P⃗ = SA[P⃗[1], maybecall!(f, x⃗[2], state, cache), P⃗[2], P⃗[3]] # only P⃗[2] is recalculated
    return LCurveCornerState(x⃗, P⃗)
end

function move_right(f, state::LCurveCornerState{T}, cache) where {T}
    @unpack x⃗, P⃗ = state
    x⃗ = SA[x⃗[2], x⃗[3], x⃗[2] + (x⃗[4] - x⃗[3]), x⃗[4]]
    P⃗ = SA[P⃗[2], P⃗[3], maybecall!(f, x⃗[3], state, cache), P⃗[4]] # only P⃗[3] is recalculated
    return LCurveCornerState(x⃗, P⃗)
end

function update_curvature!(state::LCurveCornerState{T}, cache; Pfilter = nothing, menger_curvature = true, global_search = false) where {T}
    @unpack x⃗, P⃗ = state
    Cfun(P₋, P, P₊) = menger_curvature ?
        menger(P₋, P, P₊) :
        360 - rad2deg(kahan_angle(P₋, P, P₊))
    for i in 1:4
        x, P, C = x⃗[i], P⃗[i], T(-Inf)
        if Pfilter === nothing || Pfilter(P)
            if global_search
                # Search for maximum curvature over all neighbours
                for (x₋, (P₋, _)) in cache, (x₊, (P₊, _)) in cache
                    (x₋ < x⃗[i] < x₊) && (C = max(C, Cfun(P₋, P, P₊)))
                end
            else
                # Compute curvature from nearest neighbours
                x₋, x₊ = T(-Inf), T(+Inf)
                P₋, P₊ = P, P
                for (_x, (_P, _)) in cache
                    (x₋ < _x < x ) && ((x₋, P₋) = (_x, _P))
                    (x  < _x < x₊) && ((x₊, P₊) = (_x, _P))
                end
                C = Cfun(P₋, P, P₊)
            end
        end
        cache[x] = (; P, C)
    end
    return state
end

function refine!(f, state::LCurveCornerState{T}, cache; Pfilter = nothing, analytical = false) where {T}
    # Fit spline to (negative) curvature estimates on grid and minimize over the spline
    if analytical
        x_opt, _, _ = maximize_curvature(state)
    else
        C_spl = make_spline(state.x⃗, [-cache[x].C for x in state.x⃗])
        x_opt, _ = spline_opt(C_spl)
    end
    maybecall!(f, x_opt, state, cache)
    update_curvature!(state, cache; Pfilter)
    (x_opt, (_, _)), _, _ = mapfindmax(((x, (P, C)),) -> C, collect(cache))
    return x_opt
end

function menger(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}}
    Δⱼₖ, Δₖₗ, Δₗⱼ = Pⱼ - Pₖ, Pₖ - Pₗ, Pₗ - Pⱼ
    P̄ⱼP̄ₖ, P̄ₖP̄ₗ, P̄ₗP̄ⱼ = Δⱼₖ ⋅ Δⱼₖ, Δₖₗ ⋅ Δₖₗ, Δₗⱼ ⋅ Δₗⱼ
    Cₖ = 2 * (Δⱼₖ × Δₖₗ) / √(P̄ⱼP̄ₖ * P̄ₖP̄ₗ * P̄ₗP̄ⱼ)
end

function menger(xⱼ::T, xₖ::T, xₗ::T, Pⱼ::V, Pₖ::V, Pₗ::V; interp_uniform = true, linear_deriv = true) where {T, V <: SVector{2,T}}
    if interp_uniform
        h  = min(abs(xₖ - xⱼ), abs(xₗ - xₖ)) / T(φ)
        h₋ = h₊ = h
        x₋, x₀, x₊ = xₖ - h, xₖ, xₖ + h
        P₀ = Pₖ
        P₋ = exp_interp.(x₋, xⱼ, xₖ, Pⱼ, Pₖ)
        P₊ = exp_interp.(x₊, xₖ, xₗ, Pₖ, Pₗ)
    else
        P₋, P₀, P₊ = Pⱼ, Pₖ, Pₗ
        x₋, x₀, x₊ = xⱼ, xₖ, xₗ
        h₋, h₊ = x₀ - x₋, x₊ - x₀
    end
    ξ₋, ξ₀, ξ₊ = P₋[1], P₀[1], P₊[1]
    η₋, η₀, η₊ = P₋[2], P₀[2], P₊[2]

    if linear_deriv
        ξ′  = (ξ₊ - ξ₋) / (h₊ + h₋)
        η′  = (η₊ - η₋) / (h₊ + h₋)
    else
        ξ′ = (h₋^2 * ξ₊ + (h₊ + h₋) * (h₊ - h₋) * ξ₀ - h₊^2 * ξ₋) / (h₊ * h₋ * (h₊ + h₋))
        η′ = (h₋^2 * η₊ + (h₊ + h₋) * (h₊ - h₋) * η₀ - h₊^2 * η₋) / (h₊ * h₋ * (h₊ + h₋))
    end

    ξ′′ = 2 * (h₋ * ξ₊ - (h₊ + h₋) * ξ₀ + h₊ * ξ₋) / (h₊ * h₋ * (h₊ + h₋))
    η′′ = 2 * (h₋ * η₊ - (h₊ + h₋) * η₀ + h₊ * η₋) / (h₊ * h₋ * (h₊ + h₋))

    return (ξ′ * η′′ - η′ * ξ′′) / √((ξ′^2 + η′^2)^3)
end

function menger(f; h = 1e-3)
    function menger_curvature_inner(x)
        fⱼ, fₖ, fₗ = f(x-h), f(x), f(x+h)
        Pⱼ, Pₖ, Pₗ = SA[x-h,fⱼ], SA[x,fₖ], SA[x+h,fₗ]
        menger(Pⱼ, Pₖ, Pₗ)
    end
end

function menger(x, y; h = 1e-3)
    function menger_curvature_inner(t)
        x₋, x₀, x₊ = x(t-h), x(t), x(t+h)
        y₋, y₀, y₊ = y(t-h), y(t), y(t+h)
        x′, x′′ = (x₊ - x₋) / 2h, (x₊ - 2x₀ + x₋) / h^2
        y′, y′′ = (y₊ - y₋) / 2h, (y₊ - 2y₀ + y₋) / h^2
        return (x′ * y′′ - y′ * x′′) / √((x′^2 + y′^2)^3)
    end
end

function menger(x::Dierckx.Spline1D, y::Dierckx.Spline1D)
    function menger_curvature_inner(t)
        x′  = Dierckx.derivative(x, t; nu = 1)
        x′′ = Dierckx.derivative(x, t; nu = 2)
        y′  = Dierckx.derivative(y, t; nu = 1)
        y′′ = Dierckx.derivative(y, t; nu = 2)
        return (x′ * y′′ - y′ * x′′) / √((x′^2 + y′^2)^3)
    end
end

function menger(y::Dierckx.Spline1D)
    function menger_curvature_inner(t)
        y′  = Dierckx.derivative(y, t; nu = 1)
        y′′ = Dierckx.derivative(y, t; nu = 2)
        return y′′ / √((1 + y′^2)^3)
    end
end

function maximize_curvature(state::LCurveCornerState{T}) where {T}
    # Maximize curvature and transform back from t-space to x-space
    @unpack x⃗, P⃗ = state
    t₁, t₂, t₃, t₄ = T(0), 1/T(3), 2/T(3), T(1)
    t_opt, P_opt, C_opt = maximize_curvature(P⃗...)
    x_opt =
        t₁ <= t_opt < t₂ ? lin_interp(t_opt, t₁, t₂, x⃗[1], x⃗[2]) :
        t₂ <= t_opt < t₃ ? lin_interp(t_opt, t₂, t₃, x⃗[2], x⃗[3]) :
        lin_interp(t_opt, t₃, t₄, x⃗[3], x⃗[4])
    return x_opt
end

function maximize_curvature(P₁::V, P₂::V, P₃::V, P₄::V; bezier = true) where {T, V <: SVector{2,T}}
    # Analytically maximize curvature of parametric cubic spline fit to data.
    #   see: https://cs.stackexchange.com/a/131032
    a, e = P₁[1], P₁[2]
    b, f = P₂[1], P₂[2]
    c, g = P₃[1], P₃[2]
    d, h = P₄[1], P₄[2]

    if bezier
        ξ = t -> a*(1-t)^3+3*b*(1-t)^2*t+3*c*(1-t)*t^2+d*t^3
        η = t -> e*(1-t)^3+3*f*(1-t)^2*t+3*g*(1-t)*t^2+h*t^3
        P = t -> SA{T}[ξ(t), η(t)]

        # In order to keep the equations in a more compact form, introduce the following substitutions:
        m = d-3*c+3*b-a
        n = c-2*b+a
        o = b-a
        p = h-3*g+3*f-e
        q = g-2*f+e
        r = f-e
    else
        ξ = t -> a*(t-1/3)*(t-2/3)*(t-1/1)/T((0/1-1/3)*(0/1-2/3)*(0/1-1/1)) + # (-9  * a * t^3)/2  + (+18 * a * t^2)/2 + (-11 * a * t)/2 + a
                 b*(t-0/1)*(t-2/3)*(t-1/1)/T((1/3-0/1)*(1/3-2/3)*(1/3-1/1)) + # (+27 * b * t^3)/2  + (-45 * b * t^2)/2 + (+18 * b * t)/2 +
                 c*(t-0/1)*(t-1/3)*(t-1/1)/T((2/3-0/1)*(2/3-1/3)*(2/3-1/1)) + # (-27 * c * t^3)/2  + (+36 * c * t^2)/2 + ( -9 * c * t)/2 +
                 d*(t-0/1)*(t-1/3)*(t-2/3)/T((1/1-0/1)*(1/1-1/3)*(1/1-2/3))   # (+9  * d * t^3)/2  + ( -9 * d * t^2)/2 + ( +2 * d * t)/2
        η = t -> e*(t-1/3)*(t-2/3)*(t-1/1)/T((0/1-1/3)*(0/1-2/3)*(0/1-1/1)) + # (-9  * e * t^3)/2  + (+18 * e * t^2)/2 + (-11 * e * t)/2 + e
                 f*(t-0/1)*(t-2/3)*(t-1/1)/T((1/3-0/1)*(1/3-2/3)*(1/3-1/1)) + # (+27 * f * t^3)/2  + (-45 * f * t^2)/2 + (+18 * f * t)/2 +
                 g*(t-0/1)*(t-1/3)*(t-1/1)/T((2/3-0/1)*(2/3-1/3)*(2/3-1/1)) + # (-27 * g * t^3)/2  + (+36 * g * t^2)/2 + ( -9 * g * t)/2 +
                 h*(t-0/1)*(t-1/3)*(t-2/3)/T((1/1-0/1)*(1/1-1/3)*(1/1-2/3))   # (+9  * h * t^3)/2  + ( -9 * h * t^2)/2 + ( +2 * h * t)/2
        P = t -> SA{T}[ξ(t), η(t)]

        # In order to keep the equations in a more compact form, introduce the following substitutions:
        m = ( -9 * a + 27 * b - 27 * c + 9 * d)/2
        n = (+18 * a - 45 * b + 36 * c - 9 * d)/6
        o = (-11 * a + 18 * b -  9 * c + 2 * d)/6
        p = ( -9 * e + 27 * f - 27 * g + 9 * h)/2
        q = (+18 * e - 45 * f + 36 * g - 9 * h)/6
        r = (-11 * e + 18 * f -  9 * g + 2 * h)/6
    end

    # This leads to the following simplified derivatives:
    ξ′  = t -> 3*(m*t^2+2*n*t+o)
    ξ′′ = t -> 6*(m*t+n)
    η′  = t -> 3*(p*t^2+2*q*t+r)
    η′′ = t -> 6*(p*t+q)

    # Curvature and its derivative:
    k  = t -> ((3*(m*t^2+2*n*t+o)) * (6*(p*t+q)) - (3*(p*t^2+2*q*t+r)) * (6*(m*t+n))) / ((3*(m*t^2+2*n*t+o))^2 + (3*(p*t^2+2*q*t+r))^2)^(3/2)
    k′ = t -> (-18*m*(p*t^2+2*q*t+r)+18*p*(m*t^2+2*n*t+o)-18*(m*t+n)*(2*p*t+2*q)+18*(2*m*t+2*n)*(p*t+q))/(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(3/2)-(3*(18*(p*t+q)*(m*t^2+2*n*t+o)-18*(m*t+n)*(p*t^2+2*q*t+r))*(18*(2*p*t+2*q)*(p*t^2+2*q*t+r)+18*(2*m*t+2*n)*(m*t^2+2*n*t+o)))/(2*(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(5/2))

    # Solve analytically
    coeffs = MVector{6,Complex{T}}(
        (1296*m*p^2+1296*m^3)*q-1296*n*p^3-1296*m^2*n*p,
        (1620*m*p^2+1620*m^3)*r+3240*m*p*q^2+(3240*m^2*n-3240*n*p^2)*q-1620*o*p^3+((-1620*m^2*o)-3240*m*n^2)*p,
        (5184*m*p*q+1296*n*p^2+6480*m^2*n)*r+1296*m*q^3-1296*n*p*q^2+((-6480*o*p^2)-1296*m^2*o+1296*m*n^2)*q+((-5184*m*n*o)-1296*n^3)*p,
        1296*m*p*r^2+(1944*m*q^2+6480*n*p*q-1296*o*p^2+1296*m^2*o+8424*m*n^2)*r-8424*o*p*q^2-6480*m*n*o*q+((-1296*m*o^2)-1944*n^2*o)*p,
        2592*n*p*r^2+(3888*n*q^2-2592*o*p*q+2592*m*n*o+3888*n^3)*r-3888*o*q^3+((-2592*m*o^2)-3888*n^2*o)*q,
        -324*m*r^3+(1944*n*q+324*o*p)*r^2+((-1944*o*q^2)-324*m*o^2+1944*n^2*o)*r-1944*n*o^2*q+324*o^3*p,
    )
    roots = MVector{6,Complex{T}}(undef)
    PolynomialRoots.roots5!(roots, coeffs, eps(T), true)

    # Check roots and return maximum
    t₁, t₄ = zero(T), one(T)
    k₁, k₄ = k(t₁), k(t₄)
    tmax, Pmax, kmax = k₁ > k₄ ? (t₁, P₁, k₁) : (t₄, P₄, k₄)
    for rᵢ in roots
        _t, _s = real(rᵢ), imag(rᵢ)
        !(_s ≈ 0) && continue # real roots only
        !(t₁ <= _t <= t₄) && continue # filter roots within range
        ((_k = k(_t)) > kmax) && ((tmax, Pmax, kmax) = (_t, P(_t), _k))
    end

    return tmax, Pmax, kmax
end

function directed_angle(v₁::V, v₂::V) where {T, V <: SVector{2,T}}
    α = atan(v₁[2], v₁[1]) - atan(v₂[2], v₂[1])
    α < 0 ? 2*T(π) + α : α
end
directed_angle(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}} = directed_angle(Pⱼ - Pₖ, Pₗ - Pₖ)

function kahan_angle(v₁::V, v₂::V) where {T, V <: SVector{2,T}}
    # Kahan's method for computing the angle between v₁ and v₂.
    #   see: https://scicomp.stackexchange.com/a/27694
    a, b, c = norm(v₁), norm(v₂), norm(v₁ - v₂)
    a, b = max(a,b), min(a,b)
    μ = b ≥ c ? c - (a - b) : (b - (a - c))
    num = ((a - b) + c) * max(μ, zero(μ))
    den = (a + (b + c)) * ((a - c) + b)
    α = 2 * atan(√(num / den))
    v₁ × v₂ > 0 ? 2*T(π) - α : α
end
kahan_angle(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}} = kahan_angle(Pⱼ - Pₖ, Pₗ - Pₖ)

####
#### GCV method for choosing Tikhonov regularization parameter
####

struct NNLSGCVRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    Aμ::Matrix{T}
    C_buf::Matrix{T}
    Ct_buf::Matrix{T}
    CtC_buf::Matrix{T}
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSGCVRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    Aμ = zeros(T, m, m)
    C_buf = zeros(T, m, n)
    Ct_buf = zeros(T, n, m)
    CtC_buf = zeros(T, n, n)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSGCVRegProblem(C, d, m, n, Aμ, C_buf, Ct_buf, CtC_buf, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSGCVRegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_gcv(C::AbstractMatrix, d::AbstractVector)

Returns the regularized NNLS solution, X, of the equation

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||L x||_2^2
```

where ``L`` is the identity matrix and ``\\mu`` is chosen by the Generalized Cross-Validation (GCV) method.

Details of the GCV method and L-curve theory can be can be found in:
[<NAME>., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580](https://doi.org/10.1137/1034115)

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Resulting increase in ``\\chi^2`` relative to unregularized (``\\mu = 0``) solution
"""
function lsqnonneg_gcv(C, d)
    work = lsqnonneg_gcv_work(C, d)
    lsqnonneg_gcv!(work)
end
lsqnonneg_gcv_work(C, d) = NNLSGCVRegProblem(C, d)

function lsqnonneg_gcv!(work::NNLSGCVRegProblem{T,N}) where {T,N}
    # Find μ by minimizing the function G(μ) (GCV method)
    reset_cache!(work.nnls_work_smooth_cache)

    # opt = NLopt.Opt(:LN_COBYLA, 1) # local, gradient-free, linear approximation of objective
    opt = NLopt.Opt(:LN_BOBYQA, 1) # local, gradient-free, quadratic approximation of objective
    opt.lower_bounds  = -8.0
    opt.upper_bounds  = 2.0
    opt.xtol_rel      = 0.05
    opt.min_objective = (logμ, ∇logμ) -> Float64(gcv!(work, logμ[1]))
    minf, minx, ret   = NLopt.optimize(opt, [-4.0])

    # Return the final regularized solution
    mu_final = exp(T(minx[1]))
    x_final  = solve!(work.nnls_work_smooth_cache, mu_final)
    x_unreg  = solve!(work.nnls_work)
    chi2factor_final = chi2(get_cache(work.nnls_work_smooth_cache)) / chi2(work.nnls_work)

    return (x = x_final, mu = mu_final, chi2factor = chi2factor_final)
end

# Implements equation (32) from:
# 
#   Analysis of Discrete Ill-Posed Problems by Means of the L-Curve
#   Hansen et al. 1992 (https://epubs.siam.org/doi/10.1137/1034115)
# 
# where here A = C, b = d, λ = μ, and L = identity.
function gcv!(work::NNLSGCVRegProblem, logμ; extract_subproblem = false)
    # Unpack buffers
    @unpack C, d, m, n, Aμ, C_buf, Ct_buf, CtC_buf = work

    # Solve regularized NNLS problem and record chi2 = ||C*X_reg - d||^2 which is returned
    μ = exp(logμ)
    solve!(work.nnls_work_smooth_cache, μ)
    χ² = chi2(get_cache(work.nnls_work_smooth_cache))
    x = solution(get_cache(work.nnls_work_smooth_cache))

    if extract_subproblem
        # Extract equivalent unconstrained least squares subproblem from NNLS problem
        # by extracting columns of C which correspond to nonzero components of x
        n′ = 0
        for (j,xⱼ) in enumerate(x)
            xⱼ ≈ 0 && continue
            n′ += 1
            @inbounds @simd ivdep for i in 1:m
                C_buf[i,n′] = C[i,j]
            end
        end
        C′ = reshape(uview(C_buf, 1:m*n′), m, n′)
        Ct′ = reshape(uview(Ct_buf, 1:n′*m), n′, m)
        CtC′ = reshape(uview(CtC_buf, 1:n′*n′), n′, n′)
    else
        # Use full matrix
        C′ = C
        Ct′ = Ct_buf
        CtC′ = CtC_buf
    end

    # Efficient compution of
    #   Aμ = C * (C'C + μ^2*I)^-1 * C'
    # where the matrices have sizes
    #   C: (m,n), Aμ: (m,m), Ct: (n,m), CtC: (n,n)
    mul!(CtC′, C′', C′) # C'C
    @inbounds @simd ivdep for i in 1:n
        CtC′[i,i] += μ^2 # C'C + μ^2*I
    end
    ldiv!(Ct′, cholesky!(Symmetric(CtC′)), C′') # (C'C + μ^2*I)^-1 * C'
    mul!(Aμ, C′, Ct′) # C * (C'C + μ^2*I)^-1 * C'

    # Return Generalized cross-validation. See equations 27 and 32 in
    #   <NAME>., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580
    #   https://doi.org/10.1137/1034115
    trace = m - tr(Aμ) # tr(I - Aμ) = m - tr(Aμ) for m x m matrix Aμ
    gcv = χ² / trace^2 # ||C*X_reg - d||^2 / tr(I - Aμ)^2

    return gcv
end

# Equation (27) from Hansen et al. 1992 (https://epubs.siam.org/doi/10.1137/1034115),
# specialized for L = identity:
# 
#   tr(I_m - A * (A'A + λ^2 * L'L)^-1 * A') = m - n + sum_i λ^2 / (γ_i^2 + λ^2)
# 
# where γ_i are the generalized singular values, which are equivalent to ordinary
# singular values when L = identity, and size(A) = (m,n).
function gcv_tr!(A, λ)
    m, n = size(A)
    γ  = svdvals!(A)
    γ .= λ.^2 ./ (γ.^2 .+ λ.^2)
    max(m - n, 0) + sum(γ)
end
gcv_tr(A, λ) = gcv_tr!(copy(A), λ)
gcv_tr_brute(A, λ) = tr(I - A * ((A'A + λ^2 * I) \ A'))



<|start_filename|>src/DECAES.jl<|end_filename|>
module DECAES

using LinearAlgebra, SpecialFunctions, Statistics, Random
import ArgParse, BangBang, Dierckx, DocStringExtensions, ForwardDiff, Logging, MAT, NIfTI, NLopt, Parameters, PolynomialRoots, ProgressMeter, SIMD, StaticArrays, TupleTools, UnPack, UnsafeArrays
using ArgParse: @add_arg_table!, ArgParseSettings, add_arg_group!, add_arg_table!, parse_args
using BangBang: setindex!!, setproperty!!, setproperties!!
using Base.MathConstants: φ
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDFIELDS, TYPEDSIGNATURES
using ForwardDiff: DiffResults
using Logging: ConsoleLogger, with_logger
using Parameters: @with_kw, @with_kw_noshow
using ProgressMeter: Progress, BarGlyphs, finish!, next!
using StaticArrays: FieldVector, SA, SArray, SVector, SMatrix, SizedVector, MVector
using SIMD: FloatingTypes, Vec, shufflevector
using UnPack: @unpack, @pack!
using UnsafeArrays: @uviews, uviews, uview

include("NNLS.jl")
using .NNLS

include("ParXRec.jl")
using .ParXRec

include("NormalHermiteSplines.jl")
using .NormalHermiteSplines

include("types.jl")
include("utils.jl")
include("splines.jl")
include("lsqnonneg.jl")
include("lsqnonneg_reg.jl")
include("EPGdecaycurve.jl")
include("T2mapSEcorr.jl")
include("T2partSEcorr.jl")
include("main.jl")

# Exported symbols
export MAT, NIfTI, ParXRec, load_image
export T2mapSEcorr, T2mapOptions, T2partSEcorr, T2partOptions
export EPGdecaycurve, EPGdecaycurve!, EPGdecaycurve_work
export lsqnonneg, lsqnonneg_chi2, lsqnonneg_gcv, lsqnonneg_lcurve, lcurve_corner
export main, julia_main

# include("precompile.jl")
# _precompile_()

end


<|start_filename|>src/ParXRec.jl<|end_filename|>
module ParXRec

using LightXML

export parxrec, parrec, xmlrec, SeriesInfo, ImageInfo, RecHdr

Base.@kwdef mutable struct SeriesInfo
    Patient_Name               :: String  = ""
    Examination_Name           :: String  = ""
    Protocol_Name              :: String  = ""
    Examination_Date           :: String  = ""      # TODO: Date
    Examination_Time           :: String  = ""      # TODO: Time
    Series_Data_Type           :: String  = ""
    Aquisition_Number          :: Int32   = 0
    Reconstruction_Number      :: Int32   = 0
    Scan_Duration              :: Float32 = 0
    Max_No_Phases              :: Int32   = 1
    Max_No_Echoes              :: Int32   = 1
    Max_No_Slices              :: Int32   = 1
    Max_No_Dynamics            :: Int32   = 1
    Max_No_Mixes               :: Int16   = 1
    Max_No_B_Values            :: Int32   = 1
    Max_No_Gradient_Orients    :: Int32   = 1
    No_Label_Types             :: Int32   = 0
    Patient_Position           :: String  = ""
    Preparation_Direction      :: String  = ""
    Technique                  :: String  = ""
    Scan_Resolution_X          :: Int16   = 0
    Scan_Resolution_Y          :: Int32   = 0
    Scan_Mode                  :: String  = ""
    Repetition_Times           :: NTuple{2,Float32} = (0,0)
    FOV_AP                     :: Float32 = 0
    FOV_FH                     :: Float32 = 0
    FOV_RL                     :: Float32 = 0
    Water_Fat_Shift            :: Float32 = 0
    Angulation_AP              :: Float32 = 0
    Angulation_FH              :: Float32 = 0
    Angulation_RL              :: Float32 = 0
    Off_Center_AP              :: Float32 = 0
    Off_Center_FH              :: Float32 = 0
    Off_Center_RL              :: Float32 = 0
    Flow_Compensation          :: Bool    = false
    Presaturation              :: Bool    = false
    Phase_Encoding_Velocity    :: NTuple{3,Float32} = (0,0,0)
    MTC                        :: Bool    = false
    SPIR                       :: Bool    = false
    EPI_factor                 :: Int32   = 0
    Dynamic_Scan               :: Bool    = false
    Diffusion                  :: Bool    = false
    Diffusion_Echo_Time        :: Float32 = 0
    PhotometricInterpretation  :: String  = ""
end

Base.@kwdef mutable struct ImageKey
    Slice       :: Int32  = 1
    Echo        :: Int32  = 1
    Dynamic     :: Int32  = 1
    Phase       :: Int32  = 1
    BValue      :: Int32  = 1
    Grad_Orient :: Int32  = 1
    Label_Type  :: String = ""      # Enumeration
    Type        :: String = ""      # Enumeration
    Sequence    :: String = ""      # Enumeration
    Index       :: Int32  = 0
end

Base.@kwdef mutable struct ImageInfo
    Key                        :: ImageKey
    Pixel_Size                 :: UInt16  = 0
    Scan_Percentage            :: Float64 = 0
    Resolution_X               :: UInt16  = 0
    Resolution_Y               :: UInt16  = 0
    Rescale_Intercept          :: Float64 = 0
    Rescale_Slope              :: Float64 = 0
    Scale_Slope                :: Float32 = 0
    Window_Center              :: Float64 = 0
    Window_Width               :: Float64 = 0
    Slice_Thickness            :: Float64 = 0
    Slice_Gap                  :: Float64 = 0
    Display_Orientation        :: String  = ""      # Enumeration
    fMRI_Status_Indication     :: Int16   = 0
    Image_Type_Ed_Es           :: String  = ""      # Enumeration
    Pixel_Spacing              :: NTuple{2,Float64} = (0,0)
    Echo_Time                  :: Float64 = 0
    Dyn_Scan_Begin_Time        :: Float32 = 0
    Trigger_Time               :: Float64 = 0
    Diffusion_B_Factor         :: Float32 = 0
    No_Averages                :: Float64 = 0
    Image_Flip_Angle           :: Float64 = 0
    Cardiac_Frequency          :: Int32   = 0
    Min_RR_Interval            :: Int32   = 0
    Max_RR_Interval            :: Int32   = 0
    TURBO_Factor               :: Int32   = 0
    Inversion_Delay            :: Float64 = 0
    Contrast_Type              :: String  = ""
    Diffusion_Anisotropy_Type  :: String  = ""
    Diffusion_AP               :: Float32 = 0
    Diffusion_FH               :: Float32 = 0
    Diffusion_RL               :: Float32 = 0
    Angulation_AP              :: Float64 = 0
    Angulation_FH              :: Float64 = 0
    Angulation_RL              :: Float64 = 0
    Offcenter_AP               :: Float64 = 0
    Offcenter_FH               :: Float64 = 0
    Offcenter_RL               :: Float64 = 0
    Slice_Orientation          :: String  = ""      # Enumeration
    Image_Planar_Configuration :: UInt16  = 0
    Samples_Per_Pixel          :: UInt16  = 0
end

mutable struct RecHdr
    series :: SeriesInfo
    images :: Vector{ImageInfo}
end

####
#### ParXRec
####

parxrec(path::AbstractString; load_data::Bool = true) =
    parxrec(Float64, path, load_data = load_data)

function parxrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data)
    ext = splitext(hdrpath)[2]
    hdr = lowercase(ext) == ".par" ? readpar(hdrpath) : readxml(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

parrec(path::AbstractString; load_data::Bool = true) =
    parrec(Float64, path, load_data = load_data)

function parrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data, ".PAR")
    hdr = readpar(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

xmlrec(path::AbstractString; load_data::Bool = true) =
    xmlrec(Float64, path, load_data = load_data)

function xmlrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data, ".XML")
    hdr = readxml(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

function parxrec_paths(path::AbstractString, load_data::Bool, ext = "")
    isfile(path) || error("file $path does not exist")
    path1, ext1 = splitext(path)

    if lowercase(ext1) == ".xml" || lowercase(ext1) == ".par"
        ext == "" || lowercase(ext1) == lowercase(ext) ||
            error("wrong file extension: $path (expected $ext)")

        hdrpath = path
        recpath = load_data ? (
            isfile(path1 * ".REC") ? path1 * ".REC" :
            isfile(path1 * ".rec") ? path1 * ".rec" :
            error(".REC file $path1 does not exist")
        ) : ""

    elseif lowercase(ext1) == ".rec"
        recpath = path
        hdrpath = ext == "" ? (
            isfile(path1 * ".PAR") ? path1 * ".PAR" :
            isfile(path1 * ".XML") ? path1 * ".XML" :
            isfile(path1 * ".par") ? path1 * ".par" :
            isfile(path1 * ".xml") ? path1 * ".xml" :
            error(".PAR/.XML file $path1 does not exist")

        ) : lowercase(ext) == ".par" ? (
            isfile(path1 * ".PAR") ? path1 * ".PAR" :
            isfile(path1 * ".par") ? path1 * ".par" :
            error(".PAR file $path1 does not exist")

        ) : lowercase(ext) == ".xml" ? (
            isfile(path1 * ".XML") ? path1 * ".XML" :
            isfile(path1 * ".xml") ? path1 * ".xml" :
            error(".XML file $path1 does not exist")

        ) : error("file extension $ext is not a valid Philips Rec extension")

    else
        error("file $path does not have valid Philips Rec extension")
    end

    return hdrpath, recpath
end

####
#### .REC
####

readrec(path::AbstractString, hdr::RecHdr) =
    readrec(Float64, path, hdr)

function readrec(T::DataType, path::AbstractString, hdr::RecHdr)
    validate(hdr)
    hdr = hdr.images

    prec = Dict{Int,DataType}(
        8  => UInt8,
        16 => UInt16,
        32 => Float32,
        64 => Float64,
    )

    nx       = Int(hdr[1].Resolution_X)
    ny       = Int(hdr[2].Resolution_X)
    nz       = maximum((i -> i.Key.Slice).(hdr))
    echoes   = unique((i -> i.Key.Echo).(hdr))
    dynamics = unique((i -> i.Key.Dynamic).(hdr))
    phases   = unique((i -> i.Key.Phase).(hdr))
    bvals    = unique((i -> i.Key.BValue).(hdr))
    gorients = unique((i -> i.Key.Grad_Orient).(hdr))
    types    = Dict(t => i for (i, t) in enumerate(unique((i -> i.Key.Type).(hdr))))
    labels   = Dict(l => i for (i, l) in enumerate(unique((i -> i.Key.Label_Type).(hdr))))
    seqs     = Dict(s => i for (i, s) in enumerate(unique((i -> i.Key.Sequence).(hdr))))

    if length(unique((i -> i.Pixel_Size).(hdr))) == 1
        sz = (nx * ny * hdr[1].Pixel_Size÷8 * length(hdr)) - filesize(path)
        sz < 0 && @warn "REC file larger than expected from hdr"
        sz > 0 && @warn "REC file smaller than expected from hdr"
    end

    data = zeros(T,
        nx, ny, nz, length.((echoes, dynamics, phases, bvals, gorients, labels, types, seqs))...
    )

    open(path, "r") do io
        for (i, img) in enumerate(hdr)
            ps = Int(img.Pixel_Size)

            seek(io, img.Key.Index * (nx * ny * ps÷8))
            slice = read!(io, Array{prec[ps]}(undef, nx, ny))

            ss  = inv(img.Scale_Slope)
            ri  = ss * img.Rescale_Intercept / img.Rescale_Slope

            i3  = img.Key.Slice
            i4  = img.Key.Echo
            i5  = img.Key.Dynamic
            i6  = img.Key.Phase
            i7  = img.Key.BValue
            i8  = img.Key.Grad_Orient
            i9  = labels[img.Key.Label_Type]
            i10 = types[img.Key.Type]
            i11 = seqs[img.Key.Sequence]

            @. data[:,:,i3,i4,i5,i6,i7,i8,i9,i10,i11] = ss*slice + ri
        end
    end

    return squeeze(data)
end

####
#### .XML
####

const XMLREC_HEADER           = "PRIDE_V5"
const XMLREC_SERIES_HEADER    = "Series_Info"
const XMLREC_IMAGE_ARR_HEADER = "Image_Array"
const XMLREC_IMAGE_HEADER     = "Image_Info"
const XMLREC_IMAGE_KEY_HEADER = "Key"
const XMLREC_ATTRIB_HEADER    = "Attribute"

function readxml(path::AbstractString)
    doc = parse_file(path)

    xroot = root(doc)
    name(xroot) == XMLREC_HEADER || @warn "unknown XML/REC header: $(name(xroot))"

    # Series Info
    cur = find_element(xroot, XMLREC_SERIES_HEADER)
    cur !== nothing || error("<Series_Info> tag not found")

    Ts = NamedTuple{(fieldnames(SeriesInfo)...,)}(SeriesInfo.types)
    series = SeriesInfo()

    for e in cur[XMLREC_ATTRIB_HEADER]
        key, val = _tryparse(Ts, e)
        val !== nothing ?
            setfield!(series, key, val) :
            @warn "unknown series attribute: $key. Skipping..."
    end

    # Image Info
    cur = find_element(xroot, XMLREC_IMAGE_ARR_HEADER)
    cur !== nothing || error("<Image_Array> tag not found")

    Tk = NamedTuple{(fieldnames(ImageKey)...,)}(ImageKey.types)
    Ti = NamedTuple{(fieldnames(ImageInfo)...,)}(ImageInfo.types)
    images = Vector{ImageInfo}()

    for cur in cur[XMLREC_IMAGE_HEADER]
        k = find_element(cur, XMLREC_IMAGE_KEY_HEADER)
        k !== nothing || error("<Key> tag not found")

        imagekey = ImageKey()
        for e in k[XMLREC_ATTRIB_HEADER]
            key, val = _tryparse(Tk, e)
            val !== nothing ?
                setfield!(imagekey, key, val) :
                @warn "unknown key attribute: $key. Skipping..."
        end

        image = ImageInfo(Key = imagekey)
        for e in cur[XMLREC_ATTRIB_HEADER]
            key, val = _tryparse(Ti, e)
            val !== nothing ?
                setfield!(image, key, val) :
                @warn "unknown image attribute: $key. Skipping..."
        end

        push!(images, image)
    end

    free(doc)

    return RecHdr(series, images)
end

####
#### .PAR
####

const PARSERIES = Dict{Symbol,String}(
    :Patient_Name              => "Patient name",
    :Examination_Name          => "Examination name",
    :Protocol_Name             => "Protocol name",
    :Examination_Date          => "Examination date/time",
    :Examination_Time          => "Examination date/time",
    :Series_Data_Type          => "Series Type",
    :Aquisition_Number         => "Acquisition nr",
    :Reconstruction_Number     => "Reconstruction nr",
    :Scan_Duration             => "Scan Duration",
    :Max_No_Phases             => "Max. number of cardiac phases",
    :Max_No_Echoes             => "Max. number of echoes",
    :Max_No_Slices             => "Max. number of slices",
    :Max_No_Dynamics           => "Max. number of dynamics",
    :Max_No_Mixes              => "Max. number of mixes",
    :Max_No_B_Values           => "Max. number of diffusion",
    :Max_No_Gradient_Orients   => "Max. number of gradient",
    :No_Label_Types            => "Number of label types",
    :Patient_Position          => "Patient position",
    :Preparation_Direction     => "Preparation direction",
    :Technique                 => "Technique",
    :Scan_Resolution_X         => "Scan resolution",
    :Scan_Resolution_Y         => "Scan resolution",
    :Scan_Mode                 => "Scan mode",
    :Repetition_Times          => "Repetition time",
    :FOV_AP                    => "FOV",
    :FOV_FH                    => "FOV",
    :FOV_RL                    => "FOV",
    :Water_Fat_Shift           => "Water Fat",
    :Angulation_AP             => "Angulation",
    :Angulation_FH             => "Angulation",
    :Angulation_RL             => "Angulation",
    :Off_Center_AP             => "Off Centre",
    :Off_Center_FH             => "Off Centre",
    :Off_Center_RL             => "Off Centre",
    :Flow_Compensation         => "Flow",
    :Presaturation             => "Presaturation",
    :Phase_Encoding_Velocity   => "Phase encoding",
    :MTC                       => "MTC",
    :SPIR                      => "SPIR",
    :EPI_factor                => "EPI factor",
    :Dynamic_Scan              => "Dynamic scan",
    :Diffusion                 => "Diffusion   ",
    :Diffusion_Echo_Time       => "Diffusion echo time",
)

const PARKEY = Dict{Symbol,Int}(
    :Slice       => 1,
    :Echo        => 2,
    :Dynamic     => 3,
    :Phase       => 4,
    :BValue      => 42,
    :Grad_Orient => 43,
    :Label_Type  => 49,
    :Type        => 5,
    :Sequence    => 6,
    :Index       => 7,
)

const PARIMAGE = Dict{Symbol,Int}(
    :Pixel_Size                 => 8,
    :Scan_Percentage            => 9,
    :Resolution_X               => 10,
    :Resolution_Y               => 11,
    :Rescale_Intercept          => 12,
    :Rescale_Slope              => 13,
    :Scale_Slope                => 14,
    :Window_Center              => 15,
    :Window_Width               => 16,
    :Slice_Thickness            => 23,
    :Slice_Gap                  => 24,
    :Display_Orientation        => 25,
    :fMRI_Status_Indication     => 27,
    :Image_Type_Ed_Es           => 28,
    :Pixel_Spacing              => 29,  # [29, 30]
    :Echo_Time                  => 31,
    :Dyn_Scan_Begin_Time        => 32,
    :Trigger_Time               => 33,
    :Diffusion_B_Factor         => 34,
    :No_Averages                => 35,
    :Image_Flip_Angle           => 36,
    :Cardiac_Frequency          => 37,
    :Min_RR_Interval            => 38,
    :Max_RR_Interval            => 39,
    :TURBO_Factor               => 40,
    :Inversion_Delay            => 41,
    :Contrast_Type              => 44,
    :Diffusion_Anisotropy_Type  => 45,
    :Diffusion_AP               => 46,
    :Diffusion_FH               => 47,
    :Diffusion_RL               => 48,
    :Angulation_AP              => 17,
    :Angulation_FH              => 18,
    :Angulation_RL              => 19,
    :Offcenter_AP               => 20,
    :Offcenter_FH               => 21,
    :Offcenter_RL               => 22,
    :Slice_Orientation          => 26,
)

function readpar(path::AbstractString)
    doc = read(path, String)

    # Series Info / General Information
    Ts = NamedTuple{(fieldnames(SeriesInfo)...,)}(SeriesInfo.types)
    series = SeriesInfo()

    # lines starting with `.`
    buf = [strip.(split(x.match, ": ")) for x in eachmatch(r"(?<=\n\.).*?\n", doc)];

    for (k, v) in PARSERIES
        str = buf[findfirst(a -> occursin(v, a[1]), buf)][2]
        val =
            k == :Examination_Date ?
                strip(split(str, "/")[1]) |> string :
            k == :Examination_Time ?
                strip(split(str, "/")[2]) |> string :
            k ∈ [:Scan_Resolution_X, :FOV_AP, :Angulation_AP, :Off_Center_AP] ?
                _parse(Ts[k], split(str)[1]) :
            k ∈ [:Scan_Resolution_Y, :FOV_FH, :Angulation_FH, :Off_Center_FH] ?
                _parse(Ts[k], split(str)[2]) :
            k ∈ [:FOV_RL, :Angulation_RL, :Off_Center_RL] ?
                _parse(Ts[k], split(str)[3]) :
            # else
                _parse(Ts[k], str)

        setfield!(series, k, val)
    end

    # Image Info
    Tk = NamedTuple{(fieldnames(ImageKey)...,)}(ImageKey.types)
    Ti = NamedTuple{(fieldnames(ImageInfo)...,)}(ImageInfo.types)
    images = Vector{ImageInfo}()

    # lines starting with ` ` or digit
    buf = [split(x.match) for x in eachmatch(r"(?<=\n)( |\d+).*?\n", doc)];

    for line in buf
        imagekey = ImageKey()
        for (k, v) in PARKEY
            setfield!(imagekey, k, _parse(Tk[k], line[v]))
        end
        image = ImageInfo(Key = imagekey)
        for (k, v) in PARIMAGE
            k == :Pixel_Spacing ?
                setfield!(image, k, _parse(Ti[k], line[v]*" "*line[v+1])) :
                setfield!(image, k, _parse(Ti[k], line[v]))
        end
        push!(images, image)
    end

    return RecHdr(series, images)
end

####
#### Utilities
####

function validate(hdr::RecHdr)
    series = hdr.series
    images = hdr.images

    keys = (i -> i.Key).(images)
    keymax = (v) -> maximum(getfield.(keys, v))
    keylength = (v) -> length(unique(getfield.(keys, v)))
    imglength = (v) -> length(unique(getfield.(images, v)))

    vals = [
        (series.Max_No_Slices, :Slice, "slices"),
        (series.Max_No_Echoes, :Echo, "echoes"),
        (series.Max_No_Dynamics, :Dynamic, "dynamics"),
        (series.Max_No_Phases, :Phase, "phases"),
        (series.Max_No_B_Values, :BValue, "bvalues"),
        (series.Max_No_Gradient_Orients, :Grad_Orient, "gradient orientations"),
    ]

    # Check whether data size in series info is consistent with image info
    wstr1 = (f, s, i) ->
        "maximum number of $f not matching: $s (series info) vs. $i (image info)"

    for (v1,v2,v3) in vals
        v1 == keymax(v2) || @warn wstr1(v3, v1, keymax(v2))
    end

    series.No_Label_Types + 1 == keylength(:Label_Type) ||
        @warn wstr1("label types", series.No_Label_Types, keylength(:Label_Type))

    # Non-exhaustive image info check
    wstr2 = (f, m, l) ->
        "inconsistent number of $f in image info: $m (maximum) vs. $l (length(unique))"

    for (_, v1, v2) in vals
        keymax(v1) == keylength(v1) ||
            @warn wstr2(v2, keymax(v1), keylength(v1))
    end

    keymax(:Index) + 1 == keylength(:Index) ||
        @warn wstr2("indices", keymax(:Index), keylength(:Index))

    imglength(:Pixel_Size) == 1 ||
        @warn "multiple values for `Pixel Size` found"

    imglength(:Slice_Thickness) == 1 ||
        @warn "multiple values for `Slice Thickness` found"

    imglength(:Slice_Gap) == 1 ||
        @warn "multiple values for `Slice Gap` found"

    length(unique(first.(getfield.(images, :Pixel_Spacing)))) == 1 ||
        @warn "multiple values for `Pixel Spacing` (x) found"

    length(unique(last.(getfield.(images, :Pixel_Spacing)))) == 1 ||
        @warn "multiple values for `Pixel Spacing` (y) found"

    imglength(:Resolution_X) == 1 ||
        error("multiple values for `Resolution X` found")

    imglength(:Resolution_Y) == 1 ||
        error("multiple values for `Resolution Y` found")

    return nothing
end

function squeeze(x::Array{T}) where {T}
    N = count(n -> n > 1, size(x))
    dims = ([d for d in 1:ndims(x) if size(x, d) == 1]...,)
    return dropdims(x, dims = dims) :: Array{T,N}
end

_parse(::Type{String}, val::AbstractString) = string(val)

_parse(::Type{T}, val::AbstractString) where {T<:Tuple} =
    (_parse.(T.types, split(val))...,)

_parse(::Type{Bool}, val::AbstractString) =
    (val == "1" || uppercase(val) == "Y") ? true :
    (val == "0" || uppercase(val) == "N") ? false : _parse(Bool, val)

_parse(type, str; kwargs...) = Base.parse(type, str; kwargs...) # fallback

function _tryparse(T::NamedTuple, e::XMLElement)
    key = attribute(e, "Name") |> x -> replace(strip(x), " " => "_") |> Symbol
    val = haskey(T, key) ? _parse(T[key], strip(content(e))) : nothing
    return key, val
end

_tryparse(type, str; kwargs...) = Base.tryparse(type, str; kwargs...) # fallback

end # module


<|start_filename|>docs/src/assets/logo.jl<|end_filename|>
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c2274306-49b5-45c0-ae4b-a3b8dc3a3835
begin
	using DECAES
	using Random
	import PlutoUI

	using CairoMakie
	CairoMakie.activate!(type = "png")

	using Makie.Colors
	Makie.set_theme!(resolution=(1600, 900))
end;

# ╔═╡ 9d4a6681-2b40-442f-8cc8-4ae444bc6215
md"""
ETL = $(@bind ETL_slider PlutoUI.Slider(4:32, default=32, show_value=true))

FA = $(@bind flip_angle_slider PlutoUI.Slider(range(50,180,length=100), default=150, show_value=true))

β = $(@bind refcon_slider PlutoUI.Slider(range(50,180,length=100), default=180, show_value=true))

T21 = $(@bind T21_slider PlutoUI.Slider(range(10e-3, 100e-3, length=100), default=30e-3, show_value=true))

T22 = $(@bind T22_slider PlutoUI.Slider(range(100e-3, 250e-3, length=100), default=150e-3, show_value=true))

MWF = $(@bind MWF_slider PlutoUI.NumberField(0.1:0.01:0.4, default=0.20))
"""

# ╔═╡ dfac902d-bb87-4fea-a898-d6792e9b793b
let
	sliders_md = md"""
	base height = $(@bind base_height PlutoUI.NumberField(5:5:5000, default=750))
	base width = $(@bind base_width PlutoUI.NumberField(5:5:5000, default=1610))

	E1 size = $(@bind E1_size PlutoUI.NumberField(0.01:0.01:1.0, default=0.69))
	E1 pos  = $(@bind E1_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=0.23))

	C size = $(@bind C_size PlutoUI.NumberField(0.01:0.01:1.0, default=0.56))
	C pos  = $(@bind C_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=0.41))

	A size = $(@bind A_size PlutoUI.NumberField(0.01:0.01:1.0, default=0.49))
	A pos  = $(@bind A_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=0.59))

	E2 size = $(@bind E2_size PlutoUI.NumberField(0.01:0.01:1.0, default=0.41))
	E2 pos  = $(@bind E2_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=0.75))

	S size = $(@bind S_size PlutoUI.NumberField(0.01:0.01:1.0, default=0.35))
	S pos  = $(@bind S_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=0.87))

	font = $(@bind decaes_font PlutoUI.Select(
		[
			"Dejavu Sans (default)",
			"Courier Mono",
			"Courier Bold",
			"Purisa Bold",
			"FreeMono",
			"Latin Modern Sans Demi Cond",
			"Liberation Mono",
			"Sawasdee Bold",
			"Carlito",
			"Chilanka",
			"Dyuthi",
			"Latin Modern Sans",
		];
		default = "Chilanka",
	))
	"""

	# Print current variables to terminal
	println("""
	base height = \$(@bind base_height PlutoUI.NumberField(5:5:5000, default=$(base_height)))
	base width = \$(@bind base_width PlutoUI.NumberField(5:5:5000, default=$(base_width)))

	E1 size = \$(@bind E1_size PlutoUI.NumberField(0.01:0.01:1.0, default=$(E1_size)))
	E1 pos  = \$(@bind E1_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=$(E1_pos)))

	C size = \$(@bind C_size PlutoUI.NumberField(0.01:0.01:1.0, default=$(C_size)))
	C pos  = \$(@bind C_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=$(C_pos)))

	A size = \$(@bind A_size PlutoUI.NumberField(0.01:0.01:1.0, default=$(A_size)))
	A pos  = \$(@bind A_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=$(A_pos)))

	E2 size = \$(@bind E2_size PlutoUI.NumberField(0.01:0.01:1.0, default=$(E2_size)))
	E2 pos  = \$(@bind E2_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=$(E2_pos)))

	S size = \$(@bind S_size PlutoUI.NumberField(0.01:0.01:1.0, default=$(S_size)))
	S pos  = \$(@bind S_pos  PlutoUI.NumberField(0.01:0.005:1.0, default=$(S_pos)))
	""")

	sliders_md
end

# ╔═╡ 5cf6be3c-fe58-42a8-840e-8a21a5a81584
begin
	# Global variables from sliders
	textsize = base_height .* [1.0, E1_size, C_size, A_size, E2_size, S_size]
	xposition = 20 .+ base_width .* [0, E1_pos, C_pos, A_pos, E2_pos, S_pos]
	yposition = 20 .+ [100, 50, 25, 15, 5, 0]
	textsize .= ceil.(Int, textsize)
	xposition .= ceil.(Int, xposition)
	yposition .= ceil.(Int, yposition)

	# Julia logo colors
	JLcolors = Makie.Colors.JULIA_LOGO_COLORS
end;

# ╔═╡ 15b31eb3-8ddc-4fab-9bda-41cdbe8bc710
function plot_decaes_text!(ax; textsize, xposition, yposition)
	text!(ax,
		["D", "E", "C", "A", "E", "S"],
		position = tuple.(xposition ./ 1600, yposition ./ 900),
		textsize = textsize,
		align = (:left, :baseline),
		font = decaes_font,
	)
end

# ╔═╡ 446a6ee0-b25b-47ef-b4fc-af59f30d27fd
md"""
EPG decay curve
"""

# ╔═╡ 118177e3-2b40-4cea-bba6-ab48ba492fbd
function plot_epg_curve!(ax, echo_times, epg_curve)
	epg_line = lines!(ax, echo_times, epg_curve;
		linewidth = 14,
		color = JLcolors.purple,
	)
	epg_band = band!(
		ax, echo_times, fill(0, length(echo_times)), epg_curve;
		color = (JLcolors.purple, 0.12),
	)
	return epg_line, epg_band
end;

# ╔═╡ 499d8be9-ca07-40c9-9a22-79bb99ed58fa
md"""
Inset T2 distribution figure:
"""

# ╔═╡ f217fde8-08b8-4f33-8149-69f900847bd3
function plot_t2_dist!(ax, t2maps, t2dist)
	ax.backgroundcolor = :white
	hidespines!(ax)
	hidedecorations!(ax)
	logT2 = t2maps["t2times"] .|> log10

	I = findlast(t2maps["t2times"] .<= 50e-3)
	t2_bands = Any[]
	t2_lines = Any[]
	pushinto!(x) = p -> push!(x, p)
	for (inds, color) in [
			(1:I, JLcolors.red),
			(I:length(t2dist), JLcolors.green),
		]
		band!(ax,
			logT2[inds], fill(0.0,length(logT2[inds])), t2dist[inds];
			color = (color, 0.35),
		) |> pushinto!(t2_bands)
		lines!(ax, logT2[inds], t2dist[inds],
			linewidth = 14, color = color,
		) |> pushinto!(t2_lines)
	end

	return t2_bands, t2_lines
end;

# ╔═╡ 0e7420c8-bfd6-4c91-a7bd-6c205eb24452
md"""
Compute EPG decay curve + corresponding T2 distribution:
"""

# ╔═╡ 02d229c2-e2c0-4a36-bc6a-163778be85b9
function epg_decay_curve(ETL, α, TE, T21, T22, T1, β, MWF, SNR)
	x = (0:ETL-1)./(ETL-1)
	y =    MWF  .* DECAES.EPGdecaycurve(ETL, α, TE, T21, T1, β) .+
		(1-MWF) .* DECAES.EPGdecaycurve(ETL, α, TE, T22, T1, β)
	y .*= 0.995 / maximum(y)
	y = Vector{Float64}(y)

	noise(snr, seed) = 10^(-snr/20) .* randn(MersenneTwister(seed), size(y))
	y .= sqrt.((y .+ noise(SNR, 0)).^2 .+ noise(SNR, 1).^2)

	return x, y
end;

# ╔═╡ db1d1148-7bf7-4a1d-b568-52b55a8ec1f3
function t2_distribution(y::AbstractVector, nT2, TE, T2Range, Reg)
	t2maps, t2dist = T2mapSEcorr(
		reshape(y, 1, 1, 1, :);
		nT2, TE, T2Range, Reg, Silent = true,
	)
	t2dist ./= maximum(t2dist)
	return t2maps, t2dist
end;

# ╔═╡ 7cd9afdf-6675-4416-aec7-9eba65eafbe2
function decaes_logo(;
		T1, TE, SNR, nT2, T2Range, Reg, ETL, α, β, T21, T22, MWF,
		save_figs = false,
	)
	fig = Figure(
		resolution = (1600, 900),
	)

	ax = Axis(
		fig[1, 1],
		aspect = 16/9,
		grid = false,
		ticks = false,
		backgroundcolor = :white,
	)
	xlims!(ax, 0, 1)
	ylims!(ax, 0, 1.0)
	hidespines!(ax)
	# hidespines!(ax, :t, :r) # only top and right
	hidedecorations!(ax)

	# DECAES text
	txt = plot_decaes_text!(ax; textsize, xposition, yposition)

	# Plot EPG decay curve
	echo_times, epg_curve = epg_decay_curve(
		ETL, α, TE, T21, T22, T1, β, MWF, SNR
	)
	epg_line, epg_band = plot_epg_curve!(ax, echo_times, epg_curve)

	# Inset T2 distbn
	ax_sub = Axis(fig, bbox = BBox(800, 1600, 450, 900))
	t2maps, t2dist = t2_distribution(epg_curve, nT2, TE, T2Range, Reg)
	t2_bands, t2_lines = plot_t2_dist!(ax_sub, t2maps, t2dist)

	# Save figures
	if save_figs
		save("logo.png", fig)
		save("logo.svg", fig)
	end

	(; fig, ax, ax_sub, txt, epg_line, epg_band, t2_bands, t2_lines)
end;

# ╔═╡ e55b552a-00a1-44d3-9d0b-dda4ceac1f59
let
	fig, = decaes_logo(;
		# DECAES settings
		T1 = 1.0,
		TE = 10e-3,
		SNR = 60,
		nT2 = 80,
		T2Range = (10e-3, 1.0),
		Reg = "lcurve",

		# Global variables from sliders
		ETL = ETL_slider,
		α = flip_angle_slider,
		β = refcon_slider,
		T21 = T21_slider,
		T22 = T22_slider,
		MWF = MWF_slider,

		save_figs = true,
	)

	fig
end

# ╔═╡ 1b62140b-f949-4470-a150-902c4811eb9b
#=
let
    # animation settings
    nframes = 45
    framerate = 30

	# DECAES settings
	T1 = 1.0
	TE = 10e-3
	SNR = 60
	nT2 = 80
	T2Range = (10e-3, 1.0)
	Reg = "lcurve"

    T22_iterator = cos_range(150e-3, 300e-3; length = nframes, phase = pi/2)
    FA_iterator = cos_range(150, 175; length = nframes, phase = pi)
	tsize_iterator = [
		# # Periodic grow/shrink
		# cos_range(0.950, 1.050; length = nframes, phase = pi/2);

		# Heartbeat
		smooth_bump(1.0, 1.075; width = 0.95, length = 10);
		smooth_bump(1.0, 0.975; width = 0.95, length = 5);
		smooth_bump(1.0, 1.050; width = 0.95, length = 10);
		smooth_bump(1.0, 1.000; width = 0.95, length = 20);
	]
	anim_iterator = Iterators.zip(T22_iterator, FA_iterator, tsize_iterator)

    ETL = ETL_slider
    β = refcon_slider
    T21 = T21_slider
    MWF = MWF_slider

    fig, ax, ax_sub, txt, epg_line, epg_band, t2_bands, t2_lines = decaes_logo(;
        # Global variables
        T1, TE, SNR, nT2, T2Range, Reg,

        # Local variables
        ETL, β, T21, MWF,

        # Variables for gif
        T22 = T22_iterator[1],
		α = FA_iterator[1],
    )

    record(fig, "logo.gif", anim_iterator; framerate = framerate) do (T22, α, tsize)
		# Delete current subplots
		delete!(ax, txt)
        delete!(ax, epg_line)
        delete!(ax, epg_band)
        delete!.(ax_sub, t2_lines)
        delete!.(ax_sub, t2_bands)

        # Replace subplots with updated curves
		txt = plot_decaes_text!(
			ax; textsize = tsize .* textsize, xposition, yposition)
        echo_times, epg_curve = epg_decay_curve(
			ETL, α, TE, T21, T22, T1, β, MWF, SNR)
        epg_line, epg_band = plot_epg_curve!(ax, echo_times, epg_curve)

        t2maps, t2dist = t2_distribution(epg_curve, nT2, TE, T2Range, Reg)
        t2_lines, t2_bands = plot_t2_dist!(ax_sub, t2maps, t2dist)
    end

    PlutoUI.LocalResource("logo.gif")
end
=#

# ╔═╡ 801045fa-d596-485e-84f6-f5c74629cb2d
function cos_range(lo, hi; length, phase)
	θ = range(0, 2π; length) .- phase
	x = (lo+hi)/2 .+ (hi-lo)/2 .* cos.(θ)
	return x
end;

# ╔═╡ 2fc8ef73-b5c2-455a-a1a1-f85cf45fc4ca
function smooth_bump(lo, hi; width = 1, length)
	ϵ = sqrt(eps())
	u = range(-1, 1; length)
	x = @. lo + (hi - lo) * exp(-1 / (1 - clamp(u/width, -1+ϵ, 1-ϵ)^2))
	return x
end;

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
DECAES = "d84fb938-a666-558e-89d9-d531edc6724f"
Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
CairoMakie = "~0.6.2"
DECAES = "~0.4.0"
Makie = "~0.14.2"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[Animations]]
deps = ["Colors"]
git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d"
uuid = "27a7e980-b3e6-11e9-2bcd-0b925532e340"
version = "0.4.1"

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "045ff5e1bc8c6fb1ecb28694abba0a0d55b5f4f5"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.17"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.0"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "ffabdf5297c9038973a0a3724132aa269f38c448"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.0"

[[Blosc]]
deps = ["Blosc_jll"]
git-tree-sha1 = "84cf7d0f8fd46ca6f1b3e0305b4b4a37afe50fd6"
uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
version = "0.7.0"

[[Blosc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "e747dac84f39c62aff6956651ec359686490134e"
uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
version = "1.21.0+0"

[[BufferedStreams]]
deps = ["Compat", "Test"]
git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d"
version = "1.0.0"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[Cairo]]
deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b"
uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
version = "1.0.5"

[[CairoMakie]]
deps = ["Base64", "Cairo", "Colors", "FFTW", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "SHA", "StaticArrays"]
git-tree-sha1 = "68628add03c2c7c2235834902aad96876f07756a"
uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
version = "0.6.2"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+0"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "dcc25ff085cf548bc8befad5ce048391a7c07d40"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.10.11"

[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[ColorBrewer]]
deps = ["Colors", "JSON", "Test"]
git-tree-sha1 = "61c5334f33d91e570e1d0c3eb5465835242582c4"
uuid = "a2cac450-b92f-5266-8821-25eda20663c8"
version = "0.4.0"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
git-tree-sha1 = "ed268efe58512df8c7e224d2e170afd76dd6a417"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.13.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "42a9b08d3f2f951c9b283ea427d96ed9f1f30343"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.5"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DECAES]]
deps = ["ArgParse", "Dierckx", "DocStringExtensions", "LightXML", "LinearAlgebra", "Logging", "LoggingExtras", "MAT", "NIfTI", "NLopt", "Parameters", "PolynomialRoots", "ProgressMeter", "Random", "SIMD", "SpecialFunctions", "StaticArrays", "Statistics", "TimerOutputs", "UnsafeArrays"]
git-tree-sha1 = "dceb9f969bd37e56bf889c151fb2a128ec91418f"
uuid = "d84fb938-a666-558e-89d9-d531edc6724f"
version = "0.4.0"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Dierckx]]
deps = ["Dierckx_jll"]
git-tree-sha1 = "5fefbe52e9a6e55b8f87cb89352d469bd3a3a090"
uuid = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
version = "0.5.1"

[[Dierckx_jll]]
deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
git-tree-sha1 = "a580560f526f6fc6973e8bad2b036514a4e3b013"
uuid = "cd4c43a9-7502-52ba-aa6d-59fb2a88580b"
version = "0.0.1+0"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Distributions]]
deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "d5fa76fdab7bf54615262bd639e8c2e9bf295444"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.10"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.0"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

[[FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "f985af3b9f4e278b1d24434cbb546d6092fca661"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.3"

[[FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3676abafff7e4ff07bbd2c42b3d8201f31653dcc"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.9+8"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "25b9cc23ba3303de0ad2eac03f840de9104c9253"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.0"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[FreeType]]
deps = ["CEnum", "FreeType2_jll"]
git-tree-sha1 = "cabd77ab6a6fdff49bfd24af2ebe76e6e018a2b4"
uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
version = "4.0.0"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FreeTypeAbstraction]]
deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics", "StaticArrays"]
git-tree-sha1 = "d51e69f0a2f8a3842bca4183b700cf3d9acce626"
uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
version = "0.9.1"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "4136b8a5668341e58398bb472754bff4ba0456ff"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.3.12"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "47ce50b742921377301e15005c96e979574e130b"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.1+0"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "2c1cf4df419938ece72de17f368a021ee162762e"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.0"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[GridLayoutBase]]
deps = ["GeometryBasics", "InteractiveUtils", "Match", "Observables"]
git-tree-sha1 = "97b481deb387775da2ca751cc542610750bcdfc3"
uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
version = "0.5.6"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HDF5]]
deps = ["Blosc", "Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires"]
git-tree-sha1 = "1d18a48a037b14052ca462ea9d05dee3ac607d23"
uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
version = "0.15.5"

[[HDF5_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
version = "1.12.0+1"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.12"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "8a954fed8ac097d5be04921d595f741115c1b2ad"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+0"

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[ImageCore]]
deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "75f7fea2b3601b58f24ee83617b528e57160cbfd"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.9.1"

[[ImageIO]]
deps = ["FileIO", "Netpbm", "PNGFiles", "TiffImages", "UUIDs"]
git-tree-sha1 = "d067570b4d4870a942b19d9ceacaea4fb39b69a1"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.5.6"

[[IndirectArrays]]
git-tree-sha1 = "c2a145a145dc03a7620af1444e0264ef907bd44f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "0.5.1"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "1470c80592cf1f0a35566ee5e93c5f8221ebc33a"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.3"

[[IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[Isoband]]
deps = ["isoband_jll"]
git-tree-sha1 = "f9b6d97355599074dc867318950adaa6f9946137"
uuid = "f1662d9f-8043-43de-a69a-05efc1cc6ff4"
version = "0.1.1"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JSONSchema]]
deps = ["HTTP", "JSON", "ZipFile"]
git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.3"

[[KernelDensity]]
deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.6.3"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LightXML]]
deps = ["Libdl", "XML2_jll"]
git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
version = "0.9.0"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.5"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoggingExtras]]
deps = ["Dates", "Logging"]
git-tree-sha1 = "dfeda1c1130990428720de0024d4516b1902ce98"
uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
version = "0.4.7"

[[Lz4_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
version = "1.9.3+0"

[[MAT]]
deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"]
git-tree-sha1 = "5c62992f3d46b8dce69bdd234279bb5a369db7d5"
uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
version = "0.10.1"

[[MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+1"

[[Makie]]
deps = ["Animations", "Artifacts", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Distributions", "DocStringExtensions", "FFMPEG", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "IntervalSets", "Isoband", "KernelDensity", "LinearAlgebra", "MakieCore", "Markdown", "Match", "Observables", "Packing", "PlotUtils", "PolygonOps", "Printf", "Random", "Serialization", "Showoff", "SignedDistanceFields", "SparseArrays", "StaticArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "UnicodeFun"]
git-tree-sha1 = "82f9d9de6892e5470372720d23913d14a08991fd"
uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
version = "0.14.2"

[[MakieCore]]
deps = ["Observables"]
git-tree-sha1 = "7bcc8323fb37523a6a51ade2234eee27a11114c8"
uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
version = "0.1.3"

[[MappedArrays]]
deps = ["FixedPointNumbers"]
git-tree-sha1 = "b92bd220c95a8bbe89af28f11201fd080e0e3fe7"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.3.0"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Match]]
git-tree-sha1 = "5cf525d97caf86d29307150fcba763a64eaa9cbe"
uuid = "7eb4fadd-790c-5f42-8a69-bfa0b872bfbf"
version = "1.1.0"

[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "575644e3c05b258250bb599e57cf73bbf1062901"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.9.22"

[[MathProgBase]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c"
uuid = "fdba3010-5040-5b88-9595-932c9decdf73"
version = "0.7.8"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[NIfTI]]
deps = ["Base64", "CodecZlib", "MappedArrays", "Mmap", "TranscodingStreams"]
git-tree-sha1 = "ab2dbb0dacaed807bc2b8c6f40c8dfec67a225ef"
uuid = "a3a9e032-41b5-5fc4-967a-a6b7a19844d3"
version = "0.5.6"

[[NLopt]]
deps = ["MathOptInterface", "MathProgBase", "NLopt_jll"]
git-tree-sha1 = "e7b4c6b8e10b8d1238112a21f5cc001138989dc7"
uuid = "76087f3c-5699-56af-9a33-bf431cd00edd"
version = "0.6.2"

[[NLopt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2b597c46900f5f811bec31f0dcc88b45744a2a09"
uuid = "079eb43e-fd8e-5478-9966-2cf3e3edb778"
version = "2.7.0+0"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[Netpbm]]
deps = ["FileIO", "ImageCore"]
git-tree-sha1 = "18efc06f6ec36a8b801b23f076e3c6ac7c3bf153"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.2"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.2"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "520e28d4026d16dcf7b8c8140a3041f0e20a9ca8"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.7"

[[Packing]]
deps = ["GeometryBasics"]
git-tree-sha1 = "f4049d379326c2c7aa875c702ad19346ecb2b004"
uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566"
version = "0.4.1"

[[PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "0fa5e78929aebc3f6b56e1a88cf505bb00a354c4"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.8"

[[Pango_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9bc1871464b12ed19297fbc56c4fb4ba84988b0d"
uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
version = "1.47.0+0"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "501c20a63a34ac1d015d5304da0e645f42d91c9f"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.11"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[PolygonOps]]
git-tree-sha1 = "c031d2332c9a8e1c90eca239385815dc271abb22"
uuid = "647866c9-e3ac-4575-94e7-e3d426903924"
version = "0.1.1"

[[PolynomialRoots]]
git-tree-sha1 = "5f807b5345093487f733e520a1b7395ee9324825"
uuid = "3a141323-8675-5d76-9d11-e1df1406c778"
version = "1.0.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Ratios]]
git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.0"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SIMD]]
git-tree-sha1 = "9ba33637b24341aba594a2783a502760aa0bff04"
uuid = "fdea26ae-647d-5447-a871-4b548cad5224"
version = "3.3.1"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[SignedDistanceFields]]
deps = ["Random", "Statistics", "Test"]
git-tree-sha1 = "d263a08ec505853a5ff1c1ebde2070419e3f28e9"
uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96"
version = "0.4.0"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.5.1"

[[StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "2740ea27b66a41f9d213561a04573da5d3823d4b"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.2.5"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "a43a7b58a6e7dc933b2fa2e0ca653ccf8bb8fd0e"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.6"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.8"

[[StatsFuns]]
deps = ["LogExpFunctions", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "30cd8c360c54081f806b1ee14d2eecbef3c04c49"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.8"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "Tables"]
git-tree-sha1 = "44b3afd37b17422a62aea25f04c1f7e09ce6b07f"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.5.1"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[TiffImages]]
deps = ["ColorTypes", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "OrderedCollections", "PkgVersion", "ProgressMeter"]
git-tree-sha1 = "03fb246ac6e6b7cb7abac3b3302447d55b43270e"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.4.1"

[[TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "209a8326c4f955e2442c07b56029e88bb48299c7"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.12"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.5"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[UnsafeArrays]]
git-tree-sha1 = "be71ca7d62a9d61ecaa20e2d0e1920b791623626"
uuid = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
version = "1.0.2"

[[WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.3"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[isoband_jll]]
deps = ["Libdl", "Pkg"]
git-tree-sha1 = "a1ac99674715995a536bbce674b068ec1b7d893d"
uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4"
version = "0.2.2+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"
"""

# ╔═╡ Cell order:
# ╠═c2274306-49b5-45c0-ae4b-a3b8dc3a3835
# ╠═5cf6be3c-fe58-42a8-840e-8a21a5a81584
# ╟─9d4a6681-2b40-442f-8cc8-4ae444bc6215
# ╟─dfac902d-bb87-4fea-a898-d6792e9b793b
# ╠═e55b552a-00a1-44d3-9d0b-dda4ceac1f59
# ╠═7cd9afdf-6675-4416-aec7-9eba65eafbe2
# ╠═15b31eb3-8ddc-4fab-9bda-41cdbe8bc710
# ╟─446a6ee0-b25b-47ef-b4fc-af59f30d27fd
# ╠═118177e3-2b40-4cea-bba6-ab48ba492fbd
# ╟─499d8be9-ca07-40c9-9a22-79bb99ed58fa
# ╠═f217fde8-08b8-4f33-8149-69f900847bd3
# ╟─0e7420c8-bfd6-4c91-a7bd-6c205eb24452
# ╠═02d229c2-e2c0-4a36-bc6a-163778be85b9
# ╠═db1d1148-7bf7-4a1d-b568-52b55a8ec1f3
# ╠═1b62140b-f949-4470-a150-902c4811eb9b
# ╠═801045fa-d596-485e-84f6-f5c74629cb2d
# ╠═2fc8ef73-b5c2-455a-a1a1-f85cf45fc4ca
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002
