text:
function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarAffineFunction{Float64},
   set::MOI.LessThan{Float64}
)
   check_inbounds(model, func)
   push!(model.linear_le_constraints, ConstraintInfo(func, set))
   n = length(model.linear_le_constraints)
   return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarAffineFunction{Float64},
   set::MOI.GreaterThan{Float64}
)
   check_inbounds(model, func)
   push!(model.linear_ge_constraints, ConstraintInfo(func, set))
   n = length(model.linear_ge_constraints)
   return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarAffineFunction{Float64},
   set::MOI.EqualTo{Float64}
)
   check_inbounds(model, func)
   push!(model.linear_eq_constraints, ConstraintInfo(func, set))
   n = length(model.linear_eq_constraints)
   return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarQuadraticFunction{Float64},
   set::MOI.LessThan{Float64}
)
   check_inbounds(model, func)
   push!(model.quadratic_le_constraints, ConstraintInfo(func, set))
   n = length(model.quadratic_le_constraints)
   return MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarQuadraticFunction{Float64},
   set::MOI.GreaterThan{Float64}
)
   check_inbounds(model, func)
   push!(model.quadratic_ge_constraints, ConstraintInfo(func, set))
   n = length(model.quadratic_ge_constraints)
   return MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.ScalarQuadraticFunction{Float64},
   set::MOI.EqualTo{Float64}
)
   check_inbounds(model, func)
   push!(model.quadratic_eq_constraints, ConstraintInfo(func, set))
   n = length(model.quadratic_eq_constraints)
   return MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}(n)
end

function MOI.add_constraint(
   model::Optimizer,
   func::MOI.VectorAffineFunction{Float64},
   set::MOI.Complements
)
   check_inbounds(model, func)
   push!(model.complementarity_constraints, ConstraintInfo(func, set))
   n = length(model.complementarity_constraints)
   return MathOptInterface.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Complements}(n)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
   return length(model.linear_le_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
   return length(model.linear_ge_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
   return length(model.linear_eq_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}
)
   return length(model.quadratic_le_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}
)
   return length(model.quadratic_ge_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}
)
   return length(model.quadratic_eq_constraints)
end

function MOI.get(
   model::Optimizer,
   ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64}, MOI.Complements}
)
   return length(model.complementarity_constraints)
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.SingleVariable},
   value
)
   error("Constraint names for variable bound constraints not supported.")
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.linear_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.linear_le_constraints[ci.value].name = value
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.linear_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.linear_ge_constraints[ci.value].name = value
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.linear_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.linear_eq_constraints[ci.value].name = value
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.quadratic_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.quadratic_le_constraints[ci.value].name = value
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.quadratic_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.quadratic_ge_constraints[ci.value].name = value
   return
end

