text:
struct Summary <: Output
  filename::String
  last_infections::Vector{Float32}
  num_infections::Vector{UInt32}
  peak_daily_infections::Vector{UInt32}
  peak_daily_detections::Vector{UInt32}
  Summary(filename::AbstractString, num_trajectories::Integer) = new(
    filename,
    fill(NaN32, num_trajectories),
    Vector{UInt32}(undef, num_trajectories),
    Vector{UInt32}(undef, num_trajectories),
    Vector{UInt32}(undef, num_trajectories))
end

time_or_nan(e) = MocosSim.istransmission(e) ? Float32(MocosSim.time(e)) : NaN32

function pushtrajectory!(s::Summary, trajectory_id::Integer, ::Base.AbstractLock, state::MocosSim.SimState, ::MocosSim.SimParams, cb::DetectionCallback)
  max_time = maximum(MocosSim.time, state.forest.inedges)
  s.last_infections[trajectory_id] = max_time
  s.num_infections[trajectory_id] = count(MocosSim.istransmission, state.forest.inedges)
  counts = Vector{UInt32}(undef, ceil(Int, max_time))
  s.peak_daily_infections[trajectory_id] = maximum(daily!(time_or_nan, counts, state.forest.inedges), init=0)
  s.peak_daily_detections[trajectory_id] = maximum(daily!(counts, skipmissing(cb.detection_times)), init=0)
  nothing
end

function aftertrajectories(s::Summary, ::MocosSim.SimParams)
  file = jldopen(s.filename, "w")
  try
    file["num_infections"] = s.num_infections
    file["last_infections"] = s.last_infections
    file["peak_daily_detections"] = s.peak_daily_detections
    file["peak_daily_infections"] = s.peak_daily_infections
  finally
    close(file)
  end
  nothing
end


text:
"""
ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника. Робот - в произвольной клетке поля между внешней и внутренней перегородками.
РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.
"""
function framer(r::Robot)
    steps=0
    arr_of_steps=[]
    while ((isborder(r,Sud)) && (isborder(r,West)))==false
        push!(arr_of_steps,moves!(r,West))
        push!(arr_of_steps,moves!(r,Sud))
        steps=steps+2
    end
    side=Nord
    while isborder(r,Ost)==false
        while isborder(r,side)==false
            move!(r,side)
            if isborder(r,Ost)==true #Встретилась перегородка
                break #Выход из цикла
            end
        end
        if isborder(r,Ost)==false
            move!(r,Ost)
        end
        side=inverse(side) #Смена направления на противоположное
    end
    while isborder(r,Ost)==true
        move!(r,Sud)
    end
    side=Nord
    for a ∈ 1:4
        putmarker!(r)
        move!(r,side)
        while isborder(r,next(side))==true #Установка маркеров по периметру рамки
            putmarker!(r)
            move!(r,side)
        end
        side=next(side) #Смена направления на следующее
    end
    moves!(r,Sud)
    moves!(r,West)
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for a ∈ 1:arr_of_steps[steps]
            move!(r,side)
        end
        steps=steps-1
    end
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function inverse(side::HorizonSide)
    HorizonSide(mod(Int(side)+2,4))
end
function next(side::HorizonSide)
    HorizonSide(mod(Int(side)+3,4))
end
#= Действия, описание которых начинается с 39-й строки,
выполняются аналогично тем, которые описаны в Задаче 5 =#

ext: jl

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Convenience functions for metaprogramming.
"""
module Meta

using ..CoreLogging

export quot,
       isexpr,
       isidentifier,
       isoperator,
       isunaryoperator,
       isbinaryoperator,
       ispostfixoperator,
       replace_sourceloc!,
       show_sexpr,
       @dump

using Base: isidentifier, isoperator, isunaryoperator, isbinaryoperator, ispostfixoperator
import Base: isexpr

"""
    Meta.quot(ex)::Expr

Quote expression `ex` to produce an expression with head `quote`. This can for
instance be used to represent objects of type `Expr` in the AST. See also the
manual section about [QuoteNode](@ref man-quote-node).

# Examples
```jldoctest
julia> eval(Meta.quot(:x))
:x

julia> dump(Meta.quot(:x))
Expr
  head: Symbol quote
  args: Array{Any}((1,))
    1: Symbol x

julia> eval(Meta.quot(:(1+2)))
:(1 + 2)
```
"""
quot(ex) = Expr(:quote, ex)

"""
    Meta.isexpr(ex, head[, n])::Bool

Return true if `ex` is an `Expr` with the given type `head` and optionally that
the argument list is of length `n`. `head` may be a `Symbol` or collection of
`Symbol`s. For example, to check that a macro was passed a function call
expression, you might use `isexpr(ex, :call)`.

# Examples
```jldoctest
julia> ex = :(f(x))
:(f(x))

julia> Meta.isexpr(ex, :block)
false

julia> Meta.isexpr(ex, :call)
true

julia> Meta.isexpr(ex, [:block, :call]) # multiple possible heads
true

julia> Meta.isexpr(ex, :call, 1)
false

julia> Meta.isexpr(ex, :call, 2)
true
```
"""
isexpr(@nospecialize(ex), heads) = isa(ex, Expr) && in(ex.head, heads)
isexpr(@nospecialize(ex), heads, n::Int) = isa(ex, Expr) && in(ex.head, heads) && length(ex.args) == n

"""
    replace_sourceloc!(location, expr)

Overwrite the caller source location for each macro call in `expr`, returning
the resulting AST.  This is useful when you need to wrap a macro inside a
macro, and want the inner macro to see the `__source__` location of the outer
macro.  For example:

```
macro test_is_one(ex)
    replace_sourceloc!(__source__, :(@test \$(esc(ex)) == 1))
end
@test_is_one 2
```

`@test` now reports the location of the call `@test_is_one 2` to the user,
rather than line 2 where `@test` is used as an implementation detail.
"""
function replace_sourceloc!(sourceloc, @nospecialize(ex))
    if ex isa Expr
        if ex.head == :macrocall
            ex.args[2] = sourceloc
        end
        map!(e -> replace_sourceloc!(sourceloc, e), ex.args, ex.args)
    end
    return ex
end

"""
    Meta.show_sexpr([io::IO,], ex)

Show expression `ex` as a lisp style S-expression.

# Examples
```jldoctest
julia> Meta.show_sexpr(:(f(x, g(y,z))))
(:call, :f, :x, (:call, :g, :y, :z))
```
"""
show_sexpr(ex) = show_sexpr(stdout, ex)
show_sexpr(io::IO, ex) = show_sexpr(io, ex, 0)
show_sexpr(io::IO, ex, indent::Int) = show(io, ex)

const sexpr_indent_width = 2

function show_sexpr(io::IO, ex::QuoteNode, indent::Int)
    inner = indent + sexpr_indent_width
    print(io, "(:quote, #QuoteNode\n", " "^inner)
    show_sexpr(io, ex.value, inner)
    print(io, '\n', " "^indent, ')')
end
function show_sexpr(io::IO, ex::Expr, indent::Int)
    inner = indent + sexpr_indent_width
    print(io, '(')
    show_sexpr(io, ex.head, inner)
    for arg in ex.args
        print(io, ex.head === :block ? ",\n"*" "^inner : ", ")
        show_sexpr(io, arg, inner)
    end
    if isempty(ex.args)
        print(io, ",)")
    else
        print(io, (ex.head === :block ? "\n"*" "^indent : ""), ')')
    end
end

"""
    @dump expr

Show every part of the representation of the given expression. Equivalent to
[`dump(:(expr))`](@ref dump).
"""
macro dump(expr)
    return :(dump($(QuoteNode(expr))))
end

"""
    lower(m, x)

Takes the expression `x` and returns an equivalent expression in lowered form
for executing in module `m`.
See also [`code_lowered`](@ref).
"""
lower(m::Module, @nospecialize(x)) = ccall(:jl_expand, Any, (Any, Any), x, m)

"""
    @lower [m] x

Return lowered form of the expression `x` in module `m`.
By default `m` is the module in which the macro is called.
See also [`lower`](@ref).
"""
macro lower(code)
    return :(lower($__module__, $(QuoteNode(code))))
end
macro lower(mod, code)
    return :(lower($(esc(mod)), $(QuoteNode(code))))
end

## interface to parser ##

"""
    ParseError(msg)

The expression passed to the [`parse`](@ref) function could not be interpreted as a valid Julia
expression.
"""
struct ParseError <: Exception
    msg::String
end

function _parse_string(text::AbstractString, filename::AbstractString,
                       index::Integer, options)
    if index < 1 || index > ncodeunits(text) + 1
        throw(BoundsError(text, index))
    end
    ex, offset::Int = Core._parse(text, filename, index-1, options)
    ex, offset+1
end

"""
    parse(str, start; greedy=true, raise=true, depwarn=true)

Parse the expression string and return an expression (which could later be
passed to eval for execution). `start` is the code unit index into `str` of the
first character to start parsing at (as with all string indexing, these are not
character indices). If `greedy` is `true` (default), `parse` will try to consume
as much input as it can; otherwise, it will stop as soon as it has parsed a
valid expression. Incomplete but otherwise syntactically valid expressions will
return `Expr(:incomplete, "(error message)")`. If `raise` is `true` (default),
syntax errors other than incomplete expressions will raise an error. If `raise`
is `false`, `parse` will return an expression that will raise an error upon
evaluation. If `depwarn` is `false`, deprecation warnings will be suppressed.

```jldoctest
julia> Meta.parse("(α, β) = 3, 5", 1) # start of string
(:((α, β) = (3, 5)), 16)

julia> Meta.parse("(α, β) = 3, 5", 1, greedy=false)
(:((α, β)), 9)

julia> Meta.parse("(α, β) = 3, 5", 16) # end of string
(nothing, 16)

julia> Meta.parse("(α, β) = 3, 5", 11) # index of 3
(:((3, 5)), 16)

julia> Meta.parse("(α, β) = 3, 5", 11, greedy=false)
(3, 13)
```
"""
function parse(str::AbstractString, pos::Integer; greedy::Bool=true, raise::Bool=true,
               depwarn::Bool=true)
    ex, pos = _parse_string(str, "none", pos, greedy ? :statement : :atom)
    if raise && isa(ex,Expr) && ex.head === :error
        throw(ParseError(ex.args[1]))
    end
    return ex, pos
end

"""
    parse(str; raise=true, depwarn=true)

Parse the expression string greedily, returning a single expression. An error is thrown if
there are additional characters after the first expression. If `raise` is `true` (default),
syntax errors will raise an error; otherwise, `parse` will return an expression that will
raise an error upon evaluation.  If `depwarn` is `false`, deprecation warnings will be
suppressed.

```jldoctest
julia> Meta.parse("x = 3")
:(x = 3)

julia> Meta.parse("x = ")
:($(Expr(:incomplete, "incomplete: premature end of input")))

julia> Meta.parse("1.0.2")
ERROR: Base.Meta.ParseError("invalid numeric constant \\\"1.0.\\\"")
Stacktrace:
[...]

julia> Meta.parse("1.0.2"; raise = false)
:($(Expr(:error, "invalid numeric constant \"1.0.\"")))
```
"""
function parse(str::AbstractString; raise::Bool=true, depwarn::Bool=true)
    ex, pos = parse(str, 1, greedy=true, raise=raise, depwarn=depwarn)
    if isa(ex,Expr) && ex.head === :error
        return ex
    end
    if pos <= ncodeunits(str)
        raise && throw(ParseError("extra token after end of expression"))
        return Expr(:error, "extra token after end of expression")
    end
    return ex
end

function parseatom(text::AbstractString, pos::Integer; filename="none")
    return _parse_string(text, String(filename), pos, :atom)
end

function parseall(text::AbstractString; filename="none")
    ex,_ = _parse_string(text, String(filename), 1, :all)
    return ex
end

"""
    partially_inline!(code::Vector{Any}, slot_replacements::Vector{Any},
                      type_signature::Type{<:Tuple}, static_param_values::Vector{Any},
                      slot_offset::Int, statement_offset::Int,
                      boundscheck::Symbol)

Return `code` after performing an in-place partial inlining pass on the Julia IR stored
within it.

The kind of inlining transformations performed by this function are those that are generally
possible given only a runtime type signature for a method invocation and the corresponding
method's lowered IR. Thus, this function is mainly useful when preparing Julia IR to be
emitted from a `@generated` function.

The performed transformations are:

