text:
import Base: checkindex, checkbounds_indices, DimOrInd, OneTo, Slice

@inline _similar(A::AbstractArray, ::Type{T}=eltype(A)) where {T} = similar(A, T, axes(A))
@inline _similar(A::Number, ::Type{T}=typeof(A)) where {T}        = Scalar{T}()

struct Work{T,N<:AbstractLegendreNorm,V<:AbstractArray{T}}
    norm::N
    z::V
    y¹::V
    y²::V
    pmm2::V
    pm::V
    plm2::V
    plm1::V
    pl::V
end
function Work(norm::AbstractLegendreNorm, Λ, x)
    T = promote_type(eltype(norm), eltype(Λ), eltype(x))
    z    = _similar(x, T)
    y¹   = _similar(x, T)
    y²   = _similar(x, T)
    pmm2 = _similar(x, T)
    pm   = _similar(x, T)
    plm2 = _similar(x, T)
    plm1 = _similar(x, T)
    pl   = _similar(x, T)
    return Work(norm, z, y¹, y², pmm2, pm, plm2, plm1, pl)
end

function _chkdomain(lmax, mmax)
    @noinline _chkdomain_throw_lmax(l) = throw(DomainError(l, "degree lmax must be non-negative"))
    @noinline _chkdomain_throw_mmax(m) = throw(DomainError(m, "order mmax must be non-negative and less than or equal to lmax"))

    0 ≤ lmax || _chkdomain_throw_lmax(lmax)
    0 ≤ mmax ≤ lmax || _chkdomain_throw_mmax(mmax)
    nothing
end

function _chkbounds(Λ, lmax, mmax, x)
    @noinline _chkbounds_throw_dims(M, N) = throw(DimensionMismatch(
        "Output has $M dimensions, expected $N to $(N+2)"))
    @noinline _chkbounds_throw_axes(Λ, x) = begin
        throw(DimensionMismatch(
            "Output has leading axes $(ntuple(i -> axes(Λ,i), ndims(x))), expected $(axes(x))"))
    end
    @noinline _chkbounds_throw_lmax() = throw(DimensionMismatch(
            "lmax incompatible with output array axes"))
    @noinline _chkbounds_throw_mmax() = throw(DimensionMismatch(
            "mmax incompatible with output array axes"))

    M = ndims(Λ)
    N = ndims(x)
    (0 ≤ M - N ≤ 2) || _chkbounds_throw_dims(M, N)
    # Leading dimensions of Λ are storage for the same dimensions as x
    axesΛ = axes(Λ)
    if N > 0
        axΛ = ntuple(i -> axesΛ[i], N)
        axx = axes(x)
        checkbounds_indices(Bool, axΛ, axx) || _chkbounds_throw_axes(Λ, x)
    end
    # Trailing dimensions of Λ are storage for range of ell and m, as needed
    dimΛ = ntuple(i -> axesΛ[i+N], M-N)
    if M - N > 0
        checkindex(Bool, dimΛ[1], OneTo(lmax+1)) || _chkbounds_throw_lmax()
    end
    if M - N > 1
        checkindex(Bool, dimΛ[2], OneTo(mmax+1)) || _chkbounds_throw_mmax()
    end
    nothing
end

function unsafe_legendre!(workornorm::Union{AbstractLegendreNorm,Work}, Λ, lmax, mmax, x)
    if ndims(x) > 1
        M = ndims(Λ)
        N = ndims(x)
        S = prod(size(x))
        x′ = reshape(x, S)
        Λ′ = reshape(Λ, S, ntuple(i->size(Λ,i+N), M-N)...)
    else
        x′ = x
        Λ′ = Λ
    end

    if workornorm isa AbstractLegendreNorm
        work = Work(workornorm, Λ′, x′)
        @inbounds _legendre_impl!(work, Λ′, lmax, mmax, x′)
    else
        @inbounds _legendre_impl!(workornorm, Λ′, lmax, mmax, x′)
    end
    return Λ
end

# Similar to `Base.fma` but defined for `Complex` arguments such that if all arguments are
# on the real axis, the result of `_fma(x, y, z) === _fma(real(x), real(y), real(z)) + 0im`.
#
# !!! note
#     The total operation for generic complex arguments is not fused but rather
#     proceeds via a series of FMA instructions (which include intermediate
#     rounding) --- the fusing is only strictly true for the all-real-axis case.
#
# Following definitions should be complete, but we comment out most of them since they
# are unused.
#_fma(z::Complex, w::Complex, x::Complex) =
#    Complex(_fma(real(z), real(w), -_fma(imag(z), imag(w), -real(x))),
#            _fma(real(z), imag(w),  _fma(imag(z), real(w),  imag(x))))
_fma(z::Complex, w::Complex, x::Real) =
    Complex(_fma(real(z), real(w), -_fma(imag(z), imag(w), -x)),
            _fma(real(z), imag(w), imag(z) * real(w)))
#_fma(x::Real, z::Complex, y::Number) = _fma(z, x, y)
#_fma(z::Complex, x::Real, y::Real) = Complex(_fma(real(z),x,y), imag(z)*x)
#_fma(z::Complex, x::Real, w::Complex) =
#    Complex(_fma(real(z),x,real(w)), _fma(imag(z),x,imag(w)))
#_fma(x::Real, y::Real, z::Complex) = Complex(_fma(x,y,real(z)), imag(z))
# Then cascade back to Base.fma
_fma(x, y, z) = Base.fma(x, y, z)

@propagate_inbounds function _legendre_impl!(work::Work, Λ, lmax, mmax, x)
    norm = work.norm
    z    = work.z
    y¹   = work.y¹
    y²   = work.y²
    pmm2 = work.pmm2
    pm   = work.pm
    plm2 = work.plm2
    plm1 = work.plm1
    pl   = work.pl

    local μ₁, μ₂
    T = eltype(z)
    M = ndims(x)
    N = ndims(Λ) - M
    I = CartesianIndices(map(Slice, axes(x)))

    @simd for ii in I
        z[ii]  = convert(T, x[ii])
        y²[ii] = -_fma(z[ii], z[ii], -one(real(T)))
        y¹[ii] = sqrt(y²[ii])
    end

    fill!(pm, initcond(norm, T))
    for m in 0:mmax
        @simd for ii in I
            if N == 2
                Λ[ii,m+1,m+1] = pm[ii]
            elseif N == 1 && m == mmax
                Λ[ii,m+1] = pm[ii]
            elseif N == 0 && m == lmax # == mmax
                Λ[ii] = pm[ii]
            end
        end
        m == lmax && break # exit for lmax == mmax == m

        if N == 2 || m == mmax
            # 1-term recurrence relation taking (l-1,l-1) -> (l,l-1) where l == m == m
            l = m + 1
            ν = coeff_ν(norm, real(T), l)
            @simd for ii in I
                plm1[ii] = pm[ii]
                pl[ii]   = ν * z[ii] * plm1[ii]
                if N == 2
                    Λ[ii,l+1,m+1] = pl[ii]
                elseif N == 1
                    Λ[ii,l+1] = pl[ii]
                elseif N == 0 && l == lmax
                    Λ[ii] = pl[ii]
                end
            end
            # no exit since following loop will be skipped for
            #   lmax == mmax + 1

            for l in m+2:lmax
                plm2, plm1, pl = plm1, pl, plm2
                # 2-term recurrence relation taking (l-1,m) -> (l, m)
                α = coeff_α(norm, real(T), l, m)
                β = coeff_β(norm, real(T), l, m)
                @simd for ii in I
                    pl[ii] = α * z[ii] * plm1[ii] - β * plm2[ii]
                    if N == 2
                        Λ[ii,l+1,m+1] = pl[ii]
                    elseif N == 1
                        Λ[ii,l+1] = pl[ii]
                    elseif N == 0 && l == lmax
                        Λ[ii] = pl[ii]
                    end
                end
            end
        end
        m == mmax && break # skip calculating (m,m) where m == mmax + 1

        # 1-term recurrence relation taking (m,m) -> (m+1,m+1)
        if iseven(m)
            pmm2, pm = pm, pmm2
            # Takes even m-2 to odd m-1 for following iteration where m will be odd
            μ₁ = coeff_μ(norm, real(T), m+1)
            @simd for ii in I
                pm[ii] = -μ₁ * y¹[ii] * pmm2[ii]
            end
        else
            # Takes even m-2 to even m for following iteration where m will be even again
            μ₂ = coeff_μ(norm, real(T), m+1)
            @simd for ii in I
                pm[ii] = μ₁ * μ₂ * y²[ii] * pmm2[ii]
            end
        end
    end

    return Λ
end

"""
    legendre!(norm::AbstractLegendreNorm, Λ, l::Integer, m::Integer, x)

Fills the array `Λ` with the Legendre polynomial values ``N_ℓ^m P_ℓ^m(x)``
up to/of degree(s) `l` and order(s) `m` for the normalization scheme `norm`.
`Λ` must be an array with between 0 and 2 more dimensions than `x`, with the leading
dimensions having the same shape as `x`.

- If `ndims(Λ) == ndims(x)`, then `Λ` is filled with the polynomial values at `x` for
  degree `l` and order `m`.
- If `ndims(Λ) == ndims(x) + 1`, then `l` is interpreted as `lmax`, and `Λ` filled with
  polynomial values for all degrees `0 ≤ l ≤ lmax` of order `m`.
- If `ndims(Λ) == ndims(x) + 2`, then `l` is interpreted as `lmax` and `m` as `mmax`,
  and `Λ` is filled with polynomial values for all degrees `0 ≤ l ≤ lmax` and orders
  `0 ≤ m ≤ min(mmax, l)`.
"""
function legendre!(norm::AbstractLegendreNorm, Λ, l::Integer, m::Integer, x)
    _chkdomain(l, m)
    boundscheck_hook(norm, l, m)
    _chkbounds(Λ, l, m, x)
    return unsafe_legendre!(norm, Λ, l, m, x)
end

"""
    p = legendre(norm::AbstractLegendreNorm, l::DimOrInd, m::DimOrInd, x)

Computes the associated Legendre polynomials ``N_ℓ^m P_ℓ^m(x)`` of degree(s) `l` and
order(s) `m` at `x` for the normalization scheme `norm`.

- If `l` and `m` are integers, returns a single value.
- If `l` is a range `0:lmax` and `m` an integer, returns the vector of values for order
  `m` and all degrees `0 ≤ l ≤ lmax`.
- If `l` is a range `0:lmax` and `m` is a range `0:mmax`, returns the matrix of values
  for all degrees `0 ≤ l ≤ lmax` and orders `0 ≤ m ≤ mmax`.

Note that in both the second and third cases, the ranges must have a first index of 0.
"""
function legendre end

# Scalar argument, scalar output case handled separately.
function legendre(norm::AbstractLegendreNorm, l::Integer, m::Integer, x::Number)
    _chkdomain(l, m)
    boundscheck_hook(norm, l, m)
    Λ = _similar(x)
    unsafe_legendre!(norm, Λ, l, m, x)
    return Λ[]
end

function legendre(norm::AbstractLegendreNorm, l::DimOrInd, m::DimOrInd, x)
    if l isa AbstractUnitRange
        first(l) == 0 || throw(ArgumentError("Range of degrees l must start at 0"))
    end
    if m isa AbstractUnitRange
        if !(l isa AbstractUnitRange)
            throw(ArgumentError("Range of orders m requires range of degrees l"))
        end
        first(m) == 0 || throw(ArgumentError("Range of orders m must start at 0"))
    end
    lmax = l isa AbstractUnitRange ? last(l) : l
    mmax = m isa AbstractUnitRange ? last(m) : m
    _chkdomain(lmax, mmax)
    boundscheck_hook(norm, lmax, mmax)

    Λ = zeros(eltype(x), axes(x)..., size(l)..., size(m)...)
    unsafe_legendre!(norm, Λ, lmax, mmax, x)
    return Λ
end

# Make normalizations callable with similar syntax as the legendre[!] functions
@static if VERSION < v"1.3.0-alpha"
    # Prior to julia-1.3.0, abstract types could not have methods attached to them.
    # Provide for just the defined normalization types.
    for N in (LegendreUnitNorm,LegendreSphereNorm,LegendreOrthoNorm,LegendreFourPiNorm,LegendreNormCoeff)
        @eval begin
            @inline function (norm::$N)(l, m, x)
                return legendre(norm, l, m, x)
            end
            @inline function (norm::$N)(Λ, l::Integer, m::Integer, x)
                return legendre!(norm, Λ, l, m, x)
            end
        end
    end
else
    @inline function (norm::AbstractLegendreNorm)(l, m, x)
        return legendre(norm, l, m, x)
    end
    @inline function (norm::AbstractLegendreNorm)(Λ, l::Integer, m::Integer, x)
        return legendre!(norm, Λ, l, m, x)
    end
end


text:
################################################
#
# minmax heap type and constructors
#
################################################

mutable struct BinaryMinMaxHeap{T} <: AbstractMinMaxHeap{T}
    valtree::Vector{T}

    BinaryMinMaxHeap{T}() where {T} = new{T}(Vector{T}())

    function BinaryMinMaxHeap{T}(xs::AbstractVector{T}) where {T}
        valtree = _make_binary_minmax_heap(xs)
        new{T}(valtree)
    end
end

BinaryMinMaxHeap(xs::AbstractVector{T}) where T = BinaryMinMaxHeap{T}(xs)

################################################
#
# core implementation
#
################################################

function _make_binary_minmax_heap(xs)
    valtree = copy(xs)
    for i in length(xs):-1:1
        @inbounds _minmax_heap_trickle_down!(valtree, i)
    end
    return valtree
end

