text:
####
# Shape
####
# TODO: Shape inference
typealias Shape Vector{Int}

####
# Operations
####

# TODO: Does this type hierarchy make any sense?  Think carefully about what's necessary
# Also whether it needs a hierarchy at all.  Unclear if we make use of it anywhere
# Additionally, can these be defined together with other parts of command?

abstract VarOp <: OpType
abstract CreateOp <: OpType
abstract ConstantOp <: CreateOp
abstract RandomOp <: CreateOp  # TODO: Make these!
abstract ElementWise <: OpType
abstract Activations <: ElementWise

# TODO: How to do control edges / grouping
immutable Print <: OpType end
immutable Noop <: OpType end

immutable Fill <: ConstantOp end
immutable FillLike <: ConstantOp end

# Values provided as input
immutable Placeholder <: VarOp
    shape::Shape
end

# Some value that is initialized once (on first run) and shared across runs
immutable Variable <: VarOp
    init::Node
end

immutable Constant <: ConstantOp
    value::TensorValue
end

# Specifies a mutable value
function variable(init::Node, name::AbstractString="")
    # TODO: must have shape specified / be able to infer
    Node(Variable(init), name)
end

# Specifies an input variable
function placeholder(shape::Shape, name::AbstractString="")
    Node(Placeholder(shape), name)
end

function constant(value::TensorValue, name::AbstractString="")
    Node(Constant(value), name)
end

type RandN <: RandomOp end
type Copy <: ElementWise end

type MatMul <: OpType end  # Matrix multiply
type Transpose <: OpType end
type Sum <: OpType end
type Dim <: OpType end
type Assign <: OpType end
type InPlaceAdd <: OpType end

type Mean <: OpType end
type Sum <: OpType end
type Maximum <: OpType end

# These are mainly for testing
type GetIndex <: OpType end
type SetIndex <: OpType end

#############################
# Operation creation macros #
#############################

# Each operation needs:
## 1. Type (identifier)
## 2. functions (call or infix, e.g. add and +)
## 4. Implementation (CPU/GPU)
## 3. gradients
## 4. shape inference

function gen_args(narg, typ::Type)
    assert(narg >= 0)
    args = []
    apply_args = []
    # Var names 'a'...'z'
    for var = 'a':('a' + narg - 1)
        varsym = symbol(var)
        push!(args, :($(varsym)::$(typ)))
        push!(apply_args, varsym)
    end
    args, apply_args
end

"""
@register_op  Mul (.*) 2
Expands to
function .*(a::Node, b::Node) Node(Mul(), [a, b]) end
"""
macro register_op(typ, op, narg)
    # TODO: Is there a way to interpolate an expr (like splat) into another expr with $ or similar?
    # For now, use Expr function (for which we can use splat).
    # Think it's pretty clear for now
    args, apply_args = gen_args(narg, Node)
    Expr(:function,
         Expr(:call,
              esc(op),
              args...,
              Expr(:kw,
                   Expr(:(::),
                        :name,
                        :AbstractString),
                   ""
                   )
              ),
         Expr(:call,
              :Node,
              Expr(:call, typ),
              Expr(:vect, apply_args...),
              :name))
end

"""
@register_grad Mul (b * ds) (a * ds)
 Expands to
function grad(op::Mul, out::Node, ds::Node, a::Node, b::Node)
    [b * ds, a * ds]
end
TODO: Make this take inputs[], outputs[], gradients[] explicitly as lists
"""
macro register_grad(typ, grads...)
    args, _ = gen_args(length(grads), Node)
    Expr(:function,
         Expr(:call,
              esc(:grad),
              :(op::$typ),
              :(out::Node),
              :(ds::Node),
              args...),
         Expr(:vect,
              grads...))
end

"""
@register_impl Mul 3 (a + b + c)
Expands to
function call(op::Mul, a::TensorValue, b::TensorValue, c::TensorValue)
    a + b + c
end
"""
macro register_impl(typ, narg, impl)
    args, _ = gen_args(narg, TensorValue)
    Expr(:function,
         Expr(:call,
              esc(:call),
              :(op::$typ),
              args...),
         impl)
end

"""
@register_shape Mul 3 (a + b)
Expands to
function shape(op::Mul, a::Shape, b::Shape)
    # Check compatability
    ...
    # Return new result
    a + b
end
"""
macro register_shape(typ, narg, impl)
    args, _ = gen_args(narg, Shape)
    Expr(:function,
         Expr(:call,
              esc(:shape),
              :(op::$typ),
              args...),
         impl)
end

#########################
# Operation Definitions #
#########################
# TODO: Define all parts of an operation together or keep similar parts grouped?
# e.g. @createOp Add, add, [a, b], [g, g], a + b, [shape inference]
#                type, name, args, gradients, implementation, [shape inference]
# Todo ops:
    # boolean operators
    # random

# noop basics
function noop(a::Node...)
    Node(Noop(), collect(a))
end

# Just have it return a scalar until theres some notion of ordering
function call(op::Noop, a::TensorValue...)
    1
end

# Wrapper on fill
fill(val::Float, shape::Array{Int}, name::AbstractString="") = fill(constant(val), constant(shape))
@register_op Print       print        1
@register_op Dim         dim          1
@register_op Fill        fill         2
zeros(a::Node) = fill(constant(0.0), a)
ones(a::Node) = fill(constant(1.0), a)
@register_op FillLike    fill_like    2
@register_op Copy        copy         1

@register_op MatMul      matmul       2
@register_op Transpose   t            1
@register_op Transpose   transpose    1
@register_op Assign      (.=)         2
@register_op InPlaceAdd  plusequals   2  # += doesn't work

@register_op Sum         sum          1
@register_op Sum         sum          2 # Arg 2 = axis
@register_op Mean        mean         1
@register_op Mean        mean         2

#@register_op Relu        relu         1

@register_op RandN       randn        1

@register_op Maximum     maximum      1
@register_op Maximum     maximum      2

@register_op GetIndex    getindex     2
@register_op SetIndex    setindex     3

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

@register_impl Print        1   (println(a); a)
@register_impl Constant     0   op.value
# a = scalar, b = 1d array specifying dims
@register_impl Fill         2   Base.fill(a, b...)
# a = scalar, b = array shape to copy
@register_impl FillLike     2   (dim = size(b); length(dim) == 0 ? a : fill(a, dim...))

@register_impl Dim          1   collect(Int, size(a))

@register_impl Copy         1   a

@register_impl MatMul       2   (a * b)
@register_impl Transpose    1   transpose(a)
# TODO: Run without returning anything.  Explicit NoValue
@register_impl InPlaceAdd   2   (for i in 1:length(a); a[i] += b[i] end; 0)

@register_impl Sum          1   Base.sum(a)
@register_impl Sum          2   Base.sum(a, b)
@register_impl Mean         1   Base.mean(a)
@register_impl Mean         2   Base.mean(a, b)

@register_impl RandN         1   Base.randn(a...)

@register_impl Maximum       1   maximum(a)
@register_impl Maximum       2   maximum(a, b)

@register_impl GetIndex      2   (a[b])
@register_impl SetIndex      3   (t = copy(a); t[b] = c; t)

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

@register_grad Print ds
@register_grad MatMul matmul(ds, t(b)) matmul(t(a), ds)
@register_grad Transpose t(ds)

@register_grad Sum repeatto(ds, dim(a))
@register_grad Sum repeatto(ds, dim(a)) (cg.constant(0)) # TODO: make axis nondiff

# TODO: I believe this is wrong for some edge cases
@register_grad Maximum  (eq(a, out) * ds)
@register_grad Maximum  broadcast("*", broadcast("==", a, out), ds) (b)  # 2nd one should be undefined

# Incredibly inefficient, but mostly for testing/debugging
# TODO: 0s Replace with nondiff
@register_grad GetIndex  (t = fill(cg.constant(0.0), dim(a)); setindex(t, b, ds)) (cg.constant(0))
@register_grad SetIndex  (setindex(ds, b, cg.constant(0.0))) (cg.constant(0)) (getindex(ds, b))

# TODO How to treate OnesLike etc. in gradient computations?
# TODO: Add actual GradUndefined
## TODO: May want to start grouping together like this

###########
# Scalars #
###########
abstract ScalarOp <: OpType

type Add <: ScalarOp end
type Sub <: ScalarOp end
type Mul <: ScalarOp end
type Div <: ScalarOp end
type Pow <: ScalarOp end

type Neg <: ScalarOp end
type Sign <: ScalarOp end
type Exp <: ScalarOp end
type Log <: ScalarOp end
type Sin <: ScalarOp end
type Cos <: ScalarOp end
type Abs <: ScalarOp end

type Max <: ScalarOp end
type Min <: ScalarOp end

type Eq <: ScalarOp end
type Neq <: ScalarOp end
type Le <: ScalarOp end
type Leq <: ScalarOp end
type Ge <: ScalarOp end
type Geq <: ScalarOp end

type Sigmoid <: ScalarOp end

@register_op Add     (+)          2
@register_op Sub     (-)          2
@register_op Mul     (*)          2
@register_op Div     (/)          2
@register_op Pow     (^)          2

@register_op Neg         (-)          1
@register_op Sign        sign         1
@register_op Exp         exp          1
@register_op Log         log          1
@register_op Sin         sin          1
@register_op Cos         cos          1
@register_op Abs         abs          1

@register_op Max     max          2
@register_op Min     min          2

# Would like to use == etc., but spell out for now
@register_op Eq      (eq)         2
@register_op Neq     (neq)        2
@register_op Le      (le)         2
@register_op Leq     (leq)        2
@register_op Ge      (ge)         2
@register_op Geq     (geq)        2

@register_op Sigmoid     sigmoid      1

# Basically reimplementing Base here
function elwise(op::Function, a::Real, b::Real)
    op(a, b)
end

# TODO: these are similar to broadcast in Base, but many times slower.  Why.
# TODO Similar isn't really right - should have proper type
function elwise(op::Function, a::Real, b::Array)
    out = similar(b)
    for i = eachindex(b)
        @inbounds out[i] = op(a, b[i])
    end
    reshape(out, size(b))
end
function elwise(op::Function, a::Array, b::Real)
    out = similar(a)
    for i = eachindex(a)
        @inbounds out[i] = op(a[i], b)
    end
    reshape(out, size(a))
end

function elwise(op::Function, a::Array, b::Array)
    @assert size(a) == size(b)
    out = similar(a)
    for i = eachindex(a)
        @inbounds out[i] = op(a[i], b[i])
    end
    reshape(out, size(a))
end

call(op::Add,  a::TensorValue, b::TensorValue)  = elwise(+, a, b)
call(op::Sub,  a::TensorValue, b::TensorValue)  = elwise(-, a, b)
call(op::Mul,  a::TensorValue, b::TensorValue)  = elwise(*, a, b)
call(op::Div,  a::TensorValue, b::TensorValue)  = elwise(/, a, b)
call(op::Pow,  a::TensorValue, b::TensorValue)  = elwise(^, a, b)

# These are autobroadcast
call(op::Neg,  a::TensorValue)           = -a
call(op::Sign, a::TensorValue)           = sign(a)
call(op::Exp,  a::TensorValue)           = exp(a)
call(op::Log,  a::TensorValue)           = log(a)
call(op::Sin,  a::TensorValue)           = sin(a)
call(op::Cos,  a::TensorValue)           = cos(a)
call(op::Abs,  a::TensorValue)           = abs(a)

call(op::Max,  a::TensorValue, b::TensorValue)  = elwise(max, a, b)
call(op::Min,  a::TensorValue, b::TensorValue)  = elwise(min, a, b)

call(op::Eq,   a::TensorValue, b::TensorValue)  = elwise(==, a, b)
call(op::Neq,  a::TensorValue, b::TensorValue)  = elwise(!=, a, b)
call(op::Le,   a::TensorValue, b::TensorValue)  = elwise(<, a, b)
call(op::Leq,  a::TensorValue, b::TensorValue)  = elwise(<=, a, b)
call(op::Ge,   a::TensorValue, b::TensorValue)  = elwise(>, a, b)
call(op::Geq,  a::TensorValue, b::TensorValue)  = elwise(>=, a, b)

call(op::Sigmoid, a::TensorValue) = (1.0 ./ (1.0 + exp(-a)))

# size(a) == size(b) or only one of them is a scalar
# TODO: This is absolutely definitely not a good long term solution
@register_grad Add bg(ds, a) bg(ds, b)
@register_grad Sub bg(ds, a) bg(-ds, b)
@register_grad Mul bg(b * ds, a) bg(a * ds, b)
@register_grad Div bg(ds / b, a) bg(-(ds * a) / (b * b), b)
@register_grad Pow bg(ds * b * a ^ (b - cg.constant(1.0)), a) bg(ds * log(a) * (a ^ b), b)
# TODO: make the constants the proper type

@register_grad Sign fill_like(cg.constant(0), a)
@register_grad Neg (-ds)
@register_grad Exp (exp(a) * ds)
@register_grad Log (ds / a)
@register_grad Sin (cos(a) * ds)
@register_grad Cos (-sin(a) * ds)
@register_grad Abs (sign(a) * ds)

@register_grad Max bg(eq(out, a) * ds, a) bg(eq(out, b) * ds, b)
@register_grad Min bg(eq(out, a) * ds, a) bg(eq(out, b) * ds, b)

