text:
using Base.Iterators
using SimpleMatch
using StructEquality
using ProxyInterfaces

"""
All ExprParserWithParsed have a common parse_expr method, namely that all struct fields are given directly as keyword arguments.
"""
function parse_expr(parser::ExprParserWithParsed; kw...)
  names = keys(kw)  # kw isa Base.Iterators.Pairs
  matches = [parse_expr(getproperty(parser, name), value) for (name, value) in kw]
  kw′ = NamedTuple{names}(matches)
  ExprParsed(typeof(parser))(;kw′...)
end
function parse_expr(parser::ExprParserWithParsed, any)
  throw(ParseError("$(typeof(parser)) has no `parse_expr` method defined to capture Type `$(typeof(any))`. Got: `$any`."))
end

# Overwriting show to print content by fields instead of positions.

# one-line version
function Base.show(io::IO, obj::Union{ExprParserWithParsed, ExprParsed})
  type_name = typeof(obj).name.name
  print(io, "EP.$type_name(")
  field_names = fieldnames(typeof(obj))
  for field in field_names[1:end-1]
    print(io, "$field=$(repr(getproperty(obj, field))), ")
  end
  # last field extra
  field = field_names[end]
  print(io, "$field=$(repr(getproperty(obj, field)))")
  print(io, ")")
end

# Multiline version, following https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing-1
function Base.show(io::IO, ::MIME"text/plain", obj::Union{ExprParserWithParsed, ExprParsed})
  type_name = typeof(obj).name.name
  println(io, "EP.$type_name(")
  field_lengths = length.(string.(fieldnames(typeof(obj))))
  max_field_length = maximum(field_lengths)
  for (field, field_length) in zip(fieldnames(typeof(obj)), field_lengths)
    filler = join(fill(" ", max_field_length - field_length))
    println(io, "  $field $filler= $(repr(getproperty(obj, field)))")
  end
  println(io, ")")
end

# ExprParsers
# ===========

# Expr
# ----

"""
    EP.Expr(head = EP.anything, args = EP.anything)
    EP.Expr(expr; [ignore_linenumbernodes=true])

It is the most flexible parser, but hence also the least plug-and-play.

Parses the following
```julia
Base.Expr(head, args...)
```

# Examples

```julia
julia> using ExprParsers

julia> parser = EP.Expr(head = :vect);

julia> parse_expr(parser, :([1,2,3]))
EP.Expr_Parsed(
  head = :vect
  args = Any[1, 2, 3]
)
julia> parse_expr(parser, :(f(a) = a))
ERROR: ParseError: Using default `==` comparison, but parser `:vect` ≠ value `:(=)`.
```

Also see [`Indexed`](@ref) for an example to combine `EP.Indexed` with `EP.Expr`.
"""
@exprparser struct Expr
  head = anything
  args = anything = []
end
function Expr(expr::Base.Expr; ignore_linenumbernodes=true)
  if ignore_linenumbernodes
    Expr(expr.head, Block(expr.args..., ignore_linenumbernodes=true))
  else
    Expr(expr.head, expr.args)
  end
end
parse_expr(parser::Expr, expr::Base.Expr) = parse_expr(parser, head = expr.head, args = expr.args)
function to_expr(parsed::Expr_Parsed)
  if parsed.args isa Block_Parsed
    Base.Expr(to_expr(parsed.head), to_expr(parsed.args.exprs)...)
  else
    Base.Expr(to_expr(parsed.head), to_expr(parsed.args)...)
  end
end

# Block
# -----

@doc raw"""
    EP.Block(block_expr; [ignore_linenumbernodes = true])
    EP.Block(expr1, expr2, ...; [ignore_linenumbernodes = true])
    EP.Block()

Helper to parse blocks of code (i.e. `expr.head == :block`) or a given list of expr respectively.

The main purpose is to handle linenumbernodes, otherwise it behaves similar to plain Vector of Expr.

Parses the following
```julia
quote
  any
  4
end
[:(a = 4), 42, :anyvector]
(:(a = 4), 42, :or_tuple_of_expr)
end
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Block(quote
         $(EP.anything)
         $(EP.anysymbol)
         13
       end);

julia> parse_expr(parser, [:(a = 4), :hi, 13])
EP.Block_Parsed(
  exprs = Any[:(a = 4), :hi, 13]
)
julia> parse_expr(parser, quote
         whatever(a) = a
         asymbol
         14
       end)
ERROR: ParseError: Using default `==` comparison, but parser `13` ≠ value `14`.
```
Used within [`EP.Expr`](@ref).
"""
@def_structequal struct Block <: ExprParserWithParsed
  exprs
  ignore_linenumbernodes::Bool

  function Block(expr::Base.Expr; ignore_linenumbernodes = true)
    @assert expr.head == :block "Trying to construct a Block parser from a $(expr.head). Expecting :block."
    new(expr.args, ignore_linenumbernodes)
  end
  Block(exprs...; ignore_linenumbernodes=true) = new(exprs, ignore_linenumbernodes)
  Block() = new(Iterator(repeated(anything)), true)
end
@def_structequal mutable struct Block_Parsed <: ExprParsed
  exprs
end
ExprParsed(::Base.Type{Block}) = Block_Parsed
to_expr(parsed::Block_Parsed) = Base.Expr(:block, map(to_expr, parsed.exprs)...)

function parse_expr(parser::Block, expr::Base.Expr)
  @passert expr.head == :block
  parse_expr(parser, expr.args)
end
function parse_expr(parser::Block, exprs::Union{Vector, Tuple})
  parsed_expr = if parser.ignore_linenumbernodes
    parse_expr_ignore_LineNumberNodes(parser.exprs, exprs)
  else
    parse_expr(parser.exprs, exprs)
  end
  Block_Parsed(parsed_expr)
end

parse_expr_ignore_LineNumberNodes(parser, expr) = parse_expr(parser, expr)
function parse_expr_ignore_LineNumberNodes(parsers::Union{Vector, Tuple, Iterator}, exprs::Tuple)
  parse_expr_ignore_LineNumberNodes(parsers, collect(exprs))
end
function parse_expr_ignore_LineNumberNodes(parsers::Union{Vector, Tuple, Iterator}, exprs::Vector)
  exprs_parsed = similar(exprs)
  parsers_filtered = Iterators.filter(x -> !isa(x, LineNumberNode), parsers)
  next_parser = iterate(parsers_filtered)
  for (i, expr) in enumerate(exprs)
    if expr isa LineNumberNode
      exprs_parsed[i] = expr
    else
      @passert issomething(next_parser) "Found expr `$expr` in exprs $exprs without a match in the parsers $parsers."
      parser, state = next_parser
      exprs_parsed[i] = parse_expr_ignore_LineNumberNodes(parser, expr)
      next_parser = iterate(parsers_filtered, state)
    end
  end
  if !isa(Base.IteratorSize(parsers), Base.IsInfinite)
    @passert isnothing(next_parser) "While exprs $exprs is exhausted, there is still a left over parser $(next_parser[1]) in the parsers $parsers."
  end
  exprs_parsed
end
function parse_expr_ignore_LineNumberNodes(parser::Expr, expr::Expr)
  @passert parser.head == expr.head
  Base.Expr(parser.head, parse_expr_ignore_LineNumberNodes(parser.args, expr.args)...)
end

# Module
# ------
"""
    EP.Module(no_bare = EP.Isa(Bool), name = EP.anything, body = EP.anything)

Parses the following
```julia
module A
end
```
"""
@exprparser struct Module
  "indicates whether the module is a standard module (true), or a baremodule (false)"
  no_bare = Isa(Bool) = true
  name = anything
  body = anything = quote end
end

function parse_expr(parser::Module, expr::Base.Expr)
  @passert expr.head == :module
  parse_expr(parser,
    no_bare = expr.args[1],
    name = expr.args[2],
    body = expr.args[3])
end

function to_expr(parsed::Module_Parsed)
  Base.Expr(:module, parsed.no_bare, parsed.name, parsed.body)
end

# Macro
# -----

"""
    EP.Macro(name = EP.anything, args = EP.anything, linenumber = EP.Isa(LineNumberNode))

Parses the following
```julia
@macroname arg1 arg2 ...
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Macro(name = :mymacro);

julia> parse_expr(parser, :(@mymacro 1 two))
EP.Macro_Parsed(
  name       = :mymacro
  args       = Any[1, :two]
  linenumber = :(#= none:1 =#)
)
julia> parse_expr(parser, :(@anothermacro))
ERROR: ParseError: Using default `==` comparison, but parser `:mymacro` ≠ value `:anothermacro`.
```
"""
@exprparser struct Macro
  name = anything
  args = anything = []
  linenumber = Isa(LineNumberNode) = LineNumberNode(0)
end
to_expr(parsed::Macro_Parsed) = Base.Expr(:macrocall, Base.Symbol("@", parsed.name), parsed.linenumber, parsed.args...)
function parse_expr(parser::Macro, expr::Base.Expr)
  @passert expr.head == :macrocall
  name_with_at = expr.args[1]
  name_without_at = Base.Symbol(string(name_with_at)[2:end])
  parse_expr(parser, name = name_without_at, args = expr.args[3:end], linenumber = expr.args[2])
end

# Assignment
# ----------

"""
    EP.Assignment(left = EP.anything, right = EP.anything)

Parses the following
```
left = right
```

# Examples

```jldoctest
julia> using ExprParsers

julia> parser = EP.Assignment(left = EP.anysymbol)
EP.Assignment(
  left  = ExprParsers.Isa{Symbol}()
  right = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(a = [1,2,3,4]))
EP.Assignment_Parsed(
  left  = :a
  right = :([1, 2, 3, 4])
)
julia> parse_expr(parser, :(f(a) = a))
ERROR: ParseError: Expected type `Symbol`, got `f(a)` of type `Expr`.
```
"""
@exprparser struct Assignment
  left = anything
  right = anything
end
to_expr(parsed::Assignment_Parsed) = :($(parsed.left) = $(parsed.right))
function parse_expr(parser::Assignment, expr::Base.Expr)
  @passert expr.head == :(=)  "assignment head not found"
  parse_expr(parser, left = expr.args[1], right = expr.args[2])
end

# NestedDot
# ---------

"""
    EP.NestedDot(base = EP.anything, properties = EP.anything)

Parses the following
```julia
a.b
fun(T{:hi}).b.c.d.e.f
```

# Examples
```jldoctest
julia> using ExprParsers; using Base.Iterators

julia> parser = EP.NestedDot(
         properties = EP.Iterator(repeated(
           EP.SatisfiesPredicate("It should start with 'a'.") do x
             startswith(string(x), "a")
           end
         ))
       );

julia> parse_expr(parser, :(fun(T{:hi}).aone.atwo.athree))
EP.NestedDot_Parsed(
  base       = :(fun(T{:hi}))
  properties = [:aone, :atwo, :athree]
)
julia> parse_expr(parser, :(fun(T{:hi}).aone.btwo.athree))
ERROR: ParseError: Predicate `#1` returned false on expr `btwo`. It should start with 'a'.
```
"""
@exprparser struct NestedDot
  base = anything
  properties = anything
end
function parse_expr(parser::NestedDot, expr::Base.Expr)
  @passert expr.head == :.
  exprs = _collect_nested_dots(expr)
  parse_expr(parser, base = exprs[1], properties = exprs[2:end])
end
_collect_nested_dots(expr::Base.Expr) = _collect_nested_dots(expr, Val{expr.head}())
function _collect_nested_dots(expr::Base.Expr, ::Val{:.})
  @passert length(expr.args) == 2
  @passert expr.args[2] isa QuoteNode
  [_collect_nested_dots(expr.args[1]); expr.args[2].value]
end
_collect_nested_dots(expr::Base.Expr, _::Val) = [expr]
_collect_nested_dots(any) = [any]

to_expr(parsed::NestedDot_Parsed) = foldl(parsed.properties; init=parsed.base) do expr, property
  Base.Expr(:., expr, QuoteNode(property))
end

# Reference
# ---------

"""
    EP.Reference(name = EP.anything, curlies = EP.anything)

Parses the following
```julia
a
a{b, c}
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Reference(curlies = [:A, EP.anysymbol])
EP.Reference(
  name    = ExprParsers.Isa{Any}()
  curlies = Any[:A, ExprParsers.Isa{Symbol}()]
)
julia> parse_expr(parser, :(SomeType{A, B}))
EP.Reference_Parsed(
  name    = :SomeType
  curlies = [:A, :B]
)
julia> parse_expr(parser, :(SomeType{A}))
ERROR: ParseError: length(parser) == length(values) = false
  length(parser) = 2
  parser = Any[:A, ExprParsers.Isa{Symbol}()]
  length(values) = 1
  values = Any[:A]
julia> parse_expr(parser, :(SomeType{B, C}))
ERROR: ParseError: Using default `==` comparison, but parser `:A` ≠ value `:B`.
julia> parse_expr(parser, :(SomeType{A, 1}))
ERROR: ParseError: Expected type `Symbol`, got `1` of type `Int64`.
```
"""
@exprparser struct Reference
  name = anything
  curlies = anything = []
end

