text:
""" Generic data structures for C-sets (presheaves).
"""
module CSets
export AbstractCSet, AbstractCSetType, CSet, CSetType,
  nparts, subpart, has_subpart, incident,
  add_part!, add_parts!, copy_parts!, set_subpart!, set_subparts!

using Compat
using LabelledArrays, StaticArrays

using ...Present
using ...Theories: Category, dom, codom

# C-set data types
##################

""" Abstract type for C-sets (presheaves).

The type parameters are:
- `Ob`: tuple of symbols naming objects
- `Hom`: tuple of symbols naming morphims
- `Dom`: tuple of integers giving domain of each morphism
- `Codom`: tuple of integers giving codomain of each morphism
- `Data`: tuple of symbols naming extra morphisms for data/attributes
- `DataDom`: tuple of integers giving domain of each data morphism

For example, a graph with one vertex attribute could be represented as
`AbstractCSet{(:V,:E),(:src,:tgt),(2,2),(1,1),(:vattr,),(1,)}`.

See also: [`CSet`](@ref) and [`CSetType`](@ref).
"""
abstract type AbstractCSet{Ob,Hom,Dom,Codom,Data,DataDom} end

""" Data type for C-sets (presheaves).

Instead of filling out the type parameters manually, we recommend using the
function [`CSetType`](@ref) to generate a `CSet` type from a presentation of a
category. Nevertheless, the first six type parameters are documented at
[`AbstractCSet`](@ref). The remaining type parameters are an implementation
detail and should be ignored.

Following LightGraphs.jl, the incidence vectors, stored in the `incidence`
field, are kept in sorted order. To ensure consistency, no field of the struct
should ever be mutated directly.
"""
mutable struct CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,NOb,NHom,NIndex} <:
       AbstractCSet{Ob,Hom,Dom,Codom,Data,DataDom}
  nparts::SLArray{Tuple{NOb},Int,1,NOb,Ob}
  subparts::SLArray{Tuple{NHom},Vector{Int},1,NHom,Hom}
  incident::SLArray{Tuple{NIndex},Vector{Vector{Int}},1,NIndex,Index}
  data::NamedTuple{Data}
end

function CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index}(; kw...) where
    {Ob,Hom,Dom,Codom,Data,DataDom,Index}
  CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index}((; kw...))
end
function CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index}(
    datatypes::NamedTuple{Data}) where {Ob,Hom,Dom,Codom,Data,DataDom,Index}
  NOb, NHom, NIndex = length(Ob), length(Hom), length(Index)
  @assert length(Dom) == NHom && length(Codom) == NHom
  @assert length(DataDom) == length(Data)
  CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,NOb,NHom,NIndex}(datatypes)
end
function CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,NOb,NHom,NIndex}(
    datatypes::NamedTuple{Data}) where {Ob,Hom,Dom,Codom,Data,DataDom,Index,NOb,NHom,NIndex}
  CSet{Ob,Hom,Dom,Codom,Data,DataDom,Index,NOb,NHom,NIndex}(
    SLArray{Tuple{NOb},Ob}(zeros(SVector{NOb,Int})),
    SLArray{Tuple{NHom},Hom}(Tuple(Int[] for i in 1:NHom)),
    SLArray{Tuple{NIndex},Index}(Tuple(Vector{Int}[] for i in 1:NIndex)),
    NamedTuple{Data}(T[] for T in datatypes))
end

""" Generate an abstract C-set type from a presentation of a category.

See also: [`CSetType`](@ref).
"""
function AbstractCSetType(pres::Presentation{Category}; data=())
  # Only include `Data` and `DataDom` params if there are data morphisms.
  Ts = CSetTypeParams(pres, data=data)
  AbstractCSet{(isempty(Ts[5]) ? Ts[1:4] : Ts[1:6])...}
end

""" Generate a C-set data type from a presentation of a category.

See also: [`AbstractCSetType`](@ref).
"""
function CSetType(pres::Presentation{Category}; data=(), index=())
  CSet{CSetTypeParams(pres, data=data, index=index)...}
end
function CSetTypeParams(pres::Presentation{Category}; data=(), index=())
  obs, homs = generators(pres, :Ob), generators(pres, :Hom)
  data_obs, obs = separate(ob -> nameof(ob) ∈ data, obs)
  data_homs, homs = separate(hom -> codom(hom) ∈ data_obs, homs)
  ob_index = ob -> findfirst(obs .== ob)::Int
  (Tuple(nameof.(obs)), Tuple(nameof.(homs)),
   Tuple(@. ob_index(dom(homs))), Tuple(@. ob_index(codom(homs))),
   Tuple(nameof.(data_homs)), Tuple(@. ob_index(dom(data_homs))), Tuple(index))
end
separate(f, a::AbstractArray) = (i = f.(a); (a[i], a[.!i]))

function Base.:(==)(x1::T, x2::T) where T <: CSet
  # The incidence data is redundant, so need not be compared.
  x1.nparts == x2.nparts && x1.subparts == x2.subparts && x1.data == x2.data
end

function Base.copy(cset::T) where T <: CSet
  T(cset.nparts, map(copy, cset.subparts),
    map(copy, cset.incident), map(copy, cset.data))
end

Base.empty(cset::T) where T <: CSet = T(map(eltype, cset.data))

# C-set interface
#################

""" Number of parts of given type in a C-set.
"""
nparts(cset::CSet, type::Symbol) = cset.nparts[type]

""" Get subpart of part in C-set.

Both single and vectorized access are supported.
"""
subpart(cset::CSet, part, name::Symbol) = subpart(cset, name)[part]
subpart(cset::CSet, name::Symbol) = _subpart(cset, Val(name))

@generated function _subpart(cset::T, ::Val{name}) where
    {name, obs,homs,doms,codoms,data, T <: CSet{obs,homs,doms,codoms,data}}
  if name ∈ homs
    :(cset.subparts.$name)
  elseif name ∈ data
    :(cset.data.$name)
  else
    throw(KeyError(name))
  end
end

""" Whether a C-set has a subpart with the given name.
"""
has_subpart(cset::CSet, name::Symbol) = _has_subpart(cset, Val(name))

@generated function _has_subpart(cset::T, ::Val{name}) where
    {name, obs,homs,doms,codoms,data, T<: CSet{obs,homs,doms,codoms,data}}
  name ∈ homs || name ∈ data
end

""" Get superparts incident to part in C-set.
"""
incident(cset::CSet, part, name::Symbol) = cset.incident[name][part]

""" Add part of given type to C-set, optionally setting its subparts.

Returns the ID of the added part.

See also: [`add_parts!`](@ref).
"""
add_part!(cset::CSet, type::Symbol) = only(_add_parts!(cset, Val(type), 1))

function add_part!(cset::CSet, type::Symbol, subparts)
  part = add_part!(cset, type)
  set_subparts!(cset, part, subparts)
  part
end

""" Add parts of given type to C-set, optionally setting their subparts.

Returns the range of IDs for the added parts.

See also: [`add_part!`](@ref).
"""
add_parts!(cset::CSet, type::Symbol, n::Int) = _add_parts!(cset, Val(type), n)

function add_parts!(cset::CSet, type::Symbol, n::Int, subparts)
  parts = add_parts!(cset, type, n)
  set_subparts!(cset, parts, subparts)
  parts
end

@generated function _add_parts!(cset::T, ::Val{type}, n::Int) where
    {type, obs,homs,doms,codoms,data,data_doms,indexed,
     T <: CSet{obs,homs,doms,codoms,data,data_doms,indexed}}
  obnum = findfirst(obs .== type)::Int
  in_homs = [ homs[i] for (i, codom) in enumerate(codoms) if codom == obnum ]
  out_homs = [ homs[i] for (i, dom) in enumerate(doms) if dom == obnum ]
  data_homs = [ data[i] for (i, dom) in enumerate(data_doms) if dom == obnum ]
  indexed_homs = filter(hom -> hom ∈ indexed, in_homs)
  # TODO: The three loops could (should?) be unrolled. Or is Julia's compiler
  # smart enough to do this on its own?
  quote
    if n == 0; return Int[] end
    nparts = cset.nparts.$type + n
    cset.nparts = SLVector(cset.nparts; $type=nparts)
    start = nparts - n + 1
    for name in $(Tuple(out_homs))
      sub = cset.subparts[name]
      resize!(sub, nparts)
      @inbounds sub[start:nparts] .= 0
    end
    for name in $(Tuple(indexed_homs))
      incident = cset.incident[name]
      resize!(incident, nparts)
      @inbounds for part in start:nparts; incident[part] = Int[] end
    end
    for name in $(Tuple(data_homs))
      resize!(cset.data[name], nparts)
    end
    start:nparts
  end
end

""" Copy parts from one C-set to another.

All subparts between the selected parts, including any attached data, are
preserved. Thus, if the selected parts form a sub-C-set, then the whole
sub-C-set is preserved. On the other hand, if the selected parts do *not* form a
sub-C-set, then some copied parts will have undefined subparts.
"""
copy_parts!(cset::CSet, from::CSet; kw...) = copy_parts!(cset, from, (; kw...))

copy_parts!(cset::CSet, from::CSet, type::Symbol, parts) =
  copy_parts!(cset, from, (; type => parts))

@generated function copy_parts!(cset::T, from::T, parts::NamedTuple{types}) where
    {types, obs,homs,doms,codoms, T <: CSet{obs,homs,doms,codoms}}
  obnums = [ findfirst(obs .== type)::Int for type in types ]
  in_obs, out_homs = Symbol[], Tuple{Symbol,Symbol,Symbol}[]
  for (hom, dom, codom) in zip(homs, doms, codoms)
    if dom ∈ obnums && codom ∈ obnums
      push!(in_obs, obs[codom])
      push!(out_homs, (hom, obs[dom], obs[codom]))
    end
  end
  in_obs = Tuple(unique!(in_obs))
  quote
    newparts = NamedTuple{$types}(tuple($(map(types) do type
      :(_copy_parts_data!(cset, from, Val($(QuoteNode(type))), parts.$type))
    end...)))
    partmaps = NamedTuple{$in_obs}(tuple($(map(in_obs) do type
      :(Dict{Int,Int}(zip(parts.$type, newparts.$type)))
    end...)))
    for (name, dom, codom) in $(Tuple(out_homs))
      for (p, newp) in zip(parts[dom], newparts[dom])
        q = subpart(from, p, name)
        newq = get(partmaps[codom], q, nothing)
        if !isnothing(newq)
          set_subpart!(cset, newp, name, newq)
        end
      end
    end
    newparts
  end
end

""" Copy parts of a single type from one C-set to another.

Any data attached to the parts is also copied.
"""
@generated function _copy_parts_data!(cset::T, from::T, ::Val{type}, parts) where
    {type, obs,homs,doms,codoms,data,data_doms,
     T <: CSet{obs,homs,doms,codoms,data,data_doms}}
  obnum = findfirst(obs .== type)::Int
  data_homs = [ data[i] for (i, dom) in enumerate(data_doms) if dom == obnum ]
  quote
    newparts = add_parts!(cset, $(QuoteNode(type)), length(parts))
    for name in $(Tuple(data_homs))
      cset.data[name][newparts] .= from.data[name][parts]
    end
    newparts
  end
end

""" Mutate subpart of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subparts!`](@ref).
"""
set_subpart!(cset::CSet, part::Int, name::Symbol, subpart) =
  _set_subpart!(cset, part, Val(name), subpart)

function set_subpart!(cset::CSet, part::AbstractVector{Int},
                      name::Symbol, subpart)
  broadcast(part, subpart) do part, subpart
    _set_subpart!(cset, part, Val(name), subpart)
  end
end

set_subpart!(cset::CSet, ::Colon, name::Symbol, subpart) =
  set_subpart!(cset, name, subpart)
set_subpart!(cset::CSet, name::Symbol, new_subpart) =
  set_subpart!(cset, 1:length(subpart(cset, name)), name, new_subpart)