@register_grad Sigmoid (sigmoid(a) * (cg.constant(1.0) - sigmoid(a)) * ds)

################
# Broadcasting #
################
type RepeatTo <: OpType end
@register_op RepeatTo    repeatto     2 # input array, target shp - expand singleton dims
@register_grad RepeatTo broadcastgrad(ds, a) constant(0.0) # Nondiff

# Leave scalars as scalars if dim = ()
function call(op::RepeatTo, a::Real, dim::Array)
    length(dim) == 0 ? a : fill(a, dim...)
end

# Repeats singleton dimensions of the given array to given dimensions
function call(op::RepeatTo, a::Array, dim::Array)
    cdim = collect(size(a))
    ncurdim = length(cdim)
    @assert length(size(dim)) == 1  # Is a vector
    @assert length(dim) >= length(cdim)  # target has at least as many dimensions
    if cdim == dim
        return a
    end
    repcount = zeros(dim)
    for i = 1:length(dim)
        if i > ncurdim
            # Expand all implicitly 1 dimensions
            repcount[i] = dim[i]
        elseif cdim[i] == 1
            # Expand all singleton dims
            repcount[i] = dim[i]
        elseif cdim[i] == dim[i]
            # Leave rest the same
            repcount[i] = 1
        else
            # All all nonsingleton dimensions must be the same
            @assert false
        end
    end
    repeat(a, inner=repcount)
end

# This is used for arrays that are different sizes.
# Broadcast singleton dimensions to the same size

# Calculate a common shape by expanding singleton dimensions of arguments
type BroadcastShape <: OpType end
@register_op BroadcastShape broadcastshape 2
@register_grad BroadcastShape cg.constant(0.0) cg.constant(0.0)  # Really nondiff?

function call(op::BroadcastShape, a::Array, b::Real)
    collect(size(a))
end

function call(op::BroadcastShape, a::Real, b::Array)
    collect(size(b))
end

function call(op::BroadcastShape, a::Array, b::Array)
    as = size(a)
    bs = size(b)

    if length(bs) > length(as)
        as, bs = bs, as
        a, b = b, a
    end
    out = ones(Int, length(as))
    for i = 1:length(as)
        if i > length(bs)
            out[i] = as[i]
        elseif as[i] == bs[i]
            out[i] = as[i]
        elseif as[i] == 1
            out[i] = bs[i]
        elseif bs[i] == 1
            out[i] = as[i]
        else
            throw(DimensionMismatch("Cannot broadcast arrays"))
        end
    end
    out
end

function broadcastop(op::ScalarOp, a::Node, b::Node)
    common_shape = broadcastshape(a, b)
    ba = repeatto(a, common_shape)
    bb = repeatto(b, common_shape)
    Node(op, [ba, bb])
end

function broadcast(name::AbstractString, a::Node, b::Node)
    if name== "+"
        op = Add()
    elseif name == "-"
        op = Sub()
    elseif name == "*"
        op = Mul()
    elseif name == "/"
        op = Div()
    elseif name == "^"
        op = Pow()
    elseif name == "=="
        op = Eq()
    end
    broadcastop(op, a, b)
end

# Sums along axes that were broadcast to get singleton dims back
# args: [out, arg to go back to]
# For broadcastgrad(a, b) where a = [3,5], b=[1,5], sum(a, 1) to get [1,5] output
# In place until shape inference exists?
# Some way to specialize {T <: ScalarOp}?
type BroadcastGrad <: OpType
end

@register_op BroadcastGrad broadcastgrad 2
bg = broadcastgrad

# What axes to sum to get from cur_dim to target
function get_sum_dims(curdim::Vector{Int}, targetdim::Vector{Int})
    dims = Vector{Int64}()
    nt = length(targetdim)
    nc = length(curdim)
    @assert nt <= nc
    for dim = 1:nc
        if dim > nt
            # Sum along all implicitly 1 axes of target
            if curdim[dim] > 1
                push!(dims, dim)
            end
        elseif curdim[dim] != targetdim[dim] == 1
            push!(dims, dim)
        end
    end
    dims
end

# Sum axes of "from" to match dimensions of "to"
function call(op::BroadcastGrad, from::Real, to::Real)
    from
end

function call(op::BroadcastGrad, from::Real, to::Array)
    fill(from, size(to)...)
end

function call(op::BroadcastGrad, from::Array, to::Real)
    sum(from)
end

function call(op::BroadcastGrad, from::Array, to::Array)
    fshape = size(from)
    tshape = size(to)
    if fshape == tshape
        return from
    end
    dims = get_sum_dims(collect(fshape), collect(tshape))
    return sum(from, dims)
end

#TODO question: Is {T <: ScalarOp}(op::T) more efficient than op::ScalarOp?  Look at codegen
#TODO Just calling sin(x) is faster than broadcast(sin, x) by many times.  Do this instead for most ops

###############
# Complex ops #
###############

function crossentropy(label::Node, prediction::Node)
    result = -sum(label * log(prediction))
    group_between([label, prediction], [result], string(gensym(:crossentropy)), include_in=false)
    result
end

function softmax(node::Node)
    exped = exp(node)
    summed = sum(exped, constant(1))
    div = broadcast("/", exped, summed)
    group_between([node], [div], string(gensym(:softmax)), include_in=false)
    div
end

function softmax_stable(node::Node)
    max = maximum(node, constant(1))  # Maximum columnwise
    exped = exp(broadcast("-", node, max))
    summed = sum(exped, constant(1))
    div = broadcast("/", exped, summed)
    group_between([node], [div], string(gensym(:softmax)), include_in=false)
    div
end

# TODO: Add some numeric stability optimizations so we don't need this
function softmax_crossentropy(label::Node, unnorm_prediction::Node)
    max = maximum(unnorm_prediction)
    exped = exp(unnorm_prediction - max)
    summed = sum(exped, constant(1))
    lg = unnorm_prediction - log(summed)
    result = -sum(label * lg)
    group_between([label, unnorm_prediction], [result], string(gensym(:softmax_crossentropy)), include_in=false)
    result
end

function mean_squared_error(a::Node, b::Node)
    diff = a - b
    sq = diff * diff
    mean(sq)
end

################
# Optimization #
################

# Create an optimize op and return
function sgd_optimizer(loss::Node, variables::Vector{Node}, step_size::Node)
    gradients = grad(loss, variables)
    step_sizes = map(grad -> step_size * grad, gradients)
    updates = map(vargrad -> plusequals(vargrad[1], (-step_size * vargrad[2])), zip(variables, gradients))
    noop(updates...)
end


text:
module SmallLassoDemo

using VarianceReducedSG
using ProximalOperators
using LinearAlgebra
using Random
using Test

# Square-norm gradient
sqgrad(x) = x

function setup_data(N, n, mt)
	# Generate data
	A = randn(mt, n,N)
	y = A*randn(mt,N) # No noise

	# Normalize and reformat
	norm_consts = [norm(A[i,:]) for i = 1:size(A,1)]
	data = [A[i,:]./norm_consts[i] for i = 1:size(A,1)]
	b = [y[i]/norm_consts[i] for i = 1:length(y)]

	return data, b
end

################################################################################
# SAGA
################################################################################
@testset "Lasso" begin
	println("Lasso")
	N = 100
	n = 1000
	mt = MersenneTwister(123)

	data, b = setup_data(N,n,mt)

	f(y,x,i) = (y .= data[i].*sqgrad(dot(data[i],x) - b[i]))
	fs(x,i) = sqgrad(x - b[i])
	vrg = VRGradient(f, zeros(N), n)
	vrg_lin = LinearVRG(fs, data, zeros(N))
	vrg_uni = UniformVRG(f, zeros(N), n)
	reg = NormL1(.01)

	# Solve
	niter = n*100
	alg = SAGA(vrg, 10*randn(mt,N), 1/3, mt, prox_f=reg)
	iterate!(alg, niter)

	x_star = primiterates(alg)
	fval(x,y) = norm(x - x_star)
	progress = ShowFuncVal(fval, "|| x - x^*||")

	# Solve again for different params and solvers, all should give same result
	println("Prox-SAGA - VRGradient")
	niter = n*100
	nlogs = 10
	alg = SAGA(vrg, randn(mt,N), 1/10, mt, prox_f=reg)
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-SAGA - LinearVRG - Importance")
	niter = n*100
	nlogs = 10
	alg = SAGA(
		vrg_lin, randn(mt,N), 1/4, mt, prox_f=reg, weights=0.5.+rand(mt,n))
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-SAG - LinearVRG - Importance")
	niter = n*100
	nlogs = 10
	alg = SAG(
		vrg_lin, randn(mt,N), 1/4, mt, prox_f=reg, weights=0.5.+rand(mt,n))
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-SVAG - LinearVRG")
	niter = n*100
	nlogs = 10
	alg = SVAG(vrg_lin, randn(mt,N), 1/4, n/20, mt, prox_f=reg)
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-ASVAG - LinearVRG")
	niter = n*100
	nlogs = 10
	alg = ASVAG(vrg_lin, randn(mt,N), 1/4, mt, prox_f=reg)
	iterate!(alg, niter,
		(progress,
			ShowAlgState(a->a.innoweight[], "IW"),
			ShowNewLine()
		),
		nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-SVRG - UniformVRG")
	stagelen = n/10
	niter = 50*n
	nlogs = 10
	alg = SVRG(vrg_uni, randn(mt,N), 1/8, stagelen, mt, prox_f=reg)
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-LSVRG - UniformVRG - Importance")
	q = 10/n
	niter = n*50
	nlogs = 10
	alg = LSVRG(
		vrg_uni, randn(mt,N), 1/8, q, mt, prox_f=reg,weights=0.5.+rand(mt,n))
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-ILSVRG - VRGradient - Importance")
	q = 10/n
	niter = n*50
	nlogs = 10
	alg = ILSVRG(
		vrg, randn(mt,N), 1/8, q, mt, prox_f=reg,weights=0.5.+rand(mt,n))
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-QSAGA - LinearVRG - q = 10")
	q = 10
	niter = n*50
	nlogs = 10
	alg = QSAGA(
		vrg_lin, randn(mt,N), 1/8, q, mt, prox_f=reg)
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)

	println("Prox-QSAGA - VRGadient - q = 10 - Without replacement")
	q = 10
	niter = n*50
	nlogs = 10
	alg = QSAGA(
		vrg, randn(mt,N), 1/8, q, mt, prox_f=reg, replace=false)
	iterate!(alg, niter, (progress, ShowNewLine()), nlogs)
	@test isapprox(primiterates(alg), x_star)
end

end

ext: jl


text:
"""
    println_result(msg...)

Output a result so that it stands out.
"""
println_result(msg...) = printstyled(msg..., "\n"; color=:green, bold=true)


text:
abstract type AbstractView{T,N} <: AbstractArray{T,N} end

@inline Base.size(view::AbstractView) = size(parent(view))

@inline Base.eltype(view::AbstractView{T,N}) where {T,N} = T

@inline Base.ndims(view::AbstractView{T,N}) where {T,N} = N

function Base.IndexStyle(type::Type{S}) where {S <: AbstractView}
    index_type(S)
end

@inline Base.axes(view::AbstractView{T,N}) where {T,N} = axes(parent(view))


text:
#=
fld-write.jl
write AVS .fld to file
2019-05-12, Jeff Fessler, University of Michigan
=#

export fld_write

"""
    fld_write(file, data ; kwargs...)

Write data into AVS format `.fld` file.
See README for file format.

in
- `file` name of file typically ending in `.fld`
- `data` real data array

option
- `check::Bool`         report error if file exists; default `true`
- `dir::String`         directory name to prepend file name; default `""`
- `endian::`Symbol`     `:le` little endian (default), `:be` big endian
- `head::Array{String}` comment information for file header
- `raw::Bool`           put raw data in `name.raw`, header in `name.fld`
                        where `file` = `name.fld`; default `false`
"""
function fld_write(
    file::String,
    data::AbstractArray{<:Real} ;
    check::Bool = true,
    dir::String = "",
    endian::Symbol = :le,
    head::Array{String} = empty([""]),
    warn::Bool = true,
    raw::Bool = false,
)

    data = fld_write_data_fix(data ; warn) # make data suitable for writing in .fld

    endian != :le && endian != :be && throw("endian '$endian' unknown")

    typedict = Dict([
        (Float32, endian === :le ? "float_le" : "xdr_float"),
        (Float64, endian === :le ? "double_le" : "xdr_double"),
        (UInt8, "byte"),
        (Int16, endian === :le ? "short_le" : "short_be"),
        (Int32, endian === :le ? "int_le" : "xdr_int"),
    ])

    datatype = typedict[eltype(data)] # throws error if not there

    file = joinpath(dir, file)