- replace slot numbers in the range `1:length(slot_replacements)` with the corresponding items in `slot_replacements`
- increment other slot numbers by `slot_offset`
- substitute static parameter placeholders (e.g. `Expr(:static_parameter, 1)`) with the corresponding
values in `static_param_values`
- increment any statement indices present in the IR (`GotoNode`s, `SSAValue`s, etc.) by `statement_offset`
(useful when the caller plans to prepend new statements to the IR)
- turn off boundschecking (if `boundscheck === :off`) or propagate boundschecking (if `boundscheck === :propagate`)

This function is similar to `Core.Compiler.ssa_substitute!`, but works on pre-type-inference
IR instead of the optimizer's IR.
"""
function partially_inline!(code::Vector{Any}, slot_replacements::Vector{Any},
                           @nospecialize(type_signature)#=::Type{<:Tuple}=#,
                           static_param_values::Vector{Any},
                           slot_offset::Int, statement_offset::Int,
                           boundscheck::Symbol)
    for i = 1:length(code)
        isassigned(code, i) || continue
        code[i] = _partially_inline!(code[i], slot_replacements, type_signature,
                                     static_param_values, slot_offset,
                                     statement_offset, boundscheck)
    end
    return code
end

function _partially_inline!(@nospecialize(x), slot_replacements::Vector{Any},
                            @nospecialize(type_signature), static_param_values::Vector{Any},
                            slot_offset::Int, statement_offset::Int,
                            boundscheck::Symbol)
    if isa(x, Core.SSAValue)
        return Core.SSAValue(x.id + statement_offset)
    end
    if isa(x, Core.GotoNode)
        return Core.GotoNode(x.label + statement_offset)
    end
    if isa(x, Core.SlotNumber)
        id = x.id
        if 1 <= id <= length(slot_replacements)
            return slot_replacements[id]
        end
        return Core.SlotNumber(id + slot_offset)
    end
    if isa(x, Core.NewvarNode)
        return Core.NewvarNode(_partially_inline!(x.slot, slot_replacements, type_signature,
                                                  static_param_values, slot_offset,
                                                  statement_offset, boundscheck))
    end
    if isa(x, Core.PhiNode)
        partially_inline!(x.values, slot_replacements, type_signature, static_param_values,
                          slot_offset, statement_offset, boundscheck)
        x.edges .+= slot_offset
        return x
    end
    if isa(x, Core.ReturnNode)
        return Core.ReturnNode(
            _partially_inline!(x.val, slot_replacements, type_signature, static_param_values,
                               slot_offset, statement_offset, boundscheck),
        )
    end
    if isa(x, Core.GotoIfNot)
        return Core.GotoIfNot(
            _partially_inline!(x.cond, slot_replacements, type_signature, static_param_values,
                               slot_offset, statement_offset, boundscheck),
            x.dest + statement_offset,
        )
    end
    if isa(x, Expr)
        head = x.head
        if head === :static_parameter
            return QuoteNode(static_param_values[x.args[1]])
        elseif head === :cfunction
            @assert !isa(type_signature, UnionAll) || !isempty(spvals)
            if !isa(x.args[2], QuoteNode) # very common no-op
                x.args[2] = _partially_inline!(x.args[2], slot_replacements, type_signature,
                                               static_param_values, slot_offset,
                                               statement_offset, boundscheck)
            end
            x.args[3] = _instantiate_type_in_env(x.args[3], type_signature, static_param_values)
            x.args[4] = Core.svec(Any[_instantiate_type_in_env(argt, type_signature, static_param_values) for argt in x.args[4]]...)
        elseif head === :foreigncall
            @assert !isa(type_signature, UnionAll) || !isempty(static_param_values)
            for i = 1:length(x.args)
                if i == 2
                    x.args[2] = _instantiate_type_in_env(x.args[2], type_signature, static_param_values)
                elseif i == 3
                    x.args[3] = Core.svec(Any[_instantiate_type_in_env(argt, type_signature, static_param_values) for argt in x.args[3]]...)
                elseif i == 4
                    @assert isa(x.args[4], Int)
                elseif i == 5
                    @assert isa((x.args[5]::QuoteNode).value, Symbol)
                else
                    x.args[i] = _partially_inline!(x.args[i], slot_replacements,
                                                   type_signature, static_param_values,
                                                   slot_offset, statement_offset,
                                                   boundscheck)
                end
            end
        elseif head === :boundscheck
            if boundscheck === :propagate
                return x
            elseif boundscheck === :off
                return false
            else
                return true
            end
        elseif head === :gotoifnot
            x.args[1] = _partially_inline!(x.args[1], slot_replacements, type_signature,
                                           static_param_values, slot_offset,
                                           statement_offset, boundscheck)
            x.args[2] += statement_offset
        elseif head === :enter
            x.args[1] += statement_offset
        elseif !is_meta_expr_head(head)
            partially_inline!(x.args, slot_replacements, type_signature, static_param_values,
                              slot_offset, statement_offset, boundscheck)
        end
    end
    return x
end

_instantiate_type_in_env(x, spsig, spvals) = ccall(:jl_instantiate_type_in_env, Any, (Any, Any, Ptr{Any}), x, spsig, spvals)

is_meta_expr_head(head::Symbol) = (head === :inbounds || head === :boundscheck || head === :meta || head === :loopinfo)

end # module

max_issues_repo_name: zhuangchaoxi/julia

max_forks_repo_name: zhuangchaoxi/julia


text:
function findhex(hexID::UInt16, hexIDs::Array{UInt16,1})
  k = 0
  @inbounds while k < lastindex(hexIDs)
    k += 1
    if hexID == getindex(hexIDs, k)
      return k
    end
  end
  k = -1
  return k
end

function win32_cfile!( fname::String,
                        hex_bytes::Array{UInt8,1},
                        hexIDs::Array{UInt16,1},
                        S::SeisData,
                        fc::Array{Float32,1},
                        hc::Array{Float32,1},
                        nx_new::Int64
                      )

  open(fname, "r") do cf_io
    while !eof(cf_io)
      chan_line = readline(cf_io)
      occursin(r"^\s*(?:#|$)", chan_line) && continue

      # chan_info fills S
      chan_info = String.(split(chan_line))

      # Assign identifying info to placeholder arrays
      hex2bytes!(hex_bytes, chan_info[1])
      push!(hexIDs, reinterpret(UInt16, hex_bytes)[1])
      push!(fc, one(Float32)/parse(Float32, chan_info[10]))
      push!(hc, parse(Float32, chan_info[11]))

      # Create new channel in S from chan_info
      loc           = GeoLoc()
      loc.lat       = parse(Float64, chan_info[14])
      loc.lon       = parse(Float64, chan_info[15])
      loc.el        = parse(Float64, chan_info[16])

      C = SeisChannel()
      setfield!(C, :id, string(chan_info[4], ".", chan_info[5]))
      if length(chan_info) > 18
        setfield!(C, :name, chan_info[19])
      end
      setfield!(C, :units, chan_info[9])
      setfield!(C, :gain,  Float64(parse(Float32, chan_info[13]) /
                            (parse(Float32, chan_info[8]) *
                              10.0f0^(parse(Float32, chan_info[12]) / 20.0f0))))
      setfield!(C, :loc, loc)
      setfield!(C, :x, Array{Float32,1}(undef, nx_new))
      D = getfield(C, :misc)
      D["lineDelay"]  = parse(Float32, chan_info[3]) / 1000.0f0
      D["pCorr"]      = parse(Float32, chan_info[17])
      D["sCorr"]      = parse(Float32, chan_info[18])
      push!(S, C)
    end
  end
  return nothing
end

@doc """
    S = readwin32(dfilestr, cfilestr)

Read all win32 data files matching pattern `dfilestr` into SeisData object `S`,
with channel info stored in channel files matching pattern `cfilestr`. Both
file patterns accept wild cards.

    readwin32!(S, dfilestr, cfilestr)

As above, appending data to an existing SeisData object `S`.

!!! warning

    Using multiple channel files applies no redundancy checks of any kind.