Base.@propagate_inbounds function _minmax_heap_bubble_up!(A::AbstractVector, i::Integer)
    if on_minlevel(i)
        if i > 1 && A[i] > A[hparent(i)]
            # swap to parent and bubble up max
            tmp = A[i]
            A[i] = A[hparent(i)]
            A[hparent(i)] = tmp
            _minmax_heap_bubble_up!(A, hparent(i), Reverse)
        else
            # bubble up min
            _minmax_heap_bubble_up!(A, i, Forward)
        end

    else
        # max level
        if i > 1 && A[i] < A[hparent(i)]
            # swap to parent and bubble up min
            tmp = A[i]
            A[i] = A[hparent(i)]
            A[hparent(i)] = tmp
            _minmax_heap_bubble_up!(A, hparent(i), Forward)
        else
            # bubble up max
            _minmax_heap_bubble_up!(A, i, Reverse)
        end
    end
end

Base.@propagate_inbounds function _minmax_heap_bubble_up!(A::AbstractVector, i::Integer, o::Ordering, x=A[i])
    if hasgrandparent(i)
        gparent = hparent(hparent(i))
        if lt(o, x, A[gparent])
            A[i] = A[gparent]
            A[gparent] = x
            _minmax_heap_bubble_up!(A, gparent, o)
        end
    end
end

Base.@propagate_inbounds function _minmax_heap_trickle_down!(A::AbstractVector, i::Integer)
    if on_minlevel(i)
        _minmax_heap_trickle_down!(A, i, Forward)
    else
        _minmax_heap_trickle_down!(A, i, Reverse)
    end
end

Base.@propagate_inbounds function _minmax_heap_trickle_down!(A::AbstractVector, i::Integer, o::Ordering, x=A[i])

    if haschildren(i, A)
        # get the index of the extremum (min or max) descendant
        extremum = o === Forward ? minimum : maximum
        _, m = extremum((A[j], j) for j in children_and_grandchildren(length(A), i))

        if isgrandchild(m, i)
            if lt(o, A[m], A[i])
                A[i] = A[m]
                A[m] = x
                if lt(o, A[hparent(m)], A[m])
                    t = A[m]
                    A[m] = A[hparent(m)]
                    A[hparent(m)] = t
                end
                _minmax_heap_trickle_down!(A, m, o)
            end
        else
            if lt(o, A[m], A[i])
                A[i] = A[m]
                A[m] = x
            end
        end
    end
end

################################################
#
# utilities
#
################################################

@inline level(i) = floor(Int, log2(i))
@inline lchild(i) = 2*i
@inline rchild(i) = 2*i+1
@inline children(i) = (lchild(i), rchild(i))
@inline hparent(i) = i ÷ 2
@inline on_minlevel(i) = level(i) % 2 == 0
@inline haschildren(i, A) = lchild(i) ≤ length(A)
@inline isgrandchild(j, i) = j > rchild(i)
@inline hasgrandparent(i) = i ≥ 4

"""
    children_and_grandchildren(maxlen, i)

Return the indices of all children and grandchildren of
position `i`.
"""
function children_and_grandchildren(maxlen::T, i::T) where {T <: Integer}
    left, right = children(i)
    _children_and_grandchildren = (left, children(left)..., right, children(right)...)
    return Iterators.filter(<=(maxlen), _children_and_grandchildren)
end

"""
    is_minmax_heap(h::AbstractVector) -> Bool

Return `true` if `A` is a min-max heap. A min-max heap is a
heap where the minimum element is the root and the maximum
element is a child of the root.
"""
function is_minmax_heap(A::AbstractVector)

    for i in 1:length(A)
        if on_minlevel(i)
            # check that A[i] < children A[i]
            #    and grandchildren A[i]
            for j in children_and_grandchildren(length(A), i)
                A[i] ≤ A[j] || return false
            end
        else
            # max layer
            for j in children_and_grandchildren(length(A), i)
                A[i] ≥ A[j] || return false
            end
        end
    end
    return true
end

################################################
#
# interfaces
#
################################################

Base.length(h::BinaryMinMaxHeap) = length(h.valtree)

Base.isempty(h::BinaryMinMaxHeap) = isempty(h.valtree)

"""
    pop!(h::BinaryMinMaxHeap) = popmin!(h)
"""
@inline Base.pop!(h::BinaryMinMaxHeap) = popmin!(h)

function Base.sizehint!(h::BinaryMinMaxHeap, s::Integer)
    sizehint!(h.valtree, s)
    return h
end

"""
    popmin!(h::BinaryMinMaxHeap) -> min

Remove the minimum value from the heap.
"""
function popmin!(h::BinaryMinMaxHeap)
    valtree = h.valtree
    !isempty(valtree) || throw(ArgumentError("heap must be non-empty"))
    @inbounds x = valtree[1]
    y = pop!(valtree)
    if !isempty(valtree)
        @inbounds valtree[1] = y
        @inbounds _minmax_heap_trickle_down!(valtree, 1)
    end
    return x
end

"""
    popmin!(h::BinaryMinMaxHeap, k::Integer) -> vals

Remove up to the `k` smallest values from the heap.
"""
@inline function popmin!(h::BinaryMinMaxHeap, k::Integer)
    return [popmin!(h) for _ in 1:min(length(h), k)]
end

"""
    popmax!(h::BinaryMinMaxHeap) -> max

Remove the maximum value from the heap.
"""
function popmax!(h::BinaryMinMaxHeap)
    valtree = h.valtree
    !isempty(valtree) || throw(ArgumentError("heap must be non-empty"))
    @inbounds x, i = maximum(((valtree[j], j) for j in 1:min(length(valtree), 3)))
    y = pop!(valtree)
    if !isempty(valtree) && i <= length(valtree)
        @inbounds valtree[i] = y
        @inbounds _minmax_heap_trickle_down!(valtree, i)
    end
    return x
end

"""
    popmax!(h::BinaryMinMaxHeap, k::Integer) -> vals

Remove up to the `k` largest values from the heap.
"""
@inline function popmax!(h::BinaryMinMaxHeap, k::Integer)
    return [popmax!(h) for _ in 1:min(length(h), k)]
end

function Base.push!(h::BinaryMinMaxHeap, v)
    valtree = h.valtree
    push!(valtree, v)
    @inbounds _minmax_heap_bubble_up!(valtree, length(valtree))
end

"""
    values(h::BinaryMinMaxHeap)

Returns the elements of the heap in an arbitrary order.
"""
Base.values(h::BinaryMinMaxHeap) = h.valtree

"""
    first(h::BinaryMinMaxHeap)

Get the first (minimum) of the heap.
"""
@inline Base.first(h::BinaryMinMaxHeap) = minimum(h)

@inline function Base.minimum(h::BinaryMinMaxHeap)
    valtree = h.valtree
    !isempty(h) || throw(ArgumentError("heap must be non-empty"))
    return @inbounds h.valtree[1]
end

@inline function Base.maximum(h::BinaryMinMaxHeap)
    valtree = h.valtree
    !isempty(h) || throw(ArgumentError("heap must be non-empty"))
    return @inbounds maximum(@views(valtree[1:min(end, 3)]))
end

Base.empty!(h::BinaryMinMaxHeap) = (empty!(h.valtree); h)

"""
    popall!(h::BinaryMinMaxHeap, ::Ordering = Forward)

Remove and return all the elements of `h` according to
the given ordering. Default is `Forward` (smallest to
largest).
"""
popall!(h::BinaryMinMaxHeap) = popall!(h, Forward)
popall!(h::BinaryMinMaxHeap, ::ForwardOrdering) = popmin!(h, length(h))
popall!(h::BinaryMinMaxHeap, ::ReverseOrdering) = popmax!(h, length(h))

ext: jl


text:
# text.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
module Text
using DataStructures, Devectorize, Ollam, Stage, GZip, Iterators
import Base: start, done, next, getindex, eltype

# -------------------------------------------------------------------------------------------------------------------------
# module-wide utilities
# -------------------------------------------------------------------------------------------------------------------------
getindex{K}(set :: Set{K}, key :: K) = contains(set, key) # sets don't have getindex method

include("constants.jl")
include("readers.jl")
include("tokenizers.jl")
include("features.jl")
include("models.jl")
include("tc.jl")

end


text:
module aero30
export coeffAero
function coeffAero(ma::Float64)
    if ma > 4.0
        cd0 = 0.647 * ma^(-0.5337)
    else
        cd0 = 0.0065 * ma^2 - 0.0965 * ma + 0.5901
    end
    return cd0
end
end


text:
#Pseudo-Random Bit Sequence Tests

@testset "PRBS Tests" begin #Scope for test data

_DEBUG = false

#Dumps contents of a bit sequence vector for debugging purposes:
function dbg_dumpseq(prefix, seq)
	if !_DEBUG; return; end
	print(prefix)
	for bit in seq
		if bit > 0
			print("1")
		else
			print("0")
		end
	end
	println()
end

@testset "Cyclical tests: MaxLFSR" begin
	for prbslen in 3:15 #Don't perform tests on higher-order patterns (too much time/memory)
		patlen = (2^prbslen) - 1
		genlen = patlen + prbslen #Full pattern length + "internal register length"
		regmask = (Int64(1)<<prbslen) -1 #Mask out "non-register" bits
		seed = rand(Int64) & regmask
		pattern = collect(sequence(MaxLFSR(prbslen), seed=seed, len=genlen, output=Int))

		#Make sure pattern repeats for at least as many bits as are in internal register:
		mismatch = pattern[(1:prbslen)] .!= pattern[patlen .+ (1:prbslen)]
		@test sum(convert(Vector{Int}, mismatch)) == 0

		#@show prbslen, pattern[1:prbslen], pattern[patlen .+ (1:prbslen)]
	end
end

@testset "Error detection: MaxLFSR" begin
	#TODO: Compare algorithm against known good pattern (if one can be found).
	prbslen = 15; seed = 11; len = 50
	pattern = collect(sequence(MaxLFSR(prbslen), seed=seed, len=len, output=Int))
	_errors = sequence_detecterrors(MaxLFSR(prbslen), pattern)

	dbg_dumpseq("_errors = ", _errors) #DEBUG
	@test sum(_errors) == 0

	#Test with error injected:
	pattern[prbslen+5] = 1 - pattern[prbslen+5] #Flip "bit"
	_errors = sequence_detecterrors(MaxLFSR(prbslen), pattern)

	dbg_dumpseq("_errors = ", _errors) #DEBUG
	@test sum(_errors) == 1 #Should have a single error
end

end #PRBS Tests

#Last line


text:
import Base: eltype, ==, issubset, zero, show

import LightGraphs: AbstractGraph, vertices, inneighbors, outneighbors,
add_vertex!, add_edge!, rem_vertex!, rem_edge!, edges, is_directed,
eltype, edgetype, nv, ne, has_vertex, has_edge, neighbors, src, dst

import LightGraphs.SimpleGraphs: SimpleDiGraph, SimpleEdge, Edge, fadj, badj

export BVHGraph, vertices, inneighbors, outneighbors, add_vertex!, add_edge!,
rem_vertex!, rem_edge!, props, edges, is_directed, eltype, edgetype, nv, ne,
has_vertex, has_edge, ==, zero, fadj, badj, issubset, neighbors,
name, name!, frames, frames!, frametime, frametime!, offset, offset!,
sequence, sequence!, rotations, rotations!, positions, positions!,
find, find_outneighbor, show, src, dst

Base.@kwdef mutable struct GProps
    name::AbstractString = "BVHGraph.bvh"
    offset::Vector{Float64} = zeros(Float64, 3)
    sequence::Symbol = :XYZ
    positions::Matrix{Float64} = zeros(Float64, 1, 3)
    frames::Int64 = 1
    frametime::Float64 = 0.01
end

Base.@kwdef mutable struct VProps
    name::AbstractString = "End Site"
    sequence::Symbol = :XYZ
    rotations::Matrix{Float64} = zeros(Float64, 1, 3)
    positions::Matrix{Float64} = zeros(Float64, 1, 3)
end

Base.@kwdef mutable struct EProps
    offset::Vector{Float64} = zeros(Float64, 3)
end

"""
    struct BVHGraph{T <: Integer} <: AbstractGraph{T}

Used to store motion capture data. Usually constructed by using `load` on a BVH file.
"""
struct BVHGraph{T <: Integer} <: AbstractGraph{T}
    graph::SimpleDiGraph{T}
    gprops::GProps
    vprops::Dict{T,VProps}
    eprops::Dict{SimpleEdge{T},EProps}
end

function BVHGraph(n::Integer;
    name::AbstractString = "BVHGraph.bvh",
    offset::Vector{Float64} = zeros(Float64, 3),
    sequence::Symbol = :XYZ,
    positions::Matrix{Float64} = zeros(Float64, 1, 3),
    frames::Int64 = 1,
    frametime::Float64 = 0.01)

    T = eltype(n)
    g = SimpleDiGraph(n)
    gprops = GProps(name = name,
                    offset = offset,
                    sequence = sequence,
                    positions = positions,
                    frames = frames,
                    frametime = frametime)
    vprops = Dict{T,VProps}()
    eprops = Dict{SimpleEdge{T},EProps}()

    for v in vertices(g)
        vprops[v] = VProps()
    end

    return BVHGraph(g, gprops, vprops, eprops)
end

SimpleDiGraph(g::BVHGraph) = g.graph

vertices(g::BVHGraph) = vertices(g.graph)

inneighbors(g::BVHGraph, v::Integer) = inneighbors(g.graph, v)
outneighbors(g::BVHGraph, v::Integer) = fadj(g.graph, v)

function add_vertex!(g::BVHGraph;
    name::AbstractString = "End Site",
    sequence::Symbol = :XYZ,
    rotations::Matrix{Float64} = zeros(Float64, 1, 3),
    positions::Matrix{Float64} = zeros(Float64, 1, 3))

    add_vertex!(g.graph) || return 0
    v = nv(g)
    g.vprops[v] = VProps(name, sequence, rotations, positions)
    return v
end

function add_edge!(g::BVHGraph, s::Integer, d::Integer; offset::Vector{Float64} = zeros(Float64, 3))
    add_edge!(g.graph, s, d) || return false
    e = Edge(s, d)
    g.eprops[e] = EProps(offset = offset)
    return true
end