#   @show file

    check && isfile(file) && throw("file '$file' exists")

    if raw # if writing data to separate ".raw" file, ensure it does not exist
        fileraw = file[1:(end-3)] * "raw"
        check && isfile(fileraw) && throw("file $fileraw exists")
    end

    # write header to IOBuffer
    io = IOBuffer()

    ndim = ndims(data)

    println(io, "# AVS field file ($(basename(@__FILE__)))")
    for line in head
        println(io, "# $line")
    end

    println(io, "ndim=$ndim")
    for ii=1:ndim
        println(io, "dim$ii=$(size(data,ii))")
    end
    println(io, "nspace=$ndim")
    println(io, "veclen=1")
    println(io, "data=$datatype")
    println(io, "field=uniform")

    if raw
        println(io, "variable 1 file=$fileraw filetype=binary")
    else
        write(io, "\f\f") # two form feeds: char(12)
    end

    # determine how to write the binary data
    host_is_le = () -> ENDIAN_BOM == 0x04030201
    fun = (host_is_le() == (endian === :le)) ?
        identity : # host/file same endian
        (host_is_le() && (endian === :be)) ?
        hton :
        (!host_is_le() && (endian === :le)) ?
        htol :
        throw("not done")

    # write header from IO buffer to file
    open(file, "w") do fid
        write(fid, take!(io))
    end

    # write binary data to file or fileraw
    if raw
        open(fileraw, "w") do fraw
            write(fraw, fun.(data))
        end
    else
        open(file, "a") do fid
            write(fid, fun.(data))
        end
    end

	return nothing
end

"""
    data = fld_write_data_fix(data)
Convert data to format suitable for writing to `.fld` file.
"""
function fld_write_data_fix(data::AbstractArray{BigFloat} ; warn::Bool=true)
    warn && (@warn "BigFloat downgraded to Float64")
    return Float64.(data)
end

function fld_write_data_fix(data::AbstractArray{Float16} ; warn::Bool=true)
    warn && (@warn "Float16 promoted to Float32")
    return Float32.(data)
end

function fld_write_data_fix(
    data::AbstractArray{T} ;
    warn::Bool=true,
) where {T <: Union{BigInt, Int64}}
    warn && (@warn "$T downgraded to Int32")
    return Int32.(data)
end

function fld_write_data_fix(data::AbstractArray{Bool} ; warn::Bool=true)
    warn && (@warn "Bool promoted to UInt8")
    return UInt8.(data)
end

@inline fld_write_data_fix(data::AbstractArray{T} ; warn::Bool=true) where
    {T <: Union{Float32, Float64, UInt8, Int16, Int32}} = data

fld_write_data_fix(data::AbstractArray{T} ; warn::Bool=true) where {T <: Any} =
        throw(ArgumentError("unsupported type $T"))

max_stars_repo_name: JeffFessler/AVSfldIO

max_issues_repo_name: JeffFessler/AVSfldIO


text:
# Use this script to generate grid files and save them into a given folder
using Lehmann
using Printf

# rtol = [-6, -8, -10, -12, -14]
# rtol = [-7, -9, -11, -13]
# Λ = [100, 1000, 10000, 100000, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14]
# Λ = [1e11, 1e12, 1e13, 1e14]
rtol = [-12]
Λ = [1000000,]
algorithm = :functional
# folder = "./basis/"
folder = "./"

for lambda in Λ
    for err in rtol
        if lambda <= 1e8 #the universal grid beyond 1e8 suffers from the accuracy loss, so that we turn it off
            DLRGrid(1.0, lambda, 10.0^err, true, :none, rebuild = true, algorithm = algorithm, folder = folder, verbose = true)
        end
        DLRGrid(1.0, lambda, 10.0^err, true, :ph, rebuild = true, algorithm = algorithm, folder = folder, verbose = true)
        DLRGrid(1.0, lambda, 10.0^err, true, :pha, rebuild = true, algorithm = algorithm, folder = folder, verbose = true)
    end
end


text:
using Distributed
@everywhere using MCPhylo

## Data
lsat = Dict{Symbol, Any}(
  :culm =>
    [3, 9, 11, 22, 23, 24, 27, 31, 32, 40, 40, 56, 56, 59, 61, 76, 86, 115, 129,
     210, 213, 241, 256, 336, 352, 408, 429, 602, 613, 674, 702, 1000],
  :response =>
    [0 0 0 0 0
     0 0 0 0 1
     0 0 0 1 0
     0 0 0 1 1
     0 0 1 0 0
     0 0 1 0 1
     0 0 1 1 0
     0 0 1 1 1
     0 1 0 0 0
     0 1 0 0 1
     0 1 0 1 0
     0 1 0 1 1
     0 1 1 0 0
     0 1 1 0 1
     0 1 1 1 0
     0 1 1 1 1
     1 0 0 0 0
     1 0 0 0 1
     1 0 0 1 0
     1 0 0 1 1
     1 0 1 0 0
     1 0 1 0 1
     1 0 1 1 0
     1 0 1 1 1
     1 1 0 0 0
     1 1 0 0 1
     1 1 0 1 0
     1 1 0 1 1
     1 1 1 0 0
     1 1 1 0 1
     1 1 1 1 0
     1 1 1 1 1],
  :N => 1000
)
lsat[:R] = size(lsat[:response], 1)
lsat[:T] = size(lsat[:response], 2)

n = [lsat[:culm][1]; diff(lsat[:culm])]
idx = mapreduce(i -> fill(i, n[i]), vcat, 1:length(n))
lsat[:r] = lsat[:response][idx, :]

## Model Specification

model = Model(

  r = Stochastic(2,
    (beta, theta, alpha, N, T) ->
      UnivariateDistribution[(
        p = invlogit(beta * theta[i] - alpha[j]);
        Bernoulli(p)) for i in 1:N, j in 1:T
      ],
    false
  ),

  theta = Stochastic(1,
    () -> Normal(0, 1),
    false
  ),

  alpha = Stochastic(1,
    () -> Normal(0, 100),
    false
  ),

  a = Logical(1,
    alpha -> alpha .- mean(alpha)
  ),

  beta = Stochastic(
    () -> Truncated(Flat(), 0, Inf)
  )

)

## Initial Values
inits = [
  Dict(:r => lsat[:r], :alpha => zeros(lsat[:T]), :beta => 1,
       :theta => zeros(lsat[:N])),
  Dict(:r => lsat[:r], :alpha => ones(lsat[:T]), :beta => 2,
       :theta => zeros(lsat[:N]))
]

