text:
const unary_vector_operators = [:(Base.median),
                                :(StatsBase.mad),
                                :(Base.norm),
                                :(StatsBase.skewness),
                                :(StatsBase.kurtosis)]

# TODO: dist, iqr

const binary_vector_operators = [:(Base.dot),
                                 :(Base.cor),
                                 :(Base.cov),
                                 :(StatsBase.corspearman)]

const rowwise_operators = [:rowminimums,
                           :rowmaxs,
                           :rowprods,
                           :rowsums,
                           :rowmeans,
                           :rowmedians,
                           :rowstds,
                           :rowvars,
                           :rownorms]

# Swap arguments to fname() anywhere in AST. Returns the number of
# arguments swapped
function swapargs(ast::Expr, fname::Union{Expr, Symbol})
    if ast.head == :call &&
       (ast.args[1] == fname ||
        (isa(ast.args[1], Expr) && ast.args[1].head == :curly &&
         ast.args[1].args[1] == fname)) &&
       length(ast.args) == 3

        ast.args[2], ast.args[3] = ast.args[3], ast.args[2]
        1
    elseif ast.head == :where &&
       (ast.args[1] == fname ||
        (isa(ast.args[1], Expr) && ast.args[1].head == :where &&
         ast.args[1].args[1] == fname)) &&
       length(ast.args) == 3

        ast.args[1].args[2], ast.args[1].args[3] = ast.args[1].args[3], ast.args[1].args[2]
        1
    else
        n = 0
        for arg in ast.args
            n += swapargs(arg, fname)
        end
        n
    end
end
function swapargs(ast, fname::Union{Expr, Symbol})
    ast
    0
end

# Return a block consisting of both the given function and a copy of
# the function in which arguments to the function itself and any
# 2-argument calls to a function of the same name are swapped
macro swappable(func, syms...)
    if (func.head != :function && func.head != :(=)) ||
       ((func.args[1].head != :where || length(func.args[1].args[1].args) != 3) &&
        (func.args[1].head != :call || length(func.args[1].args) != 3))
        throw(ArgumentError("@swappable may only be applied to functions of two arguments"))
    end

    func2 = deepcopy(func)
    fname = func2.args[1].args[1]
    if isa(fname, Expr) && fname.head in (:curly, :where)
        fname = fname.args[1]
    end

    for s in unique([fname, syms...])
        if swapargs(func2, s) < 1
            throw(ErrorException("No argument swapped"))
        end
    end
    esc(Expr(:block, func, func2))
end

#
# Unary operator macros for DataArrays
#

# Apply unary operator to non-missing members of a DataArray or
# AbstractDataArray
macro dataarray_unary(f, intype, outtype, N...)
    esc(quote
        function $(f)(d::$(isempty(N) ? :(DataArray{T}) : :(DataArray{T,$(N[1])}))) where {T<:$(intype)}
            data = d.data
            res = similar(data, $(outtype))
            @bitenumerate d.na i na begin
                if !na
                    @inbounds res[i] = $(f)(data[i])
                end
            end
            DataArray(res, copy(d.na))
        end
        function $(f)(adv::$(isempty(N) ? :(AbstractDataArray{T}) : :(AbstractDataArray{T,$(N[1])}))) where {T<:$(intype)}
            res = similar(adv, $(outtype))
            for i = 1:length(adv)
                res[i] = ($f)(adv[i])
            end
            res
        end
    end)
end

#
# Binary operator macros for DataArrays
#

# Binary operators with one scalar argument
macro dataarray_binary_scalar(vectorfunc, scalarfunc, outtype, swappable)
    esc(Expr(:block,
        # DataArray and AbstractDataArray with scalar
        # XXX It would be really nice to make this work with arbitrary
        # types, but doing so results in a bunch of method ambiguity
        # warnings
        Any[
            begin
                if outtype == :nothing
                    outtype = :(promote_op($scalarfunc,
                                           eltype(a), eltype(b)))
                end
                fns = Any[
                    :(function $(vectorfunc)(a::DataArray, b::$t)
                        data = a.data
                        res = similar(data, $outtype)
                        @bitenumerate a.na i na begin
                            if !na
                                @inbounds res[i] = $(scalarfunc)(data[i], b)
                            end
                        end
                        DataArray(res, copy(a.na))
                    end),
                    :(function $(vectorfunc)(a::AbstractDataArray, b::$t)
                        res = similar(a, $outtype)
                        for i = 1:length(a)
                            res[i] = $(scalarfunc)(a[i], b)
                        end
                        res
                    end)
                ]
                if swappable
                    # For /, Array/Number is valid but not Number/Array
                    # All other operators should be swappable
                    map!(x->Expr(:macrocall, Symbol("@swappable"), LineNumberNode(@__LINE__), x, scalarfunc), fns, fns)
                end
                Expr(:block, fns...)
            end
            for t in (:AbstractString, :Number)
        ]...
    ))
end

# Binary operators with two array arguments
macro dataarray_binary_array(vectorfunc, scalarfunc)
    esc(Expr(:block,
        # DataArray with other array
        Any[
            quote
                function $(vectorfunc)(a::$atype, b::$btype)
                    data1 = $(atype == :DataArray || atype == :(DataArray{Bool}) ? :(a.data) : :a)
                    data2 = $(btype == :DataArray || btype == :(DataArray{Bool}) ? :(b.data) : :b)
                    res = Array{promote_op($vectorfunc, eltype(a), eltype(b))}(
                        uninitialized, promote_shape(size(a), size(b)))
                    resna = $narule
                    @bitenumerate resna i na begin
                        if !na
                            @inbounds res[i] = $(scalarfunc)(data1[i], data2[i])
                        end
                    end
                    DataArray(res, resna)
                end
            end
            for (atype, btype, narule) in ((:(DataArray), :(AbstractRange), :(copy(a.na))),
                                           (:(AbstractRange), :(DataArray), :(copy(b.na))),
                                           (:DataArray, :DataArray, :(a.na .| b.na)),
                                           (:DataArray, :AbstractArray, :(copy(a.na))),
                                           (:AbstractArray, :DataArray, :(copy(b.na))))
        ]...,
        # AbstractDataArray with other array
        # Definitinons with DataArray necessary to avoid ambiguity
        Any[
            quote
                function $(vectorfunc)(a::$atype, b::$btype)
                    res = similar($(asim ? :a : :b),
                                  promote_op($vectorfunc, eltype(a), eltype(b)),
                                  promote_shape(size(a), size(b)))
                    for i = 1:length(a)
                        res[i] = $(scalarfunc)(a[i], b[i])
                    end
                    res
                end
            end
            for (asim, atype, btype) in ((true, :AbstractDataArray, :AbstractRange),
                                         (false, :AbstractRange, :AbstractDataArray),
                                         (true, :DataArray, :AbstractDataArray),
                                         (false, :AbstractDataArray, :DataArray),
                                         (true, :AbstractDataArray, :AbstractDataArray),
                                         (true, :AbstractDataArray, :AbstractArray),
                                         (false, :AbstractArray, :AbstractDataArray))
        ]...,
    ))
end

# Unary operators, DataArrays.
@dataarray_unary(+, Any, T)
@dataarray_unary(-, Bool, Int)
@dataarray_unary(-, Any, T)
@dataarray_unary(!, Bool, T)

# Treat ctranspose and * in a special way
for (f, elf) in ((:(Base.ctranspose), :conj), (:(Base.transpose), :identity))
    @eval begin
        function $(f)(d::DataMatrix{T}) where T
            # (c)transpose in Base uses a cache-friendly algorithm for
            # numeric arrays, which is faster than our naive algorithm,
            # but chokes on undefined values in the data array.
            # Fortunately, undefined values can only be present in
            # arrays of non-bits types.
            if isbits(T)
                DataArray($(f)(d.data), transpose(d.na))
            else
                data = d.data
                sz = (size(data, 1), size(data, 2))
                res = similar(data, size(data, 2), size(data, 1))
                @bitenumerate d.na i na begin
                    if !na
                        jnew, inew = ind2sub(sz, i)
                        @inbounds res[inew, jnew] = $(elf)(data[i])
                    end
                end
                DataArray(res, transpose(d.na))
            end
        end
    end
end

# Propagates missings
# For a dissenting view,
# http://radfordneal.wordpress.com/2011/05/21/slowing-down-matrix-multiplication-in-r/
# But we're getting 10x R while maintaining missings
for (adata, bdata) in ((true, false), (false, true), (true, true))
    @eval begin
        function (*)(a::$(adata ? :(Union{DataVector, DataMatrix}) : :(Union{Vector, Matrix})),
                           b::$(bdata ? :(Union{DataVector, DataMatrix}) : :(Union{Vector, Matrix})))
            c = $(adata ? :(a.data) : :a) * $(bdata ? :(b.data) : :b)
            res = DataArray(c, falses(size(c)))
            # Propagation can be made more efficient by storing record of corrupt
            # rows and columns, then doing fast edits.
            $(if adata
                quote
                    n1 = size(a, 1)
                    p1 = size(a, 2)
                    corrupt_rows = falses(n1)
                    for j in 1:p1, i in 1:n1
                        # Propagate missings
                        # Corrupt all rows based on i
                        corrupt_rows[i] |= a.na[i, j]
                    end
                    res.na[corrupt_rows, :] = true
                end
            end)
            $(if bdata
                quote
                    n2 = size(b, 1)
                    p2 = size(b, 2)
                    corrupt_cols = falses(p2)
                    for j in 1:p2, i in 1:n2
                        # Propagate missings
                        # Corrupt all columns based on j
                        corrupt_cols[j] |= b.na[i, j]
                    end
                    res.na[:, corrupt_cols] = true
                end
            end)
            res
        end
    end
end

#
# Elementary functions
#
# XXX: The below should be revisited once we have a way to infer what
# the proper return type of an array should be.

# One-argument elementary functions that do something different for
# Complex
for f in (:(Base.abs), :(Base.abs2))
    @eval begin
        @dataarray_unary $(f) Complex T.parameters[1]
    end
end

# One-argument elementary functions that return the same type as their
# inputs
for f in (:(Base.abs), :(Base.abs2), :(Base.conj), :(Base.sign))
    @eval begin
        @dataarray_unary $(f) Number T
    end
end

# One-argument elementary functions that always return floating points
## Base
for f in (:(Base.acos), :(Base.acosh), :(Base.asin), :(Base.asinh), :(Base.atan), :(Base.atanh),
          :(Base.sin), :(Base.sinh), :(Base.cos), :(Base.cosh), :(Base.tan), :(Base.tanh),
          :(Base.exp), :(Base.exp2), :(Base.expm1), :(Base.log), :(Base.log10), :(Base.log1p),
          :(Base.log2), :(Base.exponent), :(Base.sqrt), :(Base.gamma), :(Base.lgamma))
    @eval begin
        @dataarray_unary $(f) AbstractFloat T
        @dataarray_unary $(f) Real Float64
    end
end
## SpecialFunctions (should be a conditional module when supported)
for f in (:(SpecialFunctions.digamma), :(SpecialFunctions.erf), :(SpecialFunctions.erfc))
    @eval begin
        @dataarray_unary $(f) AbstractFloat T
        @dataarray_unary $(f) Real Float64
    end
end

# Elementary functions that take varargs
for f in (:(Base.round), :(Base.ceil), :(Base.floor), :(Base.trunc))
    @eval begin
        # ambiguity
        @dataarray_unary $(f) Real T 1
        @dataarray_unary $(f) Real T 2
        @dataarray_unary $(f) Real T

        function $(f)(d::DataArray{T}, args::Integer...) where {T<:Real}
            data = similar(d.data)
            @bitenumerate d.na i na begin
                if !na
                    @inbounds data[i] = $(f)(d[i], args...)
                end
            end
            DataArray(data, copy(d.na))
        end
        function $(f)(adv::AbstractDataArray{T}, args::Integer...) where {T<:Real}
            res = similar(adv)
            for i = 1:length(adv)
                res[i] = ($f)(adv[i], args...)
            end
            res
        end
    end
end

#
# Comparison operators
#

# This is for performance only; the definition in Base is sufficient
# for AbstractDataArrays
function Base.isequal(a::DataArray, b::DataArray)
    if size(a) != size(b) || a.na != b.na
        return false
    end
    @bitenumerate a.na i na begin
        @inbounds if !na && !isequal(a.data[i], b.data[i])
            return false
        end
    end
    return true