function rem_vertex!(g::BVHGraph, v::Integer)
    v in vertices(g) || return false

    lastv = nv(g)
    lastvprops = props(g, lastv)
    lasteoutprops = Dict(n => props(g, lastv, n) for n in outneighbors(g, lastv))
    lasteinprops = Dict(n => props(g, n, lastv) for n in inneighbors(g, lastv))

    for n in outneighbors(g, lastv)
        delete!(g.eprops, Edge(lastv, n))
    end

    for n in inneighbors(g, lastv)
        delete!(g.eprops, Edge(n, lastv))
    end

    if v != lastv

        for n in outneighbors(g, v)
            delete!(g.eprops, Edge(v, n))
        end

        for n in inneighbors(g, v)
            delete!(g.eprops, Edge(n, v))
        end

    end

    delete!(g.vprops, lastv)
    retval = rem_vertex!(g.graph, v)
    retval || return false

    if v != lastv
        g.vprops[v] = lastvprops

        for n in outneighbors(g, v)
            g.eprops[Edge(v, n)] = lasteoutprops[n]
        end

        for n in inneighbors(g, v)
            g.eprops[Edge(n, v)] = lasteinprops[n]
        end
    end

    return true
end

function rem_edge!(g::BVHGraph, s::Integer, d::Integer)
    delete!(g.eprops, Edge(s, d))
    rem_edge!(g.graph, s, d)
end

props(g::BVHGraph) = getproperty(g, :gprops)
props(g::BVHGraph, v::Integer) = get(g.vprops, v, VProps())
props(g::BVHGraph, s::Integer, d::Integer) = get(g.eprops, Edge(s, d), EProps())

edges(g::BVHGraph) = edges(g.graph)

is_directed(::Type{BVHGraph}) = true
is_directed(::Type{BVHGraph{T}}) where T = true
is_directed(g::BVHGraph) = true

eltype(g::BVHGraph) = eltype(g.graph)
edgetype(g::BVHGraph) = edgetype(g.graph)

nv(g::BVHGraph) = nv(g.graph)
ne(g::BVHGraph) = ne(g.graph)

has_vertex(g::BVHGraph, x...) = has_vertex(g.graph, x...)
@inline has_edge(g::BVHGraph, x...) = has_edge(g.graph, x...)

==(g::BVHGraph, h::BVHGraph) = g.graph == h.graph

zero(::BVHGraph) = BVHGraph(0)
zero(g::BVHGraph) = BVHGraph(0, name = g.name)

@inline fadj(g::BVHGraph, x...) = fadj(g.graph, x...)
@inline badj(g::BVHGraph, x...) = badj(g.graph, x...)

issubset(g::T, h::T) where T <: BVHGraph = issubset(g.graph, h.graph)

neighbors(g::BVHGraph, v::Integer) = outneighbors(g, v)

name(g::BVHGraph) = getproperty(props(g), :name)
name(g::BVHGraph, v::Integer) = getproperty(props(g, v), :name)
name(v::Integer) = g -> name(g, v)

name!(g::BVHGraph, name::AbstractString) = setproperty!(props(g), :name, name)
name!(name::AbstractString) = g -> name!(g, name)
name!(g::BVHGraph, v::Integer, name::AbstractString) = setproperty!(props(g, v), :name, name)
name!(v::Integer, name::AbstractString) = g -> name!(g, v, name)

frames(g::BVHGraph) = getproperty(props(g), :frames)

frames!(g::BVHGraph, f::Int64) = setproperty!(props(g), :frames, f)
frames!(f::Int64) = frames!(g, f)
# frames(g::BVHGraph) = 1:getproperty(props(g), :frames)

frametime(g::BVHGraph) = getproperty(props(g), :frametime)

frametime!(g::BVHGraph, t::Float64) = setproperty!(props(g), :frametime, t)
frametime!(t::Float64) = g -> frametime!(g, t)

offset(g::BVHGraph) = getproperty(props(g), :offset)
offset(g::BVHGraph, s::Integer, d::Integer) = getproperty(props(g, s, d), :offset)
offset(g::BVHGraph, d::Integer) = offset(g, inneighbors(g, d)[1], d)

offset!(g::BVHGraph, o::Vector{Float64}) = setproperty!(props(g), :offset, o)
offset!(g::BVHGraph, s::Integer, d::Integer, o::Vector{Float64}) = setproperty!(props(g, s, d), :offset, o)
offset!(g::BVHGraph, d::Integer, o::Vector{Float64}) = offset!(g, inneighbors(g, d)[1], d, o)

sequence(g::BVHGraph) = getproperty(props(g), :sequence)
sequence(g::BVHGraph, v::Integer) = getproperty(props(g, v), :sequence)
sequence(v::Integer) = g -> sequence(g, v)

sequence!(g::BVHGraph, sym::Symbol) = setproperty!(props(g), :sequence, sym)
sequence!(sym::Symbol) = g -> sequence!(g, sym)
sequence!(g::BVHGraph, v::Integer, sym::Symbol) = setproperty!(props(g, v), :sequence, sym)
sequence!(v::Integer, sym::Symbol) = g -> sequence!(g, v, sym)

rotations(g::BVHGraph, v::Integer) = getproperty(props(g, v), :rotations)
rotations(g::BVHGraph, v::Integer, f::Integer) = rotations(g, v)[f, :]

rotations!(g::BVHGraph, v::Integer, matrix::Matrix{Float64}) = setproperty!(props(g, v), :rotations, matrix)

positions(g::BVHGraph) = getproperty(props(g), :positions)
positions(g::BVHGraph, v::Integer) = getproperty(props(g, v), :positions)
positions(g::BVHGraph, v::Integer, f::Integer) = positions(g, v)[f, :]

positions!(g::BVHGraph, matrix::Matrix{Float64}) = setproperty!(props(g), :positions, matrix)
positions!(g::BVHGraph, v::Integer, matrix::Matrix{Float64}) = setproperty!(props(g, v), :positions, matrix)

function find(g::BVHGraph, vertexname::AbstractString)
    for v in vertices(g)
        name(g, v) |> endswith(vertexname) && return v
    end

    return 0
end

function find_outneighbor(g::BVHGraph, v, vertexname::AbstractString)
    endswith(vertexname, "End Site") && return outneighbors(g, v)[1]

    for n in outneighbors(g, v)
        name(g, n) |> endswith(vertexname) && return n
    end

    return 0
end

function show(io::IO, g::BVHGraph)
    println(io, "BVHGraph")
    println(io, "Name: $(name(g))")
    println(io, "Frames: $(frames(g))")
    println(io, "Frame Time: $(frametime(g))")
    println(io, "Vertices: $(nv(g))\n")
    println(io, get_names_formatted(g))
    return nothing
end

function get_names_formatted(g::BVHGraph, v::Integer = 1, s::AbstractString = "", spaces::Integer = 0)
    s *= "$(' '^spaces)$(v) $(name(g, v))\n"

    for n in outneighbors(g, v)
        s = get_names_formatted(g, n, s, spaces + 1)
    end

    return s
end


text:
#
# Protein - TMAO (compare new and old implementations)
#

include("../../src/MDDF.jl")
using PDBTools
using Plots

dir="/home/leandro/Drive/Alunos/Ivan/PCCP_revision"

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = PDBTools.readPDB("$dir/6Mnative.pdb")

# Set solute and solvent

protein = PDBTools.select(atoms,"protein")
solute = MDDF.Selection(protein,nmols=1 )

water = PDBTools.select(atoms,"water")
solvent = MDDF.Selection(water,natomspermol=3)

options = MDDF.Options()

# Initialize trajectroy data structure and open input stream
#R = MDDF.mddf(trajectory,options)
#MDDF.save(R,"ivan.json")

trajectory = MDDF.NamdDCD("$dir/6Mnative.dcd",solute,solvent)
#lc = MDDF.mddf_linkedcells(trajectory,options)

#trajectory = MDDF.Trajectory("$dir/6Mnative.dcd",solute,solvent)
lcP = MDDF.mddf(trajectory,options)
#lcP = MDDF.load("./ivan.json")

plot(layout=(5,1))

sp=1
plot!(ylabel="mddf",sp=sp)
plot!(lc.d,lc.mddf,linewidth=2,label="lc",subplot=sp)
scatter!(lcP.d,lcP.mddf,linewidth=2,label="lcP",subplot=sp)

sp=2
plot!(ylabel="md_count",sp=sp)
plot!(lc.d,lc.md_count,linewidth=2,label="lc",subplot=sp)
scatter!(lcP.d,lcP.md_count,linewidth=2,label="lcP",subplot=sp)

sp=3
plot!(ylabel="md_count_random",sp=sp)
plot!(lc.d,lc.md_count_random,linewidth=2,label="lc",subplot=sp)
scatter!(lcP.d,lcP.md_count_random,linewidth=2,label="lcP",subplot=sp)

sp=4
plot!(ylabel="rdf_count",sp=sp)
plot!(lc.d,lc.rdf_count,linewidth=2,label="lc",subplot=sp)
scatter!(lcP.d,lcP.rdf_count,linewidth=2,label="lcP",subplot=sp)

sp=5
plot!(ylabel="rdf_count_random",sp=sp)
plot!(lc.d,lc.rdf_count_random,linewidth=2,label="lc",subplot=sp)
scatter!(lcP.d,lcP.rdf_count_random,linewidth=2,label="lcP",subplot=sp)

plot!(size=(400,800))
savefig("./ivan.pdf")
savefig("./plots.pdf")


text:
using JuliaForHEP
using JuliaForHEP: step!
using Flux
using Flux: Data.DataLoader
using DataFrames
using Random
using LinearAlgebra
using Compat
using PyPlot

Random.seed!(123)

include("variables.jl")
include("model.jl")
include("training.jl")
include("save_model.jl")

#basedir = "/work/sschaub/h5files/ttH-processed/"
basedir = "/local/scratch/ssd/sschaub/h5files/ttH_sebastian-processed/"
h5files = joinpath.(basedir, ["ttH", "ttbb", "ttcc", "ttlf"] .* "_lola.h5")

for i in 1:10

for (scalar_features, include_jets) in zip([
    [String[]]; #=Base.vect.(Vars.scalar_features); =#
    [Vars.tier3, [Vars.tier2; Vars.tier3], Vars.scalar_features];
    [Vars.scalar_features, Vars.scalar_features]
],
    [fill(:lola, 4); :jets_as_scalars; :none]
)
    @show scalar_features, include_jets
    feature_names = if include_jets === :none
        scalar_features
    else
        [Vars.jets; Vars.tight_lepton; Vars.jet_csv; scalar_features]
    end
    inputs = extract_cols(h5files, feature_names) .|> Float32

    classes = Vars.classes
    outputs = Symbol.(extract_cols(h5files, ["class_label"]))
    outputs_onehot = Flux.onehotbatch(vec(outputs), classes)

    weights = prod(Float32, extract_cols(h5files, Vars.weights); dims=1)
    total_weights = scale_weights(weights, outputs)

    ds_train, ds_validation, ds_test = split_datasets(
        gpu,
        (inputs, outputs_onehot, total_weights, axes(inputs, 2)),
        [.75 * .8, .75 * .2, .25],
    )
    idx_train, idx_validation, idx_test = last.((ds_train, ds_validation, ds_test))
    ds_train = DataLoader(ds_train; batchsize=1024, shuffle=true)
    ds_validation = DataLoader(ds_validation; batchsize=1024, shuffle=false)
    ds_test = DataLoader(ds_test; batchsize=1024, shuffle=false)

    m = include_jets === :lola ? 15 : 0
    n_jets = include_jets === :lola ? length([Vars.jets; Vars.tight_lepton]) ÷ 4 : 0
    n_inputs = size(inputs, 1)
    n_outputs = length(classes)
    model = build_model(;
        m, n_jets, n_inputs, n_outputs,
        w_d_reducers = ((+) => 2, min => 2),
        neurons = [1024, 2048, 512, 512],
        activation = _leakyrelu,
        dropout = Dropout(0.5),
    ) |> gpu

    optimizer = ADAM(5e-5)

    penalty = l2_penalty(model)
    loss(ŷ, y; weights) = Flux.logitcrossentropy(ŷ, y; agg = x -> weighted_mean(x, weights)) + 1f-5 * penalty()

    measures = (; loss=st->st.loss, accuracy=st->accuracy(softmax(st.ŷ), st.y))

    recorded_measures = DataFrame()

    if include_jets === :lola
        # pretrain Lo+CoLa
        n_outputs_lola = (n_jets + m) * (3 + 4) + length(Vars.jet_csv)
        _model = Chain(
            model[1:end-1]...,
            #x -> view(x, 1:n_outputs_lola, :),
            x -> x[1:n_outputs_lola, :],
            build_dnn(;
                n_inputs = n_outputs_lola,
                neurons = [1024, 2048, 512, 512],
                n_outputs,
                activation = _leakyrelu,
                dropout = Dropout(0.5),
            ) |> gpu,
        )
        train!(;
            model=_model, ds_train, ds_test, ds_validation,
            loss, optimizer, measures, recorded_measures,
            max_epochs=20, min_epochs=50, early_stopping_n=10, early_stopping_percentage=2,
        )
    else
        model = model[2:end]
    end

    # main training
    train!(;
        model, ds_train, ds_test, ds_validation,
        loss, optimizer, measures, recorded_measures,
        max_epochs=120, min_epochs=50, early_stopping_n=7, early_stopping_percentage=2,
    )

    name = string(include_jets) * "+" * if isempty(scalar_features)
        "none"
    elseif length(scalar_features) == 1
        scalar_features[]
    elseif scalar_features == Vars.tier3
        "tier3"
    elseif scalar_features == [Vars.tier2; Vars.tier3]
        "tier2+3"
    elseif scalar_features == Vars.scalar_features
        "all"
    end

    save_model(;
        output_dir = "/work/sschaub/JuliaForHEP/feature_evaluation2_0131/$(name)_$i/",
        fig = nothing,
        recorded_measures,
        model,
        inputs,
        feature_names,
        n_jets,
        classes,
        outputs,
        weights,
        total_weights,
        idx_train, idx_test, idx_validation,
    )
end

end


text:
module ENA
using ManualScraper: Parser

const SYM = :ena

