text:
"""
    JavaRef is abstract parent for JavaLocalRef, JavaGlobalRef, and JavaNullRef in the JavaCall Module

    It is distinct from its parent type, JavaCall.JNI.AbstractJavaRef, since its use is defined in
    JavaCall itself rather than the JNI submodule.
"""
abstract type JavaRef <: JNI.AbstractJavaRef end

"""
    JavaLocalRef is a JavaRef that is meant to be used with local variables in a function call.
    After the function call these references may be freed and garbage collected. See note about
    JNI memory management below.

    This is the default reference type returned from the JNI.

    Use this with JNI.PushLocalFrame / JNI.PopLocalFrame for memory management.
    Also see JNI.EnsureLocalCapacity.

    The internal pointer should be deleted using JNI.DeleteLocalRef
"""
struct JavaLocalRef <: JavaRef
    ptr::Ptr{Nothing}
end

"""
    JavaGlobalRef is a JavaRef that is meant to be used with global variables that live beyond
    a single function call.
"""
struct JavaGlobalRef <: JavaRef
    ptr::Ptr{Nothing}
end

"""
    JavaNullRef is a JavaRef that serves as a placeholder to mark where references have already been deleted.

    See J_NULL
"""
struct JavaNullRef <: JavaRef
    ptr::Ptr{Nothing}
    JavaNullRef() = new(C_NULL)
end

""" Constant JavaNullRef """
const J_NULL = JavaNullRef()

Ptr(ref::JavaRef) = ref.ptr
Ptr{Nothing}(ref::JavaRef) = ref.ptr

JavaLocalRef(ref::JavaRef) = JavaLocalRef(JNI.NewLocalRef(Ptr(ref)))
JavaGlobalRef(ref::JavaRef) = JavaGlobalRef(JNI.NewGlobalRef(Ptr(ref)))

# _deleteref does local/global reference deletion without null or state checking
_deleteref(ref::JavaLocalRef ) = JNI.DeleteLocalRef( Ptr(ref))
_deleteref(ref::JavaGlobalRef) = JNI.DeleteGlobalRef(Ptr(ref))
_deleteref(ref::JavaNullRef) = nothing

"""
    deleteref deletes a JavaRef using either JNI.DeleteLocalRef or JNI.DeleteGlobalRef
"""
function deleteref(x::JavaRef)
    if x.ptr == C_NULL; return; end
    if !JNI.is_env_loaded(); return; end;
    _deleteref(x)
    return
end

"""
    jlocalframe(f, [returntype]; [capacity = 16])

    Manages java local references by using JNI's PushLocalFrame and PopLocalFrame.
    Only the local reference returned by f will be valid. Other local references
    will be freed and available for garbage collection.

    Specifying a `returntype` will allow for type stability. If `returntype` is
    specified and is not `Nothing` or `Any`, it will also be passed to the function.

    Capacity specifies the minimum number of local references that can be
    created. See the [JNI documentation](
    https://docs.oracle.com/en/java/javase/15/docs/specs/jni/functions.html#pushlocalframe
    )
    for further information.

    # Example
    ```
    julia> jlocalframe() do
               a = JObject() # Local reference created, will be GCed
               println(a)
               b = JObject() # Local reference returned
               println(b)
               b
           end

    julia> jlocalframe(JObject) do T # Specify returntype for type stability
               a = T()
               println(a)
               b = T()
               println(b)
               b
           end

    julia> jlocalframe(Nothing) do # Specify Nothing if you do want to return anything
               a = JObject()
               println(a)
           end
    ```
"""
function jlocalframe(f::Function, returntype::Type = Any; capacity = 16)
    JNI.PushLocalFrame(jint(capacity))
    result_ref = C_NULL
    return_ref = JavaLocalRef(result_ref)
    result = nothing
    try
        if returntype == Any
            result = f()
        else
            result = f(returntype)
        end
        if isa(result, JavaObject)
            result = Ptr{Nothing}(result)
        end
        if isa(result, Ptr{Nothing}) &&
           JNI.GetObjectRefType(result) == JNI.JNILocalRefType
            result_ref = result
        end
    catch err
        rethrow(err)
    finally
        return_ref = JavaLocalRef( JNI.PopLocalFrame(result_ref) )
    end

    # Return
    if returntype == Any # Not Type Stable
        if !isnull(return_ref.ptr)
            return narrow( JObject(return_ref) )
        else
            return result
        end
    elseif returntype <: JavaObject
        return returntype(return_ref)
    else
        return result::returntype
    end
end

# Closer to https://github.com/ahnlabb/BioformatsLoader.jl/commit/4d4e2d5decd87c8bfd2bfca2fdfbc4214b120977
function jlocalframe(f::Function, returntype::Type{Nothing}; capacity = 16)
    JNI.PushLocalFrame(jint(capacity))
    try
        f()
    catch err
        rethrow(err)
    finally
        JNI.PopLocalFrame(C_NULL)
    end
    return nothing
end

"""
    JavaMetaClass represents meta information about a Java class

    These are usually cached in _jmc_cache and are meant to live
    as long as the cache is valid.
"""
struct JavaMetaClass{T} <: JNI.AbstractJavaRef
    ref::JavaRef
end

#The metaclass, sort of equivalent to a the
JavaMetaClass(T, ref::JavaRef) = JavaMetaClass{T}(ref)
JavaMetaClass(T, ptr::Ptr{Nothing}) = JavaMetaClass{T}(JavaGlobalRef(ptr))

ref(mc::JavaMetaClass{T}) where T = mc.ref
Ptr(mc::JavaMetaClass{T}) where T = Ptr(mc.ref)
Ptr{Nothing}(mc::JavaMetaClass{T}) where T = Ptr(mc.ref)

"""
    JavaObject{T} is the main JavaCall type representing either an instance
    or a static class

    T is usually a symbol referring a Java class name
"""
mutable struct JavaObject{T} <: JNI.AbstractJavaRef
    ref::JavaRef

    #This below is ugly. Once we stop supporting 0.5, this can be replaced by
    # function JavaObject{T}(ptr) where T
    function JavaObject{T}(ref) where T
        j = new{T}(ref)
        finalizer(deleteref, j)
        return j
    end

    #replace with: JavaObject{T}(argtypes::Tuple, args...) where T
    JavaObject{T}(argtypes::Tuple, args...) where {T} = jnew(T, argtypes, args...)
end

# JavaObject Construction
JavaObject(ptr) = JObject(ptr)
JavaObject(T, ptr) = JavaObject{T}(ptr)
JavaObject{T}() where {T} = JavaObject{T}((),)
JavaObject{T}(ptr::Ptr{Nothing}) where {T} = JavaObject{T}(JavaLocalRef(ptr))

# JavaObject Reference Management
ref(x::JavaObject{T}) where T = x.ref
copyref(x::JavaObject{T}) where T = JavaObject{T}(JavaLocalRef(x.ref))
deleteref(x::JavaObject{T}) where T = ( deleteref(x.ref); x.ref = J_NULL )

# Obtain the underlying pointer for a JavaObject
Ptr(x::JavaObject{T}) where T = Ptr(x.ref)
Ptr{Nothing}(x::JavaObject{T}) where T = Ptr(x.ref)

"""
   jglobal(x::JavaObject) creates a new JavaGlobalRef and deletes the prior JavaRef
"""
function jglobal(x::JavaObject)
    gref = JavaGlobalRef(JNI.NewGlobalRef(Ptr(x)))
    deleteref(x.ref)
    x.ref = gref
end

"""
```
isnull(obj::JavaObject)
```
Checks if the passed JavaObject is null or not

### Args
* obj: The object of type JavaObject

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaObject) = Ptr(obj) == C_NULL
isnull(obj::Ptr{Nothing}) = obj == C_NULL

"""
```
isnull(obj::JavaMetaClass)
```
Checks if the passed JavaMetaClass is null or not

### Args
* obj: The object of type JavaMetaClass