@generated function _set_subpart!(cset::T, part::Int, ::Val{name}, subpart) where
    {name, obs,homs,doms,codoms,data,data_doms,indexed,
     T <: CSet{obs,homs,doms,codoms,data,data_doms,indexed}}
  if name ∈ indexed
    quote
      old = cset.subparts.$name[part]
      cset.subparts.$name[part] = subpart
      if old > 0
        deletesorted!(cset.incident.$name[old], part)
      end
      if subpart > 0
        insertsorted!(cset.incident.$name[subpart], part)
      end
    end
  elseif name ∈ homs
    :(cset.subparts.$name[part] = subpart)
  elseif name ∈ data
    :(cset.data.$name[part] = subpart)
  else
    throw(KeyError(name))
  end
end

""" Mutate subparts of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subpart!`](@ref).
"""
function set_subparts!(cset::CSet, part, subparts)
  for (name, subpart) in pairs(subparts)
    set_subpart!(cset, part, name, subpart)
  end
end

""" Insert into sorted vector, preserving the sorting.
"""
function insertsorted!(a::AbstractVector, x)
  insert!(a, searchsortedfirst(a, x), x)
end

""" Delete one occurrence of value from sorted vector.
"""
function deletesorted!(a::AbstractVector, x)
  i = searchsortedfirst(a, x)
  @assert i <= length(a) && a[i] == x "Value $x is not contained in $a"
  deleteat!(a, i)
end

end


text:
import StaticArrays.SVector
import StaticArrays.norm

abstract type AbstractDOF end

struct Point{d, T}<:AbstractDOF
    pos::SVector{d,T}
    id::Int
end

function distance(p1::Point{d,T}, p2::Point{d,T}) where {d, T}
    # return norm(p1.pos - p2.pos)
    out = zero(T)
    @fastmath @inbounds @simd for k = 1 : d
        out += (p1.pos[k] - p2.pos[k])^2
    end
    return sqrt(out)
end

# function that takes in a matrix and returns a n array of points corresponding to dim-slices of the matrix.
function mat2points(x::AbstractMatrix; dims::Int=1)
    @assert dims ∈ [1,2]
    if dims == 2
        x = x'
    end

    points = Vector{Point{size(x, 1), eltype(x)}}(undef, size(x, 2))
    for k = 1 : size(x, 2)
        points[k] = Point{size(x, 1), eltype(x)}(SVector{size(x, 1), eltype(x)}(x[:, k]), k)
    end

    return points
end

ext: jl


text:
## order processing (generic, not stargety-specific)

"Generate oder ID string for a new order"
function neworderid(trig::String)
  ### TODO hash order id to encrypt trigger-string
  return string(unix2datetime(time())) * "_" * string(trig)
end

"""
Prepare new order from `targ` (`(poschg::Int64, Vector[limitprice,stopprice]`)
and trigger-string `trig`.
Note: this function prepares limit and market orders for submission.
Stop-part of stoplimit orders is handled at the software level
in `orderhandling!` (even for live trading),
which calls `targ2order!` for limit order submission
if stop-price of stoplimit order is reached.
...
Overwrites `orde` and returns `Bool` request status.
"""
function targ2order!(orde::Order,
                     targ::Targ,
                     trig::String,
                     position_actual::Int64,
                     backtest::Bool)
  # this function should be called when position change is needed
  if targ[1] == 0
    warn("targ2order! should not be called when no position change is needed")
    return false
  end

  # overwrite orde components
  if length(targ[2]) == 0
    # no stop/limit conditions
    orde.price = NaN
    orde.ordertype = :market
  elseif length(targ[2]) == 1
    # limit order
    orde.price = targ[2][1]
    orde.ordertype = :limit
  else
    warn("targ2order! does not handle $targ")
    return false
  end
  orde.id = neworderid(trig)
  orde.quantity = abs(targ[1])
  orde.side = targ[1] > 0 ? :buy : :sell
  orde.status = :pending
  orde.trigger = trig

  # submit the order
  if backtest
    respid = submit_ordernew(orde, backtest)
  else
    respid = submit_ordernew(orde, position_actual)
  end
  if respid == "FAIL"
    return false
  end
  if orde.id != respid
    orde.id = respid
  end
  return true
end

"""
Order handling for backtesting and live trading.
Input:
- target `targ` as `(poschg::Int64, Vector[limitprice, stopprice]`;
- current/instantaneous price `pnow`
- current time `tnow`; for backtest, the time corresponding to `targ`
(i.e. the current OHLC step/bar time).

In-place modifies:

* `position_actual_mut` vector;
* `ordcurr` object;
* `backtestblotter` associative collection.

Returns tuple with:

* `Bool` system status;
* `Float64` current cumulative profit/loss.

NOTE: As opposed to `tradeperf` function, here total PnL is updated
at each price change time-point.
"""
function orderhandling!(targ::Targ,
                        pnow::Float64, tnow::DateTime,
                        position_actual_mut::Vector{Int64},
                        ordcurr::Order,
                        blotter::Blotter,
                        backtest::Bool
                        )
  posact = position_actual_mut[1]

  #println("enter orderhandling!: posact = $posact")
  #println(targ)

  if backtest
    if ispending(ordcurr)
      # check if previously pending order is now complete
      resp = query_orderstatus(ordcurr, pnow)
      if resp == :complete
        # assume order filled in full quantity at pnow
        poschg = getorderposchg(ordcurr)
        position_actual_mut[1] = posact + poschg
        blotter[tnow] = (poschg, pnow)
        # record different order status
        ordcurr.status = resp
        stillpending = false
        # need to re-evaluate target now that position changed
        return true, tradepnlfinal(blotter, pnow)
      elseif resp == :pending
        stillpending = true
      else
        # neither pending nor complete
        ordcurr.status = resp
        stillpending = false
      end
    else
      stillpending = false
    end
  else
    ### TODO query recent order fills
    ### update actual position based on that and rework pending order logic below
    error("yet to be implemented")

    ############################################################################
#     if ispending(ordcurr)
#     # check if previously pending order is now complete
#     if backtest
#       resp = query_orderstatus(ordcurr, pnow)
#     else
#       resp = query_orderstatus(ordcurr.id)
#     end
#     if resp == :complete
#       ### TODO partial order handling
#       ### (e.g. part of limit order execuled, part still pending)
#       ###### HIGHEST priority todo ##################################
#       ### get recent order fills and update position based on that
#       ### do that earlier in this function

#       # update actual position
#       ### for now assume order executed in full quantity
#       ### later: make opentrs length-10 or so
#       #opentrs
#       #position_actual = posact + getorderposchg(ordcurr)
#       position_actual_mut[1] = posact + getorderposchg(ordcurr)
#     #need to re-evaluate target now that position changed
#         return true, tradepnlfinal(blotter, pnow)
#       ### add new row to blotter here ###############################

#       # record different order status
#       ordcurr.status = resp
#       stillpending = false
#     elseif resp == :pending
#       stillpending = true
#     else
#       # neither pending nor complete
#       ordcurr.status = resp
#       stillpending = false
#     end
#   else
#     stillpending = false
#   end
    ##########################################################################
  end # if backtest

  # check if any position change is targeted at all
  if targ[1] == 0
    # holding the current position
    if stillpending
      # cancel pending order
      if backtest
        respstatus = submit_ordercancel(ordcurr)
      else
        respstatus = submit_ordercancel(ordcurr.id)
      end
      if respstatus
        setcancelled!(ordcurr)
        return respstatus, tradepnlfinal(blotter, pnow)
      else
        return respstatus, tradepnlfinal(blotter, pnow)
      end
    end
    # no pending orders, hold current position
    return true, tradepnlfinal(blotter, pnow)
  end

  # pending order (if any) must be in line with the target
  if stillpending
    if getorderposchg(ordcurr) == targ[1]
      if length(targ[2]) < 1 && ordcurr.ordertype == :market
        # let the order continue as is
        return true, tradepnlfinal(blotter, pnow)
      end
      if length(targ[2]) == 1 && ordcurr.ordertype == :limit
        # check if limit prices match
        if isapprox(targ[2][1], ordcurr.price)
          # let the order continue as is
          return true, tradepnlfinal(blotter, pnow)
        end
      end
      # otherwise no complete match between pending and target
      # cancel pending order not matching the target
      if backtest
        respstatus = submit_ordercancel(ordcurr)
      else
        respstatus = submit_ordercancel(ordcurr.id)
      end
      if respstatus
        setcancelled!(ordcurr)
        return respstatus, tradepnlfinal(blotter, pnow)
      else
        return respstatus, tradepnlfinal(blotter, pnow)
      end
    else
      # cancel pending order not in line with the target
      if backtest
        respstatus = submit_ordercancel(ordcurr)
      else
        respstatus = submit_ordercancel(ordcurr.id)
      end
      if respstatus
        setcancelled!(ordcurr)
        return respstatus, tradepnlfinal(blotter, pnow)
      else
        return respstatus, tradepnlfinal(blotter, pnow)
      end
    end
  end

  # no pending orders at this point
  if length(targ[2]) == 2
    # check stop condition of the stoplimit order
    if targ[1] > 0
      # buy: instantaneous price must exceed the stop-price
      if pnow >= targ[2][2]
        return (targ2order!(ordcurr,
                           (targ[1], targ[2][1:1]),
                           "buy-side software-stop",
                           posact, backtest), tradepnlfinal(blotter, pnow))
      end
    else
      # sell: instantaneous price must be below the stop-price
      if pnow <= targ[2][2]
        return (targ2order!(ordcurr,
                           (targ[1], targ[2][1:1]),
                           "sell-side software-stop",
                           posact, backtest), tradepnlfinal(blotter, pnow))
      end
    end
  elseif length(targ[2]) < 2
    # limit or market order
    return (targ2order!(ordcurr, targ, "limitormarket", posact, backtest),
            tradepnlfinal(blotter, pnow))
  else
    warn("orderhandling! does not handle $targ")
    return false, tradepnlfinal(blotter, pnow)
  end

  return true, tradepnlfinal(blotter, pnow)
end


text:
using Debug
@debug function test()
    parts = {}
    @bp
    for j=1:3
        for i=1:3
            push!(parts,"($i,$j) ")
        end
    end
    @bp
    println(parts...)
end

test()


text:
"""
    Trixi

**Trixi.jl** is a numerical simulation framework for hyperbolic conservation
laws. A key objective for the framework is to be useful to both scientists
and students. Therefore, next to having an extensible design with a fast
implementation, Trixi is focused on being easy to use for new or inexperienced
users, including the installation and postprocessing procedures.

To get started, run your first simulation with Trixi using

    trixi_include(default_example())

See also: [trixi-framework/Trixi.jl](https://github.com/trixi-framework/Trixi.jl)
"""
module Trixi

# Include other packages that are used in Trixi
# (standard library packages first, other packages next, all of them sorted alphabetically)

using LinearAlgebra: diag, dot, mul!, norm, cross, normalize
using Printf: @printf, @sprintf, println

# import @reexport now to make it available for further imports/exports
using Reexport: @reexport

import DiffEqBase: CallbackSet, DiscreteCallback,
                   ODEProblem, ODESolution, ODEFunction,
                   get_du, get_tmp_cache, u_modified!,
                   get_proposed_dt, set_proposed_dt!, terminate!, remake
using CodeTracking: code_string
@reexport using EllipsisNotation # ..
import ForwardDiff
using HDF5: h5open, attributes
using LinearMaps: LinearMap
using LoopVectorization: LoopVectorization, @turbo, indices
using LoopVectorization.ArrayInterface: static_length
import MPI
using Polyester: @batch # You know, the cheapest threads you can find...
using OffsetArrays: OffsetArray, OffsetVector
using P4est
using RecipesBase
using Requires
@reexport using StaticArrays: SVector
using StaticArrays: MVector, MArray, SMatrix
using StrideArrays: PtrArray, StrideArray, StaticInt
using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_timer!
@reexport using UnPack: @unpack
using UnPack: @pack!