parse_expr(parser::Reference, expr::Base.Symbol) = parse_expr(parser, name = expr, curlies = [])
function parse_expr(parser::Reference, expr::Base.Expr)
  if expr.head == :curly
    parse_expr(parser, name = expr.args[1], curlies = expr.args[2:end])
  elseif expr.head == :.
    parse_expr(parser, name = expr, curlies = [])
  else
    throw(ParseError("Cannot parse expr `$expr` as reference: expr.head `$(expr.head)` not in `[:curly, :.]`."))
  end
end

function to_expr(parsed::Reference_Parsed)
  if isempty(parsed.curlies)
    to_expr(parsed.name)
  else
    :( $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)} )
  end
end

# Call
# ----

@doc raw"""
    EP.Call(
      name = EP.anything,
      curlies = EP.anything,
      args = EP.anything,
      kwargs = EP.anything)

Parses the following
```julia
a()
a(b, c)
a{b, c}(d, e)
a{b, c}(d, e, f = 1; g = :two)
```

Note that all keyword arguments are collected into the `kwargs` field, also those before `;`, corresponding to standard
julia call semantics.

Note that keyword arguments are represented using the default Expr representation
`Expr(:kw, :key, "value")`.

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Call(name = :myfunc)
EP.Call(
  name    = :myfunc
  curlies = ExprParsers.Isa{Any}()
  args    = ExprParsers.Isa{Any}()
  kwargs  = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(myfunc(a, b, c = 1; d = :hi)))
EP.Call_Parsed(
  name    = :myfunc
  curlies = Any[]
  args    = Any[:a, :b]
  kwargs  = Any[:($(Expr(:kw, :c, 1))), :($(Expr(:kw, :d, :(:hi))))]
)
julia> parse_expr(parser, :(anotherfunc(a, b, c = 1; d = :hi)))
ERROR: ParseError: Using default `==` comparison, but parser `:myfunc` ≠ value `:anotherfunc`.
```
"""
@exprparser struct Call
  name = anything
  curlies = anything = []
  args = anything = []
  kwargs = anything = []
end

function parse_expr(parser::Call, expr::Base.Expr)
  @passert expr.head == :call
  called = parse_expr(Reference(), expr.args[1])
  args, kwargs = _extract_args_kwargs__collect_all_kw_into_kwargs(expr.args[2:end])
  parse_expr(parser, name = called.name, curlies = called.curlies, args = args, kwargs = kwargs)
end

function _extract_args_kwargs__collect_all_kw_into_kwargs(expr_args)
  args = []
  kwargs = []

  otherargs, parameters = if isempty(expr_args)
    [], []
  elseif isexpr(expr_args[1], :parameters)
    expr_args[2:end], expr_args[1].args
  else
    expr_args, []
  end

  for p in otherargs
    if isexpr(p, :kw)
      push!(kwargs, p)
    else
      push!(args, p)
    end
  end
  append!(kwargs, parameters)
  args, kwargs
end

to_expr(parsed::Call_Parsed) = _to_expr_Call(parsed, Val{nonempty(parsed.curlies)}(), Val{nonempty(parsed.kwargs)}())
_to_expr_Call(parsed, iscurlies::True, iswheres::True) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...),;$(to_expr(parsed.kwargs)...))
)
_to_expr_Call(parsed, iscurlies::True, iswheres::False) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...),)
)
_to_expr_Call(parsed, iscurlies::False, iswheres::True) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...),;$(to_expr(parsed.kwargs)...))
)
_to_expr_Call(parsed, iscurlies::False, iswheres::False) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...),)
)

# Signature
# ---------

@doc raw"""
    EP.Signature(
      name = EP.anything,
      curlies = EP.anything,
      args = EP.anything,
      kwargs = EP.anything,
      wheres = EP.anything)

Similar to [`EP.Call`](@ref) but with a couple of differences
- extra `wheres`
- the `name` field might stay empty
- `args` can also contain `Expr(:kw, key, value)` values
  (corresponds to default value syntax, which is only available in signatures)

Parses the following:
```julia
a(b, c::Any)
a(b::B, c) where B
(::Any, c::C) where {C <: Number}
f(::Any, c::C, d::Int=1; e=true) where {C <: Number}
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Signature()
EP.Signature(
  name    = ExprParsers.Isa{Any}()
  curlies = ExprParsers.Isa{Any}()
  args    = ExprParsers.Isa{Any}()
  kwargs  = ExprParsers.Isa{Any}()
  wheres  = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(f{T}(a::T) where T))
EP.Signature_Parsed(
  name    = :f
  curlies = Any[:T]
  args    = Any[:(a::T)]
  kwargs  = Any[]
  wheres  = Any[:T]
)
julia> parse_expr(parser, :((a, b::T) where T))
EP.Signature_Parsed(
  name    = nothing
  curlies = Any[]
  args    = Any[:a, :(b::T)]
  kwargs  = Any[]
  wheres  = Any[:T]
)
julia> parse_expr(parser, :(f(a, b::T, c::Any=3; d=true) where T))
EP.Signature_Parsed(
  name    = :f
  curlies = Any[]
  args    = Any[:a, :(b::T), :($(Expr(:kw, :(c::Any), 3)))]
  kwargs  = Any[:($(Expr(:kw, :d, true)))]
  wheres  = Any[:T]
)
julia> parse_expr(parser, :(f(a) = a))
ERROR: ParseError: expr.head in (:where, :call, :tuple) = false
  expr.head = :(=)
  expr = :(f(a) = begin
          #= none:1 =#
          a
      end)
  (:where, :call, :tuple) = (:where, :call, :tuple)
```
"""
@exprparser struct Signature
  name = anything = nothing
  curlies = anything = []
  args = anything = []
  kwargs = anything = []
  wheres = anything = []
end

function parse_expr(parser::Signature, expr::Base.Expr)
  @passert expr.head in (:where, :call, :tuple)
  # we don't reuse ParsedCallExpr because functions-signatures can actually be anonymous
  call, wheres = split_where(expr)
  @passert call isa Base.Expr  # until before here we still match assignments
  @passert call.head in (:tuple, :call)
  isanonymous = call.head == :tuple

  name, curlies, args, kwargs = if isanonymous
    _args, _kwargs = _extract_args_kwargs__no_processing(call.args)
    nothing, [], _args, _kwargs
  else
    called = parse_expr(Reference(), call.args[1])
    _args, _kwargs = _extract_args_kwargs__no_processing(call.args[2:end])
    called.name, called.curlies, _args, _kwargs
  end
  parse_expr(parser, name = name, curlies = curlies, args = args, kwargs = kwargs, wheres = wheres)
end

function _extract_args_kwargs__no_processing(expr_args)
  args, kwargs = if isempty(expr_args)
    [], []
  elseif isexpr(expr_args[1], :parameters)
    expr_args[2:end], expr_args[1].args
  else
    expr_args, []
  end
  args, kwargs
end

to_expr(parsed::Signature_Parsed) = _to_expr_Signature(parsed, parsed.name, Val{nonempty(parsed.curlies)}(), Val{nonempty(parsed.kwargs)}(), Val{nonempty(parsed.wheres)}())
_to_expr_Signature(parsed, ::Nothing, iscurlies::True, ::False, ::False) = error("Impossible Reached: Given curlies but no name. parsed.curlies = $(parsed.curlies)")
_to_expr_Signature(parsed, ::Nothing, iscurlies::True, ::False, ::True) = error("Impossible Reached: Given curlies but no name. parsed.curlies = $(parsed.curlies)")
_to_expr_Signature(parsed, ::Nothing, iscurlies::True, ::True, ::False) = error("Impossible Reached: Given curlies but no name. parsed.curlies = $(parsed.curlies)")
_to_expr_Signature(parsed, ::Nothing, iscurlies::True, ::True, ::True) = error("Impossible Reached: Given curlies but no name. parsed.curlies = $(parsed.curlies)")
_to_expr_Signature(parsed, ::Nothing, ::False, iskwargs::False, iswheres::False) = :(
  ($(to_expr(parsed.args)...),)
)
_to_expr_Signature(parsed, ::Nothing, ::False, iskwargs::True, iswheres::False) = :(
  ($(to_expr(parsed.args)...), ; $(to_expr(parsed.kwargs)...))
)
_to_expr_Signature(parsed, ::Nothing, ::False, iskwargs::False, iswheres::True) = :(
  ($(to_expr(parsed.args)...),) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, ::Nothing, ::False, iskwargs::True, iswheres::True) = :(
  ($(to_expr(parsed.args)...),; $(to_expr(parsed.kwargs)...)) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, name, iscurlies::True, iskwargs::True, iswheres::True) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...), ; $(to_expr(parsed.kwargs)...)) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, name, iscurlies::True, iskwargs::True, iswheres::False) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...), ; $(to_expr(parsed.kwargs)...))
)
_to_expr_Signature(parsed, name, iscurlies::True, iskwargs::False, iswheres::True) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...),) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, name, iscurlies::True, iskwargs::False, iswheres::False) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}($(to_expr(parsed.args)...),)
)
_to_expr_Signature(parsed, name, iscurlies::False, iskwargs::True, iswheres::True) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...), ; $(to_expr(parsed.kwargs)...)) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, name, iscurlies::False, iskwargs::True, iswheres::False) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...), ; $(to_expr(parsed.kwargs)...))
)
_to_expr_Signature(parsed, name, iscurlies::False, iskwargs::False, iswheres::True) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...),) where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Signature(parsed, name, iscurlies::False, iskwargs::False, iswheres::False) = :(
  $(to_expr(parsed.name))($(to_expr(parsed.args)...),)
)

# Function
# --------

"""
    EP.Function(
      name = EP.anything,
      curlies = EP.anything,
      args = EP.anything,
      kwargs = EP.anything,
      wheres = EP.anything,
      body = EP.anything)

Parses full functions. For instance
```julia
function a(b, c) where B
  d
end
a(b, c) = d
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.Function(
         args = [EP.anything for i in 1:3],
       )
EP.Function(
  name    = ExprParsers.Isa{Any}()
  curlies = ExprParsers.Isa{Any}()
  args    = [ExprParsers.Isa{Any}(), ExprParsers.Isa{Any}(), ExprParsers.Isa{Any}()]
  kwargs  = ExprParsers.Isa{Any}()
  wheres  = ExprParsers.Isa{Any}()
  body    = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(f(a, b, c) = a + b + c))
EP.Function_Parsed(
  name    = :f
  curlies = Any[]
  args    = [:a, :b, :c]
  kwargs  = Any[]
  wheres  = Any[]
  body    = quote
    #= none:1 =#
    a + b + c
end
)
julia> parse_expr(parser, :(
         function g(a)
           a
         end
       ))
ERROR: ParseError: length(parser) == length(values) = false
  length(parser) = 3
  parser = [ExprParsers.Isa{Any}(), ExprParsers.Isa{Any}(), ExprParsers.Isa{Any}()]
  length(values) = 1
  values = Any[:a]
julia> parse_expr(parser, :a)
ERROR: ParseError: ExprParsers.Function has no `parse_expr` method defined to capture Type `Symbol`. Got: `a`.
```
"""
@exprparser struct Function
  name = anything = nothing
  curlies = anything = []
  args = anything = []
  kwargs = anything = []
  wheres = anything = []
  body = anything = nothing
end

function parse_expr(parser::Function, expr::Base.Expr)
  @passert expr.head in (:function, :(=))
  @passert length(expr.args) == 2
  body = expr.args[2]

  signature = parse_expr(Signature(), expr.args[1])
  parse_expr(parser,
    name = signature.name,
    curlies = signature.curlies,
    args = signature.args,
    kwargs = signature.kwargs,
    wheres = signature.wheres,
    body = body)
end

to_expr(parsed::Function_Parsed) = Base.Expr(:function, to_expr(Signature_Parsed(parsed)), to_expr(parsed.body))

Signature_Parsed(parsed::Function_Parsed) = Signature_Parsed(
  name = parsed.name,
  curlies = parsed.curlies,
  args = parsed.args,
  kwargs = parsed.kwargs,
  wheres = parsed.wheres,
)

# Type
# ----

"""
    EP.Type(name = EP.anything, curlies = EP.anything, wheres = EP.anything)

Parses the following:
```julia
a
a{b, c}
a{d} where d
```

# Examples

```jldoctest
julia> using ExprParsers

julia> parser = EP.Type()
EP.Type(
  name    = ExprParsers.Isa{Any}()
  curlies = ExprParsers.Isa{Any}()
  wheres  = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(Array{T, 2} where T))
EP.Type_Parsed(
  name    = :Array
  curlies = Any[:T, 2]
  wheres  = Any[:T]
)
julia> parse_expr(parser, :(f(1,2)))
ERROR: ParseError: Cannot parse expr `f(1, 2)` as reference: expr.head `call` not in `[:curly, :.]`.
```
"""
@exprparser struct Type
  name = anything
  curlies = anything = []
  wheres = anything = []
end

parse_expr(parser::Type, expr::Base.Symbol) = parse_expr(parser, name = expr, curlies = [], wheres = [])
parse_expr(parser::Type, expr::Base.Type) = parse_expr(parser, name = expr, curlies = [], wheres = [])
function parse_expr(parser::Type, expr::Base.Expr)
  expr′, wheres = split_where(expr)
  ref = parse_expr(Reference(), expr′)
  parse_expr(parser, name = ref.name, curlies = ref.curlies, wheres = wheres)