### Returns
true if the passed object is null else false
"""
isnull(obj::JavaMetaClass) = Ptr(obj) == C_NULL

const JClass = JavaObject{Symbol("java.lang.Class")}
const JObject = JavaObject{Symbol("java.lang.Object")}
const JMethod = JavaObject{Symbol("java.lang.reflect.Method")}
const JConstructor = JavaObject{Symbol("java.lang.reflect.Constructor")}
const JField = JavaObject{Symbol("java.lang.reflect.Field")}
const JThread = JavaObject{Symbol("java.lang.Thread")}
const JClassLoader = JavaObject{Symbol("java.lang.ClassLoader")}
const JString = JavaObject{Symbol("java.lang.String")}

#JavaObject(ptr::Ptr{Nothing}) = ptr == C_NULL ? JavaObject(ptr) : JavaObject{Symbol(getclassname(getclass(ptr)))}(ptr)

function JString(str::AbstractString)
    jstring = JNI.NewStringUTF(String(str))
    if jstring == C_NULL
        geterror()
    else
        return JString(jstring)
    end
end

# jvalue(v::Integer) = int64(v) << (64-8*sizeof(v))
jvalue(v::Integer)::JNI.jvalue = JNI.jvalue(v)
jvalue(v::Float32) = jvalue(reinterpret(Int32, v))
jvalue(v::Float64) = jvalue(reinterpret(Int64, v))
jvalue(v::Ptr) = jvalue(Int(v))
jvalue(v::JavaObject) = jvalue(Ptr(v))

function _jimport(juliaclass)
    for str ∈ [" ", "(", ")"]
        juliaclass = replace(juliaclass, str=>"")
    end
    :(JavaObject{Symbol($juliaclass)})
end

macro jimport(class::Expr)
    juliaclass = sprint(Base.show_unquoted, class)
    _jimport(juliaclass)
end
macro jimport(class::Symbol)
    juliaclass = string(class)
    _jimport(juliaclass)
end
macro jimport(class::AbstractString)
    _jimport(class)
end

const primitive_names_to_types = Dict(
    :boolean => jboolean,
    :byte    => jbyte,
    :char    => jchar,
    :short   => jshort,
    :int     => jint,
    :long    => jlong,
    :float   => jfloat,
    :double  => jdouble,
    :void    => jvoid
)
jimport(juliaclass::Symbol) = juliaclass == :void ? Nothing :
    haskey(primitive_names_to_types, juliaclass) ? jimport(juliaclass, Val(true), Val(false)) : JavaObject{juliaclass}
jimport(juliaclass::Symbol, isprimitive::Val{false}, isarray::Val{false}) = jimport(juliaclass)
jimport(juliaclass::Symbol, isprimitive::Val{true},  isarray::Val{false}) = primitive_names_to_types[juliaclass]

jimport(juliaclass::String, args...) = isarray(juliaclass) ? Vector{ jimport(Symbol(juliaclass[1:end-2])) } : jimport(Symbol(juliaclass), args...)

function jimport(juliaclass::JClass)
    jimport(juliaclass, Val(isprimitive(juliaclass)), Val(isarray(juliaclass)))
end
function jimport(juliaclass::JClass, isprimitive, isarray::Val{true})
    elementType = jimport( jcall(juliaclass, "getComponentType", JClass) )
    Vector{elementType}
end
jimport(juliaclass::JClass, isprimitive, isarray::Val{false}) = jimport(getname(juliaclass), isprimitive, isarray)

isprimitive(juliaclass::JClass) = jcall(juliaclass, "isPrimitive", jboolean, ()) == 0x01
isarray(juliaclass::JClass) = jcall(juliaclass, "isArray", jboolean, ()) == 0x01
isarray(juliaclass::String) = endswith(juliaclass, "[]")

function jnew(T::Symbol, argtypes::Tuple = () , args...)
    assertroottask_or_goodenv() && assertloaded()
    sig = method_signature(Nothing, argtypes...)
    jmethodId = JNI.GetMethodID(Ptr(metaclass(T)), String("<init>"), sig)
    if jmethodId == C_NULL
        throw(JavaCallError("No constructor for $T with signature $sig"))
    end
    return  _jcall(metaclass(T), jmethodId, JNI.NewObjectA, JavaObject{T}, argtypes, args...)
end

# Call static methods
function jcall(typ::Type{JavaObject{T}}, method::AbstractString, rettype::Type, argtypes::Tuple = (),
               args... ) where T
    assertroottask_or_goodenv() && assertloaded()
    sig = method_signature(rettype, argtypes...)
    jmethodId = JNI.GetStaticMethodID(Ptr(metaclass(T)), String(method), sig)
    jmethodId==C_NULL && geterror(true)
    _jcall(metaclass(T), jmethodId, C_NULL, rettype, argtypes, args...)
end

function jcall(typ::Type{JavaObject{T}}, method::JMethod, args...) where T
    assertroottask_or_goodenv() && assertloaded()
    jmethodId = JNI.FromReflectedMethod(method)
    rettype = jimport(getreturntype(method))
    argtypes = Tuple(jimport.(getparametertypes(method)))
    jmethodId==C_NULL && geterror(true)
    _jcall(metaclass(T), jmethodId, C_NULL, rettype, argtypes, args...)
end

# Call instance methods
function jcall(obj::JavaObject, method::AbstractString, rettype::Type, argtypes::Tuple = (), args... )
    assertroottask_or_goodenv() && assertloaded()
    sig = method_signature(rettype, argtypes...)
    jmethodId = JNI.GetMethodID(Ptr(metaclass(obj)), String(method), sig)
    jmethodId==C_NULL && geterror(true)
    _jcall(obj, jmethodId, C_NULL, rettype,  argtypes, args...)
end

function jcall(obj::JavaObject, method::JMethod, args... )
    assertroottask_or_goodenv() && assertloaded()
    jmethodId = JNI.FromReflectedMethod(method)
    rettype = jimport(getreturntype(method))
    argtypes = Tuple(jimport.(getparametertypes(method)))
    jmethodId==C_NULL && geterror(true)
    _jcall(obj, jmethodId, C_NULL, rettype,  argtypes, args...)
end

# JMethod invoke
(m::JMethod)(obj, args...) = jcall(obj, m, args...)

function jfield(typ::Type{JavaObject{T}}, field::AbstractString, fieldType::Type) where T
    assertroottask_or_goodenv() && assertloaded()
    jfieldID  = JNI.GetStaticFieldID(Ptr(metaclass(T)), String(field), signature(fieldType))
    jfieldID==C_NULL && geterror(true)
    _jfield(metaclass(T), jfieldID, fieldType)
end

function jfield(obj::Type{JavaObject{T}}, field::JField) where T
    assertroottask_or_goodenv() && assertloaded()
    fieldType = jimport(gettype(field))
    jfieldID = JNI.FromReflectedField(field)
    jfieldID==C_NULL && geterror(true)
    _jfield(metaclass(T), jfieldID, fieldType)
end

function jfield(obj::JavaObject, field::AbstractString, fieldType::Type)
    assertroottask_or_goodenv() && assertloaded()
    jfieldID  = JNI.GetFieldID(Ptr(metaclass(obj)), String(field), signature(fieldType))
    jfieldID==C_NULL && geterror(true)
    _jfield(obj, jfieldID, fieldType)
end

function jfield(obj::JavaObject, field::JField)
    assertroottask_or_goodenv() && assertloaded()
    fieldType = jimport(gettype(field))
    jfieldID = JNI.FromReflectedField(field)
    jfieldID==C_NULL && geterror(true)
    _jfield(obj, jfieldID, fieldType)
end

# JField invoke
(f::JField)(obj) = jfield(obj, f)

for (x, y, z) in [(:jboolean, :(JNI.GetBooleanField), :(JNI.GetStaticBooleanField)),
                  (:jchar,    :(JNI.GetCharField),    :(JNI.GetStaticCharField))   ,
                  (:jbyte,    :(JNI.GetByteField),    :(JNI.GetStaticBypeField))   ,
                  (:jshort,   :(JNI.GetShortField),   :(JNI.GetStaticShortField))  ,
                  (:jint,     :(JNI.GetIntField),     :(JNI.GetStaticIntField))    ,
                  (:jlong,    :(JNI.GetLongField),    :(JNI.GetStaticLongField))   ,
                  (:jfloat,   :(JNI.GetFloatField),   :(JNI.GetStaticFloatField))  ,
                  (:jdouble,  :(JNI.GetDoubleField),  :(JNI.GetStaticDoubleField)) ]

    m = quote
        function _jfield(obj, jfieldID::Ptr{Nothing}, fieldType::Type{$(x)})
            callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            result = callmethod(Ptr(obj), jfieldID)
            result==C_NULL && geterror()
            return convert_result(fieldType, result)
        end
    end
    eval(m)
end

function _jfield(obj, jfieldID::Ptr{Nothing}, fieldType::Type)
    callmethod = ifelse( typeof(obj)<:JavaObject, JNI.GetObjectField , JNI.GetStaticObjectField )
    result = callmethod(Ptr(obj), jfieldID)
    result==C_NULL && geterror()
    return convert_result(fieldType, result)
end

#Generate these methods to satisfy ccall's compile time constant requirement
#_jcall for primitive and Nothing return types
for (x, y, z) in [(:jboolean, :(JNI.CallBooleanMethodA), :(JNI.CallStaticBooleanMethodA)),
                  (:jchar,    :(JNI.CallCharMethodA),    :(JNI.CallStaticCharMethodA))   ,
                  (:jbyte,    :(JNI.CallByteMethodA),    :(JNI.CallStaticByteMethodA))   ,
                  (:jshort,   :(JNI.CallShortMethodA),   :(JNI.CallStaticShortMethodA))  ,
                  (:jint,     :(JNI.CallIntMethodA),     :(JNI.CallStaticIntMethodA))    ,
                  (:jlong,    :(JNI.CallLongMethodA),    :(JNI.CallStaticLongMethodA))   ,
                  (:jfloat,   :(JNI.CallFloatMethodA),   :(JNI.CallStaticFloatMethodA))  ,
                  (:jdouble,  :(JNI.CallDoubleMethodA),  :(JNI.CallStaticDoubleMethodA)) ,
                  (:jvoid,    :(JNI.CallVoidMethodA),    :(JNI.CallStaticVoidMethodA))   ]
    m = quote
        function _jcall(obj, jmethodId::Ptr{Nothing}, callmethod::Ptr{Nothing}, rettype::Type{$(x)},
                        argtypes::Tuple, args...)
            if callmethod == C_NULL #!
                callmethod = ifelse( typeof(obj)<:JavaObject, $y , $z )
            end
            @assert callmethod != C_NULL
            @assert jmethodId != C_NULL
            isnull(obj) && throw(JavaCallError("Attempt to call method on Java NULL"))
            savedArgs, convertedArgs = convert_args(argtypes, args...)
            GC.@preserve savedArgs begin
                result = callmethod(Ptr(obj), jmethodId, Array{JNI.jvalue}(jvalue.(convertedArgs)))
            end
            deleteref.(filter(x->isa(x,JavaObject),convertedArgs))
            result==C_NULL && geterror()
            result == nothing && (return)
            return convert_result(rettype, result)
        end
    end
    eval(m)
end

#_jcall for Object return types
#obj -- receiver - Class pointer or object prointer
#jmethodId -- Java method ID
#callmethod -- the C method pointer to call
function _jcall(obj, jmethodId::Ptr{Nothing}, callmethod::Union{Function,Ptr{Nothing}}, rettype::Type, argtypes::Tuple,
                args...)
    if callmethod == C_NULL
        callmethod = ifelse(typeof(obj)<:JavaObject,
                            JNI.CallObjectMethodA  ,
                            JNI.CallStaticObjectMethodA)
    end
    @assert callmethod != C_NULL
    @assert jmethodId != C_NULL
    isnull(obj) && error("Attempt to call method on Java NULL")
    savedArgs, convertedArgs = convert_args(argtypes, args...)
    GC.@preserve savedArgs begin
        result = callmethod(Ptr(obj), jmethodId, Array{JNI.jvalue}(jvalue.(convertedArgs)))
    end
    deleteref.(filter(x->isa(x,JavaObject),convertedArgs))
    result==C_NULL && geterror()
    return convert_result(rettype, result)
end

global const _jmc_cache = [ Dict{Symbol, JavaMetaClass}() ]

function _metaclass(class::Symbol)
    jclass=javaclassname(class)
    jclassptr = JNI.FindClass(jclass)
    jclassptr == C_NULL && throw(JavaCallError("Class Not Found $jclass"))
    return JavaMetaClass(class, jclassptr)
end

function metaclass(class::Symbol)
    return _metaclass(class)
    if !haskey(_jmc_cache[ Threads.threadid() ], class)
        _jmc_cache[ Threads.threadid() ][class] = _metaclass(class)
    end
    return _jmc_cache[ Threads.threadid() ][class]
end

metaclass(::Type{JavaObject{T}}) where {T} = metaclass(T)
metaclass(::JavaObject{T}) where {T} = metaclass(T)
metaclass(::Type{T}) where T <: AbstractVector = metaclass( Symbol( JavaCall.signature(T) ) )

javaclassname(class::Symbol) = replace(string(class), "."=>"/")
javaclassname(class::AbstractString) = replace(class, "."=>"/")
javaclassname(::Type{T}) where T <: AbstractVector = JavaCall.signature(T)

function geterror(allow=false)
    isexception = JNI.ExceptionCheck()

    if isexception == JNI_TRUE
        jthrow = JNI.ExceptionOccurred()
        jthrow==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        JNI.ExceptionDescribe() #Print java stackstrace to stdout
        JNI.ExceptionClear()
        jclass = JNI.FindClass("java/lang/Throwable")
        jclass==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        jmethodId=JNI.GetMethodID(jclass, "toString", "()Ljava/lang/String;")
        jmethodId==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        res = JNI.CallObjectMethodA(jthrow, jmethodId, Int[])
        res==C_NULL && throw(JavaCallError("Java Exception thrown, but no details could be retrieved from the JVM"))
        msg = unsafe_string(JString(res))
        JNI.DeleteLocalRef(jthrow)
        throw(JavaCallError(string("Error calling Java: ",msg)))
    else
        if allow==false
            return #No exception pending, legitimate NULL returned from Java
        else
            throw(JavaCallError("Null from Java. Not known how"))
        end
    end
end

#get the JNI signature string for a method, given its
#return type and argument types
function method_signature(rettype, argtypes...)
    s=IOBuffer()
    write(s, "(")
    for arg in argtypes
        write(s, signature(arg))
    end
    write(s, ")")
    write(s, signature(rettype))
    return String(take!(s))
end

#get the JNI signature string for a given type
function signature(arg::Type)
    if arg === jboolean
        return "Z"
    elseif arg === jbyte
        return "B"
    elseif arg === jchar
        return "C"
    elseif arg === jshort
        return "S"
    elseif arg === jint
        return "I"
    elseif arg === jlong
        return "J"
    elseif arg === jfloat
        return "F"
    elseif arg === jdouble
        return "D"
    elseif arg === jvoid
        return "V"
    elseif arg <: Array
        dims = "[" ^ ndims(arg)
        return string(dims, signature(eltype(arg)))
    end
end

signature(arg::Type{JavaObject{T}}) where {T} = string("L", javaclassname(T), ";")
signature(arg::Type{JavaObject{T}}) where {T <: AbstractVector} = JavaCall.javaclassname(T)


text:
#=

Lower bound instance from

http://proceedings.mlr.press/v28/jaggi13.pdf

and

https://arxiv.org/abs/1309.5550

Example instance is to minimize || x ||^2 over the probability simplex conv(e_1, ..., e_n)

Then the primal gap is lower bounded by 1/k - 1/n in iteration k as the optimal solution has value 1/n attained by
the (1/n, ..., 1/n) vector and in iteration k we have picked up at most k vertices from the simplex lower bounding the
primal value by 1/k.

Here: slightly rewritten to consider || x - (1/n, ..., 1/n) ||^2 so that the function value becomes directly the
primal gap (squared)

Three runs are compared:
1. Frank-Wolfe with traditional step-size rule
2. Away-step Frank-Wolfe with adaptive step-size rule
3. Blended Conditional Gradients with adaptive step-size rule

NOTE:
1. ignore the timing graphs
2. as the primal gap lower bounds the dual gap we also plot the primal gap lower bound in the dual gap graph
3. AFW violates the bound in the very first round which is due to the initialization in AFW starting with two vertices
4. all methods that call an LMO at most once in an iteration are subject to this lower bound
5. the objective is strongly convex, this implies limitations also for strongly convex functions (see for a discussion https://arxiv.org/abs/1906.07867)

=#

include("activate.jl")

import LinearAlgebra

# n = Int(1e1)
n = Int(1e2)
k = Int(1e3)

xp = 1 / n * ones(n);

# definition of objective
f(x) = LinearAlgebra.norm(x - xp)^2

# definition of gradient
function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

# define LMO and do initial call to obtain starting point
lmo = FrankWolfe.ProbabilitySimplexOracle(1)
x0 = FrankWolfe.compute_extreme_point(lmo, zeros(n));

# simple benchmarking of oracles to get an idea how expensive each component is
FrankWolfe.benchmark_oracles(f, grad!, () -> rand(n), lmo; k=100)

@time x, v, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    epsilon=1e-5,
    trajectory=true,
);

@time x, v, primal, dual_gap, trajectoryAFW = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    epsilon=1e-5,
    trajectory=true,
);

@time x, v, primal, dual_gap, trajectoryBCG = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    epsilon=1e-5,
    trajectory=true,
);

# define lower bound

trajLowerbound = []
for i in 1:n-1
    push!(trajLowerbound, (i, 1 / i - 1 / n, NaN, 1 / i - 1 / n, NaN))
end

data = [trajectory, trajectoryAFW, trajectoryBCG, trajLowerbound]
label = ["FW", "AFW", "BCG", "Lowerbound"]

# ignore the timing plots - they are not relevant for this example
FrankWolfe.plot_trajectories(data, label, xscalelog=true, legend_position=:bottomleft)

ext: jl


text:
# utils for setting or finding colors

function colorwheel(n::Int)
  # RGB(1.0, 1.0, 0)
  convert(RGB, HSV((n*30)%360, 1.0,0.5))
end

function pointToColor(nm::Symbol)
  if nm == :PartialPose3XYYaw
    return RGBA(0.6,0.8,0.9,0.5)
  elseif nm == :Pose3Pose3NH
    return RGBA(1.0,1.0,0,0.5)
  else
    # println("pointToColor(..) -- warning, using default color for edges")
    return RGBA(0.0,1,0.0,0.5)
  end
end

function submapcolor(idx::Int, len::Int;
        submapcolors=SubmapColorCheat() )
  #
  n = idx%length(submapcolors.colors)+1
  smc = submapcolors.colors[n]
  return [smc for g in 1:len]
end


text:
import Base: eltype, size, length, getindex, setindex!, similar, copy

import SPTK

### Generic interface ###

abstract type SpectralParam end

function estimate(s::SpectralParam, x::AbstractArray)
    error("should provide an estimator")
end

function params(s::SpectralParam)
    error("should provide get access to parameters")
end

# subtypes of SpectralParam should have field `order`
param_order(s::SpectralParam) = s.order

### Types that characterize spectral parameters ###

abstract type FrequencyForm end
struct MelFrequency <: FrequencyForm
end
struct LinearFrequency <: FrequencyForm
end

abstract type LogForm end
struct GeneralizedLog <: LogForm
end
struct StandardLog <: LogForm
end
struct AllPoleLog <: LogForm
end

### Definitions of spectral envelope parameters ###

struct MelGeneralizedCepstrum{F<:FrequencyForm,L<:LogForm} <: SpectralParam
    order::Int
    α::Real
    γ::Real
    function (::Type{MelGeneralizedCepstrum{F,L}}){F<:FrequencyForm,L<:LogForm}(
            order::Int, α::Real, γ::Real)
        order > 0 || throw(ArgumentError("order should be larger than 0"))
        abs(α) < 1 || throw(ArgumentError("|α| < 1 is supported"))
        (-1 <= γ <= 0) || throw(ArgumentError("-1 <= γ <= 0 is supported"))
        new{F,L}(order, α, γ)
    end
end

const MelFrequencyCepstrum{L} = MelGeneralizedCepstrum{MelFrequency,L}
const LinearFrequencyCepstrum{L} = MelGeneralizedCepstrum{LinearFrequency,L}
const GeneralizedLogCepstrum{F} = MelGeneralizedCepstrum{F,GeneralizedLog}
const StandardLogCepstrum{F} = MelGeneralizedCepstrum{F,StandardLog}
const AllPoleLogCepstrum{F} = MelGeneralizedCepstrum{F,AllPoleLog}

const GeneralizedCepstrum = MelGeneralizedCepstrum{LinearFrequency,GeneralizedLog}
const MelCepstrum = MelGeneralizedCepstrum{MelFrequency,StandardLog}
const LinearCepstrum = MelGeneralizedCepstrum{LinearFrequency,StandardLog}
const AllPoleCepstrum = MelGeneralizedCepstrum{LinearFrequency,AllPoleLog}
const MelAllPoleCepstrum = MelGeneralizedCepstrum{MelFrequency,AllPoleLog}

freq_form{F<:FrequencyForm,L<:LogForm}(::Type{MelGeneralizedCepstrum{F,L}}) = F
freq_form{T<:MelGeneralizedCepstrum}(::Type{T}) = freq_form(super(T))

log_form{F<:FrequencyForm,L<:LogForm}(::Type{MelGeneralizedCepstrum{F,L}}) = L
log_form{T<:MelGeneralizedCepstrum}(::Type{T}) = freq_form(super(T))

# Generic fallback constructor, which determines type parmaters from α and γ
function MelGeneralizedCepstrum(order::Int, α::Real, γ::Real)
    F = (α == zero(α)) ? LinearFrequency : MelFrequency

    L = GeneralizedLog
    if γ == zero(γ)
        L = StandardLog
    elseif γ == -one(γ)
        L = AllPoleLog
    end

    MelGeneralizedCepstrum{F,L}(order, α, γ)
end

function MelGeneralizedCepstrum(mgc::MelGeneralizedCepstrum)
    order = param_order(mgc)
    α = allpass_alpha(mgc)
    γ = glog_gamma(mgc)
    MelGeneralizedCepstrum(order, α, γ)
end

function (::Type{GeneralizedCepstrum})(order::Int, γ::Real)
    MelGeneralizedCepstrum{LinearFrequency,GeneralizedLog}(order, 0.0, γ)
end

function (::Type{MelCepstrum})(order::Int, α::Real)
    MelGeneralizedCepstrum{MelFrequency,StandardLog}(order, α, 0.0)
end

function (::Type{LinearCepstrum})(order::Int)
    MelGeneralizedCepstrum{LinearFrequency,StandardLog}(order, 0.0, 0.0)
end

function (::Type{AllPoleCepstrum})(order::Int)
    MelGeneralizedCepstrum{LinearFrequency,AllPoleLog}(order, 0.0, -1.0)
end

function (::Type{MelAllPoleCepstrum})(order::Int, α::Real)
    MelGeneralizedCepstrum{MelFrequency,AllPoleLog}(order, α, -1.0)
end

params(c::MelGeneralizedCepstrum) = (c.order, c.α, c.γ)
params(c::GeneralizedCepstrum) = (c.order, c.γ)
params(c::MelCepstrum) = (c.order, c.α)
params(c::LinearCepstrum) = (c.order,)
params(c::AllPoleCepstrum) = (c.order,)
params(c::MelAllPoleCepstrum) = (c.order, c.α)
allpass_alpha(c::MelGeneralizedCepstrum) = c.α
glog_gamma(c::MelGeneralizedCepstrum) = c.γ

### Definitions of LPC variants ###

abstract type LinearPredictionCoefVariants <: SpectralParam end

# LPC, LSP and PARCOR
struct LinearPredictionCoef <: LinearPredictionCoefVariants
    order
end

struct LineSpectralPair <: LinearPredictionCoefVariants
    order
end

struct PartialAutoCorrelation <: LinearPredictionCoefVariants
    order
end

params(s::LinearPredictionCoefVariants) = (s.order,)

### State, which keeps actual computation results in arrays ###

abstract type AbstractParamState{T,N} <: AbstractArray{T,N} end

function paramdef(s::AbstractParamState)
    error("should provide get access for parameter definition")
end

function rawdata(s::AbstractParamState)
    error("should provide get access for raw data")
end

mutable struct SpectralParamState{S<:SpectralParam,T,N} <: AbstractParamState{T,N}
    def::S
    data::Array{T,N}

    has_loggain::Bool
    gain_normalized::Bool
    ready_to_filt::Bool  # state is explicitly converted to filter coef or not.

    function (::Type{SpectralParamState{S,T,N}}){S,T,N}(s::S, data::Array{T,N},
                                has_loggain::Bool,
                                gain_normalized::Bool,
                                ready_to_filt::Bool)
        new{S,T,N}(s, data, has_loggain, gain_normalized, ready_to_filt)
    end
end

const SpectralParamStateVector{S,T} = SpectralParamState{S,T,1}
const SpectralParamStateMatrix{S,T} = SpectralParamState{S,T,2}

function SpectralParamState{S<:SpectralParam,T,N}(s::S, data::Array{T,N},
                                                  has_loggain::Bool,
                                                  gain_normalized::Bool;
                                                  ready_to_filt::Bool=false)
    SpectralParamState{S,T,N}(s, data, has_loggain, gain_normalized,
                              ready_to_filt)
end

eltype(s::SpectralParamState) = eltype(s.data)
size(s::SpectralParamState) = size(s.data)
length(s::SpectralParamState) = length(s.data)
getindex(s::SpectralParamState, i::Real) = getindex(s.data, i)
getindex(s::SpectralParamState, i::Real...) = getindex(s.data, i...)

# colon indexing like s[:,1]
function getindex(s::SpectralParamState, i::Colon, j...)
    SpectralParamState(paramdef(s), getindex(s.data, i, j...),
                       has_loggain(s), gain_normalized(s),
                       ready_to_filt=ready_to_filt(s))
end

setindex!(s::SpectralParamState, i::Real) = setindex!(s.data, i)
setindex!(s::SpectralParamState, i::Real...) = setindex!(s.data, i...)
setindex!(s::SpectralParamState, i::Colon, j...) = setindex!(s.data, i, j...)

paramdef(s::SpectralParamState) = s.def
rawdata(s::SpectralParamState) = s.data
has_loggain(s::SpectralParamState) = s.has_loggain
gain_normalized(s::SpectralParamState) = s.gain_normalized
ready_to_filt(s::SpectralParamState) = s.ready_to_filt

function loggain!(s::SpectralParamStateVector)
    has_loggain(s) && return s
    data = rawdata(s)
    data[1] = log(data[1])
    s.has_loggain = true
    s
end

function loggain!(s::SpectralParamStateMatrix)
    has_loggain(s) && return s
    data = rawdata(s)
    for i in 1:size(s, 2)
        data[1,i] = log(data[1,i])
    end
    s.has_loggain = true
    s
end

function unloggain!(s::SpectralParamStateVector)
    !has_loggain(s) && return s
    data = rawdata(s)
    data[1] = exp(data[1])
    s.has_loggain = false
    s
end

function unloggain!(s::SpectralParamStateMatrix)
    !has_loggain(s) && return s
    data = rawdata(s)
    for i in 1:size(s, 2)
        data[1,i] = exp(data[1,i])
    end
    s.has_loggain = false
    s
end

function similar(s::SpectralParamState)
    def = paramdef(s)
    data = rawdata(s)
    SpectralParamState(def, similar(data), has_loggain(s), gain_normalized(s),
                       ready_to_filt=ready_to_filt(s))
end

function copy(s::SpectralParamState)
    def = paramdef(s)
    data = rawdata(s)
    SpectralParamState(def, copy(data), has_loggain(s), gain_normalized(s),
                       ready_to_filt=ready_to_filt(s))
end

# Force re-type mel-generalized cepstrums by its parameter values
# needed?
function retype{T<:MelGeneralizedCepstrum}(s::SpectralParamState{T})
    def = paramdef(s)
    newdef = MelGeneralizedCepstrum(def)
    SpectralParamState(newdef, rawdata(s), has_loggain(s), gain_normalized(s),
                       ready_to_filt=ready_to_filt(s))
end

### Asserts ###

function assert_ready_to_filt(state::SpectralParamState)
    ready_to_filt(state) && return nothing
    throw(ArgumentError("""ready_to_filt = $(ready_to_filt(state))
                        filter coefficient input is expected"""))
end

function assert_not_ready_to_filt(state::SpectralParamState)
    !ready_to_filt(state) && return nothing
    throw(ArgumentError("""ready_to_filt = $(ready_to_filt(state))
                        filter coefficient input (can be obtain by `mgc2b` or `mc2b`) is not allowed"""))
end

function assert_gain_normalized(state::SpectralParamState)
    gain_normalized(state) && return nothing
    throw(ArgumentError("""gain_normalized: $(gain_normalized(state))
                        gain is assumed to be normalized"""))
end

function assert_gain_unnormalized(state::SpectralParamState)
    gain_normalized(state) || return nothing
    throw(ArgumentError("""gain_normalized: $(gain_normalized(state))
                        gain is assumed to be un-normalized"""))
end


text:
#!/usr/bin/env julia
using Plots, LaTeXStrings, Measures;
pyplot()
f(x, y) = x^2 + y^2
f0(x) = f(x, 0)
f2(x) = f(x, 2)

xVals, yVals = -5:0.1:5, -5:0.1:5
plot(
  xVals,
  [f0.(xVals), f2.(xVals)],
  c = [:blue :red],
  xlims = (-5, 5),
  legend = :top,
  ylims = (-5, 25),
  ylabel = L"f(x,\cdot)",
  label = [L"f(x,0)" L"f(x,2)"],
)
p1 = annotate!(0, -0.2, text("(0,0) The minimum\n of f(x,0)", :left, :top, 10))

z = [f(x, y) for y in yVals, x in xVals]
p2 = surface(
  xVals,
  yVals,
  z,
  c = cgrad([:blue, :red]),
  legend = :none,
  ylabel = "y",
  zlabel = L"f(x,y)",
)
M = z[1:10, 1:10]
p3 = heatmap(
  M,
  c = cgrad([:blue, :red]),
  yflip = true,
  ylabel = "y",
  xticks = ([1:10;], xVals),
  yticks = ([1:10;], yVals),
)

plot(
  p1,
  p2,
  p3,
  layout = (1, 3),
  size = (1200, 400),
  xlabel = "x",
  margin = 5mm,
)
#savefig(fig1, "basicplot.pdf")


text:
using Parameters

Parameters.@with_kw mutable struct hyp
    """
    A very simple hypergraph composite type, designed to hold a node list N, an edge list E, a degree sequence D,
    """
    N::Vector{Int64}
    E::Dict{Int64, Dict}
    D::Array{Int64, 1} = Array{Int64, 1}()
end

function read_hypergraph_data(dataname::String, maxsize::Int64=25)
    labels = Int64[]
    open("../data/$dataname/node-labels-$dataname.txt") do f
        for line in eachline(f)
            push!(labels, parse(Int64, line))
        end
    end
    n = length(labels)

    E = Dict{Integer, Dict}()
    open("../data/$dataname/hyperedges-$dataname.txt") do f
        for line in eachline(f)
            edge = [parse(Int64, v) for v in split(line, ',')]
            sort!(edge)
            if length(edge) > maxsize; continue; end
            sz = length(edge)
            if !haskey(E, sz)
                E[sz] = Dict{}()
            end
            E[sz][edge] = 1
        end
    end

    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end

    N = 1:n
    H = hyp(N, E, D)

    Hin, weights = hypergraph2incidence(H)
    return Hin, weights, labels
end

max_stars_repo_name: nveldt/HypergraphHomophily

max_issues_repo_name: nveldt/HypergraphHomophily

max_forks_repo_name: nveldt/HypergraphHomophily


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

const fileIdentifierLength = 4

"""
Scalar
A Union of the Julia types `T <: Number` that are allowed in FlatBuffers schema
"""
const Scalar = Union{Bool,
Int8, Int16, Int32, Int64,
UInt8, UInt16, UInt32, UInt64,
Float32, Float64, Enum}

"""
Builder is a state machine for creating FlatBuffer objects.
Use a Builder to construct object(s) starting from leaf nodes.