# finite difference SBP operators
using SummationByPartsOperators: AbstractDerivativeOperator, DerivativeOperator, grid
import SummationByPartsOperators: integrate, left_boundary_weight, right_boundary_weight
@reexport using SummationByPartsOperators:
  SummationByPartsOperators, derivative_operator

# Define the entry points of our type hierarchy, e.g.
#     AbstractEquations, AbstractSemidiscretization etc.
# Placing them here allows us to make use of them for dispatch even for
# other stuff defined very early in our include pipeline, e.g.
#     IndicatorLöhner(semi::AbstractSemidiscretization)
include("basic_types.jl")

# Include all top-level source files
include("auxiliary/auxiliary.jl")
include("auxiliary/mpi.jl")
include("equations/equations.jl")
include("mesh/mesh.jl")
include("solvers/solvers.jl")
include("semidiscretization/semidiscretization.jl")
include("semidiscretization/semidiscretization_hyperbolic.jl")
include("callbacks_step/callbacks_step.jl")
include("callbacks_stage/callbacks_stage.jl")
include("semidiscretization/semidiscretization_euler_gravity.jl")
include("time_integration/time_integration.jl")

# `trixi_include` and special elixirs such as `convergence_test`
include("auxiliary/special_elixirs.jl")

# Plot recipes and conversion functions to visualize results with Plots.jl
include("visualization/visualization.jl")

# export types/functions that define the public API of Trixi

export AcousticPerturbationEquations2D,
       CompressibleEulerEquations1D, CompressibleEulerEquations2D, CompressibleEulerEquations3D,
       CompressibleEulerMulticomponentEquations1D, CompressibleEulerMulticomponentEquations2D,
       IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
       IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D,
       HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
       LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
       InviscidBurgersEquation1D,
       LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D

export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_godunov,
       flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang,
       flux_kennedy_gruber, flux_shima_etal, flux_ec,
       FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
       FluxLaxFriedrichs, max_abs_speed_naive,
       FluxHLL, min_max_speed_naive,
       FluxRotated

export initial_condition_constant,
       initial_condition_gauss,
       initial_condition_density_wave, initial_condition_density_pulse,
       initial_condition_isentropic_vortex,
       initial_condition_khi,
       initial_condition_weak_blast_wave, initial_condition_blast_wave,
       initial_condition_sedov_blast_wave, initial_condition_medium_sedov_blast_wave,
       initial_condition_two_interacting_blast_waves, boundary_condition_two_interacting_blast_waves,
       initial_condition_blob,
       initial_condition_orszag_tang,
       initial_condition_rotor,
       initial_condition_shock_bubble,
       initial_condition_taylor_green_vortex

export boundary_condition_periodic,
       BoundaryConditionDirichlet,
       boundary_condition_wall_noslip,
       boundary_condition_wall,
       boundary_condition_zero,
       BoundaryConditionWall,
       boundary_state_slip_wall

export initial_condition_convergence_test, source_terms_convergence_test
export initial_condition_harmonic_nonperiodic, source_terms_harmonic
export initial_condition_poisson_periodic, source_terms_poisson_periodic
export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, boundary_condition_poisson_nonperiodic
export initial_condition_briowu_shock_tube, initial_condition_torrilhon_shock_tube, initial_condition_ryujones_shock_tube,
       initial_condition_shu_osher_shock_tube, initial_condition_shu_osher_shock_tube_flipped
export initial_condition_sedov_self_gravity, boundary_condition_sedov_self_gravity
export initial_condition_eoc_test_coupled_euler_gravity, source_terms_eoc_test_coupled_euler_gravity, source_terms_eoc_test_euler
export initial_condition_lid_driven_cavity, boundary_condition_lid_driven_cavity
export initial_condition_couette_steady, initial_condition_couette_unsteady, boundary_condition_couette
export initial_condition_gauss_wall
export initial_condition_monopole, boundary_condition_monopole

export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
       cons2entropy, entropy2cons
export density, pressure, density_pressure, velocity
export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity

export TreeMesh, CurvedMesh, UnstructuredQuadMesh, P4estMesh

export DG,
       DGSEM, LobattoLegendreBasis,
       VolumeIntegralWeakForm, VolumeIntegralStrongForm,
       VolumeIntegralFluxDifferencing,
       VolumeIntegralPureLGLFiniteVolume,
       VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
       SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
       MortarL2

export nelements, nnodes, nvariables,
       eachelement, eachnode, eachvariable

export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate

export SemidiscretizationEulerGravity, ParametersEulerGravity,
       timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!

export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
       SaveRestartCallback, SaveSolutionCallback, TimeSeriesCallback, VisualizationCallback,
       AMRCallback, StepsizeCallback,
       GlmSpeedCallback, LBMCollisionCallback,
       TrivialCallback

export load_mesh, load_time

export ControllerThreeLevel, ControllerThreeLevelCombined,
       IndicatorLöhner, IndicatorLoehner, IndicatorMax

export PositivityPreservingLimiterZhangShu

export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured

export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure

# Visualization-related exports
export PlotData1D, PlotData2D, getmesh, adapt_to_mesh_level!, adapt_to_mesh_level

function __init__()
  init_mpi()

  init_p4est()

  # Enable features that depend on the availability of the Plots package
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    using .Plots: plot, plot!, savefig
  end
end

include("auxiliary/precompile.jl")
_precompile_manual_()

end


text:
using GraphQLGen

# Generate Julia code files in the same directory based on the GraphQL schema in "schema.graphql".
#
# The code fles are
#  * graphqlgen_types.jl - contains all the GraphQL types
#  * graphqlgen_functions.jl - contains all the GraphQL functions (mutations, queries, subscriptions)
GraphQLGen.generate("$(@__DIR__)", "schema.graphql")


text:
"Represents a 1D mesh with vertices and elements"
struct Mesh1D
    nvertices::Int   # Number of vertices
    vertex_coordinates::Array{Float64, 1}  # Coordinate of each vertex
    nelements::Int  # Number of elements

    # Map element ID to vertices
    element_vertices::Array{Int, 2}

    "Create a 1D mesh.
        vertex_coordinates is a 1D array of node locations
        elements is a 2D array, mapping elements to vertex indices"
    function Mesh1D(vertex_coordinates, elements)
        # Check that elements are in range
        @assert maximum(elements) == length(vertex_coordinates)
        @assert minimum(elements) == 1

        new(length(vertex_coordinates),
            vertex_coordinates,
            size(elements)[1],
            elements)
    end
end

"Generate simple equidistant grid with K elements"
function MeshGen1D(xmin, xmax, K)

    # Number of vertices
    Nv = K+1

    # Generate node coordinates
    VX = zeros(Nv)
    for i = 1:Nv
        VX[i] = (xmax-xmin)*(i-1)/(Nv-1) + xmin
    end

    # read element to node connectivity
    EToV = zeros(Int, K, 2)
    for k = 1:K
        EToV[k,1] = k
        EToV[k,2] = k+1
    end

    return Mesh1D(VX, EToV)
end


text:
export github_raw

"""
    github_raw(::URI)
Given a GitHub repository URI, return a URI for raw file
access on `master`."
"""
function github_raw(uri::URI)::URI
    @assert uri.host == "github.com"
    split = URIs.splitpath(uri)
    @assert length(split) == 2
    user = split[1]
    repo = split[2]
    URI(;
        scheme="https",
        host="raw.githubusercontent.com",
        path="/$user/$repo/master")
end

max_stars_repo_name: MarkNahabedian/WebBasedWorkspace

max_issues_repo_name: MarkNahabedian/WebBasedWorkspace

max_forks_repo_name: MarkNahabedian/WebBasedWorkspace


text:
VectorField{T,N} = SVector{3,Array{T,N}}
VectorFieldF{N} = VectorField{Float,N}

new_vectorfield(shape) = SVector(ntuple(_ -> zeros(shape), 3))

mutable struct Fields{N}
    grid::Grid{N}
    E::VectorFieldF{N} # Defined in the middle of cell edges at half time steps
    B::VectorFieldF{N} # Defined in the middle of cell faces at whole time steps
    curl_E::VectorFieldF{N} # Defined in the same place as B at whole time steps
    curl_B::VectorFieldF{N} # Defined in the same place as E at half time steps
    j::VectorFieldF{N} # Defined in the same place as E at whole time steps

    function Fields(grid::Grid{N}) where {N}
        E = new_vectorfield(size(grid))
        B = new_vectorfield(size(grid))
        curl_E = new_vectorfield(size(grid))
        curl_B = new_vectorfield(size(grid))
        j = new_vectorfield(size(grid))
        new{N}(grid, E, B, curl_E, curl_B, j)
    end
end

function advance_E!(fields::Fields{N}, Δt) where {N}
    for n = 1:N
        fields.E[n] .+=
            (Δt * SPEED_OF_LIGHT²) .* fields.curl_B[n] .-
            (Δt * VACUUM_PERMITTIVITY⁻¹) .* fields.j[n]
    end
end

function advance_B!(fields::Fields{N}, Δt) where {N}
    for n = 1:N
        fields.B[n] .-= Δt .* fields.curl_E[n]
    end
end

function update_curl_E!(fields::Fields)
    curl!(fields.curl_E, fields.E, fields.grid, Up)
end

function update_curl_B!(fields::Fields)
    curl!(fields.curl_B, fields.B, fields.grid, Down)
end

abstract type FieldType end
struct ElectricField <: FieldType end
struct MagneticField <: FieldType end

struct FieldInterpolator{N}
    interpolator::Interpolator
    grid::Grid{N}
    coordinate_references::NTuple{N,Type{<:ReferenceCoords}}

    function FieldInterpolator(
        fields::Fields{N},
        ::Type{MagneticField},
        ::Dim{D},
    ) where {N,D}
        interpolator = Interpolator(fields.B[D])
        coordinate_references = ntuple(n -> (n == D) ? LowerEdges : Centers, N)
        new{N}(interpolator, fields.grid, coordinate_references)
    end

    function FieldInterpolator(
        fields::Fields{N},
        ::Type{ElectricField},
        ::Dim{D},
    ) where {N,D}
        interpolator = Interpolator(fields.E[D])
        coordinate_references = ntuple(n -> (n == D) ? Centers : LowerEdges, N)
        new{N}(interpolator, fields.grid, coordinate_references)
    end
end

function interpolate(field_interpolator::FieldInterpolator, position)
    continuous_indices = find_indices_at(
        field_interpolator.grid,
        position,
        field_interpolator.coordinate_references,
        index_type = Continuous,
    )
    interpolate(field_interpolator.interpolator, continuous_indices)
end

struct VectorFieldInterpolator{N}
    interpolators::NTuple{3,FieldInterpolator{N}}

    function VectorFieldInterpolator(
        fields::Fields{N},
        field_type::Type{<:FieldType},
    ) where {N}
        interpolators =
            ntuple(n -> FieldInterpolator(fields, field_type, Dim{n}()), 3)
        new{N}(interpolators)
    end
end

function interpolate(
    vectorfield_interpolator::VectorFieldInterpolator,
    position,
)
    [
        interpolate(vectorfield_interpolator.interpolators[n], position)
        for n = 1:3
    ]
end

function curl!(
    dest::VectorField{<:Any,3},
    field::VectorField{<:Any,3},
    grid::Grid{3},
    direction::Type{<:SliceShift},
)
    @dbgasserteq(size(dest[1]), size(grid))
    @dbgassert dest !== field
    for n = 1:3
        dest[n][interior(grid)...] .= curlcomp(field, grid, direction, Dim{n}())
    end
end