function MOI.set(
   model::Optimizer,
   attr::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}},
   value
)
   if ! (1 <= ci.value <= length(model.quadratic_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   model.quadratic_eq_constraints[ci.value].name = value
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.linear_le_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.linear_ge_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.linear_eq_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.quadratic_le_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.quadratic_ge_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.quadratic_eq_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.LessThan(model.linear_le_constraints[ci.value].set.upper)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.GreaterThan(model.linear_ge_constraints[ci.value].set.lower)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
   if ! (1 <= ci.value <= length(model.linear_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.EqualTo(model.linear_eq_constraints[ci.value].set.value)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.LessThan(model.quadratic_le_constraints[ci.value].set.upper)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.GreaterThan(model.quadratic_ge_constraints[ci.value].set.lower)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintSet,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}
)
   if ! (1 <= ci.value <= length(model.quadratic_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return MOI.EqualTo(model.quadratic_eq_constraints[ci.value].set.value)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ConstraintName,
   ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, MOI.Complements}
)
   if ! (1 <= ci.value <= length(model.complementarity_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.complementarity_constraints[ci.value].name
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.linear_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_linear_le(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.linear_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_linear_ge(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.linear_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_linear_eq(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.quadratic_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_quadratic_le(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.quadratic_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_quadratic_ge(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintPrimal,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if ! (1 <= ci.value <= length(model.quadratic_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_quadratic_eq(model))"].level[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.linear_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return -model.sol.equ["eq$(ci.value + offset_linear_le(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.linear_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_linear_ge(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.linear_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return model.sol.equ["eq$(ci.value + offset_linear_eq(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.quadratic_le_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return -model.sol.equ["eq$(ci.value + offset_quadratic_le(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.quadratic_ge_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return -model.sol.equ["eq$(ci.value + offset_quadratic_ge(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.ConstraintDual,
   ci::MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}
)
   MOI.check_result_index_bounds(model, attr)
   if !(1 <= ci.value <= length(model.quadratic_eq_constraints))
      error("Invalid constraint index ", ci.value)
   end
   return -model.sol.equ["eq$(ci.value + offset_quadratic_eq(model))"].dual[1]
end

function MOI.get(
   model::Optimizer,
   attr::MOI.NLPBlockDual
)
   MOI.check_result_index_bounds(model, attr)
   values = zeros(model.m_nonlin)
   for i = offset_nonlin(model)+1:offset_nonlin(model)+model.m_nonlin
      values[i] = -model.sol.equ["eq$i"].dual[1]
   end
   return values
end


text:
function CC.typeinf(interp::XInterpreter, frame::InferenceState)
    ret = @invoke typeinf(interp::AbstractInterpreter, frame::InferenceState)

    return ret
end

function CC.finish!(interp::XInterpreter, caller::InferenceResult)
    ret = @invoke finish!(interp::AbstractInterpreter, caller::InferenceResult)

    return ret
end

function CC._typeinf(interp::XInterpreter, frame::InferenceState)
    ret = @invoke _typeinf(interp::AbstractInterpreter, frame::InferenceState)

    return ret
end

function CC.maybe_compress_codeinfo(interp::XInterpreter, linfo::MethodInstance, ci::CodeInfo)
    ret = @invoke maybe_compress_codeinfo(interp::AbstractInterpreter, linfo::MethodInstance, ci::CodeInfo)

    return ret
end

function CC.transform_result_for_cache(interp::XInterpreter, linfo::MethodInstance,
                                    @nospecialize(inferred_result))
    ret = @invoke transform_result_for_cache(interp::AbstractInterpreter, linfo::MethodInstance,
                                             inferred_result)

    return ret
end

function CC.cache_result!(interp::XInterpreter, result::InferenceResult, valid_worlds::WorldRange)
    ret = @invoke cache_result!(interp::AbstractInterpreter, result::InferenceResult, valid_worlds::WorldRange)

    return ret
end

function CC.finish(me::InferenceState, interp::XInterpreter)
    ret = @invoke finish(me::InferenceState, interp::AbstractInterpreter)

    return ret
end

function CC.is_same_frame(interp::XInterpreter, linfo::MethodInstance, frame::InferenceState)
    ret = @invoke is_same_frame(interp::AbstractInterpreter, linfo::MethodInstance, frame::InferenceState)

    return ret
end

function CC.resolve_call_cycle!(interp::XInterpreter, linfo::MethodInstance, parent::InferenceState)
    ret = @invoke resolve_call_cycle!(interp::AbstractInterpreter, linfo::MethodInstance, parent::InferenceState)

    return ret
end

function CC.typeinf_edge(interp::XInterpreter, method::Method, @nospecialize(atypes), sparams::SimpleVector, caller::InferenceState)
    ret = @invoke typeinf_edge(interp::AbstractInterpreter, method::Method, atypes, sparams::SimpleVector, caller::InferenceState)

    return ret
end

ext: jl


text:
const mpi_init_hooks = Any[]

macro mpi_handle(def, mpiname=nothing, extrafields...)
    name = def isa Expr ? def.args[1] : def
    _name = Symbol(:_, name)
    if mpiname == nothing
        mpiname = Symbol(:MPI_, name)
    end
    quote
        Base.@__doc__ mutable struct $(esc(def))
            val::$mpiname
            $(extrafields...)

            $(esc(name))(val::$mpiname, $(extrafields...)) = new(val, $(extrafields...))
        end

        # const initializer
        function $(esc(_name))(ival::Cint, extraargs...)
            x = $name(reinterpret($mpiname, sizeof($mpiname) == 4 ?  Int32(-1) : Int64(-1)), extraargs...)
            push!(mpi_init_hooks, () -> x.val = $mpiname(ival))
            return x
        end

        function Base.:(==)(a::$name, b::$name)
            a.val == b.val
        end

        # allows us to pass XX objects directly into MPI_XX ccall signatures
        function Base.cconvert(::Type{$mpiname}, obj::$name)
            obj.val
        end
        # allows us to pass XX objects directly into Ptr{MPI_XX} ccall signatures
        function Base.unsafe_convert(::Type{Ptr{$mpiname}}, obj::$name)
            convert(Ptr{$mpiname}, pointer_from_objref(obj))
        end
    end
end

"""
    MPI.free(obj)

Free the MPI object handle `obj`. This is typically used as the finalizer, and so need not be called directly unless otherwise noted.
"""
function free
end


text:
#Period types
value{P<:Period}(x::P) = x.value

# The default constructors for Periods work well in almost all cases
# P(x) = new((convert(Int64,x))
# The following definitions are for Period-specific safety
for p in (:Year,:Month,:Week,:Day,:Hour,:Minute,:Second,:Millisecond)
    # This ensures that we can't convert between Periods
    @eval $p(x::Period) = throw(ArgumentError("Can't convert $(typeof(x)) to $($p)"))
    # Unless the Periods are the same type
    @eval $p(x::$p) = x
    # Convenience method for show()
    @eval _units(x::$p) = " " * lowercase(string($p)) * (abs(value(x)) == 1 ? "" : "s")
    # periodisless
    @eval periodisless(x::$p,y::$p) = value(x) < value(y)
    # String parsing (mainly for IO code)
    @eval $p(x::String) = $p(parseint(Int64,x))
    # Period accessors
    @eval $p(x::TimeType) = $p($(symbol(lowercase(string(p))))(x))
end
# Now we're safe to define Period-Number conversions
# Anything an Int64 can convert to, a Period can convert to
Base.convert{T<:Number}(::Type{T},x::Period) = convert(T,value(x))
Base.convert{T<:Period}(::Type{T},x::Real) = T(int64(x))

#Print/show/traits
Base.string{P<:Period}(x::P) = string(value(x),_units(x))
Base.show(io::IO,x::Period) = print(io,string(x))
Base.zero{P<:Period}(::Union(Type{P},P)) = P(0)
Base.one{P<:Period}(::Union(Type{P},P)) = P(1)
Base.typemin{P<:Period}(::Type{P}) = P(typemin(Int64))
Base.typemax{P<:Period}(::Type{P}) = P(typemax(Int64))

# Default values (as used by TimeTypes)
default{T<:DatePeriod}(p::Union(T,Type{T})) = one(p)
default{T<:TimePeriod}(p::Union(T,Type{T})) = zero(p)

(-){P<:Period}(x::P) = P(-value(x))
Base.isless{P<:Period}(x::P,y::P) = isless(value(x),value(y))
Base.isless{R<:Real}(x::Period,y::R) = isless(int64(y),value(x))
Base.isless{R<:Real}(y::R,x::Period) = isless(int64(y),value(x))
=={P<:Period}(x::P,y::P) = ===(value(x),value(y))
=={R<:Real}(x::Period,y::R) = ==(int64(y),value(x))
=={R<:Real}(y::R,x::Period) = ==(int64(y),value(x))

Base.isless(x::Period,y::Period) = throw(ArgumentError("Can't compare $(typeof(x)) and $(typeof(y))"))
==(x::Period,y::Period) = throw(ArgumentError("Can't compare $(typeof(x)) and $(typeof(y))"))

#Period Arithmetic:
import Base.div, Base.mod, Base.gcd, Base.lcm
let vec_ops = [:.+,:.-,:.*,:.%,:div]
    for op in [:+,:-,:*,:%,:mod,:gcd,:lcm,vec_ops]
        @eval begin
        #Period-Period
        ($op){P<:Period}(x::P,y::P) = P(($op)(value(x),value(y)))
        #Period-Real
        ($op){P<:Period}(x::P,y::Real) = P(($op)(value(x),convert(Int64,y)))
        ($op){P<:Period}(x::Real,y::P) = P(($op)(convert(Int64,x),value(y)))
        end
        #Vectorized
        if op in vec_ops
            @eval begin
            ($op){P<:Period}(x::AbstractArray{P}, y::P) = reshape([$op(i,y) for i in x], size(x))
            ($op){P<:Period}(y::P, x::AbstractArray{P}) = reshape([$op(i,y) for i in x], size(x))
            end
        end
    end
end

# intfuncs
Base.gcdx{T<:Period}(a::T,b::T) = ((g,x,y)=gcdx(value(a),value(b)); return T(g),x,y)
Base.abs{T<:Period}(a::T) = T(abs(value(a)))

periodisless(::Period,::Year)        = true
periodisless(::Period,::Month)       = true
periodisless(::Year,::Month)         = false
periodisless(::Period,::Week)        = true
periodisless(::Year,::Week)          = false
periodisless(::Month,::Week)         = false
periodisless(::Period,::Day)         = true
periodisless(::Year,::Day)           = false
periodisless(::Month,::Day)          = false
periodisless(::Week,::Day)           = false
periodisless(::Period,::Hour)        = false
periodisless(::Minute,::Hour)        = true
periodisless(::Second,::Hour)        = true
periodisless(::Millisecond,::Hour)   = true
periodisless(::Period,::Minute)      = false
periodisless(::Second,::Minute)      = true
periodisless(::Millisecond,::Minute) = true
periodisless(::Period,::Second)      = false
periodisless(::Millisecond,::Second) = true
periodisless(::Period,::Millisecond) = false

# Stores multiple periods in greatest to least order by type, not values
type CompoundPeriod
    periods::Array{Period,1}
end
function Base.string(x::CompoundPeriod)
    s = ""
    for p in x.periods
        s *= ", " * string(p)
    end
    return s[3:end]
end
Base.show(io::IO,x::CompoundPeriod) = print(io,string(x))
# E.g. Year(1) + Day(1)
(+)(x::Period,y::Period) = CompoundPeriod(sort!(Period[x,y],rev=true,lt=periodisless))
(+)(x::CompoundPeriod,y::Period) = (sort!(push!(x.periods,y) ,rev=true,lt=periodisless); return x)
# E.g. Year(1) - Month(1)
(-)(x::Period,y::Period) = CompoundPeriod(sort!(Period[x,-y],rev=true,lt=periodisless))
(-)(x::CompoundPeriod,y::Period) = (sort!(push!(x.periods,-y),rev=true,lt=periodisless); return x)

# Capture TimeType+-Period methods
(+)(a::TimeType,b::Period,c::Period) = (+)(a,b+c)
(-)(a::TimeType,b::Period,c::Period) = (-)(a,b-c)
(+)(a::TimeType,b::Period,c::Period,d::Period...) = (+)((+)(a,b+c),d...)
(-)(a::TimeType,b::Period,c::Period,d::Period...) = (-)((-)(a,b-c),d...)

function (+)(x::TimeType,y::CompoundPeriod)
    for p in y.periods
        x += p
    end
    return x
end
(+)(x::CompoundPeriod,y::TimeType) = y + x

# Convert fixed value Periods to # of milliseconds
typealias FixedPeriod Union(Week,Day,Hour,Minute,Second,Millisecond)

toms(c::Millisecond) = value(c)
toms(c::Second)      = 1000*value(c)
toms(c::Minute)      = 60000*value(c)
toms(c::Hour)        = 3600000*value(c)
toms(c::Day)         = 86400000*value(c)
toms(c::Week)        = 604800000*value(c)
toms(c::Month)       = 86400000.0*30.436875*value(c)
toms(c::Year)        = 86400000.0*365.2425*value(c)

Millisecond{T<:FixedPeriod}(c::T) = Millisecond(toms(c))

days(c::Millisecond) = div(value(c),86400000)
days(c::Second)      = div(value(c),86400)
days(c::Minute)      = div(value(c),1440)
days(c::Hour)        = div(value(c),24)
days(c::Day)         = value(c)
days(c::Week)        = 7*value(c)
days(c::Year)        = 365.2425*value(c)
days(c::Month)       = 30.436875*value(c)

Day{T<:FixedPeriod}(c::T) = Day(days(c))


text:
#   Author: V. Vitale
#   Feb 2022
#   1 site and 2 sites DMRG

include("mps.jl")
include("mpo.jl")
include("contractions.jl")
using TensorOperations
using KrylovKit
using LinearAlgebra

function one_site_dmrg!(psi::MPS,
                        W::MPO,
                        sweeps::Int64)

    right_normalize!(psi)

    d = dims(psi)[2][2]
    chi = dims(psi)[2][1]
    L = construct_L(psi, W)
    R = construct_R(psi, W)

    #global psi,L,R
    Energy=0
    for sweep in 1:Int(sweeps/2)
        for i in 1:psi.N-1
            Energy,psi.data[i],psi.data[i+1] = optimize_one_site(psi.data[i], psi.data[i+1], W.data[i], L[i], R[i], "right", chi)
            L[i+1] = contract_from_left(L[i], psi.data[i], W.data[i])
        end

        for i in psi.N:-1:2
            Energy,psi.data[i-1],psi.data[i] = optimize_one_site( psi.data[i-1],psi.data[i], W.data[i], L[i], R[i], "left", chi)
            R[i-1] = contract_from_right(R[i], psi.data[i], W.data[i])
        end
    end
    println("Done! Energy= ",real(Energy),"; Variance: ",real(psi*(W*(W*psi))-(psi*(W*psi))^2))
end

function optimize_one_site(A::Array, B::Array, W::Array, E::Array, F::Array, dir::String, chi::Int )
    function H_lin(v)
        ##     +--A--+
        ##     |  |  |
        ##     L--W--R
        ##     |  |  |
        #println("E ",size(E))
        #println("MPS ",size(v))
        #println("F ",size(F))
        #println("W ",size(W))
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * W[ 2 3 ; 1 -2 ]
        return temp
    end

    if (dir == "right" )
        sA=size(A)

        en,V,info = eigsolve( H_lin , A ,  1  , :SR ; ishermitian = true  )
        en = en[1]
        V = V[1]
        sV = size(V)

        V=reshape(V,(sV[1]*sV[2],sV[3]))#s i j -> i*s, j
        A,S,V = svd(V,full=false)
        V=V'
        A = reshape( A , ( sA[1], sA[2], :) )
        #"ij,jl,slk->sik"
        S=diagm(S)
        @tensor B[:] := S[-1,1 ] * V[ 1,2 ] * B[2,-2,-3]
    else
        sB=size(B)
        en,V,info = eigsolve( H_lin , B ,  1  , :SR ; ishermitian = true  )
        en = en[1]
        V = V[1]

        sV=size(V)
        V=reshape(V,(sV[1],sV[2]*sV[3]))  #s i j -> i s j -> i, s*j

        U,S,B = svd(V,full=false)
        B=B'
        B = reshape(B,(:,sB[2],sB[3]))
        #"sij,jk,kl->sil"
        S=diagm(S)
        @tensor A[:] := A[-1,-2,2] * U[ 2,3 ] * S[ 3,-3 ]
    end

    return en, A, B

end

function two_sites_dmrg!(psi::MPS,
                        W::MPO,
                        sweeps::Int)

    right_normalize!(psi)

    d = dims(psi)[2][2]
    L = construct_L(psi, W)
    R = construct_R(psi, W)

    #global psi,L,R
    Energy=0
    for sweep in 1:Int(sweeps)
        println("Sweep: ",sweep)
        for i in 1:psi.N-1
            Energy,psi.data[i],psi.data[i+1] = two_sites_swipe_right( psi.data[i],
                                                                        psi.data[i+1],
                                                                        W.data[i],
                                                                        W.data[i+1],
                                                                        L[i],
                                                                        R[i+1])
            L[i+1] = contract_from_left(L[i], psi.data[i], W.data[i])
        end

        for i in psi.N:-1:2
            Energy,psi.data[i-1],psi.data[i] = two_sites_swipe_left(  psi.data[i-1],
                                                                        psi.data[i],
                                                                        W.data[i-1],
                                                                        W.data[i],
                                                                        L[i-1],
                                                                        R[i])
            R[i-1] = contract_from_right(R[i], psi.data[i], W.data[i])
        end

    end
    println("Done! Energy= ",real(Energy),"; Variance: ",real(psi*(W*(W*psi))-(psi*(W*psi))^2))
end

function two_sites_swipe_right(AL::Array, AR::Array, WL::Array, WR::Array, E::Array, F::Array)
    tol=1e-15

    sAL = size(AL)
    sAR = size(AR)
    @tensor A[:] := AL[-1,-2,1]*AR[1,-3,-4]
    A = reshape(A,(sAL[1],sAL[2]*sAR[2],sAR[3]))
    sWL = size(WL)
    sWR = size(WR)

    @tensor M[:] := WL[-1,1,-3,-5]*WR[1,-2,-4,-6]
    M = reshape(M,(sWL[1],sWR[2],sWL[3]*sWR[3],sWL[4]*sWR[4]))
    function H_lin(v)
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * M[ 2 3 ; 1 -2 ]
        return temp
    end

    en,V,info = eigsolve( H_lin , A ,  1  , :SR ; ishermitian = true  )
    en = en[1]
    V = V[1]

    V=reshape(V,(sAL[1]*sAL[2],sAR[2]*sAR[3]))
    U,S,V = svd(V,full=false)
    V=V'
    S=S/norm(S)
    indices = findall(1 .-cumsum(S.^2) .< tol)
    if length(indices)>0
        chi = indices[1]+1
    else
        chi = size(S)[1]
    end

    if size(S)[1] > chi
        U = U[:,1:chi]
        S = S[1:chi]
        V =  V[1:chi,:]
    end

    S /= norm(S)

    AL = reshape( U , ( sAL[1], sAL[2], :) )
    #"ij,jl,slk->sik"
    S=diagm(S)
    @tensor AR[:] := S[-1,1 ] * V[ 1,-2 ]
    AR=reshape(AR,(:,sAR[2],sAR[3]))
    return en, AL, AR
end

function two_sites_swipe_left(AL::Array, AR::Array, WL::Array, WR::Array, E::Array, F::Array)
    tol=1e-15

    sAL = size(AL)
    sAR = size(AR)

    @tensor A[:] :=AL[-1,-2,1]*AR[1,-3,-4]
    A = reshape(A,(sAL[1],sAL[2]*sAR[2],sAR[3]))

    sWL = size(WL)
    sWR = size(WR)

    @tensor M[:] := WL[-1,1,-3,-5]*WR[1,-2,-4,-6]
    M = reshape(M,(sWL[1],sWR[2],sWL[3]*sWR[3],sWL[4]*sWR[4]))

    function H_lin(v)
            @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
            @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
            @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * M[ 2 3 ; 1 -2 ]
        return temp
    end

    en,V,info = eigsolve( H_lin , A ,  1  , :SR ; ishermitian = true  )
    en = en[1]
    V = V[1]
    sV=size(V)
    V=reshape(V,(sAL[1]*sAL[2],sAR[2]*sAR[3]))

    U,S,V = svd(V,full=false)
    V=V'

    S /=norm(S)
    indices = findall(1 .-cumsum(S.^2) .< tol)
    if length(indices)>0
        chi = indices[1]+1
    else
        chi = size(S)[1]
    end
    if size(S)[1] > chi
        U = U[:,1:chi]
        S = S[1:chi]
        V =  V[1:chi,:]
    end

    S /= norm(S)

    AR = reshape(V,(:,sAR[2],sAR[3]))

    S=diagm(S)
    @tensor AL[:] :=  U[ -1,1 ] * S[ 1,-2 ]
    AL = reshape( AL , ( sAL[1], sAL[2], :) )
    return en, AL, AR
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
include("utilities/misc.jl")
include("utilities/group_actions.jl")
include("utilities/sorting_points.jl")
include("utilities/polynomials.jl")


text:
using WRDSMerger
using Documenter

Documenter.makedocs(
    modules = [WRDSMerger],
    sitename = "WRDSMerger.jl",
)


text:
## Various solver API's used in the past.  These functions are due to be standardized, and obsolete code / functions removed.

export solveTree!, solveGraph!
export fetchCliqHistoryAll!

## ==============================================================================================
## Launch the tasks/treads for cliques
## ==============================================================================================
"""
    $SIGNATURES

Start tasks (@async or Threads.@spawn threads if multithread=true) to solve the factor graph on the tree.
"""
function taskSolveTree!(dfg::AbstractDFG,
                          treel::AbstractBayesTree,
                          timeout::Union{Nothing, <:Real}=nothing;
                          oldtree::AbstractBayesTree=BayesTree(),
                          drawtree::Bool=false,
                          verbose::Bool=false,
                          verbosefid=stdout,
                          limititers::Int=-1,
                          limititercliqs::Vector{Pair{Symbol, Int}}=Pair{Symbol, Int}[],
                          downsolve::Bool=false,
                          incremental::Bool=false,
                          multithread::Bool=false,
                          skipcliqids::Vector{Symbol}=Symbol[],
                          recordcliqs::Vector{Symbol}=Symbol[],
                          delaycliqs::Vector{Symbol}=Symbol[],
                          smtasks=Task[],
                          algorithm::Symbol=:default,
                          solveKey::Symbol=algorithm )
  #
  # revert DOWNSOLVED status to INITIALIZED in preparation for new upsolve
  resetTreeCliquesForUpSolve!(treel)

  drawtree ? drawTree(treel, show=false, filepath=joinLogPath(dfg,"bt.dot")) : nothing

  cliqHistories = Dict{Int,Vector{CSMHistoryTuple}}()

  resize!(smtasks, getNumCliqs(treel))

  approx_iters = getNumCliqs(treel)*24
  solve_progressbar = verbose ? nothing : ProgressUnknown("Solve Progress: approx max $approx_iters, at iter")

  # queue all the tasks/threads
  if !isTreeSolved(treel, skipinitialized=true)
    @sync begin
      monitortask = monitorCSMs(treel, smtasks)
      # duplicate int i into async (important for concurrency)
      for i in 1:getNumCliqs(treel) # TODO, this might not always work?
        scsym = getCliqFrontalVarIds(getClique(treel, i))
        if length(intersect(scsym, skipcliqids)) == 0

          limthiscsm = filter(x -> (x[1] in scsym), limititercliqs)
          limiter = 0<length(limthiscsm) ? limthiscsm[1][2] : limititers

          if multithread
            smtasks[i] = Threads.@spawn tryCliqStateMachineSolve!(dfg, treel, i, timeout; solveKey=solveKey, algorithm=algorithm, oldtree=oldtree, verbose=verbose, verbosefid=verbosefid, drawtree=drawtree, limititers=limititers, downsolve=downsolve, incremental=incremental, delaycliqs=delaycliqs, recordcliqs=recordcliqs, solve_progressbar=solve_progressbar)
          else
            smtasks[i] = @async tryCliqStateMachineSolve!(dfg, treel, i, timeout; solveKey=solveKey, algorithm=algorithm, oldtree=oldtree, verbose=verbose, verbosefid=verbosefid, drawtree=drawtree, limititers=limiter, downsolve=downsolve, incremental=incremental, delaycliqs=delaycliqs, recordcliqs=recordcliqs, solve_progressbar=solve_progressbar)
          end
        end # if
      end # for
    end # sync
  end # if

  # if record cliques is in use, else skip computational delay
  0 == length(recordcliqs) ? nothing : fetchCliqHistoryAll!(smtasks, cliqHistories)

  !isnothing(solve_progressbar) && finish!(solve_progressbar)

  return smtasks, cliqHistories
end

function tryCliqStateMachineSolve!( dfg::G,
                                    treel::AbstractBayesTree,
                                    cliqKey::Union{Int, CliqueId},
                                    timeout::Union{Nothing, <:Real}=nothing;
                                    oldtree::AbstractBayesTree=BayesTree(),
                                    verbose::Bool=false,
                                    verbosefid=stdout,
                                    drawtree::Bool=false,
                                    limititers::Int=-1,
                                    downsolve::Bool=false,
                                    incremental::Bool=false,
                                    delaycliqs::Vector{Symbol}=Symbol[],
                                    recordcliqs::Vector{Symbol}=Symbol[],
                                    solve_progressbar=nothing,
                                    algorithm::Symbol=:default,
                                    solveKey::Symbol=algorithm ) where G <: AbstractDFG
  #
  clst = :na
  cliq = getClique(treel, cliqKey) #treel.cliques[cliqKey]
  syms = getCliqFrontalVarIds(cliq) # ids =

  oldcliq = attemptTreeSimilarClique(oldtree, getCliqueData(cliq))
  oldcliqdata = getCliqueData(oldcliq)

  opts = getSolverParams(dfg)
  # Base.rm(joinpath(opts.logpath,"logs/cliq$i"), recursive=true, force=true)
  mkpath(joinpath(opts.logpath,"logs/cliq$(cliq.id)/"))
  logger = SimpleLogger(open(joinpath(opts.logpath,"logs/cliq$(cliq.id)/log.txt"), "w+")) # NullLogger()
  # global_logger(logger)
  history = Vector{CSMHistoryTuple}()
  recordthiscliq = length(intersect(recordcliqs,syms)) > 0
  delaythiscliq = length(intersect(delaycliqs,syms)) > 0
  try
    history = initStartCliqStateMachine!( dfg, treel, cliq, timeout;
                                          oldcliqdata=oldcliqdata,
                                          drawtree=drawtree, verbose=verbose, verbosefid=verbosefid,
                                          limititers=limititers, downsolve=downsolve,
                                          recordhistory=recordthiscliq, incremental=incremental,
                                          delay=delaythiscliq, logger=logger, solve_progressbar=solve_progressbar,
                                          algorithm=algorithm, solveKey=solveKey )
    #
    # cliqHistories[cliqKey] = history
    if length(history) >= limititers && limititers != -1
      # @warn "writing logs/cliq$(cliq.id)/csm.txt"
      # @save "/tmp/cliqHistories/cliq$(cliq.id).jld2" history
      fid = open(joinpath(opts.logpath,"logs/cliq$(cliq.id)/csm.txt"), "w")
      printCliqHistorySummary(fid, history)
      close(fid)
    end
    flush(logger.stream)
    close(logger.stream)
    # clst = getCliqueStatus(cliq)
    # clst = cliqInitSolveUp!(dfg, treel, cliq, drawtree=drawtree, limititers=limititers )
  catch err
    bt = catch_backtrace()
    println()
    showerror(stderr, err, bt)
    # @warn "writing /tmp/caesar/logs/cliq$(cliq.id)/*.txt"
    fid = open(joinpath(opts.logpath,"logs/cliq$(cliq.id)/stacktrace.txt"), "w")
    showerror(fid, err, bt)
    close(fid)
    fid = open(joinpath(opts.logpath,"logs/cliq$(cliq.id)_stacktrace.txt"), "w")
    showerror(fid, err, bt)
    close(fid)
    # @save "/tmp/cliqHistories/$(cliq.label).jld2" history
    fid = open(joinpath(opts.logpath,"logs/cliq$(cliq.id)/csm.txt"), "w")
    printCliqHistorySummary(fid, history)
    close(fid)
    fid = open(joinpath(opts.logpath,"logs/cliq$(cliq.id)_csm.txt"), "w")
    printCliqHistorySummary(fid, history)
    close(fid)
    flush(logger.stream)
    close(logger.stream)
    rethrow()
  end
  # if !(clst in [UPSOLVED; DOWNSOLVED; MARGINALIZED])
  #   error("Clique $(cliq.id), initInferTreeUp! -- cliqInitSolveUp! did not arrive at the desired solution statu: $clst")
  # end
  return history
end

"""
    $SIGNATURES

Standalone state machine solution for a single clique.

Related:

initInferTreeUp!
"""
function solveCliqWithStateMachine!(dfg::G,
                                    tree::AbstractBayesTree,
                                    frontal::Symbol;
                                    iters::Int=200,
                                    downsolve::Bool=true,
                                    recordhistory::Bool=false,
                                    verbose::Bool=false,
                                    nextfnc::Function=canCliqMargRecycle_StateMachine,
                                    prevcsmc::Union{Nothing,CliqStateMachineContainer}=nothing) where G <: AbstractDFG
  #
  cliq = getClique(tree, frontal)

  children = getChildren(tree, cliq)#Graphs.out_neighbors(cliq, tree.bt)

  prnt = getParent(tree, cliq)

  destType = (G <: InMemoryDFGTypes) ? G : InMemDFGType

  csmc = isa(prevcsmc, Nothing) ? CliqStateMachineContainer(dfg, initfg(destType, solverParams=getSolverParams(dfg)), tree, cliq, prnt, children, false, true, true, downsolve, false, getSolverParams(dfg)) : prevcsmc
  statemachine = StateMachine{CliqStateMachineContainer}(next=nextfnc, name="cliq$(cliq.id)")
  while statemachine(csmc, verbose=verbose, iterlimit=iters, recordhistory=recordhistory); end
  statemachine, csmc
end

## ==============================================================================================
# Prepare CSM (based on FSM) entry points
## ==============================================================================================

"""
    $SIGNATURES

Fetch solver history from clique state machines that have completed their async Tasks and store in the `hist::Dict{Int,Tuple}` dictionary.
"""
function fetchCliqHistoryAll!(smt::Vector{Task},
                              hist::Dict{Int,Vector{CSMHistoryTuple}}=Dict{Int,Vector{CSMHistoryTuple}}() )
  #
  for i in 1:length(smt)
    sm = smt[i]
    # only fetch states that have completed processing
    if sm.state == :done
      haskey(hist, i) ? @warn("overwriting existing history key $i") : nothing
      hist[i] = fetch(sm)
    elseif !isnothing(sm.storage) && haskey(sm.storage, :statemachine)
      hist[i] = CSMHistoryTuple.(sm.storage[:statemachine].history)
    end
  end
  hist
end

## ==============================================================================================
# Nominal user interface to the solver
## ==============================================================================================

"""
    $SIGNATURES

Perform inference over the Bayes tree according to `opt::SolverParams` and keyword arguments.

Notes
- Aliased with `solveGraph!`
- Variety of options, including fixed-lag solving -- see `getSolverParams(fg)` for details.
  - See online Documentation for more details: https://juliarobotics.org/Caesar.jl/latest/
- Latest result always stored in `solvekey=:default`.
- Experimental `storeOld::Bool=true` will duplicate the current result as supersolve `:default_k`.
  - Based on `solvable==1` assumption.
- `limititercliqs` allows user to limit the number of iterations a specific CSM does.
- keywords `verbose` and `verbosefid::IOStream` can be used together to to send output to file or default `stdout`.
- keyword `recordcliqs=[:x0; :x7...]` identifies by frontals which cliques to record CSM steps.
  - See [`repeatCSMStep!`](@ref), [`printCSMHistoryLogical`](@ref), [`printCSMHistorySequential`](@ref)

DevNotes
- TODO Change keyword arguments to new @parameter `SolverOptions` type.

Example
```julia
# pass in old `tree` to enable compute recycling -- see online Documentation for more details
tree = solveTree!(fg [,tree])
```

Related

`solveGraph!`, [`solveCliqUp!`](@ref), [`solveCliqDown!`](@ref), [`buildTreeReset!`](@ref), [`repeatCSMStep`](@ref), [`printCSMHistoryLogical`](@ref)
"""
function solveTree!(dfgl::AbstractDFG,
                    oldtree::AbstractBayesTree=BayesTree();
                    timeout::Union{Nothing, <:Real}=nothing,
                    storeOld::Bool=false,
                    verbose::Bool=false,
                    verbosefid=stdout,
                    delaycliqs::Vector{Symbol}=Symbol[],
                    recordcliqs::Vector{Symbol}=Symbol[],
                    limititercliqs::Vector{Pair{Symbol, Int}}=Pair{Symbol, Int}[],
                    injectDelayBefore::Union{Nothing,Vector{<:Pair{Int,<:Pair{<:Function, <:Real}}}}=nothing,
                    skipcliqids::Vector{Symbol}=Symbol[],
                    eliminationOrder::Union{Nothing, Vector{Symbol}}=nothing,
                    eliminationConstraints::Vector{Symbol}=Symbol[],
                    smtasks::Vector{Task}=Task[],
                    dotreedraw = Int[1;],
                    runtaskmonitor::Bool=true,
                    algorithm::Symbol=:default,
                    solveKey::Symbol=:default,
                    multithread::Bool=false  )
  #
  # workaround in case isolated variables occur
  ensureSolvable!(dfgl)
  opt = getSolverParams(dfgl)

  # showtree should force drawtree
  opt.showtree && !opt.drawtree ? @info("Since .showtree=true, also bumping .drawtree=true") : nothing
  opt.drawtree |= opt.showtree

  # depcrecation
  # update worker pool incase there are more or less
  setWorkerPool!()
  if opt.multiproc && nprocs() == 1
    @info "Setting `.multiproc=false` since `Distributed.nprocs() == 1`"
    opt.multiproc = false
  end

  if opt.graphinit
    @info "Ensure variables are all initialized (graphinit)"
    initAll!(dfgl, solveKey)
    if algorithm==:parametric
      @warn "Parametric is using default graphinit (and ignoring solveKey)"
      initParametricFrom!(dfgl)
    end
  end
  # construct tree
  @info "Solving over the Bayes (Junction) tree."

  hist = Dict{Int, Vector{CSMHistoryTuple}}()

  if opt.isfixedlag
      @info "Quasi fixed-lag is enabled (a feature currently in testing, and ignoring solveKey)!"
      fifoFreeze!(dfgl)
  end

  # perhaps duplicate current value
  if storeOld || opt.dbg
    ss = listSupersolves(dfgl) .|> string
    ss_ = ss[occursin.(r"default_",ss)] .|> x->x[9:end]
    filter!(x->occursin(r"^\d+$",x), ss_)  # ss_ = ss_[occursin.(r"^\d$",ss_)]
    allk = parse.(Int,ss_)
    nextk = length(allk) == 0 ? 0 : maximum(allk)+1
    newKey = Symbol(:default_, nextk)
    deepcopySupersolve!(dfgl, newKey, :default, solvable=1)
    # foreach(x->updateVariableSolverData!(dfgl, x, getSolverData(getVariable(dfgl,x), :default), newKey, true, Symbol[]), ls(dfgl, solvable=1))
    @info "storeOld=true, previous :default deepcopied into $newKey for solvable==1 variables."
  end

  orderMethod = 0 < length(eliminationConstraints) ? :ccolamd : :qr

  # current incremental solver builds a new tree and matches against old tree for recycling.
  tree = buildTreeReset!(dfgl, eliminationOrder, drawpdf=false, show=opt.showtree,ensureSolvable=false,filepath=joinpath(opt.logpath,"bt.pdf"), eliminationConstraints=eliminationConstraints, ordering=orderMethod)

  # setAllSolveFlags!(tree, false)

  initTreeMessageChannels!(tree)

  # if desired, drawtree in a loop
  treetask, _dotreedraw = drawTreeAsyncLoop(tree, opt; dotreedraw = dotreedraw)

  @info "Do tree based init-ference"
  algorithm != :parametric ? nothing : @error("Under development, do not use, see #539")
  !storeOld ? nothing : @error("parametric storeOld keyword not wired up yet.")

  if opt.async
    @async smtasks, hist = taskSolveTree!(dfgl, tree, timeout;
                                          solveKey=solveKey, algorithm=algorithm,
                                          multithread=multithread, smtasks=smtasks, oldtree=oldtree,
                                          verbose=verbose, verbosefid=verbosefid,
                                          drawtree=opt.drawtree, recordcliqs=recordcliqs,
                                          limititers=opt.limititers, downsolve=opt.downsolve, incremental=opt.incremental,
                                          skipcliqids=skipcliqids, delaycliqs=delaycliqs, limititercliqs=limititercliqs)
  else
    smtasks, hist =        taskSolveTree!(dfgl, tree, timeout;
                                          solveKey=solveKey, algorithm=algorithm,
                                          multithread=multithread, smtasks=smtasks, oldtree=oldtree,
                                          verbose=verbose, verbosefid=verbosefid,
                                          drawtree=opt.drawtree, recordcliqs=recordcliqs,
                                          limititers=opt.limititers, downsolve=opt.downsolve, incremental=opt.incremental,
                                          skipcliqids=skipcliqids, delaycliqs=delaycliqs, limititercliqs=limititercliqs)
    @info "Finished tree based init-ference"
  end

  # NOTE copy of data from new tree in to replace outisde oldtree
  oldtree.bt = tree.bt
  oldtree.btid = tree.btid
  oldtree.cliques = tree.cliques
  oldtree.frontals = tree.frontals
  oldtree.eliminationOrder = tree.eliminationOrder
  oldtree.buildTime = tree.buildTime

  if opt.drawtree && opt.async
    @warn "due to async=true, only keeping task pointer, not stopping the drawtreerate task!  Consider not using .async together with .drawtreerate != 0"
    push!(smtasks, treetask)
  else
    dotreedraw[1] = 0
  end

  # if debugging and not async then also print the CSMHistory
  if opt.dbg && !opt.async
    hists = !opt.async ? fetchCliqHistoryAll!(smtasks) : hist
    printCSMHistorySequential(hists, joinLogPath(dfgl,"HistoryCSMAll.txt") )
  end

  return oldtree
end

"""
    solveGrapn!

Just an alias, see documentation for `solveTree!`.
"""
const solveGraph! = solveTree!

"""
    $SIGNATURES
Internal function used for solveCliqUp! to build the incoming upward message (Rx)
"""
function _buildMessagesUp(fg::AbstractDFG,
                          tree::AbstractBayesTree,
                          cliqid,
                          solveKey::Symbol;
                          status=UPSOLVED )
  #
  cliq = getClique(tree, cliqid)
  beliefMessages = Dict{Int, LikelihoodMessage}()
  for child in getChildren(tree, cliq)
    msg = prepCliqueMsgUp(fg, child, solveKey, status)
    push!(beliefMessages, child.id[]=>msg)
  end
  return beliefMessages
end

"""
    $SIGNATURES

Perform inference in the upward direction over one clique in the Bayes tree according to `opt::SolverParams`.

Example
```julia
tree = buildTreeReset!(fg)
hist, upMessageOut = solveCliqUp!(fg, tree, 2)
```

Notes
- Modifies fg with new values
- Calculates up messages from fg if not provided

DevNotes
- Test isfixedlag
- Test recordcliq

Related
[`solveTree!`](@ref), [`buildTreeReset!`](@ref), [`printCliqHistorySummary`](@ref), [`repeatCSMStep!`](@ref), `sandboxStateMachineStep`
"""
function solveCliqUp!(fg::AbstractDFG,
                      tree::AbstractBayesTree,
                      cliqid::Union{CliqueId, Int, Symbol},
                      solveKey::Symbol=:default,
                      beliefMessages::Dict{Int, LikelihoodMessage} = _buildMessagesUp(fg, tree, cliqid, solveKey); # create belief message from fg if needed
                      verbose::Bool=false,
                      recordcliq::Bool=false )
                      # cliqHistories = Dict{Int,Vector{CSMHistoryTuple}}(),
  #

  # hist = Vector{CSMHistoryTuple}()
  opt = DFG.getSolverParams(fg)

  olddown = opt.downsolve
  opt.downsolve = false
  #TODO test
  if opt.isfixedlag
      @info "Quasi fixed-lag is enabled (a feature currently in testing)!"
      fifoFreeze!(fg)
  end

  cliq = getClique(tree, cliqid)

  # TODO improve, perhaps add to constructor, sommer add all channels here regardless.
  initTreeMessageChannels!(tree)

  @debug "putting messages on up channels from $(keys(beliefMessages))"
  # put the up messages (beliefMessages) that will be used to solve this clique on the channel, the input
  for (id, msg) in pairs(beliefMessages)
    child = getClique(tree, id)
    for e in getEdgesParent(tree, child)
      @async putBeliefMessageUp!(tree, e, msg)
    end
  end

  #
  @debug "taking belief message that will be sent up"
  # take! the message that is sent up by this clique, the output
  takeUpTask = @async takeBeliefMessageUp!(tree, getEdgesParent(tree, cliq)[1])

  recordcliqs = recordcliq ?  [getFrontals(cliq)[1]] : Symbol[]

  hist = tryCliqStateMachineSolve!( fg, tree, cliq.id;
                                    solveKey=solveKey,
                                    verbose=verbose, drawtree=opt.drawtree,
                                    limititers=opt.limititers, downsolve=false,
                                    recordcliqs=recordcliqs, incremental=opt.incremental)
  #

  # post-hoc store possible state machine history in clique (without recursively saving earlier history inside state history)
  # assignTreeHistory!(tree, cliqHistories)
  beliefMessageOut = fetch(takeUpTask)
  #restore downsolve
  opt.downsolve = olddown

  return hist, beliefMessageOut
end

"""
    $SIGNATURES
Internal function used for solveCliqDown! to build the incoming downward message (Rx)
"""
function _buildMessageDown( fg::AbstractDFG,
                            tree::AbstractBayesTree,
                            cliqid,
                            solveKey::Symbol;
                            status::CliqStatus=DOWNSOLVED)
  #
  cliq = getClique(tree, cliqid)
  parent = getParent(tree, cliq)[1]
  return getCliqDownMsgsAfterDownSolve(fg, parent, solveKey; status=status)
end

function solveCliqDown!(fg::AbstractDFG,
                        tree::AbstractBayesTree,
                        cliqid::Union{CliqueId, Int, Symbol},
                        solveKey::Symbol=:default,
                        beliefMessage::LikelihoodMessage = _buildMessageDown(fg, tree, cliqid, solveKey); # create belief message from fg if needed
                        verbose::Bool=false,
                        recordcliq::Bool=false )
  #

  # hist = Vector{CSMHistoryTuple}()
  opt = DFG.getSolverParams(fg)

  upsolve = opt.upsolve

  opt.upsolve = false

  cliq = getClique(tree, cliqid)

  # TODO improve, perhaps add to constructor, sommer add all channels here regardless.
  initTreeMessageChannels!(tree)

  # Build the cliq up message to populate message factors that is needed for down
  @debug "Putting message on up channel from children"
  for (id, msg) in _buildMessagesUp(fg, tree, cliqid, solveKey)
    child = getClique(tree, id)
    for e in getEdgesParent(tree, child)
      @async putBeliefMessageUp!(tree, e, msg)
    end
  end

  # put the down message (beliefMessage) that will be used to solve this clique on the channel, the input
  @debug "putting message on down channel from parent, used by this clique"
  for e in getEdgesParent(tree, cliq)
      @async putBeliefMessageDown!(tree, e, beliefMessage)
  end

  #take! and discart the up message sent in the skip up part of the solve
  @debug "taking belief message that will be sent up"
  @async takeBeliefMessageUp!(tree, getEdgesParent(tree, cliq)[1])

  #
  @debug "taking belief message that will be sent down"
  # take! the message that is sent down by this clique, the output
  takeDownTask = @async begin
    messages = Dict{Int,LikelihoodMessage}()
    for e in getEdgesChildren(tree, cliq)
      messages[e.dst] = takeBeliefMessageDown!(tree, e)
    end
    messages
  end

  recordcliqs = recordcliq ?  [getFrontals(cliq)[1]] : Symbol[]

  hist = tryCliqStateMachineSolve!( fg, tree, cliq.id;
                                    solveKey=solveKey,
                                    verbose=verbose, drawtree=opt.drawtree, limititers=opt.limititers,
                                    recordcliqs=recordcliqs, incremental=opt.incremental)

  # fetch on down
  beliefMessageOut = fetch(takeDownTask)

  #restore
  opt.upsolve = upsolve

  return hist, beliefMessageOut
end


text:
function migrate!(agent::Ind, model::ABM)
  if model.migration_traits[agent.species] == 0
    return
  elseif model.migration_thresholds[agent.species] > get_migration_trait(agent, model)  # if migration value is lower than the threshold
    return
  end

  sites = collect(nearby_positions(agent, model, model.vision_radius[agent.species]))
  if model.check_fraction[agent.species] == 0
    destination = rand(sites)
  else
    nsites = length(sites)
    nsites_selected = ceil(Int, model.check_fraction[agent.species] * nsites)
    sites_checked = EvoDynamics.sample(sites, nsites_selected, replace=false)
    # check sites and move to the best one
    distance, site_index = pick_site(agent, sites_checked, nsites_selected, model)
    current_place = check_site(agent, agent.pos, model)
    if distance < current_place
      return
    end
    destination = sites_checked[site_index]
  end
  # NB: add migration cost? If it survives migration, then do the migration.
  move_agent!(agent, destination, model)
  # agent's base fitness is its abiotic fitness
  agent.W = abiotic_fitness(agent, model)
  adjust_fitness!(agent, model)
  return
end

"Agent evaluates the site and gives it a score"
function check_site(agent, site, model)
  phenotype = agent.abiotic_phenotype
  optimal = return_opt_phenotype(agent.species, model.step[1], site, model)
  abiotic_distance(phenotype, optimal, variance)
end

function pick_site(agent, sites, nsites, model)
  scores = Array{Float64, 1}(undef, nsites)
  for ii in 1:nsites
    scores[ii] = check_site(agent, sites[ii], model)
  end
  return findmin(scores)
end

function get_migration_trait(agent, model)
  mtrait = agent.pleiotropyMat[model.migration_traits[agent.species], :]
  sum(agent.epistasisMat[mtrait, :] * agent.q)
end


text:
using CellListMap
using StaticArrays
using Test

@testset "CellListMap.jl" begin

  if Threads.nthreads() == 1
    println("""

       WARNING: Ideally, run a multi-threaded test to check the parallel versions.

    """)
  end

  # Number of particles, sides and cutoff
  sides = @SVector [250,250,250]
  cutoff = 10.
  box = Box(sides,cutoff)

  # Particle positions
  N = 2000
  x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  # Add some pathological coordinates
  x[1] = -sides/2
  x[10] = sides
  x[100] = @SVector [sides[1]/2, -sides[2]/2, 2*sides[3]]
  y = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]

  # Initialize auxiliary linked lists
  cl = CellList(x,box)

  # Function to be evalulated for each pair: sum of displacements on x
  f(x,y,avg_dx) = avg_dx + abs(x[1] - y[1])

  naive = CellListMap.map_naive!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,x,box)
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=false) ≈ naive

  # Check if changing lcell breaks something
  box = Box(sides,cutoff,lcell=1); cl = CellList(x,box)
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
  box = Box(sides,cutoff,lcell=2); cl = CellList(x,box)
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
  box = Box(sides,cutoff,lcell=3); cl = CellList(x,box)
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive
  box = Box(sides,cutoff,lcell=5); cl = CellList(x,box)
  @test map_pairwise!((x,y,i,j,d2,avg_dx) -> f(x,y,avg_dx),0.,box,cl,parallel=true) ≈ naive

  # Function to be evalulated for each pair: build distance histogram
  function build_histogram!(d2,hist)
    d = sqrt(d2)
    ibin = floor(Int,d) + 1
    hist[ibin] += 1
    return hist
  end

  naive = CellListMap.map_naive!(
    (x,y,i,j,d2,hist) -> build_histogram!(d2,hist),
    zeros(Int,10),x,box
  )
  @test map_pairwise!(
    (x,y,i,j,d2,hist) -> build_histogram!(d2,hist),
    zeros(Int,10),box,cl,parallel=true
  ) ≈ naive
  @test map_pairwise!(
    (x,y,i,j,d2,hist) -> build_histogram!(d2,hist),
    zeros(Int,10),box,cl,parallel=false
  ) ≈ naive

  # Function to be evalulated for each pair: build distance histogram
  function potential(i,j,d2,u,mass)
    d = sqrt(d2)
    u = u - 9.8*mass[i]*mass[j]/d
    return u
  end

  # Run pairwise computation
  mass = rand(N)
  naive = CellListMap.map_naive!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,x,box)
  @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,box,cl,parallel=true) ≈ naive
  @test map_pairwise!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,box,cl,parallel=false) ≈ naive

  # Function to be evalulated for each pair: build distance histogram
  function calc_forces!(x,y,i,j,d2,mass,forces)
    G = 9.8*mass[i]*mass[j]/d2
    d = sqrt(d2)
    df = (G/d)*(x - y)
    forces[i] = forces[i] - df
    forces[j] = forces[j] + df
    return forces
  end

  # forces
  forces = [ zeros(SVector{3,Float64}) for i in 1:N ]

  # Run pairwise computation
  naive = CellListMap.map_naive!(
    (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces),
    copy(forces),x,box
  )
  @test map_pairwise!(
    (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces),
    copy(forces),box,cl,parallel=true
  ) ≈ naive
  @test map_pairwise!(
    (x,y,i,j,d2,forces) -> calc_forces!(x,y,i,j,d2,mass,forces),
    copy(forces),box,cl,parallel=false
  ) ≈ naive

  # Compute some properties of disjoint sets
  box = Box(sides,cutoff,lcell=1)
  cl = CellList(x,y,box)

  naive = CellListMap.map_naive!((x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),0.0,x,y,box)
  @test map_pairwise!(
    (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),
    0.0,box,cl,parallel=false
  ) ≈ naive
  @test map_pairwise!(
    (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),
    0.0,box,cl,parallel=true
  ) ≈ naive

  # Check different lcell
  box = Box(sides,cutoff,lcell=3)
  cl = CellList(x,y,box)
  @test map_pairwise!(
    (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),
    0.0,box,cl,parallel=false
  ) ≈ naive
  @test map_pairwise!(
    (x,y,i,j,d2,u) -> potential(i,j,d2,u,mass),
    0.0,box,cl,parallel=true
  ) ≈ naive

  # Test the examples, to check further if the parallelization didn't break something
  N = 100_000
  x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  y = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
  @test CellListMap.test1(parallel=true,x=x)[2] ≈ CellListMap.test1(parallel=false,x=x)[2]
  @test CellListMap.test2(parallel=true,x=x) ≈ CellListMap.test2(parallel=false,x=x)
  @test CellListMap.test3(parallel=true,x=x) ≈ CellListMap.test3(parallel=false,x=x)
  @test CellListMap.test4(parallel=true,x=x) ≈ CellListMap.test4(parallel=false,x=x)
  @test count(CellListMap.test5(parallel=true,x=x,y=y) .≈ CellListMap.test5(parallel=false,x=x,y=y)) == 3

  pairs1 = sort!(CellListMap.test7(parallel=true,x=x),by=x->x[3])
  pairs2 = sort!(CellListMap.test7(parallel=false,x=x),by=x->x[3])
  @test count([ count(pairs1[i] .≈ pairs2[i]) == 3 for i in 1:length(pairs1) ]) == length(pairs1)

  x = [ sides .* rand(SVector{3,Float64}) for i in 1:1_500 ]
  y = [ sides .* rand(SVector{3,Float64}) for i in 1:1_500_000 ]
  @test count(CellListMap.test6(parallel=true,x=x,y=y) .≈ CellListMap.test6(parallel=false,x=x,y=y)) == 3

  # invert x and y to test swap
  ixy = CellListMap.test6(parallel=false,x=x,y=y)
  iyx = CellListMap.test6(parallel=false,x=y,y=x)
  @test ( ixy[1] == iyx[2] && ixy[2] == iyx[1] && ixy[3] ≈ iyx[3] )
  ixy = CellListMap.test6(parallel=true,x=x,y=y)
  iyx = CellListMap.test6(parallel=true,x=y,y=x)
  @test ( ixy[1] == iyx[2] && ixy[2] == iyx[1] && ixy[3] ≈ iyx[3] )

#  # Test resizing of the cell lists
#  x = [ rand(SVector{3,Float64}) for i in 1:1000 ]
#  box = Box([0.83,0.41,0.97],0.1)
#  cl = CellList(x,box)
#  @test length(cl.cwp) == 924
#
#  box = Box([0.33,0.41,0.97],0.1)
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 924
#
#  box = Box([0.83,0.81,0.97],0.1)
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 1598
#
#  x .= 0.9*x
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 1598
#
#  x .= 1.2*x
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 1598
#
#  box = Box([0.83,0.81,0.97],0.2)
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 1598
#
#  box = Box([0.83,0.81,0.97],0.05)
#  cl = UpdateCellList!(x,box,cl)
#  @test length(cl.cwp) == 8737

end


text:
input = open("../inputs/day03.txt") do f
    read(f, String)
end

function makeMap(inp)
    lines = split(inp, "\n")
    height = length(lines)
    width = length(lines[1])
    m = fill(' ', height, width)
    for y in 1:height
        for (x, c) in enumerate(lines[y])
            m[y,x] = c
        end
    end
    m
end

function runSlope(m, dy, dx)::Int64
    width = size(m, 2)
    height = size(m, 1)
    x = 1 + dx
    y = 1 + dy
    count = 0
    while y <= height
        if m[y,x] == '#'
            count += 1
        end
        x = mod(x + dx, 1:width)
        y += dy
    end
    count
end

m = makeMap(input)
println("Part 1: ", runSlope(m, 1, 3))

slopes = [(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)]

println("Part 2: ", prod([runSlope(m, dy, dx) for (dx, dy) in slopes]))


text:
using BoardGames

struct RandomStrategy <: Strategy end

function BoardGames.getmove(board, s::RandomStrategy)
    rand(getmoves(board))
end


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 5ecc9b01-bbcf-41d3-a2e5-51c2d242ba13
begin
	using PlutoUI
	using Test
end

# ╔═╡ 420fe8b0-a0e1-11eb-3d87-8fbf9a8445a6
md"""
## Julia exploration: Channels, iterators

"""

# ╔═╡ 4cb29ef3-be1d-45e2-90c0-5cd6d33f6308
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ b5d4212a-8971-40b9-a540-0b0593f3014a
md"""
### Using channels
"""

# ╔═╡ b84302b6-6582-4d6b-9350-90a39c1c87da
## Producer
##
function fibo(n::Integer)
	Channel{Int64}(20) do ch
		a, b = 0, 1
		for _ ∈ 1:n
			a, b = b, a + b
			put!(ch, a)
		end
	end
end

# ╔═╡ 7e7b494c-0052-4192-b810-144a3c3c7e72
## Consummer
##
with_terminal() do
	for val ∈ fibo(54)
		val % 5 == 0 && println()
    	print(val, " ")
	end
end

# ╔═╡ 32674ee8-2e34-4e81-b8b8-e8fdcfcf4132
md"""
Another example.
"""

# ╔═╡ c8430355-2128-4206-ae88-319d26557f60
function generator(;n::Int=1, buffer_sz=10)::Channel
	@assert n > 0 && buffer_sz > 0

	Channel{Char}(buffer_sz) do ch
		for ix ∈ 0:n-1
			put!(ch, Char('A' + (ix % 26)))
		end
	end
end

# ╔═╡ b5155ffb-b371-45c9-8778-7f6005164022
g = generator()

# ╔═╡ 645a0289-4c3e-49c8-af2c-64555d283e25
collect(g)

# ╔═╡ 39665955-0835-459d-9476-797d50ea10a9
g₂ = generator(n=10, buffer_sz=2)

# ╔═╡ 16ab9011-cbef-4939-b3ef-69029d2cb7e7
collect(g₂)

# ╔═╡ 3dc47adb-154c-4109-bb6b-2223e0d26d9f
# we can use the splat operator
@test [g...] == collect(g)

# ╔═╡ ab40e357-7fac-4220-addb-b3add86ad982
md"""
### Details about Tasks and Channels
"""

# ╔═╡ 1aac9209-a955-430e-a898-1a25f056826c
md"""
Let us rewrite our generator function without the Channel block to explicit what is going on...
"""

# ╔═╡ 0cced5d4-77ba-46c9-b1c4-0119cf4cff46
function generator_alt(;n::Int=1, buffer_sz=10)
	@assert n > 0 && buffer_sz > 0
	ch = Channel(buffer_sz)

	# closure
	function gen_letters()
		for ix ∈ 0:n-1
			put!(ch, Char('A' + (ix % 26)))
		end
	end

	t = Task(gen_letters)    ## create a task for this coroutine
	schedule(t)              ## schedule it to run
	ch
end

# ╔═╡ f37022c3-b4bd-4904-a544-38d83ef67cbf
md"""
Channels are used to send data between coroutines.

Note we cannot have multiple coroutines running on different CPU cores. Coroutines allow for cooperative multi-tasking. The scheduler switches between coroutine, because the current one has explicitly give up control.

We specified a buffer of size 2 and when this buffer gets full, the coroutine will block and the scheduler will pick another coroutine (task) waiting for this particular channel.

Every coroutine needs a task, while a task keeps track of the execution state.
"""

# ╔═╡ 23930be8-ab35-4cae-9563-99fcd571a401
begin
	## consumer
	ch_gen = generator_alt(n=5)
	take!(ch_gen)   ## consume 1 char
end

# ╔═╡ 45046186-8f4c-44e1-9097-5b39ac678cd4
take!(ch_gen)

# ╔═╡ 4e259e35-b318-49bb-97c3-839354169b20
with_terminal() do
	for _ ∈ 1:3
		println(take!(ch_gen))
	end
	# take!(ch_gen)  ## would block forever
end

# ╔═╡ 500c5603-e0cd-40b2-80c4-6cbd6bbef7cc
md"""
Finally an alternative with the `@async` macro.
"""

# ╔═╡ d34b4627-4d7e-41be-b32d-7455a1a33105
function generator_alt₂(;n::Int=1, buffer_sz=10)
	@assert n > 0 && buffer_sz > 0
	ch = Channel(buffer_sz)

	@async begin
		for ix ∈ 0:n-1
			put!(ch, Char('A' + (ix % 26)))
		end
	end
	ch
end

# ╔═╡ 5820cab3-84e0-4c3d-a910-43a59b5ed24e
with_terminal() do
	## consumer
	n = 3
	ch_gen₂ = generator_alt₂(;n)
	for _ ∈ 1:n
		println(take!(ch_gen₂))
	end
end

# ╔═╡ 8a65122f-d236-42cc-816c-5686c8202b34
html"""
<style>
  main {
	max-width: calc(800px + 25px + 6px);
  }
  .plutoui-toc.aside {
	background: linen;
  }
  h3, h4 {
	background: wheat;
  }
</style>
"""

# ╔═╡ Cell order:
# ╟─420fe8b0-a0e1-11eb-3d87-8fbf9a8445a6
# ╠═5ecc9b01-bbcf-41d3-a2e5-51c2d242ba13
# ╟─4cb29ef3-be1d-45e2-90c0-5cd6d33f6308
# ╟─b5d4212a-8971-40b9-a540-0b0593f3014a
# ╠═b84302b6-6582-4d6b-9350-90a39c1c87da
# ╠═7e7b494c-0052-4192-b810-144a3c3c7e72
# ╟─32674ee8-2e34-4e81-b8b8-e8fdcfcf4132
# ╠═c8430355-2128-4206-ae88-319d26557f60
# ╠═b5155ffb-b371-45c9-8778-7f6005164022
# ╠═645a0289-4c3e-49c8-af2c-64555d283e25
# ╠═39665955-0835-459d-9476-797d50ea10a9
# ╠═16ab9011-cbef-4939-b3ef-69029d2cb7e7
# ╠═3dc47adb-154c-4109-bb6b-2223e0d26d9f
# ╟─ab40e357-7fac-4220-addb-b3add86ad982
# ╟─1aac9209-a955-430e-a898-1a25f056826c
# ╠═0cced5d4-77ba-46c9-b1c4-0119cf4cff46
# ╟─f37022c3-b4bd-4904-a544-38d83ef67cbf
# ╠═23930be8-ab35-4cae-9563-99fcd571a401
# ╠═45046186-8f4c-44e1-9097-5b39ac678cd4
# ╠═4e259e35-b318-49bb-97c3-839354169b20
# ╟─500c5603-e0cd-40b2-80c4-6cbd6bbef7cc
# ╠═d34b4627-4d7e-41be-b32d-7455a1a33105
# ╠═5820cab3-84e0-4c3d-a910-43a59b5ed24e
# ╟─8a65122f-d236-42cc-816c-5686c8202b34

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
module AlgebraTests

using Test

include("NonLinearOperatorMocks.jl")

@testset "LinearSolvers" begin include("LinearSolversTests.jl") end

@testset "NonLinearSolvers" begin include("NonLinearSolversTests.jl") end

@testset "NLSolvers" begin include("NLSolversTests.jl") end

end # module AlgebraTests


text:
function _is_within_range(query_start, query_length, target_start, target_length)
    target_start <= query_start <= query_start - 1 + query_length <= target_start - 1 + target_length
end

abstract type Overlap end
struct ExactMatch <: Overlap end
struct LeftNeighbor <: Overlap end
struct RightNeighbor <: Overlap end
struct Disjunct <: Overlap end
struct Within <: Overlap end
struct WithinFromStart <: Overlap end
struct WithinToEnd <: Overlap end
struct CompleteOverlap <: Overlap end
struct LeftOverlap <: Overlap end
struct RightOverlap <: Overlap end

function _determine_overlap(query_start, query_length, target_start, target_length)
    query_end = query_start - 1 + query_length
    target_end = target_start - 1 + target_length

    if query_start == target_start && query_length == target_length
        ExactMatch()
    elseif query_end == target_start - 1
        LeftNeighbor()
    elseif query_start == target_end + 1
        RightNeighbor()
    elseif query_end < target_start || target_end < query_start
        Disjunct()
    elseif target_start == query_start <= query_end < target_end
        WithinFromStart()
    elseif target_start < query_start <= query_end == target_end
        WithinToEnd()
    elseif target_start < query_start <= query_end  < target_end
        Within()
    elseif query_start < target_start <= query_end <= target_end
        LeftOverlap()
    elseif target_start <= query_start <= target_end < query_end
        RightOverlap()
    elseif query_start < target_start <= target_end  < query_end
        CompleteOverlap()
    else
        println(query_start, " ", query_length, " ", target_start, " ", target_length)
        @assert false
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
# # [Interpreting h′ and β](@id interpreting-hpbeta)
#
# This example looks at the influence of h′, β, electron collision frequency,
# and transmitter frequency on amplitude curves.
#
# ## Background
#
# [Wait and Spies (1964)](https://archive.org/details/characteristicso300wait) is a
# technical report which characterizes the influence of Earth-ionosphere waveguide
# parameters on the propagation of VLF radio waves (the attenuation rates, phase velocities,
# etc.).
# The report is probably most cited[^1] as the source of the exponential electron
# density profile parameterized by a height h′ (km) and slope β (km⁻¹).
# The profile is usually written as:
# ```math
# N(z) = 1.43\times 10^{13} \exp(-0.15 h') \exp((\beta - 0.15)(z - h'))
# ```
# for density ``N`` (m⁻³) and altitude ``z`` (km), but as described in
# [Morfitt (1977)](https://apps.dtic.mil/docs/citations/ADA047508) or
# [Thomson (1993)](https://doi.org/10.1016/0021-9169(93)90122-F), Wait and Spies (1964)
# does not explicitly introduce this form.
# Instead, Wait introduces a "conductivity parameter" ``\omega_r`` defined by
# ```math
# \omega_r = \omega_p^2/\nu
# ```
# where ``\omega_p`` is the angular plasma frequency of the electrons and ``\nu`` is the
# effective collision frequency.
# The height at which ``\omega_r = \omega``, the angular radio wave frequency is
# approximately the wave reflection height. This statement is equivalent to ``X = Z`` from
# magnetoionic theory ([Ratcliffe (1959)](https://books.google.com/books/about/The_magneto_ionic_theory_and_its_applica.html?id=lRQvAAAAIAAJ)).
#
# Citing laboratory and some observational experiments, Wait assumes an analytical
# collision frequency profile
# ```math
# \nu = 5 \times 10^6 \exp(-0.15(z - 70)).
# ```
# Using daytime measurements made by the method of partial reflection, he calculates
# the approximation
# ```math
# \omega_r = (2.5\times 10^5)\exp(0.3(z - 70)).
# ```
# To fit measurements made at other times and locations, he parameterizes the function
# as
# ```math
# \omega_r = (2.5 \times 10^5)\exp(\beta(z - h'))
# ```
# where "``\beta`` is a constant and ``h'`` is a reference height" (the height at which
# ``\omega_r = 2.5\times 10^5``).
# He immediately suggests ``\beta = 0.3`` at day and ``\beta = 0.5`` at night, which are
# largely considered typical values today.
# Wait also parameterizes the collision frequency profile as
# ```math
# \nu(z) = \nu_0\exp(-az).
# ```
# Wait fixed ``a=0.15`` for his studies in the report.
#
# [Morfitt and Shellman (1976)](http://www.dtic.mil/docs/citations/ADA032573)
# cites Wait and Spies and uses
# ```math
# \nu(z) = 1.82\times 10^{11} \exp(-0.15 z)
# ```
# from this, it is simple to derive the equation attributed to
# Wait and Spies (1964). The ``(\beta - 0.15)`` term comes from the assumption
# ``N(z) = N_0\exp(bz)`` and ``\nu(z) = \nu_0\exp(-az)`` where ``\beta = a + b`` and
# ``a = 0.15``.

using Printf, Plots, Interpolations
using LongwaveModePropagator
const LMP = LongwaveModePropagator

## Constant values
const BFIELD = BField(50e-6, π/2, 0)
const GND = GROUND[5]

const TX = Transmitter(20e3)
const RX = GroundSampler(0:5e3:3000e3, Fields.Ez)
nothing  #hide

# ## Varying h′
#
# In this section we set the frequency to 20 kHz, β is set to 0.4 km⁻¹ and we vary h′.

function varyhp(hprimes)
    amps = Vector{Vector{Float64}}(undef, length(hprimes))
    for i in eachindex(hprimes)
        electrons = Species(LMP.QE, LMP.ME, z->waitprofile(z, hprimes[i], 0.4), electroncollisionfrequency)
        waveguide = HomogeneousWaveguide(BFIELD, electrons, GND)
        E, a, p = propagate(waveguide, TX, RX)
        amps[i] = a
    end
    return amps
end

hprimes = 72:1:78
amps = varyhp(hprimes)

p = plot();
function buildplots!(p, amps)
    cmap = palette(:amp, length(hprimes)+1) # +1 allows us to use a darker lightest color

    for i in eachindex(hprimes)
        plot!(p, RX.distance/1000, amps[i];
              label=hprimes[i], color=cmap[i+1]);
    end
end
buildplots!(p, amps);
plot!(p; size=(600,400), ylims=(22, 95),
      xlabel="Range (km)", ylabel="Amplitude (dB)", legendtitle="h′")
#md savefig(p, "interpreting_hprimes.png"); nothing # hide
#md # ![](interpreting_hprimes.png)

# As h′ increases, the field strength modal interference pattern is displaced outwardly
# away from the transmitter.

# ## Varying β
#
# Now we'll fix h′ at 78 km and vary β.

function varybeta(betas)
    amps = Vector{Vector{Float64}}(undef, length(betas))
    for i in eachindex(betas)
        electrons = Species(LMP.QE, LMP.ME, z->waitprofile(z, 78, betas[i]), electroncollisionfrequency)
        waveguide = HomogeneousWaveguide(BFIELD, electrons, GND)
        E, a, p = propagate(waveguide, TX, RX)
        amps[i] = a
    end
    return amps
end

betas = [0.3, 0.4, 0.5, 0.7, 0.9, 2.0]
amps = varybeta(betas)

p = plot();
function buildplots!(p, amps)
    cmap = palette(:amp, length(betas)+1)

    for i in eachindex(betas)
        plot!(p, RX.distance/1000, amps[i];
              label=betas[i], color=cmap[i+1]);
    end
end
buildplots!(p, amps);
plot!(p; size=(600,400), ylims=(22, 95),
      xlabel="Range (km)", ylabel="Amplitude (dB)", legendtitle="β")
#md savefig(p, "interpreting_betas.png"); nothing # hide
#md # ![](interpreting_betas.png)

# The signal roughly increases with increasing β.
# Higher β profiles more closely represent a sharp reflecting boundary.
#
# It may be more accurate to say that higher β ionospheres lead to more extreme amplitudes -
# higher highs and lower lows.
# Morfitt (1977) describes this as an increase in β leads to the signal levels varying
# over a wider range in the regions of strong modal interference.

# ## Varying frequency
#
# In this section we'll look at an ionosphere described with typical daytime
# Wait profile parameters, but we'll vary the radio frequency from 5 kHz to 50 kHz.

function varyfreq(freqs)
    electrons = Species(LMP.QE, LMP.ME, z->waitprofile(z, 75, 0.35), electroncollisionfrequency)
    waveguide = HomogeneousWaveguide(BFIELD, electrons, GND)

    amps = Vector{Vector{Float64}}(undef, length(freqs))
    for i in eachindex(freqs)
        tx = Transmitter(freqs[i])
        E, a, p = propagate(waveguide, tx, RX)
        amps[i] = a
    end
    return amps
end

freqs = 5e3:5e3:50e3
amps = varyfreq(freqs)

p = plot();
function buildplots!(p, amps)
    pal = palette(:rainbow, 7)
    cmap = [pal[1]; range(pal[2], pal[3]; length=2); range(pal[4], pal[5]; length=4);
            range(pal[6], pal[7]; length=3)]

    for i in eachindex(freqs)
        fkHz = trunc(Int, freqs[i]/1000)
        λkm = trunc(LMP.C0/freqs[i]/1000; digits=1)
        plot!(p, RX.distance/1000, amps[i] .+ (10*i);
              label=string(fkHz, ",  ", λkm), color=cmap[i]);
    end
end
buildplots!(p, amps);
plot!(p; size=(600,400),
      xlabel="Range (km)", ylabel="Amplitude (dB)", legendtitle="f kHz, λ km")
#md savefig(p, "interpreting_frequencies.png"); nothing # hide
#md # ![](interpreting_frequencies.png)

# The effect of increasing frequency is similar to increasing h′.
# This makes sense because in both cases the waveguide "height" in wavelengths is
# increassing.

# ## Collision frequency
#
# So far all of the profiles generated above have used the electron collision frequency
# from Morfitt and Shellman (1976).
# Let's look at the effect of different collision frequency profiles.
#
# First we'll plot the collision frequency profiles themselves.
# On the plots we'll also mark the approximate reflection height for `TX.frequency` where
# ``\omega_r = \omega`` for a typical daytime electron density profile.

collisionfrequency(z, ν₀, a) = ν₀*exp(-a*z/1000)  # 1/1000 converts `z` to km

function reflectionheight(params, hp, β)
    species = Species(LMP.QE, LMP.ME,
        z->waitprofile(z, hp, β), z->collisionfrequency(z, params...))
    ωr = LMP.waitsparameter.(alt, (TX.frequency,), (BFIELD,), (species,))
    itp = LinearInterpolation(ωr, alt)
    eqz = itp(TX.frequency.ω)

    return eqz
end

alt = 40e3:500:110e3
params = [(1.816e11, 0.18),
          (1.816e11/5, 0.15), (1.816e11, 0.15), (1.816e11*5, 0.15),
          (1.816e11, 0.12)]

p = plot();
function buildplots!(p, params, hp, β)
    cmap = palette(:rainbow, length(params))

    for i in eachindex(params)
        eqz = reflectionheight(params[i], hp, β)

        v₀label = params[i][1]/1.816
        plot!(p, collisionfrequency.(alt, params[i]...), alt/1000;
              label=@sprintf("%.0e, %.2f", v₀label, params[i][2]), color=cmap[i]);
        hline!(p, [eqz/1000]; linestyle=:dash, color=cmap[i], linewidth=0.6, label=nothing)
    end
end
buildplots!(p, params, 75, 0.35);
annotate!(p, [(1e3, 64, text("ωᵣ = ω", 10))])
plot!(p; size=(600,400), xscale=:log10,
      xlabel="ν (s⁻¹)", ylabel="Altitude (km)", legendtitle="1.816 ν₀, a")
#md savefig(p, "interpreting_collisionprofile.png"); nothing # hide
#md # ![](interpreting_collisionprofile.png)

# Typical daytime electron density profile.

function varycollisionfreq(params, hp, β)
    amps = Vector{Vector{Float64}}(undef, length(params))
    for i in eachindex(params)
        electrons = Species(LMP.QE, LMP.ME,
            z->waitprofile(z, hp, β), z->collisionfrequency(z, params[i]...))
        waveguide = HomogeneousWaveguide(BFIELD, electrons, GND)
        E, a, p = propagate(waveguide, TX, RX)
        amps[i] = a
    end
    return amps
end
amps = varycollisionfreq(params, 75, 0.35)

p = plot();
function buildplots!(p, amps)
    cmap = palette(:rainbow, length(params))

    for i in eachindex(params)
        v₀label = params[i][1]/1.816
        plot!(p, RX.distance/1000, amps[i];
              label=@sprintf("%.0e, %.2f", v₀label, params[i][2]), color=cmap[i]);
    end
end
buildplots!(p, amps);
plot!(p; size=(600,400), ylims=(22, 95),
      xlabel="Range (km)", ylabel="Amplitude (dB)", legendtitle="1.816 ν₀, a")
#md savefig(p, "interpreting_collisionday.png"); nothing # hide
#md # ![](interpreting_collisionday.png)

# Typical nighttime electron density profile.

amps = varycollisionfreq(params, 82, 0.55)

p = plot();
function buildplots!(p, amps)
    cmap = palette(:rainbow, length(params))

    for i in eachindex(params)
        v₀label = params[i][1]/1.816
        plot!(p, RX.distance/1000, amps[i];
              label=@sprintf("%.0e, %.2f", v₀label, params[i][2]), color=cmap[i]);
    end
end
buildplots!(p, amps);
plot!(p; size=(600,400), ylims=(22, 95),
      xlabel="Range (km)", ylabel="Amplitude (dB)", legendtitle="1.816 ν₀, a")
#md savefig(p, "interpreting_collisionnight.png"); nothing # hide
#md # ![](interpreting_collisionnight.png)

# [^1]: As a government technical report and not a journal paper curated by a publishing company, it is difficult to track how many times it has been cited by others.


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LittleCMS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LittleCMS")
JLLWrappers.@generate_main_file("LittleCMS", UUID("d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"))
end  # module LittleCMS_jll


text:
using Flux
using Images
using FileIO
using Glob
using Mmap
using ImageShow
using Base.Iterators: partition
using Flux: onehotbatch,onecold, crossentropy, throttle
using Statistics:mean

# Read CIFAR10 data.(https://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz)
# Data is given in binary format where the first byte is the label of
# the first image which is a number between 0 and 9. The next 3072 bytes
# are the values of the pixels. The first 1024 are R, the next 1024 are G and
# the last 1024 are B.

function loadBatches(path)
    batch = Int.(open(read,path))
    labels = [batch[1+(n-1)*3073] for n in 1:1000]
    pixels = [reshape(batch[2+(n-1)*3073:3073+(n-1)*3073],(32,32,3)) for n in 1:1000]/255
    return (pixels,labels)
end

# Get file names
path = "./cifar-10-batches-bin/"
trainbatch = readdir(glob"data_batch_*.bin",path)
testbatch = readdir(glob"test_batch.bin",path)

@show("[INFO] Reading Data...")
# read files and prepare train and test datasets
for file in trainbatch
    if  file==trainbatch[1]
        global X_t, Y_train = loadBatches(file)
    else
        data = loadBatches(file)
        append!(X_t,data[1])
        append!(Y_train,data[2])
    end
end

# Reshape Xt
X_train =[]
for i in 1:size(X_t)[1]
    append!(X_train,X_t[i])
end
X_train = Float32.(reshape(X_train,32,32,3,:))

X_tt,Y_test = loadBatches(testbatch[1])

X_test =[]
for i in 1:size(X_tt)[1]
    append!(X_test,X_tt[i])
end
X_test = Float32.(reshape(X_test,32,32,3,:))

function make_partition_index(X,batch_size)
    idx = partition(1:Int32(length(X)/(32*32*3)),batch_size)
    indices = [(minimum(i),maximum(i)) for i in idx]
    return indices
end

function make_minibatch(X,Y,batch_size)
    indices = [i for i in make_partition_index(X,batch_size)]
    minibatch_X = [X[:,:,:,indices[i][1]:indices[i][2]] for i in 1:length(indices)]
    minibatch_Y = [Y[indices[i][1]:indices[i][2]] for i in 1:length(indices)]
    dataset = [(minibatch_X[i],onehotbatch(minibatch_Y[i],0:9)) for i in 1:length(indices)]
    return dataset
end

@show("[INFO] Creating minibatches...")
# Create minibatches
train_set = gpu.(make_minibatch(X_train,Y_train,128))
test_set = gpu.(make_minibatch(X_test,Y_test,1))
# allowscalar(false)

# VGG16
@show("[INFO] Building CNN...")
model() = Chain(
    # Size 32x32
    Conv((3,3), 3=>64,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(64),
    Dropout(0.3),
    # Size 32x32
    Conv((3,3), 64=>64,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(64),
    # Size 32x32
    x -> MaxPool((2,2))(x),
    # Size 16x16
    Conv((3,3), 64=>128,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(128),
    Dropout(0.3),
    # Size 16x16
    Conv((3,3), 128=>128,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(128),
    x -> MaxPool((2,2))(x),
    # Size 8x8
    Conv((3,3), 128=>256,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(256),
    Dropout(0.4),
    # Size 8x8
    Conv((3,3),256=>256, relu, pad=(1,1), stride=(1,1)),
    BatchNorm(256),
    Dropout(0.4),
    # Size 8x8
    Conv((3,3),256=>256, relu, pad=(1,1), stride=(1,1)),
    BatchNorm(256),
    x -> MaxPool((2,2))(x),
    # Size 4x4
    Conv((3,3), 256=>512,relu, pad=(1,1), stride=(1,1)),
    BatchNorm(512),
    Dropout(0.4),
    # Size 4x4
    Conv((3,3),512=>512, relu, pad=(1,1), stride=(1,1)),
    BatchNorm(512),
    Dropout(0.4),
    # Size 2x2
    x -> reshape(x,:,size(x,4)),
    Dense(8192,4096,relu),
    Dropout(0.5),
    Dense(4096,10),
    softmax
)|>gpu

# Make model

m = model()

loss(x,y) = crossentropy(m(x),y)

# Calculate accuracy
function accuracy(x)
    yhat = [onecold(m(x[i][1]),1:10) for i in 1:size(x)[1]]
    yact = [onecold(x[i][2],1:10) for i in 1:size(x)[1]]
    return mean(yhat .== yact)
end

evalcb = throttle(()-> @show(accuracy(test_set)),10)

opt = ADAM()

@show("[INFO] Training ...")
n_epochs = 100
for i in 1:n_epochs
    Flux.train!(loss, params(m),train_set,opt, cb=evalcb)
end


text:
using Pythia, Test

# @time @testset "Basic Methods Tests" begin include("basicMethods_test.jl") end
@time @testset "Exponential Smoothing Tests" begin include("ETS_test.jl") end


text:
convert(::Type{DoubleFloat{T}}, x::T) where {T<:IEEEFloat} = DoubleFloat{T}(x)
convert(::Type{T}, x::DoubleFloat{T}) where {T<:IEEEFloat} = HI(x)

convert(::Type{DoubleFloat{T}}, x::I) where {T<:IEEEFloat, I<:Integer} = DoubleFloat{T}(T(x))
convert(::Type{I}, x::DoubleFloat{T}) where {T<:IEEEFloat, I<:Integer} = I(round(BigInt, BigFloat(HI(x))+ BigFloat(LO(x))))

function convert(::Type{BigFloat}, x::DoubleFloat{T}) where {T<:AbstractFloat}
    res = Base.BigFloat(HI(x)) + Base.BigFloat(LO(x))
    return res
end
function convert(::Type{DoubleFloat{T}}, x::BigFloat) where {T<:AbstractFloat}
    hi = T(x)
    lo = T(x - hi)
    return DoubleFloat(hi, lo)
end

function convert(::Type{BigInt}, x::DoubleFloat{T}) where {T<:AbstractFloat}
    res = round(Base.BigInt, Base.BigFloat(HI(x)) + Base.BigFloat(LO(x)))
    return res
end
function convert(::Type{DoubleFloat{T}}, x::BigInt) where {T<:AbstractFloat}
    return convert(DoubleFloat{T}, BigFloat(x))
end

convert(::Type{DoubleFloat{T}}, x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = HI(x)
convert(::Type{DoubleFloat{DoubleFloat{T}}}, x::DoubleFloat{T}) where {T<:IEEEFloat} = DoubleFloat(x, DoubleFloat{T}(zero(T)))

convert(::Type{T}, x::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = HI(HI(x))
convert(::Type{DoubleFloat{DoubleFloat{T}}}, x::T) where {T<:IEEEFloat} = DoubleFloat(DoubleFloat{T}(x), DoubleFloat{T}(zero(T)))


text:
using MacroTools: @q

macro nograd(ex)
  isexpr(ex, :tuple) || (ex = Expr(:tuple, ex))
  blk = @q begin end
  for f in ex.args
    back = MacroTools.@q _ -> ($__source__; nothing)
    push!(blk.args, :(@inline Zygote._pullback(::Context, ::Core.Typeof($(esc(f))), args...) = $(esc(f))(args...), $back))
  end
  return blk
end

macro which(ex)
  @capture(ex, f_(args__)) || error("Zygote.@which f(args...)")
  :(InteractiveUtils.@which adjoint(Context(), $(esc(f)), $(esc.(args)...)))
end

"""

    checkpointed(f, xs...)

Use gradient checkpointing on the call `f(xs...)`. This means that
`checkpointed(f, xs...) === f(xs...)`, but when computing the derivative
intermediate results from the forward pass of `f` will not be stored. Instead the forward
pass will be repeated, when computing the derivative.
This saves memory at the cost of increasing exectution time.

!!! warning
    If `f` is not a pure function, `checkpointed` will likely give wrong results.
"""
checkpointed(f, xs...) = f(xs...)

function Zygote._pullback(ctx::Zygote.AContext, ::typeof(checkpointed), f, xs...)
    y = f(xs...)
    function pullback_checkpointed(Δy)
        y, pb = Zygote._pullback(ctx, f, xs...)
        return (nothing, pb(Δy)...)
    end
    return y, pullback_checkpointed
end

"""
    hessian(f, x)

Construct the Hessian `∂²f/∂x²`, where `x` is a real number or an array,
and `f(x)` is a real number. When `x` is an array, the result is a matrix
`H[i,j] = ∂²f/∂x[i]∂x[j]`, using linear indexing `x[i]` even if the argument
is higher-dimensional.

This uses forward over reverse, ForwardDiff over Zygote, calling `hessian_dual(f, x)`.
See [`hessian_reverse`](@ref) for an all-Zygote alternative.

# Examples

```jldoctest; setup=:(using Zygote)
julia> hessian(x -> x[1]*x[2], randn(2))
2×2 Array{Float64,2}:
 0.0  1.0
 1.0  0.0

julia> hessian(x -> sum(x.^3), [1 2; 3 4])  # uses linear indexing of x
4×4 Array{$Int,2}:
 6   0   0   0
 0  18   0   0
 0   0  12   0
 0   0   0  24

julia> hessian(sin, pi/2)
-1.0
```
"""
hessian(f, x) = hessian_dual(f, x)

hessian_dual(f, x::AbstractArray) = forward_jacobian(x -> gradient(f, x)[1], x)[2]

hessian_dual(f, x::Number) = ForwardDiff.derivative(x -> gradient(f, x)[1], x)

"""
    hessian_reverse(f, x)

This should be equivalent to [`hessian(f, x)`](@ref hessian),
but implemented using reverse over reverse mode, all Zygote.
(This is usually much slower, and more likely to find errors.)
"""
hessian_reverse(f, x::AbstractArray) = jacobian(x -> gradient(f, x)[1], x)[1]

hessian_reverse(f, x::Number) = gradient(x -> gradient(f, x)[1], x)[1]

"""
    jacobian(f, args...) -> Tuple

For each array `a ∈ args` this returns a matrix with `Ja[k,i] = ∂y[k]/∂a[i]`
where `y = f(args...)` is usually a vector.
Arrays of higher dimension are treated like `vec(a)`, or `vec(y)` for output.

For scalar `x::Number ∈ args`, the result is a vector `Jx[k] = ∂y[k]/∂x`,
while for scalar `y` all results have just one row.

With any other argument type, no result is produced, even if [`gradient`](@ref) would work.

This reverse-mode Jacobian needs to evaluate the pullback once for each element of `y`.
Doing so is usually only efficient when `length(y)` is small compared to `length(a)`,
otherwise forward mode is likely to be better.

See also [`hessian`](@ref), [`hessian_reverse`](@ref).

# Examples

```jldoctest; setup=:(using Zygote)
julia> jacobian(a -> 100*a[1:3].^2, 1:7)[1]  # first index (rows) is output
3×7 Array{$Int,2}:
 200    0    0  0  0  0  0
   0  400    0  0  0  0  0
   0    0  600  0  0  0  0

julia> jacobian((a,x) -> a.^2 .* x, [1,2,3], 1)  # scalar argument has vector jacobian
([2 0 0; 0 4 0; 0 0 6], [1, 4, 9])

julia> jacobian((a,d) -> prod(a, dims=d), [1 2; 3 4; 5 6], 2)
([2 0 … 0 0; 0 4 … 3 0; 0 0 … 0 5], [0, 0, 0])
```

!!! warning
    For arguments of any type except `Number` & `AbstractArray`, the result is `nothing`.

```jldoctest; setup=:(using Zygote)
julia> jacobian((a,s) -> a.^length(s), [1,2,3], "str")
([3 0 0; 0 12 0; 0 0 27], nothing)

julia> jacobian((a,t) -> sum(a .* t[1]) + t[2], [1,2,3], (4,5))
([4 4 4], nothing)

julia> gradient((a,t) -> sum(a .* t[1]) + t[2], [1,2,3], (4,5))  # gradient undersands the tuple
([4, 4, 4], (6, 1))
```
"""
function jacobian(f, args...)
  y, back = pullback(_jvec∘f, args...)
  out = map(args) do x
    T = promote_type(eltype(x), eltype(y))
    dx = x isa AbstractArray ? similar(x, T, length(y), length(x)) :
      x isa Number ? similar(y, T, length(y)) :
      nothing
  end
  delta = _eyelike(y)
  for k in LinearIndices(y)
    grads = back(delta[:,k])
    for (dx, grad) in zip(out, grads)
      dx isa AbstractArray || continue
      _gradcopy!(view(dx,k,:), grad)
    end
  end
  out
end

_jvec(x::AbstractArray) = vec(x)
_jvec(x::Number) = _jvec(vcat(x))
_jvec(x) = throw(ArgumentError("jacobian expected a function which returns an array, or a scalar, got $(typeof(x))"))
_jvec(x::AbstractArray{<:Complex}) = throw(ArgumentError("jacobian does not accept complex output"))

_eyelike(y::Vector) = Matrix{eltype(y)}(I, length(y), length(y))
function _eyelike(y::AbstractVector) # version which works on GPU
  out = fill!(similar(y, length(y), length(y)), 0)
  out[LinearAlgebra.diagind(out)] .= 1
  out
end

_gradcopy!(dst::AbstractArray, src::AbstractArray{<:Number}) = copyto!(dst, src)
_gradcopy!(dst::AbstractArray, src::Number) = copyto!(dst, src)
_gradcopy!(dst::AbstractArray, src::Nothing) = dst .= 0
_gradcopy!(dst::AbstractArray, src::AbstractArray) = copyto!(dst, g isa Number ? g : 0 for g in src) # e.g. Union{Nothing,Float64}

"""
    jacobian(loss, ::Params)

Like [`gradient`](@ref) with implicit parameters, this method takes a zero-argument function
and returns an `IdDict`-like object, now containing the Jacobian for each parameter.

# Examples
```jldoctest; setup=:(using Zygote)
julia> xs = [1 2; 3 4]; ys = [5,7,9];

julia> Jxy = jacobian(() -> ys[1:2] .+ sum(xs.^2), Params([xs, ys]))
Grads(...)

julia> Jxy[ys]
2×3 Array{$Int,2}:
 1  0  0
 0  1  0

julia> Jxy[xs]
2×4 Array{$Int,2}:
 2  6  4  8
 2  6  4  8
```
"""
function jacobian(f, pars::Params)
  y, back = pullback(_jvec∘f, pars)
  out = IdDict()
  for p in pars
    T = Base.promote_type(eltype(p), eltype(y))
    J = similar(y, T, length(y), length(p))
    out[p] = J
  end
  delta = _eyelike(y)
  for k in LinearIndices(y)
    grads = back(delta[:,k])
    for p in pars
      out[p] isa AbstractArray || continue
      _gradcopy!(view(out[p],k,:), grads[p])
    end
  end
  Grads(out, pars)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@generated function Carray{T, N, S}(::UndefInitializer) where {T, N, S}
	init = ((zero(UInt8) for _ in 1:S)...,)
	return :(Carray{T, N, S}(undef, $(init)))
end

(::Type{CA})(args...) where {CA<:Carray} = CA(undef, args...)
function (::Type{CA})(init::Union{NTuple{N, UInt8} where {N}, CA, Cbinding{CA}, Cptrs{CA}, UndefInitializer}, args...) where {CA<:Carray}
	result =
		init isa UndefInitializer ? bitstype(CA)(undef) :
		init isa Tuple ? bitstype(CA)(undef, init) :
		init isa Cbinding ? init[] :
		init isa Cptrs ? init[] :
		init

	for (ind, arg) in enumerate(args)
		result = arg isa Pair && first(arg) isa Integer ? initindex!(result, last(arg), first(arg)) : initindex!(result, arg, ind)
	end

	return result
end

Base.zero(::Type{CA}) where {CA<:Carray} = CA()
Base.sizeof(::Type{CA}) where {CA<:Carray} = invoke(sizeof, Tuple{Any}, bitstype(CA))
Base.size(::Type{CA}) where {T, N, CA<:Carrays{T, N}} = (N,)
Base.length(::Type{CA}) where {T, N, CA<:Carrays{T, N}} = N
Base.eltype(::Type{CA}) where {T, N, CA<:Carrays{T, N}} = T
Base.convert(::Type{CA}, t::Tuple) where {CA<:Carray} = CA(t...)
Base.convert(::Type{CA}, v::Vector) where {CA<:Carray} = CA(v...)

Base.convert(::Type{CA}, str::String) where {T<:Union{Int8, UInt8, Cconst{Int8}, Cconst{UInt8}}, CA<:Carray{T}} = CA(str...)
Base.convert(::Type{String}, ca::Carray{T}) where {T<:Union{Int8, UInt8, Cconst{Int8}, Cconst{UInt8}}} = String(map(t -> reinterpret(UInt8, convert(unqualifiedtype(T), t)), collect(ca)))
Base.String(ca::Carray{T}) where {T<:Union{Int8, UInt8, Cconst{Int8}, Cconst{UInt8}}} = convert(String, ca)
Base.show(io::IO, ca::Carray{T}) where {T<:Union{Int8, UInt8, Cconst{Int8}, Cconst{UInt8}}} = show(io, String(ca))

Base.firstindex(ca::CA) where {CA<:Carrays} = 1
Base.lastindex(ca::CA) where {CA<:Carrays} = length(ca)

Base.IndexStyle(::Type{CA}) where {CA<:Carrays} = IndexLinear()
Base.size(ca::CA) where {CA<:Carrays} = size(typeof(ca))
Base.length(ca::CA) where {CA<:Carrays} = length(typeof(ca))
Base.eltype(ca::CA) where {CA<:Carrays} = eltype(typeof(ca))

Base.keys(ca::CA) where {CA<:Carrays} = firstindex(ca):lastindex(ca)
Base.values(ca::CA) where {CA<:Carrays} = iterate(ca)
Base.iterate(ca::CA, i::Integer = 1) where {CA<:Carrays} = firstindex(ca) <= i <= lastindex(ca) ? (ca[i], i+1) : nothing

Base.getindex(ptr::Cptrs{CA}, i::Integer) where {T, N, CA<:Carrays{T, N}} = Core.Intrinsics.bitcast(Cptr{isqualifiedwith(CA, Cconst) && !isqualifiedwith(T, Cconst) ? Cconst{T} : T}, ptr) + (i-1)
Base.getindex(ca::CA, i::Integer) where {T, N, CA<:Carrays{T, N}} = load(ca, field(CA, i))

Base.setindex!(ptr::Cptrs{CA}, val, i::Integer) where {T, N, CA<:Carray{T, N}} = unsafe_store!(Core.Intrinsics.bitcast(Cptr{T}, ptr) + (i-1), val)
initindex!(ca::CA, val, i::Integer) where {T, N, CA<:Carray{T, N}} = store!(ca, field(CA, i), val)