A Builder constructs byte buffers in a last-first manner for simplicity and
performance.
"""
mutable struct Builder
    bytes::Vector{UInt8}
    minalign::Int
    vtable::Vector{UOffsetT}
    objectend::UOffsetT
    vtables::Vector{UOffsetT}
    head::UOffsetT
    nested::Bool
    finished::Bool
    sharedstrings::Dict{String, UOffsetT}
end

bytes(b::Builder) = getfield(b, :bytes)

Builder(size=0) = Builder(zeros(UInt8, size), 1, UOffsetT[], UOffsetT(0), UOffsetT[], UOffsetT(size), false, false, Dict{String, UOffsetT}())

function reset!(b::Builder)
    empty!(b.bytes)
    empty!(b.vtable)
    emtpy!(b.vtables)
    empty!(b.sharedstrings)
    b.minalign = 1
    b.nested = false
    b.finished = false
    b.head = 0
    return
end

Base.write(sink::Builder, o, x::Union{Bool,UInt8}) = sink.bytes[o+1] = UInt8(x)
function Base.write(sink::Builder, off, x::T) where {T}
    off += 1
    for (i, ind) = enumerate(off:(off + sizeof(T) - 1))
        sink.bytes[ind] = (x >> ((i-1) * 8)) % UInt8
    end
end
Base.write(b::Builder, o, x::Float32) = write(b, o, reinterpret(UInt32, x))
Base.write(b::Builder, o, x::Float64) = write(b, o, reinterpret(UInt64, x))
Base.write(b::Builder, o, x::Enum) = write(b, o, basetype(x)(x))

"""
`finishedbytes` returns a pointer to the written data in the byte buffer.
Panics if the builder is not in a finished state (which is caused by calling
`finish!()`).
"""
function finishedbytes(b::Builder)
    assertfinished(b)
    return view(b.bytes, (b.head + 1):length(b.bytes))
end

function startobject!(b::Builder, numfields)
    assertnotnested(b)
    b.nested = true
    resize!(b.vtable, numfields)
    fill!(b.vtable, 0)
    b.objectend = offset(b)
    return
end

"""
WriteVtable serializes the vtable for the current object, if applicable.