## Sampling Scheme
scheme = [AMWG(:alpha, 0.1),
          Slice(:beta, 1.0),
          Slice(:theta, 0.5)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, lsat, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
function gsreg(
        depvar,
        expvars,
        data;
        intercept=nothing,
        outsample=nothing,
        criteria=nothing,
        ttest=nothing,
        vectoroperation=nothing,
        modelavg=nothing,
        residualtest=nothing,
        time=nothing,
        summary=nothing,
        datanames=nothing,
        datatype=nothing,
        orderresults=nothing,
        onmessage=nothing,
        parallel=nothing
    )
    result = GSRegResult(
        depvar,
        expvars,
        data,
        intercept,
        outsample,
        criteria,
        ttest,
        vectoroperation,
        modelavg,
        residualtest,
        time,
        datanames,
        datatype,
        orderresults,
        onmessage,
        parallel
    )
    proc!(result)
    if summary != nothing
        f = open(summary, "w")
        write(f, to_string(result))
        close(f)
    end
    return result
end

function gsreg_single_proc_result!(
        order,
        data,
        results,
        depvar,
        expvars,
        intercept,
        outsample,
        criteria,
        ttest,
        vectoroperation,
        residualtest,
        time,
        datanames,
        datatype,
        header)

    cols = get_selected_cols(order)
    data_cols_num = size(data, 2)

    if intercept
        append!(cols, data_cols_num)
    end

    depvar_data = @view(data[1:end-outsample, 1])
    expvars_data = @view(data[1:end-outsample, cols])

    nobs = size(depvar_data, 1)
    ncoef = size(expvars_data, 2)
    qrf = qr(expvars_data)
    b = qrf \ depvar_data                   # estimate
    ŷ = expvars_data * b                    # predicted values
    er = depvar_data - ŷ                    # in-sample residuals
    er2 = er .^ 2                           # squared errors
    sse = sum(er2)                          # residual sum of squares
    df_e = nobs - ncoef                     # degrees of freedom
    rmse = sqrt(sse / nobs)                 # root mean squared error
    r2 = 1 - var(er) / var(depvar_data)     # model R-squared

    if ttest
        bstd = sqrt.( sum( (UpperTriangular(qrf.R) \ Matrix(1.0LinearAlgebra.I, ncoef, ncoef) ) .^ 2, dims=2) * (sse / df_e) ) # std deviation of coefficients
    end

    if outsample > 0
        depvar_out = @view(data[end-outsample:end, 1])
        expvars_out = @view(data[end-outsample:end, cols])
        erout = depvar_out - expvars_out * b          # out-of-sample residuals
        sseout = sum(erout .^ 2)                      # residual sum of squares
        rmseout = sqrt(sseout / outsample)            # root mean squared error
        results[order, header[:rmseout]] = rmseout
    end

    results[order, header[:index]] = order
    for (index, col) in enumerate(cols)
        header[Symbol(string(expvars[col-1], "_b"))]
        results[order, header[Symbol(string(expvars[col-1], "_b"))]] = datatype(b[index])
        if ttest
            results[order, header[Symbol(string(expvars[col-1], "_bstd"))]] = datatype(bstd[index])
        end
    end

    results[order, header[:nobs]] = nobs
    results[order, header[:ncoef]] = ncoef
    results[order, header[:sse]] = datatype(sse)
    results[order, header[:r2]] = datatype(r2)
    results[order, header[:rmse]] = datatype(rmse)
    results[order, header[:order]] = 0

    if residualtest != nothing && residualtest

        #jbtest
        x = er
        n = length(x)
        m1 = sum(x) / n
        m2 = sum((x .- m1) .^ 2) / n
        m3 = sum((x .- m1) .^ 3) / n
        m4 = sum((x .- m1) .^ 4) / n
        b1 = (m3 / m2 ^ (3 / 2)) ^ 2
        b2 = (m4 / m2 ^ 2)
        statistic = n * b1 / 6 + n * (b2 - 3) ^ 2 / 24
        d = Chisq(2.)
        jbtest = 1 .- cdf(d, statistic)

        #wtest
        regmatw = hcat((ŷ .^ 2), ŷ , ones(size(ŷ, 1)))
        qrfw = qr(regmatw)
        regcoeffw = qrfw \ er2
        residw = er2 - regmatw * regcoeffw
        rsqw = 1 - dot(residw, residw) / dot(er2, er2) # uncentered R^2
        statisticw = n * rsqw
        wtest = ccdf(Chisq(2), statisticw)

        results[order, header[:wtest]] = wtest
        results[order, header[:jbtest]] = jbtest
        if time != nothing
            e = er
            lag = 1
            xmat = expvars_data

            n = size(e,1)
            elag = zeros(Float64,n,lag)
            for ii = 1:lag  # construct lagged residuals
                elag[ii+1:end, ii] = e[1:end-ii]
            end

            offset = lag
            regmatbg = [xmat[offset+1:end,:] elag[offset+1:end,:]]
            qrfbg = qr(regmatbg)
            regcoeffbg = qrfbg \ e[offset+1:end]
            residbg = e[offset+1:end] .- regmatbg * regcoeffbg
            rsqbg = 1 - dot(residbg,residbg) / dot(e[offset+1:end], e[offset+1:end]) # uncentered R^2
            statisticbg = (n - offset) * rsqbg
            bgtest = ccdf(Chisq(lag), statisticbg)
            results[order, header[:bgtest]] = bgtest
        end
    end

    if vectoroperation == false
        if ttest
            for (index, col) in enumerate(cols)
                pos_b = header[Symbol(string(expvars[col-1], "_b"))]
                pos_bstd = header[Symbol(string(expvars[col-1], "_bstd"))]
                pos_t = header[Symbol(string(expvars[col-1], "_t"))]
                results[order, pos_t] = results[order, pos_b] / results[order, pos_bstd]
            end
        end

        if :aic in criteria
            results[order,header[:aic]] = 2 * results[order, header[:ncoef]] + results[order, header[:nobs]] * log(results[order,header[:sse]] / results[order, header[:nobs]])
        end

        if :aicc in criteria
            if :aic in criteria
                results[order, header[:aicc]] = results[order, header[:aic]] + (2(results[order, header[:ncoef]] + 1) * (results[order, header[:ncoef]]+2)) / (results[order, header[:nobs]] - (results[order, header[:ncoef]] + 1 ) - 1)
            else
                aic = 2 * results[order, header[:ncoef]] + results[order, header[:nobs]] * log(results[order, header[:sse]] / results[order, header[:nobs]])
                results[order, header[:aicc]] = aic + (2(results[order, header[:ncoef]] + 1) * (results[order, header[:ncoef]]+2)) / (results[order, header[:nobs]] - (results[order, header[:ncoef]] + 1 ) - 1)
            end
        end

        if :bic in criteria
            results[order, header[:bic]] = results[order, header[:nobs]] * log.(results[order,header[:rmse]]) + ( results[order, header[:ncoef]] - 1 ) * log.(results[order, header[:nobs]]) + results[order, header[:nobs]] + results[order, header[:nobs]] * log(2π)
        end

        if :r2adj in criteria
            results[order,header[:r2adj]] = 1 - (1 - results[order, header[:r2]]) * ((results[order, header[:nobs]] - 1) / (results[order, header[:nobs]] - results[order, header[:ncoef]]))
        end

        results[order, header[:F]] = (results[order, header[:r2]] / (results[order,header[:ncoef]] - 1)) / ((1 - results[order, header[:r2]]) / (results[order, header[:nobs]] - results[order, header[:ncoef]]))
    end
end

function gsreg_proc_result!(num_job, num_jobs, ops_by_worker, opts...)
    for j = 1:ops_by_worker
        order = (j-1) * num_jobs + num_job
        gsreg_single_proc_result!(order, opts...)
    end
end

function proc!(result::GSRegResult)
    expvars_num = size(result.expvars, 1)
    if result.intercept
        expvars_num = expvars_num-1
    end
    num_operations = 2 ^ expvars_num - 1

    result.onmessage("Creating environment")

    pdata = convert(SharedArray, result.data)
    presults = fill!(SharedArray{result.datatype}(num_operations, length(keys(result.header))),NaN)

    result.onmessage("Doing $num_operations regressions")

    if nprocs() == nworkers()
        for order = 1:num_operations
            gsreg_single_proc_result!(order, pdata, presults, result.depvar, result.expvars, result.intercept, result.outsample, result.criteria, result.ttest, result.vectoroperation,  result.residualtest, result.time, result.datanames, result.datatype, result.header)
        end
    else
        num_workers = (result.parallel != nothing) ? result.parallel : nworkers()
        ops_by_worker = div(num_operations, num_workers)
        num_jobs = (num_workers > num_operations) ? num_operations : num_workers
        remainder = num_operations - ops_by_worker * num_jobs
        jobs = []
        for num_job = 1:num_jobs
            push!(jobs, @spawnat num_job+1 gsreg_proc_result!(num_job, num_jobs, ops_by_worker, pdata, presults, result.depvar, result.expvars, result.intercept, result.outsample, result.criteria, result.ttest, result.vectoroperation, result.residualtest, result.time, result.datanames, result.datatype, result.header))
        end
        for job in jobs
            fetch(job)
        end

        if( remainder > 0 )
            for j = 1:remainder
                order = j + ops_by_worker * num_jobs
                gsreg_single_proc_result!(order, pdata, presults, result.depvar, result.expvars, result.intercept, result.outsample, result.criteria, result.ttest, result.vectoroperation, result.residualtest, result.time, result.datanames, result.datatype, result.header)
            end
        end
    end

    result.results = Array(presults)

    presult = nothing
    pdata = nothing

    if result.vectoroperation
        result.onmessage("Calculating Vector operations")
        if result.ttest
            for expvar in result.expvars
                pos_b = result.header[Symbol(string(expvar, "_b"))]
                pos_bstd = result.header[Symbol(string(expvar, "_bstd"))]
                pos_t = result.header[Symbol(string(expvar, "_t"))]
                result.results[:,pos_t] = result.results[:,pos_b] ./ result.results[:,pos_bstd]
            end
        end

        if :aic in result.criteria || :aicc in result.criteria
            result.results[:,result.header[:aic]] = 2 * result.results[:,result.header[:ncoef]] .+ result.results[:,result.header[:nobs]] .* log.(result.results[:,result.header[:sse]] ./ result.results[:,result.header[:nobs]])
        end

        if :aicc in result.criteria
            result.results[:,result.header[:aicc]] = result.results[:,result.header[:aic]] .+ (2(result.results[:,result.header[:ncoef]] .+ 1) .* (result.results[:,result.header[:ncoef]] .+ 2)) ./ (result.results[:,result.header[:nobs]] .- (result.results[:,result.header[:ncoef]] .+ 1 ) .- 1)
        end

        if :bic in result.criteria
            result.results[:,result.header[:bic]] = result.results[:,result.header[:nobs]] .* log.(result.results[:,result.header[:rmse]]) .+ ( result.results[:,result.header[:ncoef]] .- 1 ) .* log.(result.results[:,result.header[:nobs]]) .+ result.results[:,result.header[:nobs]] .+ result.results[:,result.header[:nobs]] .* log(2π)
        end

        if :r2adj in result.criteria
            result.results[:,result.header[:r2adj]] = 1 .- (1 .- result.results[:,result.header[:r2]]) .* ((result.results[:,result.header[:nobs]] .- 1) ./ (result.results[:,result.header[:nobs]] .- result.results[:,result.header[:ncoef]]))
        end

        result.results[:,result.header[:F]] = (result.results[:,result.header[:r2]] ./ (result.results[:,result.header[:ncoef]] .- 1)) ./ ((1 .- result.results[:,result.header[:r2]]) ./ (result.results[:,result.header[:nobs]] .- result.results[:,result.header[:ncoef]]))
    end

    # CP must be computed with vector operations
    if :cp in result.criteria
        result.results[:,result.header[:cp]] = (result.results[:,result.header[:nobs]] .- maximum(result.results[:,result.header[:ncoef]]) .- 2) .* (result.results[:,result.header[:rmse]] ./ minimum(result.results[:,result.header[:rmse]])) .- (result.results[:,result.header[:nobs]] .- 2 .* result.results[:,result.header[:ncoef]])
    end
    len_criteria = length(result.criteria)
    for criteria in result.criteria
        result.results[:, result.header[:order]] += AVAILABLE_CRITERIA[criteria]["index"] * (1 / len_criteria) * ( (result.results[:, result.header[criteria]] .- mean(result.results[:, result.header[criteria]]) ) ./ std(result.results[:, result.header[criteria]]) )
    end
    if result.modelavg
        result.onmessage("Generating model averaging results")
        # usar order para weight
        delta = maximum(result.results[:,result.header[:order]]) .- result.results[:,result.header[:order]]
        w1 = exp.(-delta/2)
        result.results[:,result.header[:weight]] = w1./sum(w1)
        result.average = Array{Float64}(undef, 1, length(keys(result.header)))
        weight_pos = (result.ttest) ? 4 : 2
        for expvar in result.expvars
            obs = result.results[:,result.header[Symbol(string(expvar, "_b"))]]
            if result.ttest
                obs = hcat(obs, result.results[:,result.header[Symbol(string(expvar, "_bstd"))]])
                obs = hcat(obs, result.results[:,result.header[Symbol(string(expvar, "_t"))]])
            end
            obs = hcat(obs, result.results[:,result.header[:weight]])

            #filter NaN values from selection
            obs = obs[findall(x -> !isnan(obs[x,1]), 1:size(obs,1)),:]

            #weight resizing
            obs[:, weight_pos] /= sum(obs[:, weight_pos])

            result.average[result.header[Symbol(string(expvar, "_b"))]] = sum(obs[:, 1] .* obs[:, weight_pos])
            if result.ttest
                result.average[result.header[Symbol(string(expvar, "_bstd"))]] = sum(obs[:, 2] .* obs[:, weight_pos])
                result.average[result.header[Symbol(string(expvar, "_t"))]] = sum(obs[:, 3] .* obs[:, weight_pos])
            end
        end

        for criteria in [:nobs, :r2adj, :F, :order]
            result.average[result.header[criteria]] = sum(result.results[:, result.header[criteria]] .* result.results[:, result.header[:weight]])
        end
    end

    if result.orderresults
        result.onmessage("Sorting results")
        result.results = gsregsortrows(result.results, [result.header[:order]]; rev=true)
        result.bestresult = result.results[1,:]
    else
        result.onmessage("Looking for the best result")
        max_order = result.results[1,result.header[:order]]
        best_result_index = 1
        for i = 1:num_operations
            if result.results[i,result.header[:order]] > max_order
                max_order = result.results[i,result.header[:order]]
                best_result_index = i
            end
        end
        result.bestresult = result.results[best_result_index,:]
    end
end

function to_string(result::GSRegResult)
    out = ""
    out *= @sprintf("\n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                              Best model results                              \n")
    out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf("                                     Dependent variable: %s                   \n", result.depvar)
    out *= @sprintf("                                     ─────────────────────────────────────────\n")
    out *= @sprintf("                                                                              \n")
    out *= @sprintf(" Selected covariates                 Coef.")
    if result.ttest
        out *= @sprintf("        Std.         t-test")
    end
    out *= @sprintf("\n")
    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    cols = get_selected_cols(Int64(result.bestresult[result.header[:index]]))

    data_cols_num = size(result.data, 2)
    if result.intercept
        append!(cols, data_cols_num)
    end

    for pos in cols
        varname = result.expvars[pos-1]
        out *= @sprintf(" %-35s", varname)
        out *= @sprintf(" %-10f", result.bestresult[result.header[Symbol(string(varname, "_b"))]])
        if result.ttest
            out *= @sprintf("   %-10f", result.bestresult[result.header[Symbol(string(varname, "_bstd"))]])
            out *= @sprintf("   %-10f", result.bestresult[result.header[Symbol(string(varname, "_t"))]])
        end
        out *= @sprintf("\n")
    end

    out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    out *= @sprintf(" Observations                        %-10d\n", result.bestresult[result.header[:nobs]])
    out *= @sprintf(" F-statistic                         %-10f\n", result.bestresult[result.header[:F]])
    for criteria in result.criteria
        if AVAILABLE_CRITERIA[criteria]["verbose_show"]
    out *= @sprintf(" %-30s      %-10f\n", AVAILABLE_CRITERIA[criteria]["verbose_title"], result.bestresult[result.header[criteria]])
        end
    end

    if !result.modelavg
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
    end

    if result.modelavg
        out *= @sprintf("\n")
        out *= @sprintf("\n")
        out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
        out *= @sprintf("                            Model averaging results                           \n")
        out *= @sprintf("══════════════════════════════════════════════════════════════════════════════\n")
        out *= @sprintf("                                                                              \n")
        out *= @sprintf("                                     Dependent variable: %s                   \n", result.depvar)
        out *= @sprintf("                                     ─────────────────────────────────────────\n")
        out *= @sprintf("                                                                              \n")
        out *= @sprintf(" Covariates                          Coef.")
        if result.ttest
            out *= @sprintf("        Std.         t-test")
        end
        out *= @sprintf("\n")
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

        for varname in result.expvars
            out *= @sprintf(" %-35s", varname)
            out *= @sprintf(" %-10f", result.average[result.header[Symbol(string(varname, "_b"))]])
            if result.ttest
                out *= @sprintf("   %-10f", result.average[result.header[Symbol(string(varname, "_bstd"))]])
                out *= @sprintf("   %-10f", result.average[result.header[Symbol(string(varname, "_t"))]])
            end
            out *= @sprintf("\n")
        end
        out *= @sprintf("\n")
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")
        out *= @sprintf(" Observations                        %-10d\n", result.average[result.header[:nobs]])
        out *= @sprintf(" F-statistic                         %-10f\n", result.average[result.header[:F]])
        out *= @sprintf(" Combined criteria                   %-10f\n", result.average[result.header[:order]])
        out *= @sprintf("──────────────────────────────────────────────────────────────────────────────\n")

    end
    return out
end

function to_dict(res::GSRegResult)
    result = Dict()
    for (header, index) in res.header
        push!(result, Pair(string(header), res.bestresult[index]))
    end

    dic = Dict(
        "selected_cols" => map( x -> res.datanames[x], get_selected_cols(Int64(res.bestresult[res.header[:index]]))),
        "bestresult" => result
    )

    if res.modelavg
        average = Dict()
        for (header, index) in res.header
            push!(average, Pair(string(header), res.average[index]))
        end
        push!(dic, Pair("avgresults", average))
    end

    dic
end

Base.show(io::IO, result::GSRegResult) = print(io, to_string(result))


text:
@testset "CLBLAS.L1" begin
    @testset "CLBLAS axpy" begin
        for T in testtypes
            alpha = T(2.0)
            hx = rand(T, 32)
            hy = rand(T, 32)
            x = CLArray(q, hx)
            y = CLArray(q, hy)

            cl.wait(axpy!(alpha, x, y; queue = q))
            axpy!(alpha, hx, hy)

            @test cl.to_host(y) ≈ hy
        end
    end

    @testset "CLBLAS scal" begin
        for T in [Float32,
                  Float64,
                  # Complex64,
                  Complex128]

            n = 32
            DA = T(2.0)
            hDX = rand(T, 32)
            DX = CLArray(q, hDX)
            incx = 1

            cl.wait(scal!(n, DA, DX, incx; queue=q))
            scal!(n, DA, hDX, incx)

            @test cl.to_host(DX) ≈ hDX
        end
    end
end


text:
@testset "Ambiguous nucleotide positions" begin
    function test_ns(A, seq)
        expected = Int[]
        for i in 1:length(seq)
            if seq[i] == 'N'
                push!(expected, i)
            end
        end
        bioseq = BioSequence{A}(seq)
        @test collect(ambiguous_positions(bioseq)) == expected
    end

    for len in [1, 10, 32, 1000, 10000, 100000]
        test_ns(DNAAlphabet{4}, random_dna(len))
        test_ns(RNAAlphabet{4}, random_rna(len))

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_ns(DNAAlphabet{2}, random_dna(len, probs))
        test_ns(RNAAlphabet{2}, random_rna(len, probs))
    end

    dna_seq = dna"ANANANA"
    pos = [2, 4, 6]
    for (i, p) in enumerate(ambiguous_positions(dna_seq))
        @test p == pos[i]
    end

    dna_seq = dna"NATTCGRATY"
    pos = [1, 7, 10]
    for (i, p) in enumerate(ambiguous_positions(dna_seq))
        @test p == pos[i]
    end

    rna_seq = rna"ANANANA"
    pos = [2, 4, 6]
    for (i, p) in enumerate(ambiguous_positions(rna_seq))
        @test p == pos[i]
    end

    rna_seq = rna"NAUUCGRAUY"
    pos = [1, 7, 10]
    for (i, p) in enumerate(ambiguous_positions(rna_seq))
        @test p == pos[i]
    end
end


text:
# test_local_channels2.jl

using Pkg; Pkg.activate(".")

# remove old stuff
cmd = `rm test_invoicing.txt test_invoicing_paid.txt test_ledger.txt test_journal.txt invoicenbr.txt`
run(cmd)

# enable distrbuted computing
using Distributed
@info("Enable distributed computing")

# this should be the next step
np = addprocs(4; exeflags=`--project=$(Base.active_project())`)
@info("number of processes is $(length(np))")

# activate the packages (before the processes are created)
@everywhere begin
    using AppliSales
    using AppliGeneralLedger
    using AppliAR
end;

@info("Distributed computing enabled")

# get the tasks and dispatcher
include("./api/api2.jl")

# start the dispatcher
rx = dispatcher()
@info("Dispatcher started")

# start application remote
include("client.jl")

# display aging report
using DataFrames

sleep(10)
r = DataFrame(report())
println("\nUnpaid invoices\n$("="^15)")
println(r)

println("")
include("check_general_ledger.jl")


text:
using Distributions, SpecialFunctions, FastGaussQuadrature

struct PositiveStable <: ContinuousUnivariateDistribution
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  PositiveStable(α::Real,β::Real) = ( β < -1 || β > 1 ||
    (β == -1 && α <= 1) || α <= 0 || α > 2 ) ?
    error("Parameters' requirements unmet:\n (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") :
    α == 2 ? new(2,0,0,.5) :
    new(α,β,β*(α <= 1 ? 1 :
    (α-2)/α), (1 + β*(α <= 1 ? 1 : (α-2)/α))/2)
end

import Distributions.params
function params(d::PositiveStable)
  return (d.α,d.β,d.θ,d.ρ)
end

import Distributions.minimum
function minimum(d::PositiveStable)
  return β == 1 && α == 1 ? 1 : 0
end

import Distributions.maximum
function maximum(d::PositiveStable)
  return β == 1 && α == 1 ? 1 : Inf
end

import Distributions.insupport
function insupport(d::PositiveStable,x::Real)
  return β == 1 && α == 1 ? x==1 : x >= 0
end

########################
# BASIC SIMULATION TOOLS
########################

import Distributions.rand
function rand(d::PositiveStable)
  ar = d.α*d.ρ
  if d.β == 1 && d.α == 1
    return 1.
  elseif d.β == 1 && d.α < 1
    u1 = rand(Uniform(0,pi))
    e1 = rand(Exponential())
    s1 = sin(u1)
    return (sin(d.α * u1) / s1) *
      (s1 * e1 / sin((1 -d.α) * u1)) ^ (1 - 1 /d.α)
  elseif d.β == -1 && d.α > 1
    u1 = rand(Uniform(0,pi))
    e1 = rand(Exponential())
    s1 = sin(u1)
    return ((sin(d.ρ * u1) / s1) *
      (s1 * e1 / sin((1 - d.ρ) * u1)) ^ (1 - 1 /d.ρ) ) ^ (-d.ρ)
  else
    u1 = rand(Uniform(0,pi))
    u2 = rand(Uniform(0,pi))
    e1 = rand(Exponential())
    e2 = rand(Exponential())
    s1 = sin(u1)
    s2 = sin(u2)
    return ( ((sin(ar * u1) / s1) *
      (s1 * e1 / sin((1 - ar) * u1)) ^ (1 - 1 /ar)) /
      ((sin(d.ρ * u2) / s2) *
      (s2 * e2 / sin.((1 - d.ρ) * u2)) ^ (1 - 1 /d.ρ)) ) ^ d.ρ
  end
end

function rand(d::PositiveStable,n::Int)
  ar = d.α*d.ρ
  if d.β == 1 && d.α == 1
    return ones(n)
  elseif d.β == 1 && d.α < 1
    u1 = rand(Uniform(0,pi),n)
    e1 = rand(Exponential(),n)
    s1 = sin.(u1)
    return (sin.(d.α .* u1) ./ s1) .*
      (s1 .* e1 ./ sin.((1-d.α) .* u1)) .^ (1-1/d.α)
  elseif d.β == -1 && d.α > 1
    u1 = rand(Uniform(0,pi),n)
    e1 = rand(Exponential(),n)
    s1 = sin.(u1)
    return ((sin.(d.ρ .* u1) ./ s1) .*
      (s1 .* e1 ./ sin.((1 - d.ρ) .* u1)) .^ (1 - 1/d.ρ) ) .^ (-d.ρ)
  else
    u1 = rand(Uniform(0,pi),n)
    u2 = rand(Uniform(0,pi),n)
    e1 = rand(Exponential(),n)
    e2 = rand(Exponential(),n)
    s1 = sin.(u1)
    s2 = sin.(u2)
    return ( ((sin.(ar .* u1) ./ s1) .*
      (s1 .* e1 ./ sin.((1 -ar) .* u1)) .^ (1 -1/ar)) ./
      ((sin.(d.ρ .* u2) ./ s2) .*
      (s2 .* e2 ./ sin.((1 -d.ρ) .* u2)) .^ (1 -1/d.ρ)) ) .^ d.ρ
  end
end

########################
# PDF, CDF
########################

function auxV2(x::AbstractArray,a::Real,b::Real)
  y = (pi/2).*x
  t = (pi/2)*a*b
  return ((sin.(a.*y .+ t).^a)./cos.(y)) .^ (1 /(1 -a)) .* cos.((a-1) .* y .+ t)
end

import Distributions.pdf
function pdf(d::PositiveStable,x::Real)
  if d.α <= 1 && β == 1
    return x == 1 ? 1. : 0.
  end
  if x < 0
    return 0.
  elseif x == 0
    return gamma(1 + 1 / d.α)*sin(pi*d.ρ)/(d.ρ*pi)
  end
  m = 10^3
  l = 15
  delta = 1
    #
  pir = pi*d.ρ
  pir1 = pi*(1 - d.ρ)
  #

  if d.α > 1
    a0 = d.α/(abs(d.α-1)*2)
    raa = 1 / (d.α-1)
    a1 = x == 0 ? 0 : abs(x)^raa
    a2 = -a1 .^ d.α
    s1 = d.ρ
    s2 = 1 - d.ρ
    a01 = a0*s1
    a02 = a0*s2
    if x<delta
      v = 1:Int(floor(l*d.α))
      w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v ./ d.α .+ 1) .*
        sin.(pir .* v) ./ (pi .* gamma.(v .+ 1))
      return sum(w .* abs(x) .^ (v .- 1))
    else
      nodes, weights = gausslegendre(m)
      seq1 = auxV2(nodes .* s1 .+ s2,d.α,d.θ)
      return a01*a1*(sum(seq1 .* exp.(a2 .* seq1) .* weights))
    end
  else
    if x>delta
      v = 1:Int(floor(l*d.α))
      w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v .* d.α .+ 1) .*
        sin.((pir*d.α) .* v) ./ gamma.(v .+ 1)
      return sum(w .* abs(x) .^ (-d.α .* v .- 1))/pi
    else
      nodes, weights = gausslegendre(m)
      nodes = nodes ./ 2 .+ .5
      weights = weights ./ 2
      nodes1 = nodes ./ (1 .- nodes)
      nodes2 = nodes ./ (1 .- nodes).^3
      pir = pi*d.ρ
      mat = abs(x) .* nodes1
      mat = pdf(UnilateralStable(d.α),mat)
      fC = nodes2 .* pdf(Cauchy(-cos(pir),sin(pir)),nodes1)
      return sum(fC .* mat .* weights)/d.ρ
    end
  end