# Directories to pull from inside ENA's help file.
const DIRS = map(x->"/Programming/Command_Reference/"*x,
    ["ABORT", "CALCULATE", "CONTROL", "DISPLAY", "FORMAT", "HCOPY", #"IEEE",
    "INIT", "LXI", "MEMORY", "OUTPUT", "PROGRAM", "SENSE",
    "SERVICE", "SOURCE", "STATUS", "SYSTEM", "TRIGGER"])

# Files to pull from inside ENA's help file.
const FILTER = r".htm$"

# This is necessary due to inconsistencies in the HTML files.
const INIT = [

    # Remove control chars and extended chars except CR and LF.
    (r"[^\x20-\x7e\x0a\x0d]", ""),

    # Remove certain tags that randomly interrupt command syntax
    ("
", ""),
    ("
", ""),
    ("",""),
    (r"",""),
    (r"", ""),
    (r"", ""),
    ("",""),
    ("",""),

    # Remove class attributes
    (r" style=\"[A-Za-z0-9\-\s;:]+\"", ""),
    (r" class=\"[a-z0-9]+\"", ""),

    # Contract extra white space
    (" - ", "-"),
    (" ", ""),
    (r"\s+\|", "|"),
    (r"\|\s+", "|"),

    # Fix some syntactical typos. There are a lot, thanks Keysight
    (r"COLLect\s+:ITERation", "COLLect:ITERation"),
    ("CALCulate{[1]-160}:TRACe(Tr):", "CALCulate{[1]-160}:TRACe{[1]-16}:"),
    (":CALCulate{[1]-160]", ":CALCulate{[1]-160}"),
    (" [:", "[:"),
    (" :",":"),
    (" ?", "?"),
    ("Aux1", "AUX1"),
    ("E/span>", "..."),
    (r": ([A-Z]+)", s":\1"),
    (r"([A-Z]).VOLTage.", s"\1:VOLTage"),
    ("STORe.STSource", "STORe:STSource"),
    ("DC.RANGe", "DC:RANGe"),
    ("[:STATe]{", " {"),
    ("SENSe([1]-160)", "SENSe{[1]-160}"),
    ("SENSe[1-160]", "SENSe{[1]-160}"),
    ("PORT([1]-4)", "PORT{[1]-4}"),
    ("SOLT4 1,2,3,4
","SOLT4 <numeric 1>,<numeric 2>,<numeric 3>,<numeric 4>
"),
    ("RECeiver([1]-4)", "RECeiver{[1]-4}"),
    (":[EXTernal]", "[:EXTernal]"),
    ("CALCulate[1-160]", "CALCulate{[1]-160}"),
    (r"STARt\|\s+STOP", "STARt|STOP"),
    (r"R<\/p>\s+
1\|R2\|R3", "R1|R2|R3"),
    (r"(EXTernal.SLOPe<\/h1>[\s]+
Object type<\/h2>[\s]+
Property)(
)", s"\1 (Read-Write)\2"),
    (r"(REPort\.DATA[<>\/A-Za-z0-9\s=\"\-:';.&#]+
Property)(<\/p>)", s"\1 (Read-Only)\2"),
    (r"(<string 1>)(<\^END>)", s"\1<newline>\2"),
    (r"<1/0>(<\^END>)", s"{1|0}<newline>\1"),
    ("{1/0}","{1|0}"),
    (r"<string \("n, m, l\.\.\."\)>\s*", "<string>
\r\n
"),
    (r"(<string>)\s+(:CALCulate\{\[1\]\-160\})", s"\1
\r\n
\2"),
    # ("<numeric>", "{numeric}"),
    # ("<value>", "{value}"),

    # Make parsing more convenient in a few special cases
    ("Data Array", "DataArray"),
    (r";([A-Za-z]+)\s+([0-9])&", s";\1\2&"),
    (""Sx1y1, Sx2y2, ... , Sxnxn"", "<string>"),
    ("<+ or -><value of ", "<"),
    (r"SCPI commands:<\/p>[<>a-zA-Z0-9.;\"\{\}\-?:=\s\/]+
", "SCPI commands:

"), # in :FORMAT:DATA
    (r"(CATalog\?<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Read-Only)\2"),
    (r"(:NAME<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Read-Write)\2"),
    (r"(STATe<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Read-Write)\2"),
    (r"(TRANsfer<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Read-Write)\2"),
    (r"(LOAD:PROGram<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Write-Only)\2"),
    (r"(STORe:PROGram<\/h1>\s+
)No equivalent COM Commands(<\/p>)", s"\1Property (Write-Only)\2"),

    # (r"{numeric 1},,{numeric NOP.2.\(number\s+of specified traces\)}", "{array}"),
    # ("<numeric1>,,<numeric N><^END>", "{array}<newline><^END>"),
    # ("{numeric 1},,{numeric NOP?2}", "{array}"),
    # ("<numeric1>,,<numeric NOP?2>", "{array}<newline>"),
    # ("<numeric1>,,<numeric N*2>", "{array}<newline>"),
    # (r"\{numeric 1\},\s*\.\.\.\s*,\{numeric 1\+\(N\*5\)\}","{array}"),
    # ("{numeric 1},,{numeric 1+(Nx4)}","{array}"),
    # ("{numeric 1},,{numeric 1+Nx3}","{array}"),
    # ("{numeric 1},{numeric 2},...,{numeric N}","{array}"),
    # ("{numeric 1}, ...,{numeric NOP*2}", "{array}"),
    # ("{numeric 1},...,{numeric NOP}","{array}"),
    # ("{value 1},...,{value NOP}","{array}"),
    # ("R1|R1|R2|R3}", ""),

]

# Determine if the commands are for read/write, read, or write
const GETSET_REGEX = r"((?:property)|(?:method)) \((?:(read|write|only)?[\-\s\/]*(read|write|only)?)\s?\)"i

# How to parse the ENA's HTML files for commands and args.
# l1. Preceding line must end with >\r\n
# l2. Start with
: then follow with any # of chars that may be found in a SCPI command
# l3. (variable whitespace)
# l4. One of two things follows:
#   1) Some text between < and &rt; or
#   2) Some text in curly braces
# l5. optionally have a comma and/or whitespace, then repeat to l4
# l6. optionally have whitespace, terminate with a
 tag and optional whitespace.
# const CMD_REGEX = r"(?<=>\r\n)(?:
(:[:A-Za-z\(\)\-\{\}\[\]\|?0-9\\\n\r]+)\s*(?:(?:(<[A-Z\(\)\s,a-z0-9]+>)|\{([A-Za-z0-9\-\|]+)\}),*\s*)*\s*<\/p>\s*)"
const CMD_REGEX = r"(?<=>\r\n)(?:
(:[:A-Za-z\(\)\-\{\}\[\]\|?0-9\\\n\r]+)\s*((?:(?:(?:<[A-Z\(\)\s,a-z0-9]+>)|\{(?:[A-Za-z0-9\-\|]+)\}),*\s*)+)*\s*<\/p>\s*)"

const CMD_PASS1 = [
    # Strip junk: stray HTML tag snippets
    (r"<[A-Za-z0-9]+>?", s""),
    (r"[<>]+",           s""),

    # Strip newline and form feeds
    (r"[\r\n]+", ""),

    # Strip optional segments of command
    (r"\[:[:A-Za-z\|]+\]", s"")
]

const CMD_PASS2 = [
    # Strip lowercase parts
    (r"[a-z]+",                    ""),

    # Replace infixes. Could be automated...
    (r"(ADAP)\{[\[\]0-9\-\|]+\}",  s"\1pt"),
    (r"(ARR)\{[\[\]0-9\-\|]+\}",   s"\1vnum"),
    (r"(BPOR)\{[\[\]0-9\-\|]+\}",  s"\1bp"),
    (r"(:C)\{[\[\]0-9\-\|]+\}",    s"\1el"),
    (r"(CALC)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(CHAN)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(CKIT)\{[\[\]0-9\-\|]+\}",  s"\1ckit"),
    (r"(COL)\{[\[\]0-9\-\|]+\}",   s"\1dnum"),
    (r"(DOUB)\{[\[\]0-9\-\|]+\}",  s"\1vnum"),
    (r"(FREQ)\{[\[\]0-9\-\|]+\}",  s"\1freq"),
    (r"(:G)\{[\[\]0-9\-\|]+\}",    s"\1el"),
    (r"(GRAT)\{[\[\]0-9\-\|]+\}",  s"\1gnum"),
    (r"(INCL)\{[\[\]0-9\-\|]+\}",  s"\1incl"),
    (r"(INIT)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(:L)\{[\[\]0-9\-\|]+\}",    s"\1el"),
    (r"(LIM)\{[\[\]0-9\-\|]+\}",   s"\1lnum"),
    (r"(LONG)\{[\[\]0-9\-\|]+\}",  s"\1vnum"),
    (r"(LOSS)\{[\[\]0-9\-\|]+\}",  s"\1loss"),
    (r"(MARK)\{[\[\]0-9\-\|]+\}",  s"\1m"),
    (r"(MULT)\{[\[\]0-9\-\|]+\}",  s"\1mult"),
    (r"(NETW)\{[\[\]0-9\-\|]+\}",  s"\1netw"),
    (r"(OUTP)\{[\[\]0-9\-\|]+\}",  s"\1num"),
    (r"(PAR)\{[\[\]0-9\-\|]+\}",   s"\1par"),
    (r"(PLOS)\{[\[\]0-9\-\|]+\}",  s"\1pt"),
    (r"(PORT)\{[\[\]0-9\-\|]+\}",  s"\1pt"),
    (r"(:R)\{[\[\]0-9\-\|]+\}",    s"\1el"),
    (r"(REC)\{[\[\]0-9\-\|]+\}",   s"\1pt"),
    (r"(SENS)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(STAN)\{[\[\]0-9\-\|]+\}",  s"\1std"),
    (r"(STR)\{[\[\]0-9\-\|]+\}",   s"\1vnum"),
    (r"(SOUR)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(TYPE)\{[\[\]0-9\-\|]+\}",  s"\1tr"),
    (r"(TRAC)\{[\[\]0-9\-\|]+\}",  s"\1tr"),
    (r"(WIND)\{[\[\]0-9\-\|]+\}",  s"\1ch"),
    (r"(VOLT)\{[\[\]0-9\-\|]+\}",  s"\1auxpt")
]

const ARG_REGEX = r""

const ARG_PASS = [
    (r"[\{\}\s]+", ""),
    ("ON|OFF|1|0", "v::Bool"),
    (r"<numeric[0-9]?\(?[A-Za-z]*\)?>", "v::Real"),
    (r"<value[0-9]?\(?[A-Za-z]*\)?>", "v::Real"),
    (r"<string[0-9]?>", "v::AbstractString"),
]

const RETTYPE_REGEX = r""
# const RETTYPE_REGEX = r"[Qq]uery [Rr]esponse<\/h[0-9]>\s+
(?:\{([A-Za-z0-9,_\-\?\|…\.\(\)\/\s]+)\},*\s*)*<newline"

const PARSER = Parser(INIT, GETSET_REGEX, CMD_REGEX, RETTYPE_REGEX,
    CMD_PASS1, CMD_PASS2, ARG_REGEX, ARG_PASS)

end


text:
module SeisNetwork

using SeisIO, SeisNoise, Printf, Triangle, LinearAlgebra, Statistics, Geodesy, JLD2, PyPlot, Plots
using LightXML, Dates, HDF5, Distributed

# import pre and post processing tools
include("triangulation.jl")
include("stationkml2jld.jl")
include("spatial_dvv.jl")
include("average_dvv.jl")

end # module


text:
#md# ### First generate the wordcloud on the left
using WordCloud

stwords = ["us"];
println("==Obama's==")
cs = WordCloud.randomscheme() #:Set1_8
as = WordCloud.randomangles() #(0,90,45,-45)
dens = 0.5 #not too high
wca = wordcloud(
    processtext(open(pkgdir(WordCloud)*"/res/Barack Obama's First Inaugural Address.txt"), stopwords=WordCloud.stopwords_en ∪ stwords),
    colors = cs,
    angles = as,
    density = dens,
    backgroundcolor = :maskcolor,
    ) |> generate!
#md# ### Then generate the wordcloud on the right
println("==Trump's==")
wcb = wordcloud(
    processtext(open(pkgdir(WordCloud)*"/res/Donald Trump's Inaugural Address.txt"), stopwords=WordCloud.stopwords_en ∪ stwords),
    mask = getsvgmask(wca),
    colors = cs,
    angles = as,
    density = dens,
    backgroundcolor = :maskcolor,
    maskcolor = getmaskcolor(wca),
    font = getparameter(wca, :font),
    run = x->nothing, #turn off the useless initword! and placewords! in advance
)
#md# Follow these steps to generate a wordcloud: initword! -> placewords! -> generate!
samewords = getwords(wca) ∩ getwords(wcb)
println(length(samewords), " same words")

for w in samewords
    setcolors!(wcb, w, getcolors(wca, w))
    setangles!(wcb, w, getangles(wca, w))
end
initwords!(wcb)

println("=ignore defferent words=")
keep(wcb, samewords) do
    @assert Set(wcb.words) == Set(samewords)
    centers = getpositions(wca, samewords, type=getcenter)
    setpositions!(wcb, samewords, centers, type=setcenter!) #manually initialize the position,
    setstate!(wcb, :placewords!) #and set the state flag
    generate!(wcb, 1000, teleporting=false, retry=1) #turn off the teleporting; retry=1 means no rescale
end

println("=pin same words=")
pin(wcb, samewords) do
    placewords!(wcb)
    generate!(wcb, 1000, retry=1) #allow teleport but don‘t allow rescale
end

if getstate(wcb) != :generate!
    println("=overall tuning=")
    generate!(wcb, 1000, teleporting=setdiff(getwords(wcb), samewords), retry=2) #only teleport the unique words
end

ma = paint(wca)
mb = paint(wcb)
h,w = size(ma)
space = fill(mb[1], (h, w÷20))
try mkdir("address_compare") catch end
println("results are saved in address_compare")
WordCloud.save("address_compare/compare.png", [ma space mb])
#eval# try rm("address_compare", force=true, recursive=true) catch end
gif = WordCloud.GIF("address_compare")
record(wca, "Obama", gif)
record(wcb, "Trump", gif)
WordCloud.Render.generate(gif, framerate=1)
wca, wcb
#eval# runexample(:compare)
#md# ![](address_compare/compare.png)
#md# ![](address_compare/result.gif)


text:
"""

以下のプログラムを Julia に翻訳してみる。

三角多項式グラフ
http://aoki2.si.gunma-u.ac.jp/R/Andrews.html

使用法

using RDatasets
iris = dataset("datasets", "iris")
y = Matrix(iris[:, 1:4])
andrewsgraph(y)

col = transpose(vcat(fill(1, 50), fill(2, 50), fill(3, 50)))
andrewsgraph(y, normalize=false, col=col)
"""

using Statistics, Plots

function andrewsgraph(dat; normalize=true, points=100, col=1, alpha=0.2)
    data = copy(dat)
    n, nv = size(data)
    if normalize
        means = mean(data, dims=1)
        stds = std(data, dims=1)
        for i = 1:nv
            data[:, i] = (data[:, i] .- means[i]) ./ stds[i]
        end
    end
    θ = range(-π, π, length=points)
    coef = zeros((nv, points))
    for i = 2:nv
        func = i % 2 == 0 ? sin : cos
        coef[i, :] = func.((i ÷ 2) .* θ)
    end
    coef[1, :] = repeat([1/sqrt(2)], points)
    data = transpose(data * coef)
    pyplot()
    plot(θ, data, color=col, alpha=alpha, grid=false,
            yticks=false, yshowaxis=false, label="")
    xticks!(range(-π, π, length=5), ["-π", "-π/2", "0", "π/2", "π"],
            tick_direction=:out)
end;

max_stars_repo_name: r-de-r/stats

max_issues_repo_name: r-de-r/stats

max_forks_repo_name: r-de-r/stats


text:
using Flux, Test
using Flux.CUDA
using Flux: cpu, gpu
using Statistics: mean
using LinearAlgebra: I, cholesky, Cholesky

@testset "CUDA" begin
  x = randn(5, 5)
  cx = gpu(x)
  @test cx isa CuArray

  @test Flux.onecold(gpu([1.0, 2.0, 3.0])) == 3

  x = Flux.onehotbatch([1, 2, 3], 1:3)
  cx = gpu(x)
  @test cx isa Flux.OneHotMatrix && cx.data isa CuArray
  @test (cx .+ 1) isa CuArray

  m = Chain(Dense(10, 5, tanh), Dense(5, 2), softmax)
  cm = gpu(m)

  @test all(p isa CuArray for p in params(cm))
  @test cm(gpu(rand(10, 10))) isa CuArray{Float32,2}

  xs = rand(5, 5)
  ys = Flux.onehotbatch(1:5,1:5)
  @test collect(cu(xs) .+ cu(ys)) ≈ collect(xs .+ ys)

  c = gpu(Conv((2,2),3=>4))
  x = gpu(rand(10, 10, 3, 2))
  l = c(gpu(rand(10,10,3,2)))
  @test gradient(x -> sum(c(x)), x)[1] isa CuArray

  c = gpu(CrossCor((2,2),3=>4))
  x = gpu(rand(10, 10, 3, 2))
  l = c(gpu(rand(10,10,3,2)))
  @test gradient(x -> sum(c(x)), x)[1] isa CuArray

end

@testset "onecold gpu" begin
  y = Flux.onehotbatch(ones(3), 1:10) |> gpu;
  @test Flux.onecold(y) isa CuArray
  @test y[3,:] isa CuArray
end

@testset "restructure gpu" begin
  dudt = Dense(1,1) |> gpu
  p,re = Flux.destructure(dudt)
  foo(x) = sum(re(p)(x))
  @test gradient(foo, cu(rand(1)))[1] isa CuArray
end

@testset "GPU functors" begin
  @testset "Cholesky" begin
    M = 2.0*I(10) |> collect
    Q = cholesky(M)
    Q_gpu = Q |> gpu
    @test Q_gpu isa Cholesky{<:Any,<:CuArray}
    Q_cpu = Q_gpu |> cpu
    @test Q_cpu == cholesky(eltype(Q_gpu).(M))
  end
end


text:
# Functions assisting those in run_base.jl
# get Marker object from BinaryMarker object, at start and finish indices (m: markers, i:indivs)
function get_bed_values(markers::BinMarkers, m1, m2)
    genotypes = bedvals[view(markers.binarygenotypes, :, m1:m2)]
    flip_alleles!(genotypes)
    info = markers.markerInfo[m1:m2,:]
    return Markers(info, markers.IID, Array{Float64,2}(genotypes))
end

function calc_betas(marker, y::Array{Float64,2}, n::Int64, rinv)
    m = Array{Int64,1}(marker .+ 1)
    X = sparse(1:n, m, 1.0)
    rinvX = rinv .* X
    cellcounts = sum(rinvX, dims=1)
    if cells_zero(cellcounts')
        return zeros(3,1), zeros(3,1), true
    end
    cellsums = rinvX'y
    return cellsums, cellcounts', false
end

function select_snp(submodel, i)
    return Markers(submodel.snpset.markerInfo[i:i,:], submodel.snpset.IID, submodel.snpset.genotypes[:,i:i])
end

function get_1loci_k()
    return Array{Float64,2}([1.0 0.0 0.0
    -0.5 0.0 0.5
    -0.5 1.0 -0.5])
end

function cells_zero(cells)
    return (size(cells,1) < 3 || cells[1] <= 5 || cells[2] <= 5 || cells[3] <= 5)
end

function next_phenotype_interval!(model)
    model.layerlimit += 0.5*model.losoInterval
    println("New Limit: ", model.layerlimit)
    if model.layerlimit < model.maxlayerlimit
        if model.inlayer1
            model.nlayer2 += 1
            model.phenotypesfilename = string(model.phenotypesprefix, ".layer2.n", model.nlayer2, ".tsv")
            model.phenotypes = read_yk(model.phenotypesfilename)
        else
            model.nlayer1 += 1
            model.phenotypesfilename = string(model.phenotypesprefix, ".layer1.n", model.nlayer1, ".tsv")
            model.phenotypes = read_yk(model.phenotypesfilename)
        end
        model.inlayer1 = !model.inlayer1
        get_phenotype_intersect!(model.phenotypes, model.phenoOrder)
        get_phenotype_intersect!(model.phenotypes, model.phenoOrder2)
    end
    println("Switched to new YLOSO interval: ", model.phenotypesfilename)
end
# Most variants seems to have 2 as major allele, 0 as minor, here we flip them.
function flip_alleles!(X)
    for t in eachindex(X)
        if X[t] == 0
            X[t] = 2
        elseif X[t] == 2
            X[t] = 0
        end
    end
end

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
module ReferenceTests

using Test
using ImageCore
using Distances
using FileIO
using ImageInTerminal
using SHA
using DeepDiffs
using Random

export
    @withcolor,
    @io2str,
    @test_reference,
    psnr_equality

include("utils.jl")
include("test_reference.jl")
include("fileio.jl")
include("equality_metrics.jl")
include("render.jl")

end # module


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

module ExtremeStats

using Distributions
using JuMP, Ipopt
using RecipesBase

# implement fitting methods
import Distributions: fit_mle

include("maxima.jl")
include("fitting.jl")
include("stats.jl")

# plot recipes
include("plotrecipes/return_levels.jl")
include("plotrecipes/mean_excess.jl")
include("plotrecipes/pareto_quantile.jl")

export
  # maxima types
  BlockMaxima,
  PeakOverThreshold,

  # statistics
  returnlevels,
  meanexcess

end


text:
function benchmark_carrier_replica()
    results_min = DataFrame(
        Samples = SAMPLES,
        sCPU_time = zeros(Float64,length(SAMPLES)),
        GPU_time = zeros(Float64,length(SAMPLES)),
        sGPU_time = zeros(Float64,length(SAMPLES)),
    )
    results_med = DataFrame(
        Samples = SAMPLES,
        sCPU_time = zeros(Float64,length(SAMPLES)),
        GPU_time = zeros(Float64,length(SAMPLES)),
        sGPU_time = zeros(Float64,length(SAMPLES)),
    )
    results_mean = DataFrame(
        Samples = SAMPLES,
        sCPU_time = zeros(Float64,length(SAMPLES)),
        GPU_time = zeros(Float64,length(SAMPLES)),
        sGPU_time = zeros(Float64,length(SAMPLES)),
    )
    rowpos = Int32(1)
    for N in SAMPLES
        # init signals
        cpucarrier = zeros(Complex{Float32}, N)
        gpucarrier = CuArray{ComplexF32}(cpucarrier)

        # init signals in StructArrays
        scpucarrier = StructArray{ComplexF32}(
            (real(cpucarrier),imag(cpucarrier))
        )
        sgpucarrier = StructArray{ComplexF32}(
            (real(gpucarrier),imag(gpucarrier))
        )

        # do the benchmark
        println("Benchmarking carrier_replica on CPU: StructArray{ComplexF32} ", N," samples...")
        result = @benchmark cpu_gen_carrier_replica!(
            $scpucarrier,
            1500,
            2.5e6,
            0,
            Val(7),
            1,
            $N
        )
        println(minimum(result).time)
        results_min.sCPU_time[rowpos] = minimum(result).time
        results_med.sCPU_time[rowpos] = median(result).time
        results_mean.sCPU_time[rowpos] = mean(result).time
        println("Benchmarking carrier_replica on GPU: CuArray{ComplexF32} ", N, " samples...")
        result = @benchmark CUDA.@sync gpu_gen_carrier_replica!(
            $gpucarrier,
            1500,
            2.5e6,
            0,
            Val(7),
            1,
            $N
        )
        println(minimum(result).time)
        results_min.GPU_time[rowpos] = minimum(result).time
        results_med.GPU_time[rowpos] = median(result).time
        results_mean.GPU_time[rowpos] = mean(result).time
        println("Benchmarking carrier_replica on GPU: StructArray{ComplexF32}(CuArray,CuArray) ", N, " samples...")
        result = @benchmark CUDA.@sync gpu_gen_carrier_replica!(
            $sgpucarrier,
            1500,
            2.5e6,
            0,
            Val(7),
            1,
            $N
        )
        println(minimum(result).time)
        results_min.sGPU_time[rowpos] = minimum(result).time
        results_med.sGPU_time[rowpos] = median(result).time
        results_mean.sGPU_time[rowpos] = mean(result).time
        rowpos += 1
    end
    CSV.write("data/carrierreplica_min.csv", results_min)
    CSV.write("data/carrierreplica_med.csv", results_med)
    CSV.write("data/carrierreplica_mean.csv", results_mean)
end


text:
## We need to use a library for Symbolic calculations, a very common one is SymPy from the python
## language. Luckily Julia has excellent support for it.
using Parameters
using DifferentialEquations
using ForwardDiff
using SymPy

#jl Remove the follow `#jl` to use this file standalone
include("cr_model.jl")

# We want to deal with all of these as symbolic variables, which are not the default type in Julia,
# unlike Mathematica.
@vars R C
@vars r k e a h m
# We need to make a symbolic parameter list, as `CRPar` is numeric
spar = Dict(
    :r => r,
    :k => k,
    :a => a,
    :h => h,
    :e => e,
    :m => m)

# we need to make symbolic versions of the model equations. We do this by calling the function with
# the symbolic parameters. the last parameter could be anythign as the time (`t`) argument is not
# used. I have set it to `NaN` which is a name for not a number.

## CR exp resource, type I consumer
f1, f2 = cr_expI_model([R, C], spar, NaN)
#-
SymPy.solve(f1, R)
#-
SymPy.solve(f1, C)
#-
SymPy.solve(f2, R)
#-
SymPy.solve(f2, C)
#-

## CR log resource, type I consumer
g1, g2 = cr_logI_model([R, C], spar, NaN)
#-
SymPy.solve(g1, R)
#-
SymPy.solve(g1, C)
#-
SymPy.solve(g2, R)
#-
SymPy.solve(g2, C)

## CR log resource, type II consumer
l1, l2 = cr_logII_model([R, C], spar, NaN)
#-
SymPy.solve(l1, R)
#-
SymPy.solve(l1, C)
#-
SymPy.solve(l2, R)
#-
SymPy.solve(l2, C)

## CR exp resource, type II consumer
n1, n2 = cr_expII_model([R, C], spar, NaN)
#-
SymPy.solve(n1, R)
#-
SymPy.solve(n1, C)
#-
SymPy.solve(n2, R)
#-
SymPy.solve(n2, C)

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using Base.Test
using RandomMatrices

@test Wigner{3} == GaussianHermite{3}

n = 25

for (β, T, N) in [(1, Real, n), (2, Complex, n), (4, Complex, 2n)]
    d = Wigner(β)
    A = rand(d, n)
    @test eltype(A) <: T
    @test size(A) == (N, N)

    At = tridrand(d, n)
    @test eltype(At) <: Real
    @test size(At) == (n, n)

    vals = eigvalrand(d, n)
    @test eltype(vals) <: Real
    @test length(vals) == n

    vd = RandomMatrices.VandermondeDeterminant(vals, β)
    @test isa(vd, Real)

    ed = eigvaljpdf(d, vals)
    @test isa(ed, Real)
end


text:
#####
##### Tests for flow defaults
#####

struct CircleTest <: Manifold end

struct CircleDecoratorTest <: Manifold end

base_manifold(M::CircleDecoratorTest) = CircleTest()

@traitimpl IsDecoratorManifold{CircleDecoratorTest}

function exp!(M::CircleTest, y, x, v, t::Real)
    θ = norm(v)
    u = θ == 0 ? zero(v) : v ./ θ
    α = θ * t
    y .= cos(α) * x + sin(α) * u
    return y
end

types = [Vector{Float64},
         SizedVector{2, Float64},
         MVector{2, Float64},
         Vector{Float32},
         SizedVector{2, Float32},
         MVector{2, Float32}]

@testset "geodesic_flow default" begin
    M = CircleTest()
    x₀ = normalize(randn(2))
    v₀ = normalize((I-x₀*x₀') * randn(2))

    for T in types
        @testset "Type $T" begin
            x, v = geodesic_flow(M, convert(T, x₀), convert(T, v₀), eltype(T)(π))

            @test x ≈ -x₀
            @test v ≈ -v₀
        end
    end
end

@testset "geodesic_flow decorator default" begin
    M = CircleDecoratorTest()
    x₀ = normalize(randn(2))
    v₀ = normalize((I-x₀*x₀') * randn(2))

    for T in types
        @testset "Type $T" begin
            x, v = geodesic_flow(M, convert(T, x₀), convert(T, v₀), eltype(T)(π))

            @test x ≈ -x₀
            @test v ≈ -v₀
        end
    end
end


text:
SUITE["tz_data"] = BenchmarkGroup()

SUITE["tz_data"]["parse_components"] = begin
    @benchmarkable parse_components($("2010 Europe/Warsaw"), $(DateFormat("yyyy ZZZ")))
end


text:
include("./util.jl")
include("./parsers.jl")
include("./Algorithms/TEAVAR.jl")

env = Gurobi.Env()

topology = "B4"
weibull = true
shape = 0.8
scale = 0.01
paths = "SMORE"
demand_num = 1
beta=0.90

links, capacity, link_probs, nodes = readTopology(topology)
demand, flows = readDemand("$(topology)/demand", length(nodes), demand_num, matrix=true)
T, Tf, k = parsePaths("$(topology)/paths/$(paths)", links, flows)

if weibull
  probabilities = weibullProbs(length(links), shape=shape, scale=scale)
else
  probabilities = link_probs
end
cutoff = (sum(probabilities)/length(probabilities))^2
scenarios, scenario_probs = subScenariosRecursion(probabilities, cutoff)
TEAVAR(env, links, capacity, flows, demand, beta, k,
    T, Tf, scenarios, scenario_probs, explain=true, verbose=true,
        utilization=true)


text:
include("1.1.1 Linear/.jl")
include("1.1.2 Quadratic/.jl")
include("1.1.3 General/.jl")
include("1.1.4 Improper/.jl")

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/.jl


text:
## ----------------------------------------------------------------------------
const ME_MAX_POL                = :ME_MAX_POL                 #
const ME_MAX_POL_B0             = :ME_MAX_POL_B0              #
const ME_Z_EXPECTED_G_BOUNDED   = :ME_Z_EXPECTED_G_BOUNDED    # Match ME and Dy biom average and constraint av_ug


text:
module HelperFunctionsTests

using Test
using Gridap.Helpers

for D = 0:5
  @test tfill(2,Val(D)) == tuple(fill(2,D)...)
end

@test 1 == get_val_parameter(Val{1}())

end # module


text:
# create a dictionary with all the possible models
function get_car_models(env::UrbanEnv, get_model::Function)
    d = Dict{SVector{2, LaneTag}, DriverModel}()

    #
    r1 = SVector(LaneTag(1,1), LaneTag(2,1))
    d[r1] = get_model(env, r1)

    r2 = SVector(LaneTag(1,1), LaneTag(5,1))
    d[r2] = get_model(env, r2)

    r3 = SVector(LaneTag(3,1), LaneTag(4,1))
    d[r3] = get_model(env, r3)

    r4 = SVector(LaneTag(3,1), LaneTag(5,1))
    d[r4] = get_model(env, r4)

    return d
end

function set_car_model!(mdp::PedCarMDP, s::PedCarMDPState, a::PedCarMDPAction)
    reset_hidden_state!(mdp.car_models[s.route])
    observe!(mdp.car_models[s.route], s.ego, s.car, s.ped, mdp.env.roadway)
end

function AutomotiveDrivingModels.observe!(model::UrbanDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          ped::VehicleState,
                                          roadway::Roadway)
    AutomotiveDrivingModels.observe!(model.navigator, ego, car, roadway)
    AutomotiveDrivingModels.observe!(model.intersection_driver, ego, car, roadway)
    for driver in model.crosswalk_drivers
        AutomotiveDrivingModels.observe!(driver, ego, car, ped, roadway)
    end
    a_lon_crosswalks = minimum([driver.a.a_lon for driver in model.crosswalk_drivers])
    a_lon = min(model.intersection_driver.a.a_lon, min(a_lon_crosswalks, model.navigator.a))
    # println("UrbanDriver says ", a_lon)

    model.a = LonAccelDirection(a_lon, model.navigator.dir)
    return model.a
end

# function AutomotiveDrivingModels.observe!(model::TTCIntersectionDriver,
#                                           ego::VehicleState,
#                                           car::VehicleState,
#                                           roadway::Roadway)

#     model.priority = ttc_check(model, car, ego, roadway)
#     passed = has_passed(model, car, ego, roadway)
#     a_lon_idm = model.navigator.a
#     right_of_way = model.priorities[(model.navigator.route[1].tag,model.navigator.route[end].tag)]
#     if !model.priority && !passed && engaged(model, car, ego, roadway) # emergency break
#         a_lon = -model.navigator.d_max
#     elseif !model.priority && !passed
#         a_lon = min(a_lon_idm, AutomotivePOMDPs.stop_at_end(model, car, roadway))
#     else
#         a_lon = a_lon_idm
#     end
#     # to get out of the Emergency break
#     if !model.priority && !passed && engaged(model, car, ego, roadway) && right_of_way
#         a_lon = a_lon_idm
#     end
#     model.a = LonAccelDirection(a_lon, model.navigator.dir)
# end

function AutomotiveDrivingModels.observe!(model::TTCIntersectionDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          roadway::Roadway)
    a_lon_idm = model.navigator.a
    passed =  has_passed(model, car, ego, roadway)
    is_engaged = engaged(model, car, ego, roadway)
    right_of_way = model.priorities[(model.navigator.route[1].tag,model.navigator.route[end].tag)]
    is_clogged = is_intersection_clogged(model, car, ego, roadway)
    ttc = ttc_check(model, car, ego, roadway)
    a_lon = a_lon_idm
    if isempty(model.intersection) || passed
        a_lon = a_lon_idm
    elseif !passed
        if right_of_way #&& !is_clogged
            if is_clogged && !passed && is_engaged && !model.stop && !isapprox(car.v, 0.)
                # println("Vehicle $egoid : emergency break")
                a_lon = -model.navigator.d_max
            else
                a_lon = a_lon_idm
            end
        else # left turn
            if !ttc && !is_engaged  # before left turn
                a_lon = min(a_lon_idm, AutomotivePOMDPs.stop_at_end(model, car, roadway))
                # println("!ttc && !engaged: ", a_lon)
            elseif is_clogged && !passed && is_engaged && !isapprox(car.v, 0.) #!ttc && !passed && is_engaged || (is_clogged && is_engaged)
                # println("Vehicle $egoid : emergency break")
                a_lon = -model.navigator.d_max
            elseif ttc
                a_lon = a_lon_idm
            end
        end
    end
    model.a = LonAccelDirection(a_lon, model.navigator.dir)
end

function AutomotiveDrivingModels.observe!(model::CrosswalkDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          ped::VehicleState,
                                          roadway::Roadway)
    a_lon =0.
    a_lon_idm = model.navigator.a_max
    if !model.yield
        a_lon = a_lon_idm
    else
        model.priority = !is_crossing(ped, model.crosswalk, model.conflict_lanes, roadway) ||
        AutomotivePOMDPs.has_passed(model, car, roadway)
        if !model.priority
            dist_to_cw = AutomotivePOMDPs.get_distance_to_crosswalk(model, car, roadway, -model.stop_delta)
            a_lon = min(a_lon_idm, AutomotivePOMDPs.stop_at_dist(model, car, dist_to_cw))
            # println("dist to cw $dist_to_cw, crosswalk $(model.crosswalk.tag), $(get_lane(roadway, car) ∈ model.intersection_entrances)")
            # println("priority $(model.priority), crossing $(is_crossing(ped, model.crosswalk, model.conflict_lanes, roadway)), a_lon $a_lon, v $(car.v)")
        else
            a_lon = a_lon_idm # just idm
        end
    end
    model.a = LonAccelDirection(a_lon, model.navigator.dir)
end

function AutomotiveDrivingModels.observe!(model::RouteFollowingIDM,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          roadway::Roadway)
    Δv = model.v_des - car.v
    acc = Δv*model.k_spd
    model.a = clamp(acc, -model.d_max, model.a_max)
    cur_lane = get_lane(roadway, car)
    set_direction!(model, cur_lane, roadway)
    return model
end

function is_crossing(ped::VehicleState, crosswalk::Lane, conflict_lanes::Vector{Lane{T}}, roadway::Roadway) where T<:Real
    ped_lane = get_lane(roadway, ped)
    if ped_lane.tag != crosswalk.tag
        return false
    end
    if 4.0 <= ped.posF.s <= 10.
        return true
    end
    # at this point, the pedestrian is not on the road
    # check if the pedestrian is going to cross or not
    if AutomotivePOMDPs.direction_from_center(ped, crosswalk) > 0. && get_lane(roadway, ped).tag == crosswalk.tag
        return true
    end
    return false
end

function ttc_check(model::TTCIntersectionDriver, ego::VehicleState, car::VehicleState, roadway::Roadway)
    inter_width=6.0
    posF = car.posF
    int_x, int_y, int_θ = model.intersection_pos
    Δ = (model.intersection_pos.x - car.posG.x)^2 + (model.intersection_pos.y - car.posG.y)^2
    if Δ < inter_width^2 # vehicle is in the middle
        ttc = 0.
    else
        ttc = sqrt(Δ)/car.v
    end

    if 0 <= ttc < model.ttc_threshold
        return false
    else
        return true
    end

end

function is_intersection_clogged(model::TTCIntersectionDriver, ego::VehicleState, car::VehicleState, roadway::Roadway)
    inter_width=5.0
    posF = car.posF
    int_x, int_y, int_θ = model.intersection_pos
    Δ = (model.intersection_pos.x - car.posG.x)^2 + (model.intersection_pos.y - car.posG.y)^2
    return Δ < inter_width^2 # vehicle is in the middle
end

function has_passed(model::TTCIntersectionDriver, ego::VehicleState, car::VehicleState, roadway::Roadway)
    lane = get_lane(roadway, ego)
    inter_to_car = ego.posG - model.intersection_pos
    car_vec = ego.posG + polar(VehicleDef().length/2, ego.posG.θ) - ego.posG
    has_passed = dot(inter_to_car, car_vec) > 0. || isempty(lane.exits)
    return has_passed
end

function engaged(model::TTCIntersectionDriver, ego::VehicleState, car::VehicleState, roadway::Roadway)
    inter_width = 7.5 #todo parameterized
    if normsquared(VecE2(model.intersection_pos - ego.posG)) < inter_width^2
        return true
    end
    return false
end

function is_neighbor_fore_along_lane(ego::VehicleState, car::VehicleState, roadway::Roadway, ego_def::VehicleDef = VehicleDef())
    tag_start = ego.posF.roadind.tag
    s_base = ego.posF.s + ego_def.length/2*cos(ego.posF.ϕ)
    tag_start = ego.posF.roadind.tag
    targetpoint_primary = VehicleTargetPointRear()
    targetpoint_valid = VehicleTargetPointFront()

    is_neighbor_fore_along_lane(ego, car, roadway, tag_start, s_base)
end

function is_neighbor_fore_along_lane(
    ego::VehicleState,
    veh::VehicleState,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64;
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    )

    best_ind = 0
    best_dist = max_distance_fore
    tag_target = tag_start
    lane = roadway[tag_target]

    dist_searched = 0.0
    while dist_searched < max_distance_fore

        s_adjust = NaN
        if veh.posF.roadind.tag == tag_target
            s_adjust = 0.0

        elseif is_between_segments_hi(veh.posF.roadind.ind, lane.curve) &&
                is_in_entrances(roadway[tag_target], veh.posF.roadind.tag)

            distance_between_lanes = norm(VecE2(roadway[tag_target].curve[1].pos - roadway[veh.posF.roadind.tag].curve[end].pos))
            s_adjust = -(roadway[veh.posF.roadind.tag].curve[end].s + distance_between_lanes)

        elseif is_between_segments_lo(veh.posF.roadind.ind) &&
                is_in_exits(roadway[tag_target], veh.posF.roadind.tag)

            distance_between_lanes = norm(VecE2(roadway[tag_target].curve[end].pos - roadway[veh.posF.roadind.tag].curve[1].pos))
            s_adjust = roadway[tag_target].curve[end].s + distance_between_lanes
        end

        if !isnan(s_adjust)
            s_valid = veh.posF.s + VehicleDef().length/2*cos(veh.posF.ϕ) + s_adjust
            dist_valid = s_valid - s_base + dist_searched
            if dist_valid ≥ 0.0
                s_primary = veh.posF.s - VehicleDef().length/2*cos(veh.posF.ϕ) + s_adjust
                dist = s_primary - s_base + dist_searched
                if dist < best_dist
                    best_dist = dist
                    best_ind = 1
                end
            end
        end
        if best_ind != 0
            break
        end

        if !has_next(lane) ||
           (tag_target == tag_start && dist_searched != 0.0) # exit after visiting this lane a 2nd time
            break
        end

        dist_searched += (lane.curve[end].s - s_base)
        s_base = -norm(VecE2(lane.curve[end].pos - next_lane_point(lane, roadway).pos)) # negative distance between lanes
        tag_target = next_lane(lane, roadway).tag
    end
    NeighborLongitudinalResult(best_ind, best_dist)
end


text:
#=
# Compute the node locations for the reference element
=#
include("jacobi_gauss_quad.jl");
include("triangle_nodes.jl");
include("tet_nodes.jl");

function refel_nodes!(refel, nodetype)
    if refel.dim == 0
        # 0D is a point so some of this doesn't make sense
        refel.r1d = [0];
        refel.wr1d = [1];
        refel.g1d = [0];
        refel.wg1d = [1];
        refel.r = [0];
        refel.wr = [1];
        refel.g = [0];
        refel.wg = [1];

    elseif refel.dim == 1
        # 1D has line segments
        if nodetype == UNIFORM
            refel.r1d = Array(-1:(2/(refel.Np-1)):1);
            refel.wr1d = ones(length(refel.r1d)) ./ length(refel.r1d);
        elseif nodetype == GAUSS
            (r,w) = jacobi_gauss_quad(0,0,refel.N);
            refel.r1d = r;
            refel.wr1d = w;
        elseif nodetype == LOBATTO
            if refel.N == 1
                refel.r1d = [-1; 1];
                refel.wr1d = [1; 1];
            else
                # (r,w) = jacobi_gauss_quad(1,1,refel.N-2);
                # refel.r1d = [-1; r ; 1];

                # # compute the weights
                # w = jacobi_polynomial(refel.r1d, 0, 0, refel.N)
                # adgammaN = (2*refel.N + 1) / (refel.N * (refel.N + 1))
                # w = w.*w
                # w = adgammaN./w

                # refel.wr1d = w;
                (refel.r1d, refel.wr1d) = jacobi_LGL_quad(refel.N);
            end
        end
        # Then find Gauss points
        (g,w) = jacobi_gauss_quad(0,0,refel.N);
        refel.g1d = g;
        refel.wg1d = w;

        # 1D is easy
        refel.g = refel.g1d;
        refel.wg = refel.wg1d;
        refel.r = refel.r1d;
        refel.wr = refel.wr1d;

        # surface is just points
        refel.face2local = [[1], [refel.Np]];
        refel.surf_r = [[-1], [1]];
        refel.surf_wr = [[1], [1]];
        refel.surf_g = [[-1], [1]];
        refel.surf_wg = [[1], [1]];

    elseif refel.dim == 2
        # 2D has triangles and quads
        if refel.Nfaces == 3 # triangles
            triangle_refel_nodes!(refel);
        else # quads
            if nodetype == UNIFORM
                refel.r1d = Array(-1:(2/(refel.Np-1)):1);
                refel.wr1d = ones(length(refel.r1d)) ./ length(refel.r1d);
            elseif nodetype == GAUSS
                (r,w) = jacobi_gauss_quad(0,0,refel.N);
                refel.r1d = r;
                refel.wr1d = w;
            elseif nodetype == LOBATTO
                if refel.N == 1
                    refel.r1d = [-1; 1];
                    refel.wr1d = [1; 1];
                else
                    (r,w) = jacobi_gauss_quad(1,1,refel.N-2);
                    refel.r1d = [-1; r ; 1];

                    # compute the weights
                    w = jacobi_polynomial(refel.r1d, 0, 0, refel.N);
                    adgammaN = (2*refel.N + 1) / (refel.N * (refel.N + 1));
                    w = w.*w;
                    w = adgammaN./w;

                    refel.wr1d = w;
                end
            end
            # Then find Gauss quadrature points
            (g,w) = jacobi_gauss_quad(0,0,refel.N);
            refel.g1d = g;
            refel.wg1d = w;

            # r and w
            refel.r = zeros(refel.Np,2);
            refel.wr = zeros(refel.Np);
            refel.g = zeros(refel.Np,2);
            refel.wg = zeros(refel.Np);
            n1d = length(refel.r1d);
            for j=1:n1d
                for i=1:n1d
                    k = (j-1)*n1d + i;
                    refel.r[k,:] = [refel.r1d[i]; refel.r1d[j]];
                    refel.wr[k] = refel.wr1d[i] * refel.wr1d[j];
                    refel.g[k,:] = [refel.g1d[i]; refel.g1d[j]];
                    refel.wg[k] = refel.wg1d[i] * refel.wg1d[j];
                end
            end

            ### surface lines ################
            # First get the face2local maps
            tol = 1e-8;
            qf1(x) = abs(x[1] + 1) < tol;
            qf2(x) = abs(x[2] + 1) < tol;
            qf3(x) = abs(x[1] - 1) < tol;
            qf4(x) = abs(x[2] - 1) < tol;
            refel.face2local = [get_face2local_map(refel.r, qf1),
                                get_face2local_map(refel.r, qf2),
                                get_face2local_map(refel.r, qf3),
                                get_face2local_map(refel.r, qf4)];

            tmp = zeros(2,0);
            refel.surf_r = [tmp, tmp, tmp, tmp];
            refel.surf_g = [tmp, tmp, tmp, tmp];

            for fi=1:4
                refel.surf_r[fi] = refel.r[refel.face2local[fi], :];
                refel.surf_g[fi] = refel.g[refel.face2local[fi], :];
            end
            refel.surf_wr = [refel.wr1d, refel.wr1d, refel.wr1d, refel.wr1d];
            refel.surf_wg = [refel.wg1d, refel.wg1d, refel.wg1d, refel.wg1d];
        end

    elseif refel.dim == 3
        # 3D has tets, hexs and prisms
        if refel.Nfaces == 4 # tets
            tetrahedron_refel_nodes!(refel);
        else # hexs
            if nodetype == UNIFORM
                refel.r1d = Array(-1:(2/(refel.Np-1)):1);
                refel.wr1d = ones(length(refel.r1d)) ./ length(refel.r1d);
            elseif nodetype == GAUSS
                (r,w) = jacobi_gauss_quad(0,0,refel.N);
                refel.r1d = r;
                refel.wr1d = w;
            elseif nodetype == LOBATTO
                if refel.N == 1
                    refel.r1d = [-1; 1];
                    refel.wr1d = [1; 1];
                else
                    (r,w) = jacobi_gauss_quad(1,1,refel.N-2);
                    refel.r1d = [-1; r ; 1];

                    # compute the weights
                    w = jacobi_polynomial(refel.r1d, 0, 0, refel.N);
                    adgammaN = (2*refel.N + 1) / (refel.N * (refel.N + 1));
                    w = w.*w;
                    w = adgammaN./w;

                    refel.wr1d = w;
                end
            end
            # Then find Gauss points
            (g,w) = jacobi_gauss_quad(0,0,refel.N);
            refel.g1d = g;
            refel.wg1d = w;

            # r and w
            refel.r = zeros(refel.Np,3);
            refel.wr = zeros(refel.Np);
            refel.g = zeros(refel.Np,3);
            refel.wg = zeros(refel.Np);
            n1d = length(refel.r1d);
            for k=1:n1d
                for j=1:n1d
                    for i=1:n1d
                        ind = (k-1)*n1d*n1d + (j-1)*n1d + i;
                        refel.r[ind,:] = [refel.r1d[i]; refel.r1d[j]; refel.r1d[k]];
                        refel.wr[ind] = refel.wr1d[i] * refel.wr1d[j] * refel.wr1d[k];
                        refel.g[ind,:] = [refel.g1d[i]; refel.g1d[j]; refel.g1d[k]];
                        refel.wg[ind] = refel.wg1d[i] * refel.wg1d[j] * refel.wg1d[k];
                    end
                end
            end

            ### surface quads ################
            # First get the face2local maps
            tol = 1e-8;
            hf1(x) = abs(x[1] + 1) < tol;
            hf2(x) = abs(x[2] + 1) < tol;
            hf3(x) = abs(x[3] + 1) < tol;
            hf4(x) = abs(x[1] - 1) < tol;
            hf5(x) = abs(x[2] - 1) < tol;
            hf6(x) = abs(x[3] - 1) < tol;
            refel.face2local = [get_face2local_map(refel.r, hf1),
                                get_face2local_map(refel.r, hf2),
                                get_face2local_map(refel.r, hf3),
                                get_face2local_map(refel.r, hf4),
                                get_face2local_map(refel.r, hf5),
                                get_face2local_map(refel.r, hf6)];

            tmp = zeros(3,0);
            refel.surf_r = [tmp, tmp, tmp, tmp, tmp, tmp];
            refel.surf_g = [tmp, tmp, tmp, tmp, tmp, tmp];

            for fi=1:6
                refel.surf_r[fi] = refel.r[refel.face2local[fi], :];
                refel.surf_g[fi] = refel.g[refel.face2local[fi], :];
            end

            # The weights are a little more tricky
            fwr = zeros(refel.Nfp[1]);
            fwg = zeros(refel.Nfp[1]);
            n1d = length(refel.r1d);
            for j=1:n1d
                for i=1:n1d
                    k = (j-1)*n1d + i;
                    fwr[k] = refel.wr1d[i] * refel.wr1d[j];
                    fwg[k] = refel.wg1d[i] * refel.wg1d[j];
                end
            end
            refel.surf_wr = [fwr, fwr, fwr, fwr, fwr, fwr];
            refel.surf_wg = [fwg, fwg, fwg, fwg, fwg, fwg];
        end

    else
        # Not ready for 4D
    end
end

# Determines surface nodes using a supplied function.
# Returns a face to local map. The order of which will be omnipotent.
function get_face2local_map(r, compare)
    n = size(r,1);
    map = zeros(Int,n);
    nf = 0;
    for i=1:n
        if compare(r[i,:])
            nf = nf+1;
            map[nf] = i;
        end
    end

    return map[1:nf];
end

# Returns the global node locations for elemental nodes
# Has size Np,dim
function get_node_coords(vx, elnodes)
    x = [];
    if refel.dim == 1
        hx = vx[2] - vx[1];
        x = vx[1] .+ (elnodes .+ 1) .* (hx*0.5);
    elseif refel.dim == 2
        # for now assume rectangular quads
        x = zeros(size(elnodes));
        hx = max(vx[2,1]-vx[1,1], vx[3,1]-vx[1,1]);
        hy = max(vx[2,2]-vx[1,2], vx[3,2]-vx[1,2]);
        x[:,1] = vx[1,1] .+ (elnodes[:,1] .+ 1) .*(hx*0.5);
        x[:,2] = vx[1,2] .+ (elnodes[:,2] .+ 1) .*(hy*0.5);
    elseif refel.dim == 3
        # for now assume regular hexs
        x = zeros(size(elnodes));
        hx = max(vx[2,1]-vx[1,1], vx[3,1]-vx[1,1]);
        hy = max(vx[2,2]-vx[1,2], vx[3,2]-vx[1,2]);
        x[:,1] = vx[1,1] .+ (elnodes[:,1] .+ 1) .*(hx*0.5);
        x[:,2] = vx[1,2] .+ (elnodes[:,2] .+ 1) .*(hy*0.5);
        x[:,3] = vx[1,3] .+ (elnodes[:,3] .+ 1) .*(hy*0.5);
    end

    return x;
end

# Returns the volume of an element
function get_volume(refel, vx)
    vol = 0;
    if refel.dim == 1
        vol = abs(vx[2]-vx[1]);
    elseif refel.dim == 2
        # for now assume rectangular quads
        hx = max(vx[2,1]-vx[1,1], vx[3,1]-vx[1,1]);
        hy = max(vx[2,2]-vx[1,2], vx[3,2]-vx[1,2]);
        vol = abs(hx*hy);
    elseif refel.dim == 3

    end

    return vol;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
@params struct PolynomialFit{T}
    data::Tuple{Matrix{T}, Vector{T}}
    order::Ref{Int}
    coeffs::AbstractVector{T}
    ratio::T
end
PolynomialFit{T}(n, ratio) where T = PolynomialFit{T}((zeros(T, n+1, n+1), zeros(T, n+1)), Ref(-1), zeros(T, n+1), T(ratio))

function Base.convert(::Type{PolynomialFit}, n::Int)
    PolynomialFit{Float64}(n)
end
function Base.convert(::Type{PolynomialFit{T}}, n::Int) where T
    PolynomialFit{T}((zeros(T, n+1, n+1), zeros(T, n+1)), Ref(-1),
        zeros(T, n+1), T(0.5))
end
Base.eltype(::PolynomialFit{T}) where T = T

# nth order derivatives of polynomial function f using f'''(x) notation
@inline Base.transpose(p::PolynomialFit) = Derivative(p)
@inline (poly::PolynomialFit)(x) = poly(x, 0)
@inline (poly::PolynomialFit)(x, n) = nthderivativeat(poly, n, x)
function nthderivativeat(poly::PolynomialFit{T}, n, x) where T
    o = zero(T)
    for j in (n+1):order(poly)+1
        p = n == 0 ? one(T) : prod(i->(j-i), 1:n)
        o += p*poly.coeffs[j]*x^(j-1-n)
    end
    return o
end

maxorder(p::PolynomialFit) = length(p.coeffs)-1
order(p::PolynomialFit) = p.order[]
increaseorder!(p::PolynomialFit) = p.order[]+= 1
reset!(p::PolynomialFit) = p.order[] = -1

previous(p::PolynomialFit, i) = i == 1 ? order(p) + 1 : i - 1
next(p::PolynomialFit, i) = i == order(p) + 1 ? 1 : i + 1

function newvalues!(p::PolynomialFit, X, V)
    for i in 1:length(X)
        newvalue!(p, X[i], V[i])
    end
end
function newvalue!(p::PolynomialFit{T}, x, v) where T
    if order(p) < maxorder(p)
        increaseorder!(p)
    end
    xs = [x^i for i in 0:maxorder(p)]
    p.data[1] .= (1 .- p.ratio) .* p.data[1] .+ p.ratio .* (xs * xs')
    p.data[2] .= (1 .- p.ratio) .* p.data[2] .+ p.ratio .* xs .* v
end

function newderivatives!(p::PolynomialFit, X, V)
    for i in 1:length(X)
        newderivative!(p, X[i], V[i])
    end
end
function newderivative!(p::PolynomialFit{T}, x, v) where T
    if order(p) < maxorder(p)
        increaseorder!(p)
    end
    xs = [(j-1)x^(j-2) for j in 1:maxorder(p)+1]
    p.data[1] .= (1 .- p.ratio) .* p.data[1] .+ p.ratio .* (xs * xs')
    p.data[2] .= (1 .- p.ratio) .* p.data[2] .+ p.ratio .* xs .* v
end

function solve!(p::PolynomialFit)
    n = order(p)+1
    p.coeffs[1:n] .= pinv(p.data[1][1:n, 1:n]) * p.data[2][1:n]
    @view(p.coeffs[1:n])
end

function roots(poly::Union{TP, Derivative{TP}}, min=T(0), max=T(15)) where {T, TP<:PolynomialFit{T}}
    root_ranges = IntervalRootFinding.roots(x->poly(x), min..max)
    return _filter(root_ranges)
end

Base.:/(p1::PolynomialFit, p2::PolynomialFit) = PolynomialQuotient(p1, p2)

Base.:(==)(p::PolynomialFit, n) = PolynomialEquality(p, n)

function roots(poly::PolynomialEquality{TP}, min=eltype(poly)(0), max=eltype(poly)(15)) where {TP<:Union{PolynomialFit, Derivative{<:PolynomialFit}}}
    root_ranges = IntervalRootFinding.roots(x->(poly.LHS(x) - poly.RHS), min..max)
    return _filter(root_ranges)
end

# TO-DO: derivative of a quotient using the quotient rule


text:
#!/usr/bin/env julia

using NarrativeTest
runtests()


text:
#' ---
#' title: Benders Decomposition
#' ---

#' **Originally Contributed by**: Mathieu Besançon

#' This notebook describes how to implement the Benders decomposition in JuMP using lazy constraints.
#' We keep the same notation and problem form as the first notebook *Benders decomposition*.

#' \begin{align*}
#' & \text{maximize} \quad &&c_1^T x+c_2^T v \\
#' & \text{subject to} \quad &&A_1 x+ A_2 v \preceq b \\
#' & &&x \succeq 0, x \in \mathbb{Z}^n \\
#' & &&v \succeq 0, v \in \mathbb{R}^p \\
#' \end{align*}

#' where $b \in \mathbb{R}^m$, $A_1 \in \mathbb{R}^{m \times n}$, $A_2 \in \mathbb{R}^{m \times p}$ and
#' \mathbb{Z} is the set of integers.
#' Here the symbol $\succeq$ ($\preceq$) stands for element-wise greater (less) than or equal to.
#' Any mixed integer programming problem can be written in the form above.

#' For a detailed explanation on the Benders decomposition algorithm, see the introduction notebook.

#' ### Lazy constraints

#' Some optimization solvers allow users to interact with them during the solution process
#' by providing user-defined functions which are triggered under certain conditions.
#' The generic term for these functions is **callback**.
#' In integer optimization, the main callback types are lazy constraints, user cuts and heuristic solutions.
#' See the [JuMP documentation](https://www.juliaopt.org/JuMP.jl/stable/callbacks/) for an introduction on how to use them.
#' Some callbacks define a new constraint which is only activated when necessary, i.e.
#' when a current solution does not respect them. It can avoid building an optimization model
#' with too many constraints up-front.
#' This is the case for Benders decomposition, since the sub-problem defines an exponential number of primal vertices
#' and therefore dual cuts.
#'
#' A detailed explanation on the distinction between user cuts and lazy constraints is also available on
#' [Paul Rubin's blog](https://orinanobworld.blogspot.com/2012/08/user-cuts-versus-lazy-constraints.html).

#' We use the data from the original notebook and change the solution algorithm to leverage lazy constraints:

#' **Step 1 (Initialization)**
#' **Step 2 (defining the subproblem model)**
#' **Step 3 (registering the lazy constraint of the subproblem)**

#' ### Data for the problem

#' The input data is from page 139, Integer programming by Garfinkel and Nemhauser[[1]](#c1).

c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [1 -3;
     -1 -3]
A2 = [1 -2;
     -1 -1]

M = 1000;

#' [Paul Rubin's blog on Benders Decomposition](http://orinanobworld.blogspot.ca/2011/10/benders-decomposition-then-and-now.html).

# Loading the necessary packages
#-------------------------------
using JuMP
using GLPK
using LinearAlgebra
using Test

#' Subproblem creation

function build_subproblem()
    sub_problem_model = Model(GLPK.Optimizer)
    @variable(sub_problem_model, u[1:dim_u] >= 0)
    @constraint(sub_problem_model, constr_ref_subproblem[j = 1:size(A2, 2)], dot(A2[:, j], u) >= c2[j])
    return (sub_problem_model, u)
end

# Master Problem Description
# --------------------------

master_problem_model = Model(GLPK.Optimizer)
(sub_problem_model, u) = build_subproblem()

# Variable Definition
# ----------------------------------------------------------------
@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)

# Objective Setting
# -----------------
@objective(master_problem_model, Max, t)

print(master_problem_model)

#' Track the calls to the callback

iter_num = 0

#' Define lazy constraints

function benders_lazy_constraint_callback(cb_data)
    global iter_num
    iter_num += 1
    println("Iteration number = ", iter_num)

    x_current = [callback_value(cb_data, x[j]) for j in eachindex(x)]
    fm_current = callback_value(cb_data, t)

    c_sub = b - A1 * x_current
    @objective(sub_problem_model, Min, dot(c1, x_current) + dot(c_sub, u))
    optimize!(sub_problem_model)

    print("\nThe current subproblem model is \n", sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = dot(b, u_current)

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current  ≈  fm_current # we are done
        @info("No additional constraint from the subproblem")
        # println("Optimal solution of the original problem found")
        # println("The optimal objective value t is ", fm_current)
        # println("The optimal x is ", x_current)
        # println("The optimal v is ", dual.(constr_ref_subproblem))
    end

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv = A1' * u_current - c1
        new_optimality_cons = @build_constraint(t + dot(cv, x) <= γ)
        MOI.submit(master_problem_model, MOI.LazyConstraint(cb_data), new_optimality_cons)
    end

    if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1' * u_current
        new_feasibility_cons = @build_constraint(dot(ce, x) <= γ)
        MOI.submit(master_problem_model, MOI.LazyConstraint(cb_data), new_feasibility_cons)
    end
end

MOI.set(master_problem_model, MOI.LazyConstraintCallback(), benders_lazy_constraint_callback)

optimize!(master_problem_model)

t_status = termination_status(master_problem_model)
p_status = primal_status(master_problem_model)

if p_status == MOI.INFEASIBLE_POINT
    println("The problem is infeasible :-(")
end

if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
    fm_current = M
    x_current = M * ones(dim_x)
end

if p_status == MOI.FEASIBLE_POINT
    fm_current = value(t)
    x_current = Float64[]
        for i in 1:dim_x
        push!(x_current, value(x[i]))
    end
end

println("Status of the master problem is ", t_status,
        "\nwith fm_current = ", fm_current,
        "\nx_current = ", x_current)

@test value(t) ≈ -4

#' ### References
#'
#' 1. Garfinkel, R. & Nemhauser, G. L. Integer programming. (Wiley, 1972).


text:
using Documenter
using KSVD

makedocs()


text:
#testbug_GB14
include("GB14.jl")
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

rng = Random.seed!(2019)

gbdt = GB14.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                 # Learning rate
  num_iterations = 100,                # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB14.GBLearner(gbdt, :regression)

x, y = GB14.test_data(500; p = 1000, sigma = sqrt(1))
x_test, y_test = GB14.test_data(500; p = 1000, sigma = sqrt(1))
mu = GB14.fried(x)
mu_test = GB14.fried(x_test)

# initialize the hypers
hypers = GB14.Hypers(x, y, 1, 1, 0.8)
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
hypers.α = 0.95
hypers.s =  ones(size(x, 2)) / size(x, 2)
hypers

# initialize opt
opts = GB14.Opts()
opts.num_burn = 4950# 500 in total now

# test the new function and get the rmse2 and check the Store_s
rmse2, Store_s = GB14.TBD_predict_and_rmse!(
    x,
    y,
    x_test,
    y_test,
    mu_test,
    gbdt,
    hypers,
    opts,
    rng)
a
## p = 100, σ² = 1
# num_burn = 50:
    #  rmse2 = 1.717
    #  sum(hypers.s[1:5]) = 0.457
# num_burn = 500
    # rmse2 = 1.773
    # sum(hypers.s[1:5]) = 0.447
# num_burn = 1000
    # rmse2 =
    # sum(hypers.s[1:5]) =

##### below is for testing, they work fine #####
#for l in collect(1 : opts.num_burn)

#    hypers_copy = deepcopy(hypers)

#    # fit the model
#    model = GB14.fit!(gbl, x, y, hypers_copy)

#    GB14.UpdateS!(hypers_copy) # returns the updated hypers_copy.s
#    if (opts.update_s) # true
#        hypers.s = hypers_copy.s
#    end
#    if (opts.update_varcount) # that is, the var count is going to be set to zero again
#        hypers.varcount = hypers_copy.varcount
#    end

#end
######## end of testing

hypers
sum(hypers.varcount) # 2765(burn 460)
sum(hypers.varcount[1:5]) # 1218(burn 460)
sum(hypers.s[1:5])

hypers.s

# num_burn = 50, sum(hypers.s[1:5]) = 0.49
# num_burn = 500, sum(hypers.s[1:5]) = wait..

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot; matplotlib.use("agg")
using SparseArrays

ADCME.options.sparse.auto_reorder = false

function compute_residual_and_jacobian(k_chip, S)
    # compute r and J in Jx=r for Newton's method
    # read in current step
    k_fem = k_air * constant(ones(ndof))
    k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
    # k_fem = scatter_update(k_fem, chip_fem_idx, k_chip * ones(length(chip_fem_idx))) # if k_chip is a constant
    k_fem = scatter_update(k_fem, chip_fem_idx, k_chip)

    kgauss = dof_to_gauss_points(k_fem, mesh)
    LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, mesh))

    u, v, p, T = S[1:ndof],
        S[ndof+1:2*ndof],
        S[2*ndof+1:2*ndof+nelem],
        S[2*ndof+nelem+1:end]

    grad_u = eval_grad_on_gauss_pts1(u, mesh)
    grad_v = eval_grad_on_gauss_pts1(v, mesh)

    ugauss = dof_to_gauss_points(u, mesh)
    vgauss = dof_to_gauss_points(v, mesh)
    ux, uy, vx, vy = grad_u[:,1], grad_u[:,2], grad_v[:,1], grad_v[:,2]

    # comupute jacobian
    M1 = constant(compute_fem_mass_matrix1(ux, mesh))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, mesh))

    N1 = constant(compute_fem_mass_matrix1(vy, mesh))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), mesh))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, mesh))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, mesh))

    gradT = eval_grad_on_gauss_pts1(T, mesh)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, mesh))
    DV_TY = constant(compute_fem_mass_matrix1(Ty, mesh))

    T_mat = constant(compute_fem_mass_matrix1(-buoyance_coef * constant(ones(ndof)), mesh))
    T_mat = [SparseTensor(spzeros(ndof, ndof)); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(nelem, ndof))]

    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros(ndof, nelem)) M]]

    # compute residual
    interaction = compute_interaction_term(p, mesh)
    f1 = compute_fem_source_term1(ugauss.*ux, mesh)
    f2 = compute_fem_source_term1(vgauss.*uy, mesh)
    f3 = -interaction[1:ndof]
    f4 = Laplace*u
    F = f1 + f2 + f3 + f4

    g1 = compute_fem_source_term1(ugauss.*vx, mesh)
    g2 = compute_fem_source_term1(vgauss.*vy, mesh)
    g3 = -interaction[ndof+1:end]
    g4 = Laplace*v
    T_gauss = dof_to_gauss_points(T, mesh)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, mesh)
    G = g1 + g2 + g3 + g4 + buoyance_term

    H0 = -B * [u;v]

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, mesh) * T - heat_source

    R = [F;G;H0;T0]

    return R, J