Before writing out the vtable, this checks pre-existing vtables for equality
to this one. If an equal vtable is found, point the object to the existing
vtable and return.

Because vtable values are sensitive to alignment of object data, not all
logically-equal vtables will be deduplicated.

A vtable has the following format:
<VOffsetT: size of the vtable in bytes, including this value>
<VOffsetT: size of the object in bytes, including the vtable offset>
<VOffsetT: offset for a field> * N, where N is the number of fields in
the schema for this type. Includes deprecated fields.
Thus, a vtable is made of 2 + N elements, each SizeVOffsetT bytes wide.

An object has the following format:
<SOffsetT: offset to this object's vtable (may be negative)>
<byte: data>+
"""
function writevtable!(b::Builder)
    # Prepend a zero scalar to the object. Later in this function we'll
    # write an offset here that points to the object's vtable:
    prepend!(b, SOffsetT(0))

    objectOffset = offset(b)
    existingVtable = UOffsetT(0)

    # Trim vtable of trailing zeroes.
    i = findlast(!iszero, b.vtable)
    if i !== nothing
        resize!(b.vtable, i)
    end

    # Search backwards through existing vtables, because similar vtables
    # are likely to have been recently appended. See
    # BenchmarkVtableDeduplication for a case in which this heuristic
    # saves about 30% of the time used in writing objects with duplicate
    # tables.
    for i = length(b.vtables):-1:1
        # Find the other vtable, which is associated with `i`:
        vt2Offset = b.vtables[i]
        vt2Start = length(b.bytes) - vt2Offset
        vt2Len = readbuffer(b.bytes, vt2Start, VOffsetT)

        metadata = VtableMetadataFields * sizeof(VOffsetT)
        vt2End = vt2Start + vt2Len
        vt2 = view(b.bytes, (vt2Start + metadata + 1):vt2End) #TODO: might need a +1 on the start of range here

        # Compare the other vtable to the one under consideration.
        # If they are equal, store the offset and break:
        if vtableEqual(b.vtable, objectOffset, vt2)
            existingVtable = vt2Offset
            break
        end
    end

    if existingVtable == 0
        # Did not find a vtable, so write this one to the buffer.

        # Write out the current vtable in reverse , because
        # serialization occurs in last-first order:
        for i = length(b.vtable):-1:1
            off::UOffsetT = 0
            if b.vtable[i] != 0
                # Forward reference to field;
                # use 32bit number to assert no overflow:
                off = objectOffset - b.vtable[i]
            end
            prepend!(b, VOffsetT(off))
        end

        # The two metadata fields are written last.

        # First, store the object bytesize:
        objectSize = objectOffset - b.objectend
        prepend!(b, VOffsetT(objectSize))

        # Second, store the vtable bytesize:
        vbytes = (length(b.vtable) + VtableMetadataFields) * sizeof(VOffsetT)
        prepend!(b, VOffsetT(vbytes))

        # Next, write the offset to the new vtable in the
        # already-allocated SOffsetT at the beginning of this object:
        objectStart = SOffsetT(length(b.bytes) - objectOffset)
        write(b, objectStart, SOffsetT(offset(b) - objectOffset))

        # Finally, store this vtable in memory for future
        # deduplication:
        push!(b.vtables, offset(b))
    else
        # Found a duplicate vtable.

        objectStart = SOffsetT(length(b.bytes) - objectOffset)
        b.head = objectStart

        # Write the offset to the found vtable in the
        # already-allocated SOffsetT at the beginning of this object:
        write(b, b.head, SOffsetT(existingVtable) - SOffsetT(objectOffset))
    end

    empty!(b.vtable)
    return objectOffset
end

"""
`endobject` writes data necessary to finish object construction.
"""
function endobject!(b::Builder)
    assertnested(b)
    n = writevtable!(b)
    b.nested = false
    return n
end

offset(b::Builder) = UOffsetT(length(b.bytes) - b.head)

pad!(b::Builder, n) = foreach(x->place!(b, 0x00), 1:n)

"""
`prep!` prepares to write an element of `size` after `additionalbytes`
have been written, e.g. if you write a string, you need to align such
the int length field is aligned to sizeof(Int32), and the string data follows it
directly.
If all you need to do is align, `additionalbytes` will be 0.
"""
function prep!(b::Builder, size, additionalbytes)
    # Track the biggest thing we've ever aligned to.
    if size > b.minalign
        b.minalign = size
    end
    # Find the amount of alignment needed such that `size` is properly
    # aligned after `additionalBytes`:
    alignsize = xor(Int(-1), (length(b.bytes) - b.head) + additionalbytes) + 1
    alignsize &= (size - 1)

    # Reallocate the buffer if needed:
    totalsize = alignsize + size + additionalbytes
    if b.head <= totalsize
        len = length(b.bytes)
        prepend!(b.bytes, zeros(UInt8, totalsize))
        b.head += length(b.bytes) - len
    end
    pad!(b, alignsize)
    return
end

function Base.prepend!(b::Builder, x::T) where {T}
    prep!(b, sizeof(T), 0)
    place!(b, x)
    return
end

function prependoffset!(b::Builder, off)
    prep!(b, sizeof(Int32), 0) # Ensure alignment is already done.
    if !(off <= offset(b))
        throw(ArgumentError("unreachable: $off <= $(offset(b))"))
    end
    place!(b, SOffsetT(offset(b) - off + sizeof(SOffsetT)))
    return
end

function prependoffsetslot!(b::Builder, o::Int, x::T, d) where {T}
    if x != T(d)
        prependoffset!(b, x)
        slot!(b, o)
    end
    return
end

"""
`startvector` initializes bookkeeping for writing a new vector.