end

function pdf(d::PositiveStable,X::AbstractArray{<:Real})
  if d.α == 1 && d.β == 1
    return convert(Array{Float64},X .== 1)
  end
  m = 10^3
  l = 15

  delta = 1
  x = vec(X)
  res = Float64[]
  #
  pir = pi*d.ρ
  pir1 = pi*(1 - d.ρ)
  v = 1:Int(floor(l*d.α))
  w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v ./ d.α .+ 1) .* sin.(pir .* v) ./
    ((d.ρ*pi) .* gamma.(v .+ 1))
  w1 = (-1) .^ ((v .- 1) .% 2) .* gamma.(v ./ d.α .+ 1) .* sin.(pir1 .* v) ./
    ((d.ρ*pi) .* gamma.(v .+ 1))
  #
  aux = gamma(1 + 1 / d.α)*sin(pi*d.ρ)/(d.ρ*pi)
  if d.α > 1
    a0 = d.α/((d.α-1)*2)
    raa = 1 /(d.α-1)
    a1 = abs.(x) .^ raa
    a2 = -a1 .^ d.α
    nodes, weights = gausslegendre(m)
    s1 = d.ρ
    s2 = 1 - d.ρ
    seq1 = auxV2(nodes .* s1 .+ s2,d.α,d.θ)
    weights = (a0*s1) .* weights
    for i = 1:length(x)
      if x[i] > 0
        if x[i] < delta
          push!(res,sum(w .* abs(x[i]) .^ (v .- 1)))
        else
          push!(res,a1[i]*(sum(seq1 .* exp.(a2[i] .* seq1) .* weights)))
        end
      elseif x[i] < 0
        push!(res,0.)
      else
        push!(res,aux)
      end
    end
  else
    nodes, weights = gausslegendre(m)
    nodes = nodes ./2 .+ .5
    weights = weights ./ (d.ρ * 2)
    nodes1 = nodes ./ (1 .- nodes)
    nodes2 = nodes ./ (1 .- nodes).^3
    pir = pi*d.ρ
    fC = nodes2*pdf(Cauchy(-cos(pir),sin(pir)),nodes1)
    mat = abs.(x)*transpose(nodes1)
    mat = pdf(UnilateralStable(d.α),mat)
    w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v .* d.α .+ 1) .*
      sin.((pir*d.α) .* v) ./ gamma.(v .+ 1)
    for i = 1:length(x)
      if x[i] > delta
        v = 1:Int(floor(l*d.α))
        push!(res,sum(w .* abs(x[i]) .^ (-d.α .* v .- 1))/pi)
      elseif x[i] > 0
        push!(res,sum(fC.*mat[i,:] .* weights))
      elseif x[i] < 0
        push!(res,0.)
      else
        push!(res,aux)
      end
    end
  end
  return reshape(res,size(X))
end

import Distributions.cdf
function cdf(d::PositiveStable,x::Real)
  if d.α == 1 && β == 1
    return x <= 1 ? 1. : 0.
  end
  if x <= 0
    return 0.
  end
  m = 10^3
  l = 15
  delta = 1.
  pir = pi*d.ρ
  pir1 = pi*(1-d.ρ)

  if d.α > 1
    raa = 1 / (d.α-1)
    a1 = x == 0 ? 0. : abs(x)^raa
    a2 = -a1 .^ d.α
    s1 = d.ρ
    s2 = 1 - d.ρ
    if x < delta
      v = 1:Int(floor(l*d.α))
      w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v ./ d.α .+ 1) .*
        sin.(pir .* v) ./ (v .* gamma.(v .+ 1))
      return sum(w .* abs(x) .^ v)/(d.ρ*pi)
    else
      nodes, weights = gausslegendre(m)
      seq1 = auxV2(nodes .*s1 .+ s2,d.α,d.θ)
      return 1 - s1*(sum(exp.(a2 .* seq1) .* weights))/(2 *d.ρ)
    end
  else
    if x>delta
      v = 1:Int(floor(l*d.α))
      w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v .* d.α .+ .1) .*
        sin.((pir*d.α) .* v) ./ (v .* gamma.(v .+ 1))
      return 1-sum(w .* abs(x) .^ (-d.α.*v))/(d.α*pi)
    else
      nodes, weights = gausslegendre(m)
      nodes = nodes ./ 2 .+ .5
      weights = weights ./ 2
      nodes1 = nodes ./ (1 .- nodes)
      nodes2 = 1 / (1 .- nodes) .^ 2
      pir = pi*d.ρ
      mat = abs(x)*nodes1
      mat = cdf(UnilateralStable(d.α),mat)
      fC = nodes2*pdf(Cauchy(-cos(pir),sin(pir)),nodes1)
      return sum(fC.*mat .* weights)/d.ρ
    end
  end
end

