text:
export BNFNode, Empty, Sequence, Alternatives,  NonTerminal,
    CharacterLiteral, StringLiteral
export Constructor, StringCollector
export BNFRef, recognize, logReductions, loggingReductions
export BNFGrammar, DerivationRule
export AllGrammars

using NahaJuliaLib

trace_recognize = false

"""
    BNFNode
Abstract supertype for all structs that we use to implement a grammar.
"""
abstract type BNFNode end

"""
    recognize(::BNFNode, input::String; index, finish)
Attempt to parse `input` as the specified `BNFNode`.
Return three values: whether the node matched the input,
the value represented by the matched input,
and the next index into input.

"""
@trace trace_recognize recognize(n::BNFNode, input::String; index=1, finish=length(input) + 1) =
    recognize(n, input, index, finish)

"""
    Empty()
Succeed while matching nothing.
"""
@bnfnode struct Empty <: BNFNode
end

@trace trace_recognize function recognize(n::Empty, input::String, index::Int, finish::Int)
    return true, nothing, index
end

"""
    Sequence(nodes...)
Successively match each of nodes in turn.
"""
@bnfnode struct Sequence <: BNFNode
    elements::Tuple{Vararg{<:BNFNode}}

    Sequence(elements...) = new(elements)
    Sequence(elements::Tuple) = new(elements)
end

@trace trace_recognize function recognize(n::Sequence, input::String, index::Int, finish::Int)
    collected = []
    in = index
    for n1 in n.elements
        matched, v, i = recognize(n1, input, in, finish)
        if !matched
            return false, nothing, index
        end
        in = i
        push!(collected, v)
    end
    return true, collected, in
end

"""
    Alternatives(nodes...)
Matches any one element of `nodes`.
"""
@bnfnode struct Alternatives <: BNFNode
    alternatives::Tuple{Vararg{<:BNFNode}}

    Alternatives(alternatives...) = new(alternatives)
    Alternatives(alternatives::Tuple) = new(alternatives)
end

@trace trace_recognize function recognize(n::Alternatives, input::String, index::Int, finish::Int)
    bestv= nothing
    # If one of the alternatives is Empty, we want our match to
    # succeed.
    besti = index - 1
    for n1 in n.alternatives
        matched, v, i = recognize(n1, input, index, finish)
        if matched && i > besti
            bestv = v
            besti = i
        end
    end
    if besti < index
        return false, nothing, index
    end
    return true, bestv, besti
end

"""
    CharacterLiteral(c)
Matches the single character `c`.
"""
@bnfnode struct CharacterLiteral <: BNFNode
    character::Char
end

@trace trace_recognize function recognize(n::CharacterLiteral, input::String, index::Int, finish::Int)
    if index > min(finish, length(input))
        return false, nothing, index
    end
    c = input[index]
    if c == n.character
        return true, c, index + 1
    end
    return false, nothing, index
end

"""
    StringLiteral(str)
Matches the string `str`.
"""
@bnfnode struct StringLiteral <: BNFNode
    str::AbstractString
end

@trace trace_recognize function recognize(n::StringLiteral, input::String, index::Int, finish::Int)
    e = index + length(n.str) - 1
    if e > min(finish, length(input))
        return false, nothing, index
    end
    ss = SubString(input, index, e)
    if n.str == ss
        return true, ss, e + 1
    end
    return false, nothing, index
end

"""
    Constructor(node, constructor_function)
Apply `constructor_function` to rhe result of recognizing `node`
and return that as the result.
"""
@bnfnode struct Constructor <: BNFNode
    node::BNFNode
    constructor
end

logReductions = false

function loggingReductions(f, log=true)
    global logReductions
    previous = logReductions
    try
        logReductions = true
        f()
    finally
        logReductions = previous
    end
end

@trace trace_recognize function recognize(n::Constructor, input::String, index::Int, finish::Int)
    matched, v, i = recognize(n.node, input, index, finish)
    if !matched
        return false, v, i
    end
    v2 = n.constructor(v)
    if logReductions
        @info "$(n.constructor) reduced $(typeof(v)) $v to $(typeof(v2)) $v2"
    end
    return true, v2, i
end

"""
    BNFGrammar
Represents a single grammar which can consist of a number of
`DerivationRule`s.
"""
struct BNFGrammar
    name::Symbol
    derivations # ::Dict{String, DerivationRule}

    function BNFGrammar(name::Symbol)
        g = new(name,
                # Dict{String, DerivationRule}()
                Dict())
        AllGrammars[g.name] = g
        g
    end
end

Base.haskey(grammar::BNFGrammar, rule::String) =
    haskey(grammar.derivations, rule)

Base.getindex(grammar::BNFGrammar, rule::String) =
    grammar.derivations[rule]

"""
A Dict of all defined grammars.
"""
AllGrammars = Dict{Symbol, BNFGrammar}()

function Base.getindex(grammar::BNFGrammar, nonterminal)
    grammar.derivations[nonterminal]
end

"""
    DerivationRule(grammar, rule_name, expression)
Implements a single production named `name` in the specified `grammar`.
One can include `expression` in other expressions using
`BNFRef(grammar, rule_name)`.
"""
@bnfnode struct DerivationRule <: BNFNode
    grammar_name::Symbol
    name::String
    lhs::BNFNode

    function DerivationRule(grammar::BNFGrammar, name, lhs)
        p = new(grammar.name, name, lhs)
        add_derivation(p)
        p
    end
    DerivationRule(grammar_name::Symbol, name, lhs) =
        DerivationRule(AllGrammars[grammar_name], name, lhs)
end

function Base.getproperty(p::DerivationRule, property::Symbol)
    if property in fieldnames(typeof(p))
        return getfield(p, property)
    end
    if property == :grammar
        return AllGrammars[p.grammar_name]
    end
end

function propertynames(p::DerivationRule)
    [:grammar, fieldnames(typeof(p))...]
end

@trace trace_recognize function recognize(n::DerivationRule, input::String, index::Int, finish::Int)
    recognize(n.lhs, input, index, finish)
end

function add_derivation(p::DerivationRule)
    add_derivation(p.grammar, p)
end

function add_derivation(grammar::BNFGrammar, derivation::DerivationRule)
    if haskey(grammar.derivations, derivation.name)
        throw(ErrorException(
            "The nonterminal $(derivation.name) is already defined in BNFGrammar $(grammar.name)."))
    end
    grammar.derivations[derivation.name] = derivation
end

# ????? Maybe we should merge DerivationRule and BNFRef, but it is
# really convenient for DerivationRules to automatically add
# themselves to their grammar.

"""
   BNFRef(grammar, name)
delegates to the "left hand side" of the `DerivationRule` named `name`
in `grammar`.
"""
@bnfnode struct BNFRef <:BNFNode
    grammar_name::Symbol
    name::String

    BNFRef(grammar_name::Symbol, name::String) =
        new(grammar_name, name)

    BNFRef(grammar::BNFGrammar, name::String) =
        BNFRef(grammar.name, name)
end

function lhs(n::BNFRef)
    AllGrammars[n.grammar_name][n.name].lhs
end

@trace trace_recognize function recognize(n::BNFRef, input::String, index::Int, finish::Int)
    recognize(lhs(n), input, index, finish)
end

"""
   StringCollector
StringCollector returns the entire substrring of the input that
@trace trace_recognize was recognized by its subexpression.
"""
@bnfnode struct StringCollector <: BNFNode
    node::BNFNode
end

"""
# Why is this 0?

SUBSTRING_SIZE = let
    alpha = *((('a':'z'))...)
    @allocated SubString(alpha, 3, 20)
end

# I was hoping to establish a threshhold to etermine if String or
# SubString was cheaper
"""

@trace trace_recognize function recognize(n::StringCollector, input::String, index::Int, finish::Int)
    start = index
    matched, v, i = recognize(n.node, input, index, finish)
    if !matched
        return false, v, i
    end
    return true, SubString(input, start, i - 1), i
end


text:
const ToT = UInt8
const ChannelID = UInt8
const DOMID = UInt32
const Floor = UInt8
const DU = UInt8
const HitTime = Float64
const TriggerMask = Int64

struct OMKey
    du::DU
    floor::Floor
end

struct Position{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T

    Position(x::T, y::T, z::T) where {T} = new{T}(x, y, z)
end
Position(x, y, z) = Position(promote(x, y, z)...)

struct Direction{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end
Direction(x, y, z) = Direction(promote(x, y, z)...)
Direction(ϕ, θ) = Direction(cos(ϕ)*cos(θ), sin(ϕ)*cos(θ), sin(θ))

struct MCEventInfo
    event_id::Int64
    group_id::Int64
    mc_id::Int64
    mc_time::Float64
    nanoseconds::Int64
    run_id::Int64
    timestamp::Int64
    weight_w1::Float64
    weight_w2::Float64
    weight_w3::Float64
    weight_w4::Float64
end

Base.show(io::IO, e::MCEventInfo) = begin
    print(io, "MCEventInfo: id($(e.event_id)), mc_id($(e.mc_id)), " *
          "mc_time($(e.mc_time) ($(e.nanoseconds)), ts($(e.timestamp)), run_id($(e.run_id))")
end

struct DAQEventInfo
    det_id::Int64
    event_id::Int64
    frame_index::Int64
    group_id::Int64
    mc_run_id::Int64
    mc_time::Float64
    nanoseconds::Int64
    overlays::Int64
    run_id::Int64
    timestamp::Int64
    trigger_counter::Int64
    trigger_mask::Int64
    weight_w1::Float64
    weight_w2::Float64
    weight_w3::Float64
    weight_w4::Float64
end

Base.show(io::IO, e::DAQEventInfo) = begin
    print(io, "DAQEventInfo: id($(e.event_id)), det_id($(e.det_id)), " *
          "timestamp($(e.timestamp) ($(e.nanoseconds)), run_id($(e.run_id))")
end

struct TimesliceInfo
    frame_index::UInt32
    slice_id::UInt32
    timestamp::UInt32
    nanoseconds::UInt32
    n_frames::UInt32
    group_id::UInt32
end

struct Track
    dir::Direction
    pos::Position
    time
end

# Fit
abstract type AbstractRecoTrack end

struct RecoTrack<:AbstractRecoTrack
    dir::Direction
    pos::Position
    time::Float64
end

struct NoRecoTrack<:AbstractRecoTrack end

# MC
struct MCTrack
    bjorken_y::Float64
    dir_x::Float64
    dir_y::Float64
    dir_z::Float64
    E::Float64
    group_id::Int64
    id::Int64
    interaction_channel::Int64
    is_cc::Bool
    length::Float64
    pos_x::Float64
    pos_y::Float64
    pos_z::Float64
    t::Float64
    particle_type::Int64
end

Track(t::MCTrack) = Track([t.dir_x, t.dir_y, t.dir_z], [t.pos_x, t.pos_y, t.pos_z], t.t)

Base.show(io::IO, t::MCTrack) = begin
    E = Printf.@sprintf "%0.1f" t.E
    bjorken_y = Printf.@sprintf "%0.2f" t.bjorken_y
    print(io, "MCTrack: bjorken_y($(bjorken_y)), t($(t.t)), " *
              "pos($(t.pos_x), $(t.pos_y), $(t.pos_z)), " *
              "dir($(t.dir_x), $(t.dir_y), $(t.dir_z)), " *
              "E($(E)), type($(t.particle_type))")
end

# Hardware
struct PMT
    channel_id::ChannelID
    pos::Position
    dir::Direction
end

struct DOM
    id::UInt32
    floor::Floor
    du::DU
    pmts::Vector{PMT}
end

struct Calibration
    det_id::Int32
    pos::Dict{Int32,Vector{NeRCA.Position}}
    dir::Dict{Int32,Vector{NeRCA.Direction}}
    t0::Dict{Int32,Vector{Float64}}
    du::Dict{Int32,DU}
    floor::Dict{Int32,Floor}
    omkeys::Dict{OMKey,Int32}  # maps OMKey to Module ID
    pmts::Dict{Int32,Tuple{Int32,UInt8}}
    max_z
    n_dus
end

Base.show(io::IO, c::Calibration) = begin
    print(io, "Calibration data for detector '$(c.det_id)' " *
              "with $(length(c.pos)) modules.")
end

# Signal
abstract type AbstractHit end
abstract type AbstractDAQHit<:AbstractHit end
abstract type AbstractMCHit<:AbstractHit end

Base.isless(lhs::AbstractHit, rhs::AbstractHit) = lhs.t < rhs.t

struct Hit <: AbstractDAQHit
    channel_id::ChannelID
    dom_id::DOMID
    t::HitTime
    tot::ToT
    trigger_mask::TriggerMask
end

struct SnapshotHit <: AbstractDAQHit
    channel_id::ChannelID
    dom_id::DOMID
    t::HitTime
    tot::ToT
end

struct MCHit <: AbstractHit
    a::Float32
    origin::UInt32
    pmt_id::UInt32
    t::HitTime
end

mutable struct Multiplicity
    count::Int32
    id::Int64
end

struct CalibratedHit <: AbstractDAQHit
    channel_id::ChannelID
    dom_id::UInt32
    du::DU
    floor::Floor
    t::HitTime
    tot::ToT
    pos::Position
    dir::Direction
    t0::HitTime
    trigger_mask::TriggerMask
    multiplicity::Multiplicity
end

struct TimesliceHit <: AbstractDAQHit
    channel_id::Int8
    dom_id::UInt32
    t::Int32
    tot::Int16
end

struct TriggeredHit <: AbstractDAQHit
    dom_id::Int32
    channel_id::UInt8
    t::Int32
    tot::UInt8
    trigger_mask::TriggerMask
end

struct DAQEvent
    det_id::Int32
    run_id::Int32
    timeslice_id::Int32
    timestamp::Int32
    ticks::Int32
    trigger_counter::Int64
    trigger_mask::Int64
    overlays::Int32
    n_triggered_hits::Int32
    triggered_hits::Vector{TriggeredHit}
    n_hits::Int32
    hits::Vector{Hit}
end

Base.show(io::IO, d::DAQEvent) = begin
    print(io, "DAQEvent: $(d.n_triggered_hits) triggered hits, " *
              "$(d.n_hits) snapshot hits")
end

ext: jl


text:
using DSP, JetPackDSP, Jets, Test

n1,n2,n3 = 500,5,7

@testset "Envelope^($(power)) calculation, dimension=$(dim), T=$(T)" for power in (0.25, 0.5, 1.0), T in (Float64,Float32), dim in (2,3)
    space = dim == 2 ? JetSpace(T, n3, n1) : JetSpace(T, n3, n2, n1)
    F = JopEnvelope(space, power)
    x1 = -1 .+ 2 .* rand(domain(F))
    y1 = F*x1
    y2 = (x1.^2 + imag.(hilbert(x1)).^2).^(power * 0.5)
    @test y1 ≈ y2
end

@testset "Envelope^($(power)) linearity test, dimension=$(length(N)), T=$(T)" for power in (0.25, 0.5, 1.0), T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    N=(n1,n2)
    T=Float64
    pwr = 0.25
    F = JopEnvelope(JetSpace(T,N), pwr)
    J  = jacobian!(F, -1 .+ 2 .* rand(domain(F)))

