text:
"""
This submodule provides a collection of *token* names for use in lexer and theme
definitions. The table shown below lists all the exported tokens as well as the
abbreviations used when printing stylesheets and highlighted source code.

$(TABLE)
"""
module Tokens

import ..Highlights

const __TOKENS__ = [
    :TEXT
    :WHITESPACE
    :ESCAPE
    :ERROR
    :OTHER

    :KEYWORD
    :KEYWORD_CONSTANT
    :KEYWORD_DECLARATION
    :KEYWORD_NAMESPACE
    :KEYWORD_PSEUDO
    :KEYWORD_RESERVED
    :KEYWORD_TYPE

    :NAME
    :NAME_ATTRIBUTE
    :NAME_BUILTIN
    :NAME_BUILTIN_PSEUDO
    :NAME_CLASS
    :NAME_CONSTANT
    :NAME_DECORATOR
    :NAME_ENTITY
    :NAME_EXCEPTION
    :NAME_FUNCTION
    :NAME_FUNCTION_MAGIC
    :NAME_PROPERTY
    :NAME_LABEL
    :NAME_NAMESPACE
    :NAME_OTHER
    :NAME_TAG
    :NAME_VARIABLE
    :NAME_VARIABLE_CLASS
    :NAME_VARIABLE_GLOBAL
    :NAME_VARIABLE_INSTANCE
    :NAME_VARIABLE_MAGIC

    :LITERAL
    :LITERAL_DATE

    :STRING
    :STRING_AFFIX
    :STRING_BACKTICK
    :STRING_CHAR
    :STRING_DELIMITER
    :STRING_DOC
    :STRING_DOUBLE
    :STRING_ESCAPE
    :STRING_HEREDOC
    :STRING_INTERPOL
    :STRING_OTHER
    :STRING_REGEX
    :STRING_SINGLE
    :STRING_SYMBOL

    :NUMBER
    :NUMBER_BIN
    :NUMBER_FLOAT
    :NUMBER_HEX
    :NUMBER_INTEGER
    :NUMBER_INTEGER_LONG
    :NUMBER_OCT

    :OPERATOR
    :OPERATOR_WORD

    :PUNCTUATION

    :COMMENT
    :COMMENT_HASHBANG
    :COMMENT_MULTILINE
    :COMMENT_PREPROC
    :COMMENT_PREPROCFILE
    :COMMENT_SINGLE
    :COMMENT_SPECIAL

    :GENERIC
    :GENERIC_DELETED
    :GENERIC_EMPH
    :GENERIC_ERROR
    :GENERIC_HEADING
    :GENERIC_INSERTED
    :GENERIC_OUTPUT
    :GENERIC_PROMPT
    :GENERIC_STRONG
    :GENERIC_SUBHEADING
    :GENERIC_TRACEBACK
]

const __INDICES__ = Dict([(s, n) for (n, s) in enumerate(__TOKENS__)])
const __FALLBACKS__ = zeros(Int, length(__TOKENS__))

struct TokenValue
    value::Int
end
TokenValue(s::Symbol) = TokenValue(__INDICES__[s])

function parent(s::Symbol)
    p = join(split(string(s), '_')[1:end-1], '_')
    return isempty(p) ? :TEXT : Symbol(p)
end

Base.show(io::IO, t::TokenValue) = print(io, '<', __TOKENS__[t.value], '>')

for (n, each) in enumerate(__TOKENS__)
    __FALLBACKS__[n] = __INDICES__[parent(each)]
    @eval export $each
    @eval const $each = TokenValue($n)
end

# Shortnames.

const __SHORTNAMES__ = Symbol[]

let cache = Dict{Symbol, Int}()
    for (nth, each) in enumerate(__TOKENS__)
        short = Symbol(join(map(first, split(lowercase(string(each)), '_'))))
        if haskey(cache, short)
            cache[short] += 1
            push!(__SHORTNAMES__, Symbol(string(short, Char(cache[short] + Int('A') - 1))))
        else
            cache[short] = 1
            push!(__SHORTNAMES__, short)
        end
    end
    # Check that we've actually got a unique set of shortnames.
    @assert unique(__SHORTNAMES__) == __SHORTNAMES__
end

# Make a table to include in the module's docstring.

const TABLE =
    let buffer = IOBuffer(),
        col_1_label = "Token",
        col_2_label = "Abbreviation",
        col_1 = mapreduce(s -> length(string(s)), max, __TOKENS__, init =  length(col_1_label) ) + 2,
        col_2 = mapreduce(s -> length(string(s)), max, __SHORTNAMES__, init = length(col_2_label)) + 2
        println(buffer, "| ", rpad(col_1_label, col_1), " | ", rpad(col_2_label, col_2), " |")
        println(buffer, "|:", rpad("-", col_1, "-"),    " | ", rpad("-", col_2, "-"),    ":|")
        for (long, short) in zip(__TOKENS__, __SHORTNAMES__)
            println(buffer, "| ", rpad("`$long`", col_1), " | ", rpad("`$short`", col_2), " |")
        end
        Highlights.takebuf_str(buffer)
    end

end


text:
@testset ExtendedTestSet "hausdorff_metric" begin
    @testset ExtendedTestSet "hausdorff_metric" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y2 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y2; dims=3)

        edges = find_edges(y)
        @test hausdorff_metric(edges, edges) ≈ 0
    end

    @testset ExtendedTestSet "hausdorff_metric" begin
        y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y2 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = cat(y1, y2; dims=3)

        edges = find_edges(y)
        @test hausdorff_metric(edges, edges, 50) == (0.0, 0.0)
    end
end

@testset ExtendedTestSet "mean_hausdorff_2D" begin
    @testset ExtendedTestSet "mean_hausdorff_2D" begin
        x = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        @test mean_hausdorff_2D(x, y, Distances.Euclidean(), ImageEdgeDetection.Canny()) ≈ 0
    end
end

ext: jl


text:
using CUDAdrv, CUDAnative

if Sys.iswindows()
    function hello_world()
        @cuprintf("Greetings from block %lld, thread %lld!\n", Int64(blockIdx().x), Int64(threadIdx().x))
        return
    end
else
    function hello_world()
       @cuprintf("Greetings from block %ld, thread %ld!\n", Int64(blockIdx().x), Int64(threadIdx().x))
       return
   end
end
@cuda blocks=2 threads=2 hello_world()
synchronize()


text:
struct _DUMMY_STRUCT
end

# anymissing(::_DUMMY_STRUCT) = false
nunique(::_DUMMY_STRUCT) =  false
stdze!(::_DUMMY_STRUCT) = false
stdze(::_DUMMY_STRUCT) = false
select(::_DUMMY_STRUCT) = false

byrow(ds::AbstractDataset, ::typeof(Base.sum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_sum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.sum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, sum, [col]; by = by, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.sum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_sum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.sum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, sum, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(prod), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_prod(ds, by, cols; threads = threads)
byrow(ds::AbstractDataset, ::typeof(prod), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, prod, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(count), cols::MultiColumnIndex = :; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10) = row_count(ds, by, cols; threads =threads)
byrow(ds::AbstractDataset, ::typeof(count), col::ColumnIndex; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, count, [col], by = by, threads = threads)

# byrow(ds::AbstractDataset, ::typeof(anymissing), cols::MultiColumnIndex = names(ds, Union{Missing, Number})) = row_anymissing(ds, cols)

function expand_Base_Fix(f, f2)
	if f isa Base.Fix2
		return _bool(x->f.f(f2(x), f.x))
	elseif f isa Base.Fix1
		return _bool(x->f.f(f.x, f2(x)))
	else
		return x->f(f2(x))
	end
end

function byrow(ds::AbstractDataset, ::typeof(any), cols::MultiColumnIndex = :; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10, mapformats = false)
	colsidx = multiple_getindex(index(ds), cols)
	if by isa AbstractVector
		if mapformats
			by = map((x,y)->expand_Base_Fix(x, getformat(ds, y)), by, colsidx)
		end
	else
		if mapformats
			by = map(y->expand_Base_Fix(by, getformat(ds, y)), colsidx)
		end
	end
	row_any(ds, by, colsidx, threads = threads)

end

byrow(ds::AbstractDataset, ::typeof(any), col::ColumnIndex; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10, mapformats = false) = byrow(ds, any, [col]; by = by, threads = threads, mapformats = mapformats)

function byrow(ds::AbstractDataset, ::typeof(all), cols::MultiColumnIndex = :; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10, mapformats = false)
	colsidx =  multiple_getindex(index(ds), cols)
	if by isa AbstractVector
		if mapformats
			by = map((x,y)->expand_Base_Fix(x, getformat(ds, y)), by, colsidx)
		end
	else
		if mapformats
			by = map(y->expand_Base_Fix(by, getformat(ds, y)), colsidx)
		end
	end
	row_all(ds, by, colsidx, threads = threads)
end
byrow(ds::AbstractDataset, ::typeof(all), col::ColumnIndex; by = isequal(true), threads = nrow(ds) > Threads.nthreads()*10, mapformats = false) = byrow(ds, all, [col]; by = by, threads = threads, mapformats = mapformats)

byrow(ds::AbstractDataset, ::typeof(isequal), cols::MultiColumnIndex; with = nothing, threads = nrow(ds) > Threads.nthreads()*10) = row_isequal(ds, cols, by = with, threads = threads)
byrow(ds::AbstractDataset, ::typeof(isequal), cols::ColumnIndex; with = nothing, threads = nrow(ds) > Threads.nthreads()*10) = row_isequal(ds, cols, by = with, threads = threads)

byrow(ds::AbstractDataset, ::typeof(isless), cols::MultiColumnIndex; with, threads = nrow(ds) > Threads.nthreads()*10, rev::Bool = false, lt = isless) = row_isless(ds, cols, with, threads = threads, rev = rev, lt = lt)
byrow(ds::AbstractDataset, ::typeof(isless), col::ColumnIndex; with, threads = nrow(ds) > Threads.nthreads()*10, rev::Bool = false, lt = isless) = row_isless(ds, [col], with, threads = threads, rev = rev, lt = lt)

byrow(ds::AbstractDataset, ::typeof(in), cols::MultiColumnIndex; item, threads = nrow(ds) > Threads.nthreads()*10, eq = isequal) = row_in(ds, cols, item; threads = threads, eq = eq)

byrow(ds::AbstractDataset, ::typeof(findfirst), cols::MultiColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10, item = nothing, eq = isequal) = row_findfirst(ds, by, cols; threads = threads, item = item, eq = eq)
byrow(ds::AbstractDataset, ::typeof(findlast), cols::MultiColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10, item = nothing, eq = isequal) = row_findlast(ds, by, cols; threads = threads, item = item, eq = eq)

byrow(ds::AbstractDataset, ::typeof(select), cols::MultiColumnIndex; with, threads = nrow(ds) > Threads.nthreads()*10) = row_select(ds, cols, with, threads = threads)

byrow(ds::AbstractDataset, ::typeof(fill!), cols::MultiColumnIndex; with , by = ismissing, threads = nrow(ds) > Threads.nthreads()*10, rolling = false) = row_fill!(ds, cols, with, f = by, threads = threads, rolling = rolling)
byrow(ds::AbstractDataset, ::typeof(fill!), col::ColumnIndex; with , by = ismissing, threads = nrow(ds) > Threads.nthreads()*10, rolling = false) = byrow(ds, fill!, [col], with = with, by = by, threads = threads, rolling = rolling)
byrow(ds::AbstractDataset, ::typeof(fill), cols::MultiColumnIndex; with , by = ismissing, threads = nrow(ds) > Threads.nthreads()*10, rolling = false) = row_fill!(copy(ds), cols, with, f = by, threads = threads, rolling = rolling)
byrow(ds::AbstractDataset, ::typeof(fill), col::ColumnIndex; with , by = ismissing, threads = nrow(ds) > Threads.nthreads()*10, rolling = false) = byrow(copy(ds), fill!, [col], with = with, by = by, threads = threads, rolling = rolling)

byrow(ds::AbstractDataset, ::typeof(coalesce), cols::MultiColumnIndex; threads = nrow(ds) > Threads.nthreads()*10) = row_coalesce(ds, cols; threads = threads)

byrow(ds::AbstractDataset, ::typeof(mean), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_mean(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(mean), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, mean, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.maximum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_maximum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.maximum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, maximum, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.minimum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_minimum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.minimum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, minimum, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.argmin), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_argmin(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.argmin), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, argmin, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.argmax), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_argmax(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.argmax), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, argmax, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.maximum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_maximum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.maximum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, maximum, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.minimum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_minimum(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.minimum), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, minimum, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.argmin), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_argmin(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.argmin), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, argmin, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.argmax), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, threads = nrow(ds) > Threads.nthreads()*10) = row_argmax(ds, by, cols, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.argmax), col::ColumnIndex; by = identity, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, argmax, [col]; by = by, threads = threads)