A vector has the following format:
<UOffsetT: number of elements in this vector>
<T: data>+, where T is the type of elements of this vector.
"""
function startvector!(b::Builder, elemSize, numElems, alignment)
    assertnotnested(b)
    b.nested = true
    prep!(b, sizeof(UInt32), elemSize * numElems)
    prep!(b, alignment, elemSize * numElems)
    return offset(b)
end

"""
`endvector` writes data necessary to finish vector construction.
"""
function endvector!(b::Builder, vectorNumElems)
    assertnested(b)
    place!(b, UOffsetT(vectorNumElems))
    b.nested = false
    return offset(b)
end

function createsharedstring!(b::Builder, s::AbstractString)
    get!(b.sharedstrings, s) do
        createstring!(b, s)
    end
end

"""
`createstring!` writes a null-terminated string as a vector.
"""
function createstring!(b::Builder, s::Union{AbstractString, AbstractVector{UInt8}})
    assertnotnested(b)
    b.nested = true
    s = codeunits(s)
    prep!(b, sizeof(UInt32), sizeof(s) + 1)
    place!(b, UInt8(0))

    l = sizeof(s)

    b.head -= l
    copyto!(b.bytes, b.head+1, s, 1, l)
    return endvector!(b, sizeof(s))
end

createbytevector(b::Builder, v) = createstring!(b, v)

function assertnested(b::Builder)
    # If you get this assert, you're in an object while trying to write
    # data that belongs outside of an object.
    # To fix this, write non-inline data (like vectors) before creating
    # objects.
    if !b.nested
        throw(ArgumentError("Incorrect creation order: must be inside object."))
    end
    return
end

function assertnotnested(b::Builder)
    # If you hit this, you're trying to construct a Table/Vector/String
    # during the construction of its parent table (between the MyTableBuilder
    # and builder.Finish()).
    # Move the creation of these view-objects to above the MyTableBuilder to
    # not get this assert.
    # Ignoring this assert may appear to work in simple cases, but the reason
    # it is here is that storing objects in-line may cause vtable offsets
    # to not fit anymore. It also leads to vtable duplication.
    if b.nested
        throw(ArgumentError("Incorrect creation order: object must not be nested."))
    end
    return
end

function assertfinished(b::Builder)
    # If you get this assert, you're attempting to get access a buffer
    # which hasn't been finished yet. Be sure to call builder.Finish()
    # with your root table.
    # If you really need to access an unfinished buffer, use the bytes
    # buffer directly.
    if !b.finished
        throw(ArgumentError("Incorrect use of FinishedBytes(): must call 'Finish' first."))
    end
end

"""
`prependslot!` prepends a `T` onto the object at vtable slot `o`.
If value `x` equals default `d`, then the slot will be set to zero and no
other data will be written.
"""
function prependslot!(b::Builder, o::Int, x::T, d, sh=false) where {T <: Scalar}
    if x != T(d)
        prepend!(b, x)
        slot!(b, o)
    end
    return
end

"""
`prependstructslot!` prepends a struct onto the object at vtable slot `o`.
Structs are stored inline, so nothing additional is being added.
In generated code, `d` is always 0.
"""
function prependstructslot!(b::Builder, voffset, x, d)
    if x != d
        assertnested(b)
        if x != offset(b)
            throw(ArgumentError("inline data write outside of object"))
        end
        slot!(b, voffset)
    end
    return
end

"""
`slot!` sets the vtable key `voffset` to the current location in the buffer.
"""
function slot!(b::Builder, slotnum)
    b.vtable[slotnum + 1] = offset(b)
end

# FinishWithFileIdentifier finalizes a buffer, pointing to the given `rootTable`.
# as well as applys a file identifier
function finishwithfileidentifier(b::Builder, rootTable, fid)
    if length(fid) != fileIdentifierLength
        error("incorrect file identifier length")
    end
    # In order to add a file identifier to the flatbuffer message, we need
    # to prepare an alignment and file identifier length
    prep!(b, b.minalign, sizeof(Int32) + fileIdentifierLength)
    for i = fileIdentifierLength:-1:1
        # place the file identifier
        place!(b, fid[i])
    end
    # finish
    finish!(b, rootTable)
end

"""
`finish!` finalizes a buffer, pointing to the given `rootTable`.
"""
function finish!(b::Builder, rootTable)
    assertnotnested(b)
    prep!(b, b.minalign, sizeof(UOffsetT))
    prependoffset!(b, UOffsetT(rootTable))
    b.finished = true
    return
end

"vtableEqual compares an unwritten vtable to a written vtable."
function vtableEqual(a::Vector{UOffsetT}, objectStart, b::AbstractVector{UInt8})
    if length(a) * sizeof(VOffsetT) != length(b)
        return false
    end

    for i = 0:(length(a)-1)
        x = read(IOBuffer(view(b, (i * sizeof(VOffsetT) + 1):length(b))), VOffsetT)

        # Skip vtable entries that indicate a default value.
        x == 0 && a[i+1] == 0 && continue

        y = objectStart - a[i+1]
        x != y && return false
    end
    return true
end

"""
`place!` prepends a `T` to the Builder, without checking for space.
"""
function place!(b::Builder, x::T) where {T}
    b.head -= sizeof(T)
    write(b, b.head, x)
    return
end

max_forks_repo_name: XpressAI/arrow


text:
using Plots
Random.seed!(0)

include_implicit_dynamics()
include_ddp()
include(joinpath(pwd(), "examples/implicit_dynamics/models/quadruped/model.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/models/quadruped/gait.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/models/quadruped/visuals.jl"))

# Horizon
T = 51
T2 = 2 * T - 1
Tm = 26

# Time step
tf = 1.0
tf2 = 2 * tf
h = tf2 / (T - 1)

# Generate gait
l_torso = 0.367
l_thigh = 0.2
l_leg = 0.2
w_torso = 0.267

θ1 = 0.2 * π
θ2 = θ1 + 0.25 * π
zh = l_thigh * cos(θ1) + l_leg * cos(θ2 - θ1)
θ3 = 0.275 * π
θ4 = acos((zh - l_thigh * cos(θ3)) / l_leg) + θ3

q_body_ref, pf1_ref, pf2_ref, pf3_ref, pf4_ref, contact_modes = generate_gait(θ1, θ2, θ3, θ4, T, Tm, h,
    l_torso = l_torso, l_thigh = l_thigh, l_leg = l_leg, w_torso = w_torso)
# v_body_ref = zero((q_body_ref[end] - q_body_ref[1]) / tf2)
v_body_ref = (q_body_ref[end] - q_body_ref[1]) / tf2

# Initial conditions, controls, disturbances
x1 = [q_body_ref[1]; v_body_ref[1:3]; zeros(3)]
xT = [q_body_ref[end]; v_body_ref[1:3]; zeros(3)]
u_gravity = [0.0; 0.0; model.mass * 9.81 / 4.0]
u_stand = [vcat([u_gravity for i = 1:4]...) for t = 1:T2-1]
# w_stand = [[pf1_ref[1]; pf2_ref[1]; pf3_ref[1]; pf4_ref[1]; 1; 1; 1; 1] for t = 1:T2-1]
w_stand = [[pf1_ref[t]; pf2_ref[t]; pf3_ref[t]; pf4_ref[t]; contact_modes[t]] for t = 1:T2-1]

# Rollout
x̄ = rollout(model, x1, u_stand, w_stand, h, T2)
visualize!(vis, model, x̄,
    [pf1_ref[t] for t = 1:T2],
    [pf2_ref[t] for t = 1:T2],
    [pf3_ref[t] for t = 1:T2],
    [pf4_ref[t] for t = 1:T2],
    Δt = h)

x_ref_2 = [[q_body_ref[t]; v_body_ref[1:3]; zeros(3)] for t = 1:T2]

visualize!(vis, model, x_ref_2,
    [pf1_ref[t] for t = 1:T2],
    [pf2_ref[t] for t = 1:T2],
    [pf3_ref[t] for t = 1:T2],
    [pf4_ref[t] for t = 1:T2],
    Δt = h)

# visualize!(vis, model, x̄, Δt = h)
# x̄ = linear_interpolation(x1, xT, T)
# plot(hcat(x̄...)')

# Objective
Q = [(t < T2 ? Diagonal([10.0; 10.0; 10.0; 1.0e-1 * ones(3); 1.0 * ones(3); 1.0e-1 * ones(3)])
        : Diagonal([10.0; 10.0; 10.0; 1.0e-1 * ones(3); 1.0 * ones(3); 1.0e-1 * ones(3)])) for t = 1:T2]
q = [-2.0 * Q[t] * x_ref_2[t] for t = 1:T2]

R = [0.1 * Diagonal([1.0; 1.0; 1.0e-2; 1.0; 1.0; 1.0e-2; 1.0; 1.0; 1.0e-2; 1.0; 1.0; 1.0e-2]) for t = 1:T2-1]
r = [-2.0 * R[t] * zeros(model.m)  for t = 1:T2-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T2 ? R[t] : nothing, t < T2 ? r[t] : nothing) for t = 1:T2], T2)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

prob = problem_data(model, obj, copy(x̄),
    copy(u_stand),
    # [0.001 * randn(model.m) for t = 1:T2-1],
    copy(w_stand),
    h,
    T2,
	analytical_dynamics_derivatives = true)

# Solve
@time ddp_solve!(prob,
    verbose = true,
    linesearch = :armijo,
    grad_tol = 1.0e-3,
    max_iter = 1000)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

vis = Visualizer()
render(vis)
visualize!(vis, model, x̄,
    [pf1_ref[t] for t = 1:T2],
    [pf2_ref[t] for t = 1:T2],
    [pf3_ref[t] for t = 1:T2],
    [pf4_ref[t] for t = 1:T2],
    q_body_ref = x_ref_2,
    Δt = h)

default_background!(vis)

cone_vis_scale = 0.35

for i = 1:4
    for t in [1, Tm, T, T + Tm, T2]
        foot_vis!(vis, i, t, contact_modes, tl = 0.5, r = 0.05)
    end
    for t in [1, Tm, T, T + Tm, T2-1]
        if t == T2-1
            uit = u[T2-1][(i - 1) * 3 .+ (1:3)]

            cone_vis!(vis, i, t, contact_modes, n = 25, tl = 0.1,
                h = cone_vis_scale * (uit[3] / f_max))#f_max * μ_friction * cone_vis_scale)
            force_vis!(vis, u, i, T2-1,
                f_max = f_max, μ_friction = μ_friction, cone_vis_scale = cone_vis_scale)
        end
    end
end
offset = 0.025
settransform!(vis[:world], Translation(x̄[end-1][1], 0.0, x̄[end-1][3] + offset))
default_background!(vis)

plot(tr2, hcat(ū..., ū[end])', label = "")
tr2 = range(0, stop = tf2, length = T2)
p1 = plot(tr2, hcat([[soc_projection(u[1:3], f_min, f_max, μ_friction) for u in ū]..., soc_projection(ū[end][1:3], f_min, f_max, μ_friction)]...)', linetype = :steppost, labels = "")
p2 = plot(tr2, hcat([[soc_projection(u[3 .+ (1:3)], f_min, f_max, μ_friction) for u in ū]..., soc_projection(ū[end][3 .+ (1:3)], f_min, f_max, μ_friction)]...)', linetype = :steppost, labels = "")
p3 = plot(tr2, hcat([[soc_projection(u[6 .+ (1:3)], f_min, f_max, μ_friction) for u in ū]..., soc_projection(ū[end][6 .+ (1:3)], f_min, f_max, μ_friction)]...)', linetype = :steppost, labels = "")
p4 = plot(tr2, hcat([[soc_projection(u[9 .+ (1:3)], f_min, f_max, μ_friction) for u in ū]..., soc_projection(ū[end][9 .+ (1:3)], f_min, f_max, μ_friction)]...)', linetype = :steppost, labels = "")

plot(p1, p2, p3, p4, layout = (4, 1))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Chain
using Test

@testset "1" begin
    x = [1, 2, 3]
    # one symbol
    y = @chain x begin
        sum
    end
    @test y == sum(x)

    # two expressions
    z = @chain x begin
        *(3)
        sum
    end
    @test z == sum(x .* 3)

    # interleaved expressions
    called = false
    zz = @chain x begin
        .*(3)
        @aside @assert sum(_) / length(_) == 6 # this doesn't change anything
        @aside called = true # this also doesn't do the _ insertion and doesn't change anything
        sum
    end
    @test zz == z
    @test called

    zzz = @chain x begin
        _ .* 3
        sum
    end
    @test zzz == z
end

@testset "2" begin
    x = 1:4
    y = @chain x begin
        filter(isodd, _)
        map(-, _)
        sum
        _ ^ 2
    end
    @test y == 16
end

@testset "nested begin" begin
    x = 1:5
    y = @chain x begin
        begin
            z = sum(_) + 3
            z - 7
        end
        sqrt
    end
    @test y == sqrt(sum(x) + 3 - 7)
end

@testset "no begin" begin
    x = [1, 2, 3]
    y = @chain x sum
    @test y == 6

    f() = 1
    y = @chain f() first
    @test y == 1

    y = @chain x sum(_) max(0, _) first
    @test y == 6

    y = @chain 1 (t -> t + 1)()
    @test y == 2

    y = @chain 1 (t -> t + 1)() first max(0, _)
    @test y == 2

    y = @chain 1 (==(2))
    @test y == false

    y = @chain 1 (==(2)) first (==(false))
    @test y == true

    y = @chain 1 (_ + 1)
    @test y == 2

    y = @chain 1 (_ + 1) first max(0, _)
    @test y == 2

    # the begin block will be different from the normal chain block here
    # only the last statement matters
    y = @chain x begin
        _ .+ 1
        _ .+ 2
    end sum
    @test y == sum(x .+ 2)
end

@testset "invalid invocations" begin
    # just one argument
    @test_throws LoadError eval(quote
        @chain [1, 2, 3]
    end)

    # let block
    @test_throws LoadError eval(quote
        @chain [1, 2, 3] let
            sum
        end
    end)

    # variable defined in chain block doesn't leak out
    z = @chain [1, 2, 3] begin
        @aside inside_var = 5
        @aside @test inside_var == 5
        sum(_) + inside_var
    end
    @test z == 11
    @test_throws UndefVarError inside_var
end

@testset "nested chains" begin
    x = 1:5
    local z
    y = @chain x begin
        _ * 2
        @aside @chain _ begin
            sum(_)
            _ * 2
            @aside z = _
        end
        sum
    end
    @test y == sum(x * 2)
    @test z != x * 2
end

@testset "broadcast macro symbol" begin
    x = 1:5
    y = @chain x begin
        @. sin
        sum
    end
    @test y == sum(sin.(x))

    ## leave non-symbol invocations intact
    yy = @chain x begin
        @. sin(_)
        sum
    end
    @test yy == sum(sin.(x))
end

macro sin(exp)
    :(sin($(esc(exp))))
end

macro broadcastminus(exp1, exp2)
    :(broadcast(-, $(esc(exp1)), $(esc(exp2))))
end

@testset "splicing into macro calls" begin

    x = 1
    y = @chain x begin
        @sin
    end
    @test y == sin(x)

    y = @chain x begin
        @sin()
    end
    @test y == sin(x)

    xx = [1, 2, 3, 4]
    yy = @chain xx begin
        @broadcastminus(2.5)
    end
    @test yy == broadcast(-, xx, 2.5)

    xxx = [1, 2, 3, 4]
    yyy = @chain xxx begin
        @broadcastminus(2.5, _)
    end
    @test yyy == broadcast(-, 2.5, xxx)
end

@testset "single arg version" begin
    x = [1, 2, 3]

    xx = @chain begin
        x
    end
    @test xx == x

    # this has a different internal structure (one LineNumberNode missing I think)
    @test x == @chain begin
        x
    end

    @test sum(x) == @chain begin
        x
        sum
    end

    y = @chain begin
        x
        sum
    end
    @test y == sum(x)

    z = @chain begin
        x
        @. sqrt
        sum(_)
    end
    @test z == sum(sqrt.(x))

    @test sum == @chain begin
        sum
    end
end

@testset "invalid single arg versions" begin
    # empty
    @test_throws LoadError eval(quote
        @chain begin
        end
    end)

    # rvalue _ errors
    @test_throws ErrorException eval(quote
        @chain begin
            _
        end
    end)

    @test_throws ErrorException eval(quote
        @chain begin
            sum(_)
        end
    end)
end

@testset "handling keyword argments" begin
    f(a; kwarg) = (a, kwarg)
    @test (:a, :kwarg) == @chain begin
        :a
        f(kwarg = :kwarg)
    end
    @test (:a, :kwarg) == @chain begin
        :a
        f(; kwarg = :kwarg)
    end
end

# issue 13
@testset "no argument call" begin
    x = 1
    y = @chain x begin
        sin()
    end
    @test y == sin(x)
end

# issue 13
@testset "broadcasting calls" begin

    xs = [1, 2, 3]
    ys = @chain xs begin
        sin.()
    end
    @test ys == sin.(xs)

    add(x, y) = x + y

    zs = [4, 5, 6]
    sums = @chain xs begin
        add.(zs)
    end
    @test sums == add.(xs, zs)
end

# issue 16
@testset "empty chain" begin
    a = 2
    x = @chain a + 1 begin
    end
    @test x == 3

    y = @chain begin
        a + 1
    end
    @test y == 3
end

module LocalModule
    function square(xs)
        xs .^ 2
    end

    function power(xs, pow)
        xs .^ pow
    end

    add_one(x) = x + 1

    macro sin(exp)
        :(sin($(esc(exp))))
    end

    macro broadcastminus(exp1, exp2)
        :(broadcast(-, $(esc(exp1)), $(esc(exp2))))
    end

    module SubModule
        function square(xs)
            xs .^ 2
        end

        function power(xs, pow)
            xs .^ pow
        end

        add_one(x) = x + 1

        macro sin(exp)
            :(sin($(esc(exp))))
        end

        macro broadcastminus(exp1, exp2)
            :(broadcast(-, $(esc(exp1)), $(esc(exp2))))
        end
    end
end

@testset "Module qualification" begin

    using .LocalModule

    xs = [1, 2, 3]
    pow = 4
    y = @chain xs begin
        LocalModule.square
        LocalModule.power(pow)
        Base.sum
    end
    @test y == sum(LocalModule.power(LocalModule.square(xs), pow))

    y2 = @chain xs begin
        LocalModule.SubModule.square
        LocalModule.SubModule.power(pow)
        Base.sum
    end
    @test y == sum(LocalModule.SubModule.power(LocalModule.SubModule.square(xs), pow))

    y3 = @chain xs begin
        @. LocalModule.add_one
        @. LocalModule.SubModule.add_one
    end
    @test y3 == LocalModule.SubModule.add_one.(LocalModule.add_one.(xs))

    y4 = @chain xs begin
        LocalModule.@broadcastminus(2.5)
    end
    @test y4 == LocalModule.@broadcastminus(xs, 2.5)

    y5 = @chain xs begin
        LocalModule.SubModule.@broadcastminus(2.5)
    end
    @test y5 == LocalModule.SubModule.@broadcastminus(xs, 2.5)

    y6 = @chain 3 begin
        LocalModule.@sin
    end
    @test y6 == LocalModule.@sin(3)

    y7 = @chain 3 begin
        LocalModule.SubModule.@sin
    end
    @test y7 == LocalModule.SubModule.@sin(3)
end

function kwfunc(y; x = 1)
    y * x
end

macro kwmac(exprs...)
    :(kwfunc($(esc.(exprs)...)))
end

@testset "keyword arguments" begin

    @test 6 == @chain 2 begin
        kwfunc(; x = 3)
    end

    @test 6 == @chain 2 begin
        @kwmac(; x = 3)
    end
end

@testset "@aside at the end" begin
    x = 1

    @test 1 == @chain x begin
        @aside 1 + 2
    end

    @test 2 == @chain x begin
        _ + 1
        @aside 1 + 2
    end
end

@testset "workaround for docstring parsing" begin
    @test "hi" == @chain " hi " strip
    @test "hi" == @chain " hi " begin
        strip
    end
    @test "hi" == @chain begin
        " hi "
        strip
    end
    @test "hi" == @chain begin
        "hi"
        " $_ "
        strip
    end
    @test "A" == @chain begin
        'a'
        " $_ "
        strip
        "$_"
        uppercase
    end
end

@testset "nested single line chain" begin
    @test 36 == @chain 1:3 begin
        @chain _ sum _ ^ 2
    end
end


text:
using Documenter, MLJLinearModels

makedocs(
    modules = [MLJLinearModels],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        ),
    sitename = "MLJLinearModels.jl",
    authors = "Thibaut Lienart, and contributors.",
    pages = [
        "Home" => "index.md",
        "MLJ"  => "mlj.md"
    ]
)

deploydocs(
    repo = "github.com/alan-turing-institute/MLJLinearModels.jl"
)


text:
function read2(f, block::PlainVariableBlockHeader{T}; fn=_->true, req_pts, batch_size=(req_pts[end][2]-req_pts[end][1])) where T
    if length(block.npts) != length(req_pts)
        ArgumentError("dimension mismatch")
    end

    block_size = prod(block.npts[1:end-1])
    z1, z2 = req_pts[end]
    seek(f, block.base_header.data_location + sizeof(T)*block_size*(z1-1))
    n, m = divrem(z2-z1, batch_size)
    vals = T[]
    idxs = Array{CartesianIndex{length(req_pts)},1}()
    result = Array{T, 1}(undef, block_size*batch_size)
    for i in eachindex(n)
        read_batch!(f, result, block, fn, req_pts, batch_size, block_size, vals, idxs)
    end
    if m != 0
        result = Array{T, 1}(undef, block_size*m)
        read_batch!(f, result, block, fn, req_pts, m, block_size, vals, idxs)
    end

    for (i,v) in enumerate(req_pts)
        idxs .+= Ref(cart_idx(req_pts[i][1], Val(i)))
    end

    vals, idxs
end

function read_batch!(f, result, block, fn, req_pts, batch_size, block_size, vals, idxs)
    read!(f, result)
    result = trim_block(result, block.npts, batch_size, req_pts, Val(length(req_pts)))
    idx = findall(fn, result)
    append!(vals, result[idx])
    append!(idxs, idx)
end

function trim_block(result, npts, batch_size, req_pts, ::Val{1})
    result
end

function trim_block(result, npts, batch_size, req_pts, ::Val{2})
    x1, x2 = req_pts[1]
    result = reshape(result, (npts[1], batch_size))
    result[x1:x2, :]
end

function trim_block(result, npts, batch_size, req_pts, ::Val{3})
    x1, x2 = req_pts[1]
    y1, y2 = req_pts[2]
    result = reshape(result, (npts[1], npts[2], batch_size))
    result[x1:x2, y1:y2, :]
end

cart_idx(pt, ::Val{1}) = CartesianIndex(pt-1, 0, 0)
cart_idx(pt, ::Val{2}) = CartesianIndex(0, pt-1, 0)
cart_idx(pt, ::Val{3}) = CartesianIndex(0, 0, pt-1)


text:
using NeXLSpectrum

path=raw"C:\Users\nritchie\OneDrive - National Institute of Standards and Technology (NIST)\Documents\FY2016\SE1610x\Manual\Particle 12 keV"
p = joinpath(path,"map.ptx")

props = Dict(:BeamEnergy => 12.0e3, :ProbeCurrent=>1.0, :RealTime => 0.129202)

ht1 = @time readptx(p, LinearEnergyScale(0.0,10.0), props, 1024, blocksize=8, dets=(true, false, false, false))
ht2 = @time readptx(p, LinearEnergyScale(0.0,10.0), props, 1024, blocksize=8, dets=(false, true, false, false))
ht3 = @time readptx(p, LinearEnergyScale(0.0,10.0), props, 1024, blocksize=8, dets=(false, false, true, false))
ht4 = @time readptx(p, LinearEnergyScale(0.0,10.0), props, 1024, blocksize=8, dets=(false, false, false, true))

specs = [ ht1, ht2, ht3, ht4 ]

NeXLSpectrum.checkmultispec(specs)

map(h->h[:RealTime], specs)

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
module test_julia_namedtuples

using Test

@test NamedTuple() == (;)

@test (a=1, b=2) == (a=1, :b=>2)
@test (b=2,) == (; b=2) == (; :b=>2)

b = 2
@test (; b) == (b=2,)

@test NamedTuple{(:a, :b),Tuple{Int,Int}}((1, 2)) == NamedTuple{(:a, :b)}((1, 2)) == (a=1, b=2)

d = Dict(:a=>1, :b=>2)
@test (; d...) == (a=1, b=2)
@test d == Dict(pairs((a=1, b=2)))

nt = (a = 1, b = 2)
@test NamedTuple{keys(nt)}(values(nt)) == nt

@test NamedTuple{(:a,:b), Tuple{Int, NamedTuple{(:c,)}}}((a=1, b=(c=1,))) == (a=1, b=(c=1,))

f(; kwargs...) = kwargs
@test f(; a=2, b=3) == f(; :a=>2, :b=>3)
@test  (; a=2, b=3) ==  (; :a=>2, :b=>3)
@test (; f(; a=2, b=3)...) == (; a=2, b=3) == (a=2, b=3)

name = "a"
@test NamedTuple{(Symbol(name),)}(2) == (a=2,)

@test NamedTuple{(), Tuple{}}((;)) == (;)

@test (;:).:(:)     === (:)
@test (;*).:*       === *
@test (;<:).:<:     === <:
@test (;<:, >:).:>: === >:
@test (;===).:(===) === ===
@test (;~,!,%,^,&,*,-,+,==,|,\,:,<,>,/) == (~ = ~, ! = !, % = rem, ^ = ^, (&) = &, * = *, - = -, + = +, == = ==, | = |, \ = \, (:) = Colon(), < = <, > = >, / = /)

end # module test_julia_namedtuples

module test_julia_namedtuples_structdiff

using Test

@test Base.structdiff((a=1, b=2), NamedTuple{(:a,)}) == (b=2,)
@test Base.structdiff((a=1, b=2), (a=1,))            == (b=2,)

end # module test_julia_namedtuples_structdiff


text:
using LazyArtifacts
using NPZ
using Tables
import TOML
import JSON
using Pkg
using Pkg.Artifacts

const _metadata = TOML.parsefile(joinpath(artifact"metadata", "data.toml"))
const datasets = keys(_metadata)

_check_version(version, min_version=v"1.1.0") = VersionNumber(join(version, ".")) >= min_version

struct TabularDataset{
        DataElType <: Number,
        TargetElType <: Number,
        ExtrasType <: Union{Dict{String, Any},Nothing}
    }

    input::Matrix{DataElType}
    output::Matrix{TargetElType}
    metadata::Dict{String, Any}
    extras::ExtrasType
    input_names::Vector{Symbol}
    output_names::Vector{Symbol}

    function TabularDataset(input, output, metadata, extras=nothing)
        ndims(input) != 2 && throw("input must have exactly 2 dimensions")
        ndims(output) ∉ 2 && throw("output must have 2 dimensions")
        size(input, 2) != size(output, 2) &&
            throw("input and output must have the same number of examples")

        in_type = eltype(input)
        out_type = eltype(output)
        extras_type = extras === nothing ? Nothing : Dict{String, Any}

        # check input column names
        if extras !== nothing && "column-names-attributes" in keys(extras)
            input_names = [Symbol(name) for name in extras["column-names-attributes"]]
        else
            input_names = [Symbol("A$(i)") for i in 1:size(input, 1)]
        end

        length(input_names) != size(input, 1) &&
            throw("number of attribute names does not match the number of attributes")

        # check output/target column names
        if extras !== nothing && "column-names-target" in keys(extras)
            output_names = [Symbol(name) for name in extras["column-names-target"]]
        elseif size(output, 1) > 1
            output_names = [Symbol("label$(i)") for i in 1:size(output, 1)]
        else
            output_names = [:label]
        end

        length(output_names) != size(output, 1) &&
            throw("number of output labels does not match the number of outputs")

        new{in_type, out_type, extras_type}(input, output, metadata, extras, input_names, output_names)
    end
end

function TabularDataset(name::AbstractString, split=:train)
    name = lowercase(name)
    name ∉ keys(_metadata) && throw("Unknown dataset `$name`")

    _check_version(npzread(
        joinpath(@artifact_str("$name-npz"), "$name.npz"),
        ["_version"]
    )["_version"])

    data_dict = npzread(
        joinpath(@artifact_str("$name-npz"), "$name.npz"),
        ["$split-data", "$split-labels"]
    )

    input = transpose(data_dict["$split-data"])
    output = transpose(data_dict["$split-labels"])
    metadata = _metadata[name]
    if haskey(metadata, "extras_location")
        extras = JSON.parsefile(joinpath(@artifact_str("$name-extras"), "$name.json"), dicttype=Dict, inttype=Int64)
    else
        extras = nothing
    end

    return TabularDataset(input, output, metadata, extras)
end

# "basic" getters and info
num_examples(ds::TabularDataset) = size(ds.input, 2)
num_inputs(ds::TabularDataset) = size(ds.input, 1)
num_outputs(ds::TabularDataset) = size(ds.output, 1)

task(ds::TabularDataset) = ds.metadata["task"]
num_classes(ds::TabularDataset) =
    if task(ds) != "classification"
        throw("non-classification datasets don't have a number of classes")
    else
        return ds.metadata["classes"]
    end

has_extras(ds::TabularDataset) = ds.extras !== nothing
has_extra(ds::TabularDataset, extra_name) = has_extras(ds) && extra_name in keys(ds.extras)

categorical_attributes(ds::TabularDataset) =
    has_extra(ds, "categories") ? keys(ds.extras["categories"]) : nothing

license(ds::TabularDataset) = has_extra(ds, "license") ? ds.extras["license"] : nothing
bibtex(ds::TabularDataset) = has_extra(ds, "bibtex") ? ds.extras["bibtex"] : nothing

# indexing and iteration protocol on the dataset itself
Base.getindex(ds::TabularDataset) = (ds.input, ds.output)
Base.getindex(ds::TabularDataset, i) = (ds.input[:, i], ds.output[:, i])
Base.firstindex(ds::TabularDataset) = 1
Base.lastindex(ds::TabularDataset) = length(ds)
Base.length(ds::TabularDataset) = num_examples(ds)

IndexStyle(::Type{TabularDataset}) = IndexLinear()

Base.eltype(::Type{TabularDataset{InType, OutType, ExtrasType}}) where {InType, OutType, ExtrasType} = Tuple{Vector{InType}, Vector{OutType}}

Base.iterate(ds::TabularDataset, i=1) = i <= num_examples(ds) ? (ds[i], i + 1) : nothing

# Tables.jl interface
Tables.istable(::Type{TabularDataset}) = true
Tables.rowaccess(::Type{TabularDataset}) = true
Tables.columnaccess(::Type{TabularDataset}) = true
Tables.columnnames(ds::TabularDataset) = vcat(ds.input_names, ds.output_names)

Tables.schema(ds::TabularDataset) = Tables.Schema{nothing, nothing}(
    vcat(ds.input_names, ds.output_names),
    vcat(repeat([eltype(ds.input)], num_inputs(ds)), repeat([eltype(ds.output)], num_outputs(ds)))
)

# rows interface
struct TabularDatasetRow{DatasetType <: TabularDataset}
    ds::DatasetType
    name_indices::Dict{Symbol, Int}
    index::Int
end

Tables.columnnames(row::TabularDatasetRow) = vcat(row.ds.input_names, row.ds.output_names)
function Tables.getcolumn(row::TabularDatasetRow, i::Int)
    if i > num_inputs(row.ds)
        return row.ds.output[i - num_inputs(row.ds), ds.index]
    else
        return row.ds.input[i, ds.index]
    end
end
Tables.getcolumn(row::TabularDatasetRow, nm::Symbol) = (println(row.name_indices[nm]); Tables.getcolumn(row, row.name_indices[nm]))
Tables.getcolumn(row::TabularDatasetRow, ::Type, i::Int, ::Symbol) = Tables.getcolumn(row, i)

struct TabularDatasetRows{DatasetType <: TabularDataset}
    ds::DatasetType
    name_indices::Dict{Symbol, Int}
end

function Base.iterate(iter::TabularDatasetRows, row_num=1)
    if row_num > num_examples(iter.ds)
        return nothing
    else
        return TabularDatasetRow{typeof(iter.ds)}(iter.ds, iter.name_indices, row_num), row_num + 1
    end
end

function Tables.rows(ds::TabularDataset)
    all_names = vcat(ds.input_names, ds.output_names)
    return TabularDatasetRows{typeof(ds)}(ds, Dict(name => index for (index, name) in enumerate(all_names)))
end

# column interface
Tables.columns(ds::TabularDataset) = ds
function Tables.getcolumn(ds::TabularDataset, i::Int)
    if i > num_inputs(ds)
        return ds.output[i - num_inputs(ds), :]
    else
        return ds.input[i, :]
    end
end
function Tables.getcolumn(ds::TabularDataset, nm::Symbol)
    index = findfirst(==(nm), ds.input_names)
    if index !== nothing
        return ds.input[index, :]
    else
        return ds.output[findfirst(==(nm), ds.output_names), :]
    end
end
Tables.getcolumn(ds::TabularDataset, ::Type, i::Int, ::Symbol) = Tables.getcolumn(ds, i)

# artifact pre-downloading

function ensure_downloaded(artifact_names...)
    artifacts_file = abspath(joinpath(@__DIR__, "..", "Artifacts.toml"))
    if length(artifact_names) == 0
        Pkg.Artifacts.ensure_all_artifacts_installed(artifacts_file)
    else
        for artifact_name in artifact_names
            Pkg.Artifacts.ensure_artifact_installed("$artifact_name-npz", artifacts_file)
            Pkg.Artifacts.ensure_artifact_installed("$artifact_name-extras", artifacts_file)
        end
    end
end


text:
using DynamicProgramming
using Documenter

makedocs(;
    modules=[DynamicProgramming],
    authors="Michele Fornino <fornino.michele@gmail.com>",
    repo="https://github.com/mfornino/DynamicProgramming.jl/blob/{commit}{path}#L{line}",
    sitename="DynamicProgramming.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mfornino.github.io/DynamicProgramming.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/mfornino/DynamicProgramming.jl",
)


text:
include(joinpath("Thesis model", "glu_ed_selfrep_model.jl"))
import .Glu_ED_Selfrep_Model

res_ed_selfrep = []
for var in [0.0001:0.0005:0.005;]
    push!(res_ed_selfrep, Glu_ED_Selfrep_Model.glu_ed_selfrep_model(var))
end

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

using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Thesis model", "Figures", "ED_Model")
resultspath = joinpath("Thesis model", "Results", "ED_Model")

#CSV.write(joinpath(resultspath, "Results_ED_Selfrep.csv"), res_ed[1])

x = [results_ed_selfrep["glc_ext"] for results_ed_selfrep in res_ed_selfrep]
y = [results_ed_selfrep["mu"] for results_ed_selfrep in res_ed_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter_glc = scatter!(ax, x, y)
line_glc = lines!(ax, x, y)

ax.xlabel = "External glucose concetration [mM]"
ax.ylabel = "Biomass [1/h]"

fig

FileIO.save(joinpath(imgpath, "Glu_ED_Selfrep_Growth.pdf"), fig)

## ##

x = [results_ed_selfrep["mu"] for results_ed_selfrep in res_ed_selfrep]
a = [results_ed_selfrep["g6p"] for results_ed_selfrep in res_ed_selfrep]
b = [results_ed_selfrep["glu"] for results_ed_selfrep in res_ed_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Metabolite concentration"

fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2]],
    ["Glucose-6-Phosphate", "Glutamate"],
    tellheight = false,
    tellwidth = false,
    halign = :left,
    valign = :bottom,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Glu_ED_G6P+GLU_Selfrep.pdf"), fig)

## Enzyme concentration band plot ##

x = [results_ed_selfrep["mu"] for results_ed_selfrep in res_ed_selfrep]
c = [results_ed_selfrep["pts"] for results_ed_selfrep in res_ed_selfrep]
d = [results_ed_selfrep["ed"] for results_ed_selfrep in res_ed_selfrep]
e = [results_ed_selfrep["pyk"] for results_ed_selfrep in res_ed_selfrep]
f = [results_ed_selfrep["ldh"] for results_ed_selfrep in res_ed_selfrep]
g = [results_ed_selfrep["burn"] for results_ed_selfrep in res_ed_selfrep]
h = [results_ed_selfrep["lp"] for results_ed_selfrep in res_ed_selfrep]

fracs = [c, d, e, f, g, h]
fraclabels = [
    "PTS",
    "ED",
    "Pyruvate kinase",
    "Lactate dehydrogenase",
    "ATP burning enzyme",
    "Lactate permease",
]
function plot_proteome(ax, x, fracs, fraclabels)
    d = cumsum(hcat(fracs...)'; dims = 1)
    d = [zeros(size(d, 2))'; d ./ d[end, :]']
    for i = 1:length(fraclabels)
        band!(ax, x, d[i, :], d[i+1, :], label = fraclabels[i])
    end
end

f = Figure();
ax = Axis(f[1, 1]);
plot_proteome(ax, 1:10, fracs, fraclabels)
ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Proteome fraction"
f[1, 2] = Legend(f, ax, "Protein", unique = true, framevisible = false)
f

FileIO.save(joinpath(imgpath, "Glu_ED_Selfrep_Enzymeband.pdf"), f)

## Rates ##

x = [results_ed_selfrep["mu"] for results_em_selfrep in res_ed_selfrep]
i = [results_ed_selfrep["v_pts"] for results_em_selfrep in res_ed_selfrep]
j = [results_ed_selfrep["v_pyk"] for results_em_selfrep in res_ed_selfrep]
k = [results_ed_selfrep["v_ldh"] for results_em_selfrep in res_ed_selfrep]
l = [results_ed_selfrep["v_burn"] for results_em_selfrep in res_ed_selfrep]
m = [results_ed_selfrep["v_lp"] for results_em_selfrep in res_ed_selfrep]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, i)
line1 = lines!(ax, x, i)
scatter2 = scatter!(ax, x, j)
line2 = lines!(ax, x, j)
scatter3 = scatter!(ax, x, k)
line3 = lines!(ax, x, k)
scatter4 = scatter!(ax, x, l)
line4 = lines!(ax, x, l)
scatter5 = scatter!(ax, x, m)
line5 = lines!(ax, x, m)

ax.xlabel = "Biomass [1/h]"
ax.ylabel = "Fluxes [mmol/gDW/h]"

fig[1, 1] = Legend(
    fig,
    [
        [scatter1, line1],
        [scatter2, line2],
        [scatter3, line3],
        [scatter4, line4],
        [scatter5, line5],
    ],
    ["v_pts", "v_pyk", "v_ldh", "v_burn", "v_lp"],
    tellheight = false,
    tellwidth = false,
    halign = :left,
    valign = :bottom,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Glu_ED_Selfrep_Rates.pdf"), fig)


text:
# ---
# title: JuliaRL\_MAC\_CartPole
# cover: assets/JuliaRL_MAC_CartPole.png
# description: MAC applied to CartPole
# date: 2021-05-22
# author: "[Raj Ghugare](https://github.com/RajGhugare19)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:MAC},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    N_ENV = 16
    UPDATE_FREQ = 20
    env = MultiThreadEnv([
        CartPoleEnv(; T = Float32, rng = StableRNG(hash(seed + i))) for i in 1:N_ENV
    ])
    ns, na = length(state(env[1])), length(action_space(env[1]))
    RLBase.reset!(env, is_force = true)

    agent = Agent(
        policy = QBasedPolicy(
            learner = MACLearner(
                approximator = ActorCritic(
                    actor = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 30, relu; init = glorot_uniform(rng)),
                            Dense(30, 30, relu; init = glorot_uniform(rng)),
                            Dense(30, na; init = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(1e-2),
                    ),
                    critic = NeuralNetworkApproximator(
                        model = Chain(
                            Dense(ns, 30, relu; init = glorot_uniform(rng)),
                            Dense(30, 30, relu; init = glorot_uniform(rng)),
                            Dense(30, na; init = glorot_uniform(rng)),
                        ),
                        optimizer = ADAM(3e-3),
                    ),
                ) |> cpu,
                γ = 0.99f0,
                bootstrap = true,
                update_freq = UPDATE_FREQ,
            ),
            explorer = BatchExplorer(GumbelSoftmaxExplorer()),#= seed = nothing =#
        ),
        trajectory = CircularArraySARTTrajectory(;
            capacity = UPDATE_FREQ,
            state = Matrix{Float32} => (ns, N_ENV),
            action = Vector{Int} => (N_ENV,),
            reward = Vector{Float32} => (N_ENV,),
            terminal = Vector{Bool} => (N_ENV,),
        ),
    )

    stop_condition = StopAfterStep(50_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalBatchRewardPerEpisode(N_ENV)
    Experiment(agent, env, stop_condition, hook, "# MAC with CartPole")
end

#+ tangle=false
using Plots
using Statistics
pyplot() #hide
ex = E`JuliaRL_MAC_CartPole`
run(ex)
n = minimum(map(length, ex.hook.rewards))
m = mean([@view(x[1:n]) for x in ex.hook.rewards])
s = std([@view(x[1:n]) for x in ex.hook.rewards])
plot(m,ribbon=s)
savefig("assets/JuliaRL_MAC_CartPole.png") #hide

# ![](assets/JuliaRL_MAC_CartPole.png)


text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

"""
    RectangularAperture(x, y, w, h, θ)
    RectangularAperture(position, w, h, θ)

