text:
module ExpressionExplorer
export compute_symbolreferences, try_compute_symbolreferences, compute_usings_imports, SymbolsState, FunctionName, join_funcname_parts

import ..PlutoRunner
import Markdown
import Base: union, union!, ==, push!

###
# TWO STATE OBJECTS
###

# TODO: use GlobalRef instead
FunctionName = Array{Symbol,1}

struct FunctionNameSignaturePair
    name::FunctionName
    canonicalized_head::Any
end

Base.:(==)(a::FunctionNameSignaturePair, b::FunctionNameSignaturePair) = a.name == b.name && a.canonicalized_head == b.canonicalized_head
Base.hash(a::FunctionNameSignaturePair, h::UInt) = hash(a.name, hash(a.canonicalized_head, h))

"SymbolsState trickles _down_ the ASTree: it carries referenced and defined variables from endpoints down to the root."
Base.@kwdef mutable struct SymbolsState
    references::Set{Symbol} = Set{Symbol}()
    assignments::Set{Symbol} = Set{Symbol}()
    funccalls::Set{FunctionName} = Set{FunctionName}()
    funcdefs::Dict{FunctionNameSignaturePair,SymbolsState} = Dict{FunctionNameSignaturePair,SymbolsState}()
end

"ScopeState moves _up_ the ASTree: it carries scope information up towards the endpoints."
mutable struct ScopeState
    inglobalscope::Bool
    exposedglobals::Set{Symbol}
    hiddenglobals::Set{Symbol}
    definedfuncs::Set{Symbol}
end
ScopeState() = ScopeState(true, Set{Symbol}(), Set{Symbol}(), Set{Symbol}())

# The `union` and `union!` overloads define how two `SymbolsState`s or two `ScopeState`s are combined.

function union(a::Dict{FunctionNameSignaturePair,SymbolsState}, bs::Dict{FunctionNameSignaturePair,SymbolsState}...)
    union!(Dict{FunctionNameSignaturePair,SymbolsState}(), a, bs...)
end

function union!(a::Dict{FunctionNameSignaturePair,SymbolsState}, bs::Dict{FunctionNameSignaturePair,SymbolsState}...)
    for b in bs
        for (k, v) in b
            if haskey(a, k)
                a[k] = union!(a[k], v)
            else
                a[k] = v
            end
        end
        a
    end
    return a
end

function union(a::SymbolsState, b::SymbolsState)
    SymbolsState(a.references ∪ b.references, a.assignments ∪ b.assignments, a.funccalls ∪ b.funccalls, a.funcdefs ∪ b.funcdefs)
end

function union!(a::SymbolsState, bs::SymbolsState...)
    union!(a.references, (b.references for b in bs)...)
    union!(a.assignments, (b.assignments for b in bs)...)
    union!(a.funccalls, (b.funccalls for b in bs)...)
    union!(a.funcdefs, (b.funcdefs for b in bs)...)
    return a
end

function union!(a::Tuple{FunctionName,SymbolsState}, bs::Tuple{FunctionName,SymbolsState}...)
    a[1], union!(a[2], (b[2] for b in bs)...)
end

function union(a::ScopeState, b::ScopeState)
    SymbolsState(a.inglobalscope && b.inglobalscope, a.exposedglobals ∪ b.exposedglobals, a.hiddenglobals ∪ b.hiddenglobals)
end

function union!(a::ScopeState, bs::ScopeState...)
    a.inglobalscope &= all((b.inglobalscope for b in bs)...)
    union!(a.exposedglobals, (b.exposedglobals for b in bs)...)
    union!(a.hiddenglobals, (b.hiddenglobals for b in bs)...)
    union!(a.definedfuncs, (b.definedfuncs for b in bs)...)
    return a
end

function ==(a::SymbolsState, b::SymbolsState)
    a.references == b.references && a.assignments == b.assignments && a.funccalls == b.funccalls && a.funcdefs == b.funcdefs
end

Base.push!(x::Set) = x

###
# HELPER FUNCTIONS
###

function explore_inner_scoped(ex::Expr, scopestate::ScopeState)::SymbolsState
    # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.
    innerscopestate = deepcopy(scopestate)
    innerscopestate.inglobalscope = false

    return mapfoldl(a -> explore!(a, innerscopestate), union!, ex.args, init=SymbolsState())
end