end

function solve_one_step(θ, S)
    function f(θ, S)
        r, J = compute_residual_and_jacobian(θ, S)
        J, r = impose_Dirichlet_boundary_conditions(J, r, bd, zeros(length(bd)))
        # residual_norm = norm(r)
        # op = tf.print("residual norm", residual_norm)
        # r = bind(r, op)
        return r, J
    end

    S_new = newton_raphson_with_grad(f, S, θ)
    return S_new
end

function solve_navier_stokes(S0, NT, θ)
    function condition(i, S_arr)
        i <= NT
    end

    function body(i, S_arr)
        S = read(S_arr, i)
        op = tf.print("i=",i)
        i = bind(i, op)
        S_new = solve_one_step(θ, S)
        S_arr = write(S_arr, i+1, S_new)
        return i+1, S_arr
    end

    i = constant(1, dtype=Int32)
    S_arr = TensorArray(NT+1)
    S_arr = write(S_arr, 1, S0)
    _, S = while_loop(condition, body, [i, S_arr])
    S = set_shape(stack(S), NT+1, length(S0))
end


text:
using Documenter, NiLang
using SparseArrays

using Literate
tutorialpath = joinpath(@__DIR__, "src/examples")
sourcepath = joinpath(@__DIR__, "../examples")
for jlfile in ["besselj.jl", "sparse.jl", "sharedwrite.jl", "qr.jl", "port_zygote.jl", "fib.jl", "unitary.jl", "nice.jl", "realnvp.jl", "boxmuller.jl", "lognumber.jl"]
    Literate.markdown(joinpath(sourcepath, jlfile), tutorialpath)
end

makedocs(;
    modules=[NiLang],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "What and Why" => "why.md",
        "Tutorial" => Any[
            "tutorial.md",
            "examples/port_zygote.md",
           ],
        "Examples" => Any[
            "examples/fib.md",
            "examples/besselj.md",
            "examples/sparse.md",
            "examples/lognumber.md",
            "examples/unitary.md",
            "examples/qr.md",
            "examples/nice.md",
            "examples/realnvp.md",
            "examples/boxmuller.md",
           ],
        "API & Manual" => Any[
            "instructions.md",
            "extend.md",
            "examples/sharedwrite.md",
            "api.md",
           ]
    ],
    repo="https://github.com/GiggleLiu/NiLang.jl/blob/{commit}{path}#L{line}",
    sitename="NiLang.jl",
    authors="JinGuo Liu, thautwarm",
)

deploydocs(;
    repo="github.com/GiggleLiu/NiLang.jl",
)