A rectangular aperture.

A rectangular aperture with width `w`, height `h`, and position angle `θ` in degrees.

# Examples
```jldoctest
julia> ap = RectangularAperture(0, 0, 10, 4, 0)
RectangularAperture(0, 0, w=10, h=4, θ=0°)
```
"""
struct RectangularAperture{T <: Number} <: AbstractAperture
    x::T
    y::T
    w::T
    h::T
    theta::T
    function RectangularAperture(x::T, y::T, w::T, h::T, θ::T) where T <: Number
        (w < 0  || h < 0) && error("Invalid side lengths w=$w, h=$h. Both must be positive.")
        new{T}(x, y, w, h, mod(θ, 360))
    end
end

RectangularAperture(center::AbstractVector, w, h, θ) = RectangularAperture(center..., w, h, θ)
RectangularAperture(x, y,  w, h, θ) = RectangularAperture(promote(x, y, w, h, θ)...)

function Base.show(io::IO, ap::RectangularAperture)
    print(io, "RectangularAperture($(ap.x), $(ap.y), w=$(ap.w), h=$(ap.h), θ=$(ap.theta)°)")
end

function bbox(ap::RectangularAperture{T}) where T
    w2 = ap.w / 2
    h2 = ap.h / 2
    sint, cost = sincos(deg2rad(ap.theta))

    dx1 = abs(w2 * cost - h2 * sint)
    dy1 = abs(w2 * sint + h2 * cost)
    dx2 = abs(w2 * cost + h2 * sint)
    dy2 = abs(w2 * sint - h2 * cost)

    dx = max(dx1, dx2)
    dy = max(dy1, dy2)

    xmin = ceil(Int, ap.x - dx - 0.5)
    ymin = ceil(Int, ap.y - dy - 0.5)
    xmax = ceil(Int, ap.x + dx - 0.5)
    ymax = ceil(Int, ap.y + dy - 0.5)
    return xmin, xmax, ymin, ymax
end

function mask(ap::RectangularAperture; method = :exact)
    bounds = edges(ap)
    ny, nx = size(ap)
    return rectangular_overlap(bounds..., nx, ny, ap.w, ap.h, ap.theta, method = method)
end

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

"""
    RectangularAnnulus(x, y, w_in, w_out, h_out, θ)
    RectangularAnnulus(position, w_in, w_out, h_out, θ)