function cdf(d::PositiveStable,X::AbstractArray{<:Real})
  if d.α == 1 && β == 1
    return convert(Array{Float64},X .<= 1)
  end
  m = 10^3
  l = 15
  delta = 1.
  x = vec(X)
  res = Float64[]
  pir = pi*d.ρ
  pir1 = pi*(1 - d.ρ)
  v = 1:Int(floor(l*d.α))
  w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v ./ d.α .+ 1) .*
    sin.(pir .* v) ./ ((pi*d.ρ) .* v .* gamma.(v .+ 1))
  if d.α > 1
    raa = d.α/(d.α-1)
    a2 = -abs.(x) .^ raa
    nodes, weights = gausslegendre(m)
    s1 = d.ρ
    s2 = 1 - d.ρ
    seq1 = auxV2(nodes .* s1 .+ s2,d.α,d.θ)
    s1 = s1/(2 * d.ρ)
    for i = 1:length(x)
      if x[i] > 0
        if x[i] < delta
          push!(res,sum(w .* (x[i] .^ v)))
        else
          push!(res,1 - s1*(sum(exp.(a2[i] .* seq1) .* weights)))
        end
      else
        push!(res,0.)
      end
    end
  else
    nodes, weights = gausslegendre(m)
    nodes = nodes ./ 2 .+ .5
    weights = weights ./ (2 *d.ρ)
    nodes1 = nodes./(1-nodes)
    nodes2 = 1 ./ (1 .- nodes).^2
    pir = pi*d.ρ
    fC = nodes2*pdf(Cauchy(-cos(pir),sin(pir)),nodes1)
    mat = abs(x)*transpose(nodes1)
    mat = cdf(UnilateralStable(d.α),mat)
    w = (-1) .^ ((v .- 1) .% 2) .* gamma.(v .* d.α .+ 1) .*
      sin.((pir*d.α) .* v) ./ (v .* gamma.(v .+ 1))
    for i = 1:length(x)
      if x[i] > delta
        v = 1:Int(floor(l*d.α))
        push!(res, 1 - sum(w .* abs(x[i]) .^ (-d.α .* v))/(d.α*pi))
      elseif x[i] > 0
        push!(res,aux+sum(fC .* mat[i,:] .* weights))
      else
        push!(res,0.)
      end
    end
  end
  return reshape(res,size(X))
end

import Distributions.mgf
function mgf(d::PositiveStable,x::Real)
  l = 15
  if x == 0
    return 1.
  end
  if d.α == 2
    return 2 * exp(x^2/4.)*cdf(Normal(0,sqrt(2)),x/2)
  end
  if d.β == -1 && x >= -1
    v = 0:Int(floor(l*d.α))
    w = 1 / gamma(v+1)
    return sum(w .* x .^ v)
  end

  nodes, weights = gausslegendre(m)
  nodes = nodes ./ 2 .+ .5
  weights = weights/(2 *d.ρ)
  nodes2 = 1 ./ (1 .- nodes) .^ 2
  nodes = nodes ./ (1 .- nodes)
  pir = pi*d.ρ
  fC = nodes2 .* pdf(Cauchy(-cos(pir),sin(pir)),nodes)
  if d.β == -1
    mat = exp(-(abs(x) .* nodes) .^ d.α)
    return sum(fC.*mat .* weights)
  else
    if x > 0
      return Inf
    else
      mat = exp(-(abs(x) .* nodes) .^ d.α)
      return sum(fC.*mat .* weights)
    end
  end
end

function mgf(d::PositiveStable,X::AbstractArray{<:Real})
  x = vec(X)
  l = 15
  if d.α==2
    return 2 .* exp.(x.^2 ./ 4.) .* cdf(Normal(0,sqrt(2)),x ./ 2)
  end
  nodes, weights = gausslegendre(m)
  nodes = nodes ./ 2 .+ .5
  weights = weights ./ (2 * d.ρ)
  nodes2 = 1 ./ (1 .- nodes) .^ 2
  nodes = nodes ./ (1 .- nodes)
  pir = pi*d.ρ
  fC = nodes2*pdf(Cauchy(-cos(pir),sin(pir)),nodes)
  res = Float64[]
    mat = exp.(-(abs(x)*transpose(nodes)) .^ d.α)
  if d.β == -1
    v = 0:Int(floor(18*d.α))
    w = 1 ./ gamma.(v .+ 1)
    if x[i] == 0
      push!(res,1.)
    elseif x[i] >= -1
      push!(res,sum(w .* x[i] .^ v))
    else
      push!(res,sum(fC.*mat[i,:] .* weights))
    end
  else
    for i = 1:length(x)
      if x[i] > 0
        push!(res,Inf)
      elseif x[i] < 0
        push!(res,sum(fC.*mat[i,:] .* weights))
      else
        push!(res,1.)
      end
    end
  end
  return reshape(res,size(X))
end

import Distributions.mean
function mean(d::PositiveStable)
  if d.α <= 1
    return Inf
  end
  return sin(pi*d.ρ)/(d.α*d.ρ*sin(pi/d.α)*gamma(1 + 1 / d.α))
end

import Distributions.var
function var(d::PositiveStable)
  if d.α < 2 && d.β != -1
    return Inf
  elseif d.α == 1
    if abs(d.β) != 1
      return Inf
    else
      return 0.
    end
  else
    return 2 / gamma(1 + 2 / d.α) - 1 / gamma(1 + 1 / d.α) ^ 2
  end
end

function mellin(d::PositiveStable,x::Complex)
  if (real(x) >= d.α && (d.α <= 1 || (d.α < 2 && d.β != -1))) || real(x) <= -1
    return Inf
  end
  if (d.α > 1 && d.β == -1) || d.α == 2
    return gamma(1 + x) / gamma(1 + x / d.α)
  end
  return (sin(pi * d.ρ * x) * gamma(1 + x)) /
    (d.α * d.ρ * sin(pi * x / d.α) * gamma(1 + x / d.α))
end

function mellin(d::PositiveStable,X::AbstractArray{<:Complex})
  if (d.α > 1 && d.β == -1) || d.α == 2
    return gamma.(1 .+ X) ./ gamma.(1 .+ X ./ d.α)
  end

  res = Complex{Float64}[]
  for x = X
    push!(res,
      (real(x) >= d.α && (d.α <= 1 || (d.α < 2 && d.β != -1))) || real(x) <= -1 ?
      Inf : (sin(pi * d.ρ * x) * gamma(1 + x)) / (d.α * d.ρ * sin(pi * x / d.α) * gamma(1 + x / d.α))
    )
  end
  return reshape(res,size(X))
end

function mellin(d::PositiveStable,x::Real)
  if (real(x) >= d.α && (d.α <= 1 || (d.α < 2 && d.β != -1))) || real(x) <= -1
    return Inf
  end
  if (d.α > 1 && d.β == -1) || d.α == 2
    return gamma(1 + x) / gamma(1 + x / d.α)
  end
  return (sin(pi * d.ρ * x) * gamma(1 + x)) /
  (d.α * d.ρ * sin(pi * x / d.α) * gamma(1 + x / d.α))
end

function mellin(d::PositiveStable,X::AbstractArray{<:Real})
  if (d.α > 1 && d.β == -1) || d.α == 2
    return gamma.(1 .+ X) ./ gamma.(1 .+ X ./ d.α)
  end

  res = Float64[]
  for x = X
    push!(res,
      (real(x) >= d.α && (d.α <= 1 || (d.α < 2 && d.β != -1))) || real(x) <= -1 ?
      Inf : (sin(pi * d.ρ * x) * gamma(1 + x)) / (d.α * d.ρ * sin(pi * x / d.α) * gamma(1 + x / d.α))
    )
  end
  return reshape(res,size(X))
end

export Stable, rand, minimum, maximum, insupport, pdf, cdf, mgf, mean, mellin, params


text:
 using Documenter, CellularAutomata

makedocs(sitename="CellularAutomata.jl",
pages = [
    "CellularAutomata.jl" => "index.md",
    "Examples" => [
        "One dimensional CA" => "onedim/onedimensionca.md"
        "Two dimensional CA" => "twodim/twodimensionca.md"
        ],
    "User Guide" => Any[
        "One Dimensional CA" => "user/onedimensionca.md"
        "Two Dimensial CA" => "user/twodimensionca.md"

        ]
    ])


text:
using Printf
using DelimitedFiles
using LinearAlgebra
using ForwardDiff: jacobian

const SN = V.NUM            # num of state variables
const PN = 1                # num of parameters
const VN = SN + PN          # num of variables
const MN = SN               # dim of Newton's method

const MC = 100000           # maximum of counts
const IVAL = 1e-2           # first variation
const RATE = 1e-3           # variation rate
const NEPS = 1e-12          # eps of Newton's method

# matrix transformation (large diagonal elements move to upper)
function pivoting!(s::Matrix{Float64}, pivot::Int)
    v0::Vector{Float64} = zeros(MN+1)
    v1::Vector{Float64} = zeros(MN+1)
    possess::Int = 0
    max_element::Float64 = 0.0

    for i in pivot:size(s,1)
        current_element = abs(s[i, pivot])
        if max_element <= current_element
            max_element = current_element
            possess = i
        end
    end

    for j in 1:size(s,2)
        v0[j] = s[possess, j]
        v1[j] = s[pivot, j]
    end

    for j in 1:size(s,2)
        s[possess, j] = v1[j]
        s[pivot, j] = v0[j]
    end
end

# Gaussian elimination (row reduction)
function gaussian_elimination!(s::Matrix{Float64}, e::Vector{Float64})
    for i in 1:MN
        pivoting!(s, i)
    end
    # forward
    for k in 1:size(s,1)
        w = (s[k, k] != 0.0) ? 1.0 / s[k, k] : 1.0
        for j in k:size(s,2)
            s[k, j] *= w
            for i in k:size(s,1)
                s[i, j] -= s[i, k] * s[k, j]
            end
        end
    end
    # backward
    for i in size(s,1):-1:1
        sum = 0.0
        for j in i:length(e)
            sum += s[i, j] * e[j]
        end
        e[i] = s[i, end] - sum
    end
end

# Newton's method
function newtons_method!(
        x::Vector{Float64},
        real_part::Vector{Float64},
        imaginary_part::Vector{Float64},
        fix_num::Int,
        p::Vector{Float64},
        successful::Bool)
    u::Vector{Float64} = zeros(SN)
    vx::Vector{Float64} = zeros(MN)
    s::Matrix{Float64} = zeros(MN, MN+1)

    for i in eachindex(x)
        if fix_num == i
            for j in eachindex(vx)
                idx = i + j
                if idx > length(x)
                    idx -= length(x)
                end
                vx[j] = x[idx]
            end
            break
        else
            continue
        end
    end

    # initial error
    e::Vector{Float64} = zeros(MN)
    error::Float64 = 1.0

    while error > NEPS
        for i in 1:VN
            if fix_num == i
                idx_param = VN - i
                p[BP] = (idx_param == 0) ? x[fix_num] : vx[idx_param]
                for j in eachindex(u)
                    idx = j - i
                    if idx == 0
                        u[j] = x[fix_num]
                    elseif idx < 0
                        u[j] = vx[VN+idx]
                    else
                        u[j] = vx[idx]
                    end
                end
                break
            else
                continue
            end
        end

        # initialization
        dFdx::Matrix{Float64} = jacobian(diffeq, u)
        dFdp::Vector{Float64} = get_derivatives(u, p)

        F::Vector{Float64} = diffeq(u)

        eigenvalues::Array{Complex{Float64}, 1} = eigvals(dFdx)
        for (i, eigenvalue) in enumerate(eigenvalues)
            real_part[i] = real(eigenvalue)
            imaginary_part[i] = imag(eigenvalue)
        end

        # s = [dF-F]
        for i in 1:VN
            if fix_num == i
                for k in 1:SN
                    for j in 1:SN
                        idx = i + j
                        if idx == VN
                            s[k, j] = dFdp[k]
                        elseif idx > VN
                            s[k, j] = dFdx[k, idx-VN]
                        else
                            s[k, j] = dFdx[k, idx]
                        end
                    end
                    s[k, VN] = -F[k]
                end
                break
            else
                continue
            end
        end

        gaussian_elimination!(s, e)

        # update error
        error = 0.0
        @inbounds for i in eachindex(e)
            vx[i] += e[i]
            error += e[i] * e[i]
        end
        error = sqrt(error)
        if isnan(error) || isinf(error)
            successful = false
            break
        end
    end

    for i in eachindex(x)
        if fix_num == i
            for j in eachindex(vx)
                idx = i + j
                if idx > length(x)
                    idx -= length(x)
                end
                x[idx] = vx[j]
            end
            break
        else
            continue
        end
    end
end