end

to_expr(parsed::Type_Parsed) = _to_expr_Type(parsed, Val{nonempty(parsed.curlies)}(), Val{nonempty(parsed.wheres)}())
_to_expr_Type(parsed, curlies::True, wheres::True) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)} where {$(to_expr(parsed.wheres)...)}
)
_to_expr_Type(parsed, curlies::True, wheres::False) = :(
  $(to_expr(parsed.name)){$(to_expr(parsed.curlies)...)}
)
_to_expr_Type(parsed, curlies::False, wheres::True) = :(
  $(to_expr(parsed.name)) where {$(to_expr(parsed.wheres))}
)
_to_expr_Type(parsed, curlies::False, wheres::False) = to_expr(parsed.name)

# TypeRange
# ---------

"""
    EP.TypeRange(lb = EP.anything, name = EP.anything, ub = EP.anything)

Note: Construct with `typevar = EP.anysymbol` to guarantee that only plain symbols can be used as type variable.

Parses the following.
```
TypeVar >: LowerBound
TypeVar <: UpperBound
LowerBound <: TypeVar <: UpperBound
```

# Examples
```
julia> using ExprParsers

julia> parser = EP.TypeRange(name = EP.anysymbol)
EP.TypeRange(
  lb   = ExprParsers.Isa{Any}()
  name = ExprParsers.Isa{Symbol}()
  ub   = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(Int <: T <: Number))
EP.TypeRange_Parsed(
  lb   = :Int
  name = :T
  ub   = :Number
)
julia> parse_expr(parser, :(T <: Number))
EP.TypeRange_Parsed(
  lb   = Union{}
  name = :T
  ub   = :Number
)
julia> parse_expr(parser, :(Int <: 4)) # CAUTION: when using single `<:`, the order is decisive!
EP.TypeRange_Parsed(
  lb   = Union{}
  name = :Int
  ub   = 4
)
julia> parse_expr(parser, :(4 >: Int))
ERROR: ParseError: Expected type `Symbol`, got `4` of type `Int64`.
julia> parse_expr(parser, :(4 <: Int))
ERROR: ParseError: Expected type `Symbol`, got `4` of type `Int64`.
```
"""
@exprparser struct TypeRange
  # naming is analog to Base.TypeVar
  lb = anything = Union{}
  name = anything
  ub = anything = Any
end

parse_expr(parser::TypeRange, expr::Base.Symbol) = parse_expr(parser, lb = Union{}, name = expr, ub = Any)
function parse_expr(parser::TypeRange, expr::Base.Expr)
  if expr.head == :>:
    parse_expr(parser, lb = expr.args[2], name = expr.args[1], ub = Any)
  elseif expr.head == :<:
    parse_expr(parser, lb = Union{}, name = expr.args[1], ub = expr.args[2])
  elseif expr.head == :comparison && length(expr.args) == 5 && expr.args[2] == expr.args[4] == :<:
    parse_expr(parser, lb = expr.args[1], name = expr.args[3], ub = expr.args[5])
  else
    throw(ParseError("TypeRange only parses <:, >: constructs, but got: $expr"))
  end
end

to_expr(parsed::TypeRange_Parsed) = _to_expr_TypeRange(parsed.lb, parsed.name, parsed.ub)
_to_expr_TypeRange(::Base.Type{Union{}}, name, ::Base.Type{Any}) = name
_to_expr_TypeRange(lb, name, ::Base.Type{Any}) = :($name >: $lb)
_to_expr_TypeRange(::Base.Type{Union{}}, name, ub) = :($name <: $ub)
_to_expr_TypeRange(lb, name, ub) = :($lb <: $name <: $ub)

# TypeAnnotation
# --------------

"""
    EP.TypeAnnotation(name = EP.anything, type = EP.anything)

Parses the following
```julia
a
a::B
::B
```

# Examples
```jldoctest
julia> using ExprParsers

julia> parser = EP.TypeAnnotation(type = :Int)
EP.TypeAnnotation(
  name = ExprParsers.Isa{Any}()
  type = :Int
)
julia> parse_expr(parser, :(::Int))
EP.TypeAnnotation_Parsed(
  name = nothing
  type = :Int
)
julia> parse_expr(parser, :(a::Int))
EP.TypeAnnotation_Parsed(
  name = :a
  type = :Int
)
julia> parse_expr(parser, :(a::String))
ERROR: ParseError: Using default `==` comparison, but parser `:Int` ≠ value `:String`.
```
"""
@exprparser struct TypeAnnotation
  name = anything
  type = anything
end

function parse_expr(parser::TypeAnnotation, expr::Base.Expr)
  @passert expr.head == :(::)
  @passert length(expr.args) <= 2
  if length(expr.args) == 1
    parse_expr(parser, name = nothing, type = expr.args[1])
  else
    parse_expr(parser, name = expr.args[1], type = expr.args[2])
  end
end
function to_expr(parsed::TypeAnnotation_Parsed)
  if isnothing(parsed.name)
    :(::$(to_expr(parsed.type)))
  else
    :($(to_expr(parsed.name))::$(to_expr(parsed.type)))
  end
end

# Arg
# ---

struct NoDefault end
const nodefault = NoDefault()

@doc raw"""
    EP.Arg(name = EP.anything, type = EP.anything, default = EP.anything)

A missing default value is indicated by the special variable `EP.nodefault` which is the unique instance of the
singleton type `EP.NoDefault`.

Parses the following
```julia
a
a::B
::B
a = c  # only :($(Expr(:kw, :a, :c))), not plain :(a = c)
a::B = c  # only :($(Expr(:kw, :(a::B), :c))), not plain :(a::B = c)
```

# Examples

```jldoctest
julia> using ExprParsers

julia> parser = EP.Arg()
EP.Arg(
  name    = ExprParsers.Isa{Any}()
  type    = ExprParsers.Isa{Any}()
  default = ExprParsers.Isa{Any}()
)
julia> parse_expr(parser, :(a::B))
EP.Arg_Parsed(
  name    = :a
  type    = :B
  default = ExprParsers.NoDefault()
)
julia> parse_expr(parser, Expr(:kw, :a, 3))
EP.Arg_Parsed(
  name    = :a
  type    = Any
  default = 3
)
julia> parse_expr(parser, :(a = 3))
ERROR: ParseError: AnyOf could not parse expr `a = 3` with any of the parsers `(ExprParsers.Isa{Symbol}(), EP.TypeAnnotation(name=ExprParsers.Isa{Any}(), type=ExprParsers.Isa{Any}()))`. Arg should either be a Symbol or a TypeAnnotation.
julia> parse_expr(parser, :(f(a)))
ERROR: ParseError: AnyOf could not parse expr `f(a)` with any of the parsers `(ExprParsers.Isa{Symbol}(), EP.TypeAnnotation(name=ExprParsers.Isa{Any}(), type=ExprParsers.Isa{Any}()))`. Arg should either be a Symbol or a TypeAnnotation.
```
"""
@exprparser struct Arg
  # name == nothing will be treated as anonymous type annotation
  name = anything
  type = anything = Any
  default = anything = nodefault
end

function parse_expr(parser::Arg, expr::Base.Symbol)
  parse_expr(parser, name = expr, type = Any, default = nodefault)
end

const _arg_left_handside_parser = AnyOf(anysymbol, TypeAnnotation(),
                                        errormessage = "Arg should either be a Symbol or a TypeAnnotation.")
function parse_expr(parser::Arg, expr::Base.Expr)
  if expr.head == :kw
    name, type = @match(parse_expr(_arg_left_handside_parser, expr.args[1])) do f
      f(symbol::Base.Symbol) = symbol, Any
      function f(parsed::TypeAnnotation_Parsed)
        @passert issomething(parsed.name) "need variable name if default is given"
        parsed.name, parsed.type
      end
    end
    parse_expr(parser, name = name, type = type, default = expr.args[2])
  elseif expr.head == :...
    # Vararg is not a valid Type https://github.com/JuliaLang/julia/issues/30995
    parse_expr(parser, name = expr.args[1], type = Vararg, default = nodefault)
  else
    name, type = @match(parse_expr(_arg_left_handside_parser, expr)) do f
      f(symbol::Base.Symbol) = symbol, Any
      f(parsed::TypeAnnotation_Parsed) = parsed.name, parsed.type
    end
    parse_expr(parser, name = name, type = type, default = nodefault)
  end
end

to_expr(parsed::Arg_Parsed) = _to_expr_Arg(parsed.name, parsed.type, parsed.default)
_to_expr_Arg(::Nothing, ::Base.Type{Any}, ::NoDefault) = :(::Any)
_to_expr_Arg(::Nothing, type, ::NoDefault) = :(::$(to_expr(type)))
_to_expr_Arg(name, ::Base.Type{Any}, ::NoDefault) = to_expr(name)
function _to_expr_Arg(name, type, ::NoDefault)
  if type === Vararg
    Base.Expr(:..., to_expr(name))
  else
    :($(to_expr(name))::$(to_expr(type)))
  end
end
_to_expr_Arg(name, ::Base.Type{Any}, default) = Base.Expr(:kw, to_expr(name), to_expr(default))
_to_expr_Arg(name, type, default) = Base.Expr(:kw, :($(to_expr(name))::$(to_expr(type))), to_expr(default))


text:
"""
    create_directions(ndirection, flat)

Create `ndirection` directions unformly distributed in sphere.
If `flat` is true, the directions lie in the plane instead of on the sphere.
"""
function create_directions(ndirection::Int; flat::Bool = false)
    if flat
        # Create directions (unformly distributed on unit circle)
        θ = 2π * collect(0:ndirection-1)' / ndirection
        directions = [
            cos.(θ)
            sin.(θ)
            zeros(1, ndirection)
        ]
    else
        # Create Fibonacci distributed directions
        directions = create_fibonacci_sphere(ndirection)
    end
    directions
end

ext: jl


text:
module TestElixirs

using LinearAlgebra
using Test
using Trixi

import ForwardDiff

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

cmd = string(Base.julia_cmd())
coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", cmd)