# from the source code: https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm#L9
const modifiers = [:(+=), :(-=), :(*=), :(/=), :(//=), :(^=), :(÷=), :(%=), :(<<=), :(>>=), :(>>>=), :(&=), :(⊻=), :(≔), :(⩴), :(≕)]
const modifiers_dotprefixed = [Symbol('.' * String(m)) for m in modifiers]

function will_assign_global(assignee::Symbol, scopestate::ScopeState)::Bool
    (scopestate.inglobalscope || assignee ∈ scopestate.exposedglobals) && (assignee ∉ scopestate.hiddenglobals || assignee ∈ scopestate.definedfuncs)
end

function will_assign_global(assignee::Array{Symbol,1}, scopestate::ScopeState)::Bool
    if length(assignee) == 0
        false
    elseif length(assignee) > 1
        scopestate.inglobalscope
    else
        will_assign_global(assignee[1], scopestate)
    end
end

function get_global_assignees(assignee_exprs, scopestate::ScopeState)::Set{Symbol}
    global_assignees = Set{Symbol}()
    for ae in assignee_exprs
        if isa(ae, Symbol)
            will_assign_global(ae, scopestate) && push!(global_assignees, ae)
        else
            if ae.head == :(::)
                will_assign_global(ae.args[1], scopestate) && push!(global_assignees, ae.args[1])
            else
                @warn "Unknown assignee expression" ae
            end
        end
    end
    return global_assignees
end

function get_assignees(ex::Expr)::FunctionName
    if ex.head == :tuple
        # e.g. (x, y) in the ex (x, y) = (1, 23)
        union!(Symbol[], get_assignees.(ex.args)...)
        # filter(s->s isa Symbol, ex.args)
    elseif ex.head == :(::)
        # TODO: type is referenced
        Symbol[ex.args[1]]
    elseif ex.head == :ref || ex.head == :(.)
        Symbol[]
    else
        @warn "unknown use of `=`. Assignee is unrecognised." ex
        Symbol[]
    end
end

# e.g. x = 123, but ignore _ = 456
get_assignees(ex::Symbol) = all_underscores(ex) ? Symbol[] : Symbol[ex]

# When you assign to a datatype like Int, String, or anything bad like that
# e.g. 1 = 2
# This is parsable code, so we have to treat it
get_assignees(::Any) = Symbol[]

all_underscores(s::Symbol) = all(isequal('_'), string(s))

# TODO: this should return a FunctionName, and use `split_funcname`.
"Turn :(A{T}) into :A."
function uncurly!(ex::Expr, scopestate::ScopeState)::Symbol
    @assert ex.head == :curly
    push!(scopestate.hiddenglobals, (a for a in ex.args[2:end] if a isa Symbol)...)
    Symbol(ex.args[1])
end

uncurly!(ex::Expr)::Symbol = ex.args[1]

uncurly!(s::Symbol, scopestate=nothing)::Symbol = s

"Turn `:(Base.Submodule.f)` into `[:Base, :Submodule, :f]` and `:f` into `[:f]`."
function split_funcname(funcname_ex::Expr)::FunctionName
    if funcname_ex.head == :(.)
        vcat(split_funcname.(funcname_ex.args)...)
    else
        # a call to a function that's not a global, like calling an array element: `funcs[12]()`
        # TODO: explore symstate!
        Symbol[]
    end
end

function split_funcname(funcname_ex::QuoteNode)::FunctionName
    split_funcname(funcname_ex.value)
end
function split_funcname(funcname_ex::GlobalRef)::FunctionName
    split_funcname(funcname_ex.name)
end

function split_funcname(funcname_ex::Symbol)::FunctionName
    Symbol[funcname_ex |> without_dotprefix |> without_dotsuffix]
end

function is_just_dots(ex::Expr)
    ex.head == :(.) && all(is_just_dots, ex.args)
end
is_just_dots(::Union{QuoteNode,Symbol,GlobalRef}) = true
is_just_dots(::Any) = false

# this includes GlobalRef - it's fine that we don't recognise it, because you can't assign to a globalref?
function split_funcname(::Any)::FunctionName
    Symbol[]
end

"""Turn `Symbol(".+")` into `:(+)`"""
function without_dotprefix(funcname::Symbol)::Symbol
    fn_str = String(funcname)
    if length(fn_str) > 0 && fn_str[1] == '.'
        Symbol(fn_str[2:end])
    else
        funcname
    end
end

"""Turn `Symbol("sqrt.")` into `:sqrt`"""
function without_dotsuffix(funcname::Symbol)::Symbol
    fn_str = String(funcname)
    if length(fn_str) > 0 && fn_str[end] == '.'
        Symbol(fn_str[1:end - 1])
    else
        funcname
    end
end

"""Turn `Symbol[:Module, :func]` into Symbol("Module.func").

This is **not** the same as the expression `:(Module.func)`, but is used to identify the function name using a single `Symbol` (like normal variables).
This means that it is only the inverse of `ExpressionExplorer.split_funcname` iff `length(parts) ≤ 1`."""
function join_funcname_parts(parts::FunctionName)::Symbol
	join(parts .|> String, ".") |> Symbol
end

# this is stupid -- désolé
function is_joined_funcname(joined::Symbol)
    occursin('.', String(joined))
end

assign_to_kw(e::Expr) = e.head == :(=) ? Expr(:kw, e.args...) : e
assign_to_kw(x::Any) = x

"Turn `A[i] * B[j,K[l+m]]` into `A[0] * B[0,K[0+0]]` to hide loop indices"
function strip_indexing(x, inside::Bool=false)
    if Meta.isexpr(x, :ref)
        Expr(:ref, strip_indexing(x.args[1]), strip_indexing.(x.args[2:end], true)...)
    elseif Meta.isexpr(x, :call)
        Expr(x.head, x.args[1], strip_indexing.(x.args[2:end], inside)...)
    elseif x isa Symbol && inside
        0
    else
        x
    end
end

###
# MAIN RECURSIVE FUNCTION
###

# Spaghetti code for a spaghetti problem 🍝

# Possible leaf: value
# Like: a = 1
# 1 is a value (Int64)
function explore!(value, scopestate::ScopeState)::SymbolsState
    # includes: LineNumberNode, Int64, String,
    return SymbolsState()
end

# Possible leaf: symbol
# Like a = x
# x is a symbol
# We handle the assignment separately, and explore!(:a, ...) will not be called.
# Therefore, this method only handles _references_, which are added to the symbolstate, depending on the scopestate.
function explore!(sym::Symbol, scopestate::ScopeState)::SymbolsState
    if sym ∈ scopestate.hiddenglobals
        SymbolsState()
    else
        SymbolsState(references=Set([sym]))
    end
end

"""
Returns whether or not an assignment Expr(:(=),...) is assigning to a new function
  * f(x) = ...
  * f(x)::V = ...
  * f(::T) where {T} = ...
"""
is_function_assignment(ex::Expr) = ex.args[1] isa Expr && (ex.args[1].head == :call || ex.args[1].head == :where || (ex.args[1].head == :(::) && ex.args[1].args[1] isa Expr && ex.args[1].args[1].head == :call))

# General recursive method. Is never a leaf.
# Modifies the `scopestate`.
function explore!(ex::Expr, scopestate::ScopeState)::SymbolsState
    if ex.head == :(=)
        # Does not create scope

        if is_function_assignment(ex)
            # f(x, y) = x + y
            # Rewrite to:
            # function f(x, y) x + y end
            return explore!(Expr(:function, ex.args...), scopestate)
        end

        val = ex.args[2]
        # Handle generic types assignments A{B} = C{B, Int}
        if ex.args[1] isa Expr && ex.args[1].head == :curly
            assignees, symstate = explore_funcdef!(ex.args[1], scopestate)
            innersymstate = union!(symstate, explore!(val, scopestate))
        else
            assignees = get_assignees(ex.args[1])
            symstate = innersymstate = explore!(val, scopestate)
        end

        global_assignees = get_global_assignees(assignees, scopestate)

        # If we are _not_ assigning a global variable, then this symbol hides any global definition with that name
        push!(scopestate.hiddenglobals, setdiff(assignees, global_assignees)...)
        assigneesymstate = explore!(ex.args[1], scopestate)

        push!(scopestate.hiddenglobals, global_assignees...)
        push!(symstate.assignments, global_assignees...)
        push!(symstate.references, setdiff(assigneesymstate.references, global_assignees)...)
        filter!(!all_underscores, symstate.references)  # Never record _ as a reference

        return symstate
    elseif ex.head in modifiers
        # We change: a[1] += 123
        # to:        a[1] = a[1] + 123
        # We transform the modifier back to its operator
        # for when users redefine the + function

        operator = let
            s = string(ex.head)
            Symbol(s[1:prevind(s, lastindex(s))])
        end
        expanded_expr = Expr(:(=), ex.args[1], Expr(:call, operator, ex.args[1], ex.args[2]))
        return explore!(expanded_expr, scopestate)
    elseif ex.head in modifiers_dotprefixed
        # We change: a[1] .+= 123
        # to:        a[1] .= a[1] + 123

        operator = Symbol(string(ex.head)[2:end - 1])
        expanded_expr = Expr(:(.=), ex.args[1], Expr(:call, operator, ex.args[1], ex.args[2]))
        return explore!(expanded_expr, scopestate)
    elseif ex.head == :let || ex.head == :for || ex.head == :while
        # Creates local scope
        return explore_inner_scoped(ex, scopestate)
    elseif ex.head == :filter
        # In a filter, the assignment is the second expression, the condition the first
        return mapfoldr(a -> explore!(a, scopestate), union!, ex.args, init=SymbolsState())
    elseif ex.head == :generator
        # Creates local scope

        # In a `generator`, a single expression is followed by the iterator assignments.
        # In a `for`, this expression comes at the end.

        # This is not strictly the normal form of a `for` but that's okay
        return explore!(Expr(:for, ex.args[2:end]..., ex.args[1]), scopestate)
    elseif ex.head == :macrocall
        # Does not create sccope
        new_ex = maybe_macroexpand(ex)

        newnew_ex = Meta.isexpr(new_ex, :macrocall) ? Expr(:call, new_ex.args...) : new_ex

        return explore!(newnew_ex, scopestate)
    elseif ex.head == :call
        # Does not create scope

        if is_just_dots(ex.args[1])
            funcname = ex.args[1] |> split_funcname
            symstate = if length(funcname) == 0
                explore!(ex.args[1], scopestate)
            elseif length(funcname) == 1
                if funcname[1] ∈ scopestate.hiddenglobals
                    SymbolsState()
                else
                    SymbolsState(funccalls=Set{FunctionName}([funcname]))
                end
            else
                SymbolsState(references=Set{Symbol}([funcname[1]]), funccalls=Set{FunctionName}([funcname]))
            end
            # Explore code inside function arguments:
            union!(symstate, explore!(Expr(:block, ex.args[2:end]...), scopestate))
            return symstate
        else
            return explore!(Expr(:block, ex.args...), scopestate)
        end
    elseif ex.head == :kw
        return explore!(ex.args[2], scopestate)
    elseif ex.head == :struct
        # Creates local scope

        structnameexpr = ex.args[2]
        structfields = ex.args[3].args

        equiv_func = Expr(:function, Expr(:call, structnameexpr, structfields...), Expr(:block, nothing))

        # struct should always be in Global state
        globalscopestate = deepcopy(scopestate)
        globalscopestate.inglobalscope = true

        # we register struct definitions as both a variable and a function. This is because deleting a struct is trickier than just deleting its methods.
        inner_symstate = explore!(equiv_func, globalscopestate)

        structname = first(keys(inner_symstate.funcdefs)).name |> join_funcname_parts
        push!(inner_symstate.assignments, structname)
        return inner_symstate
    elseif ex.head == :abstract
        equiv_func = Expr(:function, ex.args...)
        inner_symstate = explore!(equiv_func, scopestate)

        abstracttypename = first(keys(inner_symstate.funcdefs)).name |> join_funcname_parts
        push!(inner_symstate.assignments, abstracttypename)
        return inner_symstate
    elseif ex.head == :function || ex.head == :macro
        symstate = SymbolsState()
        # Creates local scope

        funcroot = ex.args[1]

        # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.
        innerscopestate = deepcopy(scopestate)
        innerscopestate.inglobalscope = false

        funcname, innersymstate = explore_funcdef!(funcroot, innerscopestate)
        # Macro are called using @funcname, but defined with funcname. We need to change that in our scopestate
        # (The `!= 0` is for when the function named couldn't be parsed)
        if ex.head == :macro && length(funcname) != 0
            setdiff!(innerscopestate.hiddenglobals, funcname)
            funcname = Symbol[Symbol("@$(funcname[1])")]
            push!(innerscopestate.hiddenglobals, funcname...)
        end

        union!(innersymstate, explore!(Expr(:block, ex.args[2:end]...), innerscopestate))

        funcnamesig = FunctionNameSignaturePair(funcname, canonalize(funcroot))

        if will_assign_global(funcname, scopestate)
            symstate.funcdefs[funcnamesig] = innersymstate
            if length(funcname) == 1
                push!(scopestate.definedfuncs, funcname[end])
                push!(scopestate.hiddenglobals, funcname[end])
            elseif length(funcname) > 1
                push!(symstate.references, funcname[end - 1]) # reference the module of the extended function
            end
        else
            # The function is not defined globally. However, the function can still modify the global scope or reference globals, e.g.

            # let
            #     function f(x)
            #         global z = x + a
            #     end
            #     f(2)
            # end

            # so we insert the function's inner symbol state here, as if it was a `let` block.
            symstate = innersymstate
        end

        return symstate
    elseif ex.head == :try
        symstate = SymbolsState()

        # Handle catch first
        if ex.args[3] != false
            union!(symstate, explore_inner_scoped(ex.args[3], scopestate))
            # If we catch a symbol, it could shadow a global reference, remove it
            if ex.args[2] != false
                setdiff!(symstate.references, Symbol[ex.args[2]])
            end
        end

        # Handle the try block
        union!(symstate, explore_inner_scoped(ex.args[1], scopestate))

        # Finally, handle finally
        if length(ex.args) == 4
            union!(symstate, explore_inner_scoped(ex.args[4], scopestate))
        end

        return symstate
    elseif ex.head == :(->)
        # Creates local scope

        tempname = Symbol("anon", rand(UInt64))

        # We will rewrite this to a normal function definition, with a temporary name
        funcroot = ex.args[1]
        args_ex = if funcroot isa Symbol || (funcroot isa Expr && funcroot.head == :(::))
            [funcroot]
        elseif funcroot.head == :tuple || funcroot.head == :(...) || funcroot.head == :block
            funcroot.args
        else
            @error "Unknown lambda type"
        end

        equiv_func = Expr(:function, Expr(:call, tempname, args_ex...), ex.args[2])

        return explore!(equiv_func, scopestate)
    elseif ex.head == :global
        # Does not create scope

        # We have one of:
        # global x;
        # global x = 1;
        # global x += 1;

        # where x can also be a tuple:
        # global a,b = 1,2

        globalisee = ex.args[1]

        if isa(globalisee, Symbol)
            push!(scopestate.exposedglobals, globalisee)
            return SymbolsState()
        elseif isa(globalisee, Expr)
            # temporarily set inglobalscope to true
            old = scopestate.inglobalscope
            scopestate.inglobalscope = true
            result = explore!(globalisee, scopestate)
            scopestate.inglobalscope = old
            return result
        else
            @error "unknown global use" ex
            return explore!(globalisee, scopestate)
        end

        return symstate
    elseif ex.head == :local
        # Does not create scope

        localisee = ex.args[1]

        if isa(localisee, Symbol)
            push!(scopestate.hiddenglobals, localisee)
            return SymbolsState()
        elseif isa(localisee, Expr) && (localisee.head == :(=) || localisee.head in modifiers)
            push!(scopestate.hiddenglobals, get_assignees(localisee.args[1])...)
            return explore!(localisee, scopestate)
        else
            @warn "unknown local use" ex
            return explore!(localisee, scopestate)
        end
    elseif ex.head == :tuple
        # Does not create scope

        # There are three (legal) cases:
        # 1. Creating a tuple:
        #   (a, b, c)

        # 2. Creating a named tuple:
        #   (a=1, b=2, c=3)

        # 3. Multiple assignments
        # a,b,c = 1,2,3
        # This parses to:
        # head = :tuple
        # args = [:a, :b, :(c=1), :2, :3]
        #
        # 🤔
        # we turn it into two expressions:
        #
        # (a, b) = (2, 3)
        # (c = 1)
        #
        # and explore those :)

        indexoffirstassignment = findfirst(a -> isa(a, Expr) && a.head == :(=), ex.args)
        if indexoffirstassignment !== nothing
            # we have one of two cases, see next `if`
            indexofsecondassignment = findnext(a -> isa(a, Expr) && a.head == :(=), ex.args, indexoffirstassignment + 1)

            if length(ex.args) == 1 || indexofsecondassignment !== nothing
                # 2.
                # we have a named tuple, e.g. (a=1, b=2)
                new_args = map(ex.args) do a
                    (a isa Expr && a.head == :(=)) ? a.args[2] : a
                end
                return explore!(Expr(:block, new_args...), scopestate)
            else
                # 3.
                # we have a tuple assignment, e.g. `a, (b, c) = [1, [2, 3]]`
                before = ex.args[1:indexoffirstassignment - 1]
                after = ex.args[indexoffirstassignment + 1:end]

                symstate_middle = explore!(ex.args[indexoffirstassignment], scopestate)
                symstate_outer = explore!(Expr(:(=), Expr(:tuple, before...), Expr(:block, after...)), scopestate)

                return union!(symstate_middle, symstate_outer)
            end
        else
            # 1.
            # good ol' tuple
            return explore!(Expr(:block, ex.args...), scopestate)
        end
    elseif ex.head == :(.) && ex.args[2] isa Expr && ex.args[2].head == :tuple
        # pointwise function call, e.g. sqrt.(nums)
        # we rewrite to a regular call

        return explore!(Expr(:call, ex.args[1], ex.args[2].args...), scopestate)
    elseif ex.head == :using || ex.head == :import
        imports = if ex.args[1].head == :(:)
            ex.args[1].args[2:end]
        else
            ex.args
        end

        packagenames = map(e -> e.args[end], imports)

        return SymbolsState(assignments=Set{Symbol}(packagenames))
    elseif ex.head == :quote
        # We ignore contents

        return SymbolsState()
    elseif ex.head == :module
        # We ignore contents; the module name is a definition

        return SymbolsState(assignments=Set{Symbol}([ex.args[2]]))
    else
        # fallback, includes:
        # begin, block, do, toplevel, const
        # (and hopefully much more!)

        # Does not create scope (probably)

        return mapfoldl(a -> explore!(a, scopestate), union!, ex.args, init=SymbolsState())
    end
end

"Return the function name and the SymbolsState from argument defaults. Add arguments as hidden globals to the `scopestate`.

Is also used for `struct` and `abstract`."
function explore_funcdef!(ex::Expr, scopestate::ScopeState)::Tuple{FunctionName,SymbolsState}
    if ex.head == :call
        params_to_explore = ex.args[2:end]
        # Using the keyword args syntax f(;y) the :parameters node is the first arg in the AST when it should
        # be explored last. We change from (parameters, ...) to (..., parameters)
        if length(params_to_explore) >= 2 && params_to_explore[1] isa Expr && params_to_explore[1].head == :parameters
            params_to_explore = [params_to_explore[2:end]..., params_to_explore[1]]
        end

        # Handle struct as callables, `(obj::MyType)(a, b) = ...`
        # or `function (obj::MyType)(a, b) ...; end` by rewriting it as:
        # function MyType(obj, a, b) ...; end
        funcroot = ex.args[1]
        if funcroot isa Expr && funcroot.head == :(::)
            return explore_funcdef!(Expr(:call, reverse(funcroot.args)..., params_to_explore...), scopestate)
        end

        # get the function name
        name, symstate = explore_funcdef!(funcroot, scopestate)
        # and explore the function arguments
        return mapfoldl(a -> explore_funcdef!(a, scopestate), union!, params_to_explore, init=(name, symstate))
    elseif ex.head == :(::) || ex.head == :kw || ex.head == :(=)
        # account for unnamed params, like in f(::Example) = 1
        if ex.head == :(::) && length(ex.args) == 1
            symstate = explore!(ex.args[1], scopestate)

            return Symbol[], symstate
        end

        # recurse
        name, symstate = explore_funcdef!(ex.args[1], scopestate)
        if length(ex.args) > 1
            # use `explore!` (not `explore_funcdef!`) to explore the argument's default value - these can contain arbitrary expressions
            union!(symstate, explore!(ex.args[2], scopestate))
        end
        return name, symstate

    elseif ex.head == :where
        # function(...) where {T, S <: R, U <: A.B}
        # supertypes `R` and `A.B` are referenced
        supertypes_symstate = SymbolsState()
        for a in ex.args[2:end]
            name, inner_symstate = explore_funcdef!(a, scopestate)
            if length(name) == 1
                push!(scopestate.hiddenglobals, name[1])
            end
            union!(supertypes_symstate, inner_symstate)
        end
        # recurse
        name, symstate = explore_funcdef!(ex.args[1], scopestate)
        union!(symstate, supertypes_symstate)
        return name, symstate

    elseif ex.head == :(<:)
        # for use in `struct` and `abstract`
        name = uncurly!(ex.args[1], scopestate)
        symstate = if length(ex.args) == 1
            SymbolsState()
        else
            explore!(ex.args[2], scopestate)
        end
        return Symbol[name], symstate

    elseif ex.head == :curly
        name = uncurly!(ex, scopestate)
        return Symbol[name], SymbolsState()

    elseif ex.head == :parameters || ex.head == :tuple
        return mapfoldl(a -> explore_funcdef!(a, scopestate), union!, ex.args, init=(Symbol[], SymbolsState()))

    elseif ex.head == :(.)
        return split_funcname(ex), SymbolsState()

    elseif ex.head == :(...)
        return explore_funcdef!(ex.args[1], scopestate)
    else
        return Symbol[], explore!(ex, scopestate)
    end
end

function explore_funcdef!(ex::QuoteNode, scopestate::ScopeState)::Tuple{FunctionName,SymbolsState}
    explore_funcdef!(ex.value, scopestate)
end

function explore_funcdef!(ex::Symbol, scopestate::ScopeState)::Tuple{FunctionName,SymbolsState}
    push!(scopestate.hiddenglobals, ex)
    Symbol[ex |> without_dotprefix |> without_dotsuffix], SymbolsState()
end

function explore_funcdef!(::Any, ::ScopeState)::Tuple{FunctionName,SymbolsState}
    Symbol[], SymbolsState()
end

const can_macroexpand_no_bind = Set(Symbol.(["@md_str", "Markdown.@md_str", "@gensym", "Base.@gensym", "@kwdef", "Base.@kwdef", "@enum", "Base.@enum", "@cmd"]))
const can_macroexpand = can_macroexpand_no_bind ∪ Set(Symbol.(["@bind", "PlutoRunner.@bind"]))

macro_kwargs_as_kw(ex::Expr) = Expr(:macrocall, ex.args[1:3]..., assign_to_kw.(ex.args[4:end])...)

function symbolics_mockexpand(s::Any)
    # goofy implementation of the syntax described in https://symbolics.juliasymbolics.org/dev/manual/variables/
    if Meta.isexpr(s, :ref, 2)
        :($(s.args[1]) = $(s.args[2]))
    elseif Meta.isexpr(s, :call, 2)
        second = s.args[2] === Symbol("..") ? 123 : s.args[2]
        :($(symbolics_mockexpand(s.args[1])); $(second) = 123)
    elseif s isa Symbol
        :($(s) = 123)
    else
        nothing
    end
end

is_symbolics_arg(s) = symbolics_mockexpand(s) !== nothing

maybe_untuple(es) = if length(es) == 1 && Meta.isexpr(first(es), :tuple)
    first(es).args
else
    es
end

"""
If the macro is known to Pluto, expand or 'mock expand' it, if not, return the expression.

Macros can transform the expression into anything - the best way to treat them is to `macroexpand`. The problem is that the macro is only available on the worker process, see https://github.com/fonsp/Pluto.jl/issues/196
"""
function maybe_macroexpand(ex::Expr; recursive=false, expand_bind=true)
    result = if ex.head === :macrocall
        funcname = ex.args[1] |> split_funcname
        funcname_joined = join_funcname_parts(funcname)

        args = ex.args[3:end]

        if funcname_joined ∈ (expand_bind ? can_macroexpand : can_macroexpand_no_bind)
            expanded = macroexpand(PlutoRunner, ex; recursive=false)
            Expr(:call, ex.args[1], expanded)

        elseif !isempty(args) && Meta.isexpr(args[1], :(:=))
            ex = macro_kwargs_as_kw(ex)
            # macros like @einsum C[i] := A[i,j] are assignment to C, illegal syntax without macro
            ein = args[1]
            left = if Meta.isexpr(ein.args[1], :ref)
                # assign to the symbol, and save LHS indices as fake RHS argument
                ex = Expr(ex.head, ex.args..., Expr(:ref, :Float64, ein.args[1].args[2:end]...))
                ein.args[1].args[1]
            else
                ein.args[1]  # scalar case `c := A[i,j]`
            end
            ein_done = Expr(:(=), left, strip_indexing.(ein.args[2:end])...)  # i,j etc. are local
            Expr(:call, ex.args[1:2]..., ein_done, strip_indexing.(ex.args[4:end])...)

        elseif !isempty(args) && funcname_joined === Symbol("@ode_def")
            if args[1] isa Symbol
                :($(args[1]) = @ode_def 123)
            else
                :(@ode_def)
            end
        elseif !isempty(args) && (funcname_joined === Symbol("@functor") || funcname_joined === Symbol("Flux.@functor"))
            Expr(:macrocall, ex.args[1:2]..., :($(args[1]) = 123), ex.args[4:end]...)
        elseif !isempty(args) && (funcname_joined === Symbol("@variables") || funcname_joined === Symbol("Symbolics.@variables")) && all(is_symbolics_arg, maybe_untuple(args))
            Expr(:macrocall, ex.args[1:2]..., symbolics_mockexpand.(maybe_untuple(args))...)
        # elseif length(ex.args) >= 4 && (funcname_joined === Symbol("@variable") || funcname_joined === Symbol("JuMP.@variable"))
        #     if Meta.isexpr(ex.args[4], :comparison)
        #         parts = ex.args[4].args[1:2:end]
        #         if length(parts) == 2
        #         foldl(parts) do (e,next)
        #             :($(e) = $(next))
        #         end
        #     elseif Meta.isexpr(ex.args[4], :block)

        #     end

        #     Expr(:macrocall, ex.args[1:3]..., )
            # add more macros here
        elseif length(args) ≥ 2 && ex.args[1] != GlobalRef(Core, Symbol("@doc"))
            # for macros like @test a ≈ b atol=1e-6, read assignment in 2nd & later arg as keywords
            macro_kwargs_as_kw(ex)
        else
            ex
        end
    else
        ex
    end

    if recursive && (result isa Expr)
        Expr(result.head, maybe_macroexpand.(result.args; recursive=recursive, expand_bind=expand_bind)...)
    else
        result
    end
end

maybe_macroexpand(ex::Any; kwargs...) = ex

###
# CANONICALIZE FUNCTION DEFINITIONS
###

"""
Turn a function definition expression (`Expr`) into a "canonical" form, in the sense that two methods that would evaluate to the same method signature have the same canonical form. Part of a solution to https://github.com/fonsp/Pluto.jl/issues/177. Such a canonical form cannot be achieved statically with 100% correctness (impossible), but we can make it good enough to be practical.

# Wait, "evaluate to the same method signature"?

In Pluto, you cannot do definitions of **the same global variable** in different cells. This is needed for reactivity to work, and it avoid ambiguous notebooks and stateful stuff. This rule used to also apply to functions: you had to place all methods of a function in one cell. (Go and read more about methods in Julia if you haven't already.) But this is quite annoying, especially because multiple dispatch is so important in Julia code. So we allow methods of the same function to be defined across multiple cells, but we still want to throw errors when you define **multiple methods with the same signature**, because one overrides the other. For example:
```julia
julia> f(x) = 1
f (generic function with 1 method)

julia> f(x) = 2
f (generic function with 1 method)
``

After adding the second method, the function still has only 1 method. This is because the second definition overrides the first one, instead of being added to the method table. This example should be illegal in Julia, for the same reason that `f = 1` and `f = 2` is illegal. So our problem is: how do we know that two cells will define overlapping methods?

Ideally, we would just evaluate the user's code and **count methods** afterwards, letting Julia do the work. Unfortunately, we need to know this info _before_ we run cells, otherwise we don't know in which order to run a notebook! There are ways to break this circle, but it would complicate our process quite a bit.

Instead, we will do _static analysis_ on the function definition expressions to determine whether they overlap. This is non-trivial. For example, `f(x)` and `f(y::Any)` define the same method. Trickier examples are here: https://github.com/fonsp/Pluto.jl/issues/177#issuecomment-645039993

# Wait, "function definition expressions"?
For example:

```julia
e = :(function f(x::Int, y::String)
        x + y
    end)

dump(e, maxdepth=2)

#=
gives:

Expr
  head: Symbol function
  args: Array{Any}((2,))
    1: Expr
    2: Expr
=#
```

This first arg is the function head:

```julia
e.args[1] == :(f(x::Int, y::String))
```

# Mathematics
Our problem is to find a way to compute the equivalence relation ~ on `H × H`, with `H` the set of function head expressions, defined as:

`a ~ b` iff evaluating both expressions results in a function with exactly one method.

_(More precisely, evaluating `Expr(:function, x, Expr(:block))` with `x ∈ {a, b}`.)_

The equivalence sets are isomorphic to the set of possible Julia methods.

Instead of finding a closed form algorithm for `~`, we search for a _canonical form_: a function `canonical: H -> H` that chooses one canonical expression per equivalence class. It has the property

`canonical(a) = canonical(b)` implies `a ~ b`.

We use this **canonical form** of the function's definition expression as its "signature". We compare these canonical forms when determining whether two function expressions will result in overlapping methods.

# Example
```julia
e1 = :(f(x, z::Any))
e2 = :(g(x, y))

canonalize(e1) == canonalize(e2)
```

```julia
e1 = :(f(x))
e2 = :(g(x, y))

canonalize(e1) != canonalize(e2)
```

```julia
e1 = :(f(a::X, b::wow(ie), c,      d...; e=f) where T)
e2 = :(g(z::X, z::wow(ie), z::Any, z...     ) where T)

canonalize(e1) == canonalize(e2)
```
"""
function canonalize(ex::Expr)
	if ex.head == :where
		Expr(:where, canonalize(ex.args[1]), ex.args[2:end]...)
	elseif ex.head == :call || ex.head == :tuple
		skip_index = ex.head == :call ? 2 : 1
		# ex.args[1], if ex.head == :call this is the function name, we dont want it

		interesting = filter(ex.args[skip_index:end]) do arg
			!(arg isa Expr && arg.head == :parameters)
		end

		hide_argument_name.(interesting)
    elseif ex.head == :(::)
        canonalize(ex.args[1])
    elseif ex.head == :curly || ex.head == :(<:)
        # for struct definitions, which we hackily treat as functions
        nothing
    else
		@error "Failed to canonalize this strange looking function" ex
		nothing
	end
end

# for `function g end`
canonalize(::Symbol) = nothing

function hide_argument_name(ex::Expr)
    if ex.head == :(::) && length(ex.args) > 1
        Expr(:(::), nothing, ex.args[2:end]...)
    elseif ex.head == :(...)
        Expr(:(...), hide_argument_name(ex.args[1]))
    elseif ex.head == :kw
        Expr(:kw, hide_argument_name(ex.args[1]), nothing)
    else
        ex
    end
end
hide_argument_name(::Symbol) = Expr(:(::), nothing, :Any)
hide_argument_name(x::Any) = x

###
# UTILITY FUNCTIONS
###

"Get the global references, assignment, function calls and function defintions inside an arbitrary expression."
function compute_symbolreferences(ex::Any)::SymbolsState
    symstate = explore!(ex, ScopeState())

    # We do something special to account for recursive functions:
    # If a function `f` calls a function `g`, and both are defined inside this cell, the reference to `g` inside the symstate of `f` will be deleted.
    # The motivitation is that normally, an assignment (or function definition) will add that symbol to a list of 'hidden globals' - any future references to that symbol will be ignored. i.e. the _local definition hides a global_.
    # In the case of functions, you can reference functions and variables that do not yet exist, and so they won't be in the list of hidden symbols when the function definition is analysed.
    # Of course, our method will fail if a referenced function is defined both inside the cell **and** in another cell. However, this will lead to a MultipleDefinitionError before anything bad happens.
    for (func, inner_symstate) in symstate.funcdefs
        inner_symstate.references = setdiff(inner_symstate.references, keys(symstate.funcdefs))
        inner_symstate.funccalls = setdiff(inner_symstate.funccalls, keys(symstate.funcdefs))
    end
    symstate
end

function try_compute_symbolreferences(ex::Any)::SymbolsState
	try
		compute_symbolreferences(ex)
	catch e
		@error "Expression explorer failed on: " ex
		showerror(stderr, e, stacktrace(catch_backtrace()))
		SymbolsState(references=Set{Symbol}([:fake_reference_to_prevent_it_from_looking_like_a_text_only_cell]))
	end
end

Base.@kwdef struct UsingsImports
    usings::Set{Expr}=Set{Expr}()
    imports::Set{Expr}=Set{Expr}()
end

# Performance analysis: https://gist.github.com/fonsp/280f6e883f419fb3a59231b2b1b95cab
"Preallocated version of [`compute_usings_imports`](@ref)."
function compute_usings_imports!(out::UsingsImports, ex::Any)
    if isa(ex, Expr)
        if ex.head == :using
			push!(out.usings, ex)
		elseif ex.head == :import
			push!(out.imports, ex)
        elseif ex.head != :quote
			for a in ex.args
				compute_usings_imports!(out, a)
			end
        end
    end
	out
end

"""
Given `:(using Plots, Something.Else, .LocalModule)`, return `Set([:Plots, :Something])`.
"""
function external_package_names(ex::Expr)::Set{Symbol}
	@assert ex.head == :import || ex.head == :using
	if Meta.isexpr(ex.args[1], :(:))
		external_package_names(Expr(ex.head, ex.args[1].args[1]))
	else
		out = Set{Symbol}()
		for a in ex.args
			if Meta.isexpr(a, :(.))
				if a.args[1] != :(.)
					push!(out, a.args[1])
				end
			end
		end
		out
	end
end

function external_package_names(x::UsingsImports)::Set{Symbol}
    union!(Set{Symbol}(), external_package_names.(x.usings)..., external_package_names.(x.imports)...)
end

"Get the sets of `using Module` and `import Module` subexpressions that are contained in this expression."
compute_usings_imports(ex) = compute_usings_imports!(UsingsImports(), ex)

"Return whether the expression is of the form `Expr(:toplevel, LineNumberNode(..), any)`."
function is_toplevel_expr(ex::Expr)::Bool
    (ex.head == :toplevel) && (length(ex.args) == 2) && (ex.args[1] isa LineNumberNode)
end

is_toplevel_expr(::Any)::Bool = false

"If the expression is a (simple) assignemnt at its root, return the assignee as `Symbol`, return `nothing` otherwise."
function get_rootassignee(ex::Expr, recurse::Bool=true)::Union{Symbol,Nothing}
    if is_toplevel_expr(ex) && recurse
        get_rootassignee(ex.args[2], false)
    elseif ex.head == :(=) && ex.args[1] isa Symbol
        ex.args[1]
    else
        nothing
    end
end

get_rootassignee(ex::Any, recuse::Bool=true)::Union{Symbol,Nothing} = nothing

"Is this code simple enough that we can wrap it inside a function to boost performance? Look for [`PlutoRunner.Computer`](@ref) to learn more."
function can_be_function_wrapped(x::Expr)
    if x.head === :global || # better safe than sorry
        x.head === :using ||
        x.head === :import ||
        x.head === :module ||
        x.head === :function ||
        x.head === :macro ||
        x.head === :macrocall || # we might want to get rid of this one, but that requires some work
        x.head === :struct ||
        x.head === :abstract ||
        (x.head === :(=) && is_function_assignment(x)) || # f(x) = ...
        (x.head === :call && (x.args[1] === :eval || x.args[1] === :include))
        false
    else
        all(can_be_function_wrapped, x.args)
    end

end
can_be_function_wrapped(x::Any) = true

end


text:
name(s::NMR.Spectrum) = sampleinfo(s)["Name"]
MW(s::NMR.Spectrum) = parse(Float64, sampleinfo(s)["MW"])
blank_conc(s::Spectrum) = parse(Float64, sampleinfo(s)["Blank conc. (mg/mL)"])

function protons(s::NMR.Spectrum)
    pr = get(sampleinfo(s), "Protons", "")
    parsed = eval(Meta.parse(pr))
    parsed !== nothing ? float.(parsed) : nothing
end

function quant_threshold(s::NMR.Spectrum)
    th = get(sampleinfo(s), "Quantitation threshold (%)", "0.0")
    parse(Float64, th)
end

# (internal standard) in the Name field signifies the internal standard.
isstandard(s::NMR.Spectrum) = occursin("(internal standard)", name(s))

# A spectrum is dark if its title doesn't contain the Proton field
# or if it is the internal standard.
isdark(s::NMR.Spectrum) = !haskey(sampleinfo(s), "Protons") || isstandard(s)

"""
    library(path)

Returns a tuple from a valid directory `path`: `lib`, an array of library spectra; `darks`, an array of 'dark' spectra; and `is`, the internal standard spectrum.
"""
function library(path)
    try
        expnos = [expno for expno in readdir(path) if tryparse(Int, expno) !== nothing]
        sort!(expnos, by=expno->parse(Int,expno))
        lib = [NMR.Spectrum(joinpath(path, expno),1) for expno in expnos]
        for (i,l) in enumerate(lib)
            pr = protons(l)
            if pr === nothing
                continue
            end
            l1 = length(NMR.intrng_shifts(l))
            l2 = length(pr)
            if l1 != l2
                error("Reference $(name(l)) ($i) contains $l1 integration regions != $l2 proton groups.")
            end
        end
        is = filter(isstandard, lib)
        if length(is) == 0
            error("$path: No internal standard found.")
        end
        if length(is) != 1
            error("One internal standard expected; found $(length(is)): $(name.(is))")
        end
        darks = filter(isdark, lib)
        lib, darks, is[1]
    catch e
        println(sprint(showerror, e)*"\nIgnoring directory $path.")
        return nothing
    end
end

function libraries(path)
    result = Dict( p => library(joinpath(path, p))
        for p in readdir(path) if isdir(joinpath(path, p)) )
    return filter(x -> x.second != nothing, result)
end

function getlib(s::Spectrum, libs)
    libname = get(sampleinfo(s), "Library", "D2O")
    libs[libname]
end

"""
    sampleinfo(s::NMR.Spectrum)

Returns a Dictionary of sample information obtained from the title of NMR spectrum `s`.
"""
function sampleinfo(s::NMR.Spectrum)
    result = Dict{String,String}()
    for l in split(NMR.title(s), '\n')
        if findfirst("=", l) !== nothing
            k,v = split(l, '='; limit = 2)
            result[strip(k)] = strip(v)
        end
    end
    result
end

function sampleinfovalid(s::NMR.Spectrum, libs)
    sinfo = sampleinfo(s)

    for inputfield in ["Sample weight (mg)", "Volume (mL)", "Blank conc. (mg/mL)", "Library"]
        if !haskey(sinfo, inputfield)
            error("$inputfield missing from Spectrum Title")
        end
    end

    for inputfield in ["Sample weight (mg)", "Volume (mL)", "Blank conc. (mg/mL)"]
        if tryparse(Float64, sinfo[inputfield]) == nothing
            error("Spectrum Title: Invalid $(inputfield)!")
        end
    end

    if !(sinfo["Library"] ∈ collect(keys(libs)))
        error("Spectrum Title: Invalid 'Library'")
    end

    return true
end

function statistics(percentages)
    m = mean(percentages)
    pm = (maximum(percentages)-minimum(percentages))/2
    if m > 10.0
        @sprintf "%.1f" m
    else
        @sprintf "%.2f" m
    end
end

ppm_format(δ) = @sprintf "%.2f" δ

score_format(score) = @sprintf "%.3f" score

findlibversion(lib) = get(sampleinfo(lib[1]), "Library version", "Missing")

findversion(name::String) = string(Pkg.installed()[name])

zonedtime(time::DateTime=now(), f=dateformat"yyyy-mm-dd H\hM\mS\s ZZZZ") = Dates.format(ZonedDateTime(time,localzone()), f)

function ansitohtml(text::String)
    rd = Dict(r"\e\[0m" => "",
    r"\e\[1m" => "",
    r"\e\[22m" => "",
    r"\e\[37m" => "",
    r"\e\[39m" => "",
    r"\e\[91m" => "",
    r"\e\[32m" => "",
    r"\e\[36m" => "")

    for i in rd
        text = replace(text, i)
    end
    text
end

function restorehtml(text::String)
    dict = Dict(r"<" => "<",
    r"=" => "=",
    r""" => "\"",
    r">" => ">",
    r"/" => "/")

    for i in dict
        text = replace(text, i)
    end
    text
end

ext: jl


text:
# Note: This file is executed from `Pkg` in an isolated environment.
# Warn: Since this is executed in a sandbox, the current preferences are not
#       propagated and thus no decision can be made based on them.
#       Use LazyArtifacts in this case.
#       Fixed in https://github.com/JuliaLang/Pkg.jl/pull/2920

push!(Base.LOAD_PATH, dirname(@__DIR__))
using TOML, Artifacts, Base.BinaryPlatforms
include("./platform_augmentation.jl")
artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")

# Get "target triplet" from ARGS, if given (defaulting to the host triplet otherwise)
target_triplet = get(ARGS, 1, Base.BinaryPlatforms.host_triplet())

# Augment this platform object with any special tags we require
platform = augment_platform!(HostPlatform(parse(Platform, target_triplet)))

# Select all downloadable artifacts that match that platform
artifacts = select_downloadable_artifacts(artifacts_toml; platform, include_lazy = true)

#Output the result to `stdout` as a TOML dictionary
TOML.print(stdout, artifacts)

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
using CUDA.CUSOLVER

using LinearAlgebra
using LinearAlgebra: BlasInt

m = 15
n = 10
l = 13
k = 1

@test_throws ArgumentError CUSOLVER.cusolverjob('M')

@testset "elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
    @testset "Cholesky (po)" begin
        A    = rand(elty,n,n)
        A    = A*A' #posdef
        B    = rand(elty,n,n)
        d_A  = CuArray(A)
        d_B  = CuArray(B)

        d_F  = cholesky(d_A, Val(false))
        F    = cholesky(A, Val(false))
        @test F.U   ≈ collect(d_F.U)
        @test F\(A'B) ≈ collect(d_F\(d_A'd_B))

        d_F  = cholesky(Hermitian(d_A, :L), Val(false))
        F    = cholesky(Hermitian(A, :L), Val(false))
        @test F.L   ≈ collect(d_F.L)
        @test F\(A'B) ≈ collect(d_F\(d_A'd_B))

        @test_throws DimensionMismatch LinearAlgebra.LAPACK.potrs!('U',d_A,CUDA.rand(elty,m,m))

        A    = rand(elty,m,n)
        d_A  = CuArray(A)
        @test_throws DimensionMismatch cholesky(d_A)
        @test_throws DimensionMismatch LinearAlgebra.LAPACK.potrs!('U',d_A,d_B)

        A    = zeros(elty,n,n)
        d_A  = CuArray(A)
        @test_throws LinearAlgebra.PosDefException cholesky(d_A)
    end

    CUDA.CUSOLVER.version() >= v"10.1" && @testset "Cholesky inverse (potri)" begin
        # test lower
        A    = rand(elty,n,n)
        A    = A*A' #posdef
        d_A  = CuArray(A)

        LinearAlgebra.LAPACK.potrf!('L', A)
        LinearAlgebra.LAPACK.potrf!('L', d_A)

        LinearAlgebra.LAPACK.potri!('L', A)
        LinearAlgebra.LAPACK.potri!('L', d_A)
        @test A  ≈ collect(d_A)

        # test upper
        A    = rand(elty,n,n)
        A    = A*A' #posdef
        d_A  = CuArray(A)

        LinearAlgebra.LAPACK.potrf!('U', A)
        LinearAlgebra.LAPACK.potrf!('U', d_A)
        LinearAlgebra.LAPACK.potri!('U', A)
        LinearAlgebra.LAPACK.potri!('U', d_A)
        @test A  ≈ collect(d_A)
    end

    @testset "getrf!" begin
        A          = rand(elty,m,n)
        d_A        = CuArray(A)
        d_A,d_ipiv = CUSOLVER.getrf!(d_A)
        h_A        = collect(d_A)
        h_ipiv     = collect(d_ipiv)
        alu        = LinearAlgebra.LU(h_A, convert(Vector{BlasInt},h_ipiv), zero(BlasInt))
        @test A ≈ Array(alu)
        A    = zeros(elty,n,n)
        d_A  = CuArray(A)
        @test_throws LinearAlgebra.SingularException CUSOLVER.getrf!(d_A)
    end

    @testset "getrs!" begin
        A          = rand(elty,n,n)
        d_A        = CuArray(A)
        d_A,d_ipiv = CUSOLVER.getrf!(d_A)
        B          = rand(elty,n,n)
        d_B        = CuArray(B)
        d_B        = CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
        h_B        = collect(d_B)
        @test h_B  ≈ A\B
        A          = rand(elty,m,n)
        d_A        = CuArray(A)
        @test_throws DimensionMismatch CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
        A          = rand(elty,n,n)
        d_A        = CuArray(A)
        B          = rand(elty,m,n)
        d_B        = CuArray(B)
        @test_throws DimensionMismatch CUSOLVER.getrs!('N',d_A,d_ipiv,d_B)
    end

    @testset "Matrix division" begin
        A = rand(elty,n,n)
        B = rand(elty,n,n)
        C = A \ B
        d_A = CuArray(A)
        d_B = CuArray(B)
        @test C ≈ Array(d_A \ d_B)
    end

    @testset "geqrf!" begin
        A         = rand(elty,m,n)
        d_A       = CuArray(A)
        d_A,d_tau = CUSOLVER.geqrf!(d_A)
        h_A       = collect(d_A)
        h_tau     = collect(d_tau)
        qra       = LinearAlgebra.QR(h_A, h_tau)
        @test A ≈ Array(qra)
    end

    @testset "ormqr!" begin
        A          = rand(elty, m, n)
        d_A        = CuArray(A)
        d_A, d_tau = CUSOLVER.geqrf!(d_A)
        B          = rand(elty, n, l)
        d_B        = CuArray(B)
        d_B        = CUSOLVER.ormqr!('L', 'N', d_A, d_tau, d_B)
        h_B        = collect(d_B)
        F          = qr!(A)
        @test h_B  ≈ Array(F.Q)*B
        A          = rand(elty, n, m)
        d_A        = CuArray(A)
        d_A, d_tau = CUSOLVER.geqrf!(d_A)
        B          = rand(elty, n, l)
        d_B        = CuArray(B)
        d_B        = CUSOLVER.ormqr!('L', 'N', d_A, d_tau, d_B)
        h_B        = collect(d_B)
        F          = qr!(A)
        @test h_B  ≈ Array(F.Q)*B
        A          = rand(elty, m, n)
        d_A        = CuArray(A)
        d_A, d_tau = CUSOLVER.geqrf!(d_A)
        B          = rand(elty, l, m)
        d_B        = CuArray(B)
        d_B        = CUSOLVER.ormqr!('R', 'N', d_A, d_tau, d_B)
        h_B        = collect(d_B)
        F          = qr!(A)
        @test h_B  ≈ B*Array(F.Q)
        A          = rand(elty, n, m)
        d_A        = CuArray(A)
        d_A, d_tau = CUSOLVER.geqrf!(d_A)
        B          = rand(elty, l, n)
        d_B        = CuArray(B)
        d_B        = CUSOLVER.ormqr!('R', 'N', d_A, d_tau, d_B)
        h_B        = collect(d_B)
        F          = qr!(A)
        @test h_B  ≈ B*Array(F.Q)
    end

    @testset "orgqr!" begin
        A         = rand(elty,n,m)
        d_A       = CuArray(A)
        d_A,d_tau = CUSOLVER.geqrf!(d_A)
        d_Q       = CUSOLVER.orgqr!(d_A, d_tau)
        h_Q       = collect(d_Q)
        F         = qr!(A)
        @test h_Q ≈ Array(F.Q)
        A         = rand(elty,m,n)
        d_A       = CuArray(A)
        d_A,d_tau = CUSOLVER.geqrf!(d_A)
        d_Q       = CUSOLVER.orgqr!(d_A, d_tau)
        h_Q       = collect(d_Q)
        F         = qr!(A)
        @test h_Q ≈ Array(F.Q)
    end

    @testset "sytrf!" begin
        A          = rand(elty,n,n)
        A          = A + A' #symmetric
        d_A        = CuArray(A)
        d_A,d_ipiv = CUSOLVER.sytrf!('U',d_A)
        h_A        = collect(d_A)
        h_ipiv     = collect(d_ipiv)
        A, ipiv    = LAPACK.sytrf!('U',A)
        @test ipiv == h_ipiv
        @test A ≈ h_A
        A    = rand(elty,m,n)
        d_A  = CuArray(A)
        @test_throws DimensionMismatch CUSOLVER.sytrf!('U',d_A)
        A    = zeros(elty,n,n)
        d_A  = CuArray(A)
        @test_throws LinearAlgebra.SingularException CUSOLVER.sytrf!('U',d_A)
    end

    @testset "gebrd!" begin
        A                             = rand(elty,m,n)
        d_A                           = CuArray(A)
        d_A, d_D, d_E, d_TAUQ, d_TAUP = CUSOLVER.gebrd!(d_A)
        h_A                           = collect(d_A)
        h_D                           = collect(d_D)
        h_E                           = collect(d_E)
        h_TAUQ                        = collect(d_TAUQ)
        h_TAUP                        = collect(d_TAUP)
        A,d,e,q,p                     = LAPACK.gebrd!(A)
        @test A ≈ h_A
        @test d ≈ h_D
        @test e[min(m,n)-1] ≈ h_E[min(m,n)-1]
        @test q ≈ h_TAUQ
        @test p ≈ h_TAUP
    end

    @testset "syevd!" begin
        A              = rand(elty,m,m)
        A             += A'
        d_A            = CuArray(A)
        local d_W, d_V
        if( elty <: Complex )
            d_W, d_V   = CUSOLVER.heevd!('V','U', d_A)
        else
            d_W, d_V   = CUSOLVER.syevd!('V','U', d_A)
        end
        h_W            = collect(d_W)
        h_V            = collect(d_V)
        Eig            = eigen(A)
        @test Eig.values ≈ h_W
        @test abs.(Eig.vectors'*h_V) ≈ I
        d_A            = CuArray(A)
        if( elty <: Complex )
            d_W   = CUSOLVER.heevd!('N','U', d_A)
        else
            d_W   = CUSOLVER.syevd!('N','U', d_A)
        end
        h_W            = collect(d_W)
        @test Eig.values ≈ h_W
    end

    @testset "sygvd!" begin
        A              = rand(elty,m,m)
        B              = rand(elty,m,m)
        A             *= A'
        B             *= B'
        d_A            = CuArray(A)
        d_B            = CuArray(B)
        local d_W, d_VA, d_VB
        if( elty <: Complex )
            d_W, d_VA, d_VB = CUSOLVER.hegvd!(1, 'V','U', d_A, d_B)
        else
            d_W, d_VA, d_VB = CUSOLVER.sygvd!(1, 'V','U', d_A, d_B)
        end
        h_W            = collect(d_W)
        h_VA           = collect(d_VA)
        h_VB           = collect(d_VB)
        Eig            = eigen(Hermitian(A), Hermitian(B))
        @test Eig.values ≈ h_W
        @test A*h_VA ≈ B*h_VA*Diagonal(h_W) rtol=1e-4
        # test normalization condition for eigtype 1
        @test abs.(h_VA'B*h_VA) ≈ Matrix(one(elty)*I, m, m)
        d_A            = CuArray(A)
        d_B            = CuArray(B)
        if( elty <: Complex )
            d_W   = CUSOLVER.hegvd!(1, 'N','U', d_A, d_B)
        else
            d_W   = CUSOLVER.sygvd!(1, 'N','U', d_A, d_B)
        end
        h_W            = collect(d_W)
        @test Eig.values ≈ h_W
        d_B            = CUDA.rand(elty, m+1, m+1)
        if( elty <: Complex )
            @test_throws DimensionMismatch CUSOLVER.hegvd!(1, 'N','U', d_A, d_B)
        else
            @test_throws DimensionMismatch CUSOLVER.sygvd!(1, 'N','U', d_A, d_B)
        end
    end

    @testset "syevj!" begin
        A              = rand(elty,m,m)
        B              = rand(elty,m,m)
        A             *= A'
        B             *= B'
        d_A            = CuArray(A)
        d_B            = CuArray(B)
        local d_W, d_VA, d_VB
        if( elty <: Complex )
            d_W, d_VA, d_VB = CUSOLVER.hegvj!(1, 'V','U', d_A, d_B)
        else
            d_W, d_VA, d_VB = CUSOLVER.sygvj!(1, 'V','U', d_A, d_B)
        end
        h_W            = collect(d_W)
        h_VA           = collect(d_VA)
        h_VB           = collect(d_VB)
        Eig            = eigen(Hermitian(A), Hermitian(B))
        @test Eig.values ≈ h_W
        @test A*h_VA ≈ B*h_VA*Diagonal(h_W) rtol=1e-4
        # test normalization condition for eigtype 1
        @test abs.(h_VA'B*h_VA) ≈ Matrix(one(elty)*I, m, m)
        d_A            = CuArray(A)
        d_B            = CuArray(B)
        if( elty <: Complex )
            d_W   = CUSOLVER.hegvj!(1, 'N','U', d_A, d_B)
        else
            d_W   = CUSOLVER.sygvj!(1, 'N','U', d_A, d_B)
        end
        h_W            = collect(d_W)
        @test Eig.values ≈ h_W
    end

    @testset "elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
        @testset "syevjBatched!" begin
            # Generate a random symmetric/hermitian matrix
            A = rand(elty, m,m,n)
            A += permutedims(A, (2,1,3))
            d_A = CuArray(A)

            # Run the solver
            local d_W, d_V
            if( elty <: Complex )
                d_W, d_V   = CUSOLVER.heevjBatched!('V','U', d_A)
            else
                d_W, d_V   = CUSOLVER.syevjBatched!('V','U', d_A)
            end

            # Pull it back to hardware
            h_W   = collect(d_W)
            h_V   = collect(d_V)

            # Use non-GPU blas to estimate the eigenvalues as well
            for i = 1:n
                # Get our eigenvalues
                Eig = eigen(LinearAlgebra.Hermitian(A[:,:,i]))

                # Compare to the actual ones
                @test Eig.values ≈ h_W[:,i]
                @test abs.(Eig.vectors'*h_V[:,:,i]) ≈ I
            end

            # Do it all again, but with the option to not compute eigenvectors
            d_A = CuArray(A)

            # Run the solver
            local d_W
            if( elty <: Complex )
                d_W   = CUSOLVER.heevjBatched!('N','U', d_A)
            else
                d_W   = CUSOLVER.syevjBatched!('N','U', d_A)
            end

            # Pull it back to hardware
            h_W   = collect(d_W)

            # Use non-GPU blas to estimate the eigenvalues as well
            for i = 1:n
                # Get the reference results
                Eig = eigen(LinearAlgebra.Hermitian(A[:,:,i]))

                # Compare to the actual ones
                @test Eig.values ≈ h_W[:,i]
            end
        end
    end

    @testset "svd with $method method" for
        method in (CUSOLVER.QRAlgorithm, CUSOLVER.JacobiAlgorithm),
        (_m, _n) in ((m, n), (n, m))

        A              = rand(elty, _m, _n)
        U, S, V        = svd(A, full=true)
        d_A            = CuArray(A)

        if _m > _n || method == CUSOLVER.JacobiAlgorithm
            d_U, d_S, d_V  = svd(d_A, method, full=true)
            h_S            = collect(d_S)
            h_U            = collect(d_U)
            h_V            = collect(d_V)
            @test abs.(h_U'h_U) ≈ I
            @test abs.(h_U[:,1:min(_m,_n)]'U[:,1:min(_m,_n)]) ≈ I
            @test collect(svdvals(d_A, method)) ≈ svdvals(A)
            @test abs.(h_V'*h_V) ≈ I
            @test abs.(h_V[:,1:min(_m,_n)]'*V[:,1:min(_m,_n)]) ≈ I
            @test collect(d_U'*d_A*d_V) ≈ U'*A*V
            @test collect(svd(d_A, method).V') == h_V[:,1:min(_m,_n)]'
        else
            @test_throws ArgumentError svd(d_A, method)
        end
    end
    # Check that constant propagation works
    _svd(A) = svd(A, CUSOLVER.QRAlgorithm)
    @inferred _svd(CUDA.rand(Float32, 4, 4))

    @testset "qr" begin
        tol = min(m, n)*eps(real(elty))*(1 + (elty <: Complex))

        A              = rand(elty, m, n)
        qra            = qr(A)
        d_A            = CuArray(A)
        d_F            = qr(d_A)
        d_RR           = d_F.Q'*d_A
        @test d_RR[1:n,:] ≈ d_F.R atol=tol*norm(A)
        @test norm(d_RR[n+1:end,:]) < tol*norm(A)
        @test size(d_F) == size(A)
        @test size(d_F.Q, 1) == size(A, 1)
        @test det(d_F.Q) ≈ det(collect(d_F.Q * CuMatrix{elty}(I, size(d_F.Q)))) atol=tol*norm(A)
        CUDA.@allowscalar begin
            qval = d_F.Q[1, 1]
            @test qval ≈ qra.Q[1, 1]
            qrstr = sprint(show, d_F)
            @test qrstr == "$(typeof(d_F)) with factors Q and R:\n$(sprint(show, d_F.Q))\n$(sprint(show, d_F.R))"
        end
        dQ, dR = d_F
        @test collect(dQ*dR) ≈ A
        A              = rand(elty, n, m)
        d_A            = CuArray(A)
        d_F            = qr(d_A)
        @test d_F.Q'*d_A ≈ d_F.R atol=tol*norm(A)
        @test det(d_F.Q) ≈ det(collect(d_F.Q * CuMatrix{elty}(I, size(d_F.Q)))) atol=tol*norm(A)
        A              = rand(elty, m, n)
        d_A            = CuArray(A)
        h_q, h_r       = qr(d_A)
        q, r           = qr(A)
        @test Array(h_q) ≈ Array(q)
        @test collect(CuArray(h_q)) ≈ Array(q)
        @test Array(h_r) ≈ Array(r)
        A              = rand(elty, n, m)
        d_A            = CuArray(A)
        h_q, h_r       = qr(d_A) # FixMe! Use iteration protocol when implemented
        q, r           = qr(A)
        @test Array(h_q) ≈ Array(q)
        @test Array(h_r) ≈ Array(r)
    end

    @testset "potrsBatched!" begin
        @testset "elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
            # Test lower
            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]
            bB = [rand(elty, m) for i in 1:n]

            # move to device
            bd_A = CuArray{elty, 2}[]
            bd_B = CuArray{elty, 1}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
                push!(bd_B, CuArray(bB[i]))
            end

            bd_X = CUSOLVER.potrsBatched!('L', bd_A, bd_B)
            bh_X = [collect(bd_X[i]) for i in 1:n]

            for i = 1:n
                LinearAlgebra.LAPACK.potrs!('L', bA[i], bB[i])
                @test bB[i] ≈ bh_X[i]
            end

            # Test upper
            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]
            bB = [rand(elty, m) for i in 1:n]

            # move to device
            bd_A = CuArray{elty, 2}[]
            bd_B = CuArray{elty, 1}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
                push!(bd_B, CuArray(bB[i]))
            end

            bd_X = CUSOLVER.potrsBatched!('U', bd_A, bd_B)
            bh_X = [collect(bd_X[i]) for i in 1:n]

            for i = 1:n
                LinearAlgebra.LAPACK.potrs!('U', bA[i], bB[i])
                @test bB[i] ≈ bh_X[i]
            end
            # error throwing tests
            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]
            bB = [rand(elty, m) for i in 1:n+1]

            # move to device
            bd_A = CuArray{elty, 2}[]
            bd_B = CuArray{elty, 1}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
                push!(bd_B, CuArray(bB[i]))
            end
            push!(bd_B, CuArray(bB[end]))

            @test_throws DimensionMismatch CUSOLVER.potrsBatched!('L', bd_A, bd_B)

            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]
            bB = [rand(elty, m) for i in 1:n]
            bB[1] = rand(elty, m+1)
            # move to device
            bd_A = CuArray{elty, 2}[]
            bd_B = CuArray{elty, 1}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
                push!(bd_B, CuArray(bB[i]))
            end

            @test_throws DimensionMismatch CUSOLVER.potrsBatched!('L', bd_A, bd_B)

            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]
            bB = [rand(elty, m, m) for i in 1:n]
            # move to device
            bd_A = CuArray{elty, 2}[]
            bd_B = CuArray{elty, 2}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
                push!(bd_B, CuArray(bB[i]))
            end

            @test_throws ArgumentError CUSOLVER.potrsBatched!('L', bd_A, bd_B)
        end
    end

    @testset "potrfBatched!" begin
        @testset "elty = $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
            # Test lower
            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]

            # move to device
            bd_A = CuArray{elty, 2}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
            end

            bd_A, info = CUSOLVER.potrfBatched!('L', bd_A)
            bh_A = [collect(bd_A[i]) for i in 1:n]

            for i = 1:n
                LinearAlgebra.LAPACK.potrf!('L', bA[i])
                @test bA[i] ≈ bh_A[i]
            end

            # Test upper
            bA = [rand(elty, m, m) for i in 1:n]
            bA = [bA[i]*bA[i]' for i in 1:n]

            # move to device
            bd_A = CuArray{elty, 2}[]
            for i in 1:length(bA)
                push!(bd_A, CuArray(bA[i]))
            end

            bd_A, info = CUSOLVER.potrfBatched!('U', bd_A)
            bh_A = [collect(bd_A[i]) for i in 1:n]

            for i = 1:n
                LinearAlgebra.LAPACK.potrf!('U', bA[i])
                # cuSOLVER seems to return symmetric/hermitian matrix when using 'U'
                @test Hermitian(bA[i]) ≈ bh_A[i]
            end
        end
    end

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module StratiGraphics

using GeoStatsBase

using ImageFiltering: Kernel, centered, imfilter!
using Distributions: Exponential, Uniform, wsample
using RecipesBase

import GeoStatsBase: preprocess, solve_single

include("state.jl")
include("processes.jl")
include("durations.jl")
include("environment.jl")
include("record.jl")
include("strata.jl")
include("geostats.jl")

export
  # geological environment
  Environment,
  simulate,

  # timeless processes
  GeoStatsProcess,
  SmoothingProcess,
  AnalyticalProcess,
  SequentialTimelessProcess,

  # duration processes
  ExponentialDuration,
  UniformDuration,

  # geological state
  LandState,

  # geological record
  Record,

  # stratigraphy
  Strata,
  voxelize,

  # GeoStats.jl API
  StratSim

end


text:
# Model of the Asap Dual Clocked Fifos.
@with_kw mutable struct DualClockFifo{T} <: AbstractFifo
    # Statistics
    num_reads  :: Int = 0
    num_writes :: Int = 0

    # Vector for holding data.
    buffer::Vector{T}

    # Pointers
    readside_rd_pointer :: Int = 1
    readside_wr_pointer :: Int = 1
    writeside_rd_pointer :: Int = 1
    writeside_wr_pointer :: Int = 1

    # Read side signals
    out_read_valid    :: Bool = false
    out_read_data     :: T    = zero(T)
    in_read_increment :: Bool = false

    # Write side signals
    in_write_data    :: T   = zero(T)
    in_write_request ::Bool = false
end

# Constructor
function DualClockFifo(::Type{T}, buffersize::Int) where T
    # Construct the buffer and fall-back to the default values provided
    # by the "with_kw" macro
    buffer = zeros(T, buffersize)
    return DualClockFifo{T}(buffer = buffer)
end

################################################################################
# Methods
################################################################################

function getoccupancy(w, r, buffersize)
    occupancy = w - r
    if w < r
        occupancy += buffersize
    end
    return occupancy
end

isfull(w, r, buffersize, reserve) = getoccupancy(w, r, buffersize) >= buffersize - reserve
buffersize(f::DualClockFifo) = length(f.buffer)
reserve(f::DualClockFifo) = 4

function read_occupancy(f::DualClockFifo)
    return getoccupancy(
        f.readside_wr_pointer,
        f.readside_rd_pointer,
        buffersize(f),
    )
end

function write_occupancy(f::DualClockFifo)
    return getoccupancy(
        f.writeside_wr_pointer,
        f.writeside_rd_pointer,
        buffersize(f),
    )
end

# ------- #
# Reading #
# ------- #
isreadready(f::DualClockFifo) = f.out_read_valid
isreadempty(f::DualClockFifo) = f.readside_wr_pointer == f.readside_rd_pointer

peek(f::DualClockFifo) = f.out_read_data
function read(f::DualClockFifo)
    data = peek(f)
    f.in_read_increment = true
    return data
end

function readupdate!(f::DualClockFifo)
    # Update write pointer from the write side of the fifo.
    f.readside_wr_pointer = f.writeside_wr_pointer

    if f.out_read_valid && f.in_read_increment
        f.num_reads += 1
    end

    if isreadempty(f) == false && (f.out_read_valid == false || f.in_read_increment == true)
        # Get data to read.
        f.out_read_data = f.buffer[f.readside_rd_pointer]
        f.out_read_valid = true

        # Increment read pointer.
        f.readside_rd_pointer += 1
        if f.readside_rd_pointer > buffersize(f)
            f.readside_rd_pointer = 1
        end

    elseif f.in_read_increment
        f.out_read_valid = false
    end

    f.in_read_increment = false
    return nothing
end

# ------- #
# Writing #
# ------- #

iswriteready(f::DualClockFifo) = !isfull(
    f.writeside_wr_pointer,
    f.writeside_rd_pointer,
    buffersize(f),
    reserve(f)
)
iswriteempty(f::DualClockFifo) = f.wr_wr_pointer == f.wr_rd_pointer

function write(fifo::DualClockFifo{T}, data::Integer) where T
    # Set the valid signal to true and save the data.
    fifo.in_write_request = true
    fifo.in_write_data = data
end

# Simulate a clock on the "write" side of the fifo.
function writeupdate!(f::DualClockFifo)
    # Update the read pointer from the read-side of the f.
    # TODO: Add delays to model metastability mitigation
    f.writeside_rd_pointer = f.readside_rd_pointer

    # Check if a write is pending. If so, push through the write.
    if f.in_write_request
        f.buffer[f.writeside_wr_pointer] = f.in_write_data
        # Increment the write pointer, wrapping if necessary.
        f.writeside_wr_pointer += 1
        if f.writeside_wr_pointer > buffersize(f)
            f.writeside_wr_pointer = 1
        end

        # Increment counter
        f.num_writes += 1
    end
    # Reset the pending write flag.
    f.in_write_request = false
end


text:
using Base.Test

using NURBS

include("BSplineBasis.jl")
include("NURBSBasis.jl")


text:
using Dates
using PyPlot, Colors

# don't run it in same console as Bolivia-overview plot

include("startup.jl")
include("manualMassLibrary.jl")

fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/interpolatedData_May_halfHour/"
ptr3Traces = Float64.(readdlm("$(fp)ptr3.csv")[3:end, 3:end])
ptr3unixTimes = Float64.(readdlm("$(fp)ptr3.csv")[3:end, 1])
(ptr3Peaks, peakHeader) = readdlm("$(fp)ptr3peaks.csv", header = true)
composition = Int32.(ptr3Peaks[:,3:end])
names = String.(ptr3Peaks[:,2])
masses = Float64.(ptr3Peaks[:,1])

excludeCHON = true
excludeCHO = false
showboth = !excludeCHON && !excludeCHO

# always exclude C10H12O, C10H13O, C9H11ON, C6H12O4, C10H12O2, when looking at 12.5. - 15.5. (contamination) !!!
Indices2mask = findall(x -> x in ["C10H12O", "C10H13O", "C9H11ON", "C6H12O4", "C10H12O2"], names)
ptr3Traces[:,Indices2mask] .= NaN

if excludeCHON == true
	indices2maskCHON = findall(x -> occursin("N", x), names)
	ptr3Traces[:,indices2maskCHON] .= NaN
elseif excludeCHO == true
	indices2maskCHO = findall(x -> !occursin("N", x), names)
	ptr3Traces[:,indices2maskCHO] .= NaN
end

Cnrs = [2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] 						#3 excluded, because of bullshit acetone

sumtrace = DatasetFunctions.nansum(ptr3Traces[:,findall(x -> x in Cnrs, composition[:,1])],2)

#summedTracesPlot
figure()
for n in Cnrs
	Cntrace = DatasetFunctions.nansum(ptr3Traces[:,findall(x -> x == n, composition[:,1])],2)
	plot(Dates.unix2datetime.(ptr3unixTimes), Cntrace, label = "C$(n)")
end
legend()
if excludeCHON
	title("sum over Traces from CHO peaks only with same #C and #O >= 1")
elseif excludeCHO
	title("sum over Traces from CHON peaks only with same #C and #O >= 1")
elseif showboth
	title("sum over Traces from all peaks with same #C and #O >= 1")
end

#relativeContributionPlot
n = Cnrs[1]
global Cnrtraces = DatasetFunctions.nansum(ptr3Traces[:,findall(x -> x == Cnrs[1], composition[:,1])],2)
for n in Cnrs[2:end]
	global Cnrtraces = hcat(Cnrtraces, DatasetFunctions.nansum(ptr3Traces[:,findall(x -> x == n, composition[:,1])],2))
end
cmap = get_cmap("tab20")
colArr = cmap(0:18)
figure()
stackplot(Dates.unix2datetime.(ptr3unixTimes), permutedims(Cnrtraces./sumtrace), colors = colArr)
legend(Cnrs)
title("rel. Contribution of summed peaks with different #C's")
if excludeCHON
	title("rel. Contribution of summed CHO peaks only with different #C's")
elseif excludeCHO
	title("rel. Contribution of summed CHON peaks only with different #C's")
elseif showboth
	title("rel. Contribution of all summed peaks with different #C's")
end

#comparing City and Rural Barcharts
#11. - 14.
#=
avDTimes = [DateTime(2018,5,11,12) DateTime(2018,5,11,14);	# city
		DateTime(2018,5,12,12) DateTime(2018,5,12,14);	# city
		DateTime(2018,5,13,12) DateTime(2018,5,13,14);	# rural
		DateTime(2018,5,14,12) DateTime(2018,5,14,14)]	# rural
=#

avDTimes = [DateTime(2018,5,16,12) DateTime(2018,5,16,14);	# city
		DateTime(2018,5,19,12) DateTime(2018,5,19,14);	# city
		DateTime(2018,5,17,12) DateTime(2018,5,17,14);	# rural
		DateTime(2018,5,18,12) DateTime(2018,5,18,14)]	# rural

firstday = minimum(Dates.day.(avDTimes))
lastday = maximum(Dates.day.(avDTimes))
citydays = Dates.day.(avDTimes)[1:2]
ruraldays = Dates.day.(avDTimes)[3:4]

avTimes = Dates.datetime2unix.(avDTimes)

n = 1
global meanCnrContrRel = DatasetFunctions.nanmean((Cnrtraces./sumtrace)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1)
for n in 2:size(avTimes)[1]
	global meanCnrContrRel = vcat(meanCnrContrRel, DatasetFunctions.nanmean((Cnrtraces./sumtrace)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1))
end
n = 1
global meanCnrContrAbs = DatasetFunctions.nanmean((Cnrtraces)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1)
for n in 2:size(avTimes)[1]
	global meanCnrContrAbs = vcat(meanCnrContrAbs, DatasetFunctions.nanmean((Cnrtraces)[findall(x -> avTimes[n,1] <= x <= avTimes[n,2], ptr3unixTimes),:],1))
end

# absolute
figure()
global x = collect(1:size(Cnrs)[1])  # the label locations
width = 0.35  # the width of the bars
bar(x .- width*0.5, (meanCnrContrAbs[1,:] .+ meanCnrContrAbs[2,:])./2, width, label = "mean of $(citydays[1]). and $(citydays[2]).")
bar(x .+ width*0.5, (meanCnrContrAbs[3,:] .+ meanCnrContrAbs[4,:])./2, width, label = "mean of $(ruraldays[1]). and $(ruraldays[2]).")
yscale("log")
ylabel("conc. [ppt]")
xlabel("#C")
xticks(x, Cnrs)
legend()
if excludeCHON
	title("summed peaks by Cnr (only CHO)")
	savefig("$(fp)_CnrShifts_rural_city_abs_$(firstday)-$(lastday)May_CHOonly.png")
elseif excludeCHO
	title("summed peaks by Cnr (only CHON)")
	savefig("$(fp)_CnrShifts_rural_city_abs_$(firstday)-$(lastday)May_CHONonly.png")
elseif showboth
	title("summed peaks by Cnr (all)")
	savefig("$(fp)_CnrShifts_rural_city_abs_$(firstday)-$(lastday)May.png")
end

#relative
figure()
global x = collect(1:size(Cnrs)[1])  # the label locations
width = 0.7  # the width of the bars
bar(x, 100 .* ((meanCnrContrRel[1,:] .+ meanCnrContrRel[2,:])./2 .- (meanCnrContrRel[3,:] .+ meanCnrContrRel[4,:])./2), width, label = "\"city\" - \"rural\", $(firstday).-$(lastday).")
yscale("linear")
ylabel("change of rel. contribution [%]")
xticks(x, Cnrs)
#xticklabels(Cnrs)
hlines(0,1,20)
legend()
if excludeCHON
	title("relative contribution of summed peaks by Cnr (only CHO)")
	savefig("$(fp)_CnrShifts_rural_city_rel_$(firstday)-$(lastday)May_CHOonly.png")
elseif excludeCHO
	title("relative contribution of summed peaks by Cnr (only CHON)")
	savefig("$(fp)_CnrShifts_rural_city_rel_$(firstday)-$(lastday)May_CHONonly.png")
elseif showboth
	title("relative contribution of summed peaks by Cnr (all)")
	savefig("$(fp)_CnrShifts_rural_city_rel_$(firstday)-$(lastday)May.png")
end

figure()
global x = collect(1:size(Cnrs)[1])  # the label locations
width = 0.7  # the width of the bars
bar(x, 100 .* ((meanCnrContrAbs[1,:] .+ meanCnrContrAbs[2,:])./2 .- (meanCnrContrAbs[3,:] .+ meanCnrContrAbs[4,:])./2) ./ ((meanCnrContrAbs[1,:] .+ meanCnrContrAbs[2,:])./2), width, label = "\"city\" - \"rural\", $(firstday).-$(lastday).")
yscale("linear")
ylabel("rel. change per Cnr [%]")
xticks(x, Cnrs)
#xticklabels(Cnrs)
hlines(0,1,20)
legend()
if excludeCHON
	title("relative change of summed peaks by Cnr (only CHO)")
	savefig("$(fp)_CnrShifts_rural_city_relChange_eachCnr_$(firstday)-$(lastday)May_CHOonly.png")
elseif excludeCHO
	title("relative change of summed peaks by Cnr (only CHON)")
	savefig("$(fp)_CnrShifts_rural_city_relChange_eachCnr_$(firstday)-$(lastday)May_CHONonly.png")
elseif showboth
	title("relative change of summed peaks by Cnr (all)")
	savefig("$(fp)_CnrShifts_rural_city_relChange_eachCnr_$(firstday)-$(lastday)May.png")
end

# further Informations: #########################

#nonNptr3Traces = ptr3Traces[:,findall(x -> x == 0, composition[:,5])]
#O2C = ptr3Peaks[findall(x -> x == 0, composition[:,5]),8] ./ ptr3Peaks[findall(x -> x == 0, composition[:,5]),3]
#meanO2C = DatasetFunctions.nansum(nonNptr3Traces .* transpose(O2C),2) ./ DatasetFunctions.nansum(nonNptr3Traces)
#C300K = 10 .^ ((25 .- ptr3Peaks[findall(x -> x == 0, composition[:,5]),3]) .* 0.475 .- (1.7 .* ptr3Peaks[findall(x -> x == 0, composition[:,5]),8]))
#meanlgC300K = log10.( DatasetFunctions.nansum(nonNptr3Traces .* transpose(C300K),2) ./ DatasetFunctions.nansum(nonNptr3Traces) )


text:
using NLopt

import MathProgBase
nlp = joinpath(dirname(dirname(pathof(MathProgBase))), "test", "nlp.jl")
include(nlp)

nlptest(NLoptSolver(algorithm=:LD_SLSQP))
# test derivative-free
rosenbrocktest(NLoptSolver(algorithm=:LN_PRAXIS))


text:
module RangeHelpersTests
using RangeHelpers: range
using RangeHelpers
const RH = RangeHelpers
using Test

using Test

@testset "samegrid" begin
    @test samegrid(1:2, 1:2)
    @test samegrid(1:2, 0:-1:-10)
    @test samegrid(Base.OneTo(10), 11.0:1.0:20.0)

    @test !samegrid(1:2, 1.1:2.1)
    @test !samegrid(1:2, 1.1:2.1, rtol=0.02)
    @test samegrid(1:2, 1.1:2.1, rtol=0.101)
    @test samegrid(1:2, 1.1:2.1, rtol=0.101)

    @test samegrid(1:0, 1:0)
    @test samegrid(1:0, 1:10)
    @test !samegrid(1:0, 1.1:10)
    @test samegrid(1:0, 1.1:10, rtol=0.2)
    @test !samegrid(1:0, 1.0:1.1:10)
    @test samegrid(1:0, 1.0:1.1:10, rtol=0.2)
end

@testset "subdivide" begin
    @test subdivide(1.0:2:9.0, 2) == 1:1:9
    @test subdivide(1.0:3, 10) === 1.0:0.1:3.0

    @inferred subdivide(1:3, 2)
    @inferred subdivide(1:3, 2, mode=:centers)
    @test_throws ArgumentError subdivide(1:3, 2, mode=:nonsense)
    @test_throws ArgumentError subdivide(1:3, 0)

    for _ in 1:100
        r = RH.range(start=100*randn(), step=randn(), length=rand(2:100))
        @test subdivide(subdivide(r, 2), 2) == subdivide(r, 4)
        @test subdivide(subdivide(r, 2,mode=:walls), 3,mode=:walls) == subdivide(r, 2*3, mode=:walls)

        @test subdivide(subdivide(r, 2,mode=:centers), 3,mode=:centers) ≈
            subdivide(r, 2*3, mode=:centers)
    end
end

@testset "symrange" begin
    r = @inferred symrange(step=1, length=2)
    @test r === -0.5:1.0:0.5
    r = @inferred symrange(step=1, length=3)
    @test r === -1.0:1.0:1.0

    r = @inferred symrange(center=10f0, step=2, length=3)
    @test r === 8f0:2f0:12f0

    r = symrange(start=around(-4), step=2)
    @test r === -4.0:2.0:4.0

    r = @inferred symrange(start=around(-4.1), step=2, center=1)
    @test r === -4.0:2.0:6.0

    r = @inferred symrange(stop=around(4.1), step=2, center=1)
    @test r === -2.0:2.0:4.0

    r = @inferred symrange(start=above(2), step=-1)
    @test r === 2.0:-1.0:-2.0

    r = @inferred symrange(stop=strictabove(-2), step=-1)
    @test r === 1.5:-1.0:-1.5

    r = @inferred symrange(step=2.5f0, length=106)
    @test eltype(r) === Float32
end

@testset "prolong" begin
    @testset "pre post" begin
        @inferred prolong(1:2, pre=1)
        @inferred prolong(1:2, pre=-1)
        @inferred prolong(1:2, post=2)
        @inferred prolong(1:2, post=-1)
        @inferred prolong(1:1:2, post=-1)

        @test prolong(1:1:2, pre=1) == 0:1:2
        @test prolong(1:1:2, pre=1) === 0:1:2
        @test prolong(1:2, pre=1) == 0:2
        @test prolong(1:2, pre=1) === 0:2
        @test prolong(1:2, pre=-1) == 2:2
        @test prolong(1:2, pre=-1) == 2:2
        @test prolong(1:2, pre=-1) === 2:2
        @test prolong(1:2, post=2) === 1:4
        @test prolong(1:2, post=-1) === 1:1
        @test prolong(1:3, pre=1, post=-2) === 0:1

        @test prolong(1:1, pre=1, post=-1) === 0:0
        @test prolong(1:1, pre=-1, post=1) == 2:2

        @test prolong(0:0, start=around(1), stop=around(2)) == 1:2
        @test prolong(0:0, start=around(1), stop=around(2)) === 1:2

        @test prolong(0:0, start=strictabove(10), stop=strictbelow(15)) == 11:14
        @test prolong(1e9:1e10, start=strictabove(10), stop=strictbelow(15)) == 11:14
        @test prolong(1e9:1e10, start=below(10), stop=strictbelow(15)) == 10:14
        @test prolong(1e9:1e10, start=below(10), stop=around(15.1)) == 10:15

        @test prolong(Base.OneTo(3), post=2) == Base.OneTo(5)
        @test prolong(Base.OneTo(3), post=2) === Base.OneTo(5)
        @test prolong(Base.OneTo(3), stop=around(2)) === Base.OneTo(2)
    end

    @testset "anchorrange" begin
        @inferred anchorrange(15.5, start=above(11), step=2, stop=below(15))
        @test anchorrange(15.5, start=above(11), step=2, stop=below(15)) === 11.5:2.0:13.5
    end

    @testset "start" begin
        @inferred prolong(1:10, start=around(5))
        @test prolong(1:10, start=around(5)) == 5:10
        @test prolong(1:10, start=around(5)) === 5:10

        @test prolong(1:10, start=below(-2)) == -2:10
        @test prolong(1:10, start=below(-2)) === -2:10
    end

    @testset "stop" begin
        @inferred prolong(1:10, stop=around(5))
        @test prolong(1:10, stop=around(5)) == 1:5
        @test prolong(1:10.0, stop=around(5)) === 1:5.0
        @test prolong(1:10, stop=around(5)) == 1:5
    end

end

@testset "Explict" begin
    @inferred range(1,2,length=2)
    @inferred range(start=1,stop=2,length=2)
    @inferred range(start=1,stop=2,step=0.1)
    @inferred range(stop=2,step=0.1, length=3)
    @inferred range(start=2,step=0.1, length=3)

    @testset "no directions" begin
        @test range(1, 2, length = 2) === 1.0:1.0:2.0
        @test range(1, 2, step = 0.5) == [1, 1.5, 2]
        @test range(1, stop = 2, step = 0.5) === 1.0:0.5:2.0
        @test range(start=1, stop = 2, step = 0.5) === 1.0:0.5:2.0
    end

    @testset "start" begin
        @test range(strictbelow(1), 3.1, step = 2) ≈ -0.9:2.0:3.1
        @test range(below(1), 3.1, step = 2) ≈ -0.9:2.0:3.1
        @test range(around(1), 3.1, step = 2)                     === 1.1:2.0:3.1
        @test range(above(1), 3.1, step = 2)                      === 1.1:2.0:3.1
        @test range(strictabove(1), 3.1, step = 2)                === 1.1:2.0:3.1

        @test range(strictbelow(2), 0, step=-1)                   === 1:-1:0
        @test range(below(2), 0, step=-1)                         === 2:-1:0
        @test range(around(2), 0, step=-1)                        === 2:-1:0
        @test range(above(2), 0, step=-1)                         === 2:-1:0
        @test range(strictabove(2), 0, step=-1)                   === 3:-1:0

        @test range(strictbelow(2.1), 0, step=-1)                 === 2:-1:0
        @test range(below(2.1), 0, step=-1)                       === 2:-1:0
        @test range(around(2.1), 0, step=-1)                      === 2:-1:0
        @test range(above(2.1), 0, step=-1)                       === 3:-1:0
        @test range(strictabove(2.1), 0, step=-1)                 === 3:-1:0
    end

    @testset "stop" begin
        @test range(1, strictbelow(3.1), step = 2) === 1:2:3
        @test range(1, below(3.1), step = 2)       === 1:2:3
        @test range(1, around(3.1), step = 2)      === 1:2:3
        @test range(1, above(3.1), step = 2)       === 1:2:5
        @test range(1, strictabove(3.1), step = 2) === 1:2:5

        @test range(2, strictbelow(0), step=-1) === 2:-1:-1
        @test range(2, below(0), step=-1) === 2:-1:0
        @test range(2, around(0), step=-1) === 2:-1:0
        @test range(2, above(0), step=-1) === 2:-1:0
        @test range(2, strictabove(0), step=-1) === 2:-1:1

        @test range(2.1, strictbelow(0), step=-1) === 2.1:-1.0:-0.9
        @test range(2.1, below(0), step=-1) === 2.1:-1.0:-0.9
        @test range(2.1, around(0), step=-1) === 2.1:-1.0:0.1
        @test range(2.1, above(0), step=-1) === 2.1:-1.0:0.1
        @test range(2.1, strictabove(0), step=-1) === 2.1:-1.0:0.1
    end

    @testset "step" begin
        @test range(1, 3.1, step = strictbelow(2)) === 1.0:1.05:3.1
        @test range(1, 3.1, step = below(2))       === 1.0:1.05:3.1
        @test range(1, 3.1, step = around(2))      === 1.0:2.1:3.1
        @test range(1, 3.1, step = above(2))       === 1.0:2.1:3.1
        @test range(1, 3.1, step = strictabove(2)) === 1.0:2.1:3.1

        @test range(1, 3, step=strictbelow(1)) === range(1, 3, length=4)
        @test range(1, 3, step=below(1)) === 1.0:1.0:3.0
        @test range(1, 3, step=around(1)) === 1.0:1.0:3.0
        @test range(1, 3, step=above(1)) === 1.0:1.0:3.0
        @test range(1, 3, step=strictabove(1)) === 1.0:2.0:3.0

        @test_throws ArgumentError range(1,2,step=strictabove(1))

        @test range(2, 0, step=strictbelow(-1)) === 2.0:-2.0:0.0
        @test range(2, 0, step=below(-1)) === 2.0:-1.0:0.0
        @test range(2, 0, step=around(-1)) === 2.0:-1.0:0.0
        @test range(2, 0, step=above(-1)) === 2.0:-1.0:0.0
        @test range(2, 0, step=strictabove(-1)) === range(2, stop=0, length=4)
    end
end

@testset "inferred" begin
    @inferred range(strictbelow(1), 3.1, step = 2)
    @inferred range(1, above(3.1), step = 2)
    @inferred range(1, 3.1, step = around(2))
end

# useful for creating test cases
#for dir in instances(RangeHelpers.Direction)
#    start = 2.1
#    stop = dir(0)
#    step = -1
#    r = range(start, stop, step=step)
#    println("@test range($start, $stop, step=$step) === $r")
#end

@testset "asrange" begin
    @inferred asrange([1,2,3])
    @test asrange([1,2,3]) isa AbstractRange
    @test asrange([1,2,3]) === 1:1:3
    @test asrange(1:3) === 1:3
    @test_throws ArgumentError asrange([1,2,4])
    @test asrange([1]) == range(1,stop=1,length=1)
end

@testset "binwalls, bincenters" begin
    for _ in 1:100
        start = 10randn()
        stop=10randn()
        len=rand(3:100)
        r = range(start, stop=stop, length=len)
        w = RH.binwalls(r)
        @test r ≈ RH.bincenters(RH.binwalls(r))
        @test r ≈ RH.binwalls(RH.bincenters(r))
    end
    @test RH.binwalls(1:4) == [0.5, 1.5, 2.5, 3.5, 4.5]
    @test RH.binwalls(1:4, first=false, last=true) == [1.5,2.5,3.5, 4.5]
    @test RH.bincenters(1:2:9) == [2,4,6,8]
    @inferred RH.bincenters(1:3)
    @inferred RH.binwalls(1:3)
end

@testset "searchsortedat" begin
    @test RH.searchsortedat(10:10:30, RH.around(-100)) === 1
    @test RH.searchsortedat(10:10:30, RH.around(5)) === 1
    @test RH.searchsortedat(10:10:30, RH.around(10)) === 1
    @test RH.searchsortedat(10:10:30, RH.around(14)) === 1
    @test RH.searchsortedat(10:10:30, RH.around(16.0)) === 2
    @test RH.searchsortedat(10:10:30, RH.around(20.0)) === 2
    @test RH.searchsortedat(10:10:30, RH.around(24.0)) === 2
    @test RH.searchsortedat(10:10:30, RH.around(26.0)) === 3
    @test RH.searchsortedat(10:10:30, RH.around(30)) === 3
    @test RH.searchsortedat(10:10:30, RH.around(50)) === 3
    if VERSION >= v"1.6"
        @test RH.searchsortedat(10:10:30, RH.around(-Inf)) === 1
        @test RH.searchsortedat(10:10:30, RH.around(Inf)) === 3
    end

    for _ in 1:100
        r = sort!(randn(rand(2:10)))
        x = randn()
        isbe = RH.searchsortedat(r, RH.strictbelow(x))
        ibe  = RH.searchsortedat(r, RH.below(x))
        iar  = RH.searchsortedat(r, RH.around(x))
        iab  = RH.searchsortedat(r, RH.above(x))
        isab = RH.searchsortedat(r, RH.strictabove(x))
        @test iar-1 <= isbe == ibe <= iar <= iab == isab <= isbe + 1
        @test_throws ArgumentError RH.searchsortedat(r, x)
        @test iar in eachindex(r)
        if ibe < firstindex(r)
            @test x < first(r)
        else
            @test r[ibe] < x
        end
        if iab > lastindex(r)
            @test x > last(r)
        else
            @test r[iab] > x
        end

        i = rand(eachindex(r))
        x = r[i]
        @test RH.searchsortedat(r, RH.strictbelow(x)) == i-1
        @test RH.searchsortedat(r, RH.below(x)) == i
        @test RH.searchsortedat(r, RH.around(x)) == i
        @test RH.searchsortedat(r, RH.above(x)) == i
        @test RH.searchsortedat(r, RH.strictabove(x)) == i+1
        @test RH.searchsortedat(r, x) == i
    end
    @inferred RH.searchsortedat(10:10:30, RH.around(4))
    @inferred RH.searchsortedat(10:10:30, RH.around(4.4))
    @inferred RH.searchsortedat(collect(10:10:30), RH.around(4.4))

    # double points
    coll = [1,1,1]
    @test RH.searchsortedat(coll, strictabove(1)) === 4
    @test RH.searchsortedat(coll, strictabove(0)) === 1
    @test RH.searchsortedat(coll, strictbelow(1)) === 0
    @test RH.searchsortedat(coll, strictbelow(2)) === 3
    coll = [0,1,1,1,2]
    @test RH.searchsortedat(coll, strictabove(1)) === 5
    @test RH.searchsortedat(coll, strictbelow(1)) === 1
end

@testset "doctest" begin
    import Documenter
    Documenter.doctest(RangeHelpers)
end

end#module


text:
# Turing.jl version of model at https://github.com/stan-dev/blob/master/basic_estimators/normal_loc.stan

@model nlmodel(y) = begin
  mu ~ Uniform(-10, 10)
  for n = 1:5
    y[n] ~ Normal(mu, 1.0)
  end
  mu
end


text:
@testset "Console Function Tests" begin

  @test isdefined(Julz, :console) == true

end


text:
# system/z-list,jl

include("Afft.jl")
include("linear-map/block.jl")
include("linear-map/getindex.jl")
include("linear-map/hcat.jl")


text:
# https://gridap.github.io/Tutorials/stable/pages/t009_stokes/

using Pkg
Pkg.activate("simulate")

using Gridap

n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain, partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])

V = TestFESpace(
  reffe=:QLagrangian, conformity=:H1, valuetype=VectorValue{2,Float64},
  model=model, labels=labels, order=2, dirichlet_tags=["diri0","diri1"])
Q = TestFESpace(
  reffe=:PLagrangian, conformity=:L2, valuetype=Float64,
  model=model, order=1, constraint=:zeromean)
Y = MultiFieldFESpace([V,Q])

u0 = VectorValue(0,0)
u1 = VectorValue(1,0)
U = TrialFESpace(V,[u0,u1])
P = TrialFESpace(Q)
X = MultiFieldFESpace([U,P])

trian = get_triangulation(model); degree = 2
quad = CellQuadrature(trian,degree)

function a(x,y)
  v,q = y
  u,p = x
  ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)
end
t_Ω = LinearFETerm(a,trian,quad)
op = AffineFEOperator(X,Y,t_Ω)
uh, ph = solve(op)

writevtk(trian,"results",order=2,cellfields=["uh"=>uh,"ph"=>ph])

max_stars_repo_name: stevengogogo/ProgrammableAging

max_issues_repo_name: stevengogogo/ProgrammableAging

max_forks_repo_name: stevengogogo/ProgrammableAging


text:
using MatrixDepot
using SparseArrays
using MAT
using Random
rand(0)

accepted_problems = ["least squares problem"]

# List all problems (Filter some)
list = mdlist("*/*")

# Get matrix A and vector b for all problems
prob_count = 0
md_info = 0
for problem in list
    global md_kind = mdinfo(problem).content[4].items[end][1].content[1]  # Nasty

    # Check if type is correct
    if split(md_kind)[1] == "kind:"  # Got the right kind
        kind = join(split(md_kind)[2:end], " ")
        if kind in accepted_problems
            md = mdopen(problem)
            MatrixDepot.addmetadata!(md.data)
            print("Name = $(md.data.name)\n")
            A = float(md.A)
            if issymmetric(A)
                print("Making it a generic matrix\n")
                A = SparseMatrixCSC(A)
            end
            global prob_count += 1
            print("Found $(kind) n = $(prob_count)\n")
            (m, n) = size(A)
            new_problem =  Dict("A" => A, "name" => md.data.name)
            try
                new_problem["b"] = float(md.b)
                print("Storing also b\n")
            catch e
                print("No field b. Creating random one.\n")
                s0 = randn(m)
                x0 = randn(n)
                b = A * x0 + s0
                new_problem["b"] = b
            end
            file_name = replace("$(md.data.name).mat", "/" => "_")
            try
                matwrite(file_name, new_problem)
            catch e
                print(e)
                break
            end

            print("Written to file $(file_name)\n")
        end
    else
        print("Wrong type extracted from metadata. Type = $(md_kind)\n")
    end

end

print("Total number of problems = $(prob_count)\n")


text:
"""
    ExactConvolution

Implement efficient exact integer convolution using BigInt's multiplication.
"""
module ExactConvolution

export exact_conv

include("convolution.jl")

end


text:
import Compose
import .Weave

Base.showable(m::MIME"application/svg", ctx::Context) = true
Base.showable(m::MIME"application/png", ctx::Context) = true

function Base.display(report::Weave.Report, ::MIME"image/svg+xml", ctx::Context)
    chunk = report.cur_chunk

    w = chunk.options[:fig_width]Compose.inch
    h = chunk.options[:fig_height]Compose.inch
    format = chunk.options[:fig_ext]
    dpi = chunk.options[:dpi]

    full_name, rel_name = Weave.get_figname(report, chunk, ext = format)

    push!(report.figures, rel_name)
    report.fignum += 1
    if format == ".svg"
        Compose.draw(Compose.SVG(full_name, w, h), ctx)
    elseif format == ".png"
        Compose.draw(Compose.PNG(full_name, w, h, dpi = dpi), ctx)
    elseif format == ".pdf"
        Compose.draw(Compose.PDF(full_name, w, h), ctx)
    elseif format == ".ps"
        Compose.draw(Compose.PS(full_name, w, h), ctx)
    elseif format == ".tex"
        Compose.draw(Compose.PGF(full_name, w, h, true), ctx)
    else
        @warn("Can't save figure. Unsupported format, $format")
    end
end

"""
    kill_weave_temporaries(path)

Kills temporary directories and files associated with weave...
"""
function kill_weave_temporaries(path::AbstractString)
    for fn in filter(f->!isnothing(match(r"^jl_[a-zA-Z0-9]+$", f)), readdir(path))
        rm(joinpath(path,fn);force=true, recursive=true)
    end
end

@info "Weave support loaded into NeXLSpectrum"

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
import Term: reshape_text
import Term.Consoles: clear

strings = [
    # latin alphabet
    "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua.",
    "Lorem {red}ipsum dolor sit {underline}amet, consectetur{/underline} adipiscing elit, {/red}{blue}sed do eiusmod tempor incididunt{/blue} ut labore et dolore magna aliqua.",
    "Lorem{red}ipsumdolorsit{underline}amet, consectetur{/underline} adipiscing elit, {/red}seddoeiusmo{blue}dtemporincididunt{/blue}ut labore et dolore magna aliqua.",

    # thai
    "ต้าอ่วยวาทกรรมอาว์เซี้ยว กระดี๊กระด๊า ช็อปซาดิสต์โมจิดีพาร์ตเมนต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเตอร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",
    "ต้าอ่วยวาท{red}กรรมอาว์เซี้ยว กระดี๊กระด๊า {/red}ช็อปซาดิสต์โมจิดีพาร์ตเม{blue underline}นต์ อินดอร์วิว สี่แยกมาร์กจ๊อกกี้ โซนี่บัตเต{/blue underline}อร์ฮันนีมูน ยาวีแพลนหงวนสคริปต์ แจ็กพ็อตต่อรองโทรโข่งยากูซ่ารุมบ้า บอมบ์เบอร์รีวีเจดีพาร์ทเมนท์ บอยคอตต์เฟอร์รี่บึมมาราธอน ",

    # # korean
    "국가유공자·상이군경 및 전몰군경의 유가족은 법률이 정하는 바에 의하여",
    "국[red]가유공자·상이군{bold}경 및 전{/bold}몰군경의 유{/red}가족은 법률이 정하는 바에 의하여",

    # # large chars
    "┌────────────────┬────────────────┬────────────────┬────────────────┬──────────────",
    "┌────────────abcde────┬──────────── ────┬────────abcde fghi────────┬────────────────┬──────────────",
    "┌─────────{red}───ab{/red}cde────┬──────{green}────── ────┬────────abcde fghi{/green}────────┬────────────────┬──────────────",
    "┌──────────{red}────{/red}──┬{blue bold}────────────────┬{/blue bold}──────────────┬────────────────┬──────────────end",

    # nospaces
    "."^100,
    ".{red}|||{/red}...."^10,
    ".|||...."^10,
    replace(
        """
Lorem ipsum {bold}dolor sit{/bold} amet, consectetur adipiscing elit,
ed do e{red}iusmod tempor incididunt{/red} ut {bold}labore et {underline}dolore{/underline} magna aliqua.{/bold} Ut enim ad minim
veniam, quis{green} nostrud exercitation {on_black}ullamco laboris nisi ut aliquip ex {/on_black}
ea commodo consequat.{blue} Duis aute irure dolor in{/blue} reprehenderit
in voluptate velit{/green} esse {italic}cillum dolore{/italic}{red} eu{/red}{italic green} fugiat {/italic green}nulla
pariatur. Excepteur{red} sint{/red}{blue} occaecat cupidatat {/blue}non proident,
sunt in culpa qui {italic}officia{/italic} deserunt mollit anim
id est laborum.""",
        "\n" => "",
    ),
]

print("\n"^3)
width = 33
for text in strings
    # println('_'^width)
    # println(reshape_text(text, width))
    print("\n\n")
    println('_'^width)
    tprint(reshape_text(text, width))
    println('_'^width)
    tprint(text)
end

# text = strings[end-1]
# println('_'^width)
# println(reshape_text(text, width))