function curl!(
    dest::VectorField{<:Any,2},
    field::VectorField{<:Any,2},
    grid::Grid{2},
    direction::Type{<:SliceShift},
)
    @dbgasserteq(size(dest[1]), size(grid))
    @dbgassert dest !== field
    dest[3][interior(grid)...] .= curlcomp(field, grid, direction, DIM3)
end

function curlcomp(
    field::VectorField{<:Any,3},
    grid::Grid{3},
    direction::Type{<:SliceShift},
    ::Dim{1},
)
    shifteddiff(field[3], grid, direction, DIM2) .* cell_extents⁻¹(grid)[2] .-
    shifteddiff(field[2], grid, direction, DIM3) .* cell_extents⁻¹(grid)[3]
end

function curlcomp(
    field::VectorField{<:Any,3},
    grid::Grid{3},
    direction::Type{<:SliceShift},
    ::Dim{2},
)
    shifteddiff(field[1], grid, direction, DIM3) .* cell_extents⁻¹(grid)[3] .-
    shifteddiff(field[3], grid, direction, DIM1) .* cell_extents⁻¹(grid)[1]
end

function curlcomp(
    field::VectorField{<:Any,N},
    grid::Grid{N},
    direction::Type{<:SliceShift},
    ::Dim{3},
) where {N}
    shifteddiff(field[2], grid, direction, DIM1) .* cell_extents⁻¹(grid)[1] .-
    shifteddiff(field[1], grid, direction, DIM2) .* cell_extents⁻¹(grid)[2]
end

function grad!(
    dest::VectorField{<:Any,N},
    field::Array{<:Any,N},
    grid::Grid{N},
    direction::Type{<:SliceShift},
) where {N}
    @dbgasserteq(size(dest[1]), size(grid))
    for n = 1:N
        @dbgassert dest[n] !== field
        dest[n][interior(grid)...] .=
            shifteddiff(field, grid, direction, Dim{n}()) .*
            cell_extents⁻¹(grid)[n]
    end
end

function grad_padded!(
    dest::VectorField{<:Any,N},
    padded_field::Array{<:Any,N},
    padded_grid::Grid{N},
    direction::Type{<:SliceShift},
) where {N}
    @dbgasserteq(size(dest[1]) .+ 2, size(padded_grid))
    for n = 1:N
        dest[n] .=
            shifteddiff(padded_field, padded_grid, direction, Dim{n}()) .*
            cell_extents⁻¹(padded_grid)[n]
    end
end

function div!(
    dest::Array{<:Any,3},
    field::VectorField{<:Any,3},
    grid::Grid{3},
    direction::Type{<:SliceShift},
)
    @dbgasserteq(size(dest), size(grid))
    @dbgassert (dest !== field[1]) && (dest !== field[2]) && (dest !== field[3])
    dest[interior(grid)...] .=
        shifteddiff(field[1], grid, direction, DIM1) .*
        cell_extents⁻¹(grid)[1] .+
        shifteddiff(field[2], grid, direction, DIM2) .*
        cell_extents⁻¹(grid)[2] .+
        shifteddiff(field[3], grid, direction, DIM3) .* cell_extents⁻¹(grid)[3]
end

function div!(
    dest::Array{<:Any,2},
    field::VectorField{<:Any,2},
    grid::Grid{2},
    direction::Type{<:SliceShift},
)
    @dbgasserteq(size(dest), size(grid))
    @dbgassert (dest !== field[1]) && (dest !== field[2])
    dest[interior(grid)...] .=
        shifteddiff(field[1], grid, direction, DIM1) .*
        cell_extents⁻¹(grid)[1] .+
        shifteddiff(field[2], grid, direction, DIM2) .* cell_extents⁻¹(grid)[2]
end

function div(
    field::VectorField{<:Any,N},
    grid::Grid{N},
    direction::Type{<:SliceShift},
) where {N}
    div_field =
        shifteddiff(field[1], grid, direction, DIM1) .* cell_extents⁻¹(grid)[1]
    for n = 2:N
        div_field .+=
            shifteddiff(field[n], grid, direction, Dim{n}()) .*
            cell_extents⁻¹(grid)[n]
    end
    div_field
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
"""

 This function computes the time correlation function of a vector
 defined by the center of masses of two selections.
 The trajectory might be aligned to the first frame according to the
 minimal RMSD of a third selection.

"""

using ProgressMeter

function tcf(simulation :: Simulation,
             abs_start,abs_end,emi_start,emi_end;
             lastframe=0,
             maxdt=-1.,
             align=[],
             theta=-1.,
             r0=1.,
             scaletime=1.)

  # Check if the input is correct for using or not theta

  if theta > 0.

    # This only makes sense if the starting point of the two vectors are the same

    if abs_start != emi_start
      error(" Error: If you want to use an angle (theta) to define the emission \n",
            "        vector, the initial point of the emission vector \n",
            "        must be the same as the initial point of the absorption vector. ")
    end

    # And also, this only makes sense if the vectors are not identical

    if ( abs_start == emi_start ) && ( abs_end == emi_end )
      error(" Error: If you want to compute the emission vector from a rotation theta \n",
            "        of the absorption vector, the defined absorption and emission \n",
            "        vectors cannot be equal, because they must define the plane \n",
            "        of rotation (actually, this is the only function of the definition \n",
            "        of the emission vector in this case)")
    end

    # Converting angle

    theta = 3.141592654 * theta / 180. # theta in radians
    cost = cos(theta)
    sint = sin(theta)
    r = Matrix{Float32}(undef,3,3)

  end

  nalign = length(align)
  if nalign == 0
    doalign = false
  else
    doalign = true
    xalign_ref = Matrix{Float32}(undef,nalign+4,3)
    xalign = Matrix{Float32}(undef,nalign+4,3)
  end

  if lastframe == 0
   lastframe = simulation.nframes
  end
  println(" Last frame to consider = ",lastframe)
  xabs = Matrix{Float32}(undef,lastframe,3)
  xemi = Matrix{Float32}(undef,lastframe,3)

  if maxdt < 0.
    maxstep = lastframe
  else
    maxstep = Int64(maxdt/scaletime)
  end

  # Go to first frame, to be sure
  M3GTools.firstframe(simulation)

  p = Progress(lastframe,5," Reading DCD file: ")
  for iframe in 1:lastframe

    next!(p)

    # Reading dcd data for this frame

    sides, xdcd, ydcd, zdcd = M3GTools.nextframe(simulation)

    # Computing the absorption and emission vectors

    cm_abs_start = M3GTools.cm(abs_start,simulation,xdcd,ydcd,zdcd)
    cm_abs_end = M3GTools.cm(abs_end,simulation,xdcd,ydcd,zdcd)
    cm_emi_start = M3GTools.cm(emi_start,simulation,xdcd,ydcd,zdcd)
    cm_emi_end = M3GTools.cm(emi_end,simulation,xdcd,ydcd,zdcd)

    # If align is set, move all vectors according to the alignment of the
    # selected atoms

    if doalign == true

      if iframe == 1
        for i in 1:nalign
          xalign_ref[i,1] = xdcd[align[i]]
          xalign_ref[i,2] = ydcd[align[i]]
          xalign_ref[i,3] = zdcd[align[i]]
        end
        for j in 1:3
          xalign_ref[nalign+1,j] = cm_abs_start[j]
          xalign_ref[nalign+2,j] = cm_abs_end[j]
          xalign_ref[nalign+3,j] = cm_emi_start[j]
          xalign_ref[nalign+4,j] = cm_emi_end[j]
        end
      end

      for i in 1:nalign
        xalign[i,1] = xdcd[align[i]]
        xalign[i,2] = ydcd[align[i]]
        xalign[i,3] = zdcd[align[i]]
      end
      for j in 1:3
        xalign[nalign+1,j] = cm_abs_start[j]
        xalign[nalign+2,j] = cm_abs_end[j]
        xalign[nalign+3,j] = cm_emi_start[j]
        xalign[nalign+4,j] = cm_emi_end[j]
      end

      xaligned = M3GTools.procrustes(xalign,xalign_ref)

      for j in 1:3
        cm_abs_start[j] = xaligned[nalign+1,j]
        cm_abs_end[j]   = xaligned[nalign+2,j]
        cm_emi_start[j] = xaligned[nalign+3,j]
        cm_emi_end[j]   = xaligned[nalign+4,j]
      end

    end

    for i in 1:3
      xabs[iframe,i] = cm_abs_end[i] - cm_abs_start[i]
      xemi[iframe,i] = cm_emi_end[i] - cm_emi_start[i]
    end
    vabs_norm = sqrt( xabs[iframe,1]^2 + xabs[iframe,2]^2 + xabs[iframe,3]^2 )
    vemi_norm = sqrt( xemi[iframe,1]^2 + xemi[iframe,2]^2 + xemi[iframe,3]^2 )
    for i in 1:3
      xabs[iframe,i] = xabs[iframe,i] / vabs_norm
      xemi[iframe,i] = xemi[iframe,i] / vemi_norm
    end

    # Computing emission vector from a rotation of the absoprtion vector, if theta > 0.

    if theta > 0.

      # Rotation axis is perpendicular to both vectors (external product)

      ax = xabs[iframe,2]*xemi[iframe,3] - xabs[iframe,3]*xemi[iframe,2]
      ay = xabs[iframe,3]*xemi[iframe,1] - xabs[iframe,1]*xemi[iframe,3]
      az = xabs[iframe,1]*xemi[iframe,2] - xabs[iframe,2]*xemi[iframe,1]
      vnorm = sqrt( ax^2 + ay^2 + az^2 )
      ax = ax / vnorm
      ay = ay / vnorm
      az = az / vnorm

      # Rotation matrix (quaternion derived)

      xt = 1. - cost
      r[1,1] = cost + ax*ax*xt ; r[1,2] = ax*ay*xt-az*sint; r[1,3] = ax*az*xt+ay*sint
      r[2,1] = ay*ax*xt+az*sint; r[2,2] = cost+ay*ay*xt   ; r[2,3] = ay*az*xt-ax*sint
      r[3,1] = az*ax*xt-ay*sint; r[3,2] = az*ay*xt+ax*sint; r[3,3] = cost + az*az*xt

      # Applying rotation matrix

      xt = xabs[iframe,1]
      yt = xabs[iframe,2]
      zt = xabs[iframe,3]
      xemi[iframe,1] = r[1,1]*xt + r[1,2]*yt + r[1,3]*zt
      xemi[iframe,2] = r[2,1]*xt + r[2,2]*yt + r[2,3]*zt
      xemi[iframe,3] = r[3,1]*xt + r[3,2]*yt + r[3,3]*zt

    end

  end

  # Computing the time-dependent correlation function

  tcf = zeros(maxstep)
  legendre = zeros(maxstep)
  t = Vector{Float32}(undef,maxstep)

  p = Progress(lastframe,5," Computing the tcf: ")
  for i in 1:lastframe

    next!(p)

    for j in i:min(i+maxstep-1,lastframe)

      # Computing the internal product of absoprtion and emission vectors

      int_prod = xabs[i,1]*xemi[j,1] + xabs[i,2]*xemi[j,2] + xabs[i,3]*xemi[j,3]

      # Computing the time correlation function and its second legendre polynomial

      tcf[j-i+1] = tcf[j-i+1] + int_prod
      legendre[j-i+1] = legendre[j-i+1] + 0.5*( 3. * int_prod^2 - 1. )

    end
  end

  p = Progress(maxstep,5," Final scaling: ")
  for i in 1:maxstep
    next!(p)
    tcf[i] = tcf[i] / ( lastframe - i + 1 )
    legendre[i] = legendre[i] / ( lastframe - i + 1 )
    t[i] = scaletime*(i-1)
  end

  return t, r0*legendre, tcf

end