    #  check: J m1 + J m2 - J m3 = J (m1 + m2 - m3)
    m1 = -1 .+ 2 .* rand(domain(F))
    m2 = -1 .+ 2 .* rand(domain(F))
    m3 = -1 .+ 2 .* rand(domain(F))
    @test J*m1 .+ J*m2 .- J*m3 ≈ J*(m1 .+ m2 .- m3)

    #  check: J' d1 + J' d2 -J' d3 = J' (d1 + d2 - d3)
    d1 = -1 .+ 2 .* rand(range(F))
    d2 = -1 .+ 2 .* rand(range(F))
    d3 = -1 .+ 2 .* rand(range(F))
    @test J'*d1 .+ J'*d2 .- J'*d3 ≈ J'*(d1 .+ d2 .- d3)
end

@testset "Envelope^($(power)) dot product test, dimension=$(length(N)), T=$(T)" for power in (0.25, 0.5, 1.0), T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    F = JopEnvelope(JetSpace(T,N))
    J  = jacobian!(F, -1 .+ 2 .* rand(domain(F)))
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(F)), -1 .+ 2 .* rand(range(F)))
    @test lhs ≈ rhs
end

@testset "Envelope^($(power)) dot product test with zeros, dimension=$(length(N)), T=$(T)" for power in (0.25, 0.5, 1.0), T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    op = JopEnvelope(JetSpace(T,N),power,damping=eps(T))
    m0 = -1 .+ 2 .* rand(domain(op))
    m0[1:n1] .= T(0.0)
    J  = jacobian!(op, m0)
    lhs, rhs = dot_product_test(J, -1 .+ 2 .* rand(domain(op)), -1 .+ 2 .* rand(range(op)))
    @test isapprox(lhs, rhs)
end

@testset "Envelope^($(power)) linearization test, dimension=$(length(N)), T=$(T)" for power in (0.25, 0.5, 1.0), damping in (0.0, 1e-16, 1e-8, 1e-4, 1e-2, 1, 100), T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
    T = Float32
    N =  (n1,n2,n3)
    F = JopEnvelope(JetSpace(T,N), power, damping=T(damping))
    m0 = -1 .+ 2 .* rand(domain(F))

    mu = .1 * sqrt.([1.0,1.0/2.0,1.0/4.0,1.0/8.0,1.0/16.0,1.0/32.0,1.0/64.0,1.0/128.0,1.0/256.0,1.0/512.0])
    observed, expected = linearization_test(F, m0, μ=mu)

    δ = minimum(abs, observed - expected)
    @test δ < .3
end

#TODO: Epsilon has to be large enough (~1e-4) for linearization test to pass when we have a column of zeros.
@testset "Envelope^($(power)) linearization test with zeros, dimension=$(length(N)), T=$(T)" for power in (0.25, 0.5, 1.0), damping in (1e-4, 1e-2, 1, 100), T in (Float64,Float32), N in ( (n1,n2), (n1,n2,n3) )
        op = JopEnvelope(JetSpace(T,N),power,damping=T(damping))
        m0 = -1 .+ 2 .* rand(domain(op))
        m0[1:n1] .= T(0.0)
        mu = .1 * sqrt.([1.0,1.0/2.0,1.0/4.0,1.0/8.0,1.0/16.0,1.0/32.0,1.0/64.0,1.0/128.0,1.0/256.0,1.0/512.0])

        observed, expected = linearization_test(op, m0, μ=mu)
        δ = minimum(abs, observed - expected)
        @test δ < .2
end

nothing


text:
using Similitude, Test

@test 1000normal(molarmass(Metric)) ≈ normal(molarmass(Gauss))
@test normal(luminousefficacy(Metric)) ≈ 1e7*normal(luminousefficacy(Gauss))


text:
const gfzrnx = string(@__DIR__,"/","gfzrnx",
                    isapple ? "_osx" :
                        islinux ? "_lx" :
                            throw("no gfz win yet") )


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

"""
The functions in this module are made to be able to easily simulate stress
states produced by some of the most common test machines.

Take for example the function `uniaxial_increment!`. In a push-pull machine
(with a smooth specimen), we know that the stress state is uniaxial (in the
measuring volume). Given the strain increment in the direction where the stress
is nonzero, we find a strain increment that produces zero stress in the other
directions. Similarly for the other functions.
"""
module Increments

import LinearAlgebra: norm
import Tensors: tovoigt, fromvoigt

import ..AbstractMaterial, ..integrate_material!
import ..Utilities: Symm2

export find_dstrain!, general_increment!, stress_driven_general_increment!,
       uniaxial_increment!, biaxial_increment!, stress_driven_uniaxial_increment!

"""
    find_dstrain!(material::AbstractMaterial, dstrain::AbstractVector{<:Real},
                  dt::Real, update_dstrain!::Function;
                  max_iter::Integer=50, tol::Real=1e-9)

Find a compatible strain increment for `material`.

Chances are you'll only need to call this low-level function directly if you
want to implement new kinds of strain optimizers. See `general_increment!`
and `stress_driven_general_increment!` for usage examples.

This is the skeleton of the optimizer. The individual specific optimizer
functions (`update_dstrain!)` only need to define how to update `dstrain`.
The skeleton itself isn't a Newton-Raphson root finder. It just abstracts away
the iteration loop, convergence checking and data plumbing, so it can be used,
among other kinds, to conveniently implement Newton-Raphson root finders.

The `dstrain` supplied to this function is the initial guess for the
optimization. At each iteration, it must be updated by the user-defined
corrector `update_dstrain!`, whose call signature is expected to be:

    update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T})
        where V <: AbstractVector{T} where T <: Real
      -> err::Real

`dstrain` is the current value of the strain increment, in Voigt format.
Conversion to tensor format uses `offdiagscale=2.0`. The function must update
the Voigt format `dstrain` in-place.

`dstress = stress - stress0`, where `stress` is the stress state predicted by
integrating the material for one timestep of length `dt`, using the current
value of `dstrain` as a driving strain increment, and `stress0` is the stress
state stored in `materials.variables.stress`.

`jacobian` is ∂σij/∂εkl (`material.variables_new.jacobian`), as computed by the
material implementation. In many cases, the dstrain optimization can actually be
performed by a Newton-Raphson root finder, so we pass the jacobian to facilitate
writing the update formula for such a root finder.

The return value `err` must be an error measure (Real, >= 0).

The update is iterated at most `max_iter` times, until `err` falls below `tol`.

If `max_iter` is reached and the error measure is still `tol` or greater,
`ErrorException` is thrown.

To keep features orthogonal, the timestep is **not** committed automatically.
We call `integrate_material!`, but not `update_material!`. In other words,
we only update `material.variables_new`.
"""
function find_dstrain!(material::AbstractMaterial, dstrain::AbstractVector{<:Real},
                       dt::Real, update_dstrain!::Function;
                       max_iter::Integer=50, tol::Real=1e-9)
    stress0 = tovoigt(material.variables.stress)  # observed
    T = typeof(dstrain[1])
    # @debug "---START---"
    for i=1:max_iter
        # @debug "$i, $dstrain, $stress0, $(material.variables.stress)"
        material.ddrivers.time = dt
        material.ddrivers.strain = fromvoigt(Symm2{T}, dstrain; offdiagscale=2.0)
        integrate_material!(material)
        stress = tovoigt(material.variables_new.stress)  # predicted
        dstress = stress - stress0
        jacobian = tovoigt(material.variables_new.jacobian)
        e = update_dstrain!(dstrain, dstress, jacobian)
        if e < tol
            return nothing
        end
    end
    error("No convergence in strain increment")
end

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

