text:
import Base:+,*,-,\,^,sparse,full,copy

using Base.LinAlg.BLAS:gemm!,gemm,axpy!

"""
Represents a matrix of the form A + BDBᵀ.
"""
type SymWoodbury{T,AType, BType, DType} <: AbstractMatrix{T}
  A::AType;
  B::BType;
  D::DType;
end

"""
    SymWoodbury(A, B, D)

Represents a matrix of the form A + BDBᵀ.
"""
function SymWoodbury{T}(A, B::AbstractMatrix{T}, D)
    n = size(A, 1)
    k = size(B, 2)
    if size(A, 2) != n || size(B, 1) != n || size(D,1) != k || size(D,2) != k
        throw(DimensionMismatch("Sizes of B ($(size(B))) and/or D ($(size(D))) are inconsistent with A ($(size(A)))"))
    end
    SymWoodbury{T, typeof(A),typeof(B),typeof(D)}(A,B,D)
end

function SymWoodbury{T}(A, B::AbstractVector{T}, D::T)
    n = size(A, 1)
    k = 1
    if size(A, 2) != n || length(B) != n
        throw(DimensionMismatch("Sizes of B ($(size(B))) and/or D ($(size(D))) are inconsistent with A ($(size(A)))"))
    end
    SymWoodbury{T,typeof(A),typeof(B),typeof(D)}(A,B,D)
end