text:
# GR.inline("pdf")
function plot_results_dict_bar(
      results_dict,
      n_subnets,
      subnet_array,
      idx_sorted,
      output_folder,
      plot_best_x
   )
   upscale = 2 #upscaling in resolution
   fntsm = font("serif", pointsize=round(8.0*upscale))
   fntlg = font("serif", pointsize=round(12.0*upscale))
   default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   default(size=(800*upscale,600*upscale)) #Plot canvas size

   # gr()

   # plot all results
   # println("results_dict: $results_dict")
   for (label,results) in results_dict
      if (label == "status") || (label == "subnet")
         continue
      end
      for subnet_i in 0:n_subnets
         # println("subnet_i: $subnet_i")
         br_results = Array{Float64,1}()
         br_indices = Array{String,1}()
         br_results_fixf = Array{Float64,1}()
         br_indices_fixf = Array{String,1}()
         br_results_indirPQ = Array{Float64,1}()
         br_indices_indirPQ = Array{String,1}()
         base_res = -9999.0
         if subnet_i == 0
            # Don't plot total for frequency
            if label == "frequency (Hz)"
               continue
            end
            plotlabel = "total $label"
            plotdir = "total"
         else
            plotlabel = "subnet $(subnet_array[subnet_i]) $label"
            plotdir =  "subnet$(subnet_array[subnet_i])"
         end
         for (folder,res) in results
            if results_dict["status"][folder] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED]
               if subnet_i == 0
                  subnet_result = sum(res)
               else
                  try
                     subnet_result = res[subnet_i]
                  catch
                     subnet_result = 0
                  end
               end
            else
               subnet_result = -9999.0
               println("Subnet $subnet_i, network $folder: result is -9999.0")
            end
            if occursin("br", folder)
               if occursin("_fixf", folder)
                  push!(br_results_fixf, subnet_result)
                  push!(br_indices_fixf, folder[3:(end-5)])
               elseif occursin("_indirPQ", folder)
                  push!(br_results_indirPQ, subnet_result)
                  push!(br_indices_indirPQ, folder[3:(end-8)])
               else
                  push!(br_results, subnet_result)
                  push!(br_indices, folder[3:end])
               end
            # Remove base value to plot separately
            elseif occursin("base", folder)
               base_res = subnet_result
            end
         end
         # Sort values
         # sort all by primary result instead of sorting each individually
         println("label: $label")
         perm = sortperm(br_indices)[sortperm(sortperm(idx_sorted))]
         println("perm: $perm")
         println("br_results: $br_results")
         println("br_indices: $br_indices")
         br_results_sorted = br_results[perm]
         br_indices_sorted = br_indices[perm]
         println("br_results_sorted: $br_results_sorted")
         println("br_indices_sorted: $br_indices_sorted")
         if plot_best_x > 1
            br_results_sorted = br_results_sorted[1:plot_best_x]
            br_indices_sorted = br_indices_sorted[1:plot_best_x]
         end

         # If fixed f folders are present for comparison, also sort these
         if length(br_indices_fixf) > 0
            if length(br_indices_fixf) < length(br_indices_sorted)
               println("Adding results to br_indices_fixf:")
               # println("~br_indices_fixf: $br_indices_fixf")
               # println("~br_indices_sorted: $br_indices_sorted")
               # println("~br_results_fixf: $br_results_fixf")
               # println("~br_results_sorted: $br_results_sorted")
               append!(br_indices_fixf, setdiff(br_indices_sorted, br_indices_fixf))
               # println("appending $(zeros(length(br_indices_sorted)-length(br_indices_fixf)))")
               append!(br_results_fixf, -9999.0*ones(Float64, length(br_results_sorted)-length(br_results_fixf)))

            elseif length(br_indices_fixf) > length(br_indices_sorted)
               println("Adding results to br_indices_sorted:")
               # println("br_indices_fixf: $br_indices_fixf")
               # println("br_indices_sorted: $br_indices_sorted")
               # println("setdiff(br_indices_fixf, br_indices_sorted): $(setdiff(br_indices_fixf, br_indices_sorted))")
               append!(br_indices_sorted, setdiff(br_indices_fixf, br_indices_sorted))
               append!(br_results_sorted, -9999.0*ones(Float64, length(br_results_fixf)-length(br_results_sorted)))
            end
            # println("br_indices_fixf: $br_indices_fixf")
            # println("br_indices_sorted: $br_indices_sorted")
            # println("br_results_fixf: $br_results_fixf")
            # println("br_results_sorted: $br_results_sorted")
            perm_fixf = sortperm(br_indices_fixf)[sortperm(sortperm(br_indices_sorted))]
            # println("perm_fixf: $perm_fixf")
            br_indices_fixf_sorted = br_indices_fixf[perm_fixf]
            br_results_fixf_sorted = br_results_fixf[perm_fixf]
            if plot_best_x > 1
               br_indices_fixf_sorted = br_indices_fixf_sorted[1:plot_best_x]
               br_results_fixf_sorted = br_results_fixf_sorted[1:plot_best_x]
            end
         else
            br_indices_fixf = []
            br_indices_fixf_sorted = []
            br_results_fixf_sorted = []
         end

         # If indirPQ folders are present for comparison, also sort these
         if length(br_indices_indirPQ) > 0
            if length(br_indices_indirPQ) < length(br_indices_sorted)
               append!(br_indices_indirPQ, setdiff(br_indices_sorted, br_indices_indirPQ))
               append!(br_results_indirPQ, -9999.0*ones(Float64, length(br_results_sorted)-length(br_results_indirPQ)))

            elseif length(br_indices_indirPQ) > length(br_indices_sorted)
               append!(br_indices_sorted, setdiff(br_indices_indirPQ, br_indices_sorted))
               append!(br_results_sorted, -9999.0*ones(Float64, length(br_results_indirPQ)-length(br_results_sorted)))
            end
            perm_indirPQ = sortperm(br_indices_indirPQ)[sortperm(sortperm(br_indices_sorted))]
            # println("perm_fixf: $perm_fixf")
            br_indices_indirPQ_sorted = br_indices_indirPQ[perm_indirPQ]
            br_results_indirPQ_sorted = br_results_indirPQ[perm_indirPQ]
            if plot_best_x > 1
               br_indices_indirPQ_sorted = br_indices_indirPQ_sorted[1:plot_best_x]
               br_results_indirPQ_sorted = br_results_indirPQ_sorted[1:plot_best_x]
            end
         else
            br_indices_indirPQ = []
            br_indices_indirPQ_sorted = []
            br_results_indirPQ_sorted = []
         end

         # Plot subnet values
         if br_indices_fixf_sorted == br_indices_sorted == br_indices_indirPQ_sorted
            if length(br_indices_sorted) == 0
               println("No feasible solutions were found for any of the cases. No plots will be generated.")
               return
            end
            println("Plotting variable and fixed frequency and indirect PQ values.")
            println("variable f: $br_results_sorted")
            println("fixed f: $br_results_fixf_sorted")
            println("passive PQ: $br_results_indirPQ_sorted")
            grouplabels = CategoricalArray(repeat(["passive PQ", "fixed \\omega", "variable \\omega"], inner = length(br_results)), ordered=true)
            levels!(grouplabels, ["passive PQ", "fixed \\omega", "variable \\omega"])
            nam = CategoricalArray(repeat(string.(br_indices_fixf_sorted), outer = 3), ordered=true)
            levels!(nam, string.(br_indices_fixf_sorted))
            vals = vcat(br_results_indirPQ_sorted,br_results_fixf_sorted,br_results_sorted)
            p = groupedbar(
               nam,
               vals,
               xticks=(0:(length(br_indices_sorted)-1), br_indices_sorted),
               group = grouplabels,
               bar_edges=true,
               linewidth=0, left_margin=(10*upscale)*mm,
               right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
               legend=:bottomright,
               xrotation=90
            )
            if base_res != -9999.0
               hline!([base_res], line=(2*upscale, :dash), label="base network")
               miny = min(
                  min(br_results_sorted[br_results_sorted.>-9999.0]...),
                  min(br_results_fixf_sorted[br_results_fixf_sorted.>-9999.0]...),
                  min(br_results_indirPQ_sorted[br_results_indirPQ_sorted.>-9999.0]...),
                  base_res
               )
               maxy = max(
                  max(br_results_sorted...),
                  max(br_results_fixf_sorted...),
                  max(br_results_indirPQ_sorted...),
                  base_res
               )
            else
               miny = min(
                  min(br_results_sorted[br_results_sorted.>-9999.0]...),
                  min(br_results_fixf_sorted[br_results_fixf_sorted.>-9999.0]...),
                  min(br_results_indirPQ_sorted[br_results_indirPQ_sorted.>-9999.0]...)
               )
               maxy = max(
                  max(br_results_sorted...),
                  max(br_results_fixf_sorted...),
                  max(br_results_indirPQ_sorted...)
                  )
            end
         elseif br_indices_fixf_sorted == br_indices_sorted
            if length(br_indices_sorted) == 0
               println("No feasible solutions were found for any of the cases. No plots will be generated.")
               return
            end
            println("Plotting variable and fixed frequency values.")
            grouplabels = CategoricalArray(repeat(["variable \\omega", "fixed \\omega"], inner = length(br_results)), ordered=true)
            levels!(grouplabels, ["variable \\omega", "fixed \\omega"])
            nam = CategoricalArray(repeat(string.(br_indices_fixf_sorted), outer = 2), ordered=true)
            levels!(nam, string.(br_indices_fixf_sorted))
            vals = vcat(br_results_sorted,br_results_fixf_sorted)
            p = groupedbar(
               nam,
               vals,
               xticks=(0:(length(br_indices_sorted)-1), br_indices_sorted),
               group = grouplabels,
               bar_edges=true,
               linewidth=0, left_margin=(10*upscale)*mm,
               right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
               legend=:bottomright,
               xrotation=90
            )
            if base_res != -9999.0
               hline!([base_res], line=(2*upscale, :dash), label="base network")
               miny = min(
                  min(br_results_sorted[br_results_sorted.>-9999.0]...),
                  min(br_results_fixf_sorted[br_results_fixf_sorted.>-9999.0]...),
                  base_res
               )
               maxy = max(
                  max(br_results_sorted...),
                  max(br_results_fixf_sorted...),
                  base_res
               )
            else
               miny = min(
                  min(br_results_sorted[br_results_sorted.>-9999.0]...),
                  min(br_results_fixf_sorted[br_results_fixf_sorted.>-9999.0]...)
               )
               maxy = max(
                  max(br_results_sorted...),
                  max(br_results_fixf_sorted...)
                  )
            end
         else
            # println("br_indices_fixf_sorted: $br_indices_fixf_sorted")
            # println("br_indices_sorted: $br_indices_sorted")
            # println("br_results_fixf: $br_results_fixf")
            # println("br_results_sorted: $br_results_sorted")
            println("Plotting only variable frequency values.")
            nam = string.(br_indices_sorted)
            vals = br_results_sorted
            p = bar(
               nam,
               vals,
               xticks=(0:(length(br_indices_sorted)-1), br_indices_sorted),
               bar_edges=true,
               linewidth=0, left_margin=(10*upscale)*mm,
               right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
               legend=:bottomright,
               xrotation=90
            )

            if base_res != -9999.0
               hline!([base_res], line=(2*upscale, :dash), label="base network")
               # println("br_results_sorted: $br_results_sorted")
               # println("br_results_sorted[br_results_sorted.>-9999.0]: $(br_results_sorted[br_results_sorted.>-9999.0])")
               miny = min(min(br_results_sorted[br_results_sorted.>-9999.0]...), base_res)
               maxy = max(max(br_results_sorted...), base_res)
            else
               # println("br_results_sorted: $br_results_sorted")
               # println("br_results_sorted[br_results_sorted.>-9999.0]: $(br_results_sorted[br_results_sorted.>-9999.0])")
               miny = min(br_results_sorted[br_results_sorted.>-9999.0]...)
               maxy = max(br_results_sorted...)
            end
         end
         # Set y limits
         diff = maxy - miny
         if diff > 0
            ylims!(miny-0.1*diff, maxy+0.1*diff)
         end
         xlabel!("branch index")
         ylabel!(plotlabel)
         if !isdir("$(output_folder)/plots/$plotdir")
            mkpath("$(output_folder)/plots/$plotdir")
         end
         # println("saving plot at $(output_folder)/plots/$plotdir/$plotlabel.pdf")
         savefig(p,"$(output_folder)/plots/$plotdir/$plotlabel.pdf")
      end
   end