""" readwin32
function readwin32(dfilestr::String, cfilestr::String;
                    jst::Bool=true,
                    nx_add::Int64=KW.nx_add,
                    nx_new::Int64=KW.nx_new,
                    v::Int64=KW.v
                    )
  S = SeisData()

  # Parse channel file(s)
  hex_bytes = Array{UInt8,1}(undef,2)
  hexIDs    = Array{UInt16,1}(undef, 0)
  fc        = Array{Float32,1}(undef,0)
  hc        = Array{Float32,1}(undef,0)
  if safe_isfile(cfilestr)
    win32_cfile!(cfilestr, hex_bytes, hexIDs, S, fc, hc, nx_new)
  else
    cfiles = ls(cfilestr)
    @inbounds for cfile in cfiles
      v > 1 && println("Reading channel file ", fname)
      win32_cfile!(cfile, hex_bytes, hexIDs, S, fc, hc, nx_new)
    end
  end
  L = lastindex(hexIDs)

  # Parse data files
  files = ls(dfilestr)
  nf = length(files)
  jst_const = (jst == true ? 32400000000 : 0)
  date_hex  = zeros(UInt8, 8)
  date_arr  = zeros(Int64, 6)
  checkbuf!(BUF.int32_buf, 100)
  checkbuf!(BUF.buf, 1000)
  buf       = getfield(BUF, :buf)
  x         = view(BUF.int32_buf, 1:100)

  # Preallocate arrays
  sums      = zeros(Int64, L)
  seisN     = zeros(Int64, L)
  OldTime   = zeros(Int64, L)
  locID     = Array{String,1}(undef, L)
  xi        = zeros(Int64, L)
  gapStart  = Array{Array{Int64,1},1}(undef, L)
  gapEnd    = Array{Array{Int64,1},1}(undef, L)

  @inbounds for fname in files
    v > 0 && println("Processing ", fname)
    open(fname, "r") do fid
      skip(fid, 4)
      while !eof(fid)

        # Start time
        read!(fid, date_hex)
        t_new = datehex2μs!(date_arr, date_hex)
        skip(fid, 4)

        # Bytes to read
        lsecb = Int64(ntoh(read(fid, UInt32)))
        τ = 0

        while τ < lsecb
          orgID = read(fid, UInt8)
          netID = read(fid, UInt8)
          hexID = read(fid, UInt16)
          k = findhex(hexID, hexIDs)
          V = ntoh(read(fid, UInt16))
          C = Int64(V >> 12)
          N = Int64(V & 0x0fff)
          Nh = N

          # Increment bytes read (this file), decrement N if not 4-bit
          if C == 0
            B = div(N, 2)
          else
            N -= 1
            B = C*N
          end
          τ += 10 + B
          ii = getindex(xi, k)

          # Create new channel
          if ii == 0
            nx = 60*Nh*nf
            if nx != lastindex(S.x[k])
              resize!(S.x[k], nx)
            end
            setindex!(getfield(S, :fs), Float64(Nh), k)
            t = getindex(getfield(S, :t), k)
            T = Array{Int64,2}(undef, 2, 2)
            setindex!(T, one(Int64), 1)
            setindex!(T, Int64(nx), 2)
            setindex!(T, t_new-jst_const, 3)
            setindex!(T, zero(Int64), 4)
            setindex!(getfield(S, :t), T, k)
            setindex!(gapEnd, Int64[], k)
            setindex!(gapStart, Int64[], k)
            setindex!(locID, bytes2hex([orgID & 0xff | (netID << 4) & 0xf0]), k)
            D = getindex(getfield(S, :misc), k)
            D["orgID"] = orgID
            D["netID"] = netID
            D["hexID"] = hexID
            D["locID"] = locID[k]
          end

          # Parse data
          x[1] = bswap(read(fid, Int32))
          if C == 0
            readbytes!(fid, buf, B)
            fillx_i4!(x, buf, B, 1)
          elseif C == 1
            readbytes!(fid, buf, N)
            fillx_i8!(x, buf, N, 1)
          elseif C == 2
            readbytes!(fid, buf, 2*N)
            fillx_i16_be!(x, buf, N, 1)
          elseif C == 3
            readbytes!(fid, buf, 3*N)
            fillx_i24_be!(x, buf, N, 1)
          else
            readbytes!(fid, buf, 4*N)
            fillx_i32_be!(x, buf, N, 1)
          end

          # Account for time gaps
          t_old = getindex(OldTime, k)
          gap = t_new - t_old
          if (gap > 1000000) && (t_old > 0)
            gl = div((gap - 1000000),1000000)
            (v > 0) && @warn(string("Time gap detected! (channel ", hexID, ", length ", @sprintf("%.1f",gl), "s, begin ", u2d(t_old*1.0e-6)))
            P = Nh*gl
            push!(gapStart[k], ii + 1)
            push!(gapEnd[k], ii + P)
            ii += P
          end

          y = getindex(getfield(S, :x), k)
          xa = first(x)
          j = 1
          while j < Nh
            j += 1
            xa += getindex(x,j)
            setindex!(x, xa, j)
          end
          copyto!(y, ii+1, x, 1, Nh)

          # Update counters
          OldTime[k] = t_new
          sums[k] += x[Nh]
          seisN[k] += Nh
          xi[k] = ii + Nh
        end
      end
      close(fid)
    end
  end

  # Post-process
  src = dfilestr
  κ = findall(xi.==0)
  i = 0
  @inbounds while i < S.n
    i += 1
    i in κ && continue
    χ = getindex(S.x, i)

    # Ensure we aren't overcompensating
    lastindex(χ) == getindex(xi, i) || resize!(χ, xi[i])

    # Get resp for passive velocity sensors
    fci = getindex(fc, i)
    hci = getindex(hc, i)
    resp = fctoresp(fci, hci)
    setindex!(getfield(S, :resp), resp, i)

    # There will be issues here. Japanese files use NIED or local station
    # names, which don't necessarily use international station or network codes.
    # For an example of the (lack of) correspondence see
    # http://data.sokki.jmbsc.or.jp/cdrom/seismological/catalog/appendix/apendixe.htm
    (net, sta, cha) = split(S.id[i], ".", limit=3, keepempty=true)
    # Band code
    bb = getbandcode(getindex(getfield(S, :fs), i), fc=fci)
    # Channel code
    cc = String(cha)[1:1]

    if cc == "U"
      cc = "Z"
      S.loc[i].inc = 180.0
    elseif cc == "N"
      S.loc[i].inc = 90.0
    elseif cc == "E"
      S.loc[i].az = 90.0
      S.loc[i].inc = 90.0
    end
    id = string(net, ".", sta, ".", locID[i], ".", bb, "H", cc)
    setindex!(getfield(S, :id), id, i)
    setindex!(getfield(S, :src), src, i)

    # Fill gaps with mean of data
    J = length(gapStart[i])
    if J > 0
      μ = sums[i] / seisN[i]
      gs = gapStart[i]
      ge = gapEnd[i]
      n = 0
      while n < J
        n += 1
        fill!(view(χ, gs[n]:ge[n]), μ)
      end
    end
  end
  note!(S, "+src: " * src)
  note!(S, "channel file: " * cfilestr)
  if !isempty(κ)
    deleteat!(S, κ)
    v > 0 && println("Deleted ", length(κ), " empty channels after read.")
  end
  return S
end

@doc (@doc readwin32)
function readwin32!(S::SeisData, dfilestr::String, cfilestr::String;
                    jst::Bool=true,
                    nx_add::Int64=KW.nx_add,
                    nx_new::Int64=KW.nx_new,
                    v::Int64=KW.v
                    )
  U = readwin32(dfilestr, cfilestr, jst=jst, nx_add=nx_add, nx_new=nx_new)
  append!(S,U)
  return nothing
end


text:
# Autogenerated wrapper script for Qt5Base_jll for x86_64-unknown-freebsd
export libqt5concurrent, libqt5core, libqt5dbus, libqt5gui, libqt5network, libqt5opengl, libqt5printsupport, libqt5sql, libqt5test, libqt5widgets, libqt5xml

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt5Base")
JLLWrappers.@declare_library_product(libqt5concurrent, "libQt5Concurrent.so.5")
JLLWrappers.@declare_library_product(libqt5core, "libQt5Core.so.5")
JLLWrappers.@declare_library_product(libqt5dbus, "libQt5DBus.so.5")
JLLWrappers.@declare_library_product(libqt5gui, "libQt5Gui.so.5")
JLLWrappers.@declare_library_product(libqt5network, "libQt5Network.so.5")
JLLWrappers.@declare_library_product(libqt5opengl, "libQt5OpenGL.so.5")
JLLWrappers.@declare_library_product(libqt5printsupport, "libQt5PrintSupport.so.5")
JLLWrappers.@declare_library_product(libqt5sql, "libQt5Sql.so.5")
JLLWrappers.@declare_library_product(libqt5test, "libQt5Test.so.5")
JLLWrappers.@declare_library_product(libqt5widgets, "libQt5Widgets.so.5")
JLLWrappers.@declare_library_product(libqt5xml, "libQt5Xml.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt5concurrent,
        "lib/libQt5Concurrent.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5core,
        "lib/libQt5Core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5dbus,
        "lib/libQt5DBus.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5gui,
        "lib/libQt5Gui.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5network,
        "lib/libQt5Network.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5opengl,
        "lib/libQt5OpenGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5printsupport,
        "lib/libQt5PrintSupport.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5sql,
        "lib/libQt5Sql.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5test,
        "lib/libQt5Test.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5widgets,
        "lib/libQt5Widgets.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xml,
        "lib/libQt5Xml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using GeometryPrimitives, StaticArrays, LinearAlgebra, Test
using Random: MersenneTwister
using Statistics: mean

const rtol = Base.rtoldefault(Float64)
const one⁻ = 1 - rtol  # scale factor slightly less than 1
const one⁺ = 1 + rtol  # scare factor slightly greater than 1
const one⁻⁻, one⁺⁺ = 0.9, 1.1  # (scale factor less than 1, scale factor greater than 1)

Base.isapprox(a::Tuple, b::Tuple; kws...) = all(p -> isapprox(p...; kws...), zip(a,b))
const rng = MersenneTwister(0) # test with reproducible pseudorandom numbers

randnb(lb::Real,ub::Real) = randn(rng)*(ub-lb) + 0.5*(lb+ub)
randnb(lb::AbstractVector,ub::AbstractVector) = map(randnb, lb, ub)
function inbounds(x,lb,ub)
    length(x) == length(lb) == length(lb) || throw(BoundsError())
    for i in eachindex(x)
        @inbounds lb[i] ≤ x[i] ≤ ub[i] || return false
    end
    return true
end

"check the bounding box of s with randomized trials"
function checkbounds(s::Shape{N}, ntrials=10^4) where {N}
    lb,ub = bounds(s)
    for i = 1:ntrials
        x = randnb(lb,ub)
        x ∉ s || inbounds(x,lb,ub) || return false  # return false if s - (bounding box) is nonempty
    end
    return true
end

function checktree(t::KDTree{N}, slist::Vector{<:Shape{N}}, ntrials=10^3) where {N}
    lb = SVector{N}(fill(Inf,N))
    ub = SVector{N}(fill(-Inf,N))
    for i in eachindex(slist)
        lbi,ubi = bounds(slist[i])
        lb = min.(lb,lbi)
        ub = max.(ub,ubi)
    end
    for i = 1:ntrials
        x = randnb(lb,ub)
        st = findfirst(x, t)
        sl = findfirst(x, slist)
        st == sl || return false
    end
    return true
end

@testset "GeometryPrimitives" begin

include("ball.jl")
include("cuboid.jl")
include("ellipsoid.jl")
include("polygon.jl")
include("sector.jl")
include("cylinder.jl")
include("polygonal_prism.jl")
include("sectoral_prism.jl")
include("kdtree.jl")
include("periodize.jl")
include("vxlcut.jl")

end


text:
const NO_MATCH = ([0 0], -1)
const NO_RANGE = (0:0, 0:0)

type Matches
  mesh_indices::Array{Int64, 1}
  src_points::Array
  dst_points::Array
end

"""
`GET_MAX_XC_VECTOR` - Find the cross correlation peak between two images

`[i, j, r_value], correlogram = get_max_xc_vector(A, B)`

* A: 2D array representing the first image
* B: 2D array representing the second image
* i: i displacement from the correlogram center to its peak
* j: j displacement from the correlogram center to its peak
* r_value: the value of the correlogram at its peak
* correlogram: 2D array representing the cross corelation between A & B
"""
function get_max_xc_vector(A, B)
  if std(A) == 0 || std(B) == 0
    return NO_MATCH, 0
  end
  xc = normxcorr2(A, B)
  r_max = maximum(xc)
  if isnan(r_max)
    return NO_MATCH, 0
  end
  rad = round(Int64, (size(xc, 1) - 1)/ 2)
  ind = findfirst(r_max .== xc)
  x1 = size(xc, 1)
  x2 = size(xc, 2)
  if ind == 0
    return NO_MATCH, 0
  end
  (i_max, j_max) = (rem(ind, size(xc, 1)), cld(ind, size(xc, 1)))
  if i_max == 0
    i_max = size(xc, 1)
  end
  return [i_max-1-rad j_max-1-rad], r_max, xc
end

"""
`BLOCKMATCH` - Find block matches between two images at node locations.

`matches = get_blockmatches(nodes, src_img, dst_img, src_offset, dst_offset, params)`