A rectangular annulus with inner width `w_in`, outer width `w_out`, outer height `h_out`, and position angle `θ` in degrees. `h_in` is automatically calculated from `w_in / w_out * h_out`. Note that `w_out ≥ w_in > 0`.

# Examples
```jldoctest
julia> ap = RectangularAnnulus(0, 0, 5, 10, 8, 45)
RectangularAnnulus(0.0, 0.0, w_in=5.0, w_out=10.0, h_in=4.0, h_out=8.0, θ=45.0°)
```
"""
struct RectangularAnnulus{T <: Number} <: AbstractAperture
    x::T
    y::T
    w_in::T
    w_out::T
    h_in::T
    h_out::T
    theta::T

    function RectangularAnnulus(x::T, y::T, w_in::T, w_out::T, h_out::T, θ::T) where T <: Number
        0 < w_in ≤ w_out || error("Invalid sides ($w_in, $w_out). `w_out` must be greater than or equal to `w_in` which must be greater than 0.")
        h_out ≤ 0 && error("Invalid side($h_out). `h_out` must be greater than 0.")
        P = typeof(one(T) / 1)
        new{P}(x, y, w_in, w_out, w_in / w_out * h_out, h_out, mod(θ, 360))
    end
end

RectangularAnnulus(center::AbstractVector, w_in, w_out, h_out, θ) = RectangularAnnulus(center..., w_in, w_out, h_out, θ)
RectangularAnnulus(x, y, w_in, w_out, h_out, θ) = RectangularAnnulus(promote(x, y, w_in, w_out, h_out, θ)...)

function Base.show(io::IO, ap::RectangularAnnulus)
    print(io, "RectangularAnnulus($(ap.x), $(ap.y), w_in=$(ap.w_in), w_out=$(ap.w_out), h_in=$(ap.h_in), h_out=$(ap.h_out), θ=$(ap.theta)°)")
end

function bbox(ap::RectangularAnnulus)
    w2 = ap.w_out / 2
    h2 = ap.h_out / 2
    sint, cost = sincos(deg2rad(ap.theta))

    dx1 = abs(w2 * cost - h2 * sint)
    dy1 = abs(w2 * sint + h2 * cost)
    dx2 = abs(w2 * cost + h2 * sint)
    dy2 = abs(w2 * sint - h2 * cost)

    dx = max(dx1, dx2)
    dy = max(dy1, dy2)

    xmin = ceil(Int, ap.x - dx - 0.5)
    ymin = ceil(Int, ap.y - dy - 0.5)
    xmax = ceil(Int, ap.x + dx - 0.5)
    ymax = ceil(Int, ap.y + dy - 0.5)
    return (xmin, xmax, ymin, ymax)
end

function mask(ap::RectangularAnnulus; method = :exact)
    bounds = edges(ap)
    ny, nx = size(ap)
    out = rectangular_overlap(bounds..., nx, ny, ap.w_out, ap.h_out, ap.theta, method = method)
    out .-= rectangular_overlap(bounds..., nx, ny,  ap.w_in, ap.h_in, ap.theta, method = method)
end


text:
# Function to calculate the Coappearance Matrix Q
function COAPPEARENCE(trans::Matrix{Int64})
    t_sparse = dropzeros(sparse(trans))
    Q = t_sparse'*t_sparse
    for i = 1:size(Q,1)
        Q[i,i] = 0
    end
    Q = Matrix(Q)
    Q = Int.(Q)
    return Q
end

function COAPPEARENCE(trans::SparseMatrixCSC{Float64, Int64})
    Q = trans'*trans
    for i = 1:size(Q,1)
        Q[i,i] = 0
    end
    Q = Matrix(Q)
    Q = Int.(Q)
    return Q
end

# Functions to evaluate the number of parcels necessary to fullfil all orders
## COMBINEWAREHOUSES: create all possible warehouse combinations
function COMBINEWAREHOUSES(capacity::Array{Int64,1})
    combination = [[k] for k = 1:size(capacity,1)]
    combinations_X = combinations(combination)
    combination = [c for c in combinations_X]
    return combination::Array{Array{Array{Int64,1},1},1}
end

## PARCELSSEND: number of parcels necessary to fulfill all orders
function PARCELSSEND(trans::SparseMatrixCSC{Float64, Int64},
                     X::Array{Int64,2},
                     capacity::Array{Int64,1},
                     combination::Array{Array{Array{Int64,1},1},1})
    t_sparse = trans
    if sum(capacity) == size(t_sparse,2)
        parcel = t_sparse * X
        for j = 1:size(parcel,1)
            for k = 1:size(parcel,2)
                if parcel[j,k] > 0
                    parcel[j,k] = 1
                else
                    parcel[j,k] = 0
                end
            end
        end
        parcel_out = round(Int64, sum(parcel))
    else
        warehouse_combination = Array{Int64,2}(undef,size(combination,1),size(X,1)) .= 0
        for i = 1:size(combination,1)
            for j = 1:size(combination[i],1)
                warehouse_combination[i,:] += X[:,combination[i][j]]
            end
        end
        parcel = spzeros(size(t_sparse,1),size(capacity,1))
        for i in 1:size(t_sparse,1)
            for j in 1:size(warehouse_combination,1)
                if all(x->x>=0, warehouse_combination[j,:] - t_sparse[i,:])
                    for k = 1:size(combination[j],1)
                        parcel[i,combination[j][k]] .= 1
                    end
                    break
                end
            end
        end
        parcel_out = round(Int64, sum(parcel))
    end
    return parcel_out::Int64
end

# Functions for the random allocation of SKUs to warehouses
## RANDOMALLOCONCE: allocate SKUs randomly in case each SKU can only be assigned once
function RANDOMALLOCONCE(trans::SparseMatrixCSC{Float64, Int64},
                         capacity::Array{Int64,1})
    X = Array{Int64,2}(undef,size(trans,2),size(capacity,1))
    X .= 0
    for j in 1:size(X,1)
        while sum(X[j,:]) == 0
            randomnumber = rand(1:size(capacity,1))
            if sum(X[:,randomnumber]) < capacity[randomnumber]
                X[j,randomnumber] = 1
            end
        end
    end
    X = X[shuffle(1:size(X,1)),:]
    return X::Array{Int64,2}
end

## RANDOMALLOCMULTI: allocate SKUs randomly in case each SKU can be allocated multiple times
function RANDOMALLOCMULTI(trans::SparseMatrixCSC{Float64, Int64},
                          capacity::Array{Int64,1})
    X = RANDOMALLOCONCE(trans,capacity)
    for d = 1:size(capacity,1)
        while sum(X[:,d]) < sum(capacity[d])
            randomproduct = rand(1:size(trans,2))
            if X[randomproduct,d] == 0
                X[randomproduct,d] = 1
            end
        end
    end
    return X::Array{Int64,2}
end

## RANDOMBENCH: benchmark function to evaluate multiple random allocations
function RANDOMBENCH(trans::SparseMatrixCSC{Float64, Int64},
                     capacity::Array{Int64,1},
                     iterations::Int64,
                     combination::Array{Array{Array{Int64,1},1},1})
    randomcollection = Array{Int64,1}(undef,iterations) .= 0
    Threads.@threads for r = 1:iterations
        W = RANDOMALLOCMULTI(trans,capacity)
        parcel = PARCELSSEND(trans,W,capacity,combination)
        Threads.lock(ren_lock) do
            randomcollection[r] = parcel
        end
    end
    average_parcels = round(sum(randomcollection)/iterations,digits=0)
    return average_parcels::Float64
end

## CHECKCAPACITY: function to test whether the input capacity is viable for the transactional data
function CHECKCAPACITY(trans,
                       capacity::Array{Int64,1})
    valid = 1
    if sum(capacity) < size(trans,2)
        error("Abort due to insufficient capacity!")
        valid = 0
    end
    return valid::Int64
end

# function to apply a pair-wise exchange local search on the allocation
# of all other heuristics
function LOCALSEARCH(W::Array{Int64,2},
                     Q::Array{Int64,2})
    # lock to avoid racing conditions during the local search
    ren_lock = ReentrantLock()
    iteration = 1
    coapp_sort = vec(sum(Q,dims = 2))
    coapp_sort = sortperm(coapp_sort,rev=true)
    while iteration > 0
        iteration = 0
        for i in coapp_sort
            for j in coapp_sort
                if i != j
                    for k = 2:size(W,2)
                        for g = 1:size(W,2)-1
                            if W[i,k] == 1 && W[j,g] == 1 && W[i,g] == 0 && W[j,k] == 0
                                iteration += POTENTIAL!(i,j,k,g,W,Q)
                            end
                        end
                    end
                end
            end
        end
    end
   return W::Array{Int64,2}
end

# function to calculate the potential improvement during
# the local search procedure
function POTENTIAL!(i::Int64,
                        j::Int64,
                        k::Int64,
                        g::Int64,
                        X::Array{Int64,2},
                        Q::Array{Int64,2})
    potential = 0
    impr = 0
    potential += dot(@view(X[:,g]),@view(Q[i,:]))
    potential -= dot(@view(X[:,k]),@view(Q[i,:]))
    potential += dot(@view(X[:,k]),@view(Q[j,:]))
    potential -= dot(@view(X[:,g]),@view(Q[j,:]))
    if potential > 0
        impr = 1
        X[i,k] = 0
        X[j,g] = 0
        X[i,g] = 1
        X[j,k] = 1
    end
    return impr::Int64
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
"""
    RunIf(f, rule)