end

function plot_results_dicts_bar(
      results_dict_allplots,
      n_subnets,
      subnet_array,
      idx_sorted,
      output_folder,
      plot_best_x,
      series_labels=[];
      color_palette=:tab10
   )
   upscale = 2 #upscaling in resolution
   fntsm = font("serif", pointsize=round(8.0*upscale))
   fntlg = font("serif", pointsize=round(12.0*upscale))
   default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
   default(size=(800*upscale,600*upscale)) #Plot canvas size

   # gr()

   if isa(results_dict_allplots, Dict)
      results_dict_allplots = [results_dict_allplots]
   end
   if length(series_labels) == 0
      series_labels = [[] for i in results_dict_allplots]
   end

   n_series = length(results_dict_allplots)

   function make_result_arrays!(br_results, base_results, br_results_sorted, results_dict, series_s, plotlabels, plotdirs)
      for (label,results) in results_dict
         if (label == "status") || (label == "subnet")
            continue
         end
         if !(label in keys(br_results_all))
            br_results[label] = Dict{Int64,Any}()
            br_results_sorted[label] = Dict{Int64,Any}()
            base_results[label] = Dict{Int64,Any}()
            plotlabels[label] = Dict{Int64,Any}()
            plotdirs[label] = Dict{Int64,Any}()
         end

         for subnet_i in 0:n_subnets
            if subnet_i == 0
               # Don't plot total for frequency
               if label == "frequency (Hz)"
                  continue
               end
            end
            if series_s == 1
               base_results[label][subnet_i] = NaN
            end

            if !(subnet_i in keys(br_results_all[label]))
               br_results[label][subnet_i] = Dict{Int64,Any}()
               br_results_sorted[label][subnet_i] = Dict{Int64,Any}()
            end

            if !(series_s in keys(br_results_all[label][subnet_i]))
               br_results[label][subnet_i][series_s] = Dict{String,Any}()
               br_results_sorted[label][subnet_i][series_s] = Array{Float64,1}()
            end

            for (folder,res) in results
               if !(folder in keys(br_results_all[label][subnet_i][series_s]))
                  br_results[label][subnet_i][series_s][folder] = NaN
               end
            end
         end

         for subnet_i in 0:n_subnets
            statuses = results_dict["status"]

            if subnet_i == 0
               # Don't plot total for frequency
               if label == "frequency (Hz)"
                  continue
               end
               plotlabels[label][0] = "total $label"
               plotdirs[label][0] = "total"
            # else
            #    plotlabels[label][subnet_i] = "subnet $(subnet_array[subnet_i]) $label"
            #    plotdirs[label][subnet_i] =  "subnet$(subnet_array[subnet_i])"
            end

            for (folder,res) in results
               subnet_number = 0
               if statuses[folder] in [LOCALLY_SOLVED, ALMOST_LOCALLY_SOLVED, "LOCALLY_SOLVED", "ALMOST_LOCALLY_SOLVED"]
                  if subnet_i == 0
                     # println("$folder series $series_s")
                     # println("$label res: $res")
                     subnet_result = sum(res)
                     br_results[label][0][series_s][folder] = subnet_result
                     if occursin("base", folder) && (series_s == 1)
                        base_results[label][subnet_number] = subnet_result
                        # println("set base_results[$label][$subnet_number] to $subnet_result")
                     end
                  else
                     try
                        subnet_number = results_dict["subnet"][folder][subnet_i]
                        subnet_result = res[subnet_i]
                        plotlabels[label][subnet_number] = "subnet $(subnet_number) $label"
                        plotdirs[label][subnet_number] =  "subnet$(subnet_number)"
                        if occursin("br", folder)
                           # push!(br_results[label][subnet_number][series_s][folder], subnet_result)
                           # push!(br_indices[label][subnet_number][series_s][folder], folder)
                           br_results[label][subnet_number][series_s][folder] = subnet_result
                           # println("Subnet $subnet_number, network $folder: result is $(br_results[label][subnet_number][series_s][folder])")
                           # println("Subnet_i: $subnet_i")
                           # println("results_dict[\"subnet\"][folder]=$(results_dict["subnet"][folder])")
                        # Remove base value to plot separately
                        elseif occursin("base", folder)
                           # if series_s == 1
                              base_results[label][subnet_number] = subnet_result
                              # println("set base_results[$label][$subnet_number] to $subnet_result")
                           # end
                        end
                     catch BoundsError
                        # println("BoundsError attempting to access subnet $subnet_i")
                        # println("results_dict[\"subnet\"][$folder] = $(results_dict["subnet"][folder])")
                        # println("res = $(res)")
                        # subnet_number = length(results_dict["subnet"][folder]) + 1
                        # subnet_result = -9999.0
                     end
                  end
               end
            end
         end
      end
   end

   br_results_all = Dict{String,Any}()
   base_results = Dict{String,Any}()
   br_results_all_sorted = Dict{String,Any}()
   plotlabels = Dict{String,Any}()
   plotdirs = Dict{String,Any}()
   # series_s = 1
   # Make arrays and add values from the first series
   # make_result_arrays!(br_results_all, br_indices_all, br_results_all_sorted, br_indices_all_sorted, results_dict_allplots[series_s], series_s)
   # results_string = JSON.json(br_results_all)
   # open("./br_results_temp.json", "w") do f
   #    write(f, results_string)
   # end

   # Sort values
   # Save permutation vector to sort all by primary result instead of sorting each individually
   for (series_s,results_dict) in enumerate(results_dict_allplots)
      make_result_arrays!(br_results_all, base_results, br_results_all_sorted, results_dict_allplots[series_s], series_s, plotlabels, plotdirs)
      # println("base_results: $base_results")
      for (label,results) in results_dict
         if (label == "status") || (label == "subnet")
            continue
         end
         # println("\n$label")
         # println("=================")
         for subnet_i in 0:n_subnets
            if subnet_i == 0
               # Don't plot total for frequency
               if label == "frequency (Hz)"
                  continue
               end
            end
            # Sort values
            # println("subnet $subnet_i")
            # println("series $series_s")
            # println([v for (i,v) in br_results_all[label][subnet_i][series_s]])
            # println("idx_sorted: $idx_sorted")
            # println([br_results_all[label][subnet_i][series_s][idx] for idx in idx_sorted])
            # println("br_results_all[$label][$subnet_i]: $(br_results_all[label][subnet_i])")
            # println("idx_sorted: $(idx_sorted)")
            # println("series_s: $series_s")
            br_results_all_sorted[label][subnet_i][series_s] = [br_results_all[label][subnet_i][series_s][idx] for idx in idx_sorted]

            if plot_best_x > 1
               br_results_all_sorted[label][subnet_i][series_s] = br_results_sorted[label][subnet_i][series_s][1:plot_best_x]
            end
      #       # println("br_results_all[label][subnet_i][series_s]: $(br_results_all[label][subnet_i][series_s])")
      #       # println("idx_sorted: $(idx_sorted)")
      #       # println([br_results_all[label][subnet_i][series_s][idx] for idx in idx_sorted])
      #
      #       br_results_all_sorted[label][subnet_i][series_s] = [br_results_all[label][subnet_i][series_s][idx] for idx in idx_sorted]
      #
      #       if plot_best_x > 1
      #          br_results_all_sorted[label][subnet_i][series_s] = br_results_sorted[label][subnet_i][series_s][1:plot_best_x]
      #       end
         end
      end
   end

   for (label,results) in results_dict_allplots[1]
      if (label == "status") || (label == "subnet")
         continue
      end
      # println("Plotting $label")
      for subnet_i in 0:n_subnets
         # println("Plotting subnet $subnet_i")
         if subnet_i == 0
            # Don't plot total for frequency
            if label == "frequency (Hz)"
               continue
            end
         end
         nam = CategoricalArray(repeat(string.(idx_sorted), outer = length(results_dict_allplots)), ordered=true)
         levels!(nam, string.(idx_sorted))
         # nam = repeat(string.(idx_sorted), outer = length(results_dict_allplots))
         # println("size(nam) $(size(nam))")
         vals = vcat(Any[br_results_all_sorted[label][subnet_i][series_s] for series_s in 1:length(results_dict_allplots)]...)

         if (all(vals.<=-9999.0) || all(isnan.(vals)))
            println("All values are missing. Skipping")
            continue
         end

         # println("nam: $nam")
         # println("size(vals) $(size(vals))")
         # println("vals: $vals")
         grouplabels = CategoricalArray(repeat(series_labels, inner = length(idx_sorted)), ordered=true)
         levels!(grouplabels, series_labels)
         # grouplabels = repeat(series_labels, inner = size(vals)[2])
         # println("size(grouplabels) $(size(grouplabels))")
         # println("grouplabels: $grouplabels")

         p = groupedbar(
            nam,
            vals,
            group = grouplabels,
            xticks=(0:(length(idx_sorted)-1), idx_sorted),
            bar_edges=true,
            linewidth=0, left_margin=(10*upscale)*mm,
            right_margin=(10*upscale)*mm, top_margin=(10*upscale)*mm, bottom_margin=(10*upscale)*mm,
            legend=:bottomright,
            xrotation=90,
            palette=color_palette
         )
         base_res = base_results[label][subnet_i]
         # println("base_results[$label][$subnet_i] = $base_res")
         if !isnan(base_res) && (base_res != -9999.0)
            println("plotting base = $base_res")
            hline!([base_res], line=(2*upscale, :dash), label="base network")
            miny = min(
               min(vals[.!isnan.(vals) .& (vals.>-9999.0)]...),
               base_res
            )
            maxy = max(
               max(vals[.!isnan.(vals)]...),
               base_res
            )
         else
            println("not plotting base = $base_res")
            miny = min(vals[.!isnan.(vals) .& (vals.>-9999.0)]...)
            maxy = max(vals[.!isnan.(vals)]...)
         end

         # Set y limits
         diff = maxy - miny
         if diff > 0
            ylims!(miny-0.1*diff, maxy+0.1*diff)
         end

         plotlabel = plotlabels[label][subnet_i]
         plotdir = plotdirs[label][subnet_i]
         xlabel!("branch index")
         ylabel!(plotlabel)
         if !isdir("$(output_folder)/plots/$plotdir")
            mkpath("$(output_folder)/plots/$plotdir")
         end
         println("saving plot at $(output_folder)/plots/$plotdir/$plotlabel.pdf")
         savefig(p,"$(output_folder)/plots/$plotdir/$plotlabel.pdf")
      end
   end

   # # Now repeat for all other results dictionaries, but sort according to the first dictionary
   # for (series_s,results_dict) in enumerate(results_dict_allplots[2:end])
   #    make_result_arrays!(br_results_all, br_indices_all, results_dict, series_s)
   # end
   #
   #
   # for (label,results) in results_dict_allplots[1]
   #    if (label == "status") || (label == "subnet")
   #       continue
   #    end
   #    for subnet_i in results_dict_allplots[series_s]["subnet"]
   #       br_results_s = br_results_all[label][subnet_i][series_s]
   #       br_indices_s = br_indices_all[label][subnet_i][series_s]
   #
   #       perm_s = sortperm(br_indices_i)[sortperm(sortperm(br_indices_sorted))]
   #
   #       br_indices_s_sorted = br_indices_s[perm_fixf]
   #       br_results_s_sorted = br_results_s[perm_fixf]
   #
   #       if plot_best_x > 1
   #          br_indices_s_sorted = br_indices_s_sorted[1:plot_best_x]
   #          br_results_s_sorted = br_results_s_sorted[1:plot_best_x]
   #       end
   #
   #
   #       # Plot subnet values
   #       if length(br_indices_sorted) == 0
   #          println("No feasible solutions were found for any of the cases. No plots will be generated.")
   #          return
   #       end

   #