"""
    general_increment!(material::AbstractMaterial, dstrain_knowns::AbstractVector{Union{T, Missing}},
                       dstrain::AbstractVector{Union{T, Missing}}=dstrain_knowns,
                       max_iter::Integer=50, norm_acc::T=1e-9) where T <: Real

Find a compatible strain increment for `material`.

The material state (`material.variables`) and any non-`missing` components of
the *strain* increment `dstrain_knowns` are taken as prescribed. Any `missing`
components will be solved for.

The type of the initial guess `dstrain` is `AbstractVector{Union{T, Missing}}`
only so we can make it default to `dstrain_knowns`, which has that type. Any
`missing` components in `dstrain` will be replaced by zeroes before we invoke
the solver.

This routine computes the other components of the strain increment such that the
predicted stress state matches the stored one.

See `find_dstrain!`.
"""
function general_increment!(material::AbstractMaterial,
                            dstrain_knowns::AbstractVector{<:Union{T, Missing}},
                            dt::Real,
                            dstrain::AbstractVector{<:Union{T, Missing}}=dstrain_knowns,
                            max_iter::Integer=50, norm_acc::T=1e-9) where T <: Real
    function validate_size(name::String, v::AbstractVector)
        if ndims(v) != 1 || size(v)[1] != 6
            error("""Expected a length-6 vector for $(name), got $(typeof(v)) with size $(join(size(v), "×"))""")
        end
    end
    validate_size("dstrain_knowns", dstrain_knowns)
    validate_size("dstrain", dstrain)
    dstrain_actual::AbstractVector{T} = T[((x !== missing) ? x : T(0)) for x in dstrain]
    unknowns = Integer[k for k in 1:6 if dstrain_knowns[k] === missing]
    function update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) where V <: AbstractVector{T} where T <: Real
        dstrain_correction = -jacobian[unknowns, unknowns] \ dstress[unknowns]
        dstrain[unknowns] .+= dstrain_correction
        return norm(dstrain_correction)
    end
    find_dstrain!(material, dstrain_actual, dt, update_dstrain!, max_iter=max_iter, tol=norm_acc)
    dstrain[:] = dstrain_actual
    return nothing
end

"""
    stress_driven_general_increment!(material::AbstractMaterial,
                                     dstress_knowns::AbstractVector{<:Union{T, Missing}},
                                     dt::Real,
                                     dstrain::AbstractVector{T},
                                     max_iter::Integer=50, norm_acc::T=1e-9) where T <: Real

Find a compatible strain increment for `material`.

The material state (`material.variables`) and any non-`missing` components of
the *stress* increment `dstress_knowns` are taken as prescribed. Any `missing`
components will be solved for.

This routine computes a strain increment such that the predicted stress state
matches the stored one.

See `find_dstrain!`.
"""
function stress_driven_general_increment!(material::AbstractMaterial,
                                          dstress_knowns::AbstractVector{<:Union{T, Missing}},
                                          dt::Real,
                                          dstrain::AbstractVector{T},
                                          max_iter::Integer=50, norm_acc::T=1e-9) where T <: Real
    function validate_size(name::String, v::AbstractVector)
        if ndims(v) != 1 || size(v)[1] != 6
            error("""Expected a length-6 vector for $(name), got $(typeof(v)) with size $(join(size(v), "×"))""")
        end
    end
    validate_size("dstress_knowns", dstress_knowns)
    validate_size("dstrain", dstrain)
    dstrain_actual::AbstractVector{T} = T[((x !== missing) ? x : T(0)) for x in dstrain]
    knowns = Integer[k for k in 1:6 if dstress_knowns[k] !== missing]
    function update_dstrain!(dstrain::V, dstress::V, jacobian::AbstractArray{T}) where V <: AbstractVector{T} where T <: Real
        # Mutation of `dstress` doesn't matter, since `dstress` is freshly generated at each iteration.
        dstress[knowns] -= dstress_knowns[knowns]
        dstrain_correction = -jacobian \ dstress
        dstrain .+= dstrain_correction
        return norm(dstrain_correction)
    end
    find_dstrain!(material, dstrain_actual, dt, update_dstrain!, max_iter=max_iter, tol=norm_acc)
    dstrain[:] = dstrain_actual
    return nothing
end

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

"""
    uniaxial_increment!(material::AbstractMaterial, dstrain11::Real, dt::Real;
                        dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0],
                        max_iter::Integer=50, norm_acc::Real=1e-9)

Find a compatible strain increment for `material`.

The material state (`material.variables`) and the component 11 of the *strain*
increment are taken as prescribed. This routine computes the other components of
the strain increment such that the predicted stress state matches the stored
one.

Convenience function; see `general_increment!`.

See `find_dstrain!`.
"""
function uniaxial_increment!(material::AbstractMaterial, dstrain11::Real, dt::Real;
                             dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0],
                             max_iter::Integer=50, norm_acc::Real=1e-9)
    dstrain_knowns = [dstrain11, missing, missing, missing, missing, missing]
    general_increment!(material, dstrain_knowns, dt, dstrain, max_iter, norm_acc)
    return nothing
end

"""
    biaxial_increment!(material::AbstractMaterial, dstrain11::Real, dstrain12::Real, dt::Real;
                       dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0, 0, dstrain12],
                       max_iter::Integer=50, norm_acc::Real=1e-9)

Find a compatible strain increment for `material`.

By "biaxial", we mean a stress state with one normal component and one shear
component.

The material state (`material.variables`) and the components 11 and 12 of the
*strain* increment are taken as prescribed. This routine computes the other
components of the strain increment such that the predicted stress state matches
the stored one.

Convenience function; see `general_increment!`.

See `find_dstrain!`.
"""
function biaxial_increment!(material::AbstractMaterial, dstrain11::Real, dstrain12::Real, dt::Real;
                            dstrain::AbstractVector{<:Real}=[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0, 0, dstrain12],
                            max_iter::Integer=50, norm_acc::Real=1e-9)
    dstrain_knowns = [dstrain11, missing, missing, missing, missing, dstrain12]
    general_increment!(material, dstrain_knowns, dt, dstrain, max_iter, norm_acc)
    return nothing
end

"""
    stress_driven_uniaxial_increment!(material::AbstractMaterial, dstress11::Real, dt::Real;
                                      dstrain::AbstractVector{<:Real}=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0],
                                      max_iter::Integer=50, norm_acc::Real=1e-9)

Find a compatible strain increment for `material`.

The material state (`material.variables`) and the component 11 of the *stress*
increment are taken as prescribed. This routine computes a strain increment such
that the predicted stress state matches the stored one.

Convenience function; see `stress_driven_general_increment!`.

See `find_dstrain!`.
"""
function stress_driven_uniaxial_increment!(material::AbstractMaterial, dstress11::Real, dt::Real;
                                           dstrain::AbstractVector{<:Real}=[dstress11/200e3, -0.3*dstress11/200e3, -0.3*dstress11/200e3, 0.0, 0.0, 0.0],
                                           max_iter::Integer=50, norm_acc::Real=1e-9)
    dstress_knowns = [dstress11, missing, missing, missing, missing, missing]
    stress_driven_general_increment!(material, dstress_knowns, dt, dstrain, max_iter, norm_acc)
    return nothing
end

end


text:
using StatsBase, PLINK