end

# ambiguity
@swappable (==)(a::DataArray{Bool}, b::BitArray) =
    invoke(==, Tuple{DataArray,AbstractArray}, a, b)
@swappable (==)(a::DataArray, b::BitArray) =
    invoke(==, Tuple{DataArray,AbstractArray}, a, b)
@swappable (==)(a::AbstractDataArray{Bool}, b::BitArray) =
    invoke(==, Tuple{DataArray,AbstractArray}, a, b)
@swappable (==)(a::AbstractDataArray, b::BitArray) =
    invoke(==, Tuple{DataArray,AbstractArray}, a, b)

function (==)(a::DataArray, b::DataArray)
    size(a) == size(b) || return false
    adata = a.data
    bdata = b.data
    bchunks = b.na.chunks
    has_missing = false
    @bitenumerate a.na i amissing begin
        if amissing || Base.unsafe_bitgetindex(bchunks, i)
            has_missing = true
        else
            @inbounds adata[i] == bdata[i] || return false
        end
    end
    has_missing ? missing : true
end

@swappable function (==)(a::DataArray, b::AbstractArray)
    size(a) == size(b) || return false
    adata = a.data
    has_missing = false
    @bitenumerate a.na i amissing begin
        if amissing
            has_missing = true
        else
            @inbounds adata[i] == b[i] || return false
        end
    end
    has_missing ? missing : true
end

# ambiguity
@swappable (==)(a::DataArray, b::AbstractArray{>:Missing}) =
    invoke(==, Tuple{AbstractDataArray,AbstractArray}, a, b)

#
# Binary operators
#

# Define methods for UniformScaling. Otherwise we get ambiguity
# warnings...
if isdefined(Base, :UniformScaling)

function (+)(A::DataArray{TA,2},J::UniformScaling{TJ}) where {TA,TJ}
    n = LinAlg.checksquare(A)
    B = similar(A,promote_type(TA,TJ))
    copy!(B,A)
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] += J.λ
        end
    end
    B
end
(+)(J::UniformScaling,A::DataArray{TA,2}) where {TA} = A + J

function (-)(A::DataArray{TA,2},J::UniformScaling{TJ}) where {TA,TJ<:Number}
    n = LinAlg.checksquare(A)
    B = similar(A,promote_type(TA,TJ))
    copy!(B,A)
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] -= J.λ
        end
    end
    B
end
function (-)(J::UniformScaling{TJ},A::DataArray{TA,2}) where {TA,TJ<:Number}
    n = LinAlg.checksquare(A)
    B = -A
    @inbounds for i = 1:n
        if !B.na[i,i]
            B.data[i,i] += J.λ
        end
    end
    B
end