@testset "Special elixirs" begin
  @testset "Convergence test" begin
    if !coverage
      @timed_testset "tree_2d_dgsem" begin
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"), 3, initial_refinement_level=2)
        @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05)
      end

      @timed_testset "structured_2d_dgsem" begin
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_extended.jl"), 3, cells_per_dimension=(5, 9))
        @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05)
      end

      @timed_testset "p4est_2d_dgsem" begin
        # Run convergence test on unrefined mesh
        no_refine = @cfunction((p4est, which_tree, quadrant) -> Cint(0), Cint, (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p4est_quadrant_t}))
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), 2, refine_fn_c=no_refine)
        @test isapprox(mean_convergence[:linf], [3.2, 3.2, 4.0, 3.7], rtol=0.05)
      end

      @timed_testset "structured_3d_dgsem" begin
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_3d_dgsem", "elixir_advection_basic.jl"), 2, cells_per_dimension=(7, 4, 5))
        @test isapprox(mean_convergence[:l2], [4.0], rtol=0.05)
      end

      @timed_testset "p4est_3d_dgsem" begin
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "p4est_3d_dgsem", "elixir_advection_unstructured_curved.jl"), 2, initial_refinement_level=0)
        @test isapprox(mean_convergence[:l2], [2.7], rtol=0.05)
      end

      @timed_testset "paper_self_gravitating_gas_dynamics" begin
        mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "paper_self_gravitating_gas_dynamics", "elixir_eulergravity_convergence.jl"), 2, tspan=(0.0, 0.25), initial_refinement_level=1)
        @test isapprox(mean_convergence[:l2], 4 * ones(4), atol=0.4)
      end
    else
      # Without coverage, just run simple convergence tests to cover
      # the convergence test logic
      @test_nowarn convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01))
      @test_nowarn convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"), 2, initial_refinement_level=0, tspan=(0.0, 0.1))
      @test_nowarn convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_basic.jl"), 2, tspan=(0.0, 0.01))
      @test_nowarn convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_advection_extended.jl"), 2, cells_per_dimension=(1, 1), tspan=(0.0, 0.1))
    end
  end

  @timed_testset "Test linear structure (2D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_hypdiff_lax_friedrichs.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    # check whether the user can modify `b` without changing `A`
    x = vec(ode.u0)
    Ax = A * x
    @. b = 2 * b + x
    @test A * x ≈ Ax
  end

  @testset "Test Jacobian of DG (2D)" begin
    @timed_testset "Linear advection" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_extended.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=2)
      A, _ = linear_structure(semi)

      J = jacobian_ad_forward(semi)
      @test Matrix(A) ≈ J
      λ = eigvals(J)
      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

      J = jacobian_fd(semi)
      @test Matrix(A) ≈ J
      λ = eigvals(J)
      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
    end

    @timed_testset "Compressible Euler equations" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_euler_density_wave.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=1)

      J = jacobian_ad_forward(semi)
      λ = eigvals(J)
      @test maximum(real, λ) < 7.0e-7

      J = jacobian_fd(semi)
      λ = eigvals(J)
      @test maximum(real, λ) < 7.0e-3

      # This does not work yet because of the indicators...
      @test_skip begin
        trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_euler_shockcapturing.jl"),
          tspan=(0.0, 0.0), initial_refinement_level=1)
        jacobian_ad_forward(semi)
      end

      @timed_testset "DGMulti (weak form)" begin
        gamma = 1.4
        equations = CompressibleEulerEquations2D(gamma)
        initial_condition = initial_condition_density_wave

        solver = DGMulti(polydeg = 5, element_type = Quad(), approximation_type = SBP(),
                        surface_integral = SurfaceIntegralWeakForm(flux_central),
                        volume_integral = VolumeIntegralWeakForm())

        # StartUpDG.uniform_mesh is on [-1, 1]^ndims by default
        cells_per_dimension = (2, 2)
        vertex_coordinates, EToV = StartUpDG.uniform_mesh(solver.basis.elementType, cells_per_dimension...)
        mesh = VertexMappedMesh(vertex_coordinates, EToV, solver, is_periodic=(true,true))

        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

        J = jacobian_ad_forward(semi)
        λ = eigvals(J)
        @test maximum(real, λ) < 7.0e-7
      end

      @timed_testset "DGMulti (SBP, flux differencing)" begin
        gamma = 1.4
        equations = CompressibleEulerEquations2D(gamma)
        initial_condition = initial_condition_density_wave

        solver = DGMulti(polydeg = 5, element_type = Quad(), approximation_type = SBP(),
                        surface_integral = SurfaceIntegralWeakForm(flux_central),
                        volume_integral = VolumeIntegralFluxDifferencing(flux_central))

        # StartUpDG.uniform_mesh is on [-1, 1]^ndims by default
        cells_per_dimension = (2, 2)
        vertex_coordinates, EToV = StartUpDG.uniform_mesh(solver.basis.elementType, cells_per_dimension...)
        mesh = VertexMappedMesh(vertex_coordinates, EToV, solver, is_periodic=(true,true))

        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

        J = jacobian_ad_forward(semi)
        λ = eigvals(J)
        @test maximum(real, λ) < 7.0e-7
      end
    end

    @timed_testset "MHD" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_mhd_alfven_wave.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=0)
      @test_nowarn jacobian_ad_forward(semi)
    end
  end

  @timed_testset "Test linear structure (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end

  @timed_testset "Test Jacobian of DG (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, _ = linear_structure(semi)

    J = jacobian_ad_forward(semi)
    @test Matrix(A) ≈ J

    J = jacobian_fd(semi)
    @test Matrix(A) ≈ J
  end

  @testset "AD using ForwardDiff" begin
    @timed_testset "Euler equations 1D" begin
      function entropy_at_final_time(k) # k is the wave number of the initial condition
        equations = CompressibleEulerEquations1D(1.4)
        mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level=3, n_cells_max=10^4)
        solver = DGSEM(3, flux_hll, VolumeIntegralFluxDifferencing(flux_ranocha))
        initial_condition = (x, t, equations) -> begin
          rho = 2 + sinpi(k * sum(x))
          v1  = 0.1
          p   = 10.0
          return prim2cons(SVector(rho, v1, p), equations)
        end
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                            uEltype=typeof(k))
        ode = semidiscretize(semi, (0.0, 1.0))
        summary_callback = SummaryCallback()
        analysis_interval = 100
        analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
        alive_callback = AliveCallback(analysis_interval=analysis_interval)
        callbacks = CallbackSet(
            summary_callback,
            analysis_callback,
            alive_callback
        )
        sol = solve(ode, SSPRK43(), callback=callbacks)
        Trixi.integrate(entropy, sol.u[end], semi)
      end
      ForwardDiff.derivative(entropy_at_final_time, 1.0) ≈ -0.4524664696235628
    end

    @timed_testset "Linear advection 2D" begin
      function energy_at_final_time(k) # k is the wave number of the initial condition
        equations = LinearScalarAdvectionEquation2D(0.2, -0.7)
        mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)
        solver = DGSEM(3, flux_lax_friedrichs)
        initial_condition = (x, t, equation) -> begin
          x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)
          return SVector(sinpi(k * sum(x_trans)))
        end
        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                            uEltype=typeof(k))
        ode = semidiscretize(semi, (0.0, 1.0))
        summary_callback = SummaryCallback()
        analysis_interval = 100
        analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
        alive_callback = AliveCallback(analysis_interval=analysis_interval)
        stepsize_callback = StepsizeCallback(cfl=1.6)
        callbacks = CallbackSet(
            summary_callback,
            analysis_callback,
            alive_callback,
            stepsize_callback
        )
        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), save_everystep=false, adaptive=false, dt=1.0, callback=callbacks)
        Trixi.integrate(energy_total, sol.u[end], semi)
      end
      ForwardDiff.derivative(energy_at_final_time, 1.0) ≈ 1.4388628342896945e-5
    end

    @timed_testset "elixir_euler_ad.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "special_elixirs", "elixir_euler_ad.jl"))
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
# Heuristic IRT module

function icc(θ, a, b)
    logistic(a * (θ - b))
end

function icc(x, θ, a, b)
    if x == 0
        p = 1 - icc(θ, a, b[1])
    elseif x == length(b)
        p = icc(θ, a, b[end])
    else
        p = icc(θ, a, b[x]) - icc(θ, a, b[x + 1])
    end
    return p
end

function generate_response(θ, a, b)
    X = Matrix{Int64}(undef, length(θ), length(a))
    for j in axes(X, 2)
        x = 0:1:length(b[j])
        for i in axes(X, 1)
            ps = cumsum(icc.(x, θ[i], a[j], Ref(b[j])))
            u = rand(Uniform(0, 1), 1)
            X[i, j] = count(ps .< u)
        end
    end
    return X
end

"""
    v2vv
Transform vector to vector of vector
"""
function v2vv(v)
    return [[i] for i in v]
end

# Heuristic IRT parameter estimation based on the polychoric correlations.

mutable struct HeuristicIRT
    a
    d
    b
end

function transform_FA_IRT(α::T, τ::Vector{T}) where T <: Real
    σ = √(1 - α^2)
    a = α / σ
    d = -τ ./ σ
    b = -d ./ a
    return a, d[begin+1:end-1], b[begin+1:end-1]
end

function heuristicIRT(X; method = :em)
    J = size(X, 2)
    r = Matrix{Float64}(undef, J, J)
    r[diagind(r)] .= 1.0
    τ = Vector{Vector{Float64}}(undef, J)
    for i in 1:J
        x = @view X[:, i]
        for j in i:J
            y = @view X[:, j]
            p = polyc(x, y; verbose = false)
            r[j, i] = p.ρ
            r[i, j] = r[j, i]
            if j == J
                τ[i] = p.τ₁
            end
        end
    end
    # FA
    replace_diagonal!(r)
    f = fit(FactorAnalysis, r; mean = fill(0.0, J), maxoutdim = 1, method = method)
    α = loadings(f)
    tp = transform_FA_IRT.(α, τ)
    return HeuristicIRT(getindex.(tp, 1), getindex.(tp, 2), getindex.(tp, 3))
end


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using Optim
using Random
using Distributions
using StatsBase
using Plots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*S*I
        du[2] = β*S*I - γ*I
        du[3] = γ*I
    end
    nothing
end;

tmax = 10000.0
tspan = (0.0,tmax);

cb_ss = TerminateSteadyState();

N = 1000.0;
u0 = [990.0,10.0,0.0];
p = [0.0005,0.25]; # β,γ

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

sol_ode = solve(prob_ode,Tsit5())
plot(sol_ode(0:1:40.0),
     xlabel="Time",
     ylabel="Number",
     labels=["S" "I" "R"])

lb = 0.00005
ub = 0.001
δ = 0.0000001;

θ = lb:δ:ub;

M = function(θ)
  p = prob_ode.p
  p[1] = θ
  prob = remake(prob_ode;p=p)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  ϕ = sol[end][3]/N
  ϕ
end;

using Roots
M_analytic = function(θ)
  _,γ = prob_ode.p
  β = θ
  R₀ = β*N/γ
  S0 = prob_ode.u0[1]/N
  f(ϕ) = 1-S0*exp(-R₀*ϕ)-ϕ
  ϕ = find_zero(f,(0.0,1.0))
  ϕ
end;

@benchmark M(p[1])

@benchmark M_analytic(p[1])

ϕ = M.(θ);

R₀ = θ*N/p[2];

function twiny(sp::Plots.Subplot)
    sp[:top_margin] = max(sp[:top_margin], 30Plots.px)
    plot!(sp.plt, inset = (sp[:subplot_index], bbox(0,0,1,1)))
    twinsp = sp.plt.subplots[end]
    twinsp[:xaxis][:mirror] = true
    twinsp[:background_color_inside] = RGBA{Float64}(0,0,0,0)
    Plots.link_axes!(sp[:yaxis], twinsp[:yaxis])
    twinsp
end
twiny(plt::Plots.Plot = current()) = twiny(plt[1]);

plot(θ, ϕ,
     xlabel="Infectivity parameter, β",
     ylabel="Final size, ϕ",
     legend=false)
pl = twiny()
plot!(pl, R₀, ϕ, xlabel = "Reproductive number, R₀", legend = false)

using KernelDensity
k_kd = kde_lscv(ϕ;boundary=(0.0,1.0),kernel=Normal)
q1star_kd = [pdf(k_kd, ϕᵢ) for ϕᵢ in ϕ];

using KernelDensityEstimate
k_kde = KernelDensityEstimate.kde!(ϕ)
q1star_kde = evaluateDualTree(k_kde,ϕ);

using AverageShiftedHistograms
k_ash = ash(ϕ)
q1star_ash = [AverageShiftedHistograms.pdf(k_ash,ϕᵢ) for ϕᵢ in ϕ];

using MultiKDE
k_mkde = KDEUniv(ContinuousDim(), 0.01, ϕ, MultiKDE.gaussian)
q1star_mkde = [MultiKDE.pdf(k_mkde, ϕᵢ, keep_all=false) for ϕᵢ in ϕ];

plot(ϕ,q1star_kd,label="KernelDensity",xlabel="Final size",ylabel="Density")
plot!(ϕ,q1star_kde,label="KernelDensityEstimate")
plot!(ϕ,q1star_ash,label="AverageShiftedHistograms")
plot!(ϕ,q1star_mkde,label="MultiKDE",legend=:top)

α = 0.5
q₂ = Distributions.pdf.(Uniform(0.05,0.1),ϕ)
L₁(θᵢ) = 1.0
L₂(ϕᵢ) = 1.0;

plot(ϕ,q₂,label="Prior on outputs, q₂")
plot!(ϕ,q1star_kd,label="Induced prior on outputs, q₁*")

w_kd = (q₂ ./ q1star_kd).^(1-α) .* L₁.(θ) .* L₂.(ϕ);

lw_kd = (1-α) .* log.(q₂) .- (1-α) .* log.(q1star_kd) .+ log.(L₁.(θ)) .+ log.(L₂.(ϕ));

w_kde = (q₂ ./ q1star_kde).^(1-α)
w_ash = (q₂ ./ q1star_ash).^(1-α)
w_mkde = (q₂ ./ q1star_mkde).^(1-α);

Random.seed!(123)
n = 1000
πθ = StatsBase.sample(θ, Weights(w_kd),n,replace=true);

mean(πθ),std(πθ)

histogram(πθ,
          bins=25,
          legend=false,
          xlim=(lb,ub),
          xlabel="Infectivity parameter, β",
          ylabel="Count")

w_norm = w_kd ./ sum(w_kd)
ess = 1.0/sum(w_norm .^ 2)

using ThreadsX
θ_fine = lb:1e-10:ub
# We use ThreadsX instead of ϕ_fine = M_analytic.(θ_fine)
ϕ_fine = ThreadsX.collect(M_analytic(θᵢ) for θᵢ in θ_fine)
k_kd_fine = kde_lscv(ϕ_fine;boundary=(0.0,1.0),kernel=Normal)
#k_kde_fine = KernelDensityEstimate.kde!(ϕ_fine)
k_ash_fine = ash(ϕ_fine)
k_mkde_fine = KDEUniv(ContinuousDim(), 0.01, ϕ_fine, MultiKDE.gaussian)
q1star_kd_fine = [pdf(k_kd_fine, ϕᵢ) for ϕᵢ in ϕ]
#q1star_kde_fine = evaluateDualTree(k_kde_fine,ϕ_fine)
q1star_ash_fine = [AverageShiftedHistograms.pdf(k_ash_fine,ϕᵢ) for ϕᵢ in ϕ]
q1star_mkde_fine = ThreadsX.collect(MultiKDE.pdf(k_mkde_fine, ϕᵢ, keep_all=false) for ϕᵢ in ϕ);

l = @layout [a; b; c]
pl1 = plot(ϕ,q1star_kd_fine,label="KernelDensity fine",xlabel="Final size",ylabel="Density")
plot!(pl1,ϕ,q1star_kd,label="KernelDensity coarse")
pl2=plot(ϕ,q1star_ash_fine,label="AverageShiftedHistograms fine")
plot!(pl2,ϕ,q1star_ash,label="AverageShiftedHistograms coarse")
pl3=plot(ϕ,q1star_mkde,label="MultiKDE")
plot!(pl3,ϕ,q1star_mkde_fine,label="MultiKDE fine")
plot(pl1, pl2, pl3, layout=l)