convert{W<:Woodbury}(::Type{W}, O::SymWoodbury) = Woodbury(O.A, O.B, O.D, O.B')

inv_invD_BtX(invD, B, X) = inv(invD - B'*X);
inv_invD_BtX(invD, B::AbstractVector, X) = inv(invD - vecdot(B,X));

function calc_inv(A, B, D)
  W = inv(A);
  X = W*B;
  invD = -inv(D);
  Z = inv_invD_BtX(invD, B, X);
  SymWoodbury(W,X,Z);
end

Base.inv{T<:Any, AType<:Any, BType<:AbstractVector, DType<:Real}(O::SymWoodbury{T,AType,BType,DType}) =
  calc_inv(O.A, O.B, O.D)

Base.inv{T<:Any, AType<:Any, BType<:Any, DType<:AbstractMatrix}(O::SymWoodbury{T,AType,BType,DType}) =
  calc_inv(O.A, O.B, O.D)

# D is typically small, so this is acceptable.
Base.inv{T<:Any, AType<:Any, BType<:Any, DType<:SparseMatrixCSC}(O::SymWoodbury{T,AType,BType,DType}) =
  calc_inv(O.A, O.B, full(O.D));

\(W::SymWoodbury, R::StridedVecOrMat) = inv(W)*R

"""
    partialInv(O)

Get the factors (X,Z) in W + XZXᵀ where W + XZXᵀ = inv( A + BDBᵀ )
"""
function partialInv(O::SymWoodbury)
  X = (O.A)\O.B;
  invD = -1*inv(O.D);
  Z = inv_invD_BtX(invD, O.B, X);
  return (X,Z);
end

function liftFactorVars(A,B,D)
  A  = sparse(A)
  B  = sparse(B)
  Di = sparse(inv(D))
  n  = size(A,1)
  k  = size(B,2)
  M = [A    B   ;
       B'  -Di ];
  M = lufact(M) # ldltfact once it's available.
  return x -> (M\[x; zeros(k,1)])[1:n,:];
end

function liftFactorVars(A,B,D::SparseMatrixCSC)
  liftFactorVars(A,B,full(D))
end

# This could be optimized to avoid the extra allocation of a single element matrix
function liftFactorVars(A,B,D::Real)
  liftFactorVars(A,B,ones(1,1)*D)
end

"""
    liftFactor(A)

More stable version of inv(A).  Returns a function which computs the inverse
on evaluation, i.e. `liftFactor(A)(x)` is the same as `inv(A)*x`.
"""
liftFactor(O::SymWoodbury) = liftFactorVars(O.A,O.B,O.D)

function *{T}(O::SymWoodbury{T}, x::Union{Matrix,Vector,SubArray})
  o = O.A*x;
  plusBDBtx!(o, O.B, O.D, x)
  return o
end

function plusBDBtx!(o, B, D, x)
	o[:] = o + B*(D*(B'x))
end

plusBDBtx!(o, B::AbstractVector, D, x) = plusBDBtx!(o, reshape(B,size(B,1),1),D,x)

# Optimization - use specialized BLAS package
function plusBDBtx!(o, B::Array{Float64,2}, D, x::Array{Float64,2})
  w = D*gemm('T','N',B,x);
  gemm!('N','N',1.,B,w,1., o)
end

# Minor optimization for the rank one case
function plusBDBtx!(o, B::Array{Float64,1}, d::Real, x::Union{Array{Float64,2}, SubArray})
  if size(x,2) == 1
    axpy!(vecdot(B,x)*d, B, o)
  else
    w = d*gemm('T', 'N' ,reshape(B, size(B,1), 1),x);
    gemm!('N','N',1.,B,w,1., o)
  end
end

Base.Ac_mul_B{T}(O1::SymWoodbury{T}, x::AbstractVector{T}) = O1*x
Base.Ac_mul_B{T}(O1::SymWoodbury{T}, x::AbstractMatrix{T}) = O1*x

+(O::SymWoodbury, M::SymWoodbury)    = SymWoodbury(O.A + M.A, [O.B M.B],
                                                   cat([1,2],O.D,M.D) );
*(α::Real, O::SymWoodbury)           = SymWoodbury(α*O.A, O.B, α*O.D);
*(O::SymWoodbury, α::Real)           = SymWoodbury(α*O.A, O.B, α*O.D);
+(M::AbstractMatrix, O::SymWoodbury) = SymWoodbury(O.A + M, O.B, O.D);
+(O::SymWoodbury, M::AbstractMatrix) = SymWoodbury(O.A + M, O.B, O.D);
Base.size(M::SymWoodbury)            = size(M.A);
Base.size(M::SymWoodbury, i)         = (i == 1 || i == 2) ? size(M)[1] : 1

Base.full{T}(O::SymWoodbury{T})      = full(O.A) + O.B*O.D*O.B'
Base.copy{T}(O::SymWoodbury{T})      = SymWoodbury(copy(O.A), copy(O.B), copy(O.D))

function square(O::SymWoodbury)
  A  = O.A^2
  AB = O.A*O.B
  Z  = [(AB + O.B) (AB - O.B)]
  R  = O.D*(O.B'*O.B)*O.D/4
  D  = [ O.D/2 + R  -R
        -R          -O.D/2 + R ]
  SymWoodbury(A, Z, D)
end

"""
The product of two SymWoodbury matrices will generally be a Woodbury Matrix,
except when they are the same, i.e. the user writes A'A or A*A' or A*A.

Z(A + B*D*Bᵀ) = ZA + ZB*D*Bᵀ

This package will not support support left multiplication by a generic
matrix, to keep return types consistent.
"""
function *(O1::SymWoodbury, O2::SymWoodbury)
  if (O1 === O2)
    return square(O1)
  else
    if O1.A == O2.A && O1.B == O2.B && O1.D == O2.D
      return square(O1)
    else
      throw(MethodError("ERROR: To multiply two non-identical SymWoodbury matrices, first convert to Woodbury."))
    end
  end
end

Base.Ac_mul_B(O1::SymWoodbury, O2::SymWoodbury) = O1*O2
Base.A_mul_Bc(O1::SymWoodbury, O2::SymWoodbury) = O1*O2

conjm(O::SymWoodbury, M) = SymWoodbury(M*O.A*M', M*O.B, O.D);

Base.getindex(O::SymWoodbury, I::UnitRange, I2::UnitRange) =
  SymWoodbury(O.A[I,I], O.B[I,:], O.D);

# This is a slow hack, but generally these matrices aren't sparse.
Base.sparse(O::SymWoodbury) = sparse(full(O))

# returns a pointer to the original matrix, this is consistent with the
# behavior of Symmetric in Base.
Base.ctranspose(O::SymWoodbury) = O

Base.det(W::SymWoodbury) = det(convert(Woodbury, W))


text:
println("hello world!")

sleep(5)

ext: jl

max_stars_repo_name: jwmcgettigan/project-euler-solutions

max_issues_repo_name: jwmcgettigan/project-euler-solutions

max_forks_repo_name: jwmcgettigan/project-euler-solutions


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -0.9705716506753843 <= q <= 1.6199082189427347)

                     add_NL_constraint(m, :(softplus(0.27812717705086953 + -0.4967460082986306*softplus(-0.6740257057257448 + -0.0007917455797805673*softplus(-0.08584163188401472 + 0.6683752749724832*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*softplus(-0.7769463282105753 + 0.25325339074048747*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + -0.883304587160544*softplus(0.3709718635886192 + -0.47514710764751555*softplus(-0.08584163188401472 + 0.6683752749724832*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*softplus(-0.7769463282105753 + 0.25325339074048747*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) + softplus(0.48450040123340576 + -0.32897580116872405*softplus(-0.6740257057257448 + -0.0007917455797805673*softplus(-0.08584163188401472 + 0.6683752749724832*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*softplus(-0.7769463282105753 + 0.25325339074048747*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + 0.43090905819531766*softplus(0.3709718635886192 + -0.47514710764751555*softplus(-0.08584163188401472 + 0.6683752749724832*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*softplus(-0.7769463282105753 + 0.25325339074048747*softplus(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*softplus(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# TODO: Implement q-digest
# http://www.cs.virginia.edu/~son/cs851/papers/ucsb.sensys04.pdf


text:
# Basic test script - generates the file sample_hd_traj.csv containing one hobby drone trajectory

# Include the files
include("HobbyDroneInterface.jl")
# Set the random seed
Random.seed!(1)
# Generate the trajectory file
generate_trajectory_file(1.0, "output/hd_traj.csv")

max_issues_repo_name: sisl/em-model-hobbydrone


text:
const DISTMSG = "Right-hand side of a ~ must be subtype of Distribution or a vector of " *
    "Distributions."

const INTERNALNAMES = (:__model__, :__sampler__, :__context__, :__varinfo__, :__rng__)
const DEPRECATED_INTERNALNAMES = (:_model, :_sampler, :_context, :_varinfo, :_rng)

"""
    isassumption(expr)

Return an expression that can be evaluated to check if `expr` is an assumption in the
model.

Let `expr` be `:(x[1])`. It is an assumption in the following cases:
    1. `x` is not among the input data to the model,
    2. `x` is among the input data to the model but with a value `missing`, or
    3. `x` is among the input data to the model with a value other than missing,
       but `x[1] === missing`.

When `expr` is not an expression or symbol (i.e., a literal), this expands to `false`.
"""
function isassumption(expr::Union{Symbol, Expr})
    vn = gensym(:vn)

    return quote
        let $vn = $(varname(expr))
            # This branch should compile nicely in all cases except for partial missing data
            # For example, when `expr` is `:(x[i])` and `x isa Vector{Union{Missing, Float64}}`
            if !$(DynamicPPL.inargnames)($vn, __model__) || $(DynamicPPL.inmissings)($vn, __model__)
                true
            else
                # Evaluate the LHS
                $expr === missing
            end
        end
    end
end

# failsafe: a literal is never an assumption
isassumption(expr) = :(false)

#################
# Main Compiler #
#################

"""
    @model(expr[, warn = true])

Macro to specify a probabilistic model.

If `warn` is `true`, a warning is displayed if internal variable names are used in the model
definition.

# Examples

Model definition:

```julia
@model function model(x, y = 42)
    ...
end
```

To generate a `Model`, call `model(xvalue)` or `model(xvalue, yvalue)`.
"""
macro model(expr, warn=true)
    # include `LineNumberNode` with information about the call site in the
    # generated function for easier debugging and interpretation of error messages
    esc(model(__module__, __source__, expr, warn))
end

function model(mod, linenumbernode, expr, warn)
    modelinfo = build_model_info(expr)

    # Generate main body
    modelinfo[:body] = generate_mainbody(
        mod, modelinfo[:modeldef][:body], warn
    )

    return build_output(modelinfo, linenumbernode)
end

"""
    build_model_info(input_expr)

Builds the `model_info` dictionary from the model's expression.
"""
function build_model_info(input_expr)
    # Break up the model definition and extract its name, arguments, and function body
    modeldef = MacroTools.splitdef(input_expr)

    # Print a warning if function body of the model is empty
    warn_empty(modeldef[:body])

    ## Construct model_info dictionary

    # Shortcut if the model does not have any arguments
    if !haskey(modeldef, :args) && !haskey(modeldef, :kwargs)
        modelinfo = Dict(
            :allargs_exprs => [],
            :allargs_syms => [],
            :allargs_namedtuple => NamedTuple(),
            :defaults_namedtuple => NamedTuple(),
            :modeldef => modeldef,
        )
        return modelinfo
    end

    # Extract the positional and keyword arguments from the model definition.
    allargs = vcat(modeldef[:args], modeldef[:kwargs])

    # Split the argument expressions and the default values.
    allargs_exprs_defaults = map(allargs) do arg
        MacroTools.@match arg begin
            (x_ = val_) => (x, val)
            x_ => (x, NO_DEFAULT)
        end
    end

    # Extract the expressions of the arguments, without default values.
    allargs_exprs = first.(allargs_exprs_defaults)

    # Extract the names of the arguments.
    allargs_syms = map(allargs_exprs) do arg
        MacroTools.@match arg begin
            (::Type{T_}) | (name_::Type{T_}) => T
            name_::T_ => name
            x_ => x
        end
    end

    # Build named tuple expression of the argument symbols and variables of the same name.
    allargs_namedtuple = to_namedtuple_expr(allargs_syms)

    # Extract default values of the positional and keyword arguments.
    default_syms = []
    default_vals = []
    for (sym, (expr, val)) in zip(allargs_syms, allargs_exprs_defaults)
        if val !== NO_DEFAULT
            push!(default_syms, sym)
            push!(default_vals, val)
        end
    end

    # Build named tuple expression of the argument symbols with default values.
    defaults_namedtuple = to_namedtuple_expr(default_syms, default_vals)

    modelinfo = Dict(
        :allargs_exprs => allargs_exprs,
        :allargs_syms => allargs_syms,
        :allargs_namedtuple => allargs_namedtuple,
        :defaults_namedtuple => defaults_namedtuple,
        :modeldef => modeldef,
    )

    return modelinfo
end

"""
    generate_mainbody(mod, expr, warn)

Generate the body of the main evaluation function from expression `expr` and arguments
`args`.

If `warn` is true, a warning is displayed if internal variables are used in the model
definition.
"""
generate_mainbody(mod, expr, warn) = generate_mainbody!(mod, Symbol[], expr, warn)

generate_mainbody!(mod, found, x, warn) = x
function generate_mainbody!(mod, found, sym::Symbol, warn)
    if sym in DEPRECATED_INTERNALNAMES
        newsym = Symbol(:_, sym, :__)
        Base.depwarn(
            "internal variable `$sym` is deprecated, use `$newsym` instead.",
            :generate_mainbody!,
        )
        return generate_mainbody!(mod, found, newsym, warn)
    end

    if warn && sym in INTERNALNAMES && sym ∉ found
        @warn "you are using the internal variable `$sym`"
        push!(found, sym)
    end

    return sym
end
function generate_mainbody!(mod, found, expr::Expr, warn)
    # Do not touch interpolated expressions
    expr.head === :$ && return expr.args[1]

    # If it's a macro, we expand it
    if Meta.isexpr(expr, :macrocall)
        return generate_mainbody!(mod, found, macroexpand(mod, expr; recursive=true), warn)
    end

    # Modify dotted tilde operators.
    args_dottilde = getargs_dottilde(expr)
    if args_dottilde !== nothing
        L, R = args_dottilde
        return generate_dot_tilde(
            generate_mainbody!(mod, found, L, warn),
            generate_mainbody!(mod, found, R, warn),
        ) |> Base.remove_linenums!
    end

    # Modify tilde operators.
    args_tilde = getargs_tilde(expr)
    if args_tilde !== nothing
        L, R = args_tilde
        return generate_tilde(
            generate_mainbody!(mod, found, L, warn),
            generate_mainbody!(mod, found, R, warn),
        ) |> Base.remove_linenums!
    end

    return Expr(expr.head, map(x -> generate_mainbody!(mod, found, x, warn), expr.args)...)
end

"""
    generate_tilde(left, right)

Generate an `observe` expression for data variables and `assume` expression for parameter
variables.
"""
function generate_tilde(left, right)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds isassumption
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        return quote
            $(top...)
            $isassumption = $(DynamicPPL.isassumption(left))
            if $isassumption
                $left = $(DynamicPPL.tilde_assume)(
                    __rng__, __context__, __sampler__, $tmpright, $vn, $inds, __varinfo__
                )
            else
                $(DynamicPPL.tilde_observe)(
                    __context__, __sampler__, $tmpright, $left, $vn, $inds, __varinfo__
                )
            end
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.tilde_observe)(__context__, __sampler__, $tmpright, $left, __varinfo__)
    end
end

"""
    generate_dot_tilde(left, right)

Generate the expression that replaces `left .~ right` in the model body.
"""
function generate_dot_tilde(left, right)
    @gensym tmpright
    top = [:($tmpright = $right),
           :($tmpright isa Union{$Distribution,AbstractVector{<:$Distribution}}
             || throw(ArgumentError($DISTMSG)))]

    if left isa Symbol || left isa Expr
        @gensym out vn inds isassumption
        push!(top, :($vn = $(varname(left))), :($inds = $(vinds(left))))

        return quote
            $(top...)
            $isassumption = $(DynamicPPL.isassumption(left)) || $left === missing
            if $isassumption
                $left .= $(DynamicPPL.dot_tilde_assume)(
                    __rng__, __context__, __sampler__, $tmpright, $left, $vn, $inds, __varinfo__
                )
            else
                $(DynamicPPL.dot_tilde_observe)(
                    __context__, __sampler__, $tmpright, $left, $vn, $inds, __varinfo__
                )
            end
        end
    end

    # If the LHS is a literal, it is always an observation
    return quote
        $(top...)
        $(DynamicPPL.dot_tilde_observe)(__context__, __sampler__, $tmpright, $left, __varinfo__)
    end
end

const FloatOrArrayType = Type{<:Union{AbstractFloat, AbstractArray}}
hasmissing(T::Type{<:AbstractArray{TA}}) where {TA <: AbstractArray} = hasmissing(TA)
hasmissing(T::Type{<:AbstractArray{>:Missing}}) = true
hasmissing(T::Type) = false

"""
    build_output(modelinfo, linenumbernode)

Builds the output expression.
"""
function build_output(modelinfo, linenumbernode)
    ## Build the anonymous evaluator from the user-provided model definition.

    # Remove the name.
    evaluatordef = deepcopy(modelinfo[:modeldef])
    delete!(evaluatordef, :name)

    # Add the internal arguments to the user-specified arguments (positional + keywords).
    evaluatordef[:args] = vcat(
        [
            :(__rng__::$(Random.AbstractRNG)),
            :(__model__::$(DynamicPPL.Model)),
            :(__varinfo__::$(DynamicPPL.AbstractVarInfo)),
            :(__sampler__::$(DynamicPPL.AbstractSampler)),
            :(__context__::$(DynamicPPL.AbstractContext)),
        ],
        modelinfo[:allargs_exprs],
    )

    # Delete the keyword arguments.
    evaluatordef[:kwargs] = []

    # Replace the user-provided function body with the version created by DynamicPPL.
    evaluatordef[:body] = modelinfo[:body]

    ## Build the model function.

    # Extract the named tuple expression of all arguments and the default values.
    allargs_namedtuple = modelinfo[:allargs_namedtuple]
    defaults_namedtuple = modelinfo[:defaults_namedtuple]

    # Update the function body of the user-specified model.
    # We use a name for the anonymous evaluator that does not conflict with other variables.
    modeldef = modelinfo[:modeldef]
    @gensym evaluator
    # We use `MacroTools.@q begin ... end` instead of regular `quote ... end` to ensure
    # that no new `LineNumberNode`s are added apart from the reference `linenumbernode`
    # to the call site
    modeldef[:body] = MacroTools.@q begin
        $(linenumbernode)
        $evaluator = $(MacroTools.combinedef(evaluatordef))
        return $(DynamicPPL.Model)(
            $(QuoteNode(modeldef[:name])),
            $evaluator,
            $allargs_namedtuple,
            $defaults_namedtuple,
        )
    end

    return :($(Base).@__doc__ $(MacroTools.combinedef(modeldef)))
end

function warn_empty(body)
    if all(l -> isa(l, LineNumberNode), body.args)
        @warn("Model definition seems empty, still continue.")
    end
    return
end

"""
    matchingvalue(sampler, vi, value)

Convert the `value` to the correct type for the `sampler` and the `vi` object.
"""
function matchingvalue(sampler, vi, value)
    T = typeof(value)
    if hasmissing(T)
        _value = convert(get_matching_type(sampler, vi, T), value)
        if _value === value
            return deepcopy(_value)
        else
            return _value
        end
    else
        return value
    end
end
matchingvalue(sampler, vi, value::FloatOrArrayType) = get_matching_type(sampler, vi, value)

"""
    get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T}

Get the specialized version of type `T` for sampler `spl`.

For example, if `T === Float64` and `spl::Hamiltonian`, the matching type is
`eltype(vi[spl])`.
"""
get_matching_type(spl::AbstractSampler, vi, ::Type{T}) where {T} = T
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:Union{Missing, AbstractFloat}},
)
    return Union{Missing, floatof(eltype(vi, spl))}
end
function get_matching_type(
    spl::AbstractSampler,
    vi,
    ::Type{<:AbstractFloat},
)
    return floatof(eltype(vi, spl))
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T,N}}) where {T,N}
    return Array{get_matching_type(spl, vi, T), N}
end
function get_matching_type(spl::AbstractSampler, vi, ::Type{<:Array{T}}) where T
    return Array{get_matching_type(spl, vi, T)}
end

floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
floatof(::Type) = Real # fallback if type inference failed


text:
###############################################################################
#
# run ODE on vertical layered trace gasses
#
###############################################################################
"""
    vertical_layers!(vls::VerticalLayers{FT}, t::FT) where {FT<:AbstractFloat}

Update vertical trace gas information with time, given
- `vls` [`VerticalLayers`](@ref) type struct
- `t` Time since last update
"""
function vertical_layers!(
            vls::VerticalLayers{FT},
            t::FT
) where {FT<:AbstractFloat}
    # unpack parameters
    @unpack d_CO₂, d_layer, p_CO₂, x_layer, Δ_CO₂ = vls;
    bc = DirichletBC(FT(50), FT(30));
    u0 = vls.p_CO₂;

    # define step function, move out later
    # very memory extensive for now
    step(u,p,t) = Δ_CO₂ * bc * u .* d_CO₂ ./ x_layer ./ d_layer;
    prob = ODEProblem(step, u0, (FT(0), t));
    alg  = KenCarp4();
    sol  = solve(prob, alg);

    @show sol[end];

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#=
Description
    Module for multi-agent UAV algorithm and simulations.
=#

module Multiagent

using DoubleUAVs, POMDPs, UAVs, GridInterpolations, QMDP, Interact,
      PGFPlots, JLD

export read_alpha, get_grid, viz_policy, randuavs, stress!, plot_trajs,
       maxsum, maxmin, utm, uncrd, coord, naive, distr, bulk_test

const ALPHA_FILE = "../../data/alpha.jld"
const ALPHA_VARIABLE = "alpha"

const X = 1; const XDOT = 2     # x-direction
const Y = 3; const YDOT = 4     # y-direction
const P = 5; const PDOT = 6     # heading
const B = 7; const BDOT = 8     # bank angle

const PROBLEM = DoubleUAV()

const PX = 1         # state variable indexing
const PY = 2
const PP = 3
const PV1 = 4
const PV2 = 5
const NPSTATES = size(PROBLEM.pomdp.states, 1)
const NSTATES = size(PROBLEM.pomdp.states, 2)

const COC = -1
const INDIV_ACTIONS = deg2rad([-20, -10, 0, 10, 20, rad2deg(COC)])
const NINDIV_ACTIONS = length(INDIV_ACTIONS)
actions = zeros(2, NINDIV_ACTIONS^2)
iaction = 1
for ia2 = 1:NINDIV_ACTIONS
    for ia1 = 1:NINDIV_ACTIONS
        actions[1, iaction] = INDIV_ACTIONS[ia1]
        actions[2, iaction] = INDIV_ACTIONS[ia2]
        iaction += 1
    end # for ia1
end # for ia2
const ACTIONS = actions
const NACTIONS = size(ACTIONS, 2)

const DT = 5.0              # [s]
const DTI = 1.0             # [s]
const G = 9.8               # [m/s^2]

const DIM_ACTIONS = 2
const DIM_POMDP_STATES = 5
const DIM_STATES = 8        # TYPE_SIMPLE state dimensionality
const TERM_STATE = 1e5 * ones(DIM_POMDP_STATES)

const XMIN = -2e3
const XMAX = 2e3
const YMIN = -2e3
const YMAX = 2e3
const PMIN = 0.0
const PMAX = 2 * pi
const VMIN = 10
const VMAX = 10#20

const WN = 0.20             # [rad/s]
const MIN_SEP = 500.0       # [m]
const MIN_SEP_INIT = 600.0  # [m]

const FEAS_PERCENTILE = 0.8
const CRASH_REWARD = -5e2

const NT = 100              # []
const MAX_ITER_JESP = 100

const MIN_RAD = 2000        # [m]
const MAX_RAD = 3000        # [m]
const RANGE_RAD = MAX_RAD - MIN_RAD

const MU = 0.0              # [rad] gaussian noise mean
const SIGMA = deg2rad(2)    # [rad] gaussian noise variance

const CONSENSUS_REWARD = 0.5

const SIGMA_V = 1.0         # [m/s]
const SIGMA_B = deg2rad(2)  # [rad]
const SIGMA_P = deg2rad(2)  # [rad]
const SIGMA_XY = 25         # [m]

const LARGE_NUMBER = int64(1e9)

function read_alpha(infile::ASCIIString=ALPHA_FILE)
    return load(infile, ALPHA_VARIABLE)
end # function read_alpha

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    if length(v1.rowval) < length(v2.rowval)
        for iv1 in 1:length(v1.rowval)
            sum += v1.nzval[iv1] * v2[v1.rowval[iv1]]
        end # for iv1
    else
        for iv2 in 1:length(v2.rowval)
            sum += v2.nzval[iv2] * v1[v2.rowval[iv2]]
        end # for iv2
    end # if
    return sum
end # function spdot

function spdot(v1::Vector{Float64},
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    for iv2 in 1:length(v2.rowval)
        sum += v2.nzval[iv2] * v1[v2.rowval[iv2]]
    end # for iv2
    return sum
end # function spdot

function spdot(m1::Matrix{Float64}, icol::Int64,
               v2::SparseMatrixCSC{Float64, Int64})
    sum = 0.0
    for iv2 in 1:length(v2.rowval)
        sum += v2.nzval[iv2] * m1[v2.rowval[iv2], icol]
    end # for iv2
    return sum
end # function spdot

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               v2::Vector{Float64})
    sum = 0.0
    for iv1 in 1:length(v1.rowval)
        sum += v1.nzval[iv1] * v2[v1.rowval[iv1]]
    end # for iv1
    return sum
end # function spdot

function spdot(v1::SparseMatrixCSC{Float64, Int64},
               m2::Matrix{Float64}, icol::Int64)
    sum = 0.0
    for iv1 in 1:length(v1.rowval)
        sum += v1.nzval[iv1] * m2[v1.rowval[iv1], icol]
    end # for iv1
    return sum
end # function spdot

function randuavs(nuav::Int64)
    return [UAV(VMIN + (VMAX - VMIN) * rand()) for iuav = 1:nuav]
end # function randuavs

function randactions(nuav::Int64)
    return [0.0 for iuav = 1:nuav] # for standardization
end # function randactions

function get_grid(states::Matrix{Float64})
    g = RectangleGrid(sort(unique(vec(states[PX, 1:end - 1]))),
                      sort(unique(vec(states[PY, 1:end - 1]))),
                      sort(unique(vec(states[PP, 1:end - 1]))),
                      sort(unique(vec(states[PV1, 1:end - 1]))),
                      sort(unique(vec(states[PV2, 1:end - 1]))))
    return g
end # function get_grid

function get_pomdp_state(s1::Vector{Float64}, s2::Vector{Float64})
    dx = s2[X] - s1[X]
    dy = s2[Y] - s1[Y]

    xr = dx * cos(s1[P]) + dy * sin(s1[P])
    yr = -dx * sin(s1[P]) + dy * cos(s1[P])
    pr = s2[P] - s1[P]
    v1 = sqrt(s1[XDOT]^2 + s1[YDOT]^2)
    v2 = sqrt(s2[XDOT]^2 + s2[YDOT]^2)

    if xr < XMIN || xr > XMAX ||
       yr < YMIN || yr > YMAX
        return TERM_STATE
    else
        return [xr, yr, pr, v1, v2]
    end # if
end # function get_pomdp_state

function get_belief(pstate::Vector{Float64}, grid::RectangleGrid)
    belief = spzeros(NSTATES, 1)
    if pstate == TERM_STATE
        belief[end] = 1.0
    else
        indices, weights = interpolants(grid, pstate)
        belief[indices] = weights
    end # if
    return belief
end # function get_belief

function get_belief!(pstate::Vector{Float64}, grid::RectangleGrid,
                     beliefs::SparseMatrixCSC{Float64, Int64},
                     uavIdx::Int64)
    if pstate == TERM_STATE
        beliefs[end, uavIdx] = 1.0
    else
        indices, weights = interpolants(grid, pstate)
        beliefs[indices, uavIdx] = weights
    end # if
end # function get_belief

function action_idx(a1::Float64, a2::Float64)
    idx1 = 0
    idx2 = 0
    for iaction in 1:NINDIV_ACTIONS
        if INDIV_ACTIONS[iaction] == a1
            idx1 = iaction
        end # if
        if INDIV_ACTIONS[iaction] == a2
            idx2 = iaction
        end # if
        if idx1 != 0 && idx2 != 0
            break
        end # if
    end # for iaction
    return sub2ind([NINDIV_ACTIONS, NINDIV_ACTIONS], idx1, idx2)
end # function action_idx

function and(a1::BitArray{1}, a2::BitArray{1})
    @assert length(a1) == length(a2)
    for i = 1:length(a1)
        if a2[i] == false
            a1[i] = false
        end # if
    end # for i
    return a1
end # function and

function consensus(action::Float64, suggest::Symbol)
    if (suggest == :left && action > 0) ||
       (suggest == :right && action < 0) ||
       (suggest == :straight && action == 0)
        return true
    else # suggest == :anything or no consensus
        return false
    end
end # function consensus

function maxsum(ownship::Int64, uavs::Vector{UAV},
                actions::Vector{Float64}, alpha::Matrix{Float64},
                grid::RectangleGrid, suggest::Symbol=:anything)
    #=
    Returns the max-sum ownship action and utility value based
    on utility fusion.
    =#
    nuav = length(uavs)
    utils = zeros(NINDIV_ACTIONS)

    # compute belief states for each uav pair
    beliefs = spzeros(NSTATES, nuav)
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        state = get_pomdp_state(uavs[ownship].state, uavs[iu].state)
        get_belief!(state, grid, beliefs, iu)
    end # for iu

    # find action corresponding to max-sum utility
    for iaction = 1:NINDIV_ACTIONS
        # compute sum of utility corresponding to ownship action
        for iu = 1:nuav
            if iu == ownship
                continue
            end # if
            ia = action_idx(INDIV_ACTIONS[iaction], actions[iu])
            if beliefs[end, iu] == 1
                if INDIV_ACTIONS[iaction] == COC
                    utils[iaction] += 0.0
                else
                    utils[iaction] += -0.5
                end # if
            else
                utils[iaction] += spdot(alpha, ia, beliefs[:, iu])
            end # if
        end # for iu

        # add online cost due to suggestion
        if consensus(INDIV_ACTIONS[iaction], suggest)
            utils[iaction] += CONSENSUS_REWARD * abs(utils[iaction])
        end # if
    end # for iaction

    # compute feasible range of actions
    bestAction = INDIV_ACTIONS[indmax(utils)]
    bestUtil = maximum(utils)

    return bestAction, bestUtil
end # function maxsum

function maxmin(ownship::Int64, uavs::Vector{UAV},
                actions::Vector{Float64}, alpha::Matrix{Float64},
                grid::RectangleGrid, suggest::Symbol=:anything)
    #=
    Returns the max-min ownship action and utility value based
    on utility fusion.
    =#
    nuav = length(uavs)
    utils = Inf * ones(NINDIV_ACTIONS)

    # compute belief states for each uav pair
    beliefs = spzeros(NSTATES, nuav)
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        state = get_pomdp_state(uavs[ownship].state, uavs[iu].state)
        get_belief!(state, grid, beliefs, iu)
    end # for iu

    # find action corresponding to max-min utility
    for iaction = 1:NINDIV_ACTIONS
        # compute minimum utility corresponding to ownship action
        for iu = 1:nuav
            if iu == ownship
                continue
            end # if
            aIdx = action_idx(INDIV_ACTIONS[iaction], actions[iu])
            if beliefs[end, iu] == 1
                if INDIV_ACTIONS[iaction] == COC
                    util = 0.0
                else
                    util = -0.5
                end # if
            else
                util = spdot(alpha, aIdx, beliefs[:, iu])
            end # if
            if utils[iaction] > util
                utils[iaction] = util
            end # if
        end # for iu

        # add online cost due to suggestion
        if consensus(INDIV_ACTIONS[iaction], suggest)
            utils[iaction] += CONSENSUS_REWARD * abs(utils[iaction])
        end # if
    end # for iaction

    # compute feasible range of actions
    bestAction = INDIV_ACTIONS[indmax(utils)]
    bestUtil = maximum(utils)

    return bestAction, bestUtil
end # function maxmin

function maxsum(iaction::Int64, ownship::Int64, uavs::Vector{UAV},
                actions::Vector{Float64}, alpha::Matrix{Float64},
                grid::RectangleGrid, suggest::Symbol=:anything)
    #=
    Returns the max-sum ownship utility value based on utility
    fusion.
    =#
    nuav = length(uavs)
    util = 0

    # compute belief states for each uav pair
    beliefs = spzeros(NSTATES, nuav)
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        state = get_pomdp_state(uavs[ownship].state, uavs[iu].state)
        get_belief!(state, grid, beliefs, iu)
    end # for iu

    # compute sum of utility corresponding to ownship action
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        ia = action_idx(INDIV_ACTIONS[iaction], actions[iu])
        if beliefs[end, iu] == 1
            if INDIV_ACTIONS[iaction] == COC
                util += 0.0
            else
                util += -0.5
            end # if
        else
            util += spdot(alpha, ia, beliefs[:, iu])
        end # if
    end # for iu

    # add online cost due to suggestion
    if consensus(INDIV_ACTIONS[iaction], suggest)
        util += CONSENSUS_REWARD * abs(util)
    end # if

    return util
end # function maxsum

function maxmin(iaction::Int64, ownship::Int64, uavs::Vector{UAV},
                actions::Vector{Float64}, alpha::Matrix{Float64},
                grid::RectangleGrid, suggest::Symbol=:anything)
    #=
    Returns the max-min ownship utility value based on utility
    fusion.
    =#
    nuav = length(uavs)
    min_util = Inf

    # compute belief states for each uav pair
    beliefs = spzeros(NSTATES, nuav)
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        state = get_pomdp_state(uavs[ownship].state, uavs[iu].state)
        get_belief!(state, grid, beliefs, iu)
    end # for iu

    # compute minimum utility corresponding to ownship action
    for iu = 1:nuav
        if iu == ownship
            continue
        end # if
        aIdx = action_idx(INDIV_ACTIONS[iaction], actions[iu])
        if beliefs[end, iu] == 1
            if INDIV_ACTIONS[iaction] == COC
                util = 0.0
            else
                util = -0.5
            end # if
        else
            util = spdot(alpha, aIdx, beliefs[:, iu])
        end # if
        if min_util > util
            min_util = util
        end # if
    end # for iu

    # add online cost due to suggestion
    if consensus(INDIV_ACTIONS[iaction], suggest)
        utils[iaction] += CONSENSUS_REWARD * abs(utils[iaction])
    end # if

    return min_util
end # function maxmin

function get_distance(u1::Int64, u2::Int64, uavs::Vector{UAV})
    return norm([uavs[u1].state[X] - uavs[u2].state[X],
                 uavs[u1].state[Y] - uavs[u2].state[Y]])
end # function get_distance

function get_distance(u1::Int64, u2::Int64, uavs::Vector{UAV},
                      quavs::Vector{UAV})
    return norm([uavs[u1].state[X] - quavs[u2].state[X],
                 uavs[u1].state[Y] - quavs[u2].state[Y]])
end # function get_distance

function greedy(uavs::Vector{UAV}, alpha::Matrix{Float64},
                grid::RectangleGrid, utilFn::Function)
    #=
    Returns the uav actions solution via a greedy search
    for the best policy against the closest uav, assuming
    all other intruders are white noise intruders.
    =#
    nuav = length(uavs)
    closest = zeros(Int64, nuav)
    actions = COC * ones(nuav)
    # ranges = Array(Vector{Float64}, nuav)
    utils = zeros(nuav)

    # find closest other uav
    for iu = nuav:-1:2
        distances = Inf * ones(iu - 1)
        for iu2 = 1:iu - 1
            distances[iu2] = get_distance(iu, iu2, uavs)
        end # for iu2
        closest[iu] = indmin(distances)
        closest[closest[iu]] = iu
    end # for iu

    # compute best actions
    for iu = 1:nuav
        action, util =
            utilFn(1, [uavs[iu], uavs[closest[iu]]],
                   COC * ones(2), alpha, grid)
        actions[iu] = action
        utils[iu] = util
    end # for iu
    return actions, utils
end # function greedy

function greedy(uavs::Vector{UAV}, alpha::Matrix{Float64},
                grid::RectangleGrid, utilFn::Function,
                quavs::Vector{UAV})
    #=
    Returns the uav actions solution via a greedy search
    for the best policy against the closest uav, assuming
    all other intruders are white noise intruders.
    =#
    nuav = length(uavs)
    closest = zeros(Int64, nuav)
    actions = COC * ones(nuav)
    # ranges = Array(Vector{Float64}, nuav)
    utils = zeros(nuav)

    # find closest other uav
    for iu = nuav:-1:2
        distances = Inf * ones(iu - 1)
        for iu2 = 1:iu - 1
            distances[iu2] = get_distance(iu, iu2, uavs, quavs)
        end # for iu2
        closest[iu] = indmin(distances)
        closest[closest[iu]] = iu
    end # for iu

    # compute best actions
    for iu = 1:nuav
        action, util =
            utilFn(1, [uavs[iu], quavs[closest[iu]]],
                   COC * ones(2), alpha, grid)
        actions[iu] = action
        utils[iu] = util
    end # for iu
    return actions, utils
end # function greedy

function jesp(uavs::Vector{UAV}, alpha::Matrix{Float64},
              grid::RectangleGrid, utilFn::Function)
    #=
    Returns the utility fusion uav actions solution via
    the joint equilibrium based search for policy method.
    =#
    nuav = length(uavs)

    actions = COC * ones(nuav)

    utils = -Inf * ones(nuav)

    bestUtil = -Inf
    bestActions = copy(actions)

    niter = 0
    solutionImproving = true

    while niter < MAX_ITER_JESP && solutionImproving
        solutionImproving = false

        jespOrder = randperm(nuav)
        for iu = jespOrder
            action, util =
                utilFn(iu, uavs, actions, alpha, grid)
            utils[iu] = util
            actions[iu] = action

            if bestUtil < sum(utils)
                bestUtil = sum(utils)
                bestActions = copy(actions)
                solutionImproving = true
            end # if
        end

        niter += 1
        if niter == MAX_ITER_JESP
            print("Max iter reached, check solution\n")
        end # if
    end # while

    return bestActions, bestUtil, niter
end # function jesp

function fuse_uncrd(uavs::Vector{UAV}, alpha::Matrix{Float64},
                    grid::RectangleGrid, utilFn::Function)
    #=
    Returns the utility fusion uav actions solution via
    a greedy search for policy method, assuming all other
    intruders are white noise intruders.
    =#
    nuav = length(uavs)
    dummyActions = COC * ones(nuav)
    actions = COC * ones(nuav)
    utils = zeros(nuav)

    for iu = 1:nuav
        action, util =
            utilFn(iu, uavs, dummyActions, alpha, grid)
        actions[iu] = action
        utils[iu] = util
    end # for iu
    return actions, utils
end # function fuse_uncrd

function fuse_uncrd(quavs::Vector{UAV}, alpha::Matrix{Float64},
                    grid::RectangleGrid, utilFn::Function)
    #=
    Returns the utility fusion uav actions solution via
    a greedy search for policy method, assuming all other
    intruders are white noise intruders.
    =#
    nuav = length(quavs)
    dummyActions = COC * ones(nuav)
    actions = COC * ones(nuav)
    utils = zeros(nuav)

    for iu = 1:nuav
        iuavs = [quavs[iu], quavs[1:nuav .!= iu]]
        action, util =
            utilFn(1, iuavs, dummyActions, alpha, grid)
        actions[iu] = action
        utils[iu] = util
    end # for iu
    return actions, utils
end # function fuse_uncrd

function fuse_coord(quavs::Vector{UAV}, alpha::Matrix{Float64},
                    grid::RectangleGrid, utilFn::Function)
    #=
    Returns the utility fusion uav actions solution via
    a search for the best policy based on coordination messages
    from other uavs. Note that each uav still assumes that
    all other uavs are white noise intruders.
    =#
    nuav = length(quavs)
    actions = COC * ones(nuav)
    utils = zeros(nuav)

    # compute what each uav thinks is best for all uavs
    sactions = COC * ones(nuav, nuav)
    for iu = 1:nuav
        jactions, _ = jesp(quavs, alpha, grid, utilFn)
        sactions[:, iu] = jactions
    end # for iu

    # generate what each uav guesses it should do
    gactions = diag(sactions)

    # compute best action with suggested turning sense
    for iu = 1:nuav
        # average coordination messages to get suggested sense
        jaction = sum(sactions[iu, 1:nuav .!= iu])
        if jaction > 0
            suggest = :left
        elseif jaction < 0
            suggest = :right
        else # jaction == 0
            suggest = :straight
        end # if

        # compute best action
        action, util =
            utilFn(iu, quavs, gactions, alpha, grid, suggest)
        actions[iu] = action
        utils[iu] = util
    end # for iu
    return actions, utils
end # function fuse_coord

function fuse_distr(quavs::Vector{UAV}, alpha::Matrix{Float64},
                    grid::RectangleGrid, utilFn::Function)
    nuav = length(quavs)
    actions = COC * ones(nuav)
    pactions = Inf * ones(nuav)
    utils = zeros(nuav)

    niter = 0
    while pactions != actions && niter < MAX_ITER_JESP
        niter += 1
        pactions = actions
        # compute each uav's best greedy action in parallel
        for iu = 1:nuav
            action, util =
                utilFn(iu, quavs, pactions, alpha, grid)
            actions[iu] = action
            utils[iu] = util
        end # for iu
    end # while
    return actions, utils
end # function fuse_distr

function noisify(uavs::Vector{UAV})
    #=
    Adds noise to the state variables of uavs.
    =#
    quavs = deepcopy(uavs)
    nuav = length(quavs)
    for iu = 1:nuav
        quavs[iu].state[X] += SIGMA_XY * randn()
        quavs[iu].state[Y] += SIGMA_XY * randn()
        quavs[iu].state[P] += SIGMA_P * randn()
        speed = SIGMA_V * randn() +
                sqrt(quavs[iu].state[XDOT]^2 +
                     quavs[iu].state[YDOT]^2)
        quavs[iu].state[XDOT] = speed * cos(quavs[iu].state[P])
        quavs[iu].state[YDOT] = speed * sin(quavs[iu].state[P])
    end # for iu
    return quavs
end # function noisify

function utm(uavs::Vector{UAV}, alpha::Matrix{Float64},
             grid::RectangleGrid, utilFn::Function, noise::Bool)
    #=
    Centralized controller computes the best action for all planes
    based on received measurements from each UAV and sends advisories.
    For UTM-like centralized controller, or ACAS-like coordinated
    controller with perfect information (measurements + messages).
    =#
    quavs = noisify(uavs)
    actions, _ = jesp(quavs, alpha, grid, utilFn)
    if noise
        addnoise!(actions)
    end # if
    return actions
end # function utm

function uncrd(uavs::Vector{UAV}, alpha::Matrix{Float64},
               grid::RectangleGrid, utilFn::Function, noise::Bool)
    #=
    Each UAV plans its best action based on white-noise assumption
    for all other uavs using utility fusion (no coordination). For
    ACAS-like uncoordinated controller with quantized state messages.
    =#
    quavs = noisify(uavs)
    actions, _ = fuse_uncrd(quavs, alpha, grid, utilFn)
    if noise
        addnoise!(actions)
    end # if
    return actions
end # function uncrd

function coord(uavs::Vector{UAV}, alpha::Matrix{Float64},
               grid::RectangleGrid, utilFn::Function, noise::Bool)
    #=
    Each UAV plans its best action based on coordination messages
    from all other uavs using utility fusion. For ACAS-like
    coordinated controller with quantized state messages.
    =#
    quavs = noisify(uavs)
    actions, _ = fuse_coord(quavs, alpha, grid, utilFn)
    if noise
        addnoise!(actions)
    end # if
    return actions
end # function uncrd

function naive(uavs::Vector{UAV}, alpha::Matrix{Float64},
               grid::RectangleGrid, utilFn::Function, noise::Bool)
    #=
    Each UAV plans its best action based on white-noise assumption
    for the closest other uav (no utility fusion; greedy takes
    care of this, even though the api is the same as acas and utm).
    For either UTM or ACAS-like baseline controller.
    =#
    actions, _ = greedy(uavs, alpha, grid, utilFn)
    if noise
        addnoise!(actions)
    end # if
    return actions
end # function naive

function distr(uavs::Vector{UAV}, alpha::Matrix{Float64},
               grid::RectangleGrid, utilFn::Function, noise::Bool)
    #=
    Each UAV greedily selects its own action in a distributed
    fashion at each iteration and passes its solution to server.
    =#
    quavs = noisify(uavs)
    actions, _ = fuse_distr(quavs, alpha, grid, utilFn)
    if noise
        addnoise!(actions)
    end # if
    return actions
end # function distr

function set_scenario!(iu::Int64, uavs::Vector{UAV},
                       scenario::Vector{Float64})
    uavs[iu].state[X] = scenario[1]
    uavs[iu].state[Y] = scenario[2]
    uavs[iu].state[P] = scenario[3]
    uavs[iu].state[B] = 0.0
    uavs[iu].state[XDOT] = uavs[iu].speed * cos(uavs[iu].state[P])
    uavs[iu].state[YDOT] = uavs[iu].speed * sin(uavs[iu].state[P])
    uavs[iu].state[PDOT] = 0.0
    uavs[iu].state[BDOT] = 0.0
end # function set_scenario

function viz_policy(alpha::Matrix{Float64}, grid::RectangleGrid, nbins::Int64=100)
    uavs = randuavs(3)
    set_scenario!(2, uavs, [1250, 600, deg2rad(180)])
    set_scenario!(3, uavs, [1250, -600, deg2rad(180)])

    utilFnRange = [maxmin, maxsum]
    prange = 0#:60:360

    @manipulate for iutilFn = 1:length(utilFnRange), p = prange
        utilFn = utilFnRange[iutilFn]
        function getmap1(x::Float64, y::Float64)
            set_scenario!(1, uavs, [x, y, deg2rad(p)])
            actions, _ = jesp(uavs, alpha, grid, utilFn)
            if abs(actions[1]) > 0.5
                return 2.0
            end # if
            return rad2deg(actions[1])
        end # function getmap1

        function getmap2(x::Float64, y::Float64)
            set_scenario!(1, uavs, [x, y, deg2rad(p)])
            actions, _ = jesp(uavs, alpha, grid, utilFn)
            if abs(actions[2]) > 0.5
                return 2.0
            end # if
            return rad2deg(actions[2])
        end # function getmap2

        function getmap3(x::Float64, y::Float64)
            set_scenario!(1, uavs, [x, y, deg2rad(p)])
            actions, _ = jesp(uavs, alpha, grid, utilFn)
            if abs(actions[3]) > 0.5
                return 2.0
            end # if
            return rad2deg(actions[3])
        end # function getmap3

        g = GroupPlot(3, 1, groupStyle="horizontal sep=1.5cm")
        push!(g, Axis([
            Plots.Image(getmap1, (-2500, 2500), (-2500, 2500),
                        zmin = -20, zmax = 20, xbins = nbins, ybins = nbins,
                        colormap = ColorMaps.Named("RdBu"),
                        colorbar = false),
    Plots.Node(L">", 1250, 600, style="rotate=180,font=\\huge"),
    Plots.Node(L">", 1250, -600, style="rotate=180,font=\\huge")],
                      width="9cm", height="8cm",
                      xlabel="x (m)", ylabel="y (m)",
                      title="First aircraft action",
                      style="xtick={-3000,-2000,...,3000}"))
        push!(g, Axis([
            Plots.Image(getmap2, (-2500, 2500), (-2500, 2500),
                        zmin = -20, zmax = 20, xbins = nbins, ybins = nbins,
                        colormap = ColorMaps.Named("RdBu"),
                        colorbar = false),
    Plots.Node(L">", 1250, 600, style="rotate=180,font=\\huge"),
    Plots.Node(L">", 1250, -600, style="rotate=180,font=\\huge")],
                      width="9cm", height="8cm",
                      xlabel="x (m)" ,
                      title="Second aircraft action",
                      style="xtick={-3000,-2000,...,3000}"))
        push!(g, Axis([
            Plots.Image(getmap3, (-2500, 2500), (-2500, 2500),
                        zmin = -20, zmax = 20, xbins = nbins, ybins = nbins,
                        colormap = ColorMaps.Named("RdBu")),
    Plots.Node(L">", 1250, 600, style="rotate=180,font=\\huge"),
    Plots.Node(L">", 1250, -600, style="rotate=180,font=\\huge")],
            width="9cm", height="8cm",
                      xlabel="x (m)",
                      title="Third aircraft action",
                      style="xtick={-3000,-2000,...,3000}"))
        g
    end # for utilFn, p
end # function viz_policy

function norm_angle(angle::Float64)
    return ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
end # function norm_angle

function rand_uavpos()
    # randomly generate position in polar coordinates
    r = MIN_RAD + RANGE_RAD * rand()
    p = norm_angle(2 * pi * rand())

    # convert into cartesian coordinates
    x = r * cos(p)
    y = r * sin(p)
    head = norm_angle(p + pi)
    return [x, y, head]
end # function rand_uavpos

function initconflict(x::Float64, y::Float64, uavs::Vector{UAV})
    for iu = 1:length(uavs)
        if (x - uavs[iu].state[X])^2 +
           (y - uavs[iu].state[Y])^2 <
           (MIN_SEP_INIT)^2
            return true
        end # if
    end # for iu
    return false
end # function inconflict

function inconflict(uav::UAV, uavs::Vector{UAV})
    for iu = 1:length(uavs)
        if MIN_SEP^2 > (uav.state[X] - uavs[iu].state[X])^2 +
                       (uav.state[Y] - uavs[iu].state[Y])^2
            return true
        end # if
    end # for iu
    return false
end # function inconflict

function inconflict(uav1::UAV, uav2::UAV)
    if MIN_SEP^2 > (uav1.state[X] - uav2.state[X])^2 +
                   (uav1.state[Y] - uav2.state[Y])^2
        return true
    else
        return false
    end # if
end # function inconflict

function init_test!(uavs::Vector{UAV})
    for iu = 1:length(uavs)
        newpos = rand_uavpos()
        while initconflict(newpos[1], newpos[2], uavs[1:iu - 1])
            newpos = rand_uavpos()
        end # while
        set_scenario!(iu, uavs, newpos)
    end # for iu
end # function init_test!

function plot_trajs(alg::Function, util::Function, alpha::Matrix{Float64}, g::RectangleGrid)
    trajs, n, c, t = stress!(alg, util, true, NT, randuavs(4), alpha, g)
    print("number of collisions = ", sum(c), ", average decision time = ")
    @printf("%.3e ms", t * 1000)
    PGFPlots.Axis([Plots.Linear(vec(trajs[1, X, :]), vec(trajs[1, Y, :]), mark="none", style="blue,->,smooth,thick"),
                   Plots.Linear(vec(trajs[2, X, :]), vec(trajs[2, Y, :]), mark="none", style="red,->,smooth,thick"),
                   Plots.Linear(vec(trajs[3, X, :]), vec(trajs[3, Y, :]), mark="none", style="green,->,smooth,thick"),
                   Plots.Linear(vec(trajs[4, X, :]), vec(trajs[4, Y, :]), mark="none", style="black,->,smooth,thick"),
                   Plots.Linear([trajs[1, X, 1]], [trajs[1, Y, 1]], mark="*", style="blue"),
                   Plots.Linear([trajs[2, X, 1]], [trajs[2, Y, 1]], mark="*", style="red"),
                   Plots.Linear([trajs[3, X, 1]], [trajs[3, Y, 1]], mark="*", style="green"),
                   Plots.Linear([trajs[4, X, 1]], [trajs[4, Y, 1]], mark="*", style="black")],
                  xlabel="x (m)", ylabel="y (m)", width="12cm", height="12cm",
                  style="axis equal,xtick={-4000,-2000,...,4000}, ytick={-4000,-2000,...,4000}")
end # function plot_traj

function pid_policy(bankAngle::Float64)
    function policy(t::Float64, state::Vector{Float64})
        return 2 * WN * (-state[BDOT]) + WN^2 * (bankAngle - state[B])
    end # function policy
    return policy
end # function pid_policy

function init_trajs(uavs::Vector{UAV}, nt::Int64)
    nuavs = length(uavs)
    trajs = zeros(nuavs, DIM_STATES, nt * int64(DT / DTI) + 1)
    for iu = 1:nuavs
        trajs[iu, :, 1] = uavs[iu].state
    end # for iu
    return trajs
end # function init_trajs

function gen_traj(uav::UAV, policy::Function)
    _, traj = generate_trajectory(uav, policy, 1, DTI)
    return traj
end # function gen_traj

function sim_trajs!(uavs::Vector{UAV}, actions::Vector{Float64},
                    dtcounts::Vector{Int64}, trajs::Array{Float64,3},
                    conflicts::BitArray{2})
    nuavs = length(uavs)
    nlms = 0

    # generate pid policies
    pids = Array(Function, nuavs)
    for iu = 1:nuavs
        # correct COC to 0 degrees nominal turn
        if abs(actions[iu]) > 0.5
            actions[iu] = 0.0
        end # if
        # add noise to bank angle commanded
        pids[iu] = pid_policy(actions[iu] + SIGMA_B * randn())
    end # for iu

    # simulate trajectories
    for idt = 1:int64(DT / DTI)
        # update state for each uav
        for iu = 1:nuavs
            dtcounts[iu] += 1
            traj = gen_traj(uavs[iu], pids[iu])
            trajs[iu, :, dtcounts[iu]] = traj[:, end]
            set_state!(uavs[iu], traj[:, end])
        end # for iu

        # check conflict for each uav pair
        for iu1 = 1:nuavs
            for iu2 = iu1 + 1:nuavs
                if inconflict(uavs[iu1], uavs[iu2])
                    nlms += 1
                    conflicts[iu1, iu2] = true
                end # if
            end # for iu2
        end # for iu1
    end # for idt

    return nlms
end # function sim_trajs!

function gaussnoise(mu::Float64=MU, sigma::Float64=SIGMA)
    return mu + sigma * randn()
end # function gaussnoise

function addnoise!(actions::Vector{Float64}, noise::Function=gaussnoise)
    for iu = 1:length(actions)
        actions[iu] += noise()
    end # for iu
end # function addnoise!

function compute_alerts!(alerts::BitArray{1}, actions::Vector{Float64})
    alertcount = 0
    for iaction = 1:length(actions)
        if abs(actions[iaction]) > 0.5
            alerts[iaction] = true
            alertcount += 1
        end # if
    end # for action
    return alertcount
end # function compute_alerts!

function stress!(policy::Function, utilFn::Function, noise::Bool, nt::Int64,
                 uavs::Vector{UAV}, alpha::Matrix{Float64}, grid::RectangleGrid)
    init_test!(uavs)
    nuav = length(uavs)
    trajs = init_trajs(uavs, nt)
    nlms = 0
    dtcounts = ones(Int64, nuav)
    conflicts = falses(nuav, nuav)
    avgtime = 0
    alerts = falses(nuav)
    alertcounts = 0

    for it = 1:nt
        tic()
        actions = policy(uavs, alpha, grid, utilFn, noise)
        alertcounts += compute_alerts!(alerts, actions)
        avgtime += toq()
        nlms += sim_trajs!(uavs, actions, dtcounts, trajs, conflicts)
    end # for it
    avgtime /= nt

    return trajs, nlms, conflicts, avgtime, alerts, alertcounts
end # function stress!

function bulk_test(nuavs::UnitRange{Int64}, nsim::Int64,
                   alpha::Matrix{Float64}, g::RectangleGrid,
                   nbulk::Int64=1)
    for ibulk = 1:nbulk
        print("beginning stress tests...\n")
        maxminNLMS = zeros(Int64, 5, length(nuavs))
        maxsumNLMS = zeros(Int64, 5, length(nuavs))
        maxminNLMSBool = zeros(Int64, 5, length(nuavs))
        maxsumNLMSBool = zeros(Int64, 5, length(nuavs))
        maxminTimes = zeros(Float64, 5, length(nuavs))
        maxsumTimes = zeros(Float64, 5, length(nuavs))
        maxminAlerts = zeros(Int64, 5, length(nuavs))
        maxsumAlerts = zeros(Int64, 5, length(nuavs))
        maxminAlertsCount = zeros(Int64, 5, length(nuavs))
        maxsumAlertsCount = zeros(Int64, 5, length(nuavs))

        for inu = 1:length(nuavs)
            nu = nuavs[inu]
            uavs = randuavs(nu)
            tic()
            for isim = 1:nsim
                _, n1, c1, t1, a1, as1 = stress!(utm, maxmin, true, NT, uavs, alpha, g)
                _, n2, c2, t2, a2, as2 = stress!(uncrd, maxmin, true, NT, uavs, alpha, g)
                _, n3, c3, t3, a3, as3 = stress!(coord, maxmin, true, NT, uavs, alpha, g)
                _, n4, c4, t4, a4, as4 = stress!(naive, maxmin, true, NT, uavs, alpha, g)
                _, n5, c5, t5, a5, as5 = stress!(distr, maxmin, true, NT, uavs, alpha, g)
                maxminNLMS[:, inu] += [n1, n2, n3, n4, n5]
                maxminNLMSBool[:, inu] += [sum(c1), sum(c2), sum(c3), sum(c4), sum(c5)]
                maxminTimes[:, inu] += [t1, t2, t3, t4, t5]
                maxminAlerts[:, inu] += [sum(a1), sum(a2), sum(a3), sum(a4), sum(a5)]
                maxminAlertsCount[:, inu] += [as1, as2, as3, as4, as5]
            end # for isim

            for isim = 1:nsim
                _, n1, c1, t1, a1, as1 = stress!(utm, maxsum, true, NT, uavs, alpha, g)
                _, n2, c2, t2, a2, as2 = stress!(uncrd, maxsum, true, NT, uavs, alpha, g)
                _, n3, c3, t3, a3, as3 = stress!(coord, maxsum, true, NT, uavs, alpha, g)
                _, n4, c4, t4, a4, as4 = stress!(naive, maxsum, true, NT, uavs, alpha, g)
                _, n5, c5, t5, a5, as5 = stress!(distr, maxsum, true, NT, uavs, alpha, g)
                maxsumNLMS[:, inu] += [n1, n2, n3, n4, n5]
                maxsumNLMSBool[:, inu] += [sum(c1), sum(c2), sum(c3), sum(c4), sum(c5)]
                maxsumTimes[:, inu] += [t1, t2, t3, t4, t5]
                maxsumAlerts[:, inu] += [sum(a1), sum(a2), sum(a3), sum(a4), sum(a5)]
                maxsumAlertsCount[:, inu] += [as1, as2, as3, as4, as5]
            end # for isim

            maxminTimes /= nsim
            maxsumTimes /= nsim
            @printf("nuavs = %d: cputime = %.2e sec\n", nu, toq())
            # print("maxmin nlms: ", vec(maxminNLMS[:, inu]), "\n")
            print("maxmin nlms bool: ", vec(maxminNLMSBool[:, inu]), "\n")
            print("maxmin alerts: ", vec(maxminAlertsCount[:, inu]), "\n")
            # print("maxmin alert bool: ", vec(maxminAlerts[:, inu]), "\n")
            # print("average maxmin times: ", vec(maxminTimes[:, inu]), "\n")
            # print("maxsum nlms: ", vec(maxsumNLMS[:, inu]), "\n")
            print("maxsum nlms bool: ", vec(maxsumNLMSBool[:, inu]), "\n")
            print("maxsum alerts: ", vec(maxsumAlertsCount[:, inu]), "\n")
            # print("maxsum alert bool: ", vec(maxsumAlerts[:, inu]), "\n")
            # print("average maxsum times: ", vec(maxsumTimes[:, inu]), "\n")
        end # for inu

        outfile = string("../../data/results-", ibulk, ".jld")
        save(outfile, "maxminNLMS", maxminNLMS, "maxsumNLMS", maxsumNLMS,
             "maxminNLMSBool", maxminNLMSBool, "maxsumNLMSBool", maxsumNLMSBool,
             "maxminTimes", maxminTimes, "maxsumTimes", maxsumTimes,
             "maxminAlerts", maxminAlerts, "maxsumAlerts", maxsumAlerts,
             "maxminAlertsCount", maxminAlertsCount, "maxsumAlertsCount", maxsumAlertsCount)
        println("Results saved to ", outfile)
    end # for ibulk
end # function bulk_test

function bulk_test(nuavs::UnitRange{Int64}, alphafile::ASCIIString,
                   g::RectangleGrid, tag::ASCIIString,
                   nsim::Int64=1000, nbulk::Int64=10)
    alpha = read_alpha(alphafile)
    for ibulk = 1:nbulk
        println("beginning bulk stress test no. ", ibulk)
        maxminNLMS = zeros(Int64, 5, length(nuavs))
        maxsumNLMS = zeros(Int64, 5, length(nuavs))
        maxminNLMSBool = zeros(Int64, 5, length(nuavs))
        maxsumNLMSBool = zeros(Int64, 5, length(nuavs))
        maxminTimes = zeros(Float64, 5, length(nuavs))
        maxsumTimes = zeros(Float64, 5, length(nuavs))
        maxminAlerts = zeros(Int64, 5, length(nuavs))
        maxsumAlerts = zeros(Int64, 5, length(nuavs))
        maxminAlertsCount = zeros(Int64, 5, length(nuavs))
        maxsumAlertsCount = zeros(Int64, 5, length(nuavs))

        for inu = 1:length(nuavs)
            nu = nuavs[inu]
            uavs = randuavs(nu)
            tic()
            for isim = 1:nsim
                _, n1, c1, t1, a1, as1 = stress!(utm, maxmin, true, NT, uavs, alpha, g)
                _, n2, c2, t2, a2, as2 = stress!(uncrd, maxmin, true, NT, uavs, alpha, g)
                _, n3, c3, t3, a3, as3 = stress!(coord, maxmin, true, NT, uavs, alpha, g)
                _, n4, c4, t4, a4, as4 = stress!(naive, maxmin, true, NT, uavs, alpha, g)
                _, n5, c5, t5, a5, as5 = stress!(distr, maxmin, true, NT, uavs, alpha, g)
                maxminNLMS[:, inu] += [n1, n2, n3, n4, n5]
                maxminNLMSBool[:, inu] += [sum(c1), sum(c2), sum(c3), sum(c4), sum(c5)]
                maxminTimes[:, inu] += [t1, t2, t3, t4, t5]
                maxminAlerts[:, inu] += [sum(a1), sum(a2), sum(a3), sum(a4), sum(a5)]
                maxminAlertsCount[:, inu] += [as1, as2, as3, as4, as5]
            end # for isim

            for isim = 1:nsim
                _, n1, c1, t1, a1, as1 = stress!(utm, maxsum, true, NT, uavs, alpha, g)
                _, n2, c2, t2, a2, as2 = stress!(uncrd, maxsum, true, NT, uavs, alpha, g)
                _, n3, c3, t3, a3, as3 = stress!(coord, maxsum, true, NT, uavs, alpha, g)
                _, n4, c4, t4, a4, as4 = stress!(naive, maxsum, true, NT, uavs, alpha, g)
                _, n5, c5, t5, a5, as5 = stress!(distr, maxsum, true, NT, uavs, alpha, g)
                maxsumNLMS[:, inu] += [n1, n2, n3, n4, n5]
                maxsumNLMSBool[:, inu] += [sum(c1), sum(c2), sum(c3), sum(c4), sum(c5)]
                maxsumTimes[:, inu] += [t1, t2, t3, t4, t5]
                maxsumAlerts[:, inu] += [sum(a1), sum(a2), sum(a3), sum(a4), sum(a5)]
                maxsumAlertsCount[:, inu] += [as1, as2, as3, as4, as5]
            end # for isim

            maxminTimes /= nsim
            maxsumTimes /= nsim
            @printf("nuavs = %d: cputime = %.2e sec\n", nu, toq())
            # print("maxmin nlms: ", vec(maxminNLMS[:, inu]), "\n")
            print("maxmin nlms bool: ", vec(maxminNLMSBool[:, inu]), "\n")
            print("maxmin alerts: ", vec(maxminAlertsCount[:, inu]), "\n")
            # print("maxmin alert bool: ", vec(maxminAlerts[:, inu]), "\n")
            # print("average maxmin times: ", vec(maxminTimes[:, inu]), "\n")
            # print("maxsum nlms: ", vec(maxsumNLMS[:, inu]), "\n")
            print("maxsum nlms bool: ", vec(maxsumNLMSBool[:, inu]), "\n")
            print("maxsum alerts: ", vec(maxsumAlertsCount[:, inu]), "\n")
            # print("maxsum alert bool: ", vec(maxsumAlerts[:, inu]), "\n")
            # print("average maxsum times: ", vec(maxsumTimes[:, inu]), "\n")
        end # for inu

        outfile = string("../../data/results-", tag, "-", ibulk, ".jld")
        save(outfile, "maxminNLMS", maxminNLMS, "maxsumNLMS", maxsumNLMS,
             "maxminNLMSBool", maxminNLMSBool, "maxsumNLMSBool", maxsumNLMSBool,
             "maxminTimes", maxminTimes, "maxsumTimes", maxsumTimes,
             "maxminAlerts", maxminAlerts, "maxsumAlerts", maxsumAlerts,
             "maxminAlertsCount", maxminAlertsCount, "maxsumAlertsCount", maxsumAlertsCount)
        println("Results saved to ", outfile)
    end # for ibulk
end # function bulk_test

end # module Multiagent


text:
# Practical Machine learning
# Bayesian learning - Naive Bayes example
# Chapter 9

using NaiveBayes
using RDatasets

iris = dataset("datasets", "iris")

# observations in columns and variables in rows
X = array(iris[:, 1:4])'
p, n = size(X)
# by default species is a PooledDataArray,
y = [species for species in iris[:, 5]]

# how much data use for training
train_frac = 0.9
k = int(floor(train_frac * n))
idxs = randperm(n)
train = idxs[1:k]
test = idxs[k+1:end]

model = GaussianNB(unique(y), p)
fit(model, X[:, train], y[train])

accuracy = countnz(predict(model, X[:,test]) .== y[test]) / countnz(test)

println("Accuracy: $accuracy")


text:
"""
    AgedJacobianFactors

Type containing the Jacobian Factors and the age of the Jacobian.
This allows for the Shamanskii method to not update the Jacobian at each iterate.
"""
mutable struct AgedJacobianFactors
    fac        # Jacobian factors — can be Real, Complex, Dual, or HyperDual
    age::Int # age of the Jacobian
    facage::Int # age of the factors of the Jacobian
end

"""
    \\(Jf::JacobianFactors, y)

Dispatches backslash to work with all `JacobianFactors` subtypes.
"""
function Base.:\(Jf::AgedJacobianFactors, y)
    return Jf.fac \ y
end


text:
using HTTP

export isemail

function isemail(str::AbstractString, ignoreMaxLength::Bool=false, allowIpDomain::Bool=false, domainBlacklist::Vector{String}=Vector{String}())::Bool
    maxEmailLength = 254

    if !ignoreMaxLength && length(str) > maxEmailLength
        return false
    end

    parts = split(str, "@")
    domain = pop!(parts)

    if lowercase(domain) in domainBlacklist
        return false
    end

    user = join(parts, "@")

    if !ignoreMaxLength && (length(HTTP.escapeuri(user)) > 64 || length(HTTP.escapeuri(domain)) > 254)
        return false
    end

    if !isfqdn(domain, true)
        if !allowIpDomain
            return false
        end

        if !isip(domain)
            if domain[1] !== '[' || domain[end] !== ']'
                return false
            end

            noBracketsDomain = domain[2:end-1]
            if length(noBracketsDomain) === 0 || !isip(noBracketsDomain)
                return false
            end
        end
    end

    emailUserUtf8PartReg = r"^[a-z\d!#\$%&'\*\+\-\/=\?\^_`{\|}~\u00A0-\uD7FF\uF900-\uFDCF\uFDF0-\uFFEF]+$"i
    userParts = split(user, ".")

    if !all(map((part, ) -> contains(part, emailUserUtf8PartReg), userParts))
        return false
    end

    return true
end


text:
#=
# 7. Lid-driven cavity flow
In this notebook we will simulate the flow in a square cavity with a top moving wall.
=#
using ViscousFlow
#-
using Plots

#=
## Problem specification
Take $Re=100$ for example:
=#
my_params = Dict()
my_params["Re"] = 100

#=
## Discretization
Note that the rectangle function used for making the cavity shape requires a specified half length. The
immersed boundary projection method for internal flow requires the size of the domain to
be at least a step size greater at the boundaries (i.e. `halflength + Δx`).
So, for safety we make it 2 percent greater on each side.

Also, we set the grid Reynolds number to 0.5, lower than default of 2, so we have more
accuracy. We also set the CFL number to 1.0, larger than the default of 0.5, so we
can get to steady state faster.
=#
halflength=0.5
domain_lim=1.02*halflength
xlim, ylim = (-domain_lim,domain_lim),(-domain_lim,domain_lim)
my_params["grid Re"] = 0.5
my_params["CFL"] = 1.0
g = setup_grid(xlim,ylim,my_params)

Δs = surface_point_spacing(g,my_params)

#=
## Cavity Geometry
A square cavity can be created using the `Rectangle()` function with the half length defined above.
=#
body = Rectangle(halflength,halflength,Δs)
plot(body,fillrange=nothing)

#=
## Boundary Condition at the moving wall
Assign velocity to the top boundary (in the interior of the square). To do this,
we redefine the interior boundary condition function (the "interior" key of
the bc dictionary). We don't have to touch the exterior boundary condition,
because it defaults to zero, as we desire.

Rectangles and other polygon shapes in the `RigidBodyTools.jl`
package have a `side` field that allows easy access to the index range
for each side. Here, we use that feature to set the $x$ velocity
along the top to 1. We also use the `view` function to provide access
to the part of the overall velocity vector associated with body 1. That's
not particularly important in this example, since there is only one body,
but it is useful in problems that have multiple bodies.
=#
function my_vsminus(t,base_cache,phys_params,motions)
  vsminus = zeros_surface(base_cache)
  vsu = view(vsminus.u,base_cache,1)

  top = base_cache.bl[1].side[3]
  vsu[top] .= 1.0
  return vsminus
end

bcdict = Dict("interior" => my_vsminus)

#=
## Construct the system structure
Now we provide our parameters and the boundary condition dictionary.
=#
sys = viscousflow_system(g,body,phys_params=my_params,bc=bcdict);

#=
Initialize
=#
u0 = init_sol(sys)

#=
Set up integrator
=#
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)

#=
## Solve
For demonstration purposes, we run this over 2 convective time units. (You can run it longer to get it
to steady state.)
=#
step!(integrator,2)

#=
## Examine
Plot the vorticity and streamlines. We take a little care with the
streamlines to show the small recirculation zones in the lower corners.
=#
mins, maxs = extrema(streamfunction(integrator))
slevs = vcat(range(mins,maxs,length=11),maxs-0.00438,maxs-0.00353,maxs-0.00348)
plot(
plot(vorticity(integrator),sys,title="Vorticity",clim=(-10,10),color=:turbo,linewidth=1.5,ylim=ylim,levels=-6:0.5:5),
plot(streamfunction(integrator),sys,title="Streamfunction",color=:black,ylim=ylim,levels=slevs)
)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl,md
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

using NNlib,  Interact, ImageFiltering, Images, TestImages, ImageView, Plots, OffsetArrays

# +
x=-2:.2:2
y=-2:.2:2

σ_2 = 0.5
δ = σ_2/2
K = 12;
# -

G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
Plots.wireframe(x,y,G)

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
    Plots.wireframe(x,y,D)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
    R(x,y) = relu.(D(x,y))
#     L(x,y) = -relu.(-D(x,y))
    Plots.wireframe(x,y,R)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
#     R(x,y) = relu.(D(x,y))
    L(x,y) = -relu.(-D(x,y))
    Plots.wireframe(x,y,L)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
    R(x,y) = relu.(D(x,y))
    L(x,y) = -relu.(-D(x,y))
    S_a(x,y) = R(x,y) + L(x,y)
    Plots.wireframe(x,y,S_a)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
    R(x,y) = relu.(D(x,y))
    L(x,y) = -relu.(-D(x,y))
    S_b(x,y) = R(x,y) - L(x,y)
    Plots.wireframe(x,y,S_b)
end

@manipulate for k=1:1:2K, σ_2=0.2:0.05:5
    θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
#     d_sep_rp(x,y) = relu.(d(x,y))
#     d_sep_rm(x,y) = relu.(-d(x,y))
#     A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
#     B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,d)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
    D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
#     d_sep_rp(x,y) = relu.(d(x,y))
#     d_sep_rm(x,y) = relu.(-d(x,y))
#     A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
#     B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,D_sep)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    D(x,y) = G(x+(δ*cos(θ)), y + (δ*sin(θ))) - G(x - (δ*cos(θ)), y - (δ*sin(θ)))
    Plots.wireframe(x,y,D)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
    d_sep_rp(x,y) = relu.(d(x,y))
#     d_sep_rm(x,y) = relu.(-d(x,y))
#     A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
#     B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,d_sep_rp)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
#     d_sep_rp(x,y) = relu.(d(x,y))
    d_sep_rm(x,y) = relu.(-d(x,y))
#     A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
#     B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,d_sep_rm)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
    δ = σ_2/2
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
    d_sep_rp(x,y) = relu.(d(x,y))
    d_sep_rm(x,y) = relu.(-d(x,y))
    A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
#     B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,A)
end

@manipulate for k=1:1:2K
    θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
    d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
    d_sep_rp(x,y) = relu.(d(x,y))
    d_sep_rm(x,y) = relu.(-d(x,y))
    A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
    B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
    Plots.wireframe(x,y,B)
end

x=-10:1:10
y=-10:1:10
k=1
θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
d_sep_rp(x,y) = relu.(d(x,y))
d_sep_rm(x,y) = relu.(-d(x,y))
A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
Plots.wireframe(x,y,A)

x=-10:1:10
y=-10:1:10
k=1
θ = π*(k-1)/K
#     G(x,y) = 1/(2*π*σ_2^2)*exp(-1/(2*σ_2^2)*(x^2+y^2))
d(x,y) = exp(-(x*cos(θ)+y*sin(θ))/2σ_2)-exp((x*cos(θ)+y*sin(θ))/2σ_2)
#     D_sep(x,y) = G(x,y)*exp(-1/8)*d(x,y)
d_sep_rp(x,y) = relu.(d(x,y))
d_sep_rm(x,y) = relu.(-d(x,y))
A(x,y) = d_sep_rp(x,y) - d_sep_rm(x,y)
B(x,y) = d_sep_rp(x,y) + d_sep_rm(x,y)
Plots.wireframe(x,y,B)


text:
module Showing

using AWS
using AWS.EC2
using AWS.S3

import Base.show, Base.Meta.quot

print_indented(io::IO,s,indent) = print(io," "^indent,s)
show_indented(io::IO,s::ASCIIString,indent) = print(io,s*"\n")
show_indented(io::IO,::Void,indent) = print(io,"Not set\n")
show_indented(io::IO,i,indent) = println(io,i)

function show_indented(io::IO, s::Array, indent)
    println(io,eltype(s)," Array")
    for x in s
        print_indented(io,"",indent)
        show_indented(io,x,indent+4)
    end
end

function show_indented(io::IO, s::Dict, indent)
    println(io,eltype(s)," Dict")
    for (k,v) in collect(s)
        print_indented(io,"",indent)
        show_indented(io,string(k," => ",v),indent+4)
    end
end

macro show_func(n, t)
    #n,t,id_sym = eval(n),eval(t),eval(id_sym)
    block = Expr(:block)
    for x in fieldnames(t)
        push!(block.args,quote
            $(esc(:print_indented))(io,$(string(x)*": "),indent)
            $(esc(:show_indented))(io,r.($(quot(x))),indent+4)
        end)
    end
    ret = quote
        function $(esc(:show_indented))(io::IO,r::$(t),indent)
            $(esc(:print))(io,$(string(n)))
            println(io)
            $(block)
        end
        $(esc(:show))(io::IO,r::$(t)) = $(esc(:show_indented))(io,r,0)
    end
    ret
end

# Create default show method for all types
for m in [AWS, AWS.EC2, AWS.S3]
    for n in names(m)
        t = (m).(n)
        if isa(t,Type)
            @eval @show_func $n $t
        end
    end
end

end


text:
function preprocess(::DatasetName) end
extension(dn::DatasetName) = extension(dn)
extension(::Tabular) = "csv"

"""
    function call(dataset)

Download the dataset using `DataDeps` package
and create a csv file, see [`preprocess?`](@ref)).
"""
call(dataset::DatasetName) = @datadep_str name(dataset)

"""
    preprocess(path, dataset, header, kwargs...)

Create csv file containing data from tabular dataset in "standard" format.

The format can be described as - columns represents attributes, rows instances,
attributes in a row are separated by comma. First row of file is header.
Last column contains target values and is named Target. If `dataset` has defined
`categorical` columns, appends "-C" at the beginning of column name.

#Keyword arguments
- `kwargs...`: keyword arguments that are possible in `CSV.File` function.
"""
function preprocess(
    path::AbstractString,
    dataset::Tabular;
    kwargs...
)
    # If header is Integer, it means it is in the file together with data and is used in
    # CSV.file loading.
    headers(dataset) isa Integer ? header = headers(dataset) : header = false

    df = CSV.File(
        path,
        header = header,
        transpose = transposed(dataset),
        missingstrings = ["", "NA", "?", "*", "#DIV/0!", "missing", "NaN"],
        truestrings = ["T", "t", "TRUE", "true", "y", "yes", "Y"],
        falsestrings = ["F", "f", "FALSE", "false", "n", "no", "N"];
        kwargs...
        ) |> DataFrame

    # if header is an Int or true or if header is defined in $dataset.jl
    # saves header to a file header.csv
    # does not overwrite
    if headers(dataset) isa Int || header
        save_header(names(df), path)
    elseif !isempty(headers(dataset))
        save_header(headers(dataset), path)
    end

    # place target column
    col = target(dataset)    #target column, either Int or Type{File}

    typeSplit = find_in(basename(path))

    df = place_target(col, df, typeSplit)

    categorical_cols = categorical(dataset)
    for i in categorical_cols
        rename!(df, i => names(df)[i]*"-C")
    end

    ext = extension(dataset)

    pathForSave = joinpath(dirname(path), string("data-",typeSplit,".",ext))
    CSV.write(pathForSave, df, delim=',', writeheader=true)
    rm(path)
end

function save_header(names::Vector{<:AbstractString}, path::AbstractString)
    open(joinpath(dirname(path), "header.csv"),"w") do io
        [write(io, d*"\n") for d in names]
    end
end

function rename_cols(df::AbstractDataFrame)
    return rename!(df, collect(1:(ncol(df)-1)) .=> ["Column $i" for i in 1:ncol(df)-1])
end

"""
    place_target(column::Int, df::AbstractDataFrame)

Moves column from its position to last position in `DataFrame` df.
"""
function place_target(column::Int, df::AbstractDataFrame, ::AbstractString)
    lastColIndex = ncol(df)

    if column > 0 || column < lastColIndex
        #Move target column to last position
        df.Target = df[!, column]
        df = df[!, Not(column)]
        # rename columns after target column was moved, don't if col was already the last column
        rename_cols(df)
    else
        rename!(df, lastColIndex => "Target")
    end

    return df
end

"""
    place_target(::Type{Labels}, df::AbstractDataFrame, typeSplit::AbstractString)

Pushes column with labels from a file to last position in `DataFrame` df.
"""
function place_target(::Type{Labels}, df::AbstractDataFrame, typeSplit::AbstractString)
    pathLabels = joinpath(pwd(), string("labels-", typeSplit))
    dfLabels = CSV.File(
        pathLabels,
        header = ["Target"],
        missingstrings = ["", "NA", "?", "*", "#DIV/0!"],
        truestrings = ["T", "t", "TRUE", "true", "y", "yes"],
        falsestrings = ["F", "f", "FALSE", "false", "n", "no"],
        ) |> DataFrame
    df = hcat(df, dfLabels)
    rm(pathLabels)

    return df
end

"""
    preprocess(path::AbstractString, ::Type{Labels})

Renames downloaded file with labels based on type. Filename has format labels-typeSplit.csv.
For typeSplit see [`find_in`](@ref).
"""
function preprocess(path::AbstractString, ::Type{Labels})
    typeSplit = find_in(getfilename(path))
    mv(basename(path), string("labels-", typeSplit))
end

"""
    preprocess(path::AbstractString, ::Type{Header})

Renames downloaded file with header to header.csv.
"""
function preprocess(path::AbstractString, ::Type{Header})
    mv(basename(path), string("header.csv"))
end
"""
    find_in(name::AbstractString)

Searches for strings train, test, valid in given name.

If one of the substrings is present then returns it.
If none is present then returns "train".

#Examples
```julia-repl
julia> PreprocessData.find_in("name-test.csv")
"test"

julia> PreprocessData.find_in("name.csv")
"train"
```
"""
function find_in(name::AbstractString)
    if occursin("train", name) || (!occursin("valid", name) && !occursin("test", name))
        return "train"
    elseif occursin("valid", name)
        return "valid"
    elseif occursin("test", name)
        return "test"
    end
end

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
module Runtests

import ModiaResult
using  Test

@testset "Test ModiaPlot_CairoMakie/test" begin
    ModiaResult.usePlotPackage("CairoMakie")
    include("$(ModiaResult.path)/test/runtests_withPlot.jl")
    include("test_45_Save_svg_figure.jl")
    ModiaResult.usePreviousPlotPackage()
end

end


text:
#=
common.jl

Adham Beyki, odinay@gmail.com
=#

abstract MixtureModel
abstract hMixtureModel

immutable Sent
    w::Vector{Int}
    c::Vector{Int}
    Sent(w::Vector{Int}, c::Vector{Int}) = new(w, c)
end
Sent(n::Int) = Sent(zeros(Int, n), zeros(Int, n))
Base.length(x::Sent) = length(x.w)

"""
normalizes a vector with logarithmic scale values
"""
function lognormalize!(pp::Vector{Float64})
    pp_len = length(pp)
    pp_max = maximum(pp)
    @inbounds for kk = 1:pp_len
        pp[kk] = exp(pp[kk] - pp_max)
    end
    pp_sum = sum(pp)
    @inbounds for kk = 1:pp_len
        pp[kk] /= pp_sum
    end
end

"""
draws samples from a probability vector
"""
function sample(w::Vector{Float64})
    r = rand()
    n = length(w)
    i = 1
    cw = w[1]
    while cw < r && i < n
        i += 1
        @inbounds cw += w[i]
    end
    return i
end
function sample(w::Vector{Float64}, n::Int)
    ret = zeros(Int, n)
    @inbounds for i = 1:n
        ret[i] = sample(w)
    end
    return ret
end

rand_Dirichlet(alpha::Vector{Float64}) = Distributions.rand(Distributions.Dirichlet(alpha))
rand_Dirichlet(alpha::Vector{Float64}, n) = Distributions.rand(Distributions.Dirichlet(alpha), n)

"""
deletes column kk from matrix nn and returns nn
"""
function del_column{T}(nn::Matrix{T}, kk::Int)
  KK = size(nn, 2)
  mask = 1:KK .!= kk
  return nn[:, mask]
end

"""
adds a column to the rightmost of nn and returns nn
"""
function add_column{T}(nn::Matrix{T})
  r, c = size(nn)
  mm = zeros(T, r, c+1)
  mm[:, 1:c] = nn
  return mm
end

"""
generates bar topics
"""
function gen_bars(n_bars, n_vocab, noise_level)
    KK = round(Int, sqrt(n_vocab))
    bars = zeros(Float64, n_bars, n_vocab)

    for kk in 1:round(Int, n_bars/2)
        b = zeros(KK, KK) + noise_level
        b[kk, :] = ones(KK)
        b /= sum(b)
        bars[kk, :] = b[:]
    end

    for kk in round(Int, n_bars/2)+1 : n_bars
        b = zeros(KK, KK) + noise_level
        b[:, kk-round(Int, n_bars/2)] = ones(KK)
        b /= sum(b)
        bars[kk, :] = b[:]
    end
    bars
end

"""
vertical cumulative sum in reverse order from bottom to top
"""
function reverse_cumsum{T}(nn::Matrix{T})

  J = size(nn, 1)
  mm = zeros(T, J, J)
  mm[J, :] = nn[J, :]
  for jj = 2:J
    mm[J-jj+1, :] = nn[J-jj+1, :] + mm[J-jj+2, :]
  end
  mm
end

"""
writes a topic into a CSV file
"""
function topic2csv(filename, vocab, alpha)

  idx_sorted = sortperm(alpha, rev=true)

  words = Array(Tuple{ASCIIString, Float64}, length(vocab))
  for i = 1:length(idx_sorted)
    words[i] = (vocab[idx_sorted[i]], alpha[idx_sorted[i]])
  end

  csvfile = open(filename, "w")

  for i = 1:length(vocab)
    write(csvfile, join(words[i], ","), "\n")
  end

  close(csvfile)
end

"""
writes top n topics used in a document
"""
function write_top_doctopics(pij, filename, topn=5)
  csvfile = open(filename, "w")

  for kk = 1:size(pij, 2)
    idx = sortperm(pij[:, kk], rev=true)[1:topn]
    insert!(idx, 1, kk)
    write(csvfile, join(idx, ","), "\n")
  end

  close(csvfile)
end

# macro for argument checking
macro check_args(D, cond)
    quote
        if !($cond)
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)), " is not satisfied.")))
        end
    end
end

function stick_breaking(vv::Vector{Float64})
    KK = length(vv)
    pp = ones(Float64, KK)

    pp[2:KK] = 1-vv[1:KK-1]
    pp = vv .* cumprod(pp)
    pp
end

function log_stick_breaking(vv::Vector{Float64})
  KK = length(vv)
  pp = zeros(Float64, KK)

  pp[2:KK] = log(1 - vv[1:KK-1])
  pp = log(vv) + cumsum(pp)
  pp
end

"""
vertical cumulative sum in reverse order
from bottom to top
"""
function reverse_cumsum{T}(nn::Matrix{T})

  J = size(nn, 1)
  mm = zeros(T, J, J)
  mm[J, :] = nn[J, :]
  for jj = 2:J
    mm[J-jj+1, :] = nn[J-jj+1, :] + mm[J-jj+2, :]
  end
  mm
end


text:
"""
    BroadcastTensor{T,N} <: AbstractArray{T,N}

An efficient broadcasting wrapper for `Tensor`.
Sub-types `AbstractArray` to make use of existing broadcasting utilities.
"""
struct BroadcastTensor{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
    BroadcastTensor(t::Tensor{T,N}) where {T,N} = new{T,N}(t.data)
end

# behaves like getindex(::Tensor, ...) but returns scalars instead of scalar tensors
function Base.getindex(t::BroadcastTensor, inds...)
    val = getindex(t.data, inds...)
    return typeof(val) <: AbstractArray ? Tensor(val) : val
end

# needed to allow broadcasting a BroadcastTensor
Base.axes(t::BroadcastTensor, inds...) = axes(t.data, inds...)

# define the broadcasting behaviour of tensors using BroadcastTensor
Base.broadcastable(t::Tensor) = BroadcastTensor(t)
Base.BroadcastStyle(::Type{<:Tensor}) = Broadcast.ArrayStyle{BroadcastTensor}()
Base.BroadcastStyle(::Type{<:BroadcastTensor}) = Broadcast.ArrayStyle{BroadcastTensor}()

function Base.similar(
    bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{BroadcastTensor}}, ::Type{ElType}
) where {ElType}
    return Tensor(similar(Array{ElType}, axes(bc)))
end

Base.copyto!(t::Tensor, bc::Broadcast.Broadcasted) = Tensor(copyto!(t.data, bc))


text:
mutable struct GFAssessState
    trace::GFTrace
    constraints::Any
    score::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFAssessState(constraints, params::Dict{Symbol,Any})
    GFAssessState(GFTrace(), constraints, 0., AddressVisitor(), params)
end

get_args_change(state::GFAssessState) = nothing
get_addr_change(state::GFAssessState, addr) = nothing
set_ret_change!(state::GFAssessState, value) = begin end

function addr(state::GFAssessState, dist::Distribution{T}, args, addr, delta) where {T}
    visit!(state.visitor, addr)
    retval::T = get_leaf_node(state.constraints, addr)
    score = logpdf(dist, retval, args...)
    call = CallRecord(score, retval, args)
    state.trace = assoc_primitive_call(state.trace, addr, call)
    state.trace.has_choices = true
    state.score += score
    retval
end

function addr(state::GFAssessState, gen::Generator{T,U}, args, addr, delta) where {T,U}
    visit!(state.visitor, addr)
    if has_internal_node(state.constraints, addr)
        constraints = get_internal_node(state.constraints, addr)
    else
        constraints = EmptyAssignment()
    end
    trace::U = assess(gen, args, constraints)
    call::CallRecord = get_call_record(trace)
    state.trace = assoc_subtrace(state.trace, addr, trace)
    state.trace.has_choices |= has_choices(trace)
    state.score += call.score
    call.retval::T
end

splice(state::GFAssessState, gf::GenFunction, args::Tuple) = exec(gf, state, args)

function assess(gen::GenFunction, args, constraints)
    state = Gen.GFAssessState(constraints, gen.params)
    retval = Gen.exec(gen, state, args)
    unconsumed = get_unvisited(state.visitor, constraints)
    if !isempty(unconsumed)
        error("Update did not consume all constraints: $unconsumed")
    end

    # TODO add return type annotation for gen
    call = Gen.CallRecord{Any}(state.score, retval, args)
    state.trace.call = call
    state.trace
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
function findunique(collection, item)
    for (i,x) in pairs(IndexLinear(), collection)
        if x == item
            return i
        end
    end
    error("$item not found")
end

function generate_collection_table(sets, rarity_by_set, rarity)
    Dict(map(i -> (i, zeros(Int, rarity_by_set[i][rarity])), sets))
end

function generate_collection_table(sets, rarity_by_set)
    (generate_collection_table(sets, rarity_by_set, 1),
    generate_collection_table(sets, rarity_by_set, 2),
    generate_collection_table(sets, rarity_by_set, 3),
    generate_collection_table(sets, rarity_by_set, 4))
end

function copy_to_dict(keys, x)

    Dict(map(i -> (i, deepcopy(x)), keys))

end

function make_counter(sets, dims, T = Int)
    copy_to_dict(sets, zeros(T, dims))
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
function infer_type(binding::Binding, scope, state)
    if binding isa Binding
        binding.type !== nothing && return
        if binding.val isa EXPR && CSTParser.defines_module(binding.val)
            binding.type = CoreTypes.Module
        elseif binding.val isa EXPR && typof(binding.val) === CSTParser.FunctionDef
            binding.type = CoreTypes.Function
        elseif binding.val isa EXPR && CSTParser.defines_datatype(binding.val)
            binding.type = CoreTypes.DataType
        elseif binding.val isa EXPR && typof(binding.val) === BinaryOpCall
            if kindof(binding.val[2]) === CSTParser.Tokens.EQ
                if CSTParser.is_func_call(binding.val[1])
                    binding.type = CoreTypes.Function
                elseif CSTParser.is_func_call(binding.val[3])
                    callname = CSTParser.get_name(binding.val[3])
                    if CSTParser.isidentifier(callname)
                        resolve_ref(callname, scope, state, [])
                        if hasref(callname)
                            rb = get_root_method(refof(callname), state.server)
                            if (rb isa Binding && (rb.type == CoreTypes.DataType || rb.val isa SymbolServer.DataTypeStore)) || rb isa SymbolServer.DataTypeStore
                                binding.type = rb
                            end
                        end
                    end
                elseif kindof(binding.val[3]) === CSTParser.Tokens.INTEGER
                    binding.type = CoreTypes.Int
                elseif kindof(binding.val[3]) === CSTParser.Tokens.FLOAT
                    binding.type = CoreTypes.Float64
                elseif kindof(binding.val[3]) === CSTParser.Tokens.STRING || typof(binding.val[3]) === CSTParser.Tokens.TRIPLE_STRING
                    binding.type = CoreTypes.String
                elseif typof(binding.val[3]) === IDENTIFIER && refof(binding.val[3]) isa Binding
                    binding.type = refof(binding.val[3]).type
                end
            elseif kindof(binding.val[2]) === CSTParser.Tokens.DECLARATION
                t = binding.val[3]
                if CSTParser.isidentifier(t)
                    resolve_ref(t, scope, state, [])
                end
                if typof(t) === CSTParser.Curly
                    t = t[1]
                    resolve_ref(t, scope, state, [])
                end
                if is_getfield(t) && length(t[3]) > 0
                    t = t[3][1]
                    resolve_ref(t, scope, state, [])
                end

                if refof(t) isa Binding
                    rb = get_root_method(refof(t), state.server)
                    if rb isa Binding && rb.type == CoreTypes.DataType
                        binding.type = rb
                    else
                        binding.type = refof(t)
                    end
                elseif refof(t) isa SymbolServer.DataTypeStore
                    binding.type = refof(t)
                end
            end
        end
    end
end


text:
using ImageFiltering, Base.Test

@testset "extrema_filter" begin
    function groundtruth(f, A, window::Tuple)
        Aex = copy(A)
        hshift = map(x->x>>1+1, window)
        for Ishift in CartesianRange(window)
            for I in CartesianRange(size(A))
                Aex[I] = f(Aex[I], A[map(d->clamp(I[d]+Ishift[d]-hshift[d], 1, size(A,d)), 1:ndims(A))...])
            end
        end
        Aex
    end
    groundtruth(f, A, window) = groundtruth(f, A, (window...,))
    # 1d case
    A = [0.1,0.3,0.2,0.3,0.4]
    mm = extrema_filter(A, 1)
    @test mm == [(a,a) for a in A]
    mm = extrema_filter(A, 2)
    @test mm == [(0.1,0.1),(0.1,0.3),(0.2,0.3),(0.2,0.3),(0.3,0.4)]
    mm = extrema_filter(A, 3)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.2,0.3),(0.2,0.4),(0.3,0.4)]
    mm = extrema_filter(A, 4)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.3),(0.2,0.4),(0.2,0.4)]
    mm = extrema_filter(A, 5)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.4),(0.2,0.4),(0.2,0.4)]
    mm = extrema_filter(A, 6)
    @test mm == [(0.1,0.3),(0.1,0.3),(0.1,0.4),(0.1,0.4),(0.2,0.4)]
    mm = extrema_filter(A, 7)
    @test mm == [(0.1,0.3),(0.1,0.4),(0.1,0.4),(0.1,0.4),(0.2,0.4)]
    A = [0.1,0.3,0.5,0.4,0.2]
    mm = extrema_filter(A, 1)
    @test mm == [(a,a) for a in A]
    mm = extrema_filter(A, 2)
    @test mm == [(0.1,0.1),(0.1,0.3),(0.3,0.5),(0.4,0.5),(0.2,0.4)]
    mm = extrema_filter(A, 3)
    @test mm == [(0.1,0.3),(0.1,0.5),(0.3,0.5),(0.2,0.5),(0.2,0.4)]
    mm = extrema_filter(A, 4)
    @test mm == [(0.1,0.3),(0.1,0.5),(0.1,0.5),(0.2,0.5),(0.2,0.5)]
    mm = extrema_filter(A, 5)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5),(0.2,0.5)]
    mm = extrema_filter(A, 6)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5)]
    mm = extrema_filter(A, 7)
    @test mm == [(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.1,0.5),(0.2,0.5)]
    # 2d case
    A = rand(5,5)/10
    A[2,2] = 0.8
    A[4,4] = 0.6
    for w in ((2,2), (2,3), (3,2), (3,3), (2,5))
        mm = extrema_filter(A, w)
        maxval = map(last, mm)
        Amax = groundtruth(max, A, w)
        @test maxval == Amax
        minval = map(first, mm)
        Amin = groundtruth(min, A, w)
        @test minval == Amin
    end
    # 3d case
    A = rand(5,5,5)/10
    A[2,2,2] = 0.7
    A[4,4,2] = 0.4
    A[2,2,4] = 0.5
    for w in ((2,2,2), (2,3,2), (3,2,2), (2,2,3), (3,3,3), (2,5,3))
        mm = extrema_filter(A, w)
        maxval = map(last, mm)
        Amax = groundtruth(max, A, w)
        @test maxval == Amax
        minval = map(first, mm)
        Amin = groundtruth(min, A, w)
        @test minval == Amin
    end
end

nothing


text:
using GeneralizedSDistributions
using Distributions

D = NegativeBinomial(20, 0.2)
G = GSDist(D)


text:
include("SphericalTensor.jl")

using LinearAlgebra
using SpecialFunctions

export helmholtzG, helmholtzGapprox, hologram

function sphericalbesselj(nu, x::T) where {T}
    besselj_nuhalf_x = besselj(nu + one(nu)/2, x)
    if abs(x) ≤ sqrt(eps(real(zero(besselj_nuhalf_x))))
        nu == 0 ? one(besselj_nuhalf_x) : zero(besselj_nuhalf_x)
    else
        √((float(T))(π)/2x) * besselj_nuhalf_x
    end
end

sphericalbessely(nu, x::T) where {T} = √((float(T))(π)/2x) * bessely(nu + one(nu)/2, x)

sphericalhankel1(nu, x::T) where {T} = sphericalbesselj(nu, x) + im*sphericalbessely(nu, x)

function helmholtzG(x,y;k=1)
    d = norm(x .- y)
    exp(im*k*d)/(4π*d)
end

function sphericalbesselj(nu, x::T) where {T}
    besselj_nuhalf_x = besselj(nu + one(nu)/2, x)
    if abs(x) ≤ sqrt(eps(real(zero(besselj_nuhalf_x))))
        nu == 0 ? one(besselj_nuhalf_x) : zero(besselj_nuhalf_x)
    else
        √((float(T))(π)/2x) * besselj_nuhalf_x
    end
end

function helmR(l,x;k=1) # inside solution
    conj(sphericalbesselj(l, k*norm(x)))*
        harmonic(l, x)
end

function helmS(l,x;k=1) # outside solution
    sphericalhankel1(l, k*norm(x))*
        harmonic(l, x)
end

function hologram(l,x;k=1,R=1) # projection on sphere
    # outputs a spherical tensor a, order l
    # can reconstruct it's portion of the image with harmonic(a,x)
    if norm(x) < R
        im*k*conj(helmR(l,x;k=k))*sphericalhankel1(l,k*R)
    else
        im*k*conj(sphericalbesselj(l, k*R)*helmS(l,x;k=k))
    end
end

function helmholtzGapprox(x,y;k=1,lmax=10)
    if norm(x) < norm(y)
        sum(im*k*dot(helmR(l,x;k=k),helmS(l,y;k=k)) for l in 0:lmax)
    else
        sum(im*k*dot(helmR(l,y;k=k),helmS(l,x;k=k)) for l in 0:lmax)
    end
end


text:
using SpatQSims
using Test

@testset "SpatQSims.jl" begin
    base_dir = mktempdir()
    cd(base_dir) do
        prep_sims(1, "prep.h5")
        @test "prep.h5" in readdir()

        include("qdevscaling_test.jl")
        include("sharedq_test.jl")
        include("prefintensity_test.jl")
        include("densdepq_test.jl")
        include("habq_test.jl")
    end
end