(+)(A::DataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(+, Tuple{AbstractArray{Union{Bool,Missing},2},UniformScaling{Bool}}, A, J)
(+)(J::UniformScaling{Bool},A::DataArray{Bool,2}) =
    invoke(+, Tuple{UniformScaling{Bool},AbstractArray{Union{Bool,Missing},2}}, J, A)
(-)(A::DataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(-, Tuple{AbstractArray{Union{Bool,Missing},2},UniformScaling{Bool}}, A, J)
(-)(J::UniformScaling{Bool},A::DataArray{Bool,2}) =
    invoke(-, Tuple{UniformScaling{Bool},AbstractArray{Union{Bool,Missing},2}}, J, A)

(+)(A::AbstractDataArray{TA,2},J::UniformScaling{TJ}) where {TA,TJ} =
    invoke(+, Tuple{AbstractArray{Union{TA,Missing},2},UniformScaling{TJ}}, A, J)
(+)(J::UniformScaling,A::AbstractDataArray{TA,2}) where {TA} =
    invoke(+, Tuple{UniformScaling,AbstractArray{Union{TA,Missing},2}}, J, A)
(-)(A::AbstractDataArray{TA,2},J::UniformScaling{TJ}) where {TA,TJ<:Number} =
    invoke(-, Tuple{AbstractArray{Union{TA,Missing},2},UniformScaling{TJ}}, A, J)
(-)(J::UniformScaling{TJ},A::AbstractDataArray{TA,2}) where {TA,TJ<:Number} =
    invoke(-, Tuple{UniformScaling{TJ},AbstractArray{Union{TA,Missing},2}}, J, A)

(+)(A::AbstractDataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(+, Tuple{AbstractArray{Union{Bool,Missing},2},UniformScaling{Bool}}, A, J)
(+)(J::UniformScaling{Bool},A::AbstractDataArray{Bool,2}) =
    invoke(+, Tuple{UniformScaling{Bool},AbstractArray{Union{Bool,Missing},2}}, J, A)
(-)(A::AbstractDataArray{Bool,2},J::UniformScaling{Bool}) =
    invoke(-, Tuple{AbstractArray{Union{Bool,Missing},2},UniformScaling{Bool}}, A, J)
(-)(J::UniformScaling{Bool},A::AbstractDataArray{Bool,2}) =
    invoke(-, Tuple{UniformScaling{Bool},AbstractArray{Union{Bool,Missing},2}}, J, A)

end # if isdefined(Base, :UniformScaling)

for f in (:(*), :(Base.div), :(Base.mod), :(Base.fld), :(Base.rem))
    @eval begin
        # Array with missing
        @swappable $(f)(::Missing, b::AbstractArray{T,N}) where {T, N} =
            DataArray(Array{T,N}(size(b)), trues(size(b)))

        # DataArray with scalar
        @dataarray_binary_scalar $f $f nothing true
    end
end

for f in (:(+), :(-))
    # Array with missing
    @eval @swappable $(f)(::Missing, b::AbstractArray{T,N}) where {T, N} =
        DataArray(Array{T,N}(size(b)), trues(size(b)))
end

for f in (:(+), :(-))
    @eval begin
        # Necessary to avoid ambiguity warnings
        @swappable ($f)(A::BitArray, B::AbstractDataArray{Bool}) = ($f)(Array(A), B)
        @swappable ($f)(A::BitArray, B::DataArray{Bool}) = ($f)(Array(A), B)

        @dataarray_binary_array $f $f
    end
end

# / is defined separately since it is not swappable
(/)(b::AbstractArray{T,N}, ::Missing) where {T,N} =
    DataArray(Array{T,N}(size(b)), trues(size(b)))
@dataarray_binary_scalar(/, /, nothing, false)

function Base.LinAlg.diff(dv::DataVector)
    n = length(dv)
    new_data = diff(dv.data)
    new_na = falses(n - 1)
    new_na[1] = dv.na[1]
    for i = 2:(n - 1)
        if dv.na[i]
            new_na[i - 1] = true
            new_na[i] = true
        end
    end
    new_na[n - 1] = new_na[n - 1] || dv.na[n]
    return DataArray(new_data, new_na)
end

# for f in cumulative_vector_operators
#     @eval function ($f)(dv::DataVector)
#         new_data = ($f)(dv.data)
#         new_na = falses(length(dv))
#         hitna = false
#         @bitenumerate dv.na i na begin
#             hitna |= na
#             if hitna
#                 new_na[i] = true
#             end
#         end
#         return DataArray(new_data, new_na)
#     end
# end
function Base.accumulate(f, dv::DataVector)
    new_data = accumulate(f, dv.data)
    new_na = falses(length(dv))
    hitna = false
    @bitenumerate dv.na i na begin
        hitna |= na
        if hitna
            new_na[i] = true
        end
    end
    return DataArray(new_data, new_na)
end
Base.cumsum(dv::DataVector) = accumulate(+, dv)
Base.cumprod(dv::DataVector)   = accumulate(*, dv)

for f in unary_vector_operators
    @eval ($f)(dv::DataVector) = any(dv.na) ? missing : ($f)(dv.data)
end

for f in binary_vector_operators
    @eval ($f)(dv1::DataVector, dv2::DataVector) =
            any(dv1.na) || any(dv2.na) ? missing : ($f)(dv1.data, dv2.data)
end

for f in (:(Base.minimum), :(Base.maximum), :(Base.prod), :(Base.sum),
          :(Base.mean), :(Base.median), :(Base.std), :(Base.var),
          :(Base.norm))
    colf = Symbol("col$(f)s")
    rowf = Symbol("row$(f)s")
    @eval begin
        function ($colf)(dm::AbstractDataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(p)
            for j in 1:p
                res[j] = ($f)(dm[:, j])
            end
            return res
        end
        function ($rowf)(dm::DataMatrix)
            n, p = nrow(dm), ncol(dm)
            res = datazeros(n)
            for i in 1:n
                res[i] = ($f)(DataArray(reshape(dm.data[i, :], p), reshape(dm.na[i, :], p)))
            end
            return res
        end
    end
end

#
# Boolean operators
#

function Base.all(dv::DataArray{Bool})
    data = dv.data
    hasmissings = false
    @bitenumerate dv.na i na begin
        if !na
            data[i] || return false
        else
            hasmissings = true
        end
    end
    hasmissings ? missing : true
end

function Base.all(dv::AbstractDataArray{Bool})
    hasmissings = false
    for i in 1:length(dv)
        x = dv[i]
        if !ismissing(x)
            x || return false
        else
            hasmissings = true
        end
    end
    hasmissings ? missing : true
end

function Base.any(dv::DataArray{Bool})
    hasmissings = false
    @bitenumerate dv.na i na begin
        if !na
            if dv.data[i]
                return true
            end
        else
            hasmissings = true
        end
    end
    hasmissings ? missing : false
end

function Base.any(dv::AbstractDataArray{Bool})
    hasmissings = false
    for i in 1:length(dv)
        if !ismissing(dv[i])
            if dv[i]
                return true
            end
        else
            hasmissings = true
        end
    end
    hasmissings ? missing : false
end

function rle(v::AbstractVector{T}) where T
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = similar(v, n)
    total_values = 1
    lengths = Array(Int16, n)
    total_lengths = 1
    for i in 2:n
        if v[i] == current_value
            current_length += 1
        else
            values[total_values] = current_value
            total_values += 1
            lengths[total_lengths] = current_length
            total_lengths += 1
            current_value = v[i]
            current_length = 1
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

function rle(v::AbstractDataVector{T}) where T
    n = length(v)
    current_value = v[1]
    current_length = 1
    values = DataArray(T, n)
    total_values = 1
    lengths = Vector{Int16}(uninitialized, n)
    total_lengths = 1
    for i in 2:n
        if ismissing(v[i]) || ismissing(current_value)
            if ismissing(v[i]) && ismissing(current_value)
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        else
            if v[i] == current_value
                current_length += 1
            else
                values[total_values] = current_value
                total_values += 1
                lengths[total_lengths] = current_length
                total_lengths += 1
                current_value = v[i]
                current_length = 1
            end
        end
    end
    values[total_values] = current_value
    lengths[total_lengths] = current_length
    return (values[1:total_values], lengths[1:total_lengths])
end

## inverse run-length encoding
function inverse_rle(values::AbstractVector{T},
                     lengths::Vector{I}) where {T, I <: Integer}
    total_n = Int(sum(lengths))
    pos = 0
    res = similar(values, total_n)
    n = length(values)
    for i in 1:n
        v = values[i]
        l = lengths[i]
        for j in 1:l
            pos += 1
            res[pos] = v
        end
    end
    return res
end


text:
include("derivative_example.jl")
include("example.jl")

ext: jl


text:
function test_fltype()
  Tf = Float16
  Tc = Complex{Tf}
  Tf == fltype(Tf) && Tf == fltype(Tc) && Tf == fltype((Tc, Tf)) && Tf == fltype((Tf, Tc))
end

function test_cxtype()
  Tf = Float16
  Tc = Complex{Tf}
  Tc == cxtype(Tf) && Tc == cxtype(Tc)
end

function test_innereltype(T=Complex{Float32})
  a = [zeros(T, 3), zeros(T, 2, 5)]
  T == innereltype(a)
end

function test_superzeros()
  T = Float64
  dims1 = (2, 3)
  dims2a = (7,)
  dims2b = dims1
  dims2 = (dims2a, dims2b)

  a1 = rand(T, dims1)
  a2 = [(1+im)*rand(T, dims2a), (2+im)*rand(T, dims2b)]

  a1a = zeros(T, dims1)
  a2a = [zeros(Complex{T}, dims2a), zeros(Complex{T}, dims2b)]

  a1z = superzeros(a1)
  a2z = superzeros(a2)
  a1d = superzeros(dims1)
  a2d = superzeros(Complex{T}, dims2)

  @superzeros innereltype(a1) a1 a1ma dum
  @superzeros innereltype(a2) a2 a2ma dum
  @superzeros innereltype(a1) dims1 a1md dum
  @superzeros innereltype(a2) dims2 a2md dum

  ( a1a == a1z && a1a == a1d && a1a == a1ma && a1a == a1md &&
    a2a == a2z && a2a == a2d && a2a == a2ma && a2a == a2md )
end

function test_supertuplezeros(; T1=Float64, T2=Complex{Float64}, dims1=(1,), dims2=(3, 3))
  T = (T1, T2)
  dims = (dims1, dims2)
  a = superzeros(T, dims)
  ( eltype(a[1]) == T1 && eltype(a[2]) == T2 &&
    size(a[1]) == dims1 && size(a[2]) == dims2 )
end

# This test could use some further work.
function test_radialspectrum(dev::Device, n, ahkl, ahρ; debug=false, atol=0.1, rfft=false)
  grid = TwoDGrid(dev, n, 2π)
  if rfft==true
    ah = @. ahkl(grid.kr, grid.l)
  else
    ah = @. ahkl(grid.k, grid.l)
  end
  CUDA.@allowscalar ah[1, 1] = 0.0

  ρ, ahρ_estimate = FourierFlows.radialspectrum(ah, grid; refinement=16)

  if debug
    println(sum(ahρ.(ρ) - ahρ_estimate))
    return ah, ahρ_estimate, ahρ.(ρ)
  else
    normalizeddiff = sum(abs.(ahρ.(ρ) - ahρ_estimate)) / length(ρ)
    return isapprox(normalizeddiff, 0.0, atol=atol)
  end
end

function integralsquare(func, grid::OneDGrid)
    sum(abs2.(func))*grid.dx
end

function integralsquare(func, grid::TwoDGrid)
    sum(abs2.(func))*grid.dx*grid.dy
end

function test_parsevalsum(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum = FourierFlows.parsevalsum(abs2.(funch), grid)
  isapprox(integral, parsevalsum; rtol=rtol_utils)
end

function test_parsevalsum2(func, grid; realvalued=true)
  # Compute integral in physical space
  integral = integralsquare(func, grid)
  # Compute integral in wavenumber space using Parseval's theorem
  if realvalued==true;      funch = rfft(func)
  elseif realvalued==false; funch = fft(func)
  end
  parsevalsum2 = FourierFlows.parsevalsum2(funch, grid)
  isapprox(integral, parsevalsum2; rtol=rtol_utils)
end

"""
Compute the jacobian J(a, b) and compare the result with `analytic`. Use `atol` for the comparison
to ensure validity when `analytic=0`.
"""
test_jacobian(a, b, analytic, grid) = isapprox(FourierFlows.jacobian(a, b, grid), analytic;
                                               atol=grid.nx*grid.ny*10*eps())

"""
Test the zeros macro.
"""
function test_zeros(T=Float64, dims=(13, 45))
  a1, b1 = zeros(T, dims), zeros(T, dims)
  @zeros T dims a2 b2
  a1 == a2 && b1 == b2
end

"""
Test supersize().
"""
function test_supersize()
  a = rand(16, 16)
  dimsa = size(a)

  b = [rand(1), rand(3, 34)]
  dimsb = ((1,), (3, 34))

  dimsa == supersize(a) && dimsb == supersize(b)
end

"""
Test arraytype().
"""
function test_arraytype(dev::Device)
  dev==CPU() ? ArrayType(dev)==Array : ArrayType(dev)==CuArray
end

function test_arraytypeTdim(dev::Device, T=Float64, dim=1)
  dev==CPU() ? ArrayType(dev, T, dim)<:Array{T, dim} : ArrayType(dev, T, dim)<:CuArray{T, dim}
end

function test_ongrid(dev::Device)
  nx, ny, nz = 6, 8, 10
  Lx, Ly, Lz = 2π, 2.0, 3.0

  g₁ = OneDGrid(dev, nx, Lx)
  X₁ = ArrayType(dev)(g₁.x)
  f₁(x) = x^2

  g₂ = TwoDGrid(dev, nx, Lx, ny, Ly)
  X₂, Y₂ = gridpoints(g₂)
  f₂(x, y) = x^2 - y^3

  g₃ = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)
  X₃, Y₃, Z₃ = gridpoints(g₃)
  f₃(x, y, z) = x^2 - y^3 + sin(z)

  return (FourierFlows.on_grid(f₁, g₁) ≈ f₁.(X₁) && FourierFlows.on_grid(f₂, g₂) ≈ f₂.(X₂, Y₂) && FourierFlows.on_grid(f₃, g₃) ≈ f₃.(X₃, Y₃, Z₃))
end


text:
using BlockMatching: to_cartesian
using BlockMatching: topk, topk!

@testset "to_cartesian" begin
    @test @inferred to_cartesian(2, 1, 2) == (CartesianIndex(1, 1), CartesianIndex(2, 2))
    @test @inferred to_cartesian(2, (1, 2), (3, 4)) == (CartesianIndex(1, 2), CartesianIndex(3, 4))
    @test @inferred to_cartesian(2, CartesianIndex(1, 2), CartesianIndex(3, 4)) == (CartesianIndex(1, 2), CartesianIndex(3, 4))
    @test @inferred to_cartesian(2, CartesianIndex(1, 3), 1, (4, 4)) == (CartesianIndex(1, 3), CartesianIndex(1, 1), CartesianIndex(4, 4))

    @test_throws MethodError to_cartesian(3, (1, 2), 3)
    @test_throws MethodError to_cartesian(3, CartesianIndex(1, 2), 3)
end

@testset "topk" begin
    K, N = 20, 200
    for T in (Int, Float32)
        data = T.(rand(1:N, N))
        data_sorted = sort(data)

        vals, inds = topk(data, K)
        @test isempty(setdiff(data_sorted[end-K+1:end], vals))
        @test data[inds] == vals

        inds = fill(0, K)
        vals = fill(typemin(eltype(data)), K)
        topk!(vals, inds, data, Base.Order.Forward)
        @test isempty(setdiff(data_sorted[end-K+1:end], vals))
        @test data[inds] == vals

        vals, inds = topk(data, K, Base.Order.Reverse)
        @test isempty(setdiff(data_sorted[1:K], vals))
        @test data[inds] == vals

        inds = fill(0, K)
        vals = fill(typemax(eltype(data)), K)
        topk!(vals, inds, data, Base.Order.Reverse)
        @test isempty(setdiff(data_sorted[1:K], vals))
        @test data[inds] == vals
    end
end


text:
using Primes

struct FactoredPoly
    factors::Array{Pair{Poly,Int}}

    FactoredPoly() = new(Array{Pair{Poly,Int}}[])
end

function add_factor!(f::FactoredPoly, p, power=1)
    push!(f.factors, p => power)
end

sym(f::FactoredPoly) = prod(map(v -> sym(first(v))^last(v), f.factors); init=1)
poly(f::FactoredPoly) = prod(map(v -> first(v)^last(v), f.factors); init=1)

Base.show(io::IO, f::FactoredPoly) = print(io, sym(f))

###########################################################################################

function split_poly(p::Poly)
    l, p = integer_poly(p)
    cont = gcd(extract_coef(p)...)
    pp = p ÷ cont
    cont * l, pp
end

# Lagrange interpolation algorithm for rational polynomials
function interp_poly(x, xs, ys)
    length(xs) != length(ys) && error("interp_poly needs the same number of x and y values")
    n = length(xs)
    Ps = 0
    for i = 1:n
        A = *([(x-xs[j]) for j=1:n if j!=i]...)
        B = *([(xs[i]-xs[j]) for j=1:n if j!=i]...)
        Ps += ys[i] *  A // B
    end
    return poly(Rational, Ps, x)
end

function pre_factor(p::Poly)
    l, p = integer_poly(p)
    x = p.x
    r = copy(p)
    f = Poly[]

    a₀ = p[0]
    aₙ = p[degree(p)]

    for root in find_roots(p, x)[1]
        ρ = Rational(root)
        println(ρ)
        if a₀ % numerator(ρ) == 0 && aₙ % denominator(ρ) == 0
            q = gen_poly(x, -numerator(ρ), denominator(ρ))
            r = r / q
            push!(f, q)
        end
    end

    return r, f, l
end

function list_divisors(x)
    f = factor(x)
    l = [1]

    for (p, k) in f
        r = Int[]
        for i = 0:k
            append!(r, p^i .* l)
        end
        l = r
    end
    l
end

# note k is 1-based
function generate_cross_lists(D)
    n = *(length.(D)...)
    L = zeros(Int, (length(D), n))

    for k = 1:n
        j = k - 1
        for i = 1:length(D)
            n = length(D[i])
            L[i,k] = D[i][(j % n) + 1]
            j = j ÷ n
        end
    end
    return L
end

function assemble_factors(f, r, l)
    if !iszero(r - 1)
        # push!(f, r)
        add_factor!(f, r)
    end
    if !isone(l)
        # push!(f, Poly(l))
        add_factor!(f, Poly(l))
    end
    f
end

function factor_schubert_kronecker(p::Poly)
    l, p = integer_poly(p)
    x = p.x
    r = copy(p)
    # f = Poly[]
    f = FactoredPoly()
    A = Int[]
    D = Array{Int}[]

    i = 0
    while length(D) <= degree(r)÷2
        a = ((i+1)÷2) * (isodd(i) ? -1 : +1)
        u = convert(Int, p(a))

        if iszero(u)
            while iszero(r % (x-a))
                r = r / (x-a)
                # push!(f, poly(x-a, x))
                add_factor!(f, poly(x-a, x))
            end
        else
            push!(A, a)
            divs = list_divisors(abs(u))
            push!(D, [divs; -divs])
        end
        i += 1
    end

    for d = 1:degree(r)÷2
        L = generate_cross_lists(D[1:d+1])
        xs = A[1:d+1]
        for i = 1:size(L,2)
            ys = L[:,i]
            q = interp_poly(x, xs, ys)

            while degree(q) > 0 && iszero(r % q)
                r = r / q
                # push!(f, q)
                add_factor!(f, q)
                if degree(r) < 2d
                    return assemble_factors(f, r, l)
                end
            end
        end
        if degree(r) < 2*(d+1)
            return assemble_factors(f, r, l)
        end
    end
    assemble_factors(f, r, l)
end

"""
    square-free decomposition of p
    Yun's algorithm
"""
function decompose(p::Poly)
    p = rationalize(p)
    f = FactoredPoly()
    p′ = derivative(p)
    g = gcd(p, p′)
    r = p / g
    i = 1
    c = 1

    while !iszero(g - 1) && r isa Poly
        s = gcd(g, r)
        p = r / s
        if !iszero(p - 1)
            c *= leading(p)
            #push!(f, (p / leading(p), i))
            add_factor!(f, p / leading(p), i)
        end
        i += 1
        r = s
        g = g / s
    end

    if !iszero(r - 1)
        c *= leading(r)
        # push!(f, (r / leading(r), i))
        add_factor!(f, r / leading(r), i)
    end

    if c != 1
        push!(f, (Poly(c), 0))
    end

    return f
end

function recompose(f)
    y = Poly(1)

    for (k, p) in f
        y = y * p^k
    end
    y
end

function Primes.factor(p::Poly; verbose=false)
    _, pᵢ = integer_poly(p)
    c, pₘ = to_monic(pᵢ)
    f = decompose(pₘ)
    verbose && printstyled('\t', f, '\n'; color=:red)
    h = FactoredPoly()
    p₁ = poly(Rational, 1, p.x)

    for w in f.factors
        v, k = first(w), last(w)
        if degree(v) > 0
            for w in factor_schubert_kronecker(v).factors
                w₁ = prim(from_monic(first(w), c))
                # push!(h, (w₁, k))
                add_factor!(h, w₁, k)
                p₁ = p₁ * w₁^k
            end
        end
    end

    q, r = divide_poly(p, p₁)

    if !iszero(r)
        @warn "incorrect factorization"
        println(p₁)
    end

    if !iszero(q - 1)
        add_factor!(h, q, 0)
        # push!(h, (q, 0))
    end

    return h
end


text:
using GLFW, ModernGL

# whenever the window size changed (by OS or user resize) this callback function executes
function framebuffer_size_callback(window::GLFW.Window, width::Cint, height::Cint)
    # make sure the viewport matches the new window dimensions; note that width and
    # height will be significantly larger than specified on retina displays.
	glViewport(0, 0, width, height)
end

# process all input: query GLFW whether relevant keys are pressed/released this frame and react accordingly
function processInput(window::GLFW.Window)
    GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS && GLFW.SetWindowShouldClose(window, true)
end

# misc. config
const SCR_WIDTH = 800
const SCR_HEIGHT = 600

const vertexShaderSource = """
	#version 330 core
	layout (location = 0) in vec3 aPos;
	void main(void)
	{
	    gl_Position = vec4(aPos.x, aPos.y, aPos.z, 1.0);
	}"""

const fragmentShaderSource = """
	#version 330 core
	out vec4 FragColor;
	void main(void)
	{
	    FragColor = vec4(1.0f, 0.5f, 0.2f, 1.0f);
	}"""

GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 3)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)

# create window
window = GLFW.CreateWindow(SCR_WIDTH, SCR_HEIGHT, "LearnOpenGL")
window == C_NULL && error("Failed to create GLFW window.")
GLFW.MakeContextCurrent(window)
GLFW.SetFramebufferSizeCallback(window, framebuffer_size_callback)

# build and compile shader program
vertexShader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vertexShader, 1, Ptr{GLchar}[pointer(vertexShaderSource)], C_NULL)
glCompileShader(vertexShader)
successRef = Ref{GLint}(-1)
glGetShaderiv(vertexShader, GL_COMPILE_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetShaderInfoLog(vertexShader, 512, C_NULL, infoLog)
    error("ERROR::SHADER::VERTEX::COMPILATION_FAILED\n", String(infoLog))
end

fragmentShader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(fragmentShader, 1, Ptr{GLchar}[pointer(fragmentShaderSource)], C_NULL)
glCompileShader(fragmentShader)
glGetShaderiv(vertexShader, GL_COMPILE_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetShaderInfoLog(fragmentShader, 512, C_NULL, infoLog)
    error("ERROR::SHADER::FRAGMENT::COMPILATION_FAILED\n", String(infoLog))
end

# link shaders
shaderProgram = glCreateProgram()
glAttachShader(shaderProgram, vertexShader)
glAttachShader(shaderProgram, fragmentShader)
glLinkProgram(shaderProgram)
glGetProgramiv(shaderProgram, GL_LINK_STATUS, successRef)
if successRef[] != GL_TRUE
    infoLog = Vector{GLchar}(512)
    glGetProgramInfoLog(shaderProgram, 512, C_NULL, infoLog)
    error("ERROR::SHADER::PROGRAM::LINKING_FAILED\n", String(infoLog))
end
glDeleteShader(vertexShader)
glDeleteShader(fragmentShader)

# set up vertex data (and buffer(s)) and configure vertex attributes
vertices = GLfloat[-0.5,  0.5, 0.0,   # left
                    0.5, -0.5, 0.0,   # right
                    0.0,  0.5, 0.0]   # top
vboRef = Ref{GLuint}(0)
vaoRef = Ref{GLuint}(0)
glGenVertexArrays(1, vaoRef)
glGenBuffers(1, vboRef)
vao = vaoRef[]
vbo = vboRef[]
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), Ptr{Cvoid}(0))
glEnableVertexAttribArray(0)

# note that this is allowed, the call to glVertexAttribPointer registered VBO as
# the vertex attribute's bound vertex buffer object so afterwards we can safely unbind
glBindBuffer(GL_ARRAY_BUFFER, 0);

# You can unbind the VAO afterwards so other VAO calls won't accidentally modify this VAO,
# but this rarely happens. Modifying other VAOs requires a call to glBindVertexArray anyways
# so we generally don't unbind VAOs (nor VBOs) when it's not directly necessary.
glBindVertexArray(0)

# uncomment this call to draw in wireframe polygons.
# glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)

# render loop
while !GLFW.WindowShouldClose(window)
    processInput(window)

    # render
    glClearColor(0.2, 0.3, 0.3, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)

	# draw
	glUseProgram(shaderProgram)
    glBindVertexArray(vao)
	glDrawArrays(GL_TRIANGLES, 0, 3)

    # swap buffers and poll IO events
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

# optional: de-allocate all resources once they've outlived their purpose
glDeleteVertexArrays(1, vaoRef)
glDeleteBuffers(1, vboRef)

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using SparseArrays

"""
    tgw(
        model::CAGD.Model, dim::Int; atol = 1e-7
    )::Tuple{Lar.ChainOp, Array{Array{Int,1},1}}

Compute the Topological Gift Wrapping on `model` w.r.t. `dim`-cell.

It returns the gift wrapped *chain* cycles operator along with the biconnected
components w.r.t `dim`-cells.
"""
function tgw(model, dim; atol = 1e-7)

    function select_sigma()
        for i = 1 : lo_num  if visited[i] == 1  return i  end  end
        for i = 1 : lo_num  if visited[i] == 0  return i  end  end
        return -1
    end

    function get_closer_petal(τ, stem, τ_sign)
        if !isdefined(ord_angles, Int(τ))
            ord_angles[τ] = CAGD.eval_ord_angle(model, dim - 1, τ, atol=atol)
        end
        stem_idx = findfirst(isequal(stem), ord_angles[τ])
        return ord_angles[τ][mod(stem_idx + τ_sign - 1, length(ord_angles[τ])) + 1]
    end

    # Prepearing tgw model
    chain = SparseArrays.spzeros(Int8, size(model, dim, 2), 0)
    lo_cls = model.T[dim - 1]
    lo_num, llo_num = size(lo_cls)

    # Preprocessing angles between d-1 cells w.r.t. d-2 cells
    # ord_angles has for every d-2 cell the circular ordering between d-1 cells
    ord_angles = Array{Array{Int64, 1}, 1}(undef, llo_num)

    # Either 0 (not visited), 1 (once visited) or 2 (twice visited)
    visited = zeros(Int8, lo_num)
    # If visited[σ] == 1 then other_sign[σ] is the sign for the second visit
    other_sign = zeros(Int8, lo_num)

    # Components accumulator
    components = Array{Array{Int64, 1},1}()
    comp_idx = 0
    cell_idx = 0

    # Until all dim-1 cells have been visited twice, continue
    while (σ = select_sigma()) > 0
        # Build new 3-cell
        cell_idx += 1
        c = SparseArrays.spzeros(Int8, 1, lo_num)
        if visited[σ] == 0

            # Topological invariants Check #6: Euler Characteristic
            if comp_idx > 0
                euler = [components[comp_idx]]
                push!(euler, ∪([chain[:, c].nzind for c in euler[end]]...))
                for d = dim-1 : -1 : 1
                    push!(euler, CAGD.getModelLoCell(model, d, euler[end]))
                end
                @assert sum([(-1)^(dim+1-i)*length(euler[i]) for i = 1 : dim+1]) == 0
            end

            # Set sign in representation
            c[σ] = 1
            other_sign[σ] = -1
            # Set start to a new component
            push!(components, [])
            comp_idx += 1
        else
            c[σ] = other_sign[σ]
        end
        push!(components[comp_idx], cell_idx)
        visited[σ] += 1

        # Build other faces via corolla method
        # boundary = d-2 exposed cells of currently analysed d cell
        while (boundary = (dropzeros(c * model.T[dim - 1])[1, :])).nzind != []
            # corolla are the next d-1 cells to be added
            corolla = zeros(Int8, lo_num)
            for τ in boundary.nzind
                # petals are incident d-1 cells over τ d-2 cell
                petals = model.T[dim - 1][:, τ]
                # stem is the cell generating the cell generating the petals
                #TODO it should be one only ?!?
                stem = intersect(c[1, :].nzind, petals.nzind)[1]
                # τ is travelled according the stem direction
                τ_sign = sign(boundary[τ]) # * sign(c[stem]) già fatto dalla molt
                petal = get_closer_petal(τ, stem, τ_sign)          # -sign(...) ??
                # the new petal and the stem must travel τ in opposit direction
                if model.T[dim - 1][stem, τ] == - model.T[dim - 1][petal, τ]
                    petal_sign = c[stem]
                elseif model.T[dim - 1][stem, τ] == model.T[dim - 1][petal, τ]
                    petal_sign = -c[stem]
                else
                    throw(ArgumentError("Something Weird Happened"))
                end

                # If petal was visited earlyer (is in c) then it must be visited discordely
                @assert visited[petal] < 2 "Petal $petal visited more than twice"
                @assert visited[petal] == 0 || other_sign[petal] == petal_sign "Incoherent Petal sign τ = $τ, petal = $petal"
                # If petal already is in the corolla, then it should be with the same sign
                # (otherwise it is visite twice)
                if (corolla[petal] != 0) && (corolla[petal] != petal_sign)
                    @show "Cell $petal is visited twice during $dim-d TGW."
                end

                # Update Corolla if petal not visited yet
                if corolla[petal] == 0                                                              #TODO Change for twice visiting
                    corolla[petal] += petal_sign
                end

            end

            # Add corolla to c and update visited, other_sign, componentsFaces coherently
            c += corolla'
            other_sign -= corolla
            visited += abs.(corolla)

        end

        chain = [chain c']

    end

    return chain, components
end


text:
"""
arar(y::TimeArray, h::Int, freq::DataType, max_lag::Int)

Forecasting using ARAR algorithm.

Return A matrix of forecast values and prediction intervals

# Arguments
- `y::TimeArray`: An TimeArray with only one value column.
- `h::Int`: Forecast horizon as an integer.
- `freq::DataType`: A DataType from Dates, e.g. Dates.Day or Dates.Month.
- `max_lag::Int`: An integer (>= 26) to specify maximum number of lags for sample autocovariance.

# Examples
```julia-repl
julia> arar(data, 12, Month)

```
"""
function arar(y::TimeArray, h::Int, freq::DataType, max_lag::Int=40)
  future_dates = range(maximum(timestamp(y)) + freq(1); step=freq(1), length=h)
  y = dropdims(values(y), dims = 2)
  Y = y
  @assert length(y) >= max_lag "Series is too short for arar"
  Ψ = [1]

  for k in collect(1:3)
    n = length(y)
    # memory shortening
    ϕ = map(τ -> (sum(y[(τ + 1):n] .*
    y[1:(n - τ)])/sum(y[1:(n - τ)].^2)), 1:15)
    err = map(τ -> (sum((y[(τ + 1):n] -
    ϕ[τ]*y[1:(n - τ)]).^2)/sum(y[(τ + 1):n].^2)), 1:15)
    τ = argmin(err)
    if err[τ] <= 8/n || (ϕ[τ] >= 0.93 && τ > 2)
      y = y[(τ + 1):n] - ϕ[τ] * y[1:(n - τ)]
      Ψ = [Ψ;zeros(τ)] .- ϕ[τ] .* [zeros(τ); Ψ]
    elseif ϕ[τ] >= 0.93
      A = zeros(2, 2)
      A[1, 1] = sum(y[2:(n - 1)].^2)
      A[1, 2] = sum(y[1:(n - 2)] .* y[2:(n - 1)])
      A[2, 1] = sum(y[1:(n - 2)] .* y[2:(n - 1)])
      A[2, 2] = sum(y[1:(n - 2)].^2)
      b = (sum(y[3:n] .* y[2:(n - 1)]), sum(y[3:n] .* y[1:(n - 2)]))
      ϕ = ((A' * A) ^ -1) * A'*collect(b)
      y = y[3:n] - ϕ[1] * y[2:(n - 1)] - ϕ[2] * y[1:(n - 2)]
      Ψ = (Ψ, [0], [0]) .- (ϕ[1] .* (0, Ψ, 0) ).- ϕ[2] .* (0, 0, Ψ)
    else
      break

    end
  end

S = y
Sbar = mean(S)
X = S .- Sbar
n = length(X)
xbar = Statistics.mean(X)
gamma = map(i -> sum((X[1:(n - i)] .- xbar) .* (X[(i + 1):n] .- xbar))/n, 0:max_lag)
y = Y
A = reshape(repeat([gamma[1]], 16), (4,4))
b = zeros(4)
best_σ2 = Inf
m = 26
best_lag = []
best_phi = []

for i in 2:(m-2)
    for j in (i + 1):(m-1)
        for k in (j + 1):m
            A[1, 2] = A[2, 1] = gamma[i]
            A[1, 3] = A[3, 1] = gamma[j]
            A[2, 3] = A[3, 2] = gamma[j - i + 1]
            A[1, 4] = A[4, 1] = gamma[k]
            A[2, 4] = A[4, 2] = gamma[k - i + 1]
            A[3, 4] = A[4, 3] = gamma[k - j + 1]
            b[1] = gamma[2]
            b[2] = gamma[i + 1]
            b[3] = gamma[j + 1]
            b[4] = gamma[k + 1]

            ϕ = ((A' * A) ^ -1) * A'*collect(b)

            σ2 = gamma[1] .- ϕ[1] .* gamma[2] .- ϕ[2] .*
            gamma[i + 1] .- ϕ[3] .* gamma[j + 1] .- ϕ[4] .* gamma[k + 1]

            if σ2 < best_σ2
                best_σ2 = σ2
                best_phi = ϕ
                best_lag = (1, i, j, k)
            end
        end

    end

end

i = best_lag[2]
j = best_lag[3]
k = best_lag[4]
ϕ = best_phi
σ2 = best_σ2

xi = [Ψ; zeros(k)] .-
ϕ[1] * [[0]; Ψ; zeros(k-1)] .-
ϕ[2] * [zeros(i); Ψ; zeros(k - i)] .-
ϕ[3] * [zeros(j); Ψ; zeros(k - j)] .-
ϕ[4] * [zeros(k); Ψ]

n = length(y)
k = length(xi)
y = [y; zeros(h)]
c = (1 - sum(ϕ)) .* Sbar

meanfc = map(i-> y[n + i] = - sum(xi[2:k] .* y[n + i + 1 .- (2:k)]) .+ c, 1:h)
y = y[1:n]
if h > k
  xi = append!(xi, zeros(h - k))
end
τ = zeros(h)
τ[1] = 1

if h > 1
  for j in 1:(h-1)
    τ[j + 1] = -sum((τ[1:j] .* xi[sort(((1:j) .+ 1), rev = true)]))
  end

end

se = Statistics.sqrt!(σ2 .* map(j -> sum(τ[1:j].^2), 1:h))

data = (datetime = future_dates,
Point_Forecast = meanfc,
Upper95 = meanfc + 1.96 .* se,
Upper80 = meanfc + 1.28 .* se,
Lower95 = meanfc - 1.96 .* se,
Lower80 = meanfc - 1.28 .* se)
out = TimeArray(data; timestamp = :datetime)

return out

end


text:
include("guiding_center_4d_settings.jl")

set_config(:quasi_newton_refactorize, 1)
set_config(:nls_atol, 1E-20)
set_config(:nls_rtol, 2eps())
set_config(:nls_nmax, 20)


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

add(x, y) = x + y
emit(add, Float64, Float64)

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !jlir<"typeof(Main.add)">, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = addf %arg1, %arg2 : f64
# CHECK-NEXT:     return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.fadd %arg1, %arg2  : f64
# CHECK-NEXT:     llvm.return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.add), Float64, Float64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: f64, %arg2: f64) -> f64 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.add), Float64, Float64}"(%arg0, %arg1, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, f64, f64) -> f64
# CHECK-NEXT:     llvm.return %0 : f64
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct HydroDispatch <: HydroGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        tech::TechHydro
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`, validation range: tech.reactivepowerlimits, action if invalid: warn
- `tech::TechHydro`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HydroDispatch <: HydroGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechHydro
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function HydroDispatch(name, available, bus, activepower, reactivepower, tech, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    HydroDispatch(name, available, bus, activepower, reactivepower, tech, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function HydroDispatch(; name, available, bus, activepower, reactivepower, tech, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    HydroDispatch(name, available, bus, activepower, reactivepower, tech, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function HydroDispatch(::Nothing)
    HydroDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechHydro(nothing),
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::HydroDispatch) = value.name
"""Get HydroDispatch available."""
get_available(value::HydroDispatch) = value.available
"""Get HydroDispatch bus."""
get_bus(value::HydroDispatch) = value.bus
"""Get HydroDispatch activepower."""
get_activepower(value::HydroDispatch) = value.activepower
"""Get HydroDispatch reactivepower."""
get_reactivepower(value::HydroDispatch) = value.reactivepower
"""Get HydroDispatch tech."""
get_tech(value::HydroDispatch) = value.tech
"""Get HydroDispatch services."""
get_services(value::HydroDispatch) = value.services
"""Get HydroDispatch ext."""
get_ext(value::HydroDispatch) = value.ext

InfrastructureSystems.get_forecasts(value::HydroDispatch) = value.forecasts
"""Get HydroDispatch internal."""
get_internal(value::HydroDispatch) = value.internal


text:
struct Network
    inputs::Vector{String}
    outputs::Vector{String}
    layers::Vector{Layer}
end

function Network(m::Chain, inputs::Vector{String}, outputs::Vector{String}, outputname)
    if isempty(inputs)
        inputs = repeat([""], findinputs(m))
    end
    if isempty(outputs)
        outputs = repeat([""], findoutputs(m))
    end
    layers = vcat(Layer.(m), Layer(outputname, length(outputs), 0, []))
    if layers[1].inputs == 0
        layers[1] = Layer(layers[1].name, length(inputs))
    end
    for i in eachindex(layers)[2:end - 1]
        if layers[i].inputs == 0
            layers[i] = Layer(layers[i].name, layers[i - 1].outputs)
        end
    end
    for i in eachindex(layers)[1:end - 1]
        o = layers[i].outputs
        i = layers[i + 1].inputs
        i != o && error(
            DimensionMismatch("layer A has $o outputs, layer B has $i inputs")
        )
    end
    Network(inputs, outputs, layers)
end

function findinputs(m::Chain)
    i = 1
    inputs = 0
    while inputs == 0
        inputs = Layer(m[i]).inputs
        i += 1
        if i > length(m)
            inputs = 1
        end
    end
    inputs
end

function findoutputs(m::Chain)
    i = length(m)
    outputs = 0
    while outputs == 0
        outputs = Layer(m[i]).outputs
        i -= 1
        if i < 1
            outputs = 1
        end
    end
    outputs
end


text:
function distance(p1::AbstractPoint, p2::AbstractPoint)
    return sqrt( (p2.x -p1.x)^2 + (p2.y -p1.y)^2 )
    # return norm(p2 - p1)
end

function distance(line::Line, point::AbstractPoint)
    e12 = normalize!( vector(line) )
    v1p = point - line.p1
    val = dot(v1p, v1p) - dot(e12, v1p)^2
    if val < zero(Float64)
        val = zero(Float64)
    end
    return sqrt( val )
end

function  distance(line_p1::AbstractPoint, line_p2::AbstractPoint, point::AbstractPoint)
    e12 = normalize!( line_p2 - line_p1 )
    v1p = point - line_p1
    val = dot(v1p, v1p) - dot(e12, v1p)^2
    if val < zero(Float64)
        val = zero(Float64)
    end
    return sqrt( val )
end

# ---------------------------------------------------------------------------- #
# Methods for mesh generation
# ---------------------------------------------------------------------------- #

function signed_distance_function(circle::T, p) where {T <: AbstractCircle}
    d = sqrt( (p[1] - circle.center.x)^2 + (p[2] - circle.center.y)^2 ) - circle.r
    return d
end

function signed_distance_function(rectangle::T, p) where {T <: AbstractRectangle}
    px = p[1]
    py = p[2]

    t1 = min(
        -px + rectangle.p2.x,
        px - rectangle.p1.x
    )

    t2 = min(
        -py + rectangle.p2.y,
        py - rectangle.p1.y
    )

    d = -min(t1, t2)
    return d
end

function union_signed_distance_function(d1, d2)
    d = min.(d1, d2)
    return d
end
function intersection_signed_distance_function(d1, d2)
    d = max.(d1, d2)
    return d
end
function difference_signed_distance_function(d1, d2)
    d = max.( d1, -d2 )
    return d
end

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

function point_source_node_spacing(p, src_point, h_0, h_inf, spread_factor = 1)

	dx = p[1] - src_point[1];
	dy = p[2] - src_point[2];

	norm_dp_sq 	= dx.^2 + dy.^2;
	arg 		= .- spread_factor .* norm_dp_sq;
	h 			= (h_0 .- h_inf) .* exp.(arg) .+ h_inf;

    return h
end

function uniform_node_spacing(p)
	# Uniform h(x,y) distribution
	h = 1.0
    return h
end

function left_to_right_node_spacing(p)
	# 1 + x
	# 1 + x of p
	h = 1.0 + max(p[1], 0.0);
    return h
end

function bottom_to_top_node_spacing(p)
	# 1 + y
	# 1 + y of p
	h = 1.0 + max(p[2], 0);
    return h
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function sesam3C(;name, k_N=60.0)
    # @parameters t [unit = uT]
    # D = Differential(t)

    # @parameters ϵ ϵ_tvr κ_E a_E [unit = uT^-1] m [unit = uT^-1] τ [unit = uT^-1]
    # @parameters k_L [unit = uT^-1] k_R [unit = uT^-1] k_mN [unit = uQ*uT^-1]
    # @parameters α_R0

    # # state variables
    # @variables B(t) [unit = uQ] L(t) [unit = uQ]  R(t) [unit = uQ] α_R(t) cumresp(t) [unit = uQ]

    # @variables α_L(t) tmpNoUnit(t)
    # # fluxes
    # @variables (begin
    #     syn_Enz(t),[unit = uQ/uT], r_M(t),[unit = uQ/uT], tvr_B(t),[unit = uQ/uT],
    #     dec_LPot(t),[unit = uQ/uT], dec_L(t),[unit = uQ/uT], dec_RPot(t),[unit = uQ/uT],
    #     dec_R(t),[unit = uQ/uT], u_C(t),[unit = uQ/uT],
    #     C_synBCt(t),[unit = uQ/uT], C_synBC(t),[unit = uQ/uT], r_G(t),[unit = uQ/uT],
    #     r_tvr(t),[unit = uQ/uT], syn_B(t),[unit = uQ/uT],
    #     r_B(t),[unit = uQ/uT], r_GEnz(t),[unit = uQ/uT], r_O(t),[unit = uQ/uT]
    # end)

    # #------------- N ---------------------
    # @variables L_N(t) [unit = uQ] R_N(t) [unit = uQ] I_N(t) [unit = uQ]
    # @variables (begin
    #     Φ_N(t),[unit = uQ/uT], Φ_Nu(t),[unit = uQ/uT], Φ_NB(t),[unit = uQ/uT],
    #     u_PlantN(t),[unit = uQ/uT], u_NOM(t),[unit = uQ/uT], u_immNPot(t),[unit = uQ/uT],
    #     u_NPot(t),[unit = uQ/uT], N_synBN(t),[unit = uQ/uT], M_ImbN(t),[unit = uQ/uT],
    #     β_NL(t), β_NR(t)
    # end)
    # @parameters β_NEnz β_NB l_N [unit = uT^-1] ν_N i_BN [unit = uT^-1]

    # # need to be specified by coupled system:
    # @variables (begin
    #     i_L(t),[unit = uQ/uT], β_Ni(t), i_IN(t),[unit = uQ/uT],
    #     u_PlantNmax(t),[unit = uQ/uT], k_PlantN(t),[unit = uT^-1]
    # end)

    # replace patttern to remove units: ,?\[unit[^\]]+\]
    @parameters t
    D = Differential(t)

    @parameters ϵ ϵ_tvr κ_E a_E  m  τ
    @parameters k_L  k_R  k_mN_L k_mN_R k_N=k_N
    @parameters α_R0
    @parameters ρ_CBtvr = 0.0 # proportion of carbon resorption on turnover

    @variables (begin
        B(t),  L(t),   R(t),  cumresp(t),
        dB(t), dL(t), dR(t), r_tot(t),
        syn_Enz(t), tvr_Enz(t), r_M(t), tvr_B(t),
        E_L(t), E_R(t),
        dec_LPot(t), dec_L(t), dec_RPot(t),
        dec_R(t), u_C(t),
        C_synBCt(t), C_synBC(t), r_G(t),
        r_tvr(t),
        r_B(t), r_GEnz(t), r_O(t),
        tvr_B0(t), resorp_C(t),
        # need to be defined by component across all elements
        α_L(t), α_R(t),
        # need to be specified by coupled system:
        i_L(t), syn_B(t)
    end)

    eqs = [
        D(B) ~ dB, dB ~ syn_B - tvr_B,
        D(L) ~ dL, dL ~ -dec_L + i_L,
        D(R) ~ dR, dR ~ -dec_R + ϵ_tvr*tvr_B + (1-κ_E)*syn_Enz,
        syn_Enz ~ a_E*B, tvr_Enz ~ syn_Enz,
        r_M ~ m*B,
        # element resorption changes stoichiometry: β_EB -> β_EBtvr
        tvr_B ~ tvr_B0 - resorp_C,
        tvr_B0 ~ τ*B, # turnover before resorption
        resorp_C ~ ρ_CBtvr * tvr_B0,
        dec_LPot ~ k_L * L,
        dec_L ~ dec_LPot*(α_L * syn_Enz)/(k_mN_L + α_L*syn_Enz),
        dec_RPot ~ k_R * R,
        dec_R ~ dec_RPot*(α_R * syn_Enz)/(k_mN_R + α_R*syn_Enz),
        u_C ~ dec_L + dec_R + κ_E*syn_Enz,
        C_synBCt ~ u_C - syn_Enz/ϵ - r_M,
        #C_synBC ~ IfElse.ifelse(C_synBCt > 0.0, ϵ*C_synBCt, C_synBCt),
        C_synBC ~ C_synBCt - ((C_synBCt > 0.0) * (1-ϵ)*C_synBCt),
        #syn_B ~ C_synBC # define in all-elements model
        r_O ~ u_C - (syn_Enz/ϵ +  syn_B + r_G + r_M),
        #D(cumresp) ~ r_B + r_tvr,
        r_tvr ~ (1-ϵ_tvr)*tvr_B,
        r_B ~ r_GEnz + r_G + r_M + r_O,
        r_GEnz ~ (1-ϵ)/ϵ * syn_Enz,
        #r_G ~ IfElse.ifelse(syn_B.val > 0.0, (1-ϵ)/ϵ * syn_B, 0.0),
        r_G ~ (syn_B.val > 0.0) * (1-ϵ)/ϵ * syn_B,
        r_tot ~ r_B + r_tvr,
        E_L ~ (α_L * syn_Enz)/k_N,
        E_R ~ (α_R * syn_Enz)/k_N,
        ]
    ODESystem(eqs; name)
end

function sesam3N(;name, sC = sesam3C(name=:sC))
    @parameters t
    D = Differential(t)

    sts = @variables (begin
        L_N(t),  R_N(t),  I_N(t),
        dL_N(t),  dR_N(t),  dI_N(t),
        Φ_N(t), Φ_Nu(t), Φ_NB(t), Φ_Ntvr(t),
        u_PlantN(t), u_NOM(t), u_immNPot(t),
        u_NPot(t), N_synBN(t), M_ImbN(t),
        β_NL(t), β_NR(t),
        leach_N(t),
        resorp_N(t), β_NBtvr(t),
        # need to be specified by coupled system:
        β_Ni(t), i_IN(t),
        u_PlantNmax(t), k_PlantN(t),
        SOC(t), SON(t), β_NSOM(t), BtoSOC(t)
    end)
    ps = @parameters β_NEnz β_NB l_N  ν_N i_BN ρ_NBtvr=0.0

    @unpack L, R, B, dec_L, dec_R, i_L, ϵ_tvr, tvr_B, tvr_B0, syn_B, syn_Enz, tvr_Enz, r_tvr, κ_E = sC

    eqs = [
        β_NL ~ L/L_N, β_NR ~ R/R_N,
        D(L_N) ~ dL_N, dL_N ~ -dec_L/β_NL + i_L/β_Ni,
        D(R_N) ~ dR_N, dR_N ~ -dec_R/β_NR + ϵ_tvr*tvr_B/β_NBtvr + (1-κ_E)*tvr_Enz/β_NEnz,
        D(I_N) ~ dI_N,
        u_PlantN ~ min(u_PlantNmax, k_PlantN*I_N),
        dI_N ~ i_IN - u_PlantN - leach_N + Φ_N,
        leach_N ~ l_N*I_N,
        Φ_N ~ Φ_Nu + Φ_NB + Φ_Ntvr,
        Φ_Ntvr ~ r_tvr/β_NBtvr,
        Φ_Nu ~ (1-ν_N) * u_NOM,
        u_NPot ~ ν_N * u_NOM + u_immNPot + resorp_N,
        u_NOM ~ dec_L/β_NL + dec_R/β_NR + κ_E*tvr_Enz/β_NEnz,
        u_immNPot ~ i_BN * I_N,
        N_synBN ~ u_NPot - syn_Enz/β_NEnz,
        M_ImbN ~ u_NPot - (syn_B/β_NB + syn_Enz/β_NEnz),
        Φ_NB ~ M_ImbN - u_immNPot,
        # make sure ρ_NBtvr >= 0
        resorp_N ~ ρ_NBtvr * tvr_B0/β_NB,
        β_NBtvr ~ tvr_B / (tvr_B0/β_NB - resorp_N),
        # observables for diagnostic output
        SOC ~ L + R + B,
        SON ~ L_N + R_N + B/β_NB,
        β_NSOM ~ SOC/SON,
        BtoSOC ~ B/SOC,
    ]
    extend(ODESystem(eqs, t, sts, ps; name), sC)
end

function get_revenue_eq_sesam3CN(sN)
    @parameters t
    @unpack α_L, α_R, dec_L, dec_R, β_NL, β_NR, β_NEnz, syn_Enz = sN
    @unpack u_immNPot, u_PlantN, ν_N = sN
    sts = @variables (begin
        α_LT(t), α_RT(t),
        p_uNmic(t), p_oNmic(t),
        invest_L(t), invest_R(t), return_L(t), return_R(t), revenue_L(t), revenue_R(t),
        #invest_Ln(t), invest_Rn(t), return_Ln(t), return_Rn(t),
        revenue_sum(t)
    end)
    # need to be defined in coupled component:
    @variables w_C(t), w_N(t)
    eqs = [
        p_uNmic ~ u_immNPot/(u_immNPot + u_PlantN),
        p_oNmic ~ ν_N+(1-ν_N)*p_uNmic,
        invest_L ~ α_L*syn_Enz*(w_C + w_N/β_NEnz),
        invest_R ~ α_R*syn_Enz*(w_C + w_N/β_NEnz),
        # return_L ~ dec_L * (w_C + w_N/β_NL),
        # return_R ~ dec_R * (w_C + w_N/β_NR),
        return_L ~ dec_L * (w_C + w_N/β_NL*p_oNmic),
        return_R ~ dec_R * (w_C + w_N/β_NR*p_oNmic),
        revenue_L ~ return_L / invest_L,
        revenue_R ~ return_R / invest_R,
        revenue_sum ~ revenue_L + revenue_R,
        α_LT ~ revenue_L/revenue_sum,
        α_RT ~ revenue_R/revenue_sum,
        # auxilary for plotting
        # invest_Ln ~ invest_L/(invest_L + invest_R),
        # invest_Rn ~ invest_R/(invest_L + invest_R),
        # return_Ln ~ return_L/(return_L + return_R),
        # return_Rn ~ return_R/(return_L + return_R),
        ]
    (;eqs, sts)
end

function sesam3CN(;name, δ=40.0, max_w=12, use_seam_revenue=false, sN=sesam3N(name=:sN))
    @parameters t
    D = Differential(t)
    @unpack α_L, α_R, syn_B, B, C_synBC, β_NB, N_synBN, tvr_B, τ = sN
    @unpack u_immNPot, u_PlantN, ν_N = sN
    sts = @variables (begin
        C_synBN(t),
        #dα_L(t),
        dα_R(t),
        w_C(t), w_N(t), lim_C(t), lim_N(t)
    end)
    ps = @parameters δ=δ
    eqs_rev, sts_rev = use_seam_revenue ?
        get_revenue_eq_seam(sN) : get_revenue_eq_sesam3CN(sN)
    @variables α_LT(t) α_RT(t)
    eqs = [
        C_synBN ~ β_NB*N_synBN,
        syn_B ~ min(C_synBC, C_synBN),
        # need minimum, otherwise danger of Inf and nan -> instability
        w_C ~ exp(min(max_w, -δ/tvr_B*(C_synBC - syn_B))),
        w_N ~ exp(min(max_w, -δ/tvr_B*(C_synBN - syn_B))),
        lim_C ~ w_C/(w_C + w_N), lim_N ~ w_N/(w_C + w_N), # normalized for plot
        # α_LT, α_RT by get_revenue_eq_X
        #D(α_L) ~ dα_L, dα_L ~ (α_LT - α_L)*(τ + abs(syn_B)/B),
        α_L ~ 1 - α_R,
        D(α_R) ~ dα_R, dα_R ~ (α_RT - α_R)*(τ + abs(syn_B)/B),
        ]
    extend(ODESystem(vcat(eqs,eqs_rev), t, vcat(sts, sts_rev), ps; name), sN)
end

# sesam3 now refers to sesam3CNP
#sesam3(args...;kwargs...) = sesam3CN(args...;kwargs...)

"""
R pool is a mixture of microbial turnover and enzymes.
Here the stable C/N ratio is computed based on given parameterization.
"""
function calculate_β_NR_sesam3(p,s)
    w_REnz = p[s.a_E]*(1-p[s.κ_E]) # a_E (1-κ_E) B
    w_RBtvr = p[s.ϵ_tvr]*p[s.τ]       # τ ϵ_tvr B
    # β_NBtvr may differ from β_NB because of elemental resorption
    # if ρ is not in parameter dictionary, its default is zero
    ρ_CBtvr = (hasproperty(s, :ρ_CBtvr) && haskey(p, s.ρ_CBtvr)) ? p[s.ρ_CBtvr] : 0.0
    ρ_NBtvr = (hasproperty(s, :ρ_NBtvr) && haskey(p, s.ρ_NBtvr)) ? p[s.ρ_NBtvr] : 0.0
    β_NBtvr = p[s.β_NB] * (1-ρ_CBtvr)/(1-ρ_NBtvr)
    β_NR = 1/((w_REnz/p[s.β_NEnz] + w_RBtvr/β_NBtvr)/(w_REnz + w_RBtvr))
end


text:
"Fresnel reflection, needs pol_type soon too!"
function Fᵣ(n::FT,θᵢ::FT) where FT
    #@show rad2deg(θᵢ)
    # Refractive index of air
    nᵢ = FT(1)
    # Angle of transmitted light in medium:
    θₜ = asin(nᵢ * sin(θᵢ)/n)
    #(sin(α - θₛ)^2 / sin(α + θₛ)^2 +  tan(α - θₛ)^2 / tan(α + θₛ)^2)/2
    r_perpendicular = (nᵢ * cos(θᵢ) - n * cos(θₜ)) / (nᵢ * cos(θᵢ) + n * cos(θₜ))
    r_parallel      = (nᵢ * cos(θₜ) - n * cos(θᵢ)) / (nᵢ * cos(θₜ) + n * cos(θᵢ))
    return (r_perpendicular^2 +  r_parallel^2)/2
end


text:
abstract type Renderable end

"""
Return true if an object or type is *directly renderable*, false otherwise.

New types should implement the `isrenderable(t::Type{NewType})` method.
"""
function isrenderable end

isrenderable(object) = isrenderable(typeof(object))
isrenderable(::Type{R}) where R <: Renderable = true
isrenderable(t::Type) = hasmethod(add_renderable!, Tuple{RenderModel, t})

isrenderable(t::Type{Roadway}) = true

"""
A basic drawable rectangle representing a car.
An arrow indicates the heading direction of the car.

    ArrowCar{A<:AbstractArray{Float64}, C<:Colorant} <: Renderable
    ArrowCar(pos::AbstractArray, angle::Float64=0.0; length = 4.8, width = 1.8,  color=colortheme["COLOR_CAR_OTHER"], text="", id=0)
    ArrowCar(x::Real, y::Real, angle::Float64=0.0; length = 4.8, width = 1.8,  color=colortheme["COLOR_CAR_OTHER"], text="", id=0)
"""
@with_kw struct ArrowCar{A<:AbstractArray{Float64}, C<:Colorant} <: Renderable
    pos::A         = SVector(0.0, 0.0)
    angle::Float64 = 0.0
    length::Float64 = 4.8
    width::Float64 = 1.8
    color::C       = colortheme["COLOR_CAR_OTHER"]
    text::String   = "" # some debugging text to print by the car
    id::Int        = 0
end
ArrowCar(pos::AbstractArray, angle::Float64=0.0; length = 4.8, width = 1.8,  color=colortheme["COLOR_CAR_OTHER"], text="", id=0) = ArrowCar(pos, angle, length, width, color, text, id)
ArrowCar(x::Real, y::Real, angle::Float64=0.0; length = 4.8, width = 1.8,  color=colortheme["COLOR_CAR_OTHER"], text="", id=0) = ArrowCar(SVector(x, y), angle, length, width, color, text, id)
function ArrowCar(car::E; color=id_to_color(car.id), kwargs...) where {E<:Entity}
    ArrowCar(
        posg(car)...;
        length=length(car.def), width=AutomotiveDrivingModels.width(car.def),
        color=color, kwargs...
    )
end
ArrowCar(entity::Entity, color::Color) = ArrowCar(entity; color=color)

function add_renderable!(rm::RenderModel, c::ArrowCar)
    x = c.pos[1]
    y = c.pos[2]
    add_instruction!(rm, render_vehicle, (x, y, c.angle, c.length, c.width, c.color))
    add_instruction!(rm, render_text, (c.text, x, y-c.width/2 - 2.0, 10, colorant"white"))
    return rm
end

"""
A drawable rectangle with rounded corners representing an `entity`.
"""
@with_kw struct EntityRectangle{S,D,I, C<:Colorant} <: Renderable
    entity::Entity{S,D,I}
    color::C = id_to_color(entity.id)
end

function render_entity_rectangle(ctx::CairoContext, er::EntityRectangle)
    x, y, yaw = posg(er.entity.state)
    w, h = length(er.entity.def), AutomotiveDrivingModels.width(er.entity.def)
    cr = 0.5 # [m]
    save(ctx); translate(ctx, x, y); rotate(ctx, yaw);
    color_fill = er.color
    color_line = weighted_color_mean(.4, colorant"black", color_fill)
    render_round_rect(ctx, 0, 0, w, h, 1., cr, color_fill, true, true, color_line, 1.)
    restore(ctx)
end
add_renderable!(rm::RenderModel, er::EntityRectangle) = add_instruction!(rm, render_entity_rectangle, (er,))

"""
A drawable arrow representing the current velocity vector of an `entity`.
The arrow points to the location where the vehicle will be one second in the future (assuming linear motion).
"""
@with_kw struct VelocityArrow{S,D,I, C<:Colorant} <: Renderable
    entity::Entity{S,D,I}
    color::C = colorant"white"
end

function render_velocity_arrow(ctx::CairoContext, va::VelocityArrow)
    x, y, yaw = posg(va.entity.state)
    vx, vy = velg(va.entity.state)
    save(ctx); translate(ctx, x, y);
    render_arrow(ctx, [[0.  vx];[0. vy]], va.color, .3, .8, ARROW_WIDTH_RATIO=1., ARROW_ALPHA=.12pi, ARROW_BETA=.6pi)
    restore(ctx)
end
add_renderable!(rm::RenderModel, va::VelocityArrow) = add_instruction!(rm, render_velocity_arrow, (va,))

"""
A drawable 'fancy' svg image of a race car.
The car is placed at the position of `entity` and the width and length are scaled accordingly.
The color of the car can be specified using the `color` keyword.
"""
@with_kw struct FancyCar{C<:Colorant, S, D, I} <: Renderable
    car::Entity{S, D, I}
    color::C = id_to_color(car.id)
end

function add_renderable!(rm::RenderModel, fc::FancyCar)
    x, y, yaw = posg(fc.car.state)
    l, w = length(fc.car.def), AutomotiveDrivingModels.width(fc.car.def)
    add_instruction!(rm, render_fancy_car, (x, y, yaw, l, w, fc.color))
    return rm
end

"""
A drawable 'fancy' svg image of a pedestrian.
The pedestrian is placed at the position of `entity` and the width and length of the original image are scaled accordingly.
The color of the pedestrian can be specified using the `color` keyword.
"""
@with_kw struct FancyPedestrian{C<:Colorant, S, D, I} <: Renderable
    ped::Entity{S, D, I}
    color::C = id_to_color(ped.id)
end

function add_renderable!(rm::RenderModel, fp::FancyPedestrian)
    x, y, yaw = posg(fp.ped.state)
    l, w = length(fp.ped.def), AutomotiveDrivingModels.width(fp.ped.def)
    add_instruction!(rm, render_fancy_pedestrian, (x, y, yaw, l, w, fp.color))
    return rm
end

"""
Helper function for directly rendering entities, takes care of wrapping them in renderable objects
"""
function add_renderable!(
    rendermodel::RenderModel,
    entity::E,
    color::Union{Nothing, Colorant}=nothing
) where {E<:Entity}
    if color === nothing
        color = id_to_color(entity.id)
    end
    if rendermode == :fancy
        fe = (class(entity.def) == AgentClass.PEDESTRIAN ? FancyPedestrian(ped=entity, color=color) : FancyCar(car=entity, color=color))
        add_renderable!(rendermodel, fe)
    else
        er = EntityRectangle(entity=entity, color=color)
        add_renderable!(rendermodel, er)
        va = VelocityArrow(entity=entity, color=colorant"white")
        add_renderable!(rendermodel, va)
    end
    return rendermodel
end

"""
Helper function for directly rendering frames of renderable entities
"""
function add_renderable!(rendermodel::RenderModel, scene::Scene{E}) where {E<:Entity}
    for entity in scene
        add_renderable!(rendermodel, entity)
    end
    return rendermodel
end

# render nothing by not doing anything
add_renderable!(rm::RenderModel, ::Nothing) = rm

"""
Render function for text
"""
function add_renderable!(rm::RenderModel, t::String)
    add_renderable!(rm, TextOverlay(text=split(t, '\n')))
    return rm
end

max_stars_repo_name: tawheeler/AutoViz

max_issues_repo_name: tawheeler/AutoViz

max_forks_repo_name: tawheeler/AutoViz


text:
using PyCall
@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.lines as lines
@pyimport matplotlib.animation as anim

export sptPlot,
       SeisPlot,
       waveAnim

include("sptPlot.jl")
include("SeisPlot.jl")
include("waveAnim.jl")

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
# Saturated integer arithmetic in [1..typemax]
⊕(a::T,b::T) where {T<:Unsigned}= a+b>=a ? a+b : typemax(T)
⊖(a::T,b::T) where {T<:Unsigned}= a-b<=a ? a-b : one(T)
⊕(a::T,b::T) where {T<:Signed}= a+b>=a ? a+b : typemax(T)
⊖(a::T,b::T) where {T<:Signed}= a-b>0 ? a-b : one(T)

# Side of hypercube with range φ
α(φ)= 2*round(Int,φ)+1
neighborhood(φ,Ndim)= ntuple(i->α(φ),Ndim)
mean_kernel(φ,Ndim)= ones(neighborhood(φ,Ndim)) ./ α(φ).^Ndim

# BitVector x SparseMatrixCSC{Bool}, used in overlap
import Base: *
*(z::BitVector,W::SparseMatrixCSC)= Vector(z)*W
*(z::Adjoint{Bool,BitVector},W::SparseMatrixCSC)= Vector(z.parent)'W
*(W::Adjoint{<:Any,<:SparseMatrixCSC},z::BitVector)= W*Vector(z)
*(W::SparseMatrixCSC,z::BitVector)= W*Vector(z)

# BitVector x Matrix{Bool}, used in overlap
*(z::BitVector,W::Matrix)= Vector(z)*W
*(z::Adjoint{Bool,BitVector},W::Matrix)= Vector(z.parent)'W
*(W::Adjoint{<:Any,<:Matrix},z::BitVector)= W*Vector(z)
*(W::Matrix,z::BitVector)= W*Vector(z)

#  Work around the narrow result types for SparseMatrixCSC{Bool} x Bool

# FIX SparseMatrixCSC{Bool} x Vector{Bool}
# TODO still needed?
matmul_elt_op(a,b)= a*b+a*b
*(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA<:Bool,S,Tx} =
    (T = Base.promote_op(matmul_elt_op, TA,Tx); SparseArrays.mul!(similar(x, T, A.m), A, x, one(T), zero(T)))
*(adjA::Adjoint{<:Any,<:SparseMatrixCSC{TA,S}}, x::StridedVector{Tx}) where {TA<:Bool,S,Tx} =
    (A = adjA.parent; T = Base.promote_op(matmul_elt_op, TA,Tx); SparseArrays.mul!(similar(x, T, A.n), adjoint(A), x, one(T), zero(T)))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Perl"
version = v"5.30.3"

# Collection of sources required to build perl
# with a few extra modules for polymake
sources = [
    ArchiveSource("https://www.cpan.org/src/5.0/perl-$version.tar.gz", "32e04c8bb7b1aecb2742a7f7ac0eabac100f38247352a73ad7fa104e39e7406f"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/I/IS/ISHIGAKI/JSON-4.02.tar.gz", "444a88755a89ffa2a5424ab4ed1d11dca61808ebef57e81243424619a9e8627c"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/J/JO/JOSEPHW/XML-Writer-0.625.tar.gz", "e080522c6ce050397af482665f3965a93c5d16f5e81d93f6e2fe98084ed15fbe"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/J/JS/JSTOWE/TermReadKey-2.38.tar.gz", "5a645878dc570ac33661581fbb090ff24ebce17d43ea53fd22e105a856a47290"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/H/HA/HAYASHI/Term-ReadLine-Gnu-1.36.tar.gz", "9a08f7a4013c9b865541c10dbba1210779eb9128b961250b746d26702bab6925"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-1.02.tar.gz", "4506c387043aa6a77b455f00f57409f3720aa7e553495ab2535263b4ed1ea12a"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/P/PE/PERIGRIN/XML-NamespaceSupport-1.12.tar.gz", "47e995859f8dd0413aa3f22d350c4a62da652e854267aa0586ae544ae2bae5ef"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-Base-1.09.tar.gz", "66cb355ba4ef47c10ca738bd35999723644386ac853abbeb5132841f5e8a2ad0"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/M/MA/MANWAR/SVG-2.84.tar.gz", "ec3d6ddde7a46fa507eaa616b94d217296fdc0d8fbf88741367a9821206f28af"),
    DirectorySource("./bundled")
]

# Bash recipe for building
script = raw"""
perldir=`ls -1d perl-*`
cd $WORKSPACE/srcdir/
for dir in *;
do
   [[ "$dir" == "perl-"* ]] && continue;
   [[ "$dir" == "patches" ]] && continue;
   # build extra perl modules in-tree
   # the names of the extra modules also need to appear in the
   # config.sh for all cross-compilation architectures
   sed -i '1s/^/$ENV{PERL_CORE}=0;/' $dir/Makefile.PL
   mv $dir $perldir/cpan/${dir%-*};
done

cd $perldir/

# allow combining relocation with shared library
# add patch to find binary location from shared library
atomic_patch -p1 ../patches/allow-relocate.patch

# replace some library checks that wont work in the cross-compile environment
# with the required values
atomic_patch -p1 ../patches/cross-nolibchecks.patch

if [[ $target != x86_64-linux* ]] && [[ $target != i686-linux* ]]; then
   # cross build with supplied config.sh
   # build native miniperl
   src=`pwd`
   mkdir host
   pushd host
   ../Configure -des -Dusedevel -Duserelocatableinc -Dmksymlinks -Dosname=linux -Dcc=$CC_FOR_BUILD -Dld=$LD_FOR_BUILD -Dar=$AR_FOR_BUILD -Dnm=$NM_FOR_BUILD -Dlibs=-lm
   make -j${nproc} miniperl
   make -j${nproc} generate_uudmap
   cp -p miniperl $prefix/bin/miniperl-for-build
   popd

   # copy and use prepared configure information
   cp ../patches/config-$target.sh config.sh
   ./Configure -K -S
else
   # native

   # config overrides
   if [[ $target = *-gnu ]]; then
      # disable xlocale.h usage (which was removed in recent glibc)
      cp ../patches/config.arch.gnu config.arch
   fi

   ./Configure -des -Dcc="$CC" -Dprefix=$prefix -Duserelocatableinc -Dprocselfexe -Duseshrplib -Dsysroot=/opt/$target/$target/sys-root -Dccflags="-I${prefix}/include" -Dldflags="-L${libdir} -Wl,-rpath,${libdir}" -Dlddlflags="-shared -L${libdir} -Wl,-rpath,${libdir}"
fi

make -j${nproc} depend
make -j${nproc}

make install

# put a libperl directly in lib
cd $libdir
ln -s perl5/*/*/CORE/libperl.${dlext} libperl.${dlext}

# resolve case-ambiguity:
cd $libdir/perl5/5.*.*
mv Pod/* pod
rmdir Pod

# remove sysroot and target flags from stored compiler flags:
sed -i -e "s#--sysroot[ =]\S\+##g" \
       -e "s#-target[ =]\S\+##g" \
       ${prefix}/*/perl5/*/*/Config_heavy.pl

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "macos")
    Platform("x86_64", "linux"; libc="glibc")
    Platform("i686", "linux"; libc="glibc")
    Platform("x86_64", "linux"; libc="musl")
    Platform("i686", "linux"; libc="musl")
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("perl", :perl)
    LibraryProduct("libperl", :libperl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Readline_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
import Base: isidentifier, is_id_start_char, is_id_char

const RESERVED_WORDS = Set(["begin", "while", "if", "for", "try",
    "return", "break", "continue", "function", "macro", "quote", "let",
    "local", "global", "const", "abstract", "typealias", "type", "bitstype",
    "immutable", "ccall", "do", "module", "baremodule", "using", "import",
    "export", "importall", "end", "else", "elseif", "catch", "finally"])
if VERSION >= v"0.4.0-dev+757"
    push!(RESERVED_WORDS, "stagedfunction")
end

function identifier(s::AbstractString)
    s = normalize_string(s)
    if !isidentifier(s)
        s = makeidentifier(s)
    end
    symbol(in(s, RESERVED_WORDS) ? "_"*s : s)
end

function makeidentifier(s::AbstractString)
    i = start(s)
    done(s, i) && return "x"

    res = IOBuffer(sizeof(s) + 1)

    (c, i) = next(s, i)
    under = if is_id_start_char(c)
        write(res, c)
        c == '_'
    elseif is_id_char(c)
        write(res, 'x', c)
        false
    else
        write(res, '_')
        true
    end

    while !done(s, i)
        (c, i) = next(s, i)
        if c != '_' && is_id_char(c)
            write(res, c)
            under = false
        elseif !under
            write(res, '_')
            under = true
        end
    end

    return takebuf_string(res)
end

function make_unique(names::Vector{Symbol})
    seen = Set{Symbol}()
    names = copy(names)
    dups = Int[]
    for i in 1:length(names)
        name = names[i]
        in(name, seen) ? push!(dups, i) : push!(seen, name)
    end
    for i in dups
        nm = names[i]
        k = 1
        while true
            newnm = symbol("$(nm)_$k")
            if !in(newnm, seen)
                names[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end

    return names
end

#' @description
#'
#' Generate standardized names for columns of a DataFrame. The
#' first name will be :x1, the second :x2, etc.
#'
#' @field n::Integer The number of names to generate.
#'
#' @returns names::Vector{Symbol} A vector of standardized column names.
#'
#' @examples
#'
#' DataFrames.gennames(10)
function gennames(n::Integer)
    res = Array(Symbol, n)
    for i in 1:n
        res[i] = symbol(@sprintf "x%d" i)
    end
    return res
end

#' @description
#'
#' Count the number of missing values in an Array.
#'
#' NOTE: This function always returns 0.
#'
#' @field a::Array The Array whose missing values are to be counted.
#'
#' @returns count::Int The number of missing values in `a`.
#'
#' @examples
#'
#' DataFrames.countna([1, 2, 3])
countna(a::Array) = 0

#' @description
#'
#' Count the number of missing values in a DataArray.
#'
#' @field da::DataArray The DataArray whose missing values are to be counted.
#'
#' @returns count::Int The number of missing values in `a`.
#'
#' @examples
#'
#' DataFrames.countna(@data([1, 2, 3]))
countna(da::DataArray) = sum(da.na)

#' @description
#'
#' Count the number of missing values in a PooledDataArray.
#'
#' @field pda::PooledDataArray The PooledDataArray whose missing values
#'        are to be counted.
#'
#' @returns count::Int The number of missing values in `a`.
#'
#' @examples
#'
#' DataFrames.countna(@pdata([1, 2, 3]))
function countna(da::PooledDataArray)
    res = 0
    for i in 1:length(da)
        res += da.refs[i] == 0
    end
    return res
end

function _setdiff{T}(a::AbstractVector{T}, b::AbstractVector{T})
    diff = T[]
    for val in a
        if !(val in b)
            push!(diff, val)
        end
    end
    diff
end
# because unions and parametric types don't compose, yet
function _setdiff{T}(a::AbstractVector{T}, b::T)
    diff = T[]
    for val in a
        if !(val in b)
            push!(diff, val)
        end
    end
    diff
end

function _uniqueofsorted(x::Vector)
    idx = fill(true, length(x))
    lastx = x[1]
    for i = 2:length(x)
        if lastx == x[i]
            idx[i] = false
        else
            lastx = x[i]
        end
    end
    x[idx]
end

# Gets the name of a function. Used in groupedataframe/grouping.jl
function _fnames{T<:Function}(fs::Vector{T})
    λcounter = 0
    names = map(fs) do f
        name = string(f)
        if name == "(anonymous function)" # Anonymous functions with Julia < 0.5
            λcounter += 1
            name = "λ$(λcounter)"
        end
        name
    end
    names
end


text:
doc"""
f(state) = exp{E / kT} is known as Boltzmann distribution.
The ratio of f(state1) and f(state2) is known as Boltzmann factor.
"""
function boltzmann_factor(g::Function, temperature::Float64)
    if temperature == 0.0
        return 0.0
    end
    ΔE = g()
    if ΔE < 0.0
        return 1.0
    end
    return exp(-ΔE / temperature)
end

function acceptance(sym::Type{Val{:energy}}, f::Function, states::StatePair, temperature::Float64)
    g = () -> f(states.proposal) - f(states.current)
    return boltzmann_factor(g, temperature)
end

function acceptance(sym::Type{Val{:loss}}, f::Function, states::StatePair, temperature::Float64)
    g = () -> f(states.proposal, states.current)
    return boltzmann_factor(g, temperature)
end

acceptance(sym::Symbol, f::Function, states::StatePair, temperature::Float64) = acceptance(Val{sym}, f, states, temperature)

type SAConfig
    initiate::Function
    func::Function
    func_type::Symbol
    next_state::Function
    acceptance::Function
    # markov_step=10000::Int64
    # T_0 = 10000
end

doc"""
simulated annealing algorithm
"""
function SA_algorithm()
    x = MH_init(initiate)
    cool_sch = CoolingScheduler()
    stable_counter = StabilityCounter()
    while !cool_sch.isfrozen() && !stable_counter.ishalt()
        MH_body!(x, next_state, acceptance, markov_step)
        cool_sch.cooling()
        stable_counter.update()
    end
    return x
end


text:
# The offline policy in "ilqg_pcontrol_policy_20sec_2.jld2" was computed
# using the iterative LQG algorithm in belief space presented in
# "Motion Planning under Uncertainty using Iterative Local Optimization in
# Belief Space" by Jur van den Berg, Sachin Patil, and Ron Alterovitz.
# The implementation of the algorithm that derived the policy was written in
# Julia v0.6 and is not currently included in this codebase, although we plan
# to migrate it to Julia v1.3 as well.

using SACBP

function ilqgControlUpdate(ilqgPolicy::Dict{String,Any},b::BelMvNormal{<:Real},UArray::Vector{<:MControl2D},
                           u_param_min::Vector{<:Real},u_param_max::Vector{<:Real},dto::Float64,dtexec::Float64)
    uVec,bVec,lVec,LVec = ilqgPolicy["uVec_new"],ilqgPolicy["bVec_new"],ilqgPolicy["lVec"],ilqgPolicy["LVec"];
    index = Int64(round(b.t/dto,digits=5))+1;
    if (bVec[index].t != b.t)
        error("Time Parameter Inconsistent.")
    end
    u_val = LVec[index]*(b.params - bVec[index].params) + (lVec[index] + vec(uVec[index]));
    for ii = 1:length(u_val)
        if u_val[ii] < u_param_min[ii]
            u_val[ii] = u_param_min[ii];
        elseif u_val[ii] > u_param_max[ii]
            u_val[ii] = u_param_max[ii];
        end
    end
    UArray_new = copy(UArray)
    for ii = 1:length(UArray_new)
        if UArray_new[ii].t < round(b.t + dtexec,digits=5)
            UArray_new[ii] = MControl2D(UArray_new[ii].t,u_val);
        end
    end
    return UArray_new
end


text:
# script to test new elimination tree from graph algorithm

workspace()
include("../../../src/Graph.jl")
include("../../../src/Tree.jl")

using GraphModule, TreeModule, Base.Test

A = [1 1 0 0 0; 1 1 1 0 0; 0 1 1 1 1; 0 0 1 1 1; 0 0 1 1 1]

g = Graph(sparse(A))

rng = Random.MersenneTwister(123554);
# Define number of test matrices
nn = 1000

function parentsFromOldAlgo(g)
  tic()
  t = TreeModule.createTreeFromGraph(g);
  tOld = toq()

  parOld = zeros(1:length(g.adjacencyList))
  for iii=1:length(t.nodes)
    n=t.nodes[iii]
    parOld[iii] = n.parent
  end
  return parOld, tOld
end

tNew3Arr = zeros(nn)
tNew4Arr = zeros(nn)
tOldArr = zeros(nn)

@testset "Elimination Tree Derivation from random Matrices" begin

    for iii=1:nn

        # take random dimension
        dim = rand(rng,50:200)
        density = rand(rng,0.1:0.1:0.6)
        # create random sparse matrix
        A = sprand(rng,dim,dim,density)
        A = A+A'

        # create graph from A and make Graph chordal
        g = Graph(A)
        # create Tree from Graph
        tic()
        par3 = TreeModule.etree3(g)
        tNew3 = toq()

        tic()
        par4 = TreeModule.etree4(g)
        tNew4 = toq()

        parOld,tOld = parentsFromOldAlgo(g)

        tNew3Arr[iii] = tNew3
        tNew4Arr[iii] = tNew4
        tOldArr[iii] = tOld

        @test par4 == parOld
        @test par3 == parOld
        println("$(iii) / $(nn)")
    end
end


text:
include("printing.jl")
include("base.jl")
include("circular_array_buffer.jl")
include("device.jl")
include("sum_tree.jl")


text:
using Falcons
using Documenter

DocMeta.setdocmeta!(Falcons, :DocTestSetup, :(using Falcons); recursive=true)

makedocs(;
    modules=[Falcons],
    authors="Yusuke Takase",
    repo="https://github.com/yusuke-takase/Falcons.jl/blob/{commit}{path}#{line}",
    sitename="Falcons.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://yusuke-takase.github.io/Falcons.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Scanning" => "Scanning.md",
        "Mapmake" => "Mapmake.md",
    ],
)

deploydocs(;
    repo="github.com/yusuke-takase/Falcons.jl",
)


text:
G = 6.673*10^-11
# Insira os valores da massa e do raio do planeta o qual iremos calcular a velocidade de escape
M = 3.2*10^20
R = 5*10^5
v = sqrt((2*G*M)/R)

max_stars_repo_path: Velocidade de escape.jl

max_stars_repo_name: PedroMoraesFranco/fisica-computacional

max_issues_repo_path: Velocidade de escape.jl

max_issues_repo_name: PedroMoraesFranco/fisica-computacional

max_forks_repo_path: Velocidade de escape.jl

max_forks_repo_name: PedroMoraesFranco/fisica-computacional