byrow(ds::AbstractDataset, ::typeof(var), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, dof = true, threads = nrow(ds) > Threads.nthreads()*10) = row_var(ds, by, cols; dof = dof, threads = threads)
byrow(ds::AbstractDataset, ::typeof(var), col::ColumnIndex; by = identity, dof = true, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, var, [col]; by = by, dof = dof, threads = threads)

byrow(ds::AbstractDataset, ::typeof(std), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, dof = true, threads = nrow(ds) > Threads.nthreads()*10) = row_std(ds, by, cols; dof = dof, threads = threads)
byrow(ds::AbstractDataset, ::typeof(std), col::ColumnIndex; by = identity, dof = true, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, std, [col]; by = by, dof = dof, threads = threads)

byrow(ds::AbstractDataset, ::typeof(nunique), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); by = identity, count_missing = true) = row_nunique(ds, by, cols; count_missing = count_missing)
byrow(ds::AbstractDataset, ::typeof(nunique), col::ColumnIndex; by = identity, count_missing = true) = byrow(ds, nunique, [col]; by = by, count_missing = count_missing)

byrow(ds::AbstractDataset, ::typeof(Base.cumsum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumsum(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.cumsum), col::ColumnIndex; missings = :ignore, threads = nrow(ds)> Threads.nthreads()) = byrow(ds, cumsum, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.cumprod!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumprod!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.cumprod!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumprod!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.cumprod), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumprod(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.cumprod), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumprod, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(Base.cumsum!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumsum!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(Base.cumsum!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumsum!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.cumsum), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumsum(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.cumsum), col::ColumnIndex; missings = :ignore, threads = nrow(ds)> Threads.nthreads()) = byrow(ds, cumsum, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.cumprod!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumprod!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.cumprod!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumprod!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.cumprod), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumprod(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.cumprod), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumprod, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(IMD.cumsum!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cumsum!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(IMD.cumsum!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cumsum!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(cummin!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cummin!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(cummin!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cummin!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(cummin), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cummin(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(cummin), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cummin, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(cummax!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cummax!(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(cummax!), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cummax!, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(cummax), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = row_cummax(ds, cols, missings = missings, threads = threads)
byrow(ds::AbstractDataset, ::typeof(cummax), col::ColumnIndex; missings = :ignore, threads = nrow(ds)>Threads.nthreads()*10) = byrow(ds, cummax, [col], missings = missings, threads = threads)

byrow(ds::AbstractDataset, ::typeof(sort), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); threads = true, kwargs...) = threads ? hp_row_sort(ds, cols; kwargs...) : row_sort(ds, cols; kwargs...)
byrow(ds::AbstractDataset, ::typeof(sort), col::ColumnIndex; threads = true, kwargs...) = byrow(ds, sort, [col]; threads = threads, kwargs...)

byrow(ds::AbstractDataset, ::typeof(sort!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); threads = true, kwargs...) = threads ? hp_row_sort!(ds, cols; kwargs...) : row_sort!(ds, cols; kwargs...)
# byrow(ds::AbstractDataset, ::typeof(sort!), col::ColumnIndex; threads = true, kwargs...) = byrow(ds, sort!, [col]; threads = threads, kwargs...)

byrow(ds::AbstractDataset, ::typeof(issorted), cols::MultiColumnIndex; threads = nrow(ds) > Threads.nthreads()*10, rev = false, lt = isless) = row_issorted(ds, cols; rev = rev, lt = lt, threads = threads)

byrow(ds::AbstractDataset, ::typeof(stdze), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); threads = true) = row_stdze(ds, cols, threads = threads)

byrow(ds::AbstractDataset, ::typeof(stdze!), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); threads = true) = row_stdze!(ds, cols, threads = threads)

function byrow(ds::AbstractDataset, ::typeof(hash), cols::MultiColumnIndex = :; by = identity, mapformats = false, threads = nrow(ds) > Threads.nthreads()*10)
	colsidx = multiple_getindex(index(ds), cols)
	if mapformats
		by = map(y->expand_Base_Fix(by, getformat(ds, y)), colsidx)
	end
	row_hash(ds, by, cols, threads = threads)
end

byrow(ds::AbstractDataset, ::typeof(hash), col::ColumnIndex; by = identity, mapformats = false, threads = nrow(ds) > Threads.nthreads()*10) = byrow(ds, hash, [col]; by = by, mapformats = mapformats, threads = threads)

byrow(ds::AbstractDataset, ::typeof(join), col::MultiColumnIndex; threads = nrow(ds) > Threads.nthreads()*10, delim = "", last = "") = row_join(ds, col, threads = threads, delim = delim, last = last)

byrow(ds::AbstractDataset, ::typeof(mapreduce), cols::MultiColumnIndex = names(ds, Union{Missing, Number}); op = .+, f = identity,  init = _missings(mapreduce(eltype, promote_type, view(_columns(ds),index(ds)[cols])), nrow(ds)), kwargs...) = mapreduce(f, op, eachcol(ds[!, cols]), init = init; kwargs...)

function byrow(ds::AbstractDataset, f::Function, cols::MultiColumnIndex; threads = nrow(ds)>1000)
	colsidx = multiple_getindex(index(ds), cols)
	length(colsidx) == 1 && return byrow(ds, f, colsidx[1]; threads = threads)
	threads ?  hp_row_generic(ds, f, cols) : row_generic(ds, f, cols)
end
function byrow(ds::AbstractDataset, f::Function, col::ColumnIndex; threads = nrow(ds)>1000, allowmissing::Bool = true)
	if threads
		T = Core.Compiler.return_type(f, Tuple{nonmissingtype(eltype(ds[!, col]))})
		if allowmissing
			res = Vector{Union{Missing, T}}(undef, nrow(ds))
		else
			res = Vector{T}(undef, nrow(ds))
		end
		_hp_map_a_function!(res, f, _columns(ds)[index(ds)[col]])
	else
		T = Core.Compiler.return_type(f, Tuple{nonmissingtype(eltype(ds[!, col]))})
		if allowmissing
			res = Vector{Union{Missing, T}}(undef, nrow(ds))
		else
			res = Vector{T}(undef, nrow(ds))
		end
		map!(f, res, _columns(ds)[index(ds)[col]])
	end
	res
end

# specific path for converting Any to suitable type
byrow(ds::AbstractDataset, ::typeof(identity), col::ColumnIndex) = identity.(_columns(ds)[index(ds)[col]])

# special case for converting the type of a column conveniently
byrow(ds::AbstractDataset, f::Type, col::ColumnIndex) = convert(Vector{Union{Missing, f}}, _columns(ds)[index(ds)[col]])


text:
function Base.convert{T <: GLImage}(::Type{T}, ::Type{Primitive}, image::FieldTraits.ComposableLike)
    prim = if haskey(image, Primitive)
        x = image[Primitive]
        isa(x, VertexArray) && return x
        isa(x, GeometryPrimitive) || error("Primitive needs to be of type GeometryPrimitive. Found: $x")
        x
    else
        r = get(image, Ranges)
        mins = Vec(minimum.(r))
        maxs = Vec(maximum.(r))
        FRect(mins, maxs .- mins)
    end
    resolution = get(image, MeshResolution, (2, 2))
    VertexArray(uvmesh(prim, resolution), 1)
end

"""
Supplies the resolution of a mesh. This field will be used when converting a
GeometryPrimitive to a mesh
"""
MeshResolution
@composed type GLImage
    Primitive::VertexArray
    ImageData::Texture
    SpatialOrder
end


text:
# 3a
"""
  fixef_csv(m; out="fixef.csv", intercept="days_0")

output the fixed effects of a fitted model to a csv file
Ben Bolker, 2020-08-21 (slightly modified)
https://stackoverflow.com/questions/32585319/export-linear-mixed-effects-model-outputs-in-csv-using-julia-language

"""
fixef_csv = function(m;  out="fixef.csv"; intercept="days_0")
    ct = coeftable(m);
    ct.rownms[1] = intercept;
    ct.colnms[1] = "beta";
    fixef = DataFrame(ct.cols);
    rename!(fixef, ct.colnms, makeunique = true)
    fixef[!, :term] = ct.rownms;
    CSV.write(out, fixef);
end

# Example
fixef_csv(m; outfn="fixef.csv", intercept="days_0")

# 3b
"""
    cms_csv(m, outfn="cms.csv", intercept="(Intercept)")

output the conditional modes of a fitted model to a csv file
Reinhold Kliegl, 2020-08-28
"""

cms_csv = function(m; out="cms.csv", intercept="days_0")
    ct = DataFrame(ranef(m)[:1]) .+  m.β;
    rename!(ct, subj_order[:1], makeunique = true)
    cms = DataFrame([[names(ct)]; collect.(eachrow(ct))], [:column; Symbol.(axes(ct, 1))])
    rename!(cms, ["Subj", fixef.term[1], fixef.term[2]], makeunique = true)
    CSV.write(out, cms);
end

# Example
cms_csv(m; outfn="cms.csv", intercept="days_0")

# 2
"""
    outfun(m, outfn="output.csv")

output the fixed effects of a fitted model to a csv file
Ben Bolker, 2020-08-21
https://stackoverflow.com/questions/32585319/export-linear-mixed-effects-model-outputs-in-csv-using-julia-language
"""
outfun = function(m, outfn="output.csv")
    ct = coeftable(m)
    coef_df = DataFrame(ct.cols);
    rename!(coef_df, ct.colnms, makeunique = true)
    coef_df[!, :term] = ct.rownms;
    CSV.write(outfn, coef_df);
end

# 1
# Does NOT work in VS Code as expected; possibly ok in Terminals
function restart()
    startup = """
        Base.ACTIVE_PROJECT[]=$(repr(Base.ACTIVE_PROJECT[]))
        Base.HOME_PROJECT[]=$(repr(Base.HOME_PROJECT[]))
        cd($(repr(pwd())))
        """
    cmd = `$(Base.julia_cmd()) -ie $startup`
    atexit(()->run(cmd))
    exit(0)
end


text:
##
using Pkg
Pkg.activate(".")
Pkg.instantiate()
include("src/EasyStocks.jl")

## Setup
using CairoMakie # For plotting

mp = ModPar(q=0.0,ψ=0.0,xstar=4.0,η=0.)
np = NumPar(mp,nx=500,nsav=1,nq=2501,nα=101,nh=1,ygrd=[0.0,0.0],p1cons=false)
np = NumPar(mp,nx=500,nsav=1,nq=31,nα=51,nh=1,ygrd=[0.0,0.0],p1cons=false)
np.hgrd .= 0
MS1 = ModelSolution(mp,np)
@time SolveModel!(MS1)
mp = ModPar(q=mp.q,ψ=0.025,xstar=mp.xstar,η =mp.η)
MS2 = ModelSolution(mp,np)
@time SolveModel!(MS2)

# np.πrs[1] = 0.002
# np.πrs[:] = np.πrs/sum(np.πrs)
# MS3 = ModelSolution(mp,np)
# @time SolveModel!(MS3)
#
fig = Figure()
g2 = fig[1,1] = GridLayout()
g1 = fig[2,1] = GridLayout()

ax2  = Axis(g1[1,1],xlabel = "Wealth", ylabel = "α")
ax3  = Axis(g1[2,1],xlabel = "Wealth", ylabel = "s")
v2 = Axis(g2[1,1], xlabel ="Wealth", ylabel = "Utils")
v1 = Axis(g2[1,2], xlabel ="Wealth", ylabel = "Utils")
lines!(v1,np.xgrd,MS1.V[:,1,1],linestyle=:dash,color=:gray,label="Benchmark")
lines!(v1,np.xgrd,MS2.V[:,1,1],color=:orange,label="Housing")
# ylims!(v1,(MS2.V[np.nx÷7,1,1],0))
# ylims!(v2,(MS2.V[np.nx÷7,1,2],0))

lines!(v2,np.xgrd,MS1.V[:,1,2],linestyle=:dash,color=:gray,)
lines!(v2,np.xgrd,MS2.V[:,1,2],color=:orange,)

lines!(ax3,MS2.np.xgrd, MS1.s[:,1,1]+MS1.b[:,1,1], linestyle=:dash,color=:gray)
lines!(ax3,MS2.np.xgrd, MS2.s[:,1,1]+MS2.b[:,1,1], color=:orange)
lines!(ax3,MS2.np.xgrd, MS1.s[:,1,1], linestyle=:dash,color=:blue)
lines!(ax3,MS2.np.xgrd, MS2.s[:,1,1],color=:blue)
lines!(ax3,MS2.np.xgrd, MS1.b[:,1,1], linestyle=:dash,color=:green)
lines!(ax3,MS2.np.xgrd, MS2.b[:,1,1],color=:green)

lines!(ax2,MS2.np.xgrd, MS1.α[:,1,1], linestyle=:dash,color=:gray)
lines!(ax2,MS2.np.xgrd, MS2.α[:,1,1], color=:orange)
ylims!(ax2,(-0.05,1.05))
# text!(ax2, "↑\n Gambling\n for housing", position=(3.4,0.3),align=(:right,:center))
# text!(ax2, "Stay safe \nfor housing\n↓", position=(4.,0.3),align=(:left,:center))
# text!(ax2, "Gradually \nre-enter\n↓", position=(4.9,0.8),align=(:left,:center))

# lines!(ax3,MS2.np.xgrd, MS2.α[:,1],)
# lines!(ax3,MS2.np.xgrd, MS3.α[:,1],)
# ylims!(ax3,-0.05,1.05)
#leg = Legend(g1[1, 2], ax2)
axislegend(v1,position=:rb)
Label(g2[1, 1, Top()], "Second Period Value Function", valign = :bottom,)
Label(g2[1, 2, Top()], "First Period Value Function", valign = :bottom,)

fig
## Run-code
function Prob_xstar(MS)
    tmp = fill(0.0,MS.np.nx)
    for ix in eachindex(MS.np.xgrd)
        b = MS.b[ix,1,1]
        s = MS.s[ix,1,1]
        for irs in eachindex(np.πrs)
            tmp[ix] += np.πrs[irs]*(xp(b,s,mp.r,np.rsgrd[irs]) > MS.mp.xstar)
        end
    end
    return tmp
end

Prob_xstar(MS3)
scatter!(ax2,MS2.np.xgrd, Prob_xstar(MS1), linestyle=:dash,color=:gray)
scatter!(ax2,MS2.np.xgrd, Prob_xstar(MS2), color=:orange)
fig


text:
module TestCore

using UnsafeAtomics: UnsafeAtomics, acquire, release, acq_rel
using UnsafeAtomics.Internal: OP_RMW_TABLE, inttypes
using Test

function test_default_ordering()
    @testset for T in inttypes
        test_default_ordering(T)
    end
end

function test_default_ordering(T::Type)
    xs = T[rand(T), rand(T)]
    x1 = rand(T)
    x2 = rand(T)
    @debug "xs=$(repr(xs)) x1=$(repr(x1)) x2=$(repr(x2))"

    ptr = pointer(xs, 1)
    GC.@preserve xs begin
        @test UnsafeAtomics.load(ptr) === xs[1]
        UnsafeAtomics.store!(ptr, x1)
        @test xs[1] === x1
        desired = (old = x1, success = true)
        @test UnsafeAtomics.cas!(ptr, x1, x2) === (old = x1, success = true)
        @test xs[1] === x2
        @testset for (op, name) in OP_RMW_TABLE
            xs[1] = x1
            @test UnsafeAtomics.modify!(ptr, op, x2) === (x1, op(x1, x2))
            @test xs[1] === op(x1, x2)

            rmw = getfield(UnsafeAtomics, Symbol(name, :!))
            xs[1] = x1
            @test rmw(ptr, x2) === x1
            @test xs[1] === op(x1, x2)
        end
    end
end

function test_explicit_ordering(T::Type = UInt)
    xs = T[rand(T), rand(T)]
    x1 = rand(T)
    x2 = rand(T)
    @debug "xs=$(repr(xs)) x1=$(repr(x1)) x2=$(repr(x2))"

    ptr = pointer(xs, 1)
    GC.@preserve xs begin

        @test UnsafeAtomics.load(ptr, acquire) === xs[1]
        UnsafeAtomics.store!(ptr, x1, release)
        @test xs[1] === x1
        desired = (old = x1, success = true)
        @test UnsafeAtomics.cas!(ptr, x1, x2, acq_rel, acquire) === desired
        @test xs[1] === x2
        @testset for (op, name) in OP_RMW_TABLE
            xs[1] = x1
            @test UnsafeAtomics.modify!(ptr, op, x2, acq_rel) === (x1, op(x1, x2))
            @test xs[1] === op(x1, x2)

            rmw = getfield(UnsafeAtomics, Symbol(name, :!))
            xs[1] = x1
            @test rmw(ptr, x2, acquire) === x1
            @test xs[1] === op(x1, x2)
        end
    end
end

end  # module


text:
for N in [Float64, Rational{Int}, Float32]
    # random zonotope
    rand(Zonotope)

    # 1D Zonotope
    z = Zonotope(N[0], Matrix{N}(I, 1, 1))
    # Test Dimension
    @test dim(z) == 1
    # support function
    @test ρ(N[1], z) == ρ(N[-1], z) == N(1)
    # Test Support Vector
    d = N[1]
    @test σ(d, z) == N[1]
    d = N[-1]
    @test σ(d, z) == N[-1]

    # 2D Zonotope
    z = Zonotope(N[0, 0], Matrix{N}(I, 2, 2))
    # Test Dimension
    @test dim(z) == 2
    # support function
    @test ρ(N[1, 0], z) == ρ(N[-1, 0], z) == N(1)
    # Test Support Vector
    d = N[1, 0]
    @test σ(d, z) == N[1, 1] || N[1, -1]
    d = N[-1, 0]
    @test σ(d, z) == N[-1, 1] || N[-1, -1]
    d = N[0, 1]
    @test σ(d, z) == N[1, 1] || N[-1, 1]
    d = N[0, -1]
    @test σ(d, z) == N[1, -1] || N[-1, -1]

    # 2D Zonotope not 0-centered
    z = Zonotope(N[1, 2], Matrix{N}(I, 2, 2))
    # Test Dimension
    @test dim(z) == 2
    # Test Support Vector
    d = N[1, 0]
    @test σ(d, z) == N[2, 3]
    d = N[-1, 0]
    @test σ(d, z) == N[0, 3]
    d = N[0, 1]
    @test σ(d, z) == N[2, 3]
    d = N[0, -1]
    @test σ(d, z) == N[2, 1]

    # zero column in generators
    g = zeros(N, 2, 5)
    g[:, 3] = ones(N, 2)
    g[1, 2] = N(2)
    z = Zonotope(N[1, 2], g)
    @test size(z.generators) == (2, 2)

    # boundedness
    @test isbounded(z)

    # isempty
    @test !isempty(z)

    # an_element function
    @test an_element(z) ∈ z

    # concrete operations
    gens = N[1 1; -1 1]
    Z1 = Zonotope(N[1, 1], gens)
    Z2 = Zonotope(N[-1, 1], Matrix{N}(I, 2, 2))
    A = N[0.5 1; 1 0.5]

    # translation
    @test translate(Z1, N[1, 2]) == Zonotope(N[2, 3], gens)

    # concrete Minkowski sum
    Z3 = minkowski_sum(Z1, Z2)
    @test Z3.center == N[0, 2]
    @test Z3.generators == N[1 1 1 0; -1 1 0 1]

    # concrete linear map and scale
    Z4 = linear_map(A, Z3)
    @test Z4.center == N[2, 1]
    @test Z4.generators == N[-0.5 1.5 0.5 1; 0.5 1.5 1 0.5]
    Z5 = scale(0.5, Z3)
    @test Z5.center == N[0, 1]
    @test Z5.generators == N[0.5 0.5 0.5 0; -0.5 0.5 0 0.5]

    # intersection with a hyperplane
    H1 = Hyperplane(N[1, 1], N(3))
    intersection_empty, point = is_intersection_empty(Z1, H1, true)
    @test !is_intersection_empty(Z1, H1) && !intersection_empty &&
          point ∈ Z1 && point ∈ H1
    H2 = Hyperplane(N[1, 1], N(-11))
    @test is_intersection_empty(Z1, H2) && is_intersection_empty(Z1, H2, true)[1]
    @test !is_intersection_empty(H1, Z1)
    @test is_intersection_empty(H2, Z1) && is_intersection_empty(H2, Z1, true)[1]

    # isdisjoint
    result, w = isdisjoint(Z1, Z2, true)
    @test isdisjoint(Z1, Z2) && result && w == N[]
    Z3 = Zonotope(N[2, 1], Matrix{N}(I, 2, 2))
    @test_throws ErrorException isdisjoint(Z2, Z3, true)
    @test !isdisjoint(Z2, Z3)

    # test number of generators
    Z = Zonotope(N[2, 1], N[-0.5 1.5 0.5 1; 0.5 1.5 1 0.5])
    @test ngens(Z) == 4
    @test genmat(Z) == Z.generators
    @test ispermutation(collect(generators(Z)), [genmat(Z)[:, j] for j in 1:ngens(Z)])
    # test order reduction
    Zred1 = reduce_order(Z, 1)
    @test ngens(Zred1) == 2
    @test order(Zred1) == 1
    Zred2 = reduce_order(Z, 2)
    @test ngens(Zred2) == 4
    @test order(Zred2) == 2
    Z = Zonotope(N[2, 1], N[-0.5 1.5 0.5 1 0 1; 0.5 1.5 1 0.5 1 0])
    Zred3 = reduce_order(Z, 2)
    @test ngens(Zred3) == 4
    @test order(Zred3) == 2
    Znogen = Zonotope(N[1, 2], Matrix{N}(undef, 2, 0))
    @test ngens(Znogen) == 0
    @test genmat(Znogen) == Matrix{N}(undef, 2, 0)
    @test collect(generators(Znogen)) == Vector{N}()

    # test conversion from hyperrectangular sets
    Z = convert(Zonotope, Hyperrectangle(N[2, 3], N[4, 5]))
    @test Z.center == N[2, 3] && diag(Z.generators) == N[4, 5]
    convert(Zonotope, BallInf(N[5, 3], N(2)))
    # flat hyperrectangle
    Z = convert(Zonotope, Hyperrectangle(N[2, 3], N[0, 0]))
    @test Z.center == N[2, 3] && isempty(Z.generators)

    # convert the cartesian product of two hyperrectangles to a zonotope
    h1 = Hyperrectangle(N[1/2],  N[1/2])
    h2 = Hyperrectangle(N[2.5, 4.5],  N[1/2, 1/2])
    H = convert(Hyperrectangle, h1 × h2)
    Z = convert(Zonotope, h1 × h2)
    @test Z ⊆ H && H ⊆ Z

    # same for CartesianProductArray
    Z2 = convert(Zonotope, CartesianProductArray([h1, h2]))
    @test Z == Z2

    # split a zonotope
    Z = Zonotope(N[0, 0], N[1 1; -1 1])
    Z1, Z2 = split(Z, 1) # in this case the splitting is exact
    @test Z1 ⊆ Z && Z2 ⊆ Z

    # converts the cartesian product of two zonotopes to a new zonotope
    Z1 = Zonotope(N[0], hcat(N[1]))
    Z2 = Zonotope(N[1/2], hcat(N[1/2]))
    Z = convert(Zonotope, Z1×Z2)
    @test Z isa Zonotope && Z.center == N[0, 1/2] && Matrix(Z.generators) == N[1 0; 0 1/2]

    # conversion of the lazy linear map of an abstract hyperrectangle to a zonotope
    B = BallInf(N[0], N(1))
    M = hcat(N[1])
    Z = convert(Zonotope, M * B)
    @test Z isa Zonotope && Z.center == N[0] && Z.generators == hcat(N[1])

    # conversion of the lazy linear map of the cartesian product of hyperrectangular
    # sets to a zonotope
    B = BallInf(N[0], N(1))
    M = N[1 0; 0 -1]
    Z = convert(Zonotope, M * (B × B))
    @test Z isa Zonotope && Z.center == N[0, 0] && Z.generators == M

    # same for CPA
    Z2 = convert(Zonotope, M * CartesianProductArray([B, B]))
    @test Z2 == Z

    # list of constraints
    Z = Zonotope(zeros(N, 3), Matrix(N(1)*I, 3, 3))
    B = BallInf(zeros(N, 3), N(1)) # equivalent to Z
    if N != Rational{Int} || test_suite_polyhedra
        # the rational case uses vrep => needs Polyhedra
        constraints = constraints_list(Z)
        H = HPolytope(constraints)
        @test H ⊆ B && B ⊆ H
    end

    # =============================================
    # Overapproximate the convex hull of zonotopes
    # =============================================

    # same order
    Z1 = Zonotope(zeros(N, 2), hcat(N[1, 0]))
    Z2 = Zonotope(zeros(N, 2), hcat(N[0, 1]))
    Zch = overapproximate(ConvexHull(Z1, Z2), Zonotope)
    # the result is a diamond (tight)
    @test Zch == Zonotope(N[0, 0], N[0.5 0.5; 0.5 -0.5])

    # different order
    Z1 = Zonotope(zeros(N, 2), hcat(N[1, 0]))
    Z2 = Zonotope(N[0, 0], N[1 0; 0 1])
    # the result is a box (tight)
    Zch = overapproximate(ConvexHull(Z1, Z2), Zonotope)
    @test Zch == Zonotope(N[0, 0], N[1 0; 0 1])
end

for N in [Float64, Rational{Int}]
    # conversion to HPolytope
    # 1D
    Z = Zonotope(N[0], Matrix{N}(I, 1, 1))
    P = HPolytope(constraints_list(Z))
    for d in [N[1], N[-1]]
        @test ρ(d, P) == ρ(d, Z)
    end
    # 2D
    Z = Zonotope(N[0, 0], Matrix{N}(I, 2, 2))
    P = HPolytope(constraints_list(Z))
    for d in LazySets.Approximations.BoxDiagDirections{N}(2)
        @test ρ(d, P) == ρ(d, Z)
    end
    # sparse matrix (#1468)
    constraints_list(Zonotope(N[0, 0], sparse(N[1 0 ; 0 1])))
end

for N in [Float64]
    if test_suite_polyhedra
        # constraints_list for generator matrix with a zero row
        Z = Zonotope(N[0, 0], N[2 3; 0 0])
        P = tovrep(HPolygon(constraints_list(Z)))
        @test ispermutation(vertices_list(P), [N[5, 0], [-5, 0]])
    end
end


text:
using Plots
using SankeyPlots
using ReferenceTests
using Test

names = [
    "PV",
    "Electricity Buy",
    "Battery",
    "Heat pump",
    "Biomass",
    "Electricity Sell",
    "Electricity Demand",
    "Heat Demand",
]
src = [1, 1, 1, 1, 2, 2, 2, 3, 4, 5]
dst = [6, 3, 7, 4, 3, 7, 4, 7, 8, 8]
weights = [0.1, 0.3, 0.5, 0.5, 0.2, 2.8, 1, 0.45, 4.5, 3.3]
# labels = ["", "", "0.5", "0.5", "", "2.8", "1", "", "4.5", "3.3"]
energy_colors = palette(:seaborn_colorblind)[[9, 10, 3, 5, 2, 8, 1, 4]]

@testset "SankeyPlots.jl" begin
    @testset "readme" begin
        @test_reference "refs/readme.png" sankey(src, dst, weights)
        @test_reference "refs/readme_kwargs.png" sankey(
            src, dst, weights;
            node_labels=names,
            node_colors=energy_colors,
            edge_color=:gradient,
            label_position=:bottom,
            label_size=7,
            compact=true,
        )
    end

    @testset "edge_color" begin
        for c in (:src, :dst, "#789")
            @test_reference "refs/edge_color_$c.png" sankey(src, dst, weights; edge_color=c)
        end
    end

    @testset "label_position" begin
        for p in (:node, :legend, :top, :bottom, :left, :right)
            @test_reference "refs/label_position_$p.png" sankey(
                src, dst, weights; label_position=p
            )
        end
    end
end


text:
"""
    sample_regular(shape::Tuple{Int64,Int64},redFac::Float64;kargs...)

generates a regular sampling pattern for an Array of size `shape` with a subsampling factor `redFac`.

# Arguments
* `shape::NTuple{N,Int64}` - size of the Array to be sampled
* `redFac::Float64`        - subsampling factor
"""
function sample_regular(shape::Tuple, redFac::Float64; kargs...)
  sample_regular(prod(shape), redFac)
end

function sample_regular(maxInd::Int64, redFact::Float64; kargs...)
  outSet = range(1,step=floor(Int,redFact),length=floor(Int,maxInd/redFact))
  return outSet
end


text:
# Based off of code from julia SparseArrays.jl
struct GeneralSparseMatrixCSC{Tv, Ti<:Integer, VTv<:AbstractVector{Tv}, VTi<:AbstractVector{Ti}} <: SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}
    m::Int       # Number of rows
    n::Int       # Number of columns
    colptr::VTi  # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::VTi  # Row indices of stored values
    nzval::VTv   # Stored values, typically nonzeros

    function GeneralSparseMatrixCSC{Tv, Ti, VTv, VTi}(m::Integer, n::Integer,
                                                      colptr::AbstractVector{Ti},
                                                      rowval::AbstractVector{Ti},
                                                      nzval::AbstractVector{Tv}) where {Tv, Ti<:Integer,
                                                                                VTv<:AbstractVector{Tv},
                                                                                VTi<:AbstractVector{Ti}}
        new(Int(m), Int(n), colptr, rowval, nzval)
    end
end

function GeneralSparseMatrixCSC(m::Integer, n::Integer, colptr::AbstractVector,
                                rowval::AbstractVector, nzval::AbstractVector)
    Tv = eltype(nzval)
    VTv = typeof(nzval)
    Ti = promote_type(eltype(colptr), eltype(rowval))
    VTi = promote_type(typeof(colptr), typeof(rowval))

    GeneralSparseMatrixCSC{Tv, Ti, VTv, VTi}(m, n, colptr, rowval, nzval)
end

SparseArrays.getcolptr(S::GeneralSparseMatrixCSC) = getfield(S, :colptr)
SparseArrays.rowvals(S::GeneralSparseMatrixCSC) = getfield(S, :rowval)
SparseArrays.nonzeros(S::GeneralSparseMatrixCSC) = getfield(S, :nzval)
Base.size(S::GeneralSparseMatrixCSC) = (S.m, S.n)

function GeneralSparseMatrixCSC(S::SparseArrays.AbstractSparseMatrixCSC)
    m, n = size(S)
    colptr = SparseArrays.getcolptr(S)
    rowval = rowvals(S)
    nzval = nonzeros(S)

    return GeneralSparseMatrixCSC(m, n, colptr, rowval, nzval)
end

function Adapt.adapt_structure(to, S::GeneralSparseMatrixCSC)
    m, n = size(S)
    colptr = adapt(to, SparseArrays.getcolptr(S))
    rowval = adapt(to, rowvals(S))
    nzval = adapt(to, nonzeros(S))

    return GeneralSparseMatrixCSC(m, n, colptr, rowval, nzval)
end


text:
LCR6xxx_Vertical_Dict = Dict()
# Channel enabled
LCR6xxx_Vertical_Dict["on"]  = 1
LCR6xxx_Vertical_Dict["off"] = 0
# Channel coupling
LCR6xxx_Vertical_Dict["dc50"] = 0
LCR6xxx_Vertical_Dict["gnd"]  = 1
LCR6xxx_Vertical_Dict["dc"] = 2
LCR6xxx_Vertical_Dict["ac"] = 3
# Probe attenuation
LCR6xxx_Vertical_Dict["1x"] = 1
LCR6xxx_Vertical_Dict["10x"] = 10
LCR6xxx_Vertical_Dict["100x"] = 100
# BW limit
LCR6xxx_Vertical_Dict["full"]   = 0
LCR6xxx_Vertical_Dict["200MHz"] = 1
LCR6xxx_Vertical_Dict["20MHz"]  = 2
# EnhanceRes
LCR6xxx_Vertical_Dict["0b"]   = 0
LCR6xxx_Vertical_Dict["0.5b"] = 1
LCR6xxx_Vertical_Dict["1b"]   = 2
LCR6xxx_Vertical_Dict["1.5b"] = 3
LCR6xxx_Vertical_Dict["2b"]   = 4
LCR6xxx_Vertical_Dict["2.5b"] = 5
LCR6xxx_Vertical_Dict["3b"]   = 6

# Channel enabled
function set_state(obj::T;ch::String="1",st::String="on") where {T<:LCR6xxx}
	st = LCR6xxx_Vertical_Dict[st]
	write(obj,"VBS app.Acquisition.C$ch.View=$st")
end
# Channel input impedance
#function set_input_impedance(obj::T;ch::String="1",val::String="1M") where {T<:LCR6xxx}
#	st = LCR6xxx_Vertical_Dict[st]
#	write(obj,"VBS app.Acquisition.C$ch.View=$st")
#end
# Channel coupling
function set_coupling(obj::T;ch::String="1",cpl::String="dc") where {T<:LCR6xxx}
	cpl = LCR6xxx_Vertical_Dict[cpl]
	write(obj,"VBS app.Acquisition.C$ch.Coupling=$cpl")
end
# Probe attenuation
function set_atten(obj::T;ch::String="1",att::String="10x") where {T<:LCR6xxx}
	att = LCR6xxx_Vertical_Dict[att]
	write(obj,"VBS app.Acquisition.C$ch.ProbeAttenuation = $att")
end
# BW limit
function set_bw_lim(obj::T;ch::String="1",bw::String="full") where {T<:LCR6xxx}
	bw = LCR6xxx_Vertical_Dict[bw]
	write(obj,"VBS app.Acquisition.C$ch.BandwidthLimit=$bw")
end
# Nr of averages
function set_nr_avg(obj::T;ch::String="1",navg::Int=1) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.C$ch.AverageSweeps=$navg")
end
# EnhanceRes
function set_eres(obj::T;ch::String="1",bits::String="0b") where {T<:LCR6xxx}
	bits = LCR6xxx_Vertical_Dict[bits]
	write(obj,"VBS app.Acquisition.C$ch.EnhanceResType=$bits")
end
# Degauss
function degauss(obj::T;ch::String="1",probe::String="CP030A") where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.C$ch.InputB.$probe.Degauss")
end
# Vertical range
function set_vrange(obj::T;ch::String="1",vrang=8) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.C$ch.VerScale=$(vrang/8)")
end
# Vertical scale (V/div)
function set_vscale(obj::T;ch::String="1",volt=8) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.C$ch.VerScale=$volt")
end
# Vertical offset
function set_voffs(obj::T;ch::String="1",voffs=0) where {T<:LCR6xxx}
	write(obj,"VBS app.Acquisition.C$ch.VerOffset=$voffs")
end


text:
struct StanTable{T <: AbstractMatrix} <: Tables.AbstractColumns
    names::Vector{Symbol}
    lookup::Dict{Symbol, Int}
    matrix::T
end

Tables.istable(::Type{<:StanTable}) = true

# getter methods to avoid getproperty clash
import Base: names
names(m::StanTable) = getfield(m, :names)
matrix(m::StanTable) = getfield(m, :matrix)
lookup(m::StanTable) = getfield(m, :lookup)
# schema is column names and types
Tables.schema(m::StanTable{T}) where {T} = Tables.Schema(names(m), fill(eltype(T),
    size(matrix(m), 2)))

# column interface
Tables.columnaccess(::Type{<:StanTable}) = true
Tables.columns(m::StanTable) = m
# required Tables.AbstractColumns object methods
Tables.getcolumn(m::StanTable, ::Type{T}, col::Int, nm::Symbol) where {T} = matrix(m)[:, col]
Tables.getcolumn(m::StanTable, nm::Symbol) = matrix(m)[:, lookup(m)[nm]]
Tables.getcolumn(m::StanTable, i::Int) = matrix(m)[:, i]
Tables.columnnames(m::StanTable) = names(m)

Tables.isrowtable(::Type{StanTable}) = true

function convert_a3d(a3d_array, cnames, ::Val{:table};
    chains=1:size(a3d_array, 3),
    start=1,
    return_internals=false,
    kwargs...)

    cnames = String.(cnames)

    if !return_internals
        pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
        p = filter(p -> !(p in  pi), cnames)
    else
        p = cnames
    end

    lookup_dict = Dict{Symbol, Int}()
    for (idx, name) in enumerate(p)
        lookup_dict[Symbol(name)] = idx
    end

    mats = [a3d_array[start:end, :, i] for i in chains]
    mat = vcat(mats...)
    StanTable(Symbol.(p), lookup_dict, mat)
end

function convert_a3d(a3d_array, cnames, ::Val{:tables};
    chains=1:size(a3d_array, 3),
    start=1,
    return_internals=false,
    kwargs...)

    cnames = String.(cnames)

    if !return_internals
        pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
        p = filter(p -> !(p in  pi), cnames)
    else
        p = cnames
    end

    lookup_dict = Dict{Symbol, Int}()
    for (idx, name) in enumerate(p)
        lookup_dict[Symbol(name)] = idx
    end

    mats = [a3d_array[start:end, :, i] for i in chains]
    [StanTable(Symbol.(p), lookup_dict, mats[i])  for i in chains]
end

function matrix(st::StanTable, sym::Union{Symbol, String})
    n = string.(names(st))
    syms = string(sym)
    sel = String[]
    for (i, s) in enumerate(n)
        if length(s) > length(syms) && syms == n[i][1:length(syms)] &&
            n[i][length(syms)+1] in ['[', '.']
            append!(sel, [n[i]])
        end
    end
    length(sel) == 0 && error("$syms not in $n")
    tmp = st |> TableOperations.select(sel...) |> Tables.columntable
    Tables.matrix(tmp)
end

export
    StanTable,
    matrix,
    names

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
###########
# Library #
###########
# using MatrixLM
# using Distributions, Random, Statistics, LinearAlgebra, StatsBase
# using Random
using MatrixLMnet
using Helium
using Test

####################################################
# TEST Lasso vs Elastic Net (𝛼=1) - Simulated Data #
####################################################

#=
Description:
-----------

Model: 𝐘 = 𝐗 𝛃 𝐙' + 𝜎𝜖, with 𝜖∼𝑁(0,1)

Simulate data set consisting  of 20/20/200 observations and 8 predictors.
We let 𝛽₁ = (3, 1.5, 0, 0, 2, 0, 0, 0), 𝛽₂ = (0, 1.5, 0, 3.5, 2, 0, 0 , 2) where
𝛃 = [𝛽₁, 𝛽₂] and 𝜎 = 3.
The pairwise correlation between 𝑋ᵢ and 𝑋ⱼ was set to be 𝑐𝑜𝑟(𝑖,𝑗)=(0.5)^|𝑖−𝑗|.
Here, the Z matrix is an identity matrix.
=#

# Data testing directory name
dataDir = realpath(joinpath(@__DIR__,"data"))

# Get predictors
X = Helium.readhe(joinpath(dataDir, "Xmat.he"))

# Get response
Y = Helium.readhe(joinpath(dataDir, "Ymat.he"))

# Get Z matrix
Z = Helium.readhe(joinpath(dataDir, "Zmat.he"))

# Build raw data object from MatrixLM.jl
dat = RawData(Response(Y), Predictors(X, Z));

# Hyper parameters
λ = [10.0, 5.0, 3.0]
α = [1.0]

rng = MatrixLMnet.Random.MersenneTwister(2021)

############################################
# TEST 1 Lasso vs Elastic Net (𝛼=1) - ista #
############################################

# Elastic net penalized regression
est_ista_1 = MatrixLMnet.mlmnet(dat, λ, α, method = "ista", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_ista_1 = est_ista_1.B[:, :, 3, 1];

# Elastic net penalized regression
est_ista_2 = MatrixLMnet.mlmnet(dat, λ, method = "ista",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_ista_2 = est_ista_2.B[:, :, 3, 1];

# Lasso penalized regression - ista
B_ista= Helium.readhe(joinpath(dataDir, "B_ista.he"))

println("Lasso vs Elastic Net when α=1 test 1 - ista: ", @test (B_Net_ista_1 ≈ B_ista) && (B_Net_ista_2 ≈ B_ista))

#############################################
# TEST 2 Lasso vs Elastic Net (𝛼=1) - fista #
#############################################

# Elastic net penalized regression
est_fista_1 = MatrixLMnet.mlmnet(dat, λ, α, method = "fista", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_fista_1 = est_fista_1.B[:, :, 3, 1];

# Elastic net penalized regression
est_fista_2 = MatrixLMnet.mlmnet(dat, λ, method = "fista",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_fista_2 = est_fista_2.B[:, :, 3, 1];

# Lasso penalized regression - fista
B_fista= Helium.readhe(joinpath(dataDir, "B_fista.he"))

println("Lasso vs Elastic Net when α=1 test 2 - fista: ", @test (B_Net_fista_1 ≈ B_fista) && (B_Net_fista_2 ≈ B_fista))

##########################################################
# TEST 3 Lasso vs Elastic Net (𝛼=1) - fista backtracking #
##########################################################

# Elastic net penalized regression
est_fistabt_1 = MatrixLMnet.mlmnet(dat, λ, α, method = "fista_bt", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_fistabt_1 = est_fistabt_1.B[:, :, 3, 1];

# Elastic net penalized regression
est_fistabt_2 = MatrixLMnet.mlmnet(dat, λ, method = "fista_bt",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_fistabt_2 = est_fistabt_2.B[:, :, 3, 1];

# Lasso penalized regression - fista-bt
B_fistabt = Helium.readhe(joinpath(dataDir, "B_fistabt.he"))

println("Lasso vs Elastic Net when α=1 test 3 - fista-bt: ", @test (B_Net_fistabt_1 ≈ B_fistabt) && (B_Net_fistabt_2 ≈ B_fistabt))

############################################
# TEST 4 Lasso vs Elastic Net (𝛼=1) - admm #
############################################

# Elastic net penalized regression
est_admm_1 = MatrixLMnet.mlmnet(dat, λ, α, method = "admm", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_admm_1 = est_admm_1.B[:, :, 3, 1];

# Elastic net penalized regression
est_admm_2 = MatrixLMnet.mlmnet(dat, λ, method = "admm",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_admm_2 = est_admm_2.B[:, :, 3, 1];

# Lasso penalized regression - admm
B_admm = Helium.readhe(joinpath(dataDir, "B_admm.he"))

println("Lasso vs Elastic Net when α=1 test 4 - admm: ", @test (B_Net_admm_1 ≈ B_admm) && (B_Net_admm_2 ≈ B_admm))

##########################################
# TEST 5 Lasso vs Elastic Net (𝛼=1) - cd #
##########################################

# Elastic net penalized regression
MatrixLMnet.Random.seed!(rng)
est_cd_1 = MatrixLMnet.mlmnet(dat, λ, α, method = "cd", hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_cd_1 = est_cd_1.B[:, :, 3, 1];

# Elastic net penalized regression
MatrixLMnet.Random.seed!(rng)
est_cd_2 = MatrixLMnet.mlmnet(dat, λ, method = "cd",  hasZIntercept = false, hasXIntercept = false, isVerbose = false);
B_Net_cd_2 = est_cd_2.B[:, :, 3, 1];

# Lasso penalized regression - cd
B_cd = Helium.readhe(joinpath(dataDir, "B_cd.he"))

println("Lasso vs Elastic Net when α=1 test 5 - cd: ", @test ≈(B_Net_cd_1,  B_cd; atol=1.2e-4) && ≈(B_Net_cd_2, B_cd;  atol=1.2e-4))

println("Tests mlmnet finished!")


text:
function prepvariablelist(dataobject::InfoType, datatype::Symbol, vars::Array{Symbol,1}, lmax::Real, verbose::Bool)

    if datatype == :hydro
        nvarh = dataobject.nvarh

        # manage user selected variables
        #-------------------------------------
        nvarh_list=Int[]
        ivar=1
        read_cpu = false # do not load cpu number by default
        hydrovar_buffer = copy(vars)
        used_descriptors = Dict()

        if hydrovar_buffer == [:all]
            nvarh_list=[1,2,3,4,5]

            # read_cpu = true
            if in(:cpu, hydrovar_buffer) || in(:varn1, hydrovar_buffer)
             read_cpu = true
            end

            if nvarh > 5
                for ivar=6:nvarh
                     append!(nvarh_list, [ivar])
                end
            end
            #if dataobject.descriptor.usehydro == true
            #    for (i,idvar) in enumerate(dataobject.descriptor.hydro)
            #        used_descriptors[i] = idvar
            #    end
            #end

        else
            if in(:cpu, hydrovar_buffer) || in(:varn1, hydrovar_buffer)
             read_cpu = true
             filter!(e->e≠:cpu, hydrovar_buffer)
             filter!(e->e≠:varn1, hydrovar_buffer)
            end
            if in(:rho, hydrovar_buffer) || in(:var1, hydrovar_buffer)
             append!(nvarh_list, 1)
             filter!(e->e≠:rho, hydrovar_buffer)
             filter!(e->e≠:var1, hydrovar_buffer)
            end
            if in(:vx, hydrovar_buffer) || in(:var2, hydrovar_buffer)
             append!(nvarh_list, 2)
             filter!(e->e≠:vx, hydrovar_buffer)
             filter!(e->e≠:var2, hydrovar_buffer)
            end
            if in(:vy, hydrovar_buffer) || in(:var3, hydrovar_buffer)
             append!(nvarh_list, 3)
             filter!(e->e≠:vy, hydrovar_buffer)
             filter!(e->e≠:var3, hydrovar_buffer)
            end
            if in(:vz, hydrovar_buffer) || in(:var4, hydrovar_buffer)
             append!(nvarh_list, 4)
             filter!(e->e≠:vz, hydrovar_buffer)
             filter!(e->e≠:var4, hydrovar_buffer)
            end
            if in(:p, hydrovar_buffer) || in(:var5, hydrovar_buffer)
             append!(nvarh_list, 5)
             filter!(e->e≠:p, hydrovar_buffer)
             filter!(e->e≠:var5, hydrovar_buffer)
            end

            # if dataobject.descriptor.hydro != dataobject.variable_list
            #     for (ivar,idvar) in enumerate(dataobject.descriptor.hydro)
            #         if in(idvar, vars)
            #             append!(nvarh_list, ivar)
            #             used_descriptors[ivar] = idvar
            #         end
            #     end
            # end

            if length(hydrovar_buffer)>0
                for x in hydrovar_buffer
                    if occursin("var")
                        append!( nvarh_list, parse(Int, string(x)[4:end]) )
                    end
                end
             #append!( nvarh_list, map(x->parse(Int, string(x)[4:end]), hydrovar_buffer) ) #for Symbols
             #append!( nvarh_list, map(x->parse(Int,x), hydrovar_buffer) ) #for Strings
            end
        end

        if length(nvarh_list) == 0
            error("[Mera]: Simulation vars array is empty!")
        end

        #clean for double assignment
        nvarh_list = unique(nvarh_list)

        if maximum(nvarh_list) > maximum(nvarh)
            error("[Mera]: Simulation maximum variable=$(maximum(nvarh)) < your maximum variable=$(maximum(nvarh_list))")
        end

        # create vector to use selected hydrovars
        nvarh = dataobject.nvarh
        nvarh_corr = zeros(Int, nvarh )
        nvarh_i_list=[]
        for i =1:nvarh
            if in(i,nvarh_list)
                nvarh_corr[i] = findall(x -> x == i, nvarh_list)[1]
                append!(nvarh_i_list, i)
            end
        end

        nvarh_list_strings= Symbol[]
        if read_cpu append!(nvarh_list_strings, [:cpu]) end
        for i in nvarh_list
            #if !haskey(used_descriptors, i)
                if i < 6
                    append!(nvarh_list_strings, [Symbol("$(indices_tovariables[i])")])
                elseif i > 5
                    append!(nvarh_list_strings, [Symbol("var$i")])
                end
            #else
            #    append!(nvarh_list_strings, [used_descriptors[i]])
            #end
        end

        #println("nvarh_list",nvarh_list)
        #println("nvarh_corr",nvarh_corr)

        if verbose
            if lmax != dataobject.levelmin # if AMR
                println("Key vars=(:level, :cx, :cy, :cz)")
            else # if uniform grid
                println("Key vars=(:cx, :cy, :cz)")
            end

            if read_cpu
                println("Using var(s)=$(tuple(-1, nvarh_list...)) = $(tuple(nvarh_list_strings...)) ")
            else
                println("Using var(s)=$(tuple(nvarh_list...)) = $(tuple(nvarh_list_strings...)) ")
            end
            println()
        end

        return nvarh_list, nvarh_i_list, nvarh_corr, read_cpu, used_descriptors

    elseif datatype == :particles
        nvarp = dataobject.nvarp # vx, vy, vz, mass, birth

        #:level, :x, :y, :z, :id, :cpu, :vx, :vy, :vz, :mass, :birth]
        # manage user selected variables
        #-------------------------------------
        nvarp_list=Int[]
        ivar=1
        read_cpu = false
        particlesvar_buffer = copy(vars)
        used_descriptors = Dict()
        if particlesvar_buffer == [:all]
            for invarp = 1:nvarp
                append!(nvarp_list, invarp)
            end

            # read_cpu = true
            if in(:cpu, particlesvar_buffer) || in(:varn1, particlesvar_buffer)
             read_cpu = true
            end

            if nvarp > 5 && dataobject.descriptor.pversion <= 0
                for ivar=6:nvarp
                     append!(nvarp_list, [ivar])
                end
            elseif nvarp > 8 && dataobject.descriptor.pversion > 0
                for ivar=9:nvarp
                     append!(nvarp_list, [ivar])
                end
            end
            # if dataobject.use_particles_descriptor == true
            #     for (i,idvar) in enumerate(dataobject.particles_descriptor)
            #         used_descriptors[i] = idvar
            #     end
            # end
        else
            if in(:cpu, particlesvar_buffer) || in(:varn1, particlesvar_buffer)
             read_cpu = true
             filter!(e->e≠:cpu, particlesvar_buffer)
             filter!(e->e≠:varn1, particlesvar_buffer)
            end
            if in(:vx, particlesvar_buffer) || in(:var1, particlesvar_buffer)
             append!(nvarp_list, 1)
             filter!(e->e≠:vx, particlesvar_buffer)
             filter!(e->e≠:var1, particlesvar_buffer)
            end
            if in(:vy, particlesvar_buffer) || in(:var2, particlesvar_buffer)
             append!(nvarp_list, 2)
             filter!(e->e≠:vy, particlesvar_buffer)
             filter!(e->e≠:var2, particlesvar_buffer)
            end
            if in(:vz, particlesvar_buffer) || in(:var3, particlesvar_buffer)
             append!(nvarp_list, 3)
             filter!(e->e≠:vz, particlesvar_buffer)
             filter!(e->e≠:var3, particlesvar_buffer)
            end
            if in(:mass, particlesvar_buffer) || in(:var4, particlesvar_buffer)
             append!(nvarp_list, 4)
             filter!(e->e≠:mass, particlesvar_buffer)
             filter!(e->e≠:var4, particlesvar_buffer)
            end
            if in(:birth, particlesvar_buffer) || in(:var5, particlesvar_buffer) || in(:var7, particlesvar_buffer)

                if dataobject.descriptor.pversion <= 0 || in(:var5, particlesvar_buffer)
                    append!(nvarp_list, 5)
                    filter!(e->e≠:var5, particlesvar_buffer)
                elseif dataobject.descriptor.pversion > 0 || in(:var7, particlesvar_buffer)
                    append!(nvarp_list, 7)
                    filter!(e->e≠:var7, particlesvar_buffer)
                end
                filter!(e->e≠:birth, particlesvar_buffer)
            end

            if in(:metals, particlesvar_buffer) || in(:var8, particlesvar_buffer)
                if dataobject.descriptor.pversion > 0
                    append!(nvarp_list, 8)
                    filter!(e->e≠:var8, particlesvar_buffer)
                end
                filter!(e->e≠:metals, particlesvar_buffer)
            end

            # if dataobject.particles_descriptor != dataobject.particles_variable_list
            #     for (ivar,idvar) in enumerate(dataobject.particles_descriptor)
            #         if in(idvar, vars)
            #             append!(nvarp_list, ivar)
            #             used_descriptors[ivar] = idvar
            #             filter!(e->e≠idvar, particlesvar_buffer)
            #         end
            #     end
            # end

            if length(particlesvar_buffer)>0
                for x in particlesvar_buffer
                    if occursin("var")
                        append!( nvarh_list, parse(Int, string(x)[4:end]) )
                    end
                end

             #append!( nvarp_list, map(x->parse(Int, string(x)[4:end]), particlesvar_buffer) ) #for Symbols
            end
        end

        if length(nvarp_list) == 0
            error("[Mera]: Simulation vars array is empty!")
        end

        #clean for double assignment
        nvarp_list = unique(nvarp_list)

        if maximum(nvarp_list) > maximum(nvarp)
            error("[Mera]: Simulation maximum variable=$(maximum(nvarp)) < your maximum variable=$(maximum(nvarp_list))")
        end

        # create vector to use selected particlevars
        nvarp = dataobject.nvarp
        nvarp_corr = zeros(Int, nvarp )
        nvarp_i_list=[]
        for i =1:nvarp
            if in(i,nvarp_list)
                nvarp_corr[i] = findall(x -> x == i, nvarp_list)[1]
                append!(nvarp_i_list, i)
            end
        end

        nvarp_list_strings= Symbol[]
        if read_cpu append!(nvarp_list_strings, [:cpu]) end
        for i in nvarp_list
            #if !haskey(used_descriptors, i)
            if dataobject.descriptor.pversion == 0
                if i < 6
                    append!(nvarp_list_strings, [Symbol("$(indices_toparticlevariables[i])")])
                elseif i > 5
                    append!(nvarp_list_strings, [Symbol("var$i")])
                end
            elseif dataobject.descriptor.pversion > 0
                if i < 9 && i != 5 && i != 6
                    append!(nvarp_list_strings, [Symbol("$(indices_toparticlevariables_v1[i])")])
                elseif i > 8
                    append!(nvarp_list_strings, [Symbol("var$i")])
                end
            end
            #else
            #    append!(nvarp_list_strings, [used_descriptors[i]])
            #end
        end

        if verbose
                if dataobject.descriptor.pversion == 0
                    println("Key vars=(:level, :x, :y, :z, :id)")
                elseif dataobject.descriptor.pversion > 0
                    println("Key vars=(:level, :x, :y, :z, :id, :family, :tag)")
                end

            if read_cpu
                println("Using var(s)=$(tuple(-1, nvarp_list...)) = $(tuple(nvarp_list_strings...)) ")
            else
                nvarp_i_list_buffer = nvarp_i_list
                if dataobject.descriptor.pversion > 0
                    filter!(x->x≠6,nvarp_i_list_buffer)
                    filter!(x->x≠5,nvarp_i_list_buffer)
                end
                println("Using var(s)=$(tuple(nvarp_i_list_buffer...)) = $(tuple(nvarp_list_strings...)) ")
            end
            println()
        end
        return nvarp_list, nvarp_i_list, nvarp_corr, read_cpu, used_descriptors

    end

end

# index to variable assignment
global indices_tovariables = SortedDict( -1 =>  "cpu",
                                          0 => "level",

                                          1 => "rho",
                                          2 => "vx",
                                          3 => "vy",
                                          4 => "vz",
                                          5 => "p")

global indices_toparticlevariables = SortedDict(
                                     -1 =>  "cpu",
                                      0 => "level",

                                      1 => "vx",
                                      2 => "vy",
                                      3 => "vz",
                                      4 => "mass",
                                      5 => "birth")

global indices_toparticlevariables_v1 = SortedDict(
                                   -1 =>  "cpu",
                                    0 => "level",

                                    1 => "vx",
                                    2 => "vy",
                                    3 => "vz",
                                    4 => "mass",
                                    7 => "birth",
                                    8 => "metallicity")


text:
using Test
using PDELib
import Pluto, Pkg

function runtest(t)
    include(joinpath(@__DIR__,"..","examples",t*".jl"))
    eval(Meta.parse("$(t)()"))
end

@testset "VoronoiFVM" begin
    @test runtest("fv_laplace_rect")
    @test runtest("fv_laplace_circle")
end

@testset "GradientRobustMultiPhysics" begin
    @info "GradientRobustMultiPhysics tests disabled in the moment"
#    @test runtest("fe_liddrivencavity_autonewton")
end

function testnotebook(name)
    input=joinpath(@__DIR__,"..","examples",name*".jl")
    notebook=Pluto.load_notebook_nobackup(input)
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session,input; run_async=false)
    errored=false
    for c in notebook.cells
        if c.errored
            errored=true
            @error "Error in  $(c.cell_id): $(c.output.body[:msg])\n $(c.code)"
        end
    end
    !errored
end

notebooks=["Pluto-GridsAndVisualization",
           "Pluto-MultECatGrids"]

@testset "Notebooks" begin
     for notebook in notebooks
         @test testnotebook(notebook)
     end
end


text:
using DiffEqBase, Catalyst, DiffEqJump, DiffEqProblemLibrary.JumpProblemLibrary, Plots, Statistics, DataFrames
gr()
fmt = :png
JumpProblemLibrary.importjumpproblems()

N = 256
h = 1 / N
u0 = 10*ones(Int64, N)
tf = .01
methods = (Direct(),DirectFW(),SortingDirect(),NRM(),DirectCR(),RSSA())
shortlabels = [string(leg)[12:end-2] for leg in methods]
jprob   = prob_jump_diffnetwork
prob    = DiscreteProblem(u0, (0.0, tf), [1 / (h*h)])
rn      = jprob.network(N)
ploth   = plot(reuse=false)
for (i,method) in enumerate(methods)
    println("Benchmarking method: ", method)
    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
    sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
    plot!(ploth,sol.t,sol[Int(N//2),:],label=shortlabels[i], format=fmt)
end
plot!(ploth, title="Population at middle lattice site", xlabel="time",format=fmt)

function run_benchmark!(t, jump_prob, stepper)
    sol = solve(jump_prob, stepper)
    @inbounds for i in 1:length(t)
        t[i] = @elapsed (sol = solve(jump_prob, stepper))
    end
end

nsims = 50
benchmarks = Vector{Vector{Float64}}()
for method in methods
    jump_prob = JumpProblem(rn, prob, method, save_positions=(false, false))
    stepper = SSAStepper()
    t = Vector{Float64}(undef,nsims)
    run_benchmark!(t, jump_prob, stepper)
    push!(benchmarks, t)
end

medtimes = Vector{Float64}(undef,length(methods))
stdtimes = Vector{Float64}(undef,length(methods))
avgtimes = Vector{Float64}(undef,length(methods))
for i in 1:length(methods)
    medtimes[i] = median(benchmarks[i])
    avgtimes[i] = mean(benchmarks[i])
    stdtimes[i] = std(benchmarks[i])
end

df = DataFrame(names=shortlabels,medtimes=medtimes,relmedtimes=(medtimes/medtimes[1]),
                avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)

sa = [string(round(mt,digits=4),"s") for mt in df.medtimes]
bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
ylabel!("median relative to Direct")
title!("256 Site 1D Diffusion CTRW")

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLBenchmarksOutput

max_issues_repo_name: SciML/SciMLBenchmarksOutput

max_forks_repo_name: SciML/SciMLBenchmarksOutput


text:
# generate higher order functions of the propensity to be expanded
function gen_fs(network::DiffEqBase.AbstractReactionNetwork, momentstruct::MomentStruct, k_orders::Vector{Int}, e_orders::Vector{Int}; algebra::Algebra=sympy())
	propensities = mass_rate_MA.(network.reactions, algebra=algebra)
	propensities .* α_coeff(momentstruct.species, k_orders, e_orders)
end

# expand higher order functions of the propensity
function expand_fs(momentstruct::MomentStruct, fs::Vector{T}; algebra::Algebra=sympy()) where T <: AlgebraSet
	e_of_fs = [taylor_expand(f, momentstruct, algebra=algebra) for f in fs]
end

# differentiate powers of first order raw moments w.r.t. time
function dα_dt(momentstruct::MomentStruct, dμ::Vector{T}, n_orders::Vector{Int}, k_orders::Vector{Int}; algebra::Algebra=sympy()) where T <: AlgebraSet
	dα = zero(algebra)
	for i in eachindex(n_orders)
		dα += (n_orders[i] - k_orders[i]) * (1/momentstruct.species[i]) * α_coeff(momentstruct.species, n_orders, k_orders) * dμ[i]
	end
	dα
end

# differentiate mixed raw moments w.r.t time and express in terms of central moments
function dβ_dt(network::DiffEqBase.AbstractReactionNetwork, momentstruct::MomentStruct, k_order::Vector{Int}; algebra::Algebra=sympy())
	eq = repeat([zero(algebra)], length(network.reactions))
	for e_order in Iterators.product([0:k_max for k_max in k_order]...)
		if sum(e_order) == 0
			continue
		end
		fs = gen_fs(network, momentstruct, k_order, collect(e_order), algebra=algebra)
		e_of_fs = expand_fs(momentstruct, fs, algebra=algebra)
		coeffs = [binomial_coeff(k_order, collect(e_order)) * stoich_coeff(momentstruct.species, reaction, collect(e_order)) for reaction in network.reactions]
		eq += coeffs .* e_of_fs
	end
	sum(eq)
end

# central moments. For details of calculations, see https://arxiv.org/abs/1303.5848
function dcen_dt(network::DiffEqBase.AbstractReactionNetwork, momentstruct::MomentStruct, dμ::Vector{T}; algebra::Algebra=sympy()) where T <: AlgebraSet
	dcen = repeat([zero(algebra)], length(momentstruct.cen_moments) - momentstruct.num_species - 1)
	counter = 0
	for n_order in keys(momentstruct.cen_moments)
		if sum(n_order) <= 1
			continue
		else
			counter += 1
			for k_order in Iterators.product([0:n_max for n_max in n_order]...)
				__k = collect(k_order)
				β = momentstruct.raw_moments[__k]
				α = α_coeff(momentstruct.species, n_order, __k)
				dβ = dβ_dt(network, momentstruct, __k, algebra=algebra)
				dα = dα_dt(momentstruct, dμ, n_order, __k, algebra=algebra)
				nck = binomial_coeff(n_order, __k)
				__sign = sign_coeff(n_order, __k)
				__β = raw_to_cen(momentstruct, __k)
				dcen[counter] += nck*__sign*(α*dβ+__β*dα)
			end
		end
	end
	algebra.expand.(dcen)
end


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## load packages
using LIKWID

## measure optimistic peakflops (AVX2 FMA or AVX512 FMA if available)
L1_cache_size = LIKWID.get_cpu_topology().cacheLevels[1].size ÷ 1024 # in kB

cpuinfo = LIKWID.get_cpu_info()
if occursin("AVX512", cpuinfo.features)
  likwid_bench_kernel = "peakflops_avx512_fma"
elseif occursin("AVX2", cpuinfo.features)
  likwid_bench_kernel = "peakflops_avx_fma"
else
  likwid_bench_kernel = "peakflops_sse_fma"
end

max_flops_string = read(`likwid-bench -t $likwid_bench_kernel -W N:$(L1_cache_size)kB:1`, String)
max_flops = parse(Float64, match(r"(MFlops/s:\s+)(\d+\.\d+)", max_flops_string).captures[2]) / 1024
@info "PeakPerformance" max_flops

## prepare IO
filename = joinpath(@__DIR__, "peak_performance.dat")
open(filename, "w") do io
  println(io, "# NumericalFlux\tTreeMesh\tStructuredMesh\tP4estMesh\tPeakPerformance")
end

## gather data for volume terms implemented in Trixi.jl
measured_string = read(`likwid-perfctr -C 0 -g MEM_DP -m $(Base.julia_cmd()) --check-bounds=no --threads=1 $(joinpath(@__DIR__, "measure_volume_terms.jl"))`, String)

for numerical_flux in ["flux_shima_etal_turbo", "flux_ranocha_turbo"]
  offset = findfirst("Region TreeMesh-$numerical_flux", measured_string) |> last
  m = match(r"(DP \[MFLOP/s\]\s+\|\s+)(\d+\.\d+)", measured_string, offset)
  flops_TreeMesh = parse(Float64, m.captures[2]) / 1024
  @info "TreeMesh" numerical_flux flops_TreeMesh

  offset = findfirst("Region StructuredMesh-$numerical_flux", measured_string) |> last
  m = match(r"(DP \[MFLOP/s\]\s+\|\s+)(\d+\.\d+)", measured_string, offset)
  flops_StructuredMesh = parse(Float64, m.captures[2]) / 1024
  @info "StructuredMesh" numerical_flux flops_StructuredMesh

  offset = findfirst("Region P4estMesh-$numerical_flux", measured_string) |> last
  m = match(r"(DP \[MFLOP/s\]\s+\|\s+)(\d+\.\d+)", measured_string, offset)
  flops_P4estMesh = parse(Float64, m.captures[2]) / 1024
  @info "P4estMesh" numerical_flux flops_P4estMesh

  # save results to disk
  open(filename, "a") do io
    println(io, numerical_flux, "\t",
                flops_TreeMesh, "\t",
                flops_StructuredMesh, "\t",
                flops_P4estMesh, "\t",
                max_flops)
    println(io, numerical_flux, "\t",
                flops_TreeMesh / max_flops, "\t",
                flops_StructuredMesh / max_flops, "\t",
                flops_P4estMesh / max_flops, "\t",
                max_flops / max_flops)
  end
end


text:
using BoardGamesGUI
using Test

@testset "BoardGamesGUI.jl" begin
    # Write your tests here.
end


text:
import Distributions: MvNormal, ZeroVector, PDiagMat, invcov

# NOTE: see http://qwone.com/~jason/writing/multivariateNormal.pdf
#   for derivations of gradients wrt μ and U:
#       ∇log P(z|ϕ)

"""
A transformation which is a MultivariateNormal generated using either of the equations:

dist = N(μ, σ²)  where σ is a vector of standard deviations
dist = N(μ, Σ)   where U is an upper-triangular matrix with covariance Σ = U'U

The output is sampled from this distribution: z ~ dist

`transform!` will first update the underlying distribution with input vector ϕ = (μ, U)
and then call `rand!(dist, output)`
"""
type MvNormalTransformation{T,DIST<:MvNormal} <: Transformation
    dist::DIST
    n::Int
    nμ::Int
    nU::Int
    z̄::Vector{T}         # scratch space to avoid allocation in grad!
    input::SumNode{T,1}    # the sufficient stats: ϕ = vec(μ, U)
    output::OutputNode{T,1}  # the random sample
end

function MvNormalTransformation{T}(::Type{T}, args...)
    dist = MvNormal(args...)
    n = length(dist.μ)

    # do we need inputs for μ?
    nμ = isa(dist.μ, ZeroVector) ? 0 : n

    # do we have a diagonal or upper-triangular input for Σ?
    nU = typeof(dist.Σ) <: PDiagMat ? n : div(n*(n+1), 2)

    MvNormalTransformation{T,typeof(dist)}(
        dist,
        n,
        nμ,
        nU,
        zeros(T, n),
        InputNode(T, nμ+nU),
        OutputNode(T, n)
    )
end
MvNormalTransformation(args...) = MvNormalTransformation(Float64, args...)

# update the distribution params, then randomly sample into the output node
function transform!(mv::MvNormalTransformation)
    # update the distribution using ϕ
    transform_μ!(mv, mv.dist.μ)
    transform_Σ!(mv, mv.dist.Σ)

    # output: sample from the distribution
    rand!(mv.dist, output_value(mv))

    # @show mv.dist output_value(mv)
end

# # update the gradient w.r.t. inputs ϕ = vec(μ, U)
# function grad!(mv::MvNormalTransformation)
#     Σ⁻¹ = grad_Σ!(mv, mv.dist.Σ)
#     grad_μ!(mv, mv.dist.μ, Σ⁻¹)
# end

# -----------------------------------------------------------------------
# these methods update the underlying storage of the distribution
# using the input vector

# don't do anything for zeros
function transform_μ!(mv::MvNormalTransformation, μ::ZeroVector)
    return
end

# just copy
function transform_μ!(mv::MvNormalTransformation, μ)
    copy!(μ, view(input_value(mv), 1:mv.n))
    return
end

# diagonal Σ = diagmat(σ²)
function transform_Σ!(mv::MvNormalTransformation, Σ::PDiagMat)
    input = input_value(mv)
    for i=1:mv.n
        σᵢ = compute_σ(input[mv.nμ+i])
        Σ.diag[i] = σᵢ^2
        Σ.inv_diag[i] = inv(σᵢ)^2
    end
end

# Σ = U'U   (where U is assumed upper-triangular)
# To accomplish, we update the upper triangle of the cholesky decomp of dist.Σ
function transform_Σ!(mv::MvNormalTransformation, Σ)
    input = input_value(mv)
    cf = Σ.chol.factors
    i = 1
    for r=1:mv.n, c=r:mv.n
        cf[r,c] = input[mv.nμ + i]
        i += 1
    end
end

# compute_σ(ϕᵢ) = exp(clamp(ϕᵢ, -1e1, 1e1))
compute_σ(ϕᵢ) = (ϕᵢ < 0 ? min(ϕᵢ,-1e-5) : max(ϕᵢ, 1e-5))

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

# NOTE: σ = exp(ϕ) so that we map to reasonable positive values

# NOTE: this computes the "grad log prob": ∇log P(z | ϕ)
#   and thus it will not backprop gradients.
#   It is meant to be used only for policy gradient methods for now!
# function gradlogprob!{T}(∇logP::AbstractVector, mv::MvNormalTransformation{T})
function grad!{T}(mv::MvNormalTransformation{T})
    ϕ = input_value(mv)
    ∇ϕ = input_grad(mv)
    z = output_value(mv)
    # ∇z = output_grad(mv)  # not used!!
    z̄ = mv.z̄
    nμ = mv.nμ

    # # update: z̄ = (z - μ)
    # copy!(z̄, z)
    # for i=1:nμ
    #     z̄[i] = ϕ[i] - z̄[i]
    # end
    # scalar = T(2) * (sqrt(norm(z̄)) - one(T))

    if typeof(mv.dist.Σ) <: PDiagMat

        # for i=1:mv.n
        #     # demean z
        #     z̄[i] = z[i] - (nμ > 0 ? ϕ[i] : zero(T))
        #
        #     # ensure we're not dividing by really small numbers
        #     ϕ[nμ+i] = max(1e-2, ϕ[nμ+i])
        # end

        @assert nμ > 0
        # note that μ == ϕ[1:nμ], σ == exp(ϕ[nμ+1:end])
        ϵ = 1e-6

        μ = view(ϕ, 1:nμ)
        ∇μ = view(∇ϕ, 1:nμ)
        s = view(ϕ, nμ+1:2nμ)
        ∇s = view(∇ϕ, nμ+1:2nμ)
        for i=1:nμ
            σ = compute_σ(s[i])
            ∇μ[i] = (z[i] - μ[i]) / σ^2
            ∇s[i] = -(∇μ[i] - one(T)) / σ
            # ∇σ = -(one(T) / σ - ∇μ[i]^2) / T(2)
            # ∇s[i] = ∇σ * σ
            # ∇s[i] = ∇μ[i] * (z[i] - μ[i]) - one(T)
        end
        # for i=1:nμ
        #     j = nμ+i  # index of σ
        #     ∇ϕ[i] = (z[i] - ϕ[i]) / (exp(ϕ[j])^2 + ϵ)
        #     ∇ϕ[j] = -(one(T) / (ϕ[j] + ϵ) - ∇ϕ[i]^2) / T(2)
        #     # ∇ϕ[j] = (z[i] - ϕ[i])^2 / (ϕ[j]^2 + ϵ)
        # end
        # @show nμ mv.n z z̄ ϕ ∇ϕ

        # # do update for diagonal
        # # note: diag(U) .* diag(Σ⁻¹) == 1 ./ diag(U)
        # for i=1:mv.n
        #     # if i <= nμ
        #     #     ∇ϕ[i] = T(2) *  mv.dist.Σ.inv_diag[i] * z̄[i]
        #     # end
        #     ∇ϕ[nμ+i] = scalar / ϕ[nμ+i]
        # end
        # # @show scalar z̄ ϕ ∇ϕ mv.dist.Σ.inv_diag
    # else
    #     # do update for upper-triangular
    #     Σ⁻¹ = try
    #         invcov(mv.dist)
    #     catch err
    #         warn("Error in invcov: $err")
    #         return
    #     end
    #
    #     U = UpperTriangular(mv.dist.Σ.chol.factors)
    #     # @show U Σ⁻¹ typeof(Σ⁻¹)
    #
    #     # compute gradient of μ: ∇μ = -2 Σ⁻¹ (z-μ)
    #     if nμ > 0
    #         A_mul_B!(view(∇ϕ, 1:nμ), Σ⁻¹, z̄)
    #     end
    #
    #     # compute gradient of U: ∇U = (2-2√‖z-μ‖) U Σ⁻¹
    #     i = 1
    #     for r=1:mv.n, c=r:mv.n
    #         ∇ϕᵢ = zero(T)
    #         for j=c:mv.n
    #             ∇ϕᵢ += U[r,j] * Σ⁻¹[j,c]
    #         end
    #         ∇ϕ[nμ+i] = scalar * ∇ϕᵢ
    #         i += 1
    #     end
    end

    # @show mv.n nμ ϕ ∇ϕ z z̄ scalar
    return
end


text:
struct Secant{T, Line}
    factor::T
    shift::T
    ls::Line
end
"""
    Secant(;factor, shift)

Construct a method instance of `Secant` using the two required keywords `factor`
and `shift`. The two numbers are used to construct a pertubed iterate `x'` such
that:

    x' = x*factor + shift

where x is the initial iterate. See section 1.3 in [1] for an explanation.

References
---
[1] Kelley, C. T. (2003). Solving nonlinear equations with Newton's method (Vol. 1). Siam.
"""
Secant(; factor, shift, linesearch=Static()) = Secant(factor, shift, linesearch)

function solve(prob::NLEProblem, x0, method=Secant(0.99, 1e-1), options=NLEOptions())
    if length(x0) > 1
        throw(ArgumentError("You cannot solve your problem with the Secant method, since it is a univariate method, and `length(x0)>1` where `x0` is the initial iterate."))
    end

    T = eltype(x0)
    f_abstol, f_reltol, maxiter = options.f_abstol, options.f_reltol, options.maxiter

    # Evaluate
    xnext = x0
    Fnext = value(prob, xnext)

    # define the stopping tolerance
    f_stoptol = T(f_abstol) + T(f_reltol)*Fnext

    # Check for initial convergence
    if norm(Fnext) ≤ f_stoptol
        return (F=Fnext, x=xnext, code=:success, iter=0)
    end

    # construct the shifted iterate and evaluate the value of the residual
    xcurr = method.factor*x0+method.shift
    Fcurr = value(prob, xcurr)

    for i = 1:maxiter
        # Compute the difference
        ΔF = (Fnext-Fcurr)

        # Try to rescue the approximation
        if ΔF < eps(T)
            # Try to rescue
            xcurr = xnext
            Fcurr = Fnext

            xnext = xcurr*T(method.factor) + T(method.shift)
            Fnext = value(prob, xnext)

            ΔF = Fnext-Fcurr

            if ΔF < eps(T)
                # Failed to rescue
                return (F=Fnext, x=xnext, code=:failure_ΔF, iter=i)
            end
        end

        # calculate change in x between past two iterates
        s = - Fcurr*(xnext - xcurr)/ΔF
        Δx = s

        xcurr = xnext
        Fcurr = Fnext

        xnext = xcurr + Δx
        Fnext = value(prob, xnext)

        if norm(Fnext) ≤ f_stoptol
            return (F=Fnext, x=xnext, code=:success, iter=i)
        end
    end
    (F=Fnext, x=xnext, code=:failure_itermax, iter=options.maxiter)
end


text:
XYTuple{T} = Tuple{Real, Real}

+(p1::XYTuple, p2::XYTuple) = (p1[1] + p2[1], p1[2] + p2[2])
-(p1::XYTuple, p2::XYTuple) = (p1[1] - p2[1], p1[2] - p2[2])
*(p::XYTuple, c::Real) = (c * p[1], c * p[2])
*(c::Real, p::XYTuple) = (c * p[1], c * p[2])
/(p::XYTuple, c::Real) = (p[1] / c, p[2] / c)

"""
    partition(start, stop, N)

Generates N equally spaced, centered points between start and stop.
"""
function partition(start::XYTuple, stop::XYTuple, N::Integer)
    dist = stop -start
    step = dist / N
    [start + (i - 0.5) * step for i in 1:N]
end