function new_curve!(p::Vector{Float64}; direction::Bool=false)
    count::Int = 1
    x::Vector{Float64} = zeros(VN)
    dx::Vector{Float64} = zeros(VN)

    real_part::Vector{Float64} = zeros(SN)
    imaginary_part::Vector{Float64} = zeros(SN)

    # file
    if !isdir("./Data")
        mkdir("./Data")
    else
        files::Vector{String} = readdir("./Data")
        for file in files
            rm("./Data/$file")
        end
    end
    FOUT1 = open("./Data/fp.dat", "w")  # file for fixed point
    FOUT2 = open("./Data/ev.dat", "w")  # file for eigenvalues

    # initial condition
    x[1:SN] = get_steady_state(p)
    x[end] = p[BP]  # x-axis

    # initial fixed
    fix_val::Float64 = x[end]
    fix_num::Int = VN
    x[fix_num] = fix_val

    # first Newton's method
    successful::Bool = true
    newtons_method!(x, real_part, imaginary_part, fix_num, p, successful)

    write(FOUT1, @sprintf("%d\t", count))
    for i in eachindex(x)
        write(FOUT1, @sprintf("%10.8e\t", x[i]))
    end
    write(FOUT1, @sprintf("%d\n", fix_num))
    write(FOUT2, @sprintf("%d\t", count))
    for i in 1:SN
        write(
            FOUT2, @sprintf(
                "%10.8e\t%10.8e\t", real_part[i], imaginary_part[i]
            )
        )
    end
    write(FOUT2, @sprintf("%10.8e\t%d\n", p[BP], fix_num))
    count += 1

    # keep optimums
    px::Vector{Float64} = copy(x)

    # variation
    fix_val += ifelse(direction, +IVAL, -IVAL)

    # same fixed variable
    x[fix_num] = fix_val

    while count <= MC && successful
        newtons_method!(x, real_part, imaginary_part, fix_num, p, successful)

        # maximum variation
        for (i, prev) in enumerate(px)
            @inbounds dx[i] = x[i] - prev
        end
        sum::Float64 = 0.0
        for i in eachindex(dx)
            @inbounds sum += dx[i] * dx[i]
        end
        ave::Float64 = sqrt(sum)
        for i in eachindex(dx)
            @inbounds dx[i] /= ave
        end
        px = copy(x)
        for (i, diff) in enumerate(dx)
            @inbounds x[i] += abs(RATE) * diff
        end

        # fix variable with maximum variation
        fix_num = 1
        for i in 2:length(dx)
            if abs(dx[fix_num]) < abs(dx[i])
                fix_num = i
            end
        end

        # Stop calc.
        if x[end] <= 0.0
            successful = false
        end

        write(FOUT1, @sprintf("%d\t", count))
        for i in eachindex(x)
            write(FOUT1, @sprintf("%10.8e\t", x[i]))
        end
        write(FOUT1, @sprintf("%d\n", fix_num))
        write(FOUT2, @sprintf("%d\t", count))
        for i in 1:SN
            write(
                FOUT2, @sprintf(
                    "%10.8e\t%10.8e\t", real_part[i], imaginary_part[i]
                )
            )
        end
        write(FOUT2, @sprintf("%10.8e\t%d\n", p[BP], fix_num))
        count += 1
    end

    close(FOUT1)
    close(FOUT2)
end

function get_bistable_regime(ev::Matrix{Float64})
    br::Vector{Int} = []
    for i in 1:size(ev, 1)
        if maximum(ev[i, [2j for j in 1:SN]]) > 0.0
            push!(br, i)
        end
    end
    return br
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
export to_string_dict

to_string_dict(dict::Dict) = Dict{String,Any}(string(k)=>v for (k,v) in dict)
TOML.print(io::IO,dict::Dict{Symbol,V}) where {V} = TOML.print(io,to_string_dict(dict))

export to_symbol_dict
to_symbol_dict(el) = el
to_symbol_dict(dict::Dict) = Dict{Symbol,Any}(Symbol(k)=>to_symbol_dict(v) for (k,v) in dict)
to_symbol_dict(v::Vector) = map(to_symbol_dict, v)

read_assignment(toml_dict::Dict) = toml_dict["assignment"]
read_assignment(io) = read_assignment(TOML.parsefile(io))

function TOML.parse(mat::SparseMatrixCSC{T,Int64}) where {T}
    I,J,_ = findnz(mat)
    return Dict("edge_list" => map(idx->[idx...], zip(I,J)),
        "size_i" => size(mat,1),
        "size_j" => size(mat,2)
        )
end
function read_sparse_matrix(toml_dict::Dict)
    edge_list   = toml_dict["edge_list"]
    size_i      = toml_dict["size_i"]
    size_j      = toml_dict["size_j"]
    mat = spzeros(Int,size_i,size_j)
    for e in edge_list
        mat[e[1],e[2]] = 1
    end
    mat
end
read_sparse_matrix(io) = read_sparse_matrix(TOML.parsefile(io))

export read_abstract_id

parse_abstract_id(id::A,dict=Dict{String,Any}()) where {A<:AbstractID} = merge!(dict,Dict("type"=>string(A),"id"=>get_id(id)))
TOML.parse(id::A) where {A<:AbstractID} = parse_abstract_id(id)
function read_abstract_id(dict)
    type_string = dict["type"]
    for t in [ActionID(),ObjectID(),BotID{DeliveryBot}(),BotID{CleanUpBot}()]
        if type_string == string(typeof(t))
            return typeof(t)(dict["id"])
        end
    end
    return nothing
end
function parse_predicate(a::A,dict=Dict{String,Any}()) where {A<:AbstractPlanningPredicate}
    dict["type"] = string(A)
    dict["id"] = TOML.parse()
    return dict
end
# function TOML.parse(sched::OperatingSchedule)
#
# end

export
    remap_object_id,
    remap_object_ids!

# utilities for remapping object ids
remap_object_id(x,args...) = x
remap_object_id(id::ObjectID,max_obj_id)    = ObjectID(get_id(id) + max_obj_id)
# remap_object_id(spec::PathSpec,max_obj_id)  = PathSpec(spec,object_id=spec.object_id + max_obj_id)
remap_object_id(node::OBJECT_AT,args...)    = OBJECT_AT(remap_object_id(get_object_id(node),args...), get_initial_location_id(node))
remap_object_id(node::A,args...) where {A<:Union{COLLECT,DEPOSIT}} = A(o=remap_object_id(get_object_id(node),args...),r=get_robot_id(node),x=get_location_id(node))
remap_object_id(node::A,args...) where {A<:CARRY} = A(o=remap_object_id(get_object_id(node),args...),r=node.r,x1=node.x1,x2=node.x2)
remap_object_id(node::A,args...) where {A<:TEAM_ACTION} = A(instructions=map(i->remap_object_id(i,args...),node.instructions),shape=node.shape)
remap_object_id(node::ScheduleNode,args...) = ScheduleNode(
    remap_object_id(node.id,args...),
    remap_object_id(node.node,args...),
    remap_object_id(node.spec,args...)
    )
function remap_object_ids!(vec::Vector,args...)
    map!(i->remap_object_id(i,args...),vec,vec)
end
function remap_object_ids!(dict::Dict,args...)
    for k in collect(keys(dict))
        val = remap_object_id(dict[k],args...)
        delete!(dict, k)
        dict[remap_object_id(k,args...)] = val
    end
    dict
end
function remap_object_ids!(node::Operation,args...)
    remap_object_ids!(preconditions(node),args...)
    remap_object_ids!(postconditions(node),args...)
    return node
end
remap_object_id(node::Operation,args...)    = remap_object_ids!(deepcopy(node),args...)
function remap_object_ids!(sched::OperatingSchedule,args...)
    remap_object_ids!(get_vtx_ids(sched),args...)
    remap_object_ids!(get_nodes(sched),args...)
    remap_object_ids!(get_vtx_map(sched),args...)
    @assert sanity_check(sched," after remap_object_ids!()")
    sched
end
function remap_object_ids!(new_schedule::OperatingSchedule,old_schedule::OperatingSchedule)
    max_obj_id = 0
    for id in get_vtx_ids(old_schedule)
        if typeof(id) <: ObjectID
            max_obj_id = max(get_id(id),max_obj_id)
        end
    end
    remap_object_ids!(new_schedule,max_obj_id)
end

export get_valid_robot_ids

"""
    get_valid_robot_ids(sched::OperatingSchedule,node_id::AbstractID)

Returns vector of all robot ids associated with the schedule node referenced by
node_id.
"""
function get_valid_robot_ids(sched::OperatingSchedule,node_id::A,v=get_vtx(sched,node_id)) where {A<:Union{ActionID,BotID}}
    ids = Vector{BotID}()
    node = get_node(sched,node_id).node
    if matches_template(TEAM_ACTION,node)
        for n in sub_nodes(node)
            r = get_robot_id(n)
            if get_id(r) != -1
                push!(ids,r)
            end
        end
    else
        push!(ids,get_default_robot_id(node))
    end
    return filter(CRCBS.is_valid,ids)
end
function get_valid_robot_ids(sched::OperatingSchedule,node_id::A,v=get_vtx(sched,node_id)) where {A<:Union{ObjectID,OperationID}}
    return Vector{BotID}()
end
get_valid_robot_ids(s::OperatingSchedule,v::Int) = get_valid_robot_ids(s,get_vtx_id(s,v),v)
get_valid_robot_ids(node) = BotID[]
get_valid_robot_ids(node::Union{BOT_AT,AbstractRobotAction}) = [get_robot_id(node)]
get_valid_robot_ids(node::TEAM_ACTION) = filter(CRCBS.is_valid,map(n->get_robot_id(n),sub_nodes(node)))

export robot_tip_map

"""
    robot_tip_map(sched::OperatingSchedule)

Returns a `Dict{RobotID,AbstractID}` mapping `RobotID` to the terminal node of
the `sched` corresponding to the robot's last assigned task.
"""
function robot_tip_map(sched::OperatingSchedule,vtxs=get_all_terminal_nodes(sched))
    robot_tips = Dict{BotID,AbstractID}()
    for v in vtxs
        node_id = get_vtx_id(sched,v)
        for robot_id in get_valid_robot_ids(sched,node_id,v)
            if get_id(robot_id) != -1
                @assert !haskey(robot_tips,robot_id)
                robot_tips[robot_id] = node_id
            end
        end
    end
    # @assert length(robot_tips) == length(get_robot_ICs(sched)) "length(robot_tips) == $(length(robot_tips)), but should be $(length(get_robot_ICs(sched)))"
    robot_tips
end

export construct_task_graphs_problem
"""
    `construct_task_graphs_problem`
"""
function construct_task_graphs_problem(
        project_spec::ProjectSpec,
        r0::Vector{Int},
        s0::Vector{Int},
        sF::Vector{Int},
        env;
        Δt_collect=zeros(length(s0)),
        Δt_deposit=zeros(length(sF)),
        cost_function=SumOfMakeSpans(),
        task_shapes=map(o->(1,1),s0),
        shape_dict=Dict{Int,Dict{Tuple{Int,Int},Vector{Int}}}(s=>Dict{Tuple{Int,Int},Vector{Int}}() for s in vcat(s0,sF))
        )
    N = length(r0)
    M = length(s0)
    for j in 1:M
        @assert haskey(shape_dict, s0[j]) "shape_dict has no key for s0[$j] = $(s0[j])"
        @assert haskey(shape_dict, sF[j]) "shape_dict has no key for sF[$j] = $(sF[j])"
        @assert length(task_shapes) >= j "task_shapes has no key for j = $j"
    end
    robot_ICs = [ROBOT_AT(r,x) for (r,x) in enumerate(r0)] # initial robot conditions
    object_ICs = [OBJECT_AT(j, get(shape_dict[x], s, x), s) for (j,(x,s)) in enumerate(zip(s0,task_shapes))] # initial object conditions
    object_FCs = [OBJECT_AT(j, get(shape_dict[x], s, x), s) for (j,(x,s)) in enumerate(zip(sF,task_shapes))] # initial object conditions
    new_project_spec = ProjectSpec(object_ICs,object_FCs)
    for op in get_operations(project_spec)
        add_operation!(
            new_project_spec,
            construct_operation(
                new_project_spec,
                op.station_id,
                collect(keys(preconditions(op))),
                collect(keys(postconditions(op))),
                duration(op)
            )
        )
    end

    problem_spec = ProblemSpec(
        D=env,
        cost_function=cost_function,
        Δt_collect = Dict{ObjectID,Int}(
            get_object_id(n)=>t for (n,t) in zip(object_ICs,Δt_collect)),
        Δt_deposit = Dict{ObjectID,Int}(
            get_object_id(n)=>t for (n,t) in zip(object_ICs,Δt_deposit)),
        )

    sched = construct_partial_project_schedule(new_project_spec,problem_spec,robot_ICs)
    return sched, problem_spec
end
function construct_task_graphs_problem(
    def::SimpleProblemDef,
    env::GridFactoryEnvironment;
    task_shapes = def.shapes,
    kwargs...,
)
    construct_task_graphs_problem(
        def.project_spec,
        def.r0,
        def.s0,
        def.sF,
        env;
        task_shapes = task_shapes,
        shape_dict = env.expanded_zones,
        kwargs...,
    )
end

export
    construct_random_project_spec,
    get_random_problem_instantiation,
    construct_random_task_graphs_problem

"""
    construct_random_project_spec(M::Int;max_children=1)

    Inputs:
        `M` - number of objects involved in the operation
        `max_parents` - determines the max number of inputs to any operation
        `depth_bias` ∈ [0,1] - hyperparameter for tuning depth.
            If `depth_bias` == 1.0, the project_spec graph will always be depth
            balanced (all paths through the tree will be of the same length).
            For `depth_bias` == 0.0, the graph will be as "strung out" as
            possible.
"""
function construct_random_project_spec(M::Int,object_ICs::Vector{OBJECT_AT},object_FCs::Vector{OBJECT_AT};
        max_parents=1,
        depth_bias=1.0,
        Δt_min=0,
        Δt_max=0,
        zone_map=Dict{Int,Vector{Int}}(), # maps vtx to vector of vtxs for collaborative tasks
    )
    project_spec = ProjectSpec(object_ICs,object_FCs)
    @assert length(object_ICs) == M
    object_map = map(get_object_id, initial_conditions_vector(project_spec))
    # fill with random operations going backwards
    i = M-1
    frontier = PriorityQueue{Int,Int}([M=>1])
    while i > 0
        depth = 1
        while true
            if (rand() > depth_bias) && (depth < length(frontier))
                depth += 1
            else
                break
            end
        end
        pairs = Vector{Pair{Int,Int}}()
        # pair = Pair{Int,Int}(0,0)
        for d in 1:depth
            push!(pairs, dequeue_pair!(frontier))
        end
        for p in pairs[1:end-1]
            enqueue!(frontier,p)
        end
        output_idx = pairs[end].first
        station_id = output_idx
        input_idxs = collect(max(1,1+i-rand(1:max_parents)):i)
        i = i - length(input_idxs)
        Δt=rand(Δt_min:Δt_max)
        input_ids = map(idx->object_map[idx], input_idxs)
        output_ids = map(idx->object_map[idx], [output_idx])
        add_operation!(project_spec,construct_operation(project_spec, station_id, input_ids, output_ids, Δt))
        for idx in input_idxs
            enqueue!(frontier, idx=>M-i)
        end
    end
    Δt=0
    final_idx = object_map[M]
    add_operation!(project_spec,construct_operation(project_spec, -1, [final_idx], [], Δt))
    project_spec