text:
#=
Popis
  Funkcia vypocita hodnotu ekentricity pre dany vrchol.

Parametre
  index: index vrchola :: Int64
  g: datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String

Navratova hodnota
  vert_eccentricity: hodnota ekcentricity pre dany vrchol :: Float64

Poznamky
=#

function vertex_eccentricity(index::Int64, g::GenericGraph; dist_key::AbstractString="weight")

  epi = AttributeEdgePropertyInspector{Float64}(dist_key)

  return vert_eccentricity = maximum(filter(x -> x != Inf,dijkstra_shortest_paths(g, epi, [get_vertex(index, g)]).dists))

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
using Documenter
using Distances
using Ripserer

if VERSION ≥ v"1.7-DEV" || VERSION < v"1.6-DEV" || !Sys.islinux()
    @warn "Doctests only run on Linux and Julia 1.6"
else
    DocMeta.setdocmeta!(
        Ripserer, :DocTestSetup, :(using Ripserer; using Distances); recursive=true
    )
    doctest(Ripserer; fix=true)
end


text:
module TestEmptyProxy

include("preamble.jl")

@testset begin
    @test collect(Empty(Vector)) == Union{}[]
    @test collect(Empty(Dict)) == Union{}[]
    @test collect(Empty(Set)) == Union{}[]
end

end  # module


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This file contains a description for a function, including the function's
# input/output arguments.
# This gives compiler accurate information about the function without analyzing
# it. User can easily describe his/her own functions by adding entries to the
# vector below.

# A function is assumed to have the following form
#      f(argument1, argument2, ...)
# In the description, its return result is represented by number 0,
# and its arguments by number 1, 2, and so on.
# Some arguments may be output as well, which will be described by
# the "output" field.
# ASSUMPTION: all the described is "must" information. That is, the information
# is true in all cases.
# ISSUE: depending on the specific inputs, a function "may" have different
# behavior. For example, it may update some arguments, or not update them.
# How to represent and use some "may" information?

immutable FunctionDescription
    module_name     :: AbstractString # Module of the function. It is "nothing" if it is not in a Julia module, e.g. if this is a C function
    function_name   :: AbstractString # Name of the function
    argument_types  :: Tuple          # Tuple of the function arguments' types
    output          :: Set            # The arguments updated by the function
end

const UPDATED_NONE  = Set()
const VECTOR_OR_NUM = Tuple{Vector{}, Number}
const MATRIX_OR_NUM = Tuple{AbstractMatrix{}, Number}

const element_wise_divide1_Desc = FunctionDescription(
    "Main",
    "./",
    (VECTOR_OR_NUM, Vector{}),        # The arguments must be vectors
    UPDATED_NONE,                     # No argument is updated
)