end


text:
module IntervalLinearAlgebra

# Write your package code here.
using LinearAlgebra, IntervalArithmetic, StaticArrays

export
    Krawczyk, Jacobi, GaussSeidel, GaussElimination, HansenBliekRohn,
    solve, enclose, precondition,
    comparison_matrix

include("solvers/solvers.jl")
include("utils.jl")
end


text:
l = @latexify dummyfunc(x; y=1, z=3) = x^2/y + z
@test l == raw"$\mathrm{dummyfunc}\left( x; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

@test_throws UndefVarError dummyfunc(1.)

l2 = @latexrun dummyfunc2(x; y=1, z=3) = x^2/y + z
@test l2 == raw"$\mathrm{dummyfunc2}\left( x; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

@test dummyfunc2(1.) == 4

l3 = @latexify dummyfunc2(x::Number; y=1, z=3) = x^2/y + z
@test l3 == raw"$\mathrm{dummyfunc2}\left( x::Number; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

l4 = @latexify dummyfunc2(::Number; y=1, z=3) = x^2/y + z
@test l4 == raw"$\mathrm{dummyfunc2}\left( ::Number; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"


text:
function FixedHeatFlow(; name, Q_flow=1.0, T₀=293.15, α=0.0)
    qflow, tem₀, alpha = Q_flow, T₀, α

    @parameters Q_flow T₀ α
    @named hp = HeatPort()

    eqs = [
        hp.Q_flow ~ -Q_flow * (1 + α*(hp.T - T₀))
    ]
    ODESystem(eqs, t, [], [Q_flow, T₀, α], systems=[hp], defaults=Dict(zip((Q_flow, T₀, α), (qflow, tem₀, alpha))), name=name)
end

function FixedTemperature(; name, T=0.0)
    tem = T

    @named hp = HeatPort()
    @parameters T

    eqs = [
        hp.T ~ T
    ]
    ODESystem(eqs, t, [], [T], systems=[hp], defaults=Dict(T => tem), name=name)
end


text:
"""
    interval_norm(A::Matrix{Interval})

computes the infinite norm of interval matrix A
"""
interval_norm(A::AbstractMatrix) = maximum(sum(mag.(A); dims=2)) # TODO: proper expand norm function and add 1-norm
interval_norm(v::AbstractVector) = maximum(mag.(v))
"""
Preconditions the interval system Ax = b by multipling by the (approximate) inverse of Ac,
that is the midpoint of A.
"""
function precondition(A, b)
    Ac = mid.(A)
    return Ac\A, Ac\b
end

"""
    enclose(A::Matrix{Interval}, b::Vector{Interval})

Computes an initial enclosure Σ so that x ⊆ Σ, where x is the solution of the interval
system Ax = b.
"""
function enclose(A::SMatrix{N, N, T, M}, b::SVector{N, T}) where {N, T, M}
    A1 = I - A
    e = interval_norm(b)/(1 - interval_norm(A1))
    x0 = MVector{N, T}(fill(-e..e, N))
    return x0
end
function enclose(A::MMatrix{N, N, T, M}, b::MVector{N, T}) where {N, T, M}
    A1 = I - A
    e = interval_norm(b)/(1 - interval_norm(A1))
    x0 = MVector{N, T}(fill(-e..e, N))
    return x0
end

function enclose(A, b)
    A1 = I - A
    e = interval_norm(b)/(1 - interval_norm(A1))
    x0 = fill(-e..e, length(b)  )
    return x0
end

"""
    comparison_matrix(A::Matrix{Interval})

Computes the comparison matrix ⟨A⟩ of the given matrix A according to the definition
⟨A⟩_ii = mig(A_ii)
⟨A⟩_ij = -mag(A_ij)
"""
function comparison_matrix(A)
    n = size(A, 1)
    compA = -mag.(A)
    @inbounds for (i, idx) in enumerate(diagind(A))
        compA = setindex(compA, mig(A[i, i]), idx)
    end
    return compA
end


text:
WhereTraits.BasicTraits.@overwrite_Base
using Test

# we only test for subset relationship because traits might be added during this test run
@test isbitstype(Int)
@test isiterable(Int)
@test isimmutable(Int)
@test isconcretetype(Int)

@test isbitstype(typeof(+))
@test iscallable(typeof(+))
@test isimmutable(typeof(+))
@test isconcretetype(typeof(+))

struct MyImmutableType end
@test isbitstype(MyImmutableType)
@test isimmutable(MyImmutableType)
@test isconcretetype(MyImmutableType)

mutable struct MyMutableType end
@test !isbitstype(MyMutableType)
@test ismutable(MyMutableType)
@test isconcretetype(MyMutableType)

@test !isconcretetype(Dict)
@test !isconcretetype(Any)
@test isconcretetype(Int)


text:
# Run Match tests

include(joinpath(dirname(@__FILE__), "matchtests.jl"))


text:
struct PoseUniform
    low::Vector{Float64}
    upp::Vector{Float64}
    uni::Uniform{Float64}
    function PoseUniform(xlim::Vector{Float64}, ylim::Vector{Float64})
        new([xlim[1], ylim[1], 0], [xlim[2], ylim[2], 2 * pi], Uniform())
    end
end

function uniform(mv::PoseUniform)
    c = mv.upp - mv.low
    x = [rand(mv.uni) for i = 1:3]
    return mv.low .+ (x .* c)
end

@userplot MyEllipse

@recipe function f(
    c::MyEllipse;
    n_std = 1,
    n_ellipse_vertices = 100,
    alpha = 0.0,
    color = "blue",
)
    μ, S = _covellipse_args(c.args; n_std = n_std)

    θ = range(0, 2π; length = n_ellipse_vertices)
    A = S * [cos.(θ)'; sin.(θ)']
    @series begin
        fillalpha --> 0.0
        linecolor --> color
        Shape(μ[1] .+ A[1, :], μ[2] .+ A[2, :])
    end
end

function _covellipse_args(
    (μ, Σ)::Tuple{AbstractVector{<:Real},AbstractMatrix{<:Real}};
    n_std::Real,
)
    size(μ) == (2,) && size(Σ) == (2, 2) ||
        error("covellipse requires mean of length 2 and covariance of size 2×2.")
    λ, U = eigen(Σ)
    μ, n_std * U * diagm(.√λ)
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
"""
    _SOCtoPSDaff{T}(f::MOI.VectorAffineFunction{T}, g::MOI.ScalarAffineFunction{T})

Builds a VectorAffineFunction representing the upper (or lower) triangular part of the matrix
[ f[1]     f[2:end]' ]
[ f[2:end] g * I     ]
"""
function _SOCtoPSDaff(f::MOI.VectorAffineFunction{T}, g::MOI.ScalarAffineFunction{T}) where T
    dim = MOI.output_dimension(f)
    n = div(dim * (dim+1), 2)
    # Needs to add t*I
    N0 = length(f.terms)
    Ni = length(g.terms)
    N = N0 + (dim-1) * Ni
    terms = Vector{MOI.VectorAffineTerm{T}}(undef, N)
    terms[1:N0] = MOI.VectorAffineTerm.(map(t -> trimap.(t.output_index, 1), f.terms),
                                        MOI.ScalarAffineTerm.(map(t -> t.scalar_term.coefficient, f.terms),
                                                              map(t -> t.scalar_term.variable_index, f.terms)))
    constant = [f.constants; zeros(T, n - dim)]
    cur = N0
    for i in 2:dim
        k = trimap(i, i)
        terms[cur.+(1:Ni)]  = MOI.VectorAffineTerm.(k, MOI.ScalarAffineTerm.(map(t -> t.coefficient, g.terms),
                                                                             map(t -> t.variable_index, g.terms)))
        constant[k] = g.constant
        cur += Ni
    end
    MOI.VectorAffineFunction(terms, constant)
end

"""
The `SOCtoPSDBridge` transforms the second order cone constraint ``\\lVert x \\rVert \\le t`` into the semidefinite cone constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & tI
\\end{pmatrix} \\succeq 0
```
Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  tI & \\succ 0\\\\
  t - x^\\top (tI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  t & > 0\\\\
  t^2 & > x^\\top x
\\end{align*}
```
"""
struct SOCtoPSDBridge{T} <: AbstractBridge
    dim::Int
    cr::CI{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}
end
function bridge_constraint(::Type{SOCtoPSDBridge{T}}, model, f,
                           s::MOI.SecondOrderCone) where T
    d = MOI.dimension(s)
    cr = MOI.add_constraint(model, _SOCtoPSDaff(f, T), MOI.PositiveSemidefiniteConeTriangle(d))
    SOCtoPSDBridge(d, cr)
end

_SOCtoPSDaff(f::MOI.VectorOfVariables, ::Type{T}) where T = _SOCtoPSDaff(MOI.VectorAffineFunction{T}(f), T)
_SOCtoPSDaff(f::MOI.VectorAffineFunction, ::Type) = _SOCtoPSDaff(f, MOIU.eachscalar(f)[1])

MOI.supports_constraint(::Type{SOCtoPSDBridge{T}}, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}, ::Type{MOI.SecondOrderCone}) where T = true
added_constraint_types(::Type{SOCtoPSDBridge{T}}, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}, ::Type{MOI.SecondOrderCone}) where T = [(MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle)]

function MOI.get(model::MOI.AbstractOptimizer, a::MOI.ConstraintPrimal, c::SOCtoPSDBridge)
    MOI.get(model, a, c.cr)[trimap.(1:c.dim, 1)]
end
function MOI.get(model::MOI.AbstractOptimizer, a::MOI.ConstraintDual, c::SOCtoPSDBridge)
    dual = MOI.get(model, a, c.cr)
    tdual = sum(i -> dual[trimap(i, i)], 1:c.dim)
    [tdual; dual[trimap.(2:c.dim, 1)]*2]
end

MOI.get(::SOCtoPSDBridge{T}, ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}) where T = 1
MOI.get(b::SOCtoPSDBridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}) where T = [b.cr]

function MOI.delete(model::MOI.AbstractOptimizer, c::SOCtoPSDBridge)
    MOI.delete(model, c.cr)
end

"""
The `RSOCtoPSDBridge` transforms the second order cone constraint ``\\lVert x \\rVert \\le 2tu`` with ``u \\ge 0`` into the semidefinite cone constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & 2uI
\\end{pmatrix} \\succeq 0
```
Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  uI & \\succ 0\\\\
  t - x^\\top (2uI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  u & > 0\\\\
  2tu & > x^\\top x
\\end{align*}
```
"""
struct RSOCtoPSDBridge{T} <: AbstractBridge
    dim::Int
    cr::CI{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}
end

MOI.supports_constraint(::Type{RSOCtoPSDBridge{T}}, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}, ::Type{MOI.RotatedSecondOrderCone}) where T = true
added_constraint_types(::Type{RSOCtoPSDBridge{T}}, ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}, ::Type{MOI.RotatedSecondOrderCone}) where T = [(MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle)]

function bridge_constraint(::Type{RSOCtoPSDBridge{T}}, model, f,
                           s::MOI.RotatedSecondOrderCone) where T
    d = MOI.dimension(s)-1
    cr = MOI.add_constraint(model, _RSOCtoPSDaff(f, T), MOI.PositiveSemidefiniteConeTriangle(d))
    RSOCtoPSDBridge(d, cr)
end

_RSOCtoPSDaff(f::MOI.VectorOfVariables, ::Type{T}) where T = _RSOCtoPSDaff(MOI.VectorAffineFunction{T}(f), T)
function _RSOCtoPSDaff(f::MOI.VectorAffineFunction, ::Type{T}) where T
    n = MOI.output_dimension(f)
    f_scalars = MOIU.eachscalar(f)
    g = MOIU.operate!(*, T, f_scalars[2], convert(T, 2))
    _SOCtoPSDaff(f_scalars[[1; 3:n]], g)