end
# function construct_random_project_spec(M::Int,s0::Vector{Int},sF::Vector{Int};
function construct_random_project_spec(M::Int,s0::Vector{V},sF::Vector{V};kwargs...) where {V<:Union{Int,Vector{Int}}}
    object_ICs = Vector{OBJECT_AT}([OBJECT_AT(o,s) for (o,s) in enumerate(s0)])
    object_FCs = Vector{OBJECT_AT}([OBJECT_AT(o,s) for (o,s) in enumerate(sF)])
    construct_random_project_spec(M,object_ICs,object_FCs;
        kwargs...)
end

export
    choose_random_object_sizes,
    choose_and_shuffle_object_sizes

"""
    Tool for randomly selecting how many robots (and in what configuration)
        should deliver each task.
"""
function choose_random_object_sizes(M,probs::Dict{Tuple{Int,Int},Float64})
    k = sort(collect(keys(probs)))
    v = cumsum(map(key->probs[key], k))
    v = v / sum(collect(values(probs)))
    sizes = Vector{Tuple{Int,Int}}()
    for j in 1:M
        n = rand()
        i = 1
        while n > v[i]
            i += 1
        end
        push!(sizes,k[i])
    end
    sizes
end
function choose_and_shuffle_object_sizes(M,probs::Dict{Tuple{Int,Int},Float64})
    k = sort(collect(keys(probs)),by=i->i[1]*i[2])
    v = map(key->probs[key], k)
    v = M * v / sum(v)
    v_rounded = Int.(round.(v))
    while sum(v_rounded) < sum(v)
        idx = sort(collect(1:length(v)),by=i->abs(v[i]-v_rounded[i]))[end]
        v_rounded[idx] = v_rounded[idx] + 1
    end
    while sum(v_rounded) > sum(v)
        idx = sort(collect(1:length(v)),by=i->abs(v[i]-v_rounded[i]))[end]
        v_rounded[idx] = v_rounded[idx] - 1
    end
    sizes = Vector{Tuple{Int,Int}}()
    for (idx,n) in enumerate(v_rounded)
        for i in 1:n
            push!(sizes,k[idx])
        end
    end
    shuffle(sizes)
end
function choose_random_object_sizes(M,probs::Dict{Int,Float64}=Dict(1=>1.0,2=>0.0,4=>0.0),choices=[(1,1),(2,1),(1,2),(2,2)])
    k = sort(collect(keys(probs)))
    size_probs = Dict{Tuple{Int,Int},Float64}(s=>probs[s[1]*s[2]] for s in choices)
    for ksize in k
        counts = sum(map(s->s[1]*s[2]==ksize, choices))
        for s in choices
            size_probs[s] = size_probs[s] / counts
        end
    end
    choose_and_shuffle_object_sizes(M,size_probs)
end

"""
    `get_random_problem_instantiation`

    Args:
    - `N`: number of robots
    - `M`: number of delivery tasks
    - `robot_zones`: list of possible start locations for robots
    - `pickup_zones`: list of possible start locations for objects
    - `dropoff_zones`: list of possible destinations for objects
"""
function get_random_problem_instantiation(N::Int,M::Int,pickup_zones,dropoff_zones,robot_zones)
    ##### Random Problem Initialization #####
    r0 = robot_zones[sortperm(rand(length(robot_zones)))][1:N]
    s0 = pickup_zones[sortperm(rand(length(pickup_zones)))][1:M]
    sF = dropoff_zones[sortperm(rand(length(dropoff_zones)))][1:M]
    return r0,s0,sF
end

"""
    `construct_randomd_task_graphs_problem`
"""
function construct_random_task_graphs_problem(N::Int,M::Int,
    pickup_vtxs::Vector{Int},dropoff_vtxs::Vector{Int},free_vtxs::Vector{Int},dist_matrix,
    Δt_collect::Vector{Float64}=zeros(M),
    Δt_deposit::Vector{Float64}=zeros(M)
    )
    # select subset of pickup, dropoff and free locations to instantiate objects and robots
    r0,s0,sF        = get_random_problem_instantiation(N,M,pickup_vtxs,dropoff_vtxs,free_vtxs)
    project_spec    = construct_random_project_spec(M,s0,sF;max_parents=3,depth_bias=1.0,Δt_min=0,Δt_max=0)

    construct_task_graphs_problem(project_spec,r0,s0,sF,dist_matrix;Δt_collect=Δt_collect,Δt_deposit=Δt_deposit)
end

################################################################################
################################### Rendering ##################################
################################################################################

export
    get_display_metagraph

function get_display_metagraph(project_schedule::OperatingSchedule;
    verbose=true,
    f=(v,p)->title_string(p,verbose),
    object_color="orange",
    robot_color="lime",
    clean_up_bot_color="yellow",
    action_color="cyan",
    operation_color="red",
    remove_leaf_robots=false
    )
    graph = MetaDiGraph(deepcopy(project_schedule.graph))
    for (id,pred) in get_object_ICs(project_schedule)
        v = get_vtx(project_schedule, get_object_id(pred))
        set_prop!(graph, v, :vtype, :object_ic)
        set_prop!(graph, v, :text, f(v,pred))
        set_prop!(graph, v, :color, object_color)
        set_prop!(graph, v, :vtx_id, v)
    end
    for (id,op) in get_operations(project_schedule)
        v = get_vtx(project_schedule, OperationID(id))
        set_prop!(graph, v, :vtype, :operation)
        set_prop!(graph, v, :text, f(v,op))
        set_prop!(graph, v, :color, operation_color)
        set_prop!(graph, v, :vtx_id, v)
    end
    for (id,a) in get_actions(project_schedule)
        v = get_vtx(project_schedule, ActionID(id))
        set_prop!(graph, v, :vtype, :action)
        set_prop!(graph, v, :text, f(v,a))
        set_prop!(graph, v, :color, action_color)
        set_prop!(graph, v, :vtx_id, v)
    end
    for (id,pred) in get_robot_ICs(project_schedule)
        v = get_vtx(project_schedule, get_robot_id(pred))
        set_prop!(graph, v, :vtype, :robot_ic)
        set_prop!(graph, v, :text, f(v,pred))
        if isa(id,CleanUpBotID)
            set_prop!(graph, v, :color, clean_up_bot_color)
        else
            set_prop!(graph, v, :color, robot_color)
        end
        set_prop!(graph, v, :vtx_id, v)
    end
    if remove_leaf_robots == true
        for v in reverse(vertices(graph))
            # if get_prop(graph,v,:vtype) ==:robot_ic
            if typeof(get_vtx_id(project_schedule,v)) <: RobotID
                if length(outneighbors(graph,v)) == 0
                    rem_vertex!(graph,v)
                end
            end
        end
    end
    graph
end

################################################################################
######################### COMPARISONS AND VERIFICATION #########################
################################################################################
function Base.:(==)(o1::OBJECT_AT,o2::OBJECT_AT)
    try
        @assert(o1.o == o2.o)
        @assert(o1.x == o2.x)
    catch e
        if isa(e,AssertionError)
            # println(e.msg)
        else
            throw(e)
        end
        return false
    end
    return true
end
function Base.:(==)(op1::Operation,op2::Operation)
    try
        @assert(preconditions(op1) == preconditions(op2))
        @assert(postconditions(op1) == postconditions(op2))
        @assert(op1.Δt == op2.Δt)
        @assert(op1.station_id == op2.station_id)
    catch e
        if isa(e,AssertionError)
            # println(e.msg)
        else
            throw(e)
        end
        return false
    end
    return true
end
function Base.:(==)(spec1::ProjectSpec,spec2::ProjectSpec)
    try
        @assert(spec1.initial_conditions == spec2.initial_conditions)
        @assert(spec1.final_conditions == spec2.final_conditions)
    catch e
        if isa(e,AssertionError)
            # println(e.msg)
        else
            throw(e)
        end
        return false
    end
    return true
end

export
    get_object_paths

function get_object_paths(solution,sched)
    robot_paths = convert_to_vertex_lists(solution)
    tF = maximum(map(length, robot_paths))
    object_paths = Vector{Vector{Int}}()
    object_intervals = Vector{Vector{Int}}()
    object_ids = Int[]
    path_idxs = Int[]
    for v in vertices(sched.graph)
        node = get_node_from_id(sched,get_vtx_id(sched,v))
        if isa(node, Union{CARRY,TEAM_CARRY})
            if isa(node, CARRY)
                object_id = get_object_id(node)
                agent_id_list = [get_id(get_robot_id(node))]
                s0_list = [get_id(get_initial_location_id(node))]
                sF_list = [get_id(get_destination_location_id(node))]
            else
                object_id = get_object_id(node.instructions[1])
                agent_id_list = map(n->get_id(get_robot_id(n)), node.instructions)
                s0_list = map(n->get_id(get_initial_location_id(n)), node.instructions)
                sF_list = map(n->get_id(get_destination_location_id(n)), node.instructions)
            end
            object_vtx = get_vtx(sched,object_id)
            for (idx,(agent_id,s0,sF)) in enumerate(zip(agent_id_list,s0_list,sF_list))
                if get_id(agent_id) != -1
                    push!(object_paths,[
                        map(t->s0,0:get_t0(sched,v)-1)...,
                        map(t->robot_paths[agent_id][t],min(get_t0(sched,v)+1,tF):min(get_tF(sched,v)+1,tF,length(robot_paths[agent_id])))...,
                        map(t->sF,min(get_tF(sched,v)+1,tF):tF)...
                    ])
                    push!(object_intervals,[get_t0(sched,object_vtx),get_tF(sched,v)+1])
                    push!(object_ids, get_id(object_id))
                    push!(path_idxs, idx)
                end
            end
        end
    end
    object_paths, object_intervals, object_ids, path_idxs
end
# function get_object_paths(solution,env::SearchEnv)
#     get_object_paths(solution,get_schedule(env))
# end

export
    fill_object_path_dicts!,
    convert_to_path_vectors

function fill_object_path_dicts!(solution,project_schedule,cache,
        object_path_dict = Dict{Int,Vector{Vector{Int}}}(),
        object_interval_dict = Dict{Int,Vector{Int}}()
    )
    object_paths, object_intervals, object_ids, path_idxs = get_object_paths(solution,project_schedule)
    for (path,interval,id,idx) in zip(object_paths, object_intervals, object_ids, path_idxs)
        if !haskey(object_path_dict,id)
            object_path_dict[id] = Vector{Vector{Int}}()
        end
        if length(object_path_dict[id]) >= idx
            object_path_dict[id][idx] = path
        else
            push!(object_path_dict[id], path)
        end
        object_interval_dict[id] = interval
    end
    object_path_dict, object_interval_dict
end

function convert_to_path_vectors(object_path_dict, object_interval_dict)
    object_paths = Vector{Vector{Int}}()
    object_intervals = Vector{Vector{Int}}()
    for id in sort(collect(keys(object_path_dict)))
        paths = object_path_dict[id]
        for path in paths
            push!(object_intervals, object_interval_dict[id])
            push!(object_paths, path)
        end
    end
    object_paths, object_intervals
end


text:
function loadBPP(instance::Symbol)
    return load(instance, :BPP)
end

function loadBPP(file_name::String)
    return load(file_name, :BPP)
end

function loadCSP(instance::Symbol)
    data = load(instance, :BPP)
    return data === nothing ? data : toCSP(data)
end

function loadCSP(file_name::String)
    return load(file_name, :CSP)
end

function load(instance::Symbol, type::Symbol)
    file_name = joinpath(data_path, string(instance) * ".zip")
    if !isfile(file_name)
        println("File $(string(instance)) not found!")
        return nothing
    end
    name = splitext(basename(file_name))[1]

    file = ZipFile.Reader(file_name)
    values = parse.(Int64, split(read(file.files[1], String)))
    return load(values, name, type)
end

function load(file_name::String, type::Symbol)
    if !isfile(file_name)
        println("File $file_name not found!")
        return nothing
    end

    name = splitext(basename(file_name))[1]
    values = parse.(Int64, split(read(file_name, String)))
    return load(values, name, type)
end

function load(values::Array{Int64}, name::String, type::Symbol)
    capacity = values[2]
    lb, ub = loadBounds(name)
    if type == :BPP
        weights = values[3:end]
        return BPPData(name, capacity, weights, lb, ub)
    else
        weights = values[3:2:end]
        demands = values[4:2:end]
        return CSPData(name, capacity, weights, demands, lb, ub)
    end
end