const star_Desc = FunctionDescription(
    "Main",
    "*",                              # *(A::SparseMatrixCSC, x::Vector)
    (SparseMatrixCSC, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const star1_Desc = FunctionDescription(
    "Main",
    "*",
    (Number, SparseMatrixCSC, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const star2_Desc = FunctionDescription(
    "Main",
    "*",
    (Number, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const dot_Desc = FunctionDescription(
    "Main",
    "dot",                            # Dot(x::Vector, y::Vector)
    (Vector{}, Vector{}),             # The arguments must be vectors
    UPDATED_NONE,                     # No argument is updated
)

const copy_Desc = FunctionDescription(
    "Main",
    "copy",
    (Vector{}, ),                     # The arguments must be a vector
    UPDATED_NONE,                     # No argument is updated
)

const add_vector_Desc = FunctionDescription(
    "Main",
    "+",
    (VECTOR_OR_NUM, VECTOR_OR_NUM),
    UPDATED_NONE,                     # No argument is updated
)

const add_matrix_Desc = FunctionDescription(
    "Main",
    "+",
    (MATRIX_OR_NUM, MATRIX_OR_NUM),
    UPDATED_NONE,                     # No argument is updated
)

const sub_vector_Desc = FunctionDescription(
    "Main",
    "-",
    (VECTOR_OR_NUM, VECTOR_OR_NUM),
    UPDATED_NONE,                     # No argument is updated
)

const norm_Desc = FunctionDescription(
    "Main",
    "norm",
    (AbstractArray, ),
    UPDATED_NONE,                     # No argument is updated
)

const spones_Desc = FunctionDescription(
    "Main",
    "spones",
    (AbstractSparseMatrix{}, ),
    UPDATED_NONE,                     # No argument is updated
)

const size_Desc = FunctionDescription(
    "Main",
    "size",
    (AbstractArray, Number),
    UPDATED_NONE,                     # No argument is updated
)

const Array_Desc = FunctionDescription(
    "Main",
    "Array",                          # Array(Any, Integer)
    (Any, Integer),
    UPDATED_NONE,                     # No argument is updated
)

const fill!_Desc = FunctionDescription(
    "Main",
    "fill!",
    (AbstractArray, Number),
    Set(1),                           # argument 1 is updated
)

const max_Desc = FunctionDescription(
    "Main",
    "max",
    (Vector{}, Number),
    UPDATED_NONE,                     # No argument is updated
)

const scale_Desc = FunctionDescription(
    "Main",
    "scale",
    (AbstractMatrix{}, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const sum_Desc = FunctionDescription(
    "Main",
    "sum",
    (AbstractArray, ),
    UPDATED_NONE,                     # No argument is updated
)

const convert_Desc = FunctionDescription(
    "Main",
    "convert",
    (Array{Float64,1}, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const eltype_Desc = FunctionDescription(
    "Main",
    "eltype",
    (AbstractArray, ),
    UPDATED_NONE,                     # No argument is updated
)

const inverse_divide_Desc = FunctionDescription(
    "Main",
    "\\",
    (AbstractSparseMatrix{}, Vector{}),
    UPDATED_NONE,                     # No argument is updated
)

const fwdTriSolve!_Desc = FunctionDescription(
    "Base.SparseMatrix",
    "fwdTriSolve!",
    (AbstractSparseMatrix{}, Vector{}),
    Set(2),                           # Argument 2 (the vector) is updated
)

const bwdTriSolve!_Desc = FunctionDescription(
    "Base.SparseMatrix",
    "bwdTriSolve!",
    (AbstractSparseMatrix{}, Vector{}),
    Set(2),                           # Argument 2 (the vector) is updated
)

function_descriptions  = [
    element_wise_divide1_Desc,
    star_Desc,
    star1_Desc,
    star2_Desc,
    dot_Desc,
    copy_Desc,
    add_vector_Desc,
    add_matrix_Desc,
    sub_vector_Desc,
    norm_Desc,
    spones_Desc,
    size_Desc,
    Array_Desc,
    fill!_Desc,
    max_Desc,
    scale_Desc,
    sum_Desc,
    convert_Desc,
    eltype_Desc,
    inverse_divide_Desc,
    fwdTriSolve!_Desc,
    bwdTriSolve!_Desc
]

"""
Takes a module and a function both as Strings. Looks up the specified module as
part of the "Main" module and then looks and returns the Function object
corresponding to the "func" String in that module.
"""
function get_function_from_string(mod :: AbstractString, func :: AbstractString)
    # A module string may have submodules like "Base.SparseMatrix". We need to
    # get module object level by level
    modobj  = eval(:Main)
    modules = split(mod, '.')
    for m in modules
        msym = Symbol(m)
        modobj = eval(:($modobj.$msym))
    end

    fsym   = Symbol(func)
    return eval(:($modobj.$fsym))
end

"""
Convert the function_descriptions table into a dictionary that can be passed to
LivenessAnalysis to indicate which args are unmodified by which functions.
"""
function create_unmodified_args_dict()
    res = Dict{Tuple{Any,Array{DataType,1}}, Array{Int64,1}}()

    for desc in function_descriptions
        num_args    = length(desc.argument_types)   # Get the number of arguments to the functions.
        arg_type_array = collect(desc.argument_types)
        #println(desc.function_name, " arg_type_array = ", arg_type_array, " type = ", typeof(arg_type_array))
        unmodifieds = ones(Int64, num_args)         # desc.output contains "modifies" so we default to true and then turn off based on desc.output.
        for j in desc.output
            unmodifieds[j] = 0
        end
        res[(get_function_from_string(desc.module_name, desc.function_name),
             arg_type_array)] = unmodifieds
    end

    return res
end


text:
using WoodburyMatrices

d  = 2+rand(5)
dl = -rand(4)
du = -rand(4)
M = Tridiagonal(dl, d, du)
F = lufact(M)
U = sprand(5,2,0.2)
V = sprand(2,5,0.2)
C = rand(2,2)
W = Woodbury(F, U, C, V)

src = rand(5, 8)
@test_approx_eq W\src AxisAlgorithms.A_ldiv_B_md(W, src, 1)
src = rand(5, 5, 5)
for dim = 1:3
    dest1 = mapslices(x->W\x, copy(src), dim)
    dest2 = similar(src)
    AxisAlgorithms.A_ldiv_B_md!(dest2, W, src, dim)
    @test_approx_eq dest1 dest2
end


text:
@testset "evaluation tests" begin
    x = MPC.Variable("x", Complex)
    y = MPC.Variable("y", Complex)
    z = MPC.Variable("z", Real)
    b = MPC.Variable("b", Bool)

    ## Polynomial algebra
    pt = Point(SortedDict(x=>2, y=>1+im, z=>0+im, b=>2.1))
    # print(pt)

    p1 = x^2 + 3*y + conj(x) + conj(z) + b
    evaluate(p1, pt) == 10 + 3im

    p2 = y^6 - y^6 + (-y*x*b) / 4 + π*conj(b)

    ## Poly operators
    @test evaluate(x, pt) == 2
    @test evaluate(y, pt) == 1+im
    @test evaluate(conj(y), pt) == 1-im
    @test evaluate(z, pt) == 0
    @test evaluate(b, pt) == 1

    @test evaluate(x*y, pt) == 2+2im
    @test evaluate(p2, pt) == (-1-im)/2 + π

    @test evaluate(real(y), pt) == 1
    @test evaluate(imag(y), pt) == 1

    @test evaluate(x*y*conj(y), pt) == 2*(1+im)*(1-im)
    @test evaluate(x*y + conj(y), pt) == 3+im
    @test evaluate(x*y + 1 + 1im, pt) == 3+3im

    ## Point algebra
    pt = Point(SortedDict(x=>2, y=>1+im, z=>0+im, b=>2.1))
    pt1 = Point(SortedDict(z=>0+im))
    pt2 = Point(SortedDict(x=>3, y=>-1+2im, b=>-1))
    @test merge(pt, pt1, pt2) == Point(SortedDict(x=>5, y=>3im, z=>0+im, b=>1))

    @test pt + pt1 - 3*pt2 == Point(SortedDict(x=>-7, y=>4-5im, z=>0+im, b=>1))

    @test norm(pt2, Inf) == 3
    @test norm(pt2, 1) == 6
    @test norm(pt2, 2) == √(14)

end


text:
# This file is part of FlagSets.jl project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

function flagset_flags end
function flag_bit_map end
function flags end
function flagkeys end

"""
    BitMask(bitmask::Integer)

Structure used internally by `FlagSets` for default constructor of subtypes of `FlagSet`.
"""
struct BitMask{T<:Integer}
    bitmask::T
end

basetype(::Type{<:FlagSet{F,B}}) where {F,B<:Integer} = B

Base.isvalid(::Type{T}, x::Integer) where {T<:FlagSet} = (x & basetype(T)(typemax(T)) == x)
Base.isvalid(::Type{T}, x::BitMask) where {T<:FlagSet} = (x.bitmask & basetype(T)(typemax(T)) == x.bitmask)
Base.isvalid(::Type{T}, x::T) where {T<:FlagSet} = true
Base.isvalid(x::T) where {T<:FlagSet} = true

function Base.isvalid(::Type{T}, x::S) where {S,T<:FlagSet}
    Base.isiterable(S) && all(elt ∈ flags(T) for elt ∈ x)
end

function Base.isvalid(::Type{T}, x::Symbol) where {T<:FlagSet{Symbol}}
    x ∈ flags(T)
end

(::Type{I})(x::FlagSet) where {I<:Integer} = I(x.bitflags)::I
Base.cconvert(::Type{I}, x::FlagSet) where {I<:Integer} = I(x.bitflags)
Base.write(io::IO, x::T) where {T<:FlagSet} = write(io, x.bitflags)
Base.read(io::IO, ::Type{T}) where {T<:FlagSet} = T(BitMask(read(io, basetype(T))))
#Base.isless(x::FlagSet{T}, y::FlagSet{T}) where {T<:Integer} = isless(T(x), T(y))
Base.:|(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags | y.bitflags))
Base.:&(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags & y.bitflags))
Base.:⊻(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags ⊻ y.bitflags))
Base.:~(x::T) where {T<:FlagSet} = setdiff(typemax(T), x)
Base.iszero(x::FlagSet) = isempty(x)

# Iterator interface

function Base.iterate(x::T, xi::Integer = x.bitflags) where {T<:FlagSet}
    iszero(xi) && (return nothing)
    fs = flagset_flags(T)
    nbit = trailing_zeros(xi)
    xi ⊻= typeof(xi)(1) << nbit
    return (fs[nbit+1], xi)
end

Base.length(x::FlagSet) = count_ones(x.bitflags)
Base.isempty(x::FlagSet) = iszero(x.bitflags)

"""
    flags(x::FlagSet)

Return the set of flags in `x` as a `Tuple`.

# Examples
```julia
julia> @symbol_flagset FontFlags bold=1 italic=2 large=4

julia> flags(FontFlags(3))
(:bold, :italic)
```
"""
flags(x::T) where {T<:FlagSet} = NTuple{length(x),eltype(T)}(x)

"""
    flags(T)

Return the set of flags of type `T` as a `Tuple`.

# Examples
```julia
julia> @flagset RoundingFlags RoundDown RoundUp RoundNearest

julia> flags(RoundingFlags)
(RoundingMode{:Down}(), RoundingMode{:Up}(), RoundingMode{:Nearest}())
```
"""
flags(::Type{<:FlagSet})

# Set manipulation
Base.union(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags | y.bitflags))
Base.intersect(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags & y.bitflags))
Base.setdiff(x::T, y::T) where {T<:FlagSet} = T(BitMask(x.bitflags & ~y.bitflags))
Base.issubset(x::T, y::T) where {T<:FlagSet} = (x.bitflags & y.bitflags) == x.bitflags
Base.in(elt, x::T) where {T<:FlagSet} = !iszero(get_flag_bit(T, elt, zero(basetype(T))) & x.bitflags)
Base.:⊊(x::T, y::T) where {T<:FlagSet} = x != y && (x.bitflags & y.bitflags) == x.bitflags

Base.empty(s::FlagSet{B,F}, ::Type{F}) where {B,F} = typeof(s)()
Base.empty(s::FlagSet) = typeof(s)()

# Printing
function Base.show(io::IO, x::T) where {T<:FlagSet}
    compact = get(io, :compact, false)::Bool
    type_p = (get(io, :typeinfo, false) == T)::Bool
    malformed = !isvalid(x)
    if compact || malformed
        xi = Integer(x)
        type_p &= malformed
        type_p || show(io, T)
        type_p || print(io, "(")
        show(io, xi)
        type_p || print(io, malformed ? " #= Invalid code =#)" : ")")
    else
        show(io, T)
        print(io, "([")
        join(io, repr.(x), ", ")
        print(io, "])")
    end
end

function Base.show(io::IO, mime::MIME"text/plain", x::FlagSet)
    if get(io, :compact, false)
        show(io, x)
    else
        invoke(show, Tuple{typeof(io),typeof(mime),AbstractSet}, io, mime,  x)
    end
end

# Printing FlagSet type
function Base.show(io::IO, mime::MIME"text/plain", type::Type{<:FlagSet})
    if isconcretetype(type) && !(get(io, :compact, false))
        print(io, "FlagSet ")
        Base.show_datatype(io, type)
        print(io, ":")
        keys = map(key -> isnothing(key) ? "" : String(key), flagkeys(type))
        if !all(isempty, keys)
            klen = maximum(length, keys)
            keys = map(k -> isempty(k) ? ' '^(klen+2) : k * ' '^(klen-length(k)) * " = ", keys)
        end
        for (flag, key) ∈ zip(flags(type), keys)
            bit = get_flag_bit(type, flag)
            print(io, "\n ", key, repr(bit), " --> ", repr(flag))
        end
    else
        invoke(show, Tuple{typeof(io),typeof(mime),Type}, io, mime, type)
    end
end

function flagset_argument_error(typename, x)
    throw(ArgumentError(string("invalid value for FlagSet $(typename): ", repr(x))))
end

function get_flag_bit(::Type{T}, flag, default) where {T<:FlagSet}
    get(flag_bit_map(T), flag, default)
end

function get_flag_bit(::Type{T}, flag) where {T<:FlagSet}
    not_found = basetype(T)(0x7f)
    val = get(flag_bit_map(T), flag, not_found)
    val == not_found && flagset_argument_error(T, flag)
    val
end

# Default constructors

function (::Type{T})(flag::Symbol, flags::Symbol...) where {T<:FlagSet{Symbol}}
    T((flag, flags...))
end

function (::Type{T})(itr) where {T<:FlagSet}
    Base.isiterable(typeof(itr)) || flagset_argument_error(T, itr)
    bitmask::basetype(T) = 0
    for flag ∈ itr
        bitmask |= get_flag_bit(T, flag)
    end
    T(BitMask(bitmask))
end

function (::Type{T})(bitmask::Integer) where {T<:FlagSet}
    T(BitMask(bitmask))
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# All types and constructors are defined here.

# Structs to define a figure

## Term configuration
Base.@kwdef mutable struct TermConf
    font::String = ""
    size::String = ""
    background::String = ""
end

Base.@kwdef mutable struct PrintConf
    print_term::String       = config[:print][:print_term]
    print_font::String       = config[:print][:print_font]
    print_size::String       = config[:print][:print_size]
    print_linewidth::String  = config[:print][:print_linewidth]
    print_background::String = config[:print][:print_background]
end

## Coordinates
const Coord = Union{AbstractRange{T},Array{T}} where T <: Real
Coord() = Float64[] # return an empty coordinate
const ComplexCoord = Array{T} where T <: Complex

# storage for financial coordinates
struct FinancialCoords
    open::Coord
    low::Coord
    high::Coord
    close::Coord
end
const FCuN = Union{FinancialCoords, Nothing}
FinancialCoords() = nothing

# storage for error coordinates
struct ErrorCoords
    ylow::Coord
    yhigh::Coord
end
const ECuN = Union{ErrorCoords, Nothing}
ErrorCoords() = nothing
ErrorCoords(l::Coord) = ErrorCoords(l,Float64[])

## Curves
# Curve configuration
Base.@kwdef struct CurveConf
    legend::String    = ""
    plotstyle::String = config[:curve][:plotstyle]
    linecolor::String = config[:curve][:linecolor]
    linewidth::String = config[:curve][:linewidth]
    linestyle::String = config[:curve][:linestyle]
    pointtype::String = config[:curve][:pointtype]
    pointsize::String = config[:curve][:pointsize]
    fillstyle::String = config[:curve][:fillstyle]
    fillcolor::String = config[:curve][:pointsize]
end

# A curve is a set of coordinates and a configuration
Base.@kwdef mutable struct Curve
    x::Coord = Coord()
    y::Coord = Coord()
    Z::Coord = Coord()
    F::FCuN = FinancialCoords()
    E::ECuN = ErrorCoords()
    conf::CurveConf = CurveConf()
end

# Storage for axes configuration
Base.@kwdef mutable struct AxesConf
    title::String      = ""
    xlabel::String     = ""
    ylabel::String     = ""
    zlabel::String     = ""
    fillstyle::String  = config[:axes][:fillstyle]
    grid::String       = config[:axes][:grid]
    boxwidth::String   = config[:axes][:boxwidth]
    keyoptions::String = config[:axes][:keyoptions]
    axis::String       = config[:axes][:axis]
    xrange::String     = config[:axes][:xrange]
    yrange::String     = config[:axes][:yrange]
    zrange::String     = config[:axes][:zrange]
    xzeroaxis::String  = config[:axes][:xzeroaxis]
    yzeroaxis::String  = config[:axes][:yzeroaxis]
    zzeroaxis::String  = config[:axes][:zzeroaxis]
    palette::String    = config[:axes][:palette]
    output::String     = config[:axes][:output]
end

# At the top level, a figure is a handle, an axes configuration, and a
# set of curves.
const Handle = Union{Int,Nothing}  # handle type

mutable struct Figure
    handle::Handle               # each figure has a unique handle
    term::TermConf               # term options
    print::PrintConf             # print optinos
    axes::AxesConf               # figure configuration
    curves::Union{Nothing,Vector{Curve}} # a vector of curves
    svg::String          # SVG data returned by gnuplot (used in IJulia)
    gpcom::String        # a gnuplot command to run before plotting
end
# Construct an empty figure with given handle
Figure(handle) = Figure(handle,TermConf(),PrintConf(),AxesConf(),nothing,"","")

# We need a global variable to keep track of gnuplot's state
mutable struct GnuplotState
    current                      # current figure -- "pointer" to figs
    process                      # the gnuplot process
    gp_stdout::String            # last data read from gnuplot's stdout
    gp_stderr::String            # last data read from  gnuplot's stderr
    gp_lasterror::String         # gnuplot's last error output
    gp_error::Bool               # true if last command resulted in gp error
    figs::Vector{Figure}          # storage for all figures
end

GnuplotState() = GnuplotState(nothing,nothing,"","","",false,Figure[])

# Types for dealing with arguments to `plot`
Base.@kwdef mutable struct PlotArgs
    legend::String          = ""
    title::String           = ""
    xlabel::String          = ""
    ylabel::String          = ""
    zlabel::String          = ""
    plotstyle::String       = config[:curve][:plotstyle]
    linecolor::String       = config[:curve][:linecolor]
    linewidth::String       = config[:curve][:linewidth]
    linestyle::String       = config[:curve][:linestyle]
    pointtype::String       = config[:curve][:pointtype]
    pointsize::String       = config[:curve][:pointsize]
    fillcolor::String       = config[:curve][:fillcolor]
    fillstyle::String       = config[:curve][:fillstyle]
    grid::String            = config[:axes][:grid]
    boxwidth::String        = config[:axes][:boxwidth]
    keyoptions::String      = config[:axes][:keyoptions]
    axis::String            = config[:axes][:axis]
    xrange::String          = config[:axes][:xrange]
    yrange::String          = config[:axes][:yrange]
    zrange::String          = config[:axes][:zrange]
    xzeroaxis::String       = config[:axes][:xzeroaxis]
    yzeroaxis::String       = config[:axes][:yzeroaxis]
    zzeroaxis::String       = config[:axes][:zzeroaxis]
    palette::String         = config[:axes][:palette]
    output::String          = config[:axes][:output]
    font::String            = config[:term][:font]
    size::String            = config[:term][:size]
    background::String      = config[:term][:background]
    printfont::String       = config[:print][:print_font]
    printsize::String       = config[:print][:print_size]
    printlinewidth::String  = config[:print][:print_linewidth]
    printbackground::String = config[:print][:print_background]
    gpcom::String           = ""
end


text:
using BenchmarkTools
using CSV
using DataFrames
using PSFModels: Gaussian, AiryDisk, Moffat
using PyCall

@info "Starting PSF evaluation benchmark"

ap = pyimport("astropy")
apm = pyimport("astropy.modeling.models")
@info "Using astropy version $(ap.__version__)"

julia_models = [
    Gaussian(fwhm=(5.3, 4.7)),
    AiryDisk(fwhm=12.7),
    Moffat(fwhm=8.1)
]
astropy_models = [
    apm.Gaussian2D(x_stddev=5.3, y_stddev=4.7),
    apm.AiryDisk2D(radius=12.7),
    apm.Moffat2D(gamma=8.1, alpha=-1)
]

names = ["Gaussian", "AiryDisk", "Moffat"]
jl_times = []
py_times = []

for (jl_mod, py_mod, name) in zip(julia_models, astropy_models, names)
    jl_time = @belapsed $jl_mod(0.5, 0.7)
    push!(jl_times, jl_time)

    py_time = @belapsed $py_mod(0.5, 0.7)
    push!(py_times, py_time)

    @info "$name" PSFModels=jl_time astropy=py_time
end

filename = joinpath(@__DIR__, "results.csv")
DataFrame(name=names, psfmodels=jl_times, astropy=py_times) |> CSV.write(filename)

@info "Results saved to $filename"

@info "Finished PSF evaluation benchmark"

nothing


text:
import Base.Random: rand, srand
import RNG: gen_seed

@inline rand(r::MT19937, ::Type{UInt32}) = mt_get(r)

srand(r::MT19937, seed::Integer) = mt_set!(r, seed % UInt32)
srand(r::MT19937, seed::NTuple{N, UInt32}=gen_seed(UInt32, N)) = mt_set!(r, seed)


text:
"""
    fitonecycle!(learner, nepochs[, lrmax])

Fit `learner` for `nepochs` using a one-cycle learning rate schedule.
The learning rate starts at `lrmax/div` for `pct_start*nepochs` epochs, rising
to `lrmax` and then goes down to `lrmax/div_final` over the remaining duration.

## Keyword arguments

- `wd = 0`: weight decay
- `pct_start = 0.25`: Percentage of time spent raising the learning rate
- `div = 25`: Starting learning rate is `lr_max/div`
- `div_final = 1e5`: Ending learning rate is `lr_max/div_final`
"""
function fitonecycle!(
        learner::Learner, nepochs::Int, maxlr=0.1;
        dataiters=(learner.data.training, learner.data.validation),
        wd=0.,
        kwargs...)

    nsteps = length(learner.data.training)
    scheduler = Scheduler(LearningRate => onecycle(
        nepochs * nsteps,
        maxlr;
        kwargs...))

    wdoptim = wd > 0 ? decay_optim(learner.optimizer, wd) : learner.optimizer
    withfields(learner, optimizer=wdoptim) do
        withcallbacks(learner, scheduler) do
            fit!(learner, nepochs, dataiters)
        end
    end
end

"""
    decay_optim(optim, wd)

Add [`WeightDecay`](#) with value `wd` to optimizer `optim`.
"""
decay_optim(optim, wd) = Optimiser(WeightDecay(wd), optim)
function decay_optim(optim::Optimiser, wd)
    # change weight decay if present, else add it
    i = findfirst(o -> o isa WeightDecay, optim.os)
    if isnothing(i)
        return Optimiser(WeightDecay(wd), optim.os...)
    else
        return Optimiser(WeightDecay(wd), optim.os[1:i - 1]..., optim.os[i + 1:end]...)
    end
end


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

"""
    SimulationSolution

A solution to a spatial simulation problem.
"""
struct SimulationSolution{D<:AbstractDomain}
  domain::D
  realizations::Dict{Symbol,Vector{<:AbstractVector}}
  nreals::Int

  function SimulationSolution{D}(domain, realizations) where {D}
    n = [length(r) for (var, r) in realizations]
    @assert length(unique(n)) == 1 "number of realizations must be unique"
    new(domain, realizations, n[1])
  end
end

SimulationSolution(domain, realizations) =
  SimulationSolution{typeof(domain)}(domain, realizations)

"""
    getindex(solution, var)

Return simulation solution for specific variable `var`
as a vector of realizations.
"""
function Base.getindex(solution::SimulationSolution, var::Symbol)
  solution.realizations[var]
end

function Base.getindex(solution::SimulationSolution{<:RegularGrid}, var::Symbol)
  sz = size(solution.domain)
  [reshape(real, sz) for real in solution.realizations[var]]
end

#---------------
# ITERATOR API
#---------------

"""
    iterate(solution, state=1)

Iterate over realizations in simulation `solution`.
"""
Base.iterate(solution::SimulationSolution, state=1) =
  state > solution.nreals ? nothing : (solution[state], state + 1)

"""
    length(solution)

Return the number of realizations in simulation `solution`.
"""
Base.length(solution::SimulationSolution) = solution.nreals

#----------------
# INDEXABLE API
#----------------

"""
    getindex(solution, ind)

Return the `ind`-th realization of simulation `solution`.
"""
function Base.getindex(solution::SimulationSolution, ind::Int)
  sdomain = solution.domain
  sreals  = solution.realizations
  idata   = OrderedDict([(var, reals[ind]) for (var, reals) in sreals])
  georef(idata, sdomain)
end

Base.getindex(solution::SimulationSolution, inds::AbstractVector{Int}) =
  [getindex(solution, ind) for ind in inds]

"""
    firstindex(solution)

Return the first index of simulation `solution`.
"""
Base.firstindex(solution::SimulationSolution) = 1

"""
    lastindex(solution)

Return the last index of simulation `solution`.
"""
Base.lastindex(solution::SimulationSolution) = length(solution)

# ------------
# IO methods
# ------------
function Base.show(io::IO, solution::SimulationSolution)
  dim = ndims(solution.domain)
  print(io, "$(dim)D SimulationSolution")
end

function Base.show(io::IO, ::MIME"text/plain", solution::SimulationSolution)
  println(io, solution)
  println(io, "  domain: ", solution.domain)
  println(io, "  variables: ", join(keys(solution.realizations), ", ", " and "))
  print(  io, "  N° reals:  ", solution.nreals)
end


text:
using CSV
using DataFrames
using StatsBase

@testset "LooCompare" begin
    using Random
    Random.seed!(129111)

    df = CSV.read(joinpath("data", "WaffleDivorce.csv"), DataFrame)
    df.D = zscore(df.Divorce)
    df.M = zscore(df.Marriage)
    df.A = zscore(df.MedianAgeMarriage)
    data = (D=df.D, A=df.A)

    function lin(a, b, c, x...)
        result = @. a + b * c
        for i in 1:2:length(x)
            @. result += x[i] * x[i + 1]
        end
        return result
    end

    @model function m5_1t(A, D)
        a ~ Normal(0, 0.2)
        bA ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = lin(a, A[i], bA)
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_1t = sample(m5_1t(df.A, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 4)

    @model function m5_2t(M, D)
        a ~ Normal(0, 0.2)
        bM ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = lin(a, M[i], bM)
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_2t = sample(m5_2t(df.M, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 4)

    @model function m5_3t(A, M, D)
        a ~ Normal(0, 0.2)
        bA ~ Normal(0, 0.5)
        bM ~ Normal(0, 0.5)
        σ ~ Exponential(1)
        for i in eachindex(D)
            μ = a + M[i] * bM + A[i] * bA
            D[i] ~ Normal(μ, σ)
        end
    end

    chn5_3t = sample(m5_3t(df.A, df.M, df.D), NUTS(1000, 0.9), MCMCThreads(), 1000, 4)
    pw_lls5_1t = pointwise_log_likelihoods(m5_1t(df.A, df.D), chn5_1t)
    psis_loo_output5_1t = psis_loo(m5_1t(df.A, df.D), chn5_1t)
    psis_loo_output5_1t |> display

    pw_lls5_2t = pointwise_log_likelihoods(m5_2t(df.M, df.D), chn5_2t)
    psis_loo_output5_2t = psis_loo(m5_2t(df.M, df.D), chn5_2t)
    psis_loo_output5_2t |> display

    pw_lls5_3t = pointwise_log_likelihoods(m5_3t(df.A, df.M, df.D), chn5_3t)
    psis_loo_output5_3t = psis_loo(m5_3t(df.A, df.M, df.D), chn5_3t)
    psis_loo_output5_3t |> display

    n_tuple = (
        m5_1t=psis_loo_output5_1t,
        m5_2t=psis_loo_output5_2t,
        m5_3t=psis_loo_output5_3t,
    )
    comps = loo_compare(n_tuple)
    comps |> display

    @test comps.estimates(:m5_1t, :cv_est) ≈ 0.00 atol = .01
    @test comps.estimates(:m5_1t, :se_cv_est) ≈ 0.00 atol = .01
    @test comps.estimates(:m5_1t, :weight) ≈ 0.67 atol = .01
    @test comps.estimates(:m5_2t, :cv_est) ≈ -6.68 atol = .01
    @test comps.estimates(:m5_2t, :se_cv_est) ≈ 4.79 atol = .01
    @test comps.estimates(:m5_2t, :weight) ≈ 0.00 atol = .01
    @test comps.estimates(:m5_3t, :cv_est) ≈ -0.69 atol = .01
    @test comps.estimates(:m5_3t, :se_cv_est) ≈ 0.42 atol = .01
    @test comps.estimates(:m5_3t, :weight) ≈ 0.33 atol = .01
    @test sum(comps.estimates(:, :weight, :)) ≈ 1
    total = NamedDims.unname(sum(comps.pointwise(:, :cv_est, :); dims=:data))
    @test reshape(total, 3) ≈ comps.estimates(:, :cv_est) atol=.001

end


text:
const deck_list = [
  "ARC", "DEMO",
  "Steady", "Pulsed"
]

for deck in deck_list
  include("$(lowercase(deck)).jl")
  FusionSystems.@eval export $(Symbol(deck))
end


text:
using LinearAlgebra
using Nystrom
using StaticArrays
using Random
using Plots
plotlyjs()
Random.seed!(1)

## general definitions
xout = SVector(3,3)
ops = (
    Laplace(;dim=2),
    Helmholtz(;dim=2,k=1.2),
    Elastostatic(;dim=2,μ=2,λ=3)
)
methods = (
    Nystrom.assemble_dim,
    Nystrom.assemble_gk
)

pde      = ops[1]     # choose the pde
assemble = methods[1] # choose the singular integration method

T    = Nystrom.default_density_eltype(pde)
c    = rand(T)

## place to store dofs and error
ndofs = []
ee = []

npatches = [2,4,8,16,32,64,128,256,512]
qorder   =  5

for n in npatches
    ## generate mesh
    Geometry.clear_entities!()
    Ω   = ParametricSurfaces.Circle() |> Geometry.Domain
    Γ   = boundary(Ω)
    M   = meshgen(Γ,(n,))
    mesh = NystromMesh(view(M,Γ),order=qorder)
    push!(ndofs,length(Nystrom.dofs(mesh)))
    ## create exact solution
    u    = (qnode) -> SingleLayerKernel(pde)(xout,qnode)*c
    dudn = (qnode) -> transpose(DoubleLayerKernel(pde)(xout,qnode))*c
    γ₀u   = Density(u,mesh)
    γ₁u   = Density(dudn,mesh)
    γ₀u_norm = norm(norm.(γ₀u,Inf),Inf)
    # define integral operators
    Sop = SingleLayerOperator(pde,mesh)
    Dop = DoubleLayerOperator(pde,mesh)
    ## compute error
    S = assemble(Sop)
    D = assemble(Dop)
    e1    = Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)/γ₀u_norm
    push!(ee,norm(e1,Inf))
end

fig = plot(ndofs,ee,xscale=:log10,yscale=:log10,m=:o,label="error",lc=:black);
plot!(fig,xlabel="n",ylabel="error");
for p in 1:5
    c = ee[end]*ndofs[end]^p
    plot!(fig,ndofs,c ./ ndofs.^(p),xscale=:log10,yscale=:log10,label="h^$p",ls=:dash);
end
display(fig)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
#Tests des packages
include("MPCCProblems.jl")
include("../ALAS/MPCCPlot2D.jl")

#Est-ce que l'on veut les figures :
plot=false

r0=1.0
s0=0.5
#t0=sqrt(0.1) #t=o(r) normalement
t0=0.5
sig_r=0.1
sig_s=0.1
sig_t=0.05

mpcc=ex3 #ex1, ex2, ex3

@time xb,fb,orb,nb_eval = MPCCsolve.solve(mpcc,r0,sig_r,s0,sig_s,t0,sig_t)

try
 oa=orb.inner_output_alas[1]
 @printf("%s	",orb.solve_message)
end
@show nb_eval

nbc=mpcc.nb_comp;n=length(xb)-2*nbc

IterationsPlot2D(1,2,n+1,n+2,orb)


text:
"""
    steadystate_bicg(H, J, l=2; [log=false], kwargs...) -> rho[, log]

Compute the steady state density matrix of a Hamiltonian and a set of jump
operators by solving `L rho = 0` via the stabilized biconjugate gradient method.

# Arguments
* `H`: dense Hamiltonian.
* `J`: array of dense jump operators.
* `l`: number of GMRES steps per iteration.
* `kwargs...`: further arguments are passed on to the iterative solver.

See also: [`steadystate_bicg!`](@ref)
"""
function steadystate_bicg(H::O, J::Vector{O}, l::Int=2; log::Bool=false, kwargs...) where {B<:Basis,O<:AbstractOperator{B,B}}
    ρ0 = O<:DenseOperator ? DenseOperator(H.basis_l) : SparseOperator(H.basis_l)
    ρ0.data[1,1] = ComplexF64(1.0)
    return steadystate_iterative!(ρ0,H,J,:bicgstabl!,l;log=log,kwargs...)
end

"""
    steadystate_bicg!(rho0, H, J, l=2; [log=false], kwargs...) -> rho[, log]

Compute the steady state density matrix of a Hamiltonian and a set of jump
operators by solving `L rho = 0` via the stabilized biconjugate gradient method.

# Arguments
* `rho0`: Initial guess.
* `H`: dense Hamiltonian.
* `J`: array of dense jump operators.
* `l`: number of GMRES steps per iteration.
* `kwargs...`: further arguments are passed on to the iterative solver.

See also: [`steadystate_bicg`](@ref)
"""
steadystate_bicg!(ρ0::DenseOperator{B,B}, H::DenseOperator{B,B}, J::Vector{O}, l::Int=2; log::Bool=false, tol::Float64 = sqrt(eps(real(ComplexF64))), kwargs...) where {B<:Basis,O<:DenseOperator{B,B}} = steadystate_iterative!(ρ0,H,J,:bicgstabl!,l;log=log,kwargs...)
steadystate_bicg!(ρ0::SparseOperator{B,B}, H::SparseOperator{B,B}, J::Vector{O}, l::Int=2; log::Bool=false, tol::Float64 = sqrt(eps(real(ComplexF64))), kwargs...) where {B<:Basis,O<:SparseOperator{B,B}} = steadystate_iterative!(ρ0,H,J,:bicgstabl!,l;log=log,kwargs...)


text:
using Symbolics

include_model("quadrotor")

function fd(model::Quadrotor{Midpoint, FixedTime}, x, u, w, h, t)
    x + h[1] * f(model, x + 0.5 * h[1] * f(model, x, u, w), u, w)
end

n = model.n
m = model.m
d = model.d

@variables x_sym[1:n]
@variables u_sym[1:m]
@variables w_sym[1:d]
@variables h_sym[1:1]
@variables t_sym[1:1]
@variables P_sym[1:n, 1:n]
@variables p_sym[1:n]

fd_sym = fd(model, x_sym, u_sym, w_sym, h_sym, t_sym)

Q = 10.0 * Diagonal(ones(n))
R = 1.0 * Diagonal(ones(m))

function gt(x, u, t)
    transpose(x) * Q * x + transpose(u) * R * u
end

gt_sym = gt(x_sym, u_sym, t_sym)

function Vt(x, P, p)
    0.5 * transpose(x) * P * x + transpose(x) * p
end

V_sym = Vt(x_sym, P_sym, p_sym)

function Qt(x, u, P⁺, p⁺, h, t)
    x⁺ = fd(model, x, u, zeros(model.d), h, t)
    gt(x, u, t) + Vt(x⁺, P⁺, p⁺)
end

Qt_sym = Qt(x_sym, u_sym, P_sym, p_sym, h_sym, t_sym)

Qt_grad = Symbolics.gradient(Qt_sym, [x_sym; u_sym], simplify = true)
Qt_hess = Symbolics.hessian(Qt_sym, [x_sym; u_sym], simplify = true)
Qt_hess_sp = Symbolics.sparsehessian(Qt_sym, [x_sym; u_sym], simplify = true)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
struct IHT{T} <: AbstractSparseCoder
    x
    Φ
    Φᵀ
    iterations
    nonzeros::T

    absTol
    filter_domain::Symbol

    IHT(x, Φ, Φᵀ; iterations=nothing, nonzeros::T, absTol=1e2*eps(), filter_domain=:convolution) where {T<:Union{Integer,Tuple}} = new{T}(x, Φ, Φᵀ, iterations, nonzeros, absTol, filter_domain)
end

# argmin_{y} || x - Φy ||_2^2/2 s.t. ||y||_0 ≤ K
# Φ     : synthesis operator
# Φᵀ    : adjoint of Φ
#
function (iht::IHT)(y0; iterations=iht.iterations, isverbose::Bool=false)
    len = length(y0)
    yₖ = y0
    εx = Inf
    εy = Inf

    x̃ = iht.Φ(yₖ)
    for k = 1:iht.iterations
        yₖ₋₁ = yₖ

        yₖ = hardshrink(yₖ₋₁ + iht.Φᵀ(iht.x - x̃), iht.nonzeros)
        x̃ = iht.Φ(yₖ)

        εx = norm(iht.x - x̃)^2 / 2
        εy = norm(yₖ - yₖ₋₁)

        if isverbose
            println("number of Iterations $k: err = $εx, ||Δy|| = $εy.")
        end

        if εy <= iht.absTol
            break
        end
    end
    (yₖ, εx)
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
## Complex-to-complex transforms.

"""
    FFT()

Complex-to-complex FFT.

See also
[`AbstractFFTs.fft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.fft).
"""
struct FFT <: AbstractTransform end

"""
    FFT!()

In-place version of [`FFT`](@ref).

See also
[`AbstractFFTs.fft!`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.fft!).
"""
struct FFT! <: AbstractTransform end

"""
    BFFT()

Unnormalised backward complex-to-complex FFT.

Like `AbstractFFTs.bfft`, this transform is not normalised.
To obtain the inverse transform, divide the output by the length of the
transformed dimension.

See also
[`AbstractFFTs.bfft`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.bfft).
"""
struct BFFT <: AbstractTransform end

"""
    BFFT()

In-place version of [`BFFT`](@ref).

See also
[`AbstractFFTs.bfft!`](https://juliamath.github.io/AbstractFFTs.jl/stable/api/#AbstractFFTs.bfft!).
"""
struct BFFT! <: AbstractTransform end

const TransformC2C = Union{FFT, FFT!, BFFT, BFFT!}

length_output(::TransformC2C, length_in::Integer) = length_in
eltype_output(::TransformC2C, ::Type{Complex{T}}) where {T <: FFTReal} = Complex{T}
eltype_input(::TransformC2C, ::Type{T}) where {T <: FFTReal} = Complex{T}

plan(::FFT, args...; kwargs...) = FFTW.plan_fft(args...; kwargs...)
plan(::FFT!, args...; kwargs...) = FFTW.plan_fft!(args...; kwargs...)
plan(::BFFT, args...; kwargs...) = FFTW.plan_bfft(args...; kwargs...)
plan(::BFFT!, args...; kwargs...) = FFTW.plan_bfft!(args...; kwargs...)

binv(::FFT, d) = BFFT()
binv(::FFT!, d) = BFFT!()
binv(::BFFT, d) = FFT()
binv(::BFFT!, d) = FFT!()

is_inplace(::Union{FFT, BFFT}) = false
is_inplace(::Union{FFT!, BFFT!}) = true

_intprod(x::Int, y::Int...) = x * _intprod(y...)
_intprod() = one(Int)
_prod_dims(s::Dims, dims) = _intprod((s[i] for i in dims)...)
_prod_dims(A::AbstractArray, dims) = _prod_dims(size(A), dims)

scale_factor(::TransformC2C, A, dims) = _prod_dims(A, dims)


text:
const _attributes_string =
"""
Protected
Orderless
Flat
OneIdentity
Listable
Constant
NumericFunction
Locked
ReadProtected
HoldFirst
HoldRest
HoldAll
HoldAllComplete
NHoldFirst
NHoldRest
NHoldAll
SequenceHold
Temporary
Stub
"""

const possible_attributes = map(Symbol,split(_attributes_string))

# For Heads that are not symbols
get_attribute(args...) = false

# Return true if sj has attribute attr
function get_attribute(sj::SJSymbol, attr::Symbol)
    get(getssym(sj).attr,attr,false)
end

_set_attribute(s,attr) = getssym(Symbol(s)).attr[Symbol(attr)] = true

symattr(s::SJSymbol) = getssym(s).attr
get_attributes(sj::SymString) = ( ks = sort!(collect(Any, keys(symattr(Symbol(sj))))) )

## We will get rid of Qsym
set_attribute(sj::Qsym, attr::Symbol) = (getssym(sj).attr[attr] = true)

"""
    unset_attribute(sj::SJSymbol, attr::Symbol)

unset a single attribute `attr` for `sj`.
"""
unset_attribute(sj::SJSymbol, attr::Symbol) = delete!(getssym(sj).attr, attr)

"""
    clear_attributes(sja::SymString...)

clear all attributes set for the symbols `sja`.
"""
clear_attributes(sja::SymString...) =  foreach( x -> empty!(getssym(x).attr), sja)

## to be removed
clear_attributes(sj::Qsym) =  empty!(getssym(sj).attr)

### All access to data structures is above this line, or in mxpr.jl
### i.e. the interface is above this line

# Return true if head of mx has attribute attr
function get_attribute(mx::Mxpr{T}, attr::Symbol) where T
    get_attribute(T,attr)
end

## need a version that warns for users and asserts for developing ?
function set_attribute1(sj::SymString, attr::SymString)
    (ssj, sattr) = (Symbol(sj),Symbol(attr))
    if sattr in possible_attributes
        _set_attribute(sj,attr)
        return nothing
    end
    @symwarn("'$attr' is not an attribute")
    return nothing
end

"""
    set_attribute(s::SymString,attr::SymString)

set the `attr` attribute for `s`. Arguments may be `Symbol` or `String`.
"""
set_attribute(s::SymString,attr::SymString)= set_attribute1(Symbol(s),Symbol(attr))
set_attribute(sa::AbstractArray, attr::SymString) = foreach( x -> set_attribute1(x,attr), sa)
set_attribute(s::SymString, attra::AbstractArray) = foreach( x -> set_attribute1(s,x), attra)
set_attribute(sa::AbstractArray, attra::AbstractArray) = foreach( x -> set_attribute(sa,x), attra)

"""
    set_pattribute(s::SymString, attr::SymString)

set the attribute `attr` for symbol `s`. Also set the `Protected` attribute.
"""
function set_pattribute(s::SymString, attr::SymString)
    protect(s)
    set_attribute(s,attr)
end

"""
    set_pattribute(sa::AbstractArray, attr::SymString)

call `set_pattribute` for each of the elements in `sa`.
"""
set_pattribute(sa::AbstractArray, attr::SymString) = foreach( x -> set_pattribute(x,attr), sa)

"""
    set_pattribute(sa::AbstractArray, attr::SymString)

call `set_pattribute` for each of the elements in `attra`.
"""
function set_pattribute(s::SymString, attra::AbstractArray)
    protect(s)
    set_attribute(s,attra)
end

"""
    set_pattribute(sa::AbstractArray, attra::AbstractArray)

set all of the attributes in `attra` for all of the symbols in `sa`. Also
set the attribute `Protected` for each of the `sa`.
"""
function set_pattribute(sa::AbstractArray, attra::AbstractArray)
    foreach(protect, sa)
    set_attribute(sa,attra)
end

"""
    set_sysattributes(sym, attr)

register `sym` as a system symbol and set attributes. Attribute
`Protected` is set even it is not specified. Both `sym` and `attr`
can a strings or a symbol or arrays of these. In case both are arrays,
all symbols get all attributes.
"""
function set_sysattributes(sym, attr)
    register_system_symbol(sym)
    set_pattribute(sym,attr)
end

function set_sysattributes(syms::AbstractArray, attr)
    foreach(register_system_symbol, syms)
    set_pattribute(syms,attr)
end

function set_sysattributes(sym)
    register_system_symbol(sym)
    protect(sym)
end

function set_sysattributes(syms::AbstractArray)
    foreach(register_system_symbol, sym)
    foreach(protect,sym)
end

function protectedsymbols_strings()
    symstrings = Array{String}(undef, 0)
    for s in keys(CurrentContext.symtab)
        if get_attribute(s,:Protected) && s != :ans
            push!(symstrings,string(getsym(s))) end
    end
    sort!(symstrings)
end

function protectedsymbols()
    args = newargs()
    for s in keys(CurrentContext.symtab)
        if get_attribute(s,:Protected) && s != :ans
            push!(args,getsym(s)) end
    end
    mx = mxpra(:List, sort!(args))
end

# Related code in predicates.jl and attributes.jl
unprotect(sj::SJSym) = unset_attribute(sj,:Protected)

protect(sj::SymString) = set_attribute(sj,:Protected)

## hmmm, we could do this for all attributes
for s in (:HoldFirst,:HoldAll,:HoldRest,:HoldAllComplete, :SequenceHold, :Flat, :Orderless,
          :Listable, :Protected, :Constant)
    sf = Symbol("is",s)
    @eval ($sf)(x) = get_attribute(x,$(QuoteNode(s)))
end


text:
module BouncingSphere3D_2

using Modia
using ModiaLang.StaticArrays
@usingModiaPlot

BouncingSphere = Model3D(
    boxHeigth = 0.05,
    groundMaterial = VisualMaterial(color="DarkGreen", transparency=0.5),
    sphereMaterial = VisualMaterial(color="Red"),
    gravField = UniformGravityField(g=9.81, n=[0, -1, 0]),
    world = Object3D(feature=Scene(gravityField=:gravField,
                                   visualizeFrames=false, defaultFrameLength=0.2, gap=0.2,
                                   enableContactDetection=true, visualizeContactPoints=false)),
    worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
    ground = Object3D(parent=:world, translation=:[0.0, -boxHeigth/2, 0.0],
                      feature=Solid(shape=Box(lengthX=0.8, lengthY=:boxHeigth, lengthZ=0.5),
                                    visualMaterial=:groundMaterial, solidMaterial="Steel",
                                    collision=true)),
    bottom = Object3D(parent=:world, translation=:[0.0, -0.4-boxHeigth/2, 0.5],
                      feature=Solid(shape=Box(lengthX=0.8, lengthY=:boxHeigth, lengthZ=0.5),
                                    visualMaterial=:groundMaterial, solidMaterial="Steel",
                                    collision=true)),
    wall = Object3D(parent=:world, translation=:[0.0, -0.2, 0.75+boxHeigth/2],
                    feature=Solid(shape=Box(lengthX=0.8, lengthY=0.4, lengthZ=:boxHeigth),
                                  visualMaterial=:groundMaterial, solidMaterial="Steel",
                                  collision=true)),
    sphere = Object3D(feature=Solid(shape=Sphere(diameter=0.2),
                                    visualMaterial=:sphereMaterial, solidMaterial="Steel",
                                    collision=true)),
    free = FreeMotion(obj1=:world, obj2=:sphere, r=Var(init=SVector{3,Float64}(0.0, 1.0, 0.0)), w=Var(init=SVector{3,Float64}(10.0, 0.0, -5.0)))
)

bouncingSphere = @instantiateModel(BouncingSphere, unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 2.7
tolerance = 1e-8
requiredFinalStates = [0.29302652789657424, -1.2392524758992525, -0.13171402091799628, 0.11357140027502943, -4.417276610857313, -0.8729857092038694, 2.3810648360145152, -0.426312211044484 , 0.13656197405807574, -7.335181005993924, 0.812877179123953, 4.870682494334355]
simulate!(bouncingSphere, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true,
          requiredFinalStates_rtol = 0.2, requiredFinalStates_atol = 0.2, requiredFinalStates=requiredFinalStates)

plot(bouncingSphere, ["free.r" "free.rot"; "free.v" "free.w"], figure=1)

end