* nodes: Nx2 array of nodes for locations in an image where blockmatches will be
    conducted. Nodes are expected to be in global coordinates (not local to the
    image's pixel space).
* src_img: 2D array representing the image that will be moving when warped
* dst_img: 2D array representing the image that will stay fixed (aka template)
* src_offset: 2-element array representing the location of src_img's [0,0] pixel
    in the global space.
* dst_offset: 2-element array representing the location of dst_img's [0,0] pixel
    in the global space.
* params: Dict object containing elements defining:
  * block_size: radius from a src node that will be sliced from src_img and used
      in the cross correlation
  * search_r: additional radius beyond block_size that will be sliced from
      dst_img and used in the cross correlation
  * min_r: the minimum threshold of a cross correlation peak for a blockmatch to
      be accepted.

This method is parallelized. Make sure to start Julia session with additional
processors:

`julia -p n`

"""
function blockmatch(nodes, src_img, dst_img, src_offset, dst_offset, params)
  n = size(nodes, 1)
  displacements = Array{Tuple{Array{Float64, 2}, Float64}, 1}(n)
  src_ranges = Array{Tuple{UnitRange{Int64}, UnitRange{Int64}}, 1}(n)
  dst_ranges = Array{Tuple{UnitRange{Int64}, UnitRange{Int64}}, 1}(n)

  block_size = params["block_size"]
  search_r = params["search_r"]
  min_r = params["min_r"]
  b_rad = block_size + search_r

  src_sz = size(src_img)
  dst_sz = size(dst_img)

  for i in 1:n
    pt = nodes[i]
    src_ranges[i] = get_range(src_sz, pt, src_offset, block_size)
    dst_ranges[i] = get_range(dst_sz, pt, dst_offset, b_rad)
  end

  num_procs = nprocs()
  k = 1
  nexti() = (i=k; k+=1; i)
  @sync begin
  for p in 1:num_procs
    if p != myid() || num_procs == 1
      @async begin
      while true
        i = nexti()
        if i > n
          break
        end
        if src_ranges[i] == NO_RANGE || dst_ranges[i] == NO_RANGE
          displacements[i] = NO_MATCH
          continue
        end

        src = src_img[src_ranges[i]...]
        dst = dst_img[dst_ranges[i]...]

        max_vect_xc = remotecall_fetch(get_max_xc_vector, p, src,  dst)
        displacements[i] = max_vect_xc[1:2]
      end
      end
    end
  end
  end

  mesh_indices = []
  src_points = []
  dst_points = []
  for i in 1:n
    v = displacements[i]
    if v == NO_MATCH
      continue
    end
    if v[2] < min_r
      continue
    end
    push!(mesh_indices, i)
    push!(src_points, nodes[i, :])
    push!(dst_points, nodes[i, :] + v[1])
  end
  if length(dst_points) == 0
    return Void
  end
  return Matches(mesh_indices, vcat(src_points...), vcat(dst_points...))
end

"""
Is a point contained within an image perimeter, less a certain radius?
"""
function is_internal(sz, pt, radius)
  above_min_i = pt[1] > radius
  above_min_j = pt[2] > radius
  below_max_i = pt[1] <= (sz[1]-radius)
  below_max_j = pt[2] <= (sz[2]-radius)
  return above_min_i && above_min_j && below_max_i && below_max_j
end

"""
From global point and offset, calculate range of pixels within radius to slice
"""
function get_range(sz, pt, offset, radius)
  if !is_internal(sz, pt-offset, radius)
    return NO_RANGE
  end
  pt_local = pt-offset
  i = round(Int64, ceil(pt_local[1]))
  j = round(Int64, ceil(pt_local[2]))
  i_range = i-radius:i+radius
  j_range = j-radius:j+radius
  return i_range, j_range
end


text:
module PolytopesTests

using Test
using LinearAlgebra
using Gridap.Helpers
using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields
using Gridap.ReferenceFEs

import Gridap.ReferenceFEs: Polytope
import Gridap.ReferenceFEs: get_faces
import Gridap.ReferenceFEs: get_dimranges
import Gridap.ReferenceFEs: get_vertex_coordinates
import Base: ==

struct MockVertex <: Polytope{0} end

get_faces(p::MockVertex) = [[1]]

get_dimranges(p::MockVertex) = [1:1]

function Polytope{0}(p::MockVertex,faceid::Integer)
  @assert faceid == 1
  p
end

function (==)(a::MockVertex,b::MockVertex)
  true
end

get_vertex_coordinates(p::MockVertex) = [zero(Point{0,Float64})]

struct MockSegment <: Polytope{1} end

get_faces(p::MockSegment) = [[1],[2],[1,2,3]]

get_dimranges(p::MockSegment) = [1:2,3:3]

function Polytope{0}(p::MockSegment,faceid::Integer)
  @assert faceid in 1:2
  MockVertex()
end

function Polytope{1}(p::MockSegment,faceid::Integer)
  @assert faceid == 1
  p
end

function (==)(a::MockSegment,b::MockSegment)
  true
end

get_vertex_coordinates(p::MockSegment) = Point{1,Float64}[(0),(1)]

struct MockQuad <: Polytope{2} end

function get_faces(p::MockQuad)
  [
    [1],[2],[3],[4],
    [1,2,5],[3,4,6],[1,3,7],[2,4,8],
    [1,2,3,4,5,6,7,8,9]
  ]
end

function get_dimranges(p::MockQuad)
  [1:4,5:8,9:9]
end

function Polytope{0}(p::MockQuad,faceid::Integer)
  @assert faceid in 1:4
  MockVertex()
end

function Polytope{1}(p::MockQuad,faceid::Integer)
  @assert faceid in 1:4
  MockSegment()
end

function Polytope{2}(p::MockQuad,faceid::Integer)
  @assert faceid == 1
  p
end

function (==)(a::MockQuad,b::MockQuad)
  true
end

get_vertex_coordinates(p::MockQuad) = Point{2,Float64}[(0,0),(1,0),(0,1),(1,1)]

v = MockVertex()
test_polytope(v)
@test get_faces(v,0,0) == [[1]]
@test get_facedims(v) == [0]
@test get_offsets(v) == [0]
@test num_facets(v) == 0
@test num_edges(v) == 0

s = MockSegment()
test_polytope(s)
@test get_faces(s,0,0) == [[1], [2]]
@test get_faces(s,1,0) == [[1, 2]]
@test get_faces(s,0,1) == [[1], [1]]
@test get_faces(s,1,1) == [[1]]
@test get_facedims(s) == [0,0,1]
@test get_offsets(s) == [0,2]
@test num_facets(s) == 2
@test num_edges(s) == 1

q = MockQuad()
test_polytope(q)
@test get_faces(q,0,0) == [[1], [2], [3], [4]]
@test get_faces(q,1,0) == [[1, 2], [3, 4], [1, 3], [2, 4]]
@test get_faces(q,2,0) == [[1, 2, 3, 4]]
@test get_faces(q,0,1) == [[1, 3], [1, 4], [2, 3], [2, 4]]
@test get_faces(q,1,1) == [[1], [2], [3], [4]]
@test get_faces(q,2,1) == [[1, 2, 3, 4]]
@test get_faces(q,0,2) == [[1], [1], [1], [1]]
@test get_faces(q,1,2) == [[1], [1], [1], [1]]
@test get_faces(q,2,2) == [[1]]
@test get_facedims(q) == [0,0,0,0,1,1,1,1,2]
@test get_offsets(q) == [0,4,8]
@test num_facets(q) == 4
@test num_edges(q) == 4

@test get_face_vertices(q,1) == [[1, 2], [3, 4], [1, 3], [2, 4]]
@test get_face_vertices(q) == [[1], [2], [3], [4], [1, 2], [3, 4], [1, 3], [2, 4], [1, 2, 3, 4]]

d = 1
@test get_dimrange(q,d) == get_dimranges(q)[d+1]

pmin, pmax = get_bounding_box(HEX)
@test pmin == Point(0,0,0)
@test pmax == Point(1,1,1)

end # module


text:
using GeoData, Test, Dates, Statistics
using GeoData: Start, Center, End,
      formatdims, dims, aggregate, upsample, downsample, span

@testset "upsample" begin
    @test upsample(1, 2) == 1
    @test upsample(2, 2) == 3
    @test upsample(3, 2) == 5
    @test upsample(1, 3) == 1
    @test upsample(2, 3) == 4
    @test upsample(3, 3) == 7
end

@testset "downsample" begin
    @test downsample(1, 2) == 1
    @test downsample(3, 2) == 2
    @test downsample(5, 2) == 3
    @test downsample(1, 3) == 1
    @test downsample(4, 3) == 2
    @test downsample(7, 3) == 3
end

data1 = [ 1  2  3  4  5  6 -1
          7  8  9 10 11 12 -1
         13 14 15 16 17 18 -1]
data2 = 2 * data1
data3 = 3 * data1
data4 = 4 * data1
dimz = Lon([30., 40., 50.]; mode=Sampled(Ordered(), Regular(10.0), Points())),
       Lat(LinRange(-10., 20., 7); mode=Sampled(Ordered(), Regular(5.0), Points()))
array1 = GeoArray(data1, dimz)
array2 = GeoArray(data2, dimz)
array1a = GeoArray(data3, dimz)
array2a = GeoArray(data4, dimz)
stack1 = GeoStack(array1, array2; keys=(:array1, :array2))
stack2 = GeoStack(array1a, array2a; keys=(:array1, :array2))
dates = [DateTime(2017), DateTime(2018)]
series = GeoSeries([stack1, stack2], (Ti(dates),));

@testset "Aggregate a dimension" begin
    lat = Lat(LinRange(3, 13, 6);
              mode=Sampled(Ordered(), Regular(2.0), Intervals(Start())))
    aglat = aggregate(Start(), lat, 3)
    @test span(mode(aglat)) == Regular(6.0)
    @test disaggregate(Start(), aglat, 3) == lat

    aglon = aggregate(Start(), dimz[1], 3)
    @test step(mode(aglon)) === 30.0
    @test val(aglon) == [30.0]
    disaglon = disaggregate(Start(), aglon, 3)
    @test val(disaglon) == val(dimz[1])
    @test step(disaglon) == step(dimz[1])
    @test mode(disaglon) == mode(dimz[1])

    aglat = aggregate(Start(), dimz[2], 3)
    @test step(mode(aglat)) === 15.0
    @test val(aglat) == LinRange(-10.0, 5.0, 2)
    disaglat = disaggregate(Start(), aglat, 3)
    # The last item is lost due to rounding in `aggregate`
    @test val(disaglat) != val(dimz[2])
    @test val(disaglat) === LinRange(-10.0, 15.0, 6)
    @test step(disaglat) == step(dimz[2])
    @test mode(disaglat) == mode(dimz[2])
end

@testset "aggregate and disaggregate at a locus" begin
    @testset "single scale single locus" begin
        scale = 3
        array1_aggstart = aggregate(Start(), array1, scale)
        @test array1_aggstart == [1 4]
        @test length.(dims(array1_aggstart)) == size(array1_aggstart)
        array1_disagstart = disaggregate(Start(), array1_aggstart, scale)
        disaggstart =
            [1 1 1 4 4 4
             1 1 1 4 4 4
             1 1 1 4 4 4]
        @test array1_disagstart == disaggstart
        array1_aggcenter = aggregate(Center(), array1, scale)
        @test array1_aggcenter == [8 11]
        array1_disagcenter = disaggregate(Center(), array1_aggcenter, scale)
        disaggcenter =
            [8 8 8 11 11 11
             8 8 8 11 11 11
             8 8 8 11 11 11]
        @test array1_disagcenter == disaggcenter
        array1_aggend = aggregate(End(), array1, scale)
        @test array1_aggend == [15 18]
        array1_disaggend = disaggregate(End(), array1_aggend, scale)
        @test array1_disaggend ==
            [15 15 15 18 18 18
             15 15 15 18 18 18
             15 15 15 18 18 18]

        stack1_aggstart = aggregate(Start(), stack1, scale)
        @test stack1_aggstart[:array2] == [2 8]
        stack1_disaggstart = disaggregate(Start(), stack1_aggstart, scale)
        @test stack1_disaggstart[:array1] == disaggstart
        @test aggregate(Center(), stack1, scale)[:array2] == [16 22]

        series_aggcenter = aggregate(Center(), series, scale)
        @test series_aggcenter[2][:array2] == [32 44]
        series_disaggcenter = disaggregate(Center(), series_aggcenter, scale)
        @test series_disaggcenter[2][:array2] ==
            [32 32 32 44 44 44
             32 32 32 44 44 44
             32 32 32 44 44 44]
        @test typeof(aggregate(Start(), series, scale)) <: GeoSeries
    end

    @testset "mixed scales" begin
        scale = (3, 2)
        @test aggregate(Start(), array1, scale) == [1 3 5]
        scale = (Lat(2), Lon(3))
        @test aggregate(Center(), array1, scale) == [8 10 12]
        @test aggregate(End(), array1, scale) == [14 16 18]
        A = aggregate(Start(), array1, scale)
        @test length.(dims(A)) == size(A)
    end

    @testset "mixed locus" begin
        scale = 3
        @test aggregate((End(), Start()), array1, 3) == [13 16]
        @test aggregate((End(), Start()), array1, (3, 2)) == [13 15 17]
        A = aggregate((End(), Start()), array1, scale)
        @test length.(dims(A)) == size(A)
    end

    @testset "dim scale" begin
        agg = aggregate(Start(), array1, (Lat(3), Lon(1)))
        @test agg == aggregate(Start(), array1, (1, 3))
        disagg = disaggregate(Start(), agg, (Lat(3), Lon(1)))
        @test disagg ==
            [1  1  1  4  4  4
             7  7  7 10 10 10
            13 13 13 16 16 16]
        agg = aggregate(Start(), array1, (Lon(1), Lat(Near(-4))))
        @test agg == aggregate(Start(), array1, (1, 2))
        @testset "scale 1 dims are unchanged" begin
            @test dims(agg, Lon) === dims(array1, Lon)
        end
    end

end

@testset "Aggregate with a function" begin
    @test aggregate(sum, array1, 3) == [72 99]
    @test aggregate(median, array1, 3) == [8 11]
    @test aggregate(sum, array1, (3, 2)) == [45 57 69]
    A = aggregate(sum, array1, (3, 2))
    @test length.(dims(A)) == size(A)
end

@testset "Aggregate different index modes" begin
    dimz = Band(1:3), Dim{:category}([:a, :b, :c]), X([10, 20, 30, 40])
    a1 = [1 2 3; 4 5 6; 7 8 9]
    A = cat(a1, a1 .+ 10, a1 .+ 20, a1 .+ 30, dims=3)
    da = DimensionalArray(A, dimz)
    @test vec(aggregate(sum, da, (3, 2, 2))) == [114, 354]
end


text:
struct FormatError <: Exception
    msg::String
end

function Base.showerror(out::IO, e::FormatError)
    print(out, "FormatError: ", e.msg)
end

function parse(fmt::String)
    list = Union{String, Field}[]
    auto = 0  # automatic numbering
    str = IOBuffer()
    last = lastindex(fmt)
    i = firstindex(fmt)
    while i ≤ last
        c = fmt[i]
        if c == '{'
            while i + 1 ≤ last && fmt[i] == '{' && fmt[i+1] == '{'
                write(str, '{')
                i += 2
            end
            if i ≤ last && fmt[i] == '{'
                i == last && throw(FormatError("single '{' is not allowed; use '{{' instead"))
                str.size > 0 && push!(list, String(take!(str)))
                field, i, auto = parse_field(fmt, i + 1, auto)
                push!(list, field)
            end
        elseif c == '}'
            while i + 1 ≤ last && fmt[i] == '}' && fmt[i+1] == '}'
                write(str, '}')
                i += 2
            end
            if i ≤ last && fmt[i] == '}' && !(i + 1 ≤ last && fmt[i+1] == '}')
                throw(FormatError("single '}' is not allowed; use '}}' instead"))
            end
        else
            write(str, c)
            i = nextind(fmt, i)
        end
    end
    str.size > 0 && push!(list, String(take!(str)))

    # check fields
    interpolated = missing
    for f in list
        f isa Field || continue
        if ismissing(interpolated)
            interpolated = f.argument isa Expr
        elseif interpolated != (f.argument isa Expr)
            throw(FormatError("mixing interpolated and non-interpolated fields is not allowed"))
        end
    end

    if isempty(list)
        push!(list, "")
    end
    return list
end

function parse_field(fmt::String, i::Int, auto::Int)
    incomplete_field() = throw(FormatError("incomplete field"))
    last = lastindex(fmt)
    arg, i, auto = parse_argument(fmt, i, auto)
    i ≤ last || incomplete_field()
    conv = CONV_UNSPECIFIED
    if fmt[i] == '/'
        i + 1 ≤ last || incomplete_field()
        conv, i = parse_conv(fmt, i + 1)
        i ≤ last || incomplete_field()
    end
    spec = SPEC_DEFAULT
    if fmt[i] == ':'
        i + 1 ≤ last || incomplete_field()
        spec, i, auto = parse_spec(fmt, i + 1, auto)
        i ≤ last || incomplete_field()
    end
    fmt[i] == '}' || throw(FormatError("invalid character $(repr(fmt[i]))"))
    # check consistency of arguments
    if arg isa Expr
        spec.fill      isa Union{Char,         Expr} &&
        spec.width     isa Union{Int, Nothing, Expr} &&
        spec.precision isa Union{Int, Nothing, Expr} ||
        throw(FormatError("inconsistent interpolation of arguments"))
    else
        spec.fill      isa Union{Char,         Positional, Keyword} &&
        spec.width     isa Union{Int, Nothing, Positional, Keyword} &&
        spec.precision isa Union{Int, Nothing, Positional, Keyword} ||
        throw(FormatError("inconsistent interpolation of arguments"))
    end
    return Field(arg, conv, spec), i + 1, auto
end

function parse_argument(s::String, i::Int, auto::Int)
    c = s[i]  # the first character after '{'
    if c == '$'
        i < lastindex(s) && (Base.is_id_start_char(s[i+1]) || s[i+1] == '(') ||
            throw(FormatError("identifier or '(' is expected after '\$'"))
        ast, i = Meta.parse(s, i + 1, greedy = false)
        arg = ast isa Expr ? ast : Expr(:block, ast)
    elseif isdigit(c)
        num, i = parse_digits(s, i)
        num == 0 && throw(FormatError("argument 0 is not allowed; use 1 or above"))
        arg = Positional(num)
    elseif Base.is_id_start_char(c)
        name, i = Meta.parse(s, i, greedy = false)
        arg = Keyword(name)
    else
        auto += 1
        arg = Positional(auto)
    end
    return arg, i, auto
end

function parse_conv(fmt::String, i::Int)
    c = fmt[i]
    if c == 'r'
        return CONV_REPR, i + 1
    elseif c == 's'
        return CONV_STRING, i + 1
    else
        throw(FormatError("invalid conversion character $(repr(c))"))
    end
end

function parse_spec(fmt::String, i::Int, auto::Int)
    # default
    fill = FILL_DEFAULT
    align = ALIGN_UNSPECIFIED
    sign = SIGN_DEFAULT
    altform = ALTFORM_DEFAULT
    zero = ZERO_DEFAULT
    width = WIDTH_UNSPECIFIED
    grouping = GROUPING_UNSPECIFIED
    precision = PRECISION_UNSPECIFIED
    type = TYPE_UNSPECIFIED

    incomplete_argument() = throw(FormatError("incomplete argument"))
    char2align(c) = c == '<' ? ALIGN_LEFT :
                    c == '^' ? ALIGN_CENTER :
                    c == '>' ? ALIGN_RIGHT : @assert false

    # align
    last = lastindex(fmt)
    if fmt[i] == '{'
        # dynamic fill or dynamic width?
        _arg, _i, _auto = parse_argument(fmt, i + 1, auto)
        _i ≤ last && fmt[_i] == '}' || incomplete_argument()
        if _i + 1 ≤ last && fmt[_i+1] ∈ "<^>"
            # it was a dynamic fill
            fill = _arg
            align = char2align(fmt[_i+1])
            auto = _auto
            i = _i + 2
            i ≤ last || @goto END
        end
    elseif fmt[i] != '}' && nextind(fmt, i) ≤ last && fmt[nextind(fmt, i)] ∈ "<^>"
        # fill + align
        fill = fmt[i]
        i = nextind(fmt, i)
        align = char2align(fmt[i])
        i += 1
        i ≤ last || @goto END
    elseif fmt[i] ∈ "<^>"
        # align only
        align = char2align(fmt[i])
        i += 1
        i ≤ last || @goto END
    end

    # sign
    if fmt[i] ∈ "-+ "
        sign = fmt[i] == '-' ? SIGN_MINUS : fmt[i] == '+' ? SIGN_PLUS : SIGN_SPACE
        i += 1
        i ≤ last || @goto END
    end

    # alternative form
    if fmt[i] == '#'
        altform = true
        i += 1
        i ≤ last || @goto END
    end

    # width
    if fmt[i] == '{'
        width, i, auto = parse_argument(fmt, i + 1, auto)
        i ≤ last && fmt[i] == '}' || incomplete_argument()
        i += 1
        i ≤ last || @goto END
    elseif isdigit(fmt[i])
        if fmt[i] == '0' && i + 1 ≤ last && isdigit(fmt[i+1])
            # preceded by zero
            zero = true
            i += 1
        end
        width, i = parse_digits(fmt, i)
        i ≤ last || @goto END
    end

    # grouping
    if fmt[i] == ','
        grouping = GROUPING_COMMA
        i += 1
        i ≤ last || @goto END
    elseif fmt[i] == '_'
        grouping = GROUPING_UNDERSCORE
        i += 1
        i ≤ last || @goto END
    end

    # precision
    if fmt[i] == '.'
        i += 1
        i ≤ last || @goto END
        if fmt[i] == '{'
            precision, i, auto = parse_argument(fmt, i + 1, auto)
            i ≤ last && fmt[i] == '}' || incomplete_argument()
            i += 1
            i ≤ last || @goto END
        elseif isdigit(fmt[i])
            precision, i = parse_digits(fmt, i)
            i ≤ last || @goto END
        else
            throw(FormatError("unexpected $(repr(fmt[i])) after '.'"))
        end
    end

    # type
    if fmt[i] ∈ "dXxoBbcpsFfEeGgAa%"
        type = fmt[i]
        i += 1
        i ≤ last || @goto END
    end

    @label END
    return Spec(fill, align, sign, altform, zero, width, grouping, precision, type), i, auto
end

function parse_digits(s::String, i::Int)
    n = 0
    last = lastindex(s)
    while i ≤ last && isdigit(s[i])
        n′ = 10n + Int(s[i] - '0')
        n′ ≥ n || throw(FormatError("number overflows"))
        n = n′
        i += 1
    end
    return n, i
end


text:
# Sample script for mesh plotting.
#
# Represent cell centers as markers.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot

file = "test/data/bulk.amr.vlsv"
meta = load(file)

fig = plt.figure()

ax = fig.add_subplot(projection="3d")

plotmesh(meta, marker="+")


text:
    #
    # Sin TESTEAR. TTT-D
    #

    function likelihood_teams_draw(g::Game)
        r = g.result
        o = sortperm(r)
        t = [team_messages(performance(g,o[e]), Ninf, Ninf, Ninf) for e in 1:length(g)]
        u = [draw_messages(draw_performance(g,o[e]), draw_performance(g,o[e]) + t[e].prior, Ninf, Ninf) for e in 1:length(g)]
        tie = [r[o[e]]==r[o[e+1]] for e in 1:length(g)-1]
        d = [(diff_messages(Ninf, Ninf), diff_messages(Ninf, Ninf),) for e in 1:length(tie) ]
        step = (Inf, Inf)::Tuple{Float64,Float64}; iter = 0::Int64

        while (step > 1e-6) & (iter < 20)
            step = (0., 0.)
            for e in 1:length(d)#e=2
                if !tie[e]
                    #TODO: crear par\'ametros por defecto para trunc()
                    d[e][1].prior = posterior_win(t[e]) - posterior_lose(u[e+1])
                    d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                    u[e+1].likelihood_lose =  posterior_win(t[e]) - d[e][1].likelihood
                else
                    d[e][1].prior = posterior_win(u[e]) - posterior_lose(t[e+1])
                    d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                    t[e+1].likelihood_lose = posterior_win(u[e]) - d[e][1].likelihood
                    d[e][2].prior = posterior_win(u[e+1]) - posterior_lose(t[e])
                    d[e][2].likelihood = trunc(d[e][2].prior,0.,false)/d[e][2].prior
                    u[e+1].likelihood_win = posterior_lose(t[e]) + d[e][2].likelihood
                end
                t[e+1].likelihood_draw = likelihood(u[e+1]) - u[e+1].prior
            end
            d21_likelihood = d[2][1].likelihood
            for e in length(d):-1:1
                if !tie[e]
                    d[e][1].prior = posterior_win(t[e]) - posterior_lose(u[e+1])
                    d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                    t[e].likelihood_win = posterior_lose(u[e+1]) + d[e][1].likelihood
                else
                    d[e][1].prior = posterior_win(u[e]) - posterior_lose(t[e+1])
                    d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                    u[e].likelihood_win = posterior_lose(t[e+1]) + d[e][1].likelihood
                    d[e][2].prior = posterior_win(u[e+1]) - posterior_lose(t[e])
                    d[e][2].likelihood = trunc(d[e][2].prior,0.,false)/d[e][2].prior
                    t[e].likelihood_lose = posterior_win(u[e+1]) - d[e][2].likelihood
                end
                u[e].prior_team = posterior_draw(t[e]) + u[e].prior
            end
            step = max(step,delta(d[2][1].likelihood,d21_likelihood))
            iter += 1
        end
        if length(d)==1
            e=1
            if !tie[e]
                d[e][1].prior = posterior_win(t[e]) - posterior_lose(u[e+1])
                d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                u[e+1].likelihood_lose =  posterior_win(t[e]) - d[e][1].likelihood
                t[e+1].likelihood_draw = likelihood(u[e+1]) - u[e+1].prior
                t[e].likelihood_win = posterior_lose(u[e+1]) + d[e][1].likelihood
            else
                while (step > 1e-6) & (iter < 10)
                    d11_likelihood = d[e][1].likelihood

                    u[e].prior_team = posterior_draw(t[e]) + u[e].prior

                    d[e][1].prior = posterior_win(u[e]) - posterior_lose(t[e+1])
                    d[e][1].likelihood = trunc(d[e][1].prior,0.,false)/d[e][1].prior
                    u[e].likelihood_win = posterior_lose(t[e+1]) + d[e][1].likelihood
                    t[e+1].likelihood_lose = posterior_win(u[e]) - d[e][1].likelihood

                    d[e][2].prior = posterior_win(u[e+1]) - posterior_lose(t[e])
                    d[e][2].likelihood = trunc(d[e][2].prior,0.,false)/d[e][2].prior
                    u[e+1].likelihood_win = posterior_lose(t[e]) + d[e][2].likelihood
                    t[e].likelihood_lose = posterior_win(u[e+1]) - d[e][2].likelihood

                    t[e+1].likelihood_draw = likelihood(u[e+1]) - u[e+1].prior
                    t[e].likelihood_draw = likelihood(u[e]) - u[e].prior

                    step = delta(d[1][1].likelihood_win,d11_likelihood)
                end
            end
        end
        return [ likelihood(t[o[e]]) for e in 1:length(t)]
    end


text:
#################################################################################
#
#   Factor.jl : Factorization
#
#################################################################################

export Fac, factors, unit

################################################################################
#
#   Type
#
################################################################################

mutable struct Fac{T <: RingElem}
   unit::T
   fac::Dict{T, Int}

   function Fac{T}() where {T}
     f = new()
     f.fac = Dict{T, Int}()
     return f
   end
end

function Fac(u::T, d::Dict{T, Int}) where {T}
   f = Fac{T}()
   f.unit = u
   f.fac = d
   return f
end

doc"""
    unit(a::Fac{T}) -> T

> Returns the unit of the factorization.
"""
unit(a::Fac) = a.unit

#primes(a::Fac) = collect(keys(a.fac))

################################################################################
#
#   Syntax sugar
#
################################################################################

doc"""
    in(a::T, b::Fac{T})

> Test whether `a` is a factor of `b`.
"""
function Base.in(a::T, b::Fac{T}) where {T}
   a in keys(b.fac)
end

doc"""
    getindex(a::Fac{T}, b::T) -> Int

> If `b` is a factor of `a`, the corresponding exponent is returned. Otherwise
> an error is thrown.
"""
function getindex(a::Fac{T}, b::T) where {T}
  if haskey(a.fac, b)
    return a.fac[b]
  else
    error("$b is not a factor of $a")
  end
end

doc"""
    setindex!(a::Fac{T}, c::Int, b::T)

> If `b` is a factor of `a`, the corresponding entry is set to c.
"""
function setindex!(a::Fac{T}, c::Int, b::T) where {T}
  if haskey(a.fac, b)
    error("$b is already set (to $(a[b]))")
  else
    setindex!(a.fac, c, b)
  end
end

################################################################################
#
#   String I/O
#
################################################################################

function Base.show(io::IO, a::Fac)
  if isdefined(a, :unit)
    print(io, "$(a.unit)")
  else
    print(io, "unit not set")
  end
  for (p, e) in a.fac
    bracket = needs_parentheses(p)
    print(io, " * ")
    bracket ? print(io, "($p)") : print(io, "$p")
    a[p] == 1 ? nothing : print(io, "^$e")
  end
end

################################################################################
#
#   Make Factor objects iterable
#
################################################################################

Base.start(a::Fac) = Base.start((a.fac))

Base.done(a::Fac, state) = Base.done((a.fac), state)

Base.next(a::Fac, state) = Base.next((a.fac), state)

doc"""
    length(a::Fac) -> Int

> Returns the number of factors of `a`, not including the unit.
"""
Base.length(a::Fac) = Base.length(a.fac)

################################################################################
#
#   Special functions for Fac{fmpz}
#
################################################################################

Base.in(b::Integer, a::Fac{fmpz}) = Base.in(fmpz(b), a)

Base.getindex(a::Fac{fmpz}, b::Integer) = getindex(a, fmpz(b))


text:
module JungleHelperThinCrumbleBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Crumble Blocks (Thin) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.CrumbleBlock,
        "rectangle",
        Dict{String, Any}(
          "texture" => "JungleHelper/thin"
        )
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
struct NonDecomposedModel <: AbstractNLPEvaluator
    problem::ProxALProblem
    modelinfo::ModelInfo
    algparams::AlgParams
    opfdata::OPFData
    rawdata::RawData
    space::AbstractBackend
end

"""
    NonDecomposedModel(
        case_file::String,
        load_file::String,
        modelinfo::ModelInfo,
        algparams::AlgParams,
        space::AbstractBackend=JuMPBackend(),
        time_horizon_start = 1,
    )

Instantiate non-decomposed multi-period ACOPF instance
specified in `case_file` with loads in `load_file` with model parameters
`modelinfo` and algorithm parameters `algparams`.
The problem is defined over the horizon
[`time_horizon_start`, `modelinfo.num_time_periods`]

"""
function NonDecomposedModel(
    case_file::String, load_file::String,
    modelinfo::ModelInfo,
    algparams::AlgParams,
    space::AbstractBackend=JuMPBackend(),
    time_horizon_start=1
)
    rawdata = RawData(case_file, load_file)
    opfdata = opf_loaddata(
        rawdata;
        time_horizon_start = time_horizon_start,
        time_horizon_end = modelinfo.num_time_periods,
        load_scale = modelinfo.load_scale,
        ramp_scale = modelinfo.ramp_scale,
        corr_scale = modelinfo.corr_scale
    )

    # ctgs_arr = deepcopy(rawdata.ctgs_arr)
    problem = ProxALProblem(opfdata, rawdata, modelinfo, algparams, space, nothing)
    return NonDecomposedModel(problem, modelinfo, algparams, opfdata, rawdata, space)
end

"""
    optimize!(nlp::NonDecomposedModel)

Solve problem using the `nlp` evaluator
of the nondecomposed model.

"""
function optimize!(nlp::NonDecomposedModel)

    opfmodel = opf_model_nondecomposed(nlp.opfdata, nlp.rawdata, nlp.modelinfo, nlp.algparams)
    return opf_solve_nondecomposed(opfmodel, nlp.opfdata, nlp.modelinfo, nlp.algparams)
end

function opf_model_nondecomposed(opfdata::OPFData, rawdata::RawData, modelinfo::ModelInfo, algparams::AlgParams)
    opfmodel = JuMP.Model(algparams.optimizer)
    opf_model_add_variables(opfmodel, opfdata, modelinfo, algparams)
    opf_model_add_block_constraints(opfmodel, opfdata, rawdata, modelinfo)
    obj_expr = compute_objective_function(opfmodel, opfdata, modelinfo, algparams)

    if algparams.mode == :lyapunov_bound
        lyapunov_expr = compute_quadratic_penalty(opfmodel, opfdata,  modelinfo, algparams)
        if !algparams.decompCtgs
            opf_model_add_ctgs_linking_constraints(opfmodel, opfdata, modelinfo)
        end
    else
        lyapunov_expr = 0
        opf_model_add_time_linking_constraints(opfmodel, opfdata, modelinfo)
        opf_model_add_ctgs_linking_constraints(opfmodel, opfdata, modelinfo)
    end

    @objective(opfmodel, Min, obj_expr + lyapunov_expr)

    return opfmodel
end

function opf_solve_nondecomposed(opfmodel::JuMP.Model, opfdata::OPFData,
                                 modelinfo::ModelInfo,
                                 algparams::AlgParams)
    JuMP.optimize!(opfmodel)
    status = termination_status(opfmodel)
    if status ∉ MOI_OPTIMAL_STATUSES
        (algparams.verbose > 0) &&
            @warn("$(algparams.mode) model not solved to optimality. status: $status")
        return nothing
    end

    x = OPFPrimalSolution(opfdata, modelinfo)
    x.Pg .= value.(opfmodel[:Pg])
    x.Qg .= value.(opfmodel[:Qg])
    x.Vm .= value.(opfmodel[:Vm])
    x.Va .= value.(opfmodel[:Va])
    x.ωt .= value.(opfmodel[:ωt])
    x.St .= value.(opfmodel[:St])
    x.Zt .= value.(opfmodel[:Zt])
    x.Sk .= value.(opfmodel[:Sk])
    x.Zk .= value.(opfmodel[:Zk])
    if modelinfo.allow_constr_infeas
        x.sigma_real .= value.(opfmodel[:sigma_real])
        x.sigma_imag .= value.(opfmodel[:sigma_imag])
        x.sigma_lineFr .= value.(opfmodel[:sigma_lineFr])
        x.sigma_lineTo .= value.(opfmodel[:sigma_lineTo])
    end

    result = Dict()
    result["primal"] = x
    result["objective_value_" * String(algparams.mode)] = objective_value(opfmodel)
    result["solve_time"] = JuMP.solve_time(opfmodel)

    return result
end

max_stars_repo_name: exanauts/proxALM

max_issues_repo_name: exanauts/proxALM

max_forks_repo_name: exanauts/proxALM


text:
module TestInterface
using Compat.Test

# Checking interface through Diversity.jl
using Diversity
using SpatialEcology
using EcoBase

@testset "SpatialEcology-Diversity interface" begin
    numspecies = 10
    numcommunities = 8
    manyweights = rand(numspecies, numcommunities)
    manyweights /= sum(manyweights)

    species = map(n -> "Species $n", 1:numspecies)
    communities = map(n -> "SC $n", 1:numcommunities)
    #ut = UniqueTypes(species)
    #sc = Subcommunities(communities)
    mc = ComMatrix(manyweights, specnames = species, sitenames = communities)

    @test all(EcoBase.thingnames(mc) .== species)
    @test all(EcoBase.placenames(mc) .== communities)
    @test all(EcoBase.occurrences(mc) .≈ manyweights)
    @test all(EcoBase.richness(mc) .== repeat([numspecies], inner=numcommunities))
    @test all(EcoBase.occupancy(mc) .== repeat([numcommunities], inner=numspecies))
    fewerweights = deepcopy(manyweights)
    fewerweights[1, 1] = 0
    fewerweights /= sum(fewerweights)
    fmc = ComMatrix(fewerweights, specnames = species, sitenames = communities)

    @test EcoBase.noccupied(fmc) == numcommunities
    @test EcoBase.noccurring(fmc) == numspecies
    @test EcoBase.noccupied(fmc, 1) == numcommunities - 1
    @test EcoBase.noccurring(fmc, 1) == numspecies - 1
    @test EcoBase.nthings(fmc) == numspecies
    @test EcoBase.nplaces(fmc) == numcommunities
end

@testset "EcoBase-Diversity interface" begin
    numspecies = 10
    numcommunities = 8
    manyweights = rand(numspecies, numcommunities)
    manyweights /= sum(manyweights)

    species = map(n -> "Species $n", 1:numspecies)
    communities = map(n -> "SC $n", 1:numcommunities)
    ut = UniqueTypes(species)
    sc = Subcommunities(communities)
    mc = Metacommunity(manyweights, ut, sc)

    @test all(thingnames(mc) .== species)
    @test all(placenames(mc) .== communities)
    @test all(occurrences(mc) .≈ manyweights)
    @test all(richness(mc) .== repeat([numspecies], inner=numcommunities))
    @test all(occupancy(mc) .== repeat([numcommunities], inner=numspecies))
    fewerweights = deepcopy(manyweights)
    fewerweights[1, 1] = 0
    fewerweights /= sum(fewerweights)
    fmc = Metacommunity(fewerweights, ut, sc)

    @test noccupied(fmc) == numcommunities
    @test noccurring(fmc) == numspecies
    @test noccupied(fmc, 1) == numcommunities - 1
    @test noccurring(fmc, 1) == numspecies - 1
    @test nthings(fmc) == numspecies
    @test nplaces(fmc) == numcommunities
end

end


text:
mutable struct TreeModelXMLElement <: MyXMLElement
    el::XMLOrNothing
    node_traits::Vector{String}
    trait_dims::Vector{Int}
    param_names::Vector{String}
    newick_el::NewickXMLElement
    fixed_tree::Bool
    id::String

    function TreeModelXMLElement(newick_el::NewickXMLElement,
                dim::Int;
                trait::String = bn.DEFAULT_TRAIT_NAME)
        param_name = bn.LEAF_TRAITS
        return new(nothing, [trait], [dim], [param_name], newick_el, true,
                    bn.TREE_MODEL)
    end

    function TreeModelXMLElement(newick_el::NewickXMLElement,
                node_traits::Vector{String},
                trait_dims::Vector{Int},
                param_names::Vector{String}
                )
        return new(nothing, node_traits, trait_dims, param_names, newick_el, true)
    end

    function TreeModelXMLElement(newick_el::NewickXMLElement,
                                 node_traits::Vector{String},
                                 trait_dims::Vector{Int}
                                )
        param_names = ["$trait.leafTrait" for trait in node_traits]
        return TreeModelXMLElement(newick_el, node_traits, trait_dims, param_names)
    end
end

function name(::TreeModelXMLElement)
    return bn.TREE_MODEL
end

function get_id(tm::TreeModelXMLElement)
    return tm.id
end

function make_xml(tl::TreeModelXMLElement)
    make_xml(tl.newick_el)
    tl.el =  make_treeModel(tl.newick_el.el,
                    tl.node_traits,
                    tl.trait_dims,
                    tl.param_names,
                    id=tl.id)
end

function make_treeModel(newick_el::XMLElement,
        trait_names::Vector{String},
        trait_dims::Vector{Int},
        param_names::Vector{String};
        fixed_tree::Bool = false,
        id::String = bn.TREE_MODEL)

    @assert name(newick_el) == bn.NEWICK

    el = new_element(bn.TREE_MODEL)
    set_attribute(el, bn.ID, id)
    set_attribute(el, bn.FIX_HEIGHTS, bn.TRUE)

    if fixed_tree
        set_attribute(el, bn.FIXED_TREE, bn.TRUE)
    end

    tree_el = new_element(bn.TREE)
    newick_id = attribute(newick_el, bn.ID)
    set_attribute(tree_el, bn.IDREF, newick_id)
    add_child(el, tree_el)

    rh_el = new_element(bn.ROOT_HEIGHT)
    add_parameter_id(rh_el, "$(bn.TREE_MODEL).$(bn.ROOT_HEIGHT)")
    add_child(el, rh_el)

    add_nodeHeights(el, "$(bn.TREE_MODEL).$(bn.INTERNAL_NODE_HEIGHTS)",
            internal = true)
    add_nodeHeights(el, "$(bn.TREE_MODEL).$(bn.ALL_INTENAL_NODE_HEIGHTS)",
            internal = true, root = true)

    for i = 1:length(trait_names)
        nt_el = new_element(bn.NODE_TRAITS)
        attrs = [(bn.ROOT_NODE, bn.FALSE),
                (bn.INTERNAL_NODES, bn.FALSE),
                (bn.LEAF_NODES, bn.TRUE),
                (bn.AS_MATRIX, bn.TRUE),
                (bn.TRAIT_DIM, string(trait_dims[i])),
                (bn.NAME, trait_names[i])]
        set_attributes(nt_el, attrs)
        add_parameter_id(nt_el, param_names[i])
        add_child(el, nt_el)
    end
    return el
end

function add_nodeHeights(pel::XMLElement, param_id::String;
        internal::Bool = false, root::Bool = false)

    el = new_element(bn.NODE_HEIGHTS)
    if internal
        set_attribute(el, bn.INTERNAL_NODES, bn.TRUE)
    end
    if root
        set_attribute(el, bn.ROOT_NODE, bn.TRUE)
    end
    add_parameter_id(el, param_id)
    add_child(pel, el)
    return el
end

function add_leaf_param!(txl::TreeModelXMLElement, trait_name::String, trait_dim::Int)
    push!(txl.node_traits, trait_name)
    push!(txl.trait_dims, trait_dim)
    push!(txl.param_names, "$(trait_name).leaf")
end


text:
import ExperimentalAnalysis
using Base.Test
using DataFrames
using GLM
# write your own tests here

df = DataFrame(a = randn(10),b = randn(10),c = randn(10),d = randn(10))
df[2,2] = NA

@test ExperimentalAnalysis.perform_example_analysis() == nothing
ExperimentalAnalysis.scattermatrix(randn(10,3))
ExperimentalAnalysis.scattermatrix(df)
ExperimentalAnalysis.scattermatrix(df, reglines=true)
ExperimentalAnalysis.scattermatrix(df, f ~ a + b)
ExperimentalAnalysis.scattermatrix(df, c ~ a + b)
ExperimentalAnalysis.scattermatrix(randn(10,3),["x1", "x2", "x3"])

modela = lm(b ~ a + c, df)
modelb = lm(d ~ a + c, df)
ExperimentalAnalysis.scattermatrix([modela, modelb])
ExperimentalAnalysis.modelheatmap([modela, modelb])
ExperimentalAnalysis.modelheatmap(["a", "b"], modela, modelb)
ExperimentalAnalysis.modelheatmap(["a", "b"], [modela, modelb])


text:
module RunTests

using Test

@time @testset "Helpers" begin include("HelpersTests/runtests.jl") end

@time @testset "Inference" begin include("InferenceTests/runtests.jl") end

@time @testset "TensorValues" begin include("TensorValuesTests/runtests.jl") end

@time @testset "Arrays" begin include("ArraysTests/runtests.jl") end

@time @testset "Fields" begin include("FieldsTests/runtests.jl") end

@time @testset "Polynomials" begin include("PolynomialsTests/runtests.jl") end

@time @testset "ReferenceFEs" begin include("ReferenceFEsTests/runtests.jl") end

include("../bench/runbenchs.jl")

end # module


text:
using JuliaConSubmission
using Test

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

max_stars_repo_name: baggepinnen/Soss-probprog

max_issues_repo_name: baggepinnen/Soss-probprog

max_forks_repo_name: baggepinnen/Soss-probprog


text:
# New factor interface, something perhaps like this

export calcFactorResidualTemporary

getFactorOperationalMemoryType(dfg::SolverParams) = CommonConvWrapper
getFactorOperationalMemoryType(dfg::NoSolverParams) = CommonConvWrapper

# Helper function to construct CF from a CCW
CalcFactor(ccwl::CommonConvWrapper) = CalcFactor( ccwl.usrfnc!, _getFMdThread(ccwl), 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)

"""
    $SIGNATURES

Sample the factor stochastic model `N::Int` times and store the samples in the preallocated `ccw.measurement` container.

DevNotes
- Use in place operations where possible and remember `measurement` is a `::Tuple`.
- TODO only works on `.threadid()==1` at present, see #1094
- Also see, JuliaRobotics/RoME.jl#465
"""
sampleFactor(cf::CalcFactor{<:AbstractFactor}, N::Int=1  ) = [getSample(cf) for _=1:N]

function Base.show(io::IO, x::CalcFactor)
  println(io, )
  printstyled(io, " CalcFactor:\n", color=:blue)
  println(io, "  .factor: ", typeof(x.factor))
end

Base.show(io::IO, ::MIME"text/plain", x::CalcFactor) = show(io, x)

"""
    $SIGNATURES

Function to calculate measurement dimension from factor sampling.

Notes
- Will not work in all situations, but good enough so far.
  - # TODO standardize via domain or manifold definition...??
"""
function calcZDim(cf::CalcFactor{T}) where {T <: AbstractFactor}
  #
  try
    M = getManifold(T)
    return manifold_dimension(M)
  catch
    try
      M = getManifold(cf.factor)
      return manifold_dimension(M)
    catch
      @warn "no method getManifold(::$(string(T))), calcZDim will attempt legacy length(sample) method instead"
    end
  end

  # NOTE try to make sure we get matrix back (not a vector)
  smpls = sampleFactor(cf, 2)[1]
  return length(smpls)
end

calcZDim(ccw::CommonConvWrapper) = calcZDim(CalcFactor(ccw))

calcZDim(cf::CalcFactor{<:GenericMarginal}) = 0

calcZDim(cf::CalcFactor{<:ManifoldPrior}) = manifold_dimension(cf.factor.M)

"""
    $SIGNATURES

Helper function for evaluating factor residual functions, by adding necessary `CalcFactor` wrapper.

Notes
- Factor must already be in a factor graph to work
- Will not yet properly support all multihypo nuances, more a function for testing
- Useful for debugging a factor.

Example
```julia
fg = generateCanonicalFG_Kaess()

residual = calcFactorResidual(fg, :x1x2f1, [1.0], [0.0], [0.0])
```

Related

[`calcFactorResidualTemporary`](@ref), [`_evalFactorTemporary!`](@ref), [`evalFactor`](@ref), [`approxConv`](@ref)
"""
calcFactorResidual(dfgfct::DFGFactor, args...; ccw::CommonConvWrapper=IIF._getCCW(dfgfct)) = CalcFactor(ccw)(args...)
calcFactorResidual(dfg::AbstractDFG, fctsym::Symbol, args...) = calcFactorResidual(getFactor(dfg, fctsym), args...)

"""
    $SIGNATURES

Evaluate the residual function for a single sample.

Notes
- Binary factors only at this stage, and `multihypo` does not have to be considered in this test
- Assumes calculation is for a single particle, so `meas::Tuple{Z,other}` is only a single particles value.

Example
```julia
residual = calcFactorResidualTemporary(Pose2Pose2(...), (RoME.Pose2,RoME.Pose2), (z_i,), (x1, x2))
```

Related

[`calcFactorResidual`](@ref), [`CalcResidual`](@ref), [`_evalFactorTemporary!`](@ref), [`approxConv`](@ref), [`_buildGraphByFactorAndTypes!`](@ref)
"""
function calcFactorResidualTemporary( fct::AbstractRelative,
                                      varTypes::Tuple,
                                      measurement,
                                      pts::Tuple;
                                      tfg::AbstractDFG = initfg(),
                                      _blockRecursion::Bool=false )
  #

  # build a new temporary graph
  _, _dfgfct = _buildGraphByFactorAndTypes!(fct, varTypes, pts, dfg=tfg, _blockRecursion=_blockRecursion)

  # get a fresh measurement if needed
  _measurement = if measurement != [] #length(measurement) != 0
    measurement
  else
    # now use the CommonConvWrapper object in `_dfgfct`
    cfo = CalcFactor(_getCCW(_dfgfct))
    sampleFactor(cfo, 1)[1]
  end

  # assume a single sample point is being run
  return calcFactorResidual(_dfgfct, _measurement, pts...)
end

## =============================================================================================
## FactorOperationalMemory helper constructors
## =============================================================================================

function ConvPerThread( X::AbstractVector{P},
                        zDim::Int,
                        factormetadata::FactorMetadata;
                        particleidx::Int=1,
                        activehypo= 1:length(params),
                        perturb=zeros(zDim),
                        res=zeros(zDim),
                        thrid_ = 0  ) where P
  #
  return ConvPerThread{typeof(res), typeof(factormetadata), Any}( thrid_,
                        particleidx,
                        factormetadata,
                        Int[activehypo;],
                        perturb,
                        X,
                        res )
end

function CommonConvWrapper( fnc::T,
                            X::AbstractVector{P},
                            zDim::Int,
                            params::AbstractVector{<:AbstractVector{Q}},
                            factormetadata::FactorMetadata;
                            partial::Bool=false,
                            hypotheses::H=nothing,
                            certainhypo=nothing,
                            activehypo= 1:length(params),
                            nullhypo::Real=0,
                            varidx::Int=1,
                            measurement::AbstractVector=Vector(Vector{Float64}(),),
                            particleidx::Int=1,
                            xDim::Int=size(X,1),
                            partialDims::AbstractVector{<:Integer}=1:length(X),
                            perturb=zeros(zDim),
                            res::AbstractVector{<:Real}=zeros(zDim),
                            threadmodel::Type{<:_AbstractThreadModel}=MultiThreaded,
                            inflation::Real=3.0,
                            vartypes=typeof.(getVariableType.(factormetadata.fullvariables)),
                            gradients=nothing) where {T<:AbstractFactor,P,H,Q}
  #
  # FIXME, deprecate params::Vector throughout codebase
  tup = tuple(params...)
  nms = tuple(factormetadata.variablelist...)
  p_ntup = NamedTuple{nms,typeof(tup)}(tup)

  return  CommonConvWrapper(fnc,
                            xDim,
                            zDim,
                            partial,
                            hypotheses,
                            certainhypo,
                            Float64(nullhypo),
                            p_ntup, # params,
                            varidx,
                            measurement,
                            threadmodel,
                            (i->ConvPerThread(X, zDim,factormetadata, particleidx=particleidx,
                                              activehypo=activehypo,
                                              perturb=perturb, res=res )).(1:Threads.nthreads()),
                            inflation,
                            partialDims,
                            DataType[vartypes...],
                            gradients)
end

function _resizePointsVector!(vecP::AbstractVector{P}, mkd::ManifoldKernelDensity, N::Int) where P
  #
  pN = length(vecP)
  resize!(vecP, N)
  for j in pN:N
    smp = AMP.sample(mkd, 1)[1]
    # @show j, smp, typeof(smp), typeof(vecP[j])
    vecP[j] = smp[1]
  end

  vecP
end

"""
    $(SIGNATURES)

Prepare the particle arrays `ARR` to be used for approximate convolution.
This function ensures that ARR has te same dimensions among all the parameters.
Function returns with ARR[sfidx] pointing at newly allocated deepcopy of the
existing values in getVal(Xi[.label==solvefor]).

Notes
- Return values `sfidx` is the element in ARR where `Xi.label==solvefor` and
- `maxlen` is length of all (possibly resampled) `ARR` contained particles.
- `Xi` is order sensitive.
- for initialization, solveFor = Nothing.
- `P = getPointType(<:InferenceVariable)`

DevNotes
- FIXME ARR internally should become a NamedTuple
"""
function prepareparamsarray(Xi::Vector{<:DFGVariable},
                            solvefor::Union{Nothing, Symbol},
                            N::Int=0;
                            solveKey::Symbol=:default  ) where P
  #
  # FIXME ON FIRE, refactor to new NamedTuple instead
  varParamsAll = Vector{Vector{Any}}()

  LEN = Int[]
  maxlen = N # FIXME see #105
  count = 0
  sfidx = 0

  for xi in Xi
    vecP = getVal(xi, solveKey=solveKey)
    push!(varParamsAll, vecP)
    LEN = length.(varParamsAll)
    maxlen = maximum([N; LEN])
    count += 1
    if xi.label == solvefor
      sfidx = count #xi.index
    end
  end

  # resample variables with too few kernels (manifolds points)
  SAMP = LEN .< maxlen
  for i in 1:count
    if SAMP[i]
      Pr = getBelief(Xi[i], solveKey)
      _resizePointsVector!(varParamsAll[i], Pr, maxlen)
    end
  end

  # TODO --rather define reusable memory for the proposal
  # we are generating a proposal distribution, not direct replacement for existing memory and hence the deepcopy.
  if sfidx > 0
    varParamsAll[sfidx] = deepcopy(varParamsAll[sfidx])
  end

  # get solvefor manifolds
  # FIXME deprecate use of (:null,)
  mani = length(Xi)==0 || sfidx==0 ? (:null,) : getManifold(Xi[sfidx])

  # FIXME, forcing maxlen to N results in errors (see test/testVariousNSolveSize.jl) see #105
  # maxlen = N == 0 ? maxlen : N
  return varParamsAll, maxlen, sfidx, mani
end

"""
    $SIGNATURES
Internal method to set which dimensions should be used as the decision variables for later numerical optimization.
"""
function _setCCWDecisionDimsConv!(ccwl::Union{CommonConvWrapper{F},
                                              CommonConvWrapper{Mixture{N_,F,S,T}}} ) where {N_,F<:Union{AbstractManifoldMinimize, AbstractRelativeMinimize, AbstractRelativeRoots, AbstractPrior},S,T}
  #

  # NOTE should only be done in the constructor
  ccwl.partialDims = if ccwl.partial
    Int[ccwl.usrfnc!.partial...]
  else
    Int[1:ccwl.xDim...]
  end

  nothing
end

"""
    $(SIGNATURES)

Prepare a common functor computation object `prepareCommonConvWrapper{T}` containing
the user factor functor along with additional variables and information using during
approximate convolution computations.

DevNotes
- TODO consolidate with others, see https://github.com/JuliaRobotics/IncrementalInference.jl/projects/6
"""
function prepareCommonConvWrapper!( F_::Type{<:AbstractRelative},
                                    ccwl::CommonConvWrapper{F},
                                    Xi::AbstractVector{<:DFGVariable},
                                    solvefor::Symbol,
                                    N::Int;
                                    needFreshMeasurements::Bool=true,
                                    solveKey::Symbol=:default  ) where {F <: AbstractFactor}
  #

  # FIXME, order of fmd ccwl cf are a little weird and should be revised.
  pttypes = getVariableType.(Xi) .|> getPointType
  PointType = 0 < length(pttypes) ? pttypes[1] : Vector{Float64}

  #FIXME, see #1321
  # vecPtsArr = Vector{Vector{Any}}()

  #TODO some better consolidate is needed
  ccwl.vartypes = typeof.(getVariableType.(Xi))

  # FIXME maxlen should parrot N (barring multi-/nullhypo issues)
  vecPtsArr, maxlen, sfidx, mani = prepareparamsarray( Xi, solvefor, N, solveKey=solveKey)

  # FIXME ON FIRE, what happens if this is a partial dimension factor?  See #1246
  ccwl.xDim = getDimension(getVariableType(Xi[sfidx]))
  # ccwl.xDim = length(vecPtsArr[sfidx][1])
  # TODO should be selecting for the correct multihypothesis mode

    # setup the partial or complete decision variable dimensions for this ccwl object
    # NOTE perhaps deconv has changed the decision variable list, so placed here during consolidation phase
    # TODO, should this not be part of `prepareCommonConvWrapper` -- only here do we look for .partial
    _setCCWDecisionDimsConv!(ccwl)

  # SHOULD WE SLICE ARR DOWN BY PARTIAL DIMS HERE (OR LATER)?
  # FIXME refactor new type higher up.
  tup = tuple(vecPtsArr...)
  nms = tuple(getLabel.(Xi)...)
  ntp = NamedTuple{nms,typeof(tup)}(tup)
  ccwl.params = ntp

  # get factor metadata -- TODO, populate, also see #784
  fmd = FactorMetadata(Xi, getLabel.(Xi), ccwl.params, solvefor, nothing)

  # TODO consolidate with ccwl??
  # FIXME do not divert Mixture for sampling
  # cf = _buildCalcFactorMixture(ccwl, fmd, 1, ccwl.measurement, ccwl.params) # TODO perhaps 0 is safer
  cf = CalcFactor( ccwl.usrfnc!, fmd, 0, length(ccwl.measurement), ccwl.measurement, ccwl.params)

  #  get variable node data
  vnds = Xi

  # option to disable fresh samples
  if needFreshMeasurements
    # TODO refactor
    ccwl.measurement = sampleFactor(cf, maxlen)
    # sampleFactor!(ccwl, maxlen, fmd, vnds)
  end

  # cache the measurement dimension
  ccwl.zDim = calcZDim(CalcFactor(ccwl))
  # set the 'solvefor' variable index -- i.e. which connected variable of the factor is being computed in this convolution.
  ccwl.varidx = sfidx

  # set each CPT
  for thrid in 1:Threads.nthreads()
    cpt_ = ccwl.cpt[thrid]
    cpt_.X = ccwl.params[sfidx]

    # used in ccw functor for AbstractRelativeMinimize
    # TODO JT - Confirm it should be updated here. Testing in _prepCCW
    resize!(cpt_.res, ccwl.zDim)
    fill!(cpt_.res, 0.0)
  end

  # calculate new gradients perhaps
  # J = ccwl.gradients(measurement..., pts...)

  return sfidx, maxlen, mani
end

function prepareCommonConvWrapper!( ccwl::Union{CommonConvWrapper{F},
                                                CommonConvWrapper{Mixture{N_,F,S,T}}},
                                    Xi::AbstractVector{<:DFGVariable},
                                    solvefor::Symbol,
                                    N::Int;
                                    kw...  ) where {N_,F<:AbstractRelative,S,T}
  #
  prepareCommonConvWrapper!(F, ccwl, Xi, solvefor, N; kw...)
end

#


text:
#=
suitable univariate convex functions defined on ℝ₊₊ and associated oracles
- h_val evaluates (sum of) h
- h_conj_dom_pos is true if domain of conjugate of h is ℝ₊₊, else domain is ℝ
- h_conj evaluates the convex conjugate of (sum of) h
- h_der1/h_der2/h_der3 evaluate 1st/2nd/3rd derivatives of h pointwise
- get_initial_point gives initial point for u,v,λ_i

TODO derive central initial points for get_initial_point and use real type
=#

"""
$(TYPEDEF)

The negative logarithm function ``x \\to - \\log(x)``.
"""
struct NegLogSSF <: SepSpectralFun end

h_val(xs::Vector{T}, ::NegLogSSF) where {T <: Real} =
    -sum(log, xs)

h_conj_dom_pos(::NegLogSSF) = true

h_conj(xs::Vector{T}, ::NegLogSSF) where {T <: Real} =
    -length(xs) - sum(log, xs)

h_der1(ds::Vector{T}, xs::Vector{T}, ::NegLogSSF) where {T <: Real} =
    (@. ds = -inv(xs))
h_der2(ds::Vector{T}, xs::Vector{T}, ::NegLogSSF) where {T <: Real} =
    (@. ds = xs^-2)
h_der3(ds::Vector{T}, xs::Vector{T}, ::NegLogSSF) where {T <: Real} =
    (@. ds = -2 * xs^-3)

function get_initial_point(d::Int, ::NegLogSSF)
    # TODO initial central point
    return (1, 1, 1)
end

"""
$(TYPEDEF)

The negative entropy function ``x \\to x \\log(x)``.
"""
struct NegEntropySSF <: SepSpectralFun end

h_val(xs::Vector{T}, ::NegEntropySSF) where {T <: Real} =
    sum(x * log(x) for x in xs)

h_conj_dom_pos(::NegEntropySSF) = false

h_conj(xs::Vector{T}, ::NegEntropySSF) where {T <: Real} =
    sum(exp(-x - 1) for x in xs)

h_der1(ds::Vector{T}, xs::Vector{T}, ::NegEntropySSF) where {T <: Real} =
    (@. ds = 1 + log(xs))
h_der2(ds::Vector{T}, xs::Vector{T}, ::NegEntropySSF) where {T <: Real} =
    (@. ds = inv(xs))
h_der3(ds::Vector{T}, xs::Vector{T}, ::NegEntropySSF) where {T <: Real} =
    (@. ds = -xs^-2)

function get_initial_point(d::Int, ::NegEntropySSF)
    # TODO initial central point
    return (1, 1, 1)
end

"""
$(TYPEDEF)

The negative square root function ``x \\to -x^{1/2}``. Note this is a special
case of the negative power: `NegPower01SSF(0.5)`.
"""
struct NegSqrtSSF <: SepSpectralFun end

h_val(xs::Vector{T}, ::NegSqrtSSF) where {T <: Real} =
    -sum(sqrt, xs)

h_conj_dom_pos(::NegSqrtSSF) = true

h_conj(xs::Vector{T}, ::NegSqrtSSF) where {T <: Real} =
    sum(inv, xs) / 4

h_der1(ds::Vector{T}, xs::Vector{T}, ::NegSqrtSSF) where {T <: Real} =
    (@. ds = inv(-2 * sqrt(xs)))
h_der2(ds::Vector{T}, xs::Vector{T}, ::NegSqrtSSF) where {T <: Real} =
    (@. ds = inv(4 * xs * sqrt(xs)))
h_der3(ds::Vector{T}, xs::Vector{T}, ::NegSqrtSSF) where {T <: Real} =
    (@. ds = -3 / (8 * abs2(xs) * sqrt(xs)))

function get_initial_point(d::Int, ::NegSqrtSSF)
    # TODO initial central point
    return (0, 1, 1)
end

"""
$(TYPEDEF)

The negative power function ``x \\to -x^p`` parametrized by ``p \\in (0, 1)``
"""
struct NegPower01SSF <: SepSpectralFun
    p::Real
    NegPower01SSF(p::Real) = (@assert 0 < p < 1; new(p))
end

h_val(xs::Vector{T}, h::NegPower01SSF) where {T <: Real} =
    -sum(x^T(h.p) for x in xs)

h_conj_dom_pos(::NegPower01SSF) = true

# -(p - 1) * sum((x / p)^q for x in xs)
function h_conj(xs::Vector{T}, h::NegPower01SSF) where {T <: Real}
    p = T(h.p)
    q = p / (p - 1)
    return (1 - p) * p^-q * sum(x^q for x in xs)
end

function h_der1(ds::Vector{T}, xs::Vector{T}, h::NegPower01SSF) where {T <: Real}
    p = T(h.p)
    pm1 = p - 1
    @. ds = -p * xs^pm1
    return ds
end
function h_der2(ds::Vector{T}, xs::Vector{T}, h::NegPower01SSF) where {T <: Real}
    p = T(h.p)
    pm2 = p - 2
    coef = -p * (p - 1)
    @. ds = coef * xs^pm2
    return ds
end
function h_der3(ds::Vector{T}, xs::Vector{T}, h::NegPower01SSF) where {T <: Real}
    p = T(h.p)
    pm3 = p - 3
    coef = -p * (p - 1) * (p - 2)
    @. ds = coef * xs^pm3
    return ds
end

function get_initial_point(d::Int, h::NegPower01SSF)
    # TODO initial central point
    return (0, 1, 1)
end

"""
$(TYPEDEF)

The power function ``x \\to x^p`` parametrized by ``p \\in (1, 2]``. Note for
``p = 2``, it is more efficient to use [`EpiPerSquare`](@ref).
"""
struct Power12SSF <: SepSpectralFun
    p::Real
    Power12SSF(p::Real) = (@assert 1 < p <= 2; new(p))
end

h_val(xs::Vector{T}, h::Power12SSF) where {T <: Real} =
    sum(x^T(h.p) for x in xs)

h_conj_dom_pos(::Power12SSF) = false

# (p - 1) * sum((x >= 0 ? zero(x) : (abs(x) / p)^q) for x in xs)
function h_conj(xs::Vector{T}, h::Power12SSF) where {T <: Real}
    p = T(h.p)
    q = p / (p - 1)
    val = zero(T)
    for x in xs
        if x < 0
            val += abs(x)^q
        end
    end
    return (p - 1) * p^-q * val
end

function h_der1(ds::Vector{T}, xs::Vector{T}, h::Power12SSF) where {T <: Real}
    p = T(h.p)
    pm1 = p - 1
    @. ds = p * xs^pm1
    return ds
end
function h_der2(ds::Vector{T}, xs::Vector{T}, h::Power12SSF) where {T <: Real}
    p = T(h.p)
    pm2 = p - 2
    coef = p * (p - 1)
    @. ds = coef * xs^pm2
    return ds
end
function h_der3(ds::Vector{T}, xs::Vector{T}, h::Power12SSF) where {T <: Real}
    p = T(h.p)
    pm3 = p - 3
    coef = p * (p - 1) * (p - 2)
    @. ds = coef * xs^pm3
    return ds
end

function get_initial_point(d::Int, h::Power12SSF)
    # TODO initial central point
    return (2d, 1, 1)
end


text:
module IBFS
# Required packages

    using DataArrays
    using YAML
    include("DataTypes.jl")
    include("Grids.jl")
    include("Solvers.jl")
    include("Simulation.jl")

    export Simulation

end # module