`RunIf`s allows wrapping a rule in a condition, passed the `SimData` object and the cell state and index.

```julia
RunIf(dispersal) do data, state, I
    state >= oneunit(state)
end
`` `
"""
struct RunIf{R,W,F,T<:Rule{R,W}} <: RuleWrapper{R,W}
    f::F
    rule::T
end

rule(runif::RunIf) = runif.rule
# Forward ruletype, radius and neighborhoodkey to the contained rule
ruletype(runif::RunIf) = ruletype(rule(runif))
radius(runif::RunIf) = radius(rule(runif))
neighborhoodkey(runif::RunIf) = neighborhoodkey(rule(runif))
neighborhood(runif::RunIf) = neighborhood(rule(runif))
neighbors(runif::RunIf) = neighbors(rule(runif))

@inline function setwindow(runif::RunIf{R,W}, win) where {R,W}
    f = runif.f
    r = setwindow(rule(runif), win)
    RunIf{R,W,typeof(f),typeof(r)}(f, r)
end

# We have to hook into cell_kernel! to handle the option of no return value
@inline function cell_kernel!(
    simdata, ruletype::Val{<:Rule}, condition::RunIf, rkeys, wkeys, I...
)
    readval = _readcell(simdata, rkeys, I...)
    if condition.f(simdata, readval, I)
        writeval = applyrule(simdata, rule(condition), readval, I)
        _writecell!(simdata, ruletype, wkeys, writeval, I...)
    else
        # Otherwise copy source to dest without change
        _writecell!(simdata, ruletype, wkeys, _readcell(simdata, wkeys, I...), I...)
    end
    return nothing
end
# We have to hook into cell_kernel! to handle the option of no return value
@inline function cell_kernel!(
    simdata, ::Type{<:SetRule}, condition::RunIf, rkeys, wkeys, I...
)
    readval = _readcell(simdata, rkeys, I...)
    if condition.f(data, readval, I)
        applyrule!(simdata, rule(condition), readval, I)
    end
    return nothing
end

"""
    RunAt(rules...)
    RunAt(rules::Tuple)

`RunAt`s allow running a `Rule` or multiple `Rule`s at a lower frequeny
than the main simulation, using a `range` matching the main `tspan` but with a larger
span, or specific events - by using a vector of arbitrary times in `tspan`.
"""
struct RunAt{R,W,Ru<:Tuple,Ti<:AbstractVector} <: RuleWrapper{R,W}
    rules::Ru
    times::Ti
end
RunAt(rules...; times) = RunAt(rules, times)
RunAt(rules::Tuple; times) = RunAt(rules, times)
function RunAt(rules::Tuple, times)
    rkeys = Tuple{union(map(k -> _asiterable(_readkeys(k)), rules)...)...}
    wkeys = Tuple{union(map(k -> _asiterable(_writekeys(k)), rules)...)...}
    RunAt{rkeys,wkeys,typeof(rules),typeof(times)}(rules, times)
end

rules(runat::RunAt) = runat.rules
# Only the first rule in runat can be a NeighborhoodRule, but this seems annoying...
radius(runat::RunAt) = radius(first(rules(runat)))
# Forward ruletype to the contained rule
ruletype(runat::RunAt) = ruletype(first(rules(runat)))
neighborhoodkey(runat::RunAt) = neighborhoodkey(first(rules(runat)))
neighborhood(runat::RunAt) = neighborhood(neighborhood(rules(runat)))
neighbors(runat::RunAt) = neighbors(rules(runat))

function sequencerules!(simdata::AbstractSimData, rules::Tuple{<:RunAt,Vararg})
    runat = rules[1]
    if currenttime(simdata) in runat.times
        # Run the sequenced rule
        simdata = sequencerules!(simdata, DynamicGrids.rules(runat))
    end
    # Run the rest of the rules recursively
    sequencerules!(simdata, tail(rules))
end

function Base.tail(runat::RunAt{R,W}) where {R,W}
    rulestail = tail(rules(runat))
    RunAt{R,W,typeof(rulestail),typeof(runat.times)}(rulestail, runat.times)
end
Base.getindex(runat::RunAt, i) = getindex(rules(runat), i)
Base.iterate(runat::RunAt) = iterate(rules(runat))
Base.iterate(runat::RunAt, nothing) = iterate(rules(runat), nothing)
Base.length(runat::RunAt) = length(rules(runat))
Base.firstindex(runat::RunAt) = firstindex(rules(runat))
Base.lastindex(runat::RunAt) = lastindex(rules(runat))