end

function MOI.get(model::MOI.AbstractOptimizer, a::MOI.ConstraintPrimal, c::RSOCtoPSDBridge)
    x = MOI.get(model, MOI.ConstraintPrimal(), c.cr)[[trimap(1, 1); trimap(2, 2); trimap.(2:c.dim, 1)]]
    x[2] /= 2 # It is (2u*I)[1,1] so it needs to be divided by 2 to get u
    x
end
function MOI.get(model::MOI.AbstractOptimizer, a::MOI.ConstraintDual, c::RSOCtoPSDBridge)
    dual = MOI.get(model, MOI.ConstraintDual(), c.cr)
    udual = sum(i -> dual[trimap(i, i)], 2:c.dim)
    [dual[1]; 2udual; dual[trimap.(2:c.dim, 1)]*2]
end

MOI.get(::RSOCtoPSDBridge{T}, ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}) where T = 1
MOI.get(b::RSOCtoPSDBridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, MOI.PositiveSemidefiniteConeTriangle}) where T = [b.cr]

function MOI.delete(model::MOI.AbstractOptimizer, c::RSOCtoPSDBridge)
    MOI.delete(model, c.cr)
end


text:
import FFTW

# returns the lengths of the bounding rectangle in reciprocal space
# that encloses the sphere of radius Gmax
function bounding_rectangle(lattice::AbstractMatrix{T}, Gmax; tol=sqrt(eps(T))) where {T}
    # If |B G| ≤ Gmax, then
    # |Gi| = |e_i^T B^-1 B G| ≤ |B^-T e_i| Gmax = |A_i| Gmax
    # with B the reciprocal lattice matrix, e_i the i-th canonical
    # basis vector and A_i the i-th column of the lattice matrix
    Glims = [norm(lattice[:, i]) / 2T(π) * Gmax for i in 1:3]

    # Round up, unless exactly zero (in which case keep it zero in
    # order to just have one G vector for 1D or 2D systems)
    Glims = [Glim == 0 ? 0 : ceil(Int, Glim .- tol) for Glim in Glims]
    Glims
end
function diameter(lattice)
    diam = zero(eltype(lattice))
    # brute force search
    for vec in Vec3.(Iterators.product(-1:1, -1:1, -1:1))
        diam = max(diam, norm(lattice*vec))
    end
    diam
end

const _smallprimes = [2, 3, 5]

@doc raw"""
Determine the minimal grid size for the cubic basis set to be able to
represent product of orbitals (with the default `supersampling=2`).

Optionally optimize the grid afterwards for the FFT procedure by
ensuring factorization into small primes.

The function will determine the smallest parallelepiped containing the wave vectors
 ``|G|^2/2 \leq E_\text{cut} ⋅ \text{supersampling}^2``.
For an exact representation of the density resulting from wave functions
represented in the spherical basis sets, `supersampling` should be at least `2`.
"""
function compute_fft_size(lattice::AbstractMatrix{T}, Ecut;
                          supersampling=2,
                          tol=sqrt(eps(T)),
                          ensure_smallprimes=true) where T
    Gmax = supersampling * sqrt(2 * Ecut)
    Glims = bounding_rectangle(lattice, Gmax; tol=tol)

    fft_size = Vec3(2 .* Glims .+ 1)
    # Optimize FFT grid size: Make sure the number factorises in small primes only
    if ensure_smallprimes
        fft_size = nextprod.(Ref(_smallprimes), fft_size)
    end
    fft_size
end
function compute_fft_size(model::Model, Ecut; kwargs...)
    compute_fft_size(model.lattice, Ecut; kwargs...)
end

# This uses a more precise and slower algorithm than the one above,
# simply enumerating all G vectors and seeing where their difference
# is. It needs the kpoints to do so.
# TODO This function is strange ... it should only depend on the kcoords
#      It should be merged with build_kpoints somehow
function compute_fft_size_precise(lattice::AbstractMatrix{T}, Ecut, kpoints;
                                  supersampling=2, ensure_smallprimes=true) where T
    recip_lattice  = compute_recip_lattice(lattice)
    recip_diameter = diameter(recip_lattice)
    Glims = [0, 0, 0]
    # get the bounding rectangle that contains all G-G' vectors
    # (and therefore densities and potentials)
    # This handles the case `supersampling=2`
    for kpt in kpoints
        for G in G_vectors(kpt)
            if norm(recip_lattice * (G + kpt.coordinate)) ≤ sqrt(2Ecut) - recip_diameter
                # each of the 8 neighbors (in ∞-norm) also belongs to the grid
                # so we can safely skip the search knowing at least one of them
                # will have larger |G-Gp|.
                # Savings with this trick are surprisingly small :
                # for silicon, 50% at Ecut 30, 70% at Ecut 100
                continue
            end
            for Gp in G_vectors(kpt)
                for i = 1:3
                    @inbounds Glims[i] = max(Glims[i], abs(G[i] - Gp[i]))
                end
            end
        end
    end
    if supersampling != 2
        # no guarantees there, we just do our best to satisfy the
        # target supersampling ratio
        Glims = round.(Int, supersampling ./ 2 .* Glims)
    end

    fft_size = Vec3(2 .* Glims .+ 1)
    # Optimize FFT grid size: Make sure the number factorises in small primes only
    if ensure_smallprimes
        fft_size = nextprod.(Ref(_smallprimes), fft_size)
    end
    fft_size
end

# Main entry point from pwbasis. Uses the above functions to find out
# the correct fft_size according to user specification
function compute_fft_size(model::Model{T}, Ecut, supersampling,
                          optimize_fft_size, kcoords) where {T}
    fft_size = compute_fft_size(model, Ecut; supersampling=supersampling)
    if optimize_fft_size
        # We build a temporary set of kpoints here
        # We don't reuse this kpoint construction for the pwbasis
        # because build_kpoints builds index mapping from the
        # k-point-specific basis to the global basis and thus the
        # fft_size needs to be final at kpoint construction time
        fft_size = Tuple{Int, Int, Int}(fft_size)
        kpoints_temp = build_kpoints(model, fft_size, kcoords, Ecut)
        fft_size = compute_fft_size_precise(model.lattice, Ecut, kpoints_temp;
                                            supersampling=supersampling)
    end

    # TODO generic FFT is kind of broken for some fft sizes
    #      ... temporary workaround, see more details in workarounds/fft_generic.jl
    fft_size = next_working_fft_size(T, fft_size)
    Tuple{Int, Int, Int}(fft_size)
end

# For Float32 there are issues with aligned FFTW plans, so we
# fall back to unaligned FFTW plans (which are generally discouraged).
_fftw_flags(::Type{Float32}) = FFTW.MEASURE | FFTW.UNALIGNED
_fftw_flags(::Type{Float64}) = FFTW.MEASURE

"""
Plan a FFT of type `T` and size `fft_size`, spending some time on finding an
optimal algorithm. (Inplace, out-of-place) x (forward, backward) FFT plans are returned.
"""
function build_fft_plans(T::Union{Type{Float32}, Type{Float64}}, fft_size)
    tmp = Array{Complex{T}}(undef, fft_size...)
    ipFFT = FFTW.plan_fft!(tmp, flags=_fftw_flags(T))
    opFFT = FFTW.plan_fft(tmp, flags=_fftw_flags(T))
    # backward by inverting and stripping off normalizations
    ipFFT, opFFT, inv(ipFFT).p, inv(opFFT).p
end

# TODO Some grid sizes are broken in the generic FFT implementation
# in FourierTransforms, for more details see workarounds/fft_generic.jl
# This function is needed to provide a noop fallback for grid adjustment for
# for floating-point types natively supported by FFTW
next_working_fft_size(::Type{Float32}, size::Int) = size
next_working_fft_size(::Type{Float64}, size::Int) = size
next_working_fft_size(T, sizes::Union{Tuple, AbstractArray}) = next_working_fft_size.(T, sizes)


text:
import POMDPToolbox: obs_weight
import POMDPs: observation

struct P <: POMDP{Void, Void, Void} end

@test !@implemented obs_weight(::P, ::Void, ::Void, ::Void, ::Void)
@test !@implemented obs_weight(::P, ::Void, ::Void, ::Void)
@test !@implemented obs_weight(::P, ::Void, ::Void)

obs_weight(::P, ::Void, ::Void, ::Void) = 1.0
@test @implemented obs_weight(::P, ::Void, ::Void, ::Void)
@test @implemented obs_weight(::P, ::Void, ::Void, ::Void, ::Void)
@test !@implemented obs_weight(::P, ::Void, ::Void)

@test obs_weight(P(), nothing, nothing, nothing, nothing) == 1.0

observation(::P, ::Void) = nothing
@test @implemented obs_weight(::P, ::Void, ::Void)


text:
using Images
using ImageView

function imresize!(resized, original)
    scale1 = (size(original,1)-1)/(size(resized,1)-0.999f0)
    scale2 = (size(original,2)-1)/(size(resized,2)-0.999f0)
    for jr = 0:size(resized,2)-1
        jo = scale2*jr
        ijo = itrunc(jo)
        fjo = jo - oftype(jo, ijo)
        @inbounds for ir = 0:size(resized,1)-1
            io = scale1*ir
            iio = itrunc(io)
            fio = io - oftype(io, iio)
            tmp = (1-fio)*((1-fjo)*original[iio+1,ijo+1] +
                           fjo*original[iio+1,ijo+2]) +
            fio*((1-fjo)*original[iio+2,ijo+1] +
                   fjo*original[iio+2,ijo+2])
            resized[ir+1,jr+1] = convertsafely(eltype(resized), tmp)
        end
    end
    resized
end
imresize(original, new_size) =
    imresize!(similar(original,
                      new_size), original)

convertsafely{T<:FloatingPoint}(::Type{T}, val) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::Integer) = convert(T, val)
convertsafely{T<:Integer}(::Type{T}, val::FloatingPoint) =
    itrunc(T, val+oftype(val, 0.5))

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
function fresnelS0(x::T) where {T}
    x2 = x*x
    x4 = x2*x2
    n = @evalpoly(x4, T(3.18016297876567817986E11),
                      -T(4.42979518059697779103E10),
                      T(2.54890880573376359104E9),
                      -T(6.29741486205862506537E7),
                      T(7.08840045257738576863E5),
                      -T(2.99181919401019853726E3))
    d = @evalpoly(x4, T(6.07366389490084639049E11),
                      T(2.24411795645340920940E10),
                      T(4.19320245898111231129E8),
                      T(5.17343888770096400730E6),
                      T(4.55847810806532581675E4),
                      T(2.81376268889994315696E2),
                      T(1))
    x*x2*n/d
end
function fresnelC0(x::T) where {T}
    x2 = x*x
    x4 = x2*x2
    n = @evalpoly(x4, T(9.99999999999999998822E-1),
                      -T(2.05525900955013891793E-1),
                      T(1.88843319396703850064E-2),
                      -T(6.45191435683965050962E-4),
                      T(9.50428062829859605134E-6),
                      -T(4.98843114573573548651E-8))
    d = @evalpoly(x4, T(1.00000000000000000118E0),
                      T(4.12142090722199792936E-2),
                      T(8.68029542941784300606E-4),
                      T(1.22262789024179030997E-5),
                      T(1.25001862479598821474E-7),
                      T(9.15439215774657478799E-10),
                      T(3.99982968972495980367E-12))
    x*n/d
end
fresnelS(x::T) where {T} = x < T(1.6) ? fresnelS0(x) : imag((1 + im)*erf(x*sqrt(T(pi))*(1 - im)/2)/2)
fresnelC(x::T) where {T} = x < T(1.6) ? fresnelC0(x) : real((1 + im)*erf(x*sqrt(T(pi))*(1 - im)/2)/2)