# generate random binary genotype data
n = 5000
p = 24000
srand(2016)
x_float = sample([0.0, 1.0, 2.0], (n,p))
y = PLINK.compress(x_float)
yt = PLINK.compress(x_float')
yt[3] = PLINK.ZERO8 # indicate row-major ordering

# save data to file
write(open("x_test.bed", "w"), y)
write(open("xt_test.bed", "w"), yt)

# need to make BIM, FAM files
# make BIM file first
BIM = open("x_test.bim", "w")
#for i = 1:(p-1)
for i = 1:p
    @printf(BIM, "%d\t%s\t%d\t%d\t%s\t%s\n", 1, "rs" * "$i", 0, i, "A", "T")
end
#@printf(BIM, "%d\t%s\t%d\t%d\t%s\t%s", 1, "rs" * "$p", 0, p, "A", "T")
close(BIM)

# now make FAM
FAM = open("x_test.fam", "w")
#for i = 1:(n-1)
for i = 1:n
    @printf(FAM, "%d\t%d\t%d\t%d\t%d\t%d\n", i, i, 0, 0, 1, 1)
end
#@printf(FAM, "%d\t%d\t%d\t%d\t%d\t%d", n, n, 0, 0, 1, 1)
close(FAM)

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
# See src/jd_vars : def_GLOBAL_LXDEFS

"""
$SIGNATURES

Internal function to resolve a `\\output{rpath}` (finds the output and shows it).
The `reprocess` argument indicates whether the output should be processed as judoc markdown or
inserted as is.
"""
function resolve_lx_output(lxc::LxCom, reprocess::Bool=false)::String
    rpath = strip(content(lxc.braces[1])) # [assets]/subpath/script{.jl}
    return resolve_lx_input_plainoutput(rpath, reprocess)
end

"""
$SIGNATURES

Internal function to resolve a `\\textoutput{rpath}` (finds the output and shows it as markdown).
Note: the output is re-processed (if there's any markdown in the output, it will be parsed).
"""
resolve_lx_textoutput(lxc::LxCom) = resolve_lx_output(lxc, true)

"""
$SIGNATURES

Internal function to resolve a `\\textinput{rpath}` (reads the file and show it after processing).
"""
function resolve_lx_textinput(lxc::LxCom, reprocess::Bool=false)::String
    rpath = strip(content(lxc.braces[1]))
    return resolve_lx_input_textfile(rpath)
end

"""
$SIGNATURES

Internal function to resolve a `\\figalt{alt}{rpath}` (finds a plot and includes it with alt).
"""
function resolve_lx_figalt(lxc::LxCom)::String
    rpath = strip(content(lxc.braces[2]))
    alt = strip(content(lxc.braces[1]))
    path = resolve_assets_rpath(rpath; canonical=false)
    fdir, fext = splitext(path)
    # there are several cases
    # A. a path with no extension --> guess extension
    # B. a path with extension --> use that
    # then in both cases there can be a relative path set but the user may mean
    # that it's in the subfolder /output/ (if generated by code) so should look
    # both in the relpath and if not found and if /output/ not already the last subdir
    candext = ifelse(isempty(fext), (".png", ".jpeg", ".jpg", ".svg", ".gif"), (fext,))
    for ext ∈ candext
        candpath = fdir * ext
        syspath  = joinpath(PATHS[:folder], split(candpath, '/')...)
        isfile(syspath) && return html_img(candpath, alt)
    end
    # now try in the output dir just in case (provided we weren't already looking there)
    p1, p2 = splitdir(fdir)
    if splitdir(p1)[2] != "output"
        for ext ∈ candext
            candpath = joinpath(p1, "output", p2 * ext)
            syspath  = joinpath(PATHS[:folder], split(candpath, '/')...)
            isfile(syspath) && return html_img(candpath, alt)
        end
    end
    return html_err("image matching '$path' not found")
end

"""
$SIGNATURES

Internal function to resolve a `\\tableinput{header}{rpath}` (finds a csv and includes it with header).
"""
function resolve_lx_tableinput(lxc::LxCom)::String
    rpath = strip(content(lxc.braces[2]))
    header = strip(content(lxc.braces[1]))
    path = resolve_assets_rpath(rpath; canonical=false)
    fdir, fext = splitext(path)
    # copy-paste from resolve_lx_figalt()
    # A. a path with extension --> use that
    # there can be a relative path set but the user may mean
    # that it's in the subfolder /output/ (if generated by code) so should look
    # both in the relpath and if not found and if /output/ not already the last subdir
    syspath  = joinpath(PATHS[:folder], split(path, '/')...)
    if isfile(syspath)
        return csv2html(syspath, header)
    end
    # now try in the output dir just in case (provided we weren't already looking there)
    p1, p2 = splitdir(fdir)
    if splitdir(p1)[2] != "output"
        candpath = joinpath(p1, "output", p2 * fext)
        syspath  = joinpath(PATHS[:folder], split(candpath, '/')...)
        isfile(syspath) && return csv2html(syspath, header)
    end
    return html_err("table matching '$path' not found")
end

"""
$SIGNATURES

Internal function to process an array of strings to markdown table (in one single string).
If header is empty, the first row of the file will be used for header.
"""
function csv2html(path, header)::String
    csvcontent = readdlm(path, ',', String, header=false)
    nrows, ncols = size(csvcontent)
    io = IOBuffer()
    # writing the header
    if ! isempty(header)
        # header provided
        newheader = split(header, ",")
        hs = size(newheader,1)
        hs != ncols && return html_err("header size ($hs) and number of columns ($ncols) do not match")
        write(io, prod("| " * h * " " for h in newheader))
        rowrange = 1:nrows
    else
        # header from csv file
        write(io, prod("| " * csvcontent[1, i] * " " for i in 1:ncols))
        rowrange = 2:nrows
    end
    # writing end of header & header separator
    write(io, "|\n|", repeat( " ----- |", ncols), "\n")
    # writing content
    for i in rowrange
        for j in 1:ncols
            write(io, "| ", csvcontent[i,j], " ")
        end
        write(io, "|\n")
    end
    return md2html(String(take!(io)))
end

"""
$SIGNATURES

Internal function to resolve a `\\file{name}{rpath}` (finds a file and includes it with link name).
Note that while `\\figalt` is tolerant to extensions not being specified, this one is not.
"""
function resolve_lx_file(lxc::LxCom)::String
    rpath = strip(content(lxc.braces[2]))
    name = strip(content(lxc.braces[1]))
    path = resolve_assets_rpath(rpath; canonical=false)
    syspath = joinpath(PATHS[:folder], split(path, '/')...)
    isfile(syspath) && return html_ahref(path, name)
    return html_err("file matching '$path' not found")
end

"""
$SIGNATURES

Dictionary of functions to use for different default latex-like commands. The output of these
commands is inserted "as-is" (without re-processing) in the HTML.
"""
const LXCOM_SIMPLE = Dict{String, Function}(
    "\\output" => resolve_lx_output,  # include plain output generated by code
    "\\figalt" => resolve_lx_figalt,  # include a figure (may or may not have been generated)
    "\\file"   => resolve_lx_file,    # include a file
    "\\tableinput" => resolve_lx_tableinput,    # include table from a csv file
    )

"""
$SIGNATURES

Same as [`LXCOM_SIMPLE`](@ref) except the output is re-processed before being inserted in the HTML.
"""
const LXCOM_SIMPLE_REPROCESS = Dict{String, Function}(
    "\\textoutput" => resolve_lx_textoutput, # include output generated by code and reproc
    "\\textinput"  => resolve_lx_textinput,
    )


text:
push!(LOAD_PATH, "../src/")
using Documenter, Healpix

makedocs(
    modules = [Healpix],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename = "Healpix.jl",
    pages = [
        "Introduction"  => "index.md",
        "Working with resolutions" => "resolutions.md",
        "Pixel functions" => "pixelfunc.md",
        "Map functions" => "mapfunc.md",
        "Spherical harmonics" => "alm.md",
        "Visualization" => "visualization.md",
        "Miscellanea" => "misc.md",
    ])

deploydocs(
    repo = "github.com/ziotom78/Healpix.jl.git",
)


text:
type InputLayer
    inputActivation::Vol
    outputActivation::Vol

    function InputLayer(dimensions)
        new(zeros(dimensions), zeros(dimensions))
    end
end

function foward(self::InputLayer, inputActivation::Vol)
    self.outputActivation = copy(self.inputActivation)
    self.outputActivation.w = copy(self.inputActivation.w)
    self.outputActivation.dw = copy(self.inputActivation.dw)
end

function backward(self::InputLayer)
    self.inputActivation.dw = copy(self.inputActivation.dw)
end

gradientStates(self::InputLayer) = []


text:
"""
    required_halo_size(tendency_term)

Returns the required size of halos for a term appearing
in the tendency for a velocity field or tracer field.

Example
=======

```jldoctest
using Oceananigans.Advection: CenteredFourthOrder
using Oceananigans.Utils: required_halo_size

required_halo_size(CenteredFourthOrder())

# output
2
"""
required_halo_size(anything) = 1

function inflate_halo_size(Hx, Hy, Hz, tendency_terms...)
    for term in tendency_terms
        H = required_halo_size(term)

        Hx = max(Hx, H)
        Hy = max(Hy, H)
        Hz = max(Hz, H)
    end

    return Hx, Hy, Hz
end


text:
# Code used to generate gradient waveforms that are moment-compensated
# Sequence optimization for diffusion motion-compensation

using MRIsim, JuMP, Ipopt #, Gtk
using MRIsim: get_Bmatrix, get_SRmatrix, get_Mmatrix,
              dur, get_bvalue, write_diff_fwf, delay
## Parameters
dwell_time = 6.4e-6
Gmax =  21e-3 #21e-3 45e-3 T/m
Smax = 50 #52.5 mT/m/ms
plots = true
N1 = 250 #floor(Int64, τ * 1e3 * 15625 / 100) + 2 # Δt = 6.4e-6 #dwell-time
#Path were to write the waveforms
#path = "/media/ccp/Samsung_T5/"
path_file = "/home/ccp/Documents/MRIsim.jl/"
k = 2
sym = false
maxwell = true
moment_nulled, bvalue0 = false, 50
seq_name = maxwell ? "MX_MC$(k)_" : "MC$(k)_"
seq_name = sym ? seq_name*"sym_" : seq_name
seq_name = seq_name * "b" * string(bvalue0)
# Read files
f = readlines(path_file*"qte_timings.txt")
timings = tryparse.(Float64, split(f[2]," "))
δ1, rf180, δ2 = timings[timings .≠ nothing]*1e-3
if sym
    rf180 += δ1 - δ2
    δ1 = δ2
end
#Grads# - Pre-defined RF waveforms.
τ = (δ1 + rf180 + δ2) # τ/Nt = Δt => Nt = τ/Δt
N2 = floor(Int, N1 * δ2 / δ1)
DIF = Sequence([Grad(x -> 1, δ1 - dwell_time, N1) delay(rf180) Grad(x -> 1, δ2 - dwell_time, N2)])
idx180 = N1 + 1
_, N = size(DIF.GR)
#Opt matrices
B =  get_Bmatrix(DIF)  #B-value
SR = get_SRmatrix(DIF) #Slew-rate
M =  get_Mmatrix(DIF; order=0); #Moments
M0v, M1v, M2v = M[1,:], M[2,:], M[3,:]
t = DIF.GR.T
T = [sum(t[1:i])  for i=1:N] #Acumulated time
# Optimazation
Mm = M[2:k+1,:] ./ 10 .^(2:k+1)
println("###")
model = Model(); set_optimizer(model, Ipopt.Optimizer); set_silent(model)
@variable(model, -Gmax <= x[1:N] <= Gmax, start=1, start=Gmax); #max-grads

if moment_nulled
    @objective(model, Min, -x'*B*x); #b-value
    @constraint(model, moments, Mm*x .== 0); #moments
else
    @variable(model, t1)
    # @variable(model, t2)
    @objective(model, Min, t1);
    @constraint(model, bvalue, x'*B*x == bvalue0); #b-value
    @constraint(model, moments, [t1; Mm*x] in MOI.NormOneCone(k+1))
    # @constraint(model, [t2; 1e-7*SR*x] in MOI.NormOneCone(N+2))
end
@constraint(model, M0, M0v'*x == 0) #diffusion condition
@constraint(model, start_seq,  x[1] .== 0); #seq
@constraint(model, end_seq,  x[end] .== 0); #seq
@constraint(model, start_RF180,  x[idx180] .== 0); #rf
@constraint(model, end_RF180,  x[idx180 + 1] .== 0); #rf
@constraint(model, slewrate, -Smax .<= SR*x .<= Smax); #slew rate
if maxwell
    @constraint(model, concomitant, sum(x[1:idx180-1].^2) == sum(x[idx180+1:end].^2)); #concomitant
end
optimize!(model);
gx = value.(x) #retrieving solution
#Seq object
setproperty!.(DIF.GR, :A, [gx'; 0*gx'; 0*gx'])
bmax = gx'*B*gx
@info "λ0 = $(abs(round(M0v'*gx/Gmax,digits=1))), λ1 = $(abs(round(M1v'*gx/Gmax,digits=1))), λ2 = $(abs(round(M2v'*gx/Gmax,digits=1)))"
@info "b-value: $(round(bmax, digits=2)) s/mm2"
@info round(δ1*1e3,digits=4), round(rf180*1e3,digits=4), round(δ2*1e3,digits=4)
@info seq_name
if termination_status(model) == MOI.LOCALLY_SOLVED
    @info "Solved!"
else
    @info "NOT Solved :("
end
#TO SCANNER
write_diff_fwf(DIF,idx180,Gmax,floor(Int64,bmax); filename=path_file*"QTE_Waveforms/$(seq_name).txt", name=seq_name)
write_diff_fwf(DIF,idx180,Gmax,floor(Int64,bmax); filename=path_file*"QTE_Waveforms/qte_vectors_input.txt", name=seq_name)
# Plots
# if plots
MRIsim.plot_grads_moments(DIF,title="ODTI, b=$(round(get_bvalue(DIF)*1e-6, digits=2)) s/mm2, λ0 = $(abs(round(1e3*M0v'*gx,digits=1))), λ1 = $(abs(round(1e3*M1v'*gx,digits=1))), λ2 = $(abs(round(1e3*M2v'*gx,digits=1)))")
# end
# end


text:
"""
    Ant <: Environment

    four-legged insect-like robot, based on https://gym.openai.com/envs/Ant-v2/
"""
struct Ant end

function ant(;
    representation=:minimal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.5,
    spring=0.0,
    damper=0.0,
    parse_damper=true,
    seed=1,
    contact_feet=true,
    contact_body=true,
    limits::Bool=true,
    info=nothing,
    vis=Visualizer(),
    name=:robot,
    opts_step=SolverOptions(),
    opts_grad=SolverOptions(),
    T=Float64)

    mechanism = get_ant(;
        timestep,
        gravity,
        friction_coefficient,
        spring,
        damper,
        parse_damper,
        contact_feet,
        contact_body,
        limits)

    initialize_ant!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end

    nu = 8
    no = nx + length(mechanism.contacts)

    aspace = BoxSpace(nu,
        low=(-ones(nu)),
        high=(ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = [zeros(8, 6) I(nu)]
    control_scaling = Diagonal(timestep * 150.0 * ones(nu))

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Ant, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev,
        control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)
    return env
end

function Base.step(env::Environment{Ant}, x, u;
    gradients=false,
    attitude_decompress=false)

    mechanism = env.mechanism
    timestep= mechanism.timestep

    # x position (before)
    xposbefore = x[1]

    x0 = copy(x)
    env.input_previous .= u  # for rendering in Gym
	u_scaled = env.control_map * u

    z0 = env.representation == :minimal ? minimal_to_maximal(mechanism, x0) : x0
    z1 = step!(mechanism, z0, u_scaled; opts=env.opts_step)
    env.state .= env.representation == :minimal ? maximal_to_minimal(mechanism, z1) : z1

    # x position (after)
    xposafter = env.state[1]

    # forward reward
    forward_reward = 2.0 * (xposafter - xposbefore) / timestep

    # control cost
	ctrl_cost = (0.05 * u' * u)[1]

    # contact cost
    contact_cost = 0.0

    for contact in mechanism.contacts
        contact_cost += 0.5 * 1.0e-3 * max(-1.0, min(1.0, contact.impulses[2][1]))^2.0
    end

	# survive_reward = 1.0
    survive_reward = 0.05

    # total reward
    reward = forward_reward - ctrl_cost - contact_cost + survive_reward

    # done ?
    done = !(all(isfinite.(env.state)) && (env.state[3] >= 0.2) && (env.state[3] <= 1.0))

    # Gradients
    if gradients
        if env.representation == :minimal
            fx, fu = get_minimal_gradients!(env.mechanism, z0, u_scaled, opts=env.opts_grad)
        elseif env.representation == :maximal
            fx, fu = get_maximal_gradients!(env.mechanism, z0, u_scaled, opts=env.opts_grad)
            if attitude_decompress
                if attitude_decompress
                    A0 = attitude_jacobian(z0, length(env.mechanism.bodies))
                    A1 = attitude_jacobian(z1, length(env.mechanism.bodies))
                    fx = A1 * fx * A0'
                    fu = A1 * fu
                end
            end
        end
        env.dynamics_jacobian_state .= fx
        env.dynamics_jacobian_input .= fu * env.control_map
    end

    info = Dict()

    return get_observation(env), reward, done, info
end

# TODO add random noise
function Base.reset(env::Environment{Ant};
    x=nothing)

    initialize!(env.mechanism, type2symbol(Ant))

    if x != nothing
        env.state .= x
    else
        x = get_minimal_state(env.mechanism)
        if env.representation == :minimal
            set_maximal_state!(env.mechanism, minimal_to_maximal(env.mechanism, x))
            env.state .= x
        elseif env.representation == :maximal
            z = minimal_to_maximal(env.mechanism, x)
            set_maximal_state!(env.mechanism, z)
            env.state .= z
        end
        env.input_previous .= 0.0
    end

    return get_observation(env)
end

function get_observation(env::Environment{Ant,T}) where T
    contact_force = T[]
    for contact in env.mechanism.contacts
        push!(contact_force, max(-1.0, min(1.0, contact.impulses[2][1])))
    end
    return [env.state; contact_force]
end


text:
using IntervalArithmetic, IntervalRootFinding
using ForwardDiff
using Base.Test

const D = IntervalRootFinding.derivative

setprecision(Interval, Float64)
float_pi = @interval(pi)

setprecision(Interval, 10000)
big_pi = @interval(pi)
# Using precision "only" 256 leads to overestimation of the true roots for `cos`
# i.e the Newton method gives more accurate results!

half_pi = big_pi / 2
three_halves_pi = 3*big_pi/2

@testset "Testing newton1d" begin

    f(x) = exp(x^2) - cos(x)
    f′(x) = 2*x*exp(x^2) + sin(x)
    f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24
    f1′(x) = 4x^3 - 30x^2 + 70x - 50

    for autodiff in (false, true)
        if autodiff
            rts1 = newton1d(sin, -5..5)
            rts2 = newton1d(f, -∞..∞)
            rts3 = newton1d(f1, -10..10)

        else
            rts1 = newton1d(sin, cos, -5..5)
            rts2 = newton1d(f, f′, -∞..∞)
            rts3 = newton1d(f1, f1′, -10..10)
        end

        @test length(rts1) == 3
        L = [(-π.. -π), (0..0), (π..π)]
        for i = 1:length(rts1)
            @test L[i] == rts1[i].interval && :unique == rts1[i].status
        end

        @test length(rts2) == 1
        @test (0..0) == rts2[1].interval && :unique == rts2[1].status

        @test length(rts3) == 4
        L = [1, 2, 3, 4]
        for i = 1:length(rts3)
            @test L[i] in rts3[i].interval && :unique == rts3[i].status
        end
    end
end


text:
@testset "523.continuous-subarray-sum.jl" begin
    @test check_subarray_sum([23, 2, 4, 6, 7], 6) == true
    @test check_subarray_sum([23, 2, 6, 4, 7], 6) == true
    @test check_subarray_sum([23, 2, 6, 4, 7], 13) == false
end


text:
using OnlineAI
using MLPlots
plotly()
default(leg=false)

F = Float32
S = SkanSynapse{F}
params = SkanParams{F}(
    ramp_step_adjustment = 0,
    ramp_step_initfunc = () -> F(rand(1:2)),
    weight_init = 2
  )

txt = " ABC ACD ADA CDB"
letters = sort(unique(txt))
# letteridx = Dict([(l,i) for (i,l) in enumerate(letters)])
I = length(letters)

# init the neurons
numinputs = I
inputs = [InputNeuron() for i in 1:numinputs]
numoutputs = I
outputs = [SkanNeuron(S[],
                      false,
                      false,
                      F(0),
                      params.weight_init,
                      F(0),
                      F(0)
                      # F(20numinputs),
                      # F(100numinputs)
                     ) for i in 1:numoutputs]
neurons = vcat(inputs, outputs)

# init the synapses
synapses = SkanSynapse[]
for input in inputs, output in outputs
    synapse = SkanSynapse(input,
                          output,
                          params.weight_init,
                          zero(F),
                          zero(F),
                          params.ramp_step_initfunc(),
                          zero(Int8))
    push!(output.incoming_synapses, synapse)
    push!(synapses, synapse)
end

# do a sim
n = 100

synapse_fields = [:ramp,:ramp_step,:ramp_flag]
sdict = Dict(zip(synapse_fields, [zeros(n, length(synapses)) for i=1:length(synapse_fields)]))

ispikes = SpikeTrains(numinputs, c=:red, title="Input Spikes", ms=10)
ospikes = SpikeTrains(numoutputs, title="Output Spikes", ms=10)
opulses = SpikeTrains(numoutputs, title="Output Pulses", ms=10)

# neuron_fields = [:pulse, :spike, :potential, :threshold]
neuron_fields = [:potential, :threshold]
ndict = Dict(zip(neuron_fields, [zeros(n, numoutputs) for i=1:length(neuron_fields)]))
# opulses, spikes, potentials, thresholds = [zeros(n, numoutputs) for i=1:4]

# inputspikes = zeros(n, numinputs)

for t in 1:n

    # create random spiking activity
    # for input in inputs
    #     input.spike = rand() < 0.01
    # end

    # only the input neuron corresponding to the current letter spikes
    l = txt[mod1(round(Int,t), length(txt))]
    for (i,input) in enumerate(inputs)
      input.spike = l == letters[i]
    end

    # update the network
    step!(params, neurons, synapses)

    # store the network state for plotting
    for (j,s) in enumerate(synapses), (k,v) in sdict
        v[t,j] = getfield(s, k)
    end
    for (j,s) in enumerate(outputs), (k,v) in ndict
        v[t,j] = getfield(s, k)
    end
    for (j,s) in enumerate(outputs)
        s.pulse && push!(opulses, j, t)
        s.spike && push!(ospikes, j, t)
    end
    for (j,s) in enumerate(inputs)
        # inputspikes[t,j] = s.spike
        s.spike && push!(ispikes, j, t)
    end

end

splts = [plot(sdict[k], title=string(k)) for k in synapse_fields]
# nplts = [plot(ndict[k], title=k) for k in neuron_fields]
nplt = plot(ndict[:potential], title="Membrane Potential")
plot!(ndict[:threshold])
# inplt = plot(inputspikes,title="Inputs Spikes")

subplot(splts..., nplt, ispikes.plt, ospikes.plt, opulses.plt, size=(1000,1200), nc=1, link=true)


text: const ASSET_FINGERPRINT = ""


text:
using GeoData, ArchGDAL, NCDatasets, Plots, Statistics

geturl(url, filename=splitdir(url)[2]) = begin
    isfile(filename) || download(url, filename)
    filename
end

# Load some layers from NetCDF #############################################

ncurl = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc"
ncfilename = geturl(ncurl, "tos_O1_2001-2002.nc")
stack = NCDstack(ncfilename)

# Load the sea surface temperature layer
A = stack[:tos]

# Plot the 1st, 4th, 7th and 10th months
A[Ti(1:3:12)] |> plot
savefig("tos_4.png")

# Plot the Australia region
A[Ti(Contains(DateTime360Day(2001, 01, 17))), Lat(Between(0, -50)), Lon(Between(100, 160))] |> plot
savefig("tos_jan_australia.png")

# Plot the mean
mean(A; dims=Ti) |> plot
savefig("mean_tos.png")

# Save the mean
write("mean_tos.ncd", NCDarray, mean(A; dims=Ti))

# Plo a transect
A[Lat(Contains(20)), Ti(1)] |> plot
savefig("tos_20deg_lat.png")

# Create a GeoSeries from multiple files.
# This uses the same data three times to avoid downloads, you
# would really use a series of datasets matching the time dimension dates.
timedim = Ti([DateTime360Day(2001, 01, 1), DateTime360Day(2001, 02, 1), DateTime360Day(2001, 03, 1)])
filenames = [ncfilename, ncfilename, ncfilename]
series = GeoSeries(filenames, (timedim,); childtype=NCDstack)

# Get a single array from the series
A = series[Near<|DateTime360Day(2001, 01, 1)]["tos"][Lon<|Between(50, 200)]

# Plot single slices
view(A, Ti(4)) |> plot

# Or reduce over some dimension
mean(A; dims=Ti) |> plot
std(A; dims=Ti) |> plot

# Replace missing values and to reduce over a dimension
minimum(replace_missing(A, NaN); dims=Ti) |> plot
maximum(replace_missing(A, NaN); dims=Ti) |> plot
reduce(+, A; dims=Ti) |> plot

# Plot the mean sea surface temperature for australia in the second half of 2002
stack = NCDstack(ncfilename)
t = Ti<|Between(DateTime360Day(2002, 07, 1), DateTime360Day(2002, 012, 30))
stack["tos"][t, Lat<|Between(-45, 0.5), Lon<|Between(110, 160)] |> x->mean(x; dims=Ti) |> plot

# Permute the dimensions in the underlying data
# It stil plots the right way up
permutedims(A, (Lat, Lon, Ti))[Ti(1:3:12)] |> plot

# Line plots have (kind of) useful labels
A = series[1]["tos"]
replace(A[Lat(80)], missing=>NaN) |> plot
A[Lon(170), Ti(10)] |> plot
A[Lat(1:80), Lon(170), Ti(10)] |> plot

# Load a tif with GDAL ######################################################
gdal_url = "https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif"
tif_filename = geturl(gdal_url)

array = GDALarray(tif_filename)
array[Band(1)] |> plot

# Make a stack (just using the same file three times for demonstration)
filepaths = [tif_filename, tif_filename, tif_filename]
diskstack = GDALstack(filepaths, (:one, :two, :three))
memstack = GeoStack(diskstack)
# Plot a section of the tif file
diskstack[:two][Band(1), Lat(1:100), Lon(1:150)] |> plot

# See the test suit for more examples...


text:
using Test
using CompScienceMeshes

mesh = meshrectangle(1.0, 1.0, 1.0)
faces = skeleton(mesh, 2)
verts = vertices(mesh, faces.faces[1])
p = simplex(verts)

@test p.vertices == [
    point(0.0, 0.0, 0.0),
    point(0.0, 1.0, 0.0),
    point(1.0, 0.0, 0.0)]
@test p.tangents == [
    point(-1.0, 0.0, 0.0),
    point(-1.0, 1.0, 0.0)]
@test p.normals == [
    point(0.0, 0.0, -1.0)]
@test p.volume == 0.5


text:
isbanded(A::AbstractTriangular) = isbanded(parent(A))
bandwidths(A::Union{UpperTriangular,UnitUpperTriangular}) =
    (min(0,bandwidth(parent(A),1)), bandwidth(parent(A),2))
bandwidths(A::Union{LowerTriangular,UnitLowerTriangular}) =
    (bandwidth(parent(A),1), min(0,bandwidth(parent(A),2)))

triangularlayout(::Type{Tri}, ML::BandedColumnMajor) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::BandedRowMajor) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::ConjLayout{BandedRowMajor}) where {Tri} = Tri(ML)

function tribandeddata(::TriangularLayout{'U'}, A)
    B = triangulardata(A)
    u = bandwidth(B,2)
    D = bandeddata(B)
    view(D, 1:u+1, :)
end

function tribandeddata(::TriangularLayout{'L'}, A)
    B = triangulardata(A)
    l,u = bandwidths(B)
    D = bandeddata(B)
    view(D, u+1:l+u+1, :)
end

tribandeddata(A) = tribandeddata(MemoryLayout(A), A)

Base.replace_in_print_matrix(A::Union{UpperTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitUpperTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

Base.replace_in_print_matrix(A::Union{LowerTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitLowerTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

# Mul
@lazylmul UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatMulVec{T, <:TriangularLayout{'U',UNIT,BandedColumnMajor},
                                   <:AbstractStridedLayout}) where {UNIT,T <: BlasFloat}
    A,x = M.A, M.B
    x ≡ dest || copyto!(dest, x)
    tbmv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatMulVec{T, <:TriangularLayout{'L',UNIT,BandedColumnMajor},
                                   <:AbstractStridedLayout}) where {UNIT,T <: BlasFloat}
    A,x = M.A, M.B
    x ≡ dest || copyto!(dest, x)
    tbmv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatMulVec{T, <:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.A, M.B
    x ≡ dest || copyto!(dest, x)
    tbmv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatMulVec{T, <:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.A, M.B
    x ≡ dest || copyto!(dest, x)
    tbmv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end

# Ldiv
@lazyldiv UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{T, <:TriangularLayout{'U',UNIT,BandedColumnMajor},
                                   <:AbstractStridedLayout}) where {UNIT,T <: BlasFloat}
    A,x = inv(M.A), M.B
    x ≡ dest || copyto!(dest, x)
    tbsv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{T, <:TriangularLayout{'L',UNIT,BandedColumnMajor},
                                   <:AbstractStridedLayout}) where {UNIT,T <: BlasFloat}
    A,x = inv(M.A), M.B
    x ≡ dest || copyto!(dest, x)
    tbsv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{T, <:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = inv(M.A), M.B
    x ≡ dest || copyto!(dest, x)
    tbsv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{T, <:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                   <:AbstractStridedLayout}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = inv(M.A), M.B
    x ≡ dest || copyto!(dest, x)
    tbsv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end


text:
using Missings
using SimpleTraits
using Unitful

newempty(::Type{Data}) where Data = Data()

struct LinkBranch{RT, NL, Data, LenUnits} <: AbstractBranch{RT, NL}
    name::Int
    inout::Tuple{AbstractNode{RT, NL}, AbstractNode{RT, NL}}
    length::Union{Missing, LenUnits}
    data::Data
end

mutable struct LinkNode{RT, NL, Data,
                        B <: AbstractBranch{RT, NL}} <: AbstractNode{RT, NL}
    name::NL
    inbound::Union{B, Missing}
    other::Vector{B}
    data::Data
    LinkNode{RT, NL, Data, B}(name::NL, data::Data = newempty(Data)) where
    {RT, NL, Data, B} = new{RT, NL, Data, B}(name, missing, Vector{B}(), data)

end

import Phylo.API: _prefernodeobjects
_prefernodeobjects(::Type{<:LinkNode}) = true

mutable struct LinkTree{RT, NL, N <: LinkNode{RT, NL},
                        B <: LinkBranch{RT, NL}, TD} <:
                            AbstractTree{OneTree, RT, NL, N, B}
    name::Union{String, Missing}
    nodedict::Dict{NL, Int}
    roots::Vector{N}
    nodes::Vector{Union{N, Missing}}
    branches::Vector{Union{B, Missing}}
    data::Dict{String, Any}
    tipdata::TD
    rootheight
    isvalid::Union{Bool, Missing}
    cache::Dict{TraversalOrder, Vector{N}}

    function LinkTree{RT, NL, N,
                      B, TD}(tipnames::Vector{NL} = NL[];
                             treename::Union{String, Missing} = missing,
                             tipdata::TD = newempty(TD),
                             rootheight = missing) where {RT, NL, N, B, TD}
        tree = new{RT, NL, N, B, TD}(treename, Dict{NL, N}(), Vector{N}(),
                                     Vector{Union{N, Missing}}(), Vector{B}(),
                                     Dict{String, Any}(),
                                     tipdata, rootheight, missing,
                                     Dict{TraversalOrder, Vector{N}}())
        if !isempty(tipnames)
            createnodes!(tree, tipnames)
        elseif !isempty(tipdata)
            tips = unique(keys(tipdata))
            createnodes!(tree, tips)
        end
        return tree
    end
end
function LinkTree{RT, NL, N, B, TD}(leafinfos::TD) where {RT, NL, N, B, TD}
    leafnames = unique(info[1] for info in getiterator(leafinfos))
    return LinkTree{RT, NL, N, B, TD}(leafnames; tipdata = leafinfos)
end

import Phylo.API: _branchdims
_branchdims(T::Type{<:LinkTree}) = _branchdims(branchtype(T))

const LB{RT, LenUnits} = LinkBranch{RT, String, Dict{String, Any}, LenUnits}
const LN{RT, LenUnits} = LinkNode{RT, String, Dict{String, Any}, LB{RT, LenUnits}}
const LT{RT, TD, LenUnits} = LinkTree{RT, String, LN{RT, LenUnits}, LB{RT, LenUnits}, TD}
const LTD{RT, LenUnits} = LT{RT, Dict{String, Any}, LenUnits}
const RootedTree = LTD{OneRoot, Float64}
const ManyRootTree = LTD{ManyRoots, Float64}
const UnrootedTree = LTD{Unrooted, Float64}

# LinkBranch methods
function LinkBranch(name::Int,
                    from::LinkNode{RT, NL},
                    to::LinkNode{RT, NL},
                    len = missing,
                    data::Data = nothing) where {RT, NL, Data}
    ismissing(len) || ustrip(len) >= 0.0 ||
        error("Branch length must be positive or missing (no recorded length), not $len")
    return ismissing(len) ?
        LinkBranch{RT, NL, Data, Float64}(name, (from, to), len, data) :
        LinkBranch{RT, NL, Data, typeof(len)}(name, (from, to), len, data)
end

import Phylo.API: _preferbranchobjects
_preferbranchobjects(::Type{<:LinkBranch}) = true
_branchdims(::Type{B}) where {RT, NL, Data, LenUnits,
                               B <: LinkBranch{RT, NL, Data, LenUnits}} =
    dimension(LenUnits)

import Phylo.API: _src
_src(::AbstractTree, branch::LinkBranch{<:Rooted}) = branch.inout[1]
import Phylo.API: _dst
_dst(::AbstractTree, branch::LinkBranch{<:Rooted}) = branch.inout[2]
import Phylo.API: _conn
function _conn(::AbstractTree, branch::LinkBranch{RT, NL, D},
               exclude::AbstractNode{RT, NL}) where {RT, NL, D}
    return exclude ≡ branch.inout[1] ? branch.inout[2] :
        (exclude ≡ branch.inout[2] ? branch.inout[1] :
         error("Branch $(branch.name) not connected to $(exclude.name)"))
end
import Phylo.API: _getlength
_getlength(::AbstractTree, branch::LinkBranch) = branch.length
import Phylo.API: _getbranchname
_getbranchname(tree::LinkTree, id::Int) = id
function _getbranchname(tree::LinkTree{RT, NL, N, B}, branch::B) where
    {RT, NL, N, B}
    id = branch.name
    branch ≡ tree.branches[id] || error("Branch $id not in tree")
    return id
end

import Phylo.API: _branchdatatype
_branchdatatype(::Type{<: LinkBranch{RT, NL, Data}}) where {RT, NL, Data} = Data
import Phylo.API: _getbranchdata
_getbranchdata(::AbstractTree, branch::LinkBranch) = branch.data
import Phylo.API: _setbranchdata!
_setbranchdata!(::AbstractTree, branch::LinkBranch{RT, NL, Data},
                data::Data) where {RT, NL, Data} = (branch.data = data)

# LinkNode methods
function LinkNode(tree::T, name::NL,
                  data::Dict{String, Any} = Dict{String, Any}()) where
    {RT, NL, N <: LinkNode, B, T <: AbstractTree{OneTree, RT, NL, N, B}}
    return N(name, missing, Vector{B}(), data)
end

import Phylo.API: _nodedatatype
_nodedatatype(::Type{LinkNode{RT, NL, Data, B}}) where {RT, NL, Data, B} = Data

import Phylo.API: _getnodedata
_getnodedata(::LinkTree, node::LinkNode) = node.data

import Phylo.API: _setnodedata!
_setnodedata!(::LinkTree, node::LinkNode{RT, NL, Data, B},
              data::Data) where {RT, NL, B, Data} = (node.data = data)

import Phylo.API: _hasinbound
_hasinbound(::LinkTree, node::LinkNode{<: Rooted}) = !ismissing(node.inbound)

import Phylo.API: _degree
_degree(::LinkTree, node::LinkNode{Unrooted}) = length(node.other)

import Phylo.API: _getinbound
_getinbound(::LinkTree, node::LinkNode) = node.inbound

import Phylo.API: _addinbound!
function _addinbound!(tree::LinkTree,
                      node::LinkNode{<: Rooted, NL, Data, B},
                      branch::B) where {NL, Data, B}
    _hasinbound(tree, node) &&
        error("LinkNode already has an inbound connection")
    node.inbound = branch
    filter!(n -> n ≠ node, tree.roots)
    tree.isvalid = missing
end

import Phylo.API: _removeinbound!
function _removeinbound!(tree::LinkTree, node::LinkNode{<: Rooted, NL, Data, B},
                         branch::B) where {NL, Data, B}
    _hasinbound(tree, node) || error("Node has no inbound connection")
    node.inbound ≡ branch ||
        error("Node has no inbound connection from branch $inbound")
    node.inbound = missing
    push!(tree.roots, node)
    tree.isvalid = missing
end

import Phylo.API: _getoutbounds
_getoutbounds(::LinkTree, node::LinkNode{<: Rooted}) = node.other

import Phylo.API: _addoutbound!
_addoutbound!(::LinkTree,
              node::LinkNode{<: Rooted, NL, Data, B},
              branch::B) where {NL, Data, B} = push!(node.other, branch)

import Phylo.API: _removeoutbound!
function _removeoutbound!(tree::LinkTree,
                          node::LinkNode{<: Rooted, NL, Data, B},
                          branch::B) where {NL, Data, B}
    branch ∈ _getoutbounds(tree, node) ? filter!(p -> p ≢ branch, node.other) :
        error("Node does not have outbound connection to branch $branch")
end

import Phylo.API: _getconnections
_getconnections(::AbstractTree, node::LinkNode{Unrooted}) = node.other

import Phylo.API: _addconnection!
function _addconnection!(::AbstractTree, node::LinkNode{Unrooted, NL, Data, B},
                         branch::B) where {NL, Data, B}
    push!(node.other, branch)
end

import Phylo.API: _removeconnection!
function _removeconnection!(tree::AbstractTree,
                            node::LinkNode{Unrooted, NL}, branch::B) where
    {NL, B <: LinkBranch{Unrooted, NL}}
    outbound ∈ _getoutbounds(tree, node) ?
        filter!(p -> p ≢ branch, node.other) :
        error("Node does not have a connection to branch $branch")
end

# LinkTree methods
const TREENAME = "Tree"
const NODENAME = "Node"

import Phylo.API: _validate!
function _validate!(tree::LinkTree{RT, NL, N, B, TD}) where {RT, NL, N, B, TD}
    tree.isvalid = true
    if TD <: Dict
        if !isempty(tree.tipdata)
            tree.isvalid &= (Set(keys(tree.tipdata)) == Set(getleafnames(tree)))
        end
    end
    nr = nroots(tree)
    if RT == OneRoot
        if nr != 1
            @warn "Wrong number of roots for $RT tree ($nr)"
            tree.isvalid = false
        end
    elseif RT == ManyRoots
        if nr < 1
            @warn "Wrong number of roots for $RT tree ($nr)"
            tree.isvalid = false
        end
    end
    return tree.isvalid
end

import Phylo.API: _treenametype
_treenametype(::Type{<: LinkTree}) = String
import Phylo.API: _gettreename
_gettreename(tree::LinkTree) = ismissing(tree.name) ? TREENAME : tree.name
import Phylo.API: _newnodelabel
function _newnodelabel(tree::LinkTree{RT, String}) where RT
    id = length(tree.nodes) + 1
    name = NODENAME * " $id"
    return haskey(tree.nodedict, name) ?
        _newlabel(collect(getnodenames(tree)), NODENAME) : name
end
import Phylo.API: _newbranchlabel
_newbranchlabel(tree::LinkTree) = length(tree.branches) + 1
import Phylo.API: _deletenode!
function _deletenode!(tree::LinkTree{RT, NL, N}, node::N) where {RT, NL, N}
    (haskey(tree.nodedict, node.name) &&
     tree.nodes[tree.nodedict[node.name]] ≡ node) ||
     error("Node $(node.name) is not in tree $(_gettreename(tree)), cannot be deleted")
    # Does nothing for unrooted tree, as inbound is never set
    !ismissing(node.inbound) && _deletebranch!(tree, node.inbound)
    # Delete outbound connections of rooted tree, all connections of unrooted
    while _degree(tree, node) > 0
        _deletebranch!(tree, first(node.other))
    end
    tree.nodes[tree.nodedict[node.name]] = missing
    delete!(tree.nodedict, node.name)
    filter!(n -> n ≢ node, tree.roots)
    tree.isvalid = missing
    return true
end

import Phylo.API: _getroots
_getroots(tree::LinkTree{<: Rooted}) = tree.roots
_getroots(tree::LinkTree{Unrooted}) = error("Unrooted trees do not have roots")

import Phylo.API: _getnodes
_getnodes(tree::LinkTree) = skipmissing(tree.nodes)

import Phylo.API: _getnodenames
_getnodenames(tree::LinkTree) = keys(tree.nodedict)

import Phylo.API: _getnodename
_getnodename(::LinkTree{RT, NL, N}, node::N) where {RT, NL, N} = node.name

import Phylo.API: _hasnode
_hasnode(tree::LinkTree{RT, NL}, name::NL) where {RT, NL} =
    haskey(tree.nodedict, name)
_hasnode(tree::LinkTree{RT, NL, N}, node::N) where {RT, NL, N} =
    haskey(tree.nodedict, node.name)

import Phylo.API: _getnode
@traitfn _getnode(tree::T,
                  node::NL) where {RT, NL,
                                   T <: LinkTree{RT, NL};
                                   PreferNodeObjects{T}} =
                                       tree.nodes[tree.nodedict[node]]

import Phylo.API: _getbranches
_getbranches(tree::LinkTree) = skipmissing(tree.branches)

import Phylo.API: _getbranchnames
_getbranchnames(tree::LinkTree) =
    [getbranchname(tree, b) for b in tree.branches if !ismissing(b)]

import Phylo.API: _hasbranch
_hasbranch(tree::LinkTree, id::Int) =
    1 ≤ id ≤ length(tree.branches) && !ismissing(tree.branches[id])
_hasbranch(tree::LinkTree{RT, NL, N, B}, branch::B) where
{RT, NL, N <: LinkNode{RT, NL}, B <: LinkBranch{RT, NL}} =
    branch ∈ skipmissing(tree.branches)

import Phylo.API: _getbranch
_getbranch(tree::LinkTree, id::Int) = tree.branches[id]

import Phylo.API: _createbranch!
function _createbranch!(tree::LinkTree{RT, NL, N, B},
                        from::N, to::N, len::Union{LenUnits, Missing} = missing,
                        name::Union{Int, Missing} = missing,
                        data::Data = newempty(Data)) where
    {RT <: Rooted, NL, N, Data, LenUnits,
     B <: LinkBranch{RT, NL, Data, LenUnits}}
    id = ismissing(name) ? length(tree.branches) + 1 : name
    branch = LinkBranch(id, from, to, len, data)
    if ismissing(name)
        push!(tree.branches, branch)
    else
        l = length(tree.branches)
        if id > l
            resize!(tree.branches, id)
            tree.branches[l+1:id-1] .= missing
        else
            ismissing(tree.branches[id]) || error("Branch $name already exists")
        end
        tree.branches[id] = branch
    end
    _addoutbound!(tree, from, branch)
    _addinbound!(tree, to, branch)
    tree.isvalid = missing
    return branch
end
function _createbranch!(tree::LinkTree{Unrooted, NL, N, B},
                        from::N, to::N, len::Union{LenUnits, Missing} = missing,
                        name::Union{Int, Missing} = missing,
                        data::Data = newempty(Data)) where
    {NL, N, Data, LenUnits,
     B <: LinkBranch{Unrooted, NL, Data, LenUnits}}
    id = ismissing(name) ? length(tree.branches) + 1 : name
    branch = LinkBranch(id, tree, from, to, len, data)
    push!(tree.branches, branch)
    _addconnection!(tree, from, branch)
    _addconnection!(tree, to, branch)
    tree.isvalid = missing
    return branch
end

import Phylo.API: _deletebranch!
function _deletebranch!(tree::LinkTree{<: Rooted}, id::Int)
    branch = tree.branches[id]
    _removeoutbound!(tree, branch.inout[1], branch)
    _removeinbound!(tree, branch.inout[2], branch)
    tree.branches[id] = missing
    tree.isvalid = missing
    return true
end
function _deletebranch!(tree::LinkTree{Unrooted}, id::Int)
    branch = tree.branches[id]
    _removeconnection!(tree, branch.inout[1], branch)
    _removeconnection!(tree, branch.inout[2], branch)
    tree.branches[id] = missing
    tree.isvalid = missing
    return true
end
function _deletebranch!(tree::LinkTree{<: Rooted, NL, N, B}, branch::B) where
    {NL, N, B}
    _removeoutbound!(tree, branch.inout[1], branch)
    _removeinbound!(tree, branch.inout[2], branch)
    id = _getbranchname(tree, branch)
    tree.branches[id] = missing
    tree.isvalid = missing
    return true
end
function _deletebranch!(tree::LinkTree{Unrooted, NL, N, B}, branch::B) where
    {NL, N, B}
    _removeconnection!(tree, branch.inout[1], branch)
    _removeconnection!(tree, branch.inout[2], branch)
    id = _getbranchname(tree, branch)
    tree.branches[id] = missing
    tree.isvalid = missing
    return true
end

import Phylo.API: _createnode!
function _createnode!(tree::LinkTree{RT, NL, N, B}, name::Union{NL, Missing},
                      data::Data = newempty(Data)) where
    {RT <: Rooted, NL, Data, B, N <: LinkNode{RT, NL, Data, B}}
    nodename = ismissing(name) ? _newnodelabel(tree) : name
    _hasnode(tree, nodename) && error("Node $nodename already exists in tree.")
    node = LinkNode{RT, NL, Data, B}(nodename, data)
    id = length(tree.nodes) + 1
    push!(tree.nodes, node)
    tree.nodedict[nodename] = id
    push!(tree.roots, node)
    return node
end
function _createnode!(tree::LinkTree{RT, NL, N, B}, name::Union{NL, Missing},
                      data::Data = newempty(Data)) where
    {RT <: Unrooted, NL, Data, B, N <: LinkNode{RT, NL, Data, B}}
    nodename = ismissing(name) ? _newnodelabel(tree) : name
    _hasnode(tree, nodename) && error("Node $nodename already exists in tree.")
    node = LinkNode{RT, NL, Data, B}(nodename, data)
    id = length(tree.nodes) + 1
    push!(tree.nodes, node)
    tree.nodedict[nodename] = id
    return node
end

import Phylo.API: _leafinfotype
_leafinfotype(::Type{<:LinkTree{RT, NL, N, B, TD}}) where {RT, NL, N, B, TD} =
    TD

import Phylo.API: _getleafinfo
_getleafinfo(tree::LinkTree) = tree.tipdata

import Phylo.API: _setleafinfo!
_setleafinfo!(tree::LinkTree{RT, NL, N, B, TD}, leafinfo::TD) where
{RT, NL, N, B, TD} = tree.tipdata = leafinfo

import Phylo.API: _nodedatatype
_nodedatatype(::Type{<:LinkTree{RT, NL, N}}) where {RT, NL, N} =
    _nodedatatype(N)

import Phylo.API: _branchdatatype
_branchdatatype(::Type{<:LinkTree{RT, NL, N, B}}) where {RT, NL, N, B} =
    _branchdatatype(B)

import Phylo.API: _nnodes
_nnodes(tree::LinkTree) = count((!)∘ismissing, tree.nodes)
import Phylo.API: _nbranches
_nbranches(tree::LinkTree) = count((!)∘ismissing, tree.branches)

import Base.show
function show(io::IO, node::LinkNode{Unrooted})
    print(io, "Unrooted LinkNode '$(node.name)', with $(length(node.other)) connection(s).")
    if length(node.other) == 0
        println(io, "a node with no connections.")
    elseif length(node.other) == 1
        println(io, "a node with 1 connection (branch $(node.other[1].name))")
    else
        print(io, "a node with $(length(node.other)) outbound connections (branches $(node.other[1].name)")
        for i in 2:length(node.other)
            print(io, ", $(node.other[i].name)")
        end
        println(io, ")")
    end
end

function show(io::IO, node::LinkNode)
    print(io, "LinkNode $(node.name), ")
    if ismissing(node.inbound)
        if length(node.other) == 0
            println(io, "an isolated node with no connections.")
        elseif length(node.other) == 1
            println(io, "a root node with 1 outbound connection (branch $(node.other[1].name))")
        else
            print(io, "a root node with $(length(node.other)) outbound connections (branches $(node.other[1].name)")
            for i in 2:length(node.other)
                print(io, ", $(node.other[i].name)")
            end
            println(io, ")")
        end
    else
        if length(node.other) == 0
            println(io, "a tip of the tree with an incoming connection (branch $(node.inbound.name)).")
        elseif length(node.other) == 1
            println(io, "an internal node with 1 inbound and 1 outbound connection (branches $(node.inbound.name) and $(node.other[1].name))")
        else
            print(io, "an internal node with 1 inbound and $(length(node.other)) outbound connections (branches $(node.inbound.name) and $(node.other[1].name)")
            for i in 2:length(node.other)
                print(io, ", $(node.other[i].name)")
            end
            println(io, ")")
        end
    end
end

function show(io::IO, branch::LinkBranch{Unrooted})
    println(io, "Unrooted LinkBranch $(branch.name), connecting nodes $(branch.inout[1].name) and $(branch.inout[2].name)$(ismissing(branch.length) ? "" : " (length $(branch.length))").")
end

function show(io::IO, branch::LinkBranch)
    println(io, "LinkBranch $(branch.name), from node $(branch.inout[1].name) to node $(branch.inout[2].name)$(ismissing(branch.length) ? "" : " (length $(branch.length))").")
end


text:
"""
    append_hs!(model::JuMP.Model, x, nh::Int)

Append equalty constraints (auto-generated)

# Arguments
    `model::JuMP.Model`: optimization model
    `x::JuMP.variables`: optimization model variables
    `nh::Int`: number of equality constraints
"""
function append_hs!(model::JuMP.Model, x, nh::Int)
    for ih = 1:nh
        if ih == 1
            @NLconstraint(model, h1(x...) == 0)
        elseif ih == 2
            @NLconstraint(model, h2(x...) == 0)
        elseif ih == 3
            @NLconstraint(model, h3(x...) == 0)
        elseif ih == 4
            @NLconstraint(model, h4(x...) == 0)
        elseif ih == 5
            @NLconstraint(model, h5(x...) == 0)
        elseif ih == 6
            @NLconstraint(model, h6(x...) == 0)
        elseif ih == 7
            @NLconstraint(model, h7(x...) == 0)
        elseif ih == 8
            @NLconstraint(model, h8(x...) == 0)
        elseif ih == 9
            @NLconstraint(model, h9(x...) == 0)
        elseif ih == 10
            @NLconstraint(model, h10(x...) == 0)
        elseif ih == 11
            @NLconstraint(model, h11(x...) == 0)
        elseif ih == 12
            @NLconstraint(model, h12(x...) == 0)
        elseif ih == 13
            @NLconstraint(model, h13(x...) == 0)
        elseif ih == 14
            @NLconstraint(model, h14(x...) == 0)
        elseif ih == 15
            @NLconstraint(model, h15(x...) == 0)
        elseif ih == 16
            @NLconstraint(model, h16(x...) == 0)
        elseif ih == 17
            @NLconstraint(model, h17(x...) == 0)
        elseif ih == 18
            @NLconstraint(model, h18(x...) == 0)
        elseif ih == 19
            @NLconstraint(model, h19(x...) == 0)
        elseif ih == 20
            @NLconstraint(model, h20(x...) == 0)
        end
    end
end

max_stars_repo_name: Yuricst/NLPSaUT

max_issues_repo_name: Yuricst/NLPSaUT

max_forks_repo_name: Yuricst/NLPSaUT


text:
module Rand
    export Crand

    mutable struct Crand
        iw::Array{Int64,1}
        jr::Int64
        kr::Int64
    end

    mutable struct Random_LCGs{T}
        crand_public::Crand
        Random_LCGs(ndelay;julian_random_number=false) = new{julian_random_number}(cint3(ndelay))
    end

    function (g::Random_LCGs{false})()
        #return rand()
        return ranf!(g.crand_public)
    end

    function (g::Random_LCGs{true})()
        return rand()
        #return ranf!(g.crand_public)
    end

    const IP = 521
    const IQ = 32
    const MACRM = 40
    const MACRI = 1
    const FNORM=Float32(0.465661e-9)

    const intmax32 = Int32(2147483647)

    function cint3(ndelay)
        println("ndelay : $ndelay")
        #ccall((:readseed,"./randfortran.so"),Nothing,(),)
        crand_public = init3()
        delay3!(crand_public,ndelay)
        #display(crand_public.iw)

        return crand_public
    end

    function make_value_Int32_again(ix)
        while ix < -intmax32 - 1
            ix += 2*intmax32 + 1
        end
        while ix > intmax32
            ix -= 2*(intmax32 +1)
        end
        return ix
    end

    function init3()

        #println(2147483647+1-2*(intmax32+1))
        ib = zeros(Int64,IP) #INTEGER ib(0:IP-1)
        iw = zeros(Int64,IP) #INTEGER iw(0:IP-1)

        ix = Int32(MACRI)
        for i=0:IP-1
            ix = ix*69069
            ix = make_value_Int32_again(ix)

            #while ix < -intmax32-1 || intmax32 < ix
            #    ix -= Int32(Int32(2)*(intmax32+Int32(1)))
            #end
            #println(Int32(ix) >>> -(-31))
            ib[i+1] = Int32(ix) >>> -(-31)
            #println("$i ",ix," ",typeof(ix))
            #println(ix,"\t",ib[i+1])
            #ib[i+1] = ishift(ix,-31)
        end

        jr = 0
        kr = IP-IQ
        for j=0:IP-1
            iwork = 0
            for i=0:31
                iwork=iwork*2+ib[jr+1]
                #while iwork < -intmax32-1 || intmax32 < iwork
                #    iwork -= Int32(Int32(2)*(intmax32+Int32(1)))
                #end
                iwork = make_value_Int32_again(iwork)

                ib[jr+1] = Int32(ib[jr+1]) ⊻ Int32(ib[kr+1])
                jr += 1
                if jr == IP
                    jr = 0
                end
                kr += 1
                if kr == IP
                    kr = 0
                end
            end
            iw[j+1] = Int32(iwork) >>> -(-1)
            #println(j,"\t",iw[j+1])
        end

        return Crand(iw,jr,kr)

#        exit()

    end

    function delay3!(crand_public::Crand,lambda)

        iwk = zeros(Int32,2*IP-2+1)
        c = zeros(Int32,2*IP-1+1)
        ib = zeros(Int32,IP+31+1)

        mu = MACRM
        for i=0:IP-1
            iwk[i+1] = crand_public.iw[i+1]
        end

        for i=IP:2*IP-2
            iwk[i+1] = Int32(iwk[i-IP+1]) ⊻ Int32(iwk[i-IQ+1])
        end

        for i=0:mu-1
            ib[i+1] = 0
        end

        m = lambda
        nb = mu-1

        #count = 0
        while m>IP-1
            nb += 1
            ib[nb+1] = m % 2
            m = div(m,2)
            #count += 1
            #if count > 100000000
            #    error("something is wrong in rand.jl!")
            #end
        end

        for i=0:IP-1
            c[i+1] = 0
        end

        c[m+1] = 1

        for j=nb:-1:0
            for i=IP-1:-1:0
                c[2*i+ib[j+1]+1] = c[i+1]
                c[2*i+1-ib[j+1]+1] = 0
            end

            for i=2*IP-1:-1:IP
                c[i-IP+1] = Int32(c[i-IP+1]) ⊻ Int32(c[i+1])
                c[i-IQ+1] = Int32(c[i-IQ+1]) ⊻ Int32(c[i+1])
            end
        end

        for j=0:IP-1
            iwork = 0
            for i=0:IP-1
                #while iwork < -intmax32-1 || intmax32 < iwork
                #    iwork -= Int32(Int32(2)*(intmax32+Int32(1)))
                #end
                iwork = make_value_Int32_again(iwork)
                iwork = Int32(iwork) ⊻ Int32(c[i+1]*iwk[j+i+1])
            end
            crand_public.iw[j+1] = iwork
        end

    end

    function ranf!(crand_public)
        #ranfs = zeros(Float64,1)
        #ccall((:call_random,"./randfortran.so"),Nothing,(Ref{Float64},),ranfs)
        #return ranfs[1]

        ir = crand_public.iw
        j = crand_public.jr
        k = crand_public.kr
        crand_public.iw[j+1] = Int32(crand_public.iw[j+1]) ⊻ Int32(crand_public.iw[k+1])
        irnd = crand_public.iw[j+1]
        #println("irnd = $irnd ", FNORM,"\t",irnd*FNORM)
        #ranf = Float32(irnd*FNORM)
        #granf = Float64(irnd)*FNORM
        frnd = Float32(irnd*FNORM)
        #println(frnd,"\t",irnd*FNORM,"\t",Float64(irnd),"\t",granf)
        #exit()
        ranf = Float32(frnd)

        crand_public.jr = crand_public.jr + 1
        if crand_public.jr >= IP
            crand_public.jr = crand_public.jr - IP
        end
        crand_public.kr = crand_public.kr + 1
        if crand_public.kr >= IP
            crand_public.kr = crand_public.kr - IP
        end
        return ranf
    end

    function rgauss1!(crand_public,nup)
        xr = zeros(Float64,nup)

        nup2 = div(nup,2)
        if nup != 2*nup2
            error(
                """
                !!!!  Attention nup in rgauss should be even nup = $nup
                """
            )
        end

        for i=1:nup2
            i2 = 2*i
            i1 = i2 - 1
            v1 = sqrt(-log(ranf!(crand_public)+1e-10)*2)
            v2 = 2pi*ranf!(crand_public)

            xr[i1] = v1*cos(v2)
            xr[i2] = v1*sin(v2)

        end

        return xr
    end
end


text:
#    Problem : GROUP A
#    *********
#	 This problem is a sum of n+1 least-squares groups, the first n of
#    which have only a linear element.
#    It Hessian matrix is dense.
#
#    Origonal SIF Source: Problem 23 in
#    J.J. More', B.S. Garbow and K.E. Hillstrom,
#    "Testing Unconstrained Optimization Software",
#    ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    PENALTY1.jl.SIF classification SUR2-AN-V-0
#
#    Number of variables is variable

f = (x) -> begin
	a  = 1e-5
	fx = tail = 0.0
	for i in firstindex(x):lastindex(x)
		item  = x[i] - 1
      	fx   += item^2*a
      	tail += x[i]^2
	end
	fx += (tail - 0.25)^2
    return fx
end

g! = (g, x) -> begin
	a  = 1e-5
	tail = 0.0
	for i in firstindex(x):lastindex(x)
		item  = x[i] - 1
      	g[i]  = 2a*item
      	tail += x[i]^2
	end
	for i in firstindex(x):lastindex(x)
		g[i] += 4*(tail-0.25)x[i]
	end
    return g
end

fg! = (g, x) -> begin
	a  = 1e-5
	tail = 0.0
	for i in firstindex(x):lastindex(x)
		item  = x[i] - 1
      	fx   += item^2*a
      	g[i]  = 2a*item
      	tail += x[i]^2
	end
	fx += (tail - 0.25)^2
	for i in firstindex(x):lastindex(x)
		g[i] += 4*(tail-0.25)x[i]
	end
    return fx, g
end

init = (n::Int=1000) -> begin
	x0 = Vector(1.0:1.0:n)
    return n, x0
end

TestSet["PENALTY1"] = UncProgram("PENALTY1", f, g!, fg!, init)


text:
@testset "evaluation.jl" begin

    @testset "SameInstancesEvaluator constructor" begin
        generator = SeaPearl.HomogenousGraphColoringGenerator(10, 0.1)

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        valueheuristic = SeaPearl.BasicHeuristic(my_heuristic)
        valueSelectionArray = [valueheuristic]
        evaluator = SeaPearl.SameInstancesEvaluator(valueSelectionArray,generator; seed=nothing, evalFreq = 50, nbInstances = 2)

        @test evaluator.nbInstances == 2
        @test evaluator.evalFreq == 50
        @test length(evaluator.instances) == 2
        @test length(values(evaluator.instances[1].variables)) == 11
        @test length(values(evaluator.instances[2].variables)) == 11
    end

    @testset "evaluate(::SameInstancesEvaluator)" begin

        generator = SeaPearl.HomogenousGraphColoringGenerator(10, 0.1)

        my_heuristic(x::SeaPearl.IntVar; cpmodel=nothing) = minimum(x.domain)
        valueheuristic = SeaPearl.BasicHeuristic(my_heuristic)
        valueSelectionArray = [valueheuristic]
        evaluator = SeaPearl.SameInstancesEvaluator(valueSelectionArray,generator; seed=1, evalFreq = 50, nbInstances = 2)

        variableheuristic = SeaPearl.MinDomainVariableSelection{false}()

        SeaPearl.evaluate(evaluator, variableheuristic, SeaPearl.DFSearch())
        @test evaluator.metrics[1,1].nodeVisited[1] == [12, 23]
    end
end


text:
@mkCppFun IsRoleAllowedAsHand
@mkCppFun ButtonMaskFromId
@mkCppFun VRToken
@mkCppFun OpenVRInternal_ModuleContext
@mkCppFun VRSystem
@mkCppFun VRChaperone
@mkCppFun VRChaperoneSetup
@mkCppFun VRCompositor
@mkCppFun VROverlay
@mkCppFun VRScreenshots
@mkCppFun VRRenderModels
@mkCppFun VRApplications
@mkCppFun VRSettings
@mkCppFun VRResources
@mkCppFun VRExtendedDisplay
@mkCppFun VRTrackedCamera
@mkCppFun VRDriverManager
@mkCppFun VRInput
@mkCppFun VRIOBuffer
@mkCppFun VRSpatialAnchors
# # @mkCppFun VRNotifications
@mkCppFun VR_Init
@mkCppFun VR_Shutdown


text:
function parse_elem!(file::AbstractString, ex::Expr, dest::Set{String})
    if ex.head == :call && ex.args[1] == :include && ex.args[2] isa String
        dir = dirname(file)
        include_file = normpath(joinpath(dir, ex.args[2]))
        parse_includes!(include_file, dest)
        return
    end
    for arg in ex.args
        parse_elem!(file, arg, dest)
    end
end
function parse_elem!(file::AbstractString, ex, dest::Set{String}) end

function parse_includes!(file::AbstractString, dest::Set{String})
    !isfile(file) && return
    file in dest && return
    push!(dest, abspath(file))

    ex = Base.parse_input_line(read(file, String); filename=file)
    parse_elem!(file, ex, dest)

end

function parse_includes(file::AbstractString)
    result = Set{String}()
    parse_includes!(file, result)
    return result
end


text:
function f_spl(msap_prop, rho_0, c_0)

    # Compute SPL
    spl = 10*log10.(msap_prop) .+ 20*log10.(rho_0 .* c_0.^2) #.* ones(eltype(msap_prop), (1, settings.N_f))

    # Remove all SPL below 0
    spl = clamp.(spl, 0, Inf)

    return spl
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA