text:
module ParallelIncrements

export atomicref

using Base: gc_alignment, llvmcall
using Base.Threads: ArithmeticTypes, WORD_SIZE, atomictypes, inttype, llvmtypes
using BenchmarkTools

struct AtomicRef{T,A} <: Ref{T}
    pointer::Ptr{T}
    x::A
end

atomicref(x, args...) = AtomicRef(pointer(x, args...), x)

include("atomics.jl")

function seq_unsafe_incrementat!(dest, indices)
    @inbounds for i in indices
        dest[i] += one(eltype(dest))
    end
    return dest
end

function unsafe_incrementat!(dest, indices)
    for i in indices
        r = atomicref(dest, i)
        Threads.atomic_add!(r, one(eltype(dest)))
    end
    return dest
end

function benchsuite()
    suite = BenchmarkGroup()
    let s1 = suite["single"] = BenchmarkGroup()
        s1["atomic"] = @benchmarkable(
            begin
                fill!(dest, 0)
                for _ in 1:1000
                    r = atomicref(dest, 1)
                    Threads.atomic_add!(r, one(eltype(dest)))
                end
            end,
            setup = (dest = [0]),
        )
        s1["nonatomic"] = @benchmarkable(begin
            fill!(dest, 0)
            for _ in 1:1000
                @inbounds dest[1] += 1
            end
        end, setup = (dest = [0]))
    end
    for n in [1_000]
        s1 = suite["n=$n"] = BenchmarkGroup()
        dest = zeros(Int, n)
        for m in [1_000_000]
            indices = rand(eachindex(dest), m)
            s2 = s1["m=$m"] = BenchmarkGroup()
            s2["nonatomic"] = @benchmarkable(
                seq_unsafe_incrementat!($dest, $indices),
                setup = (fill!($dest, 0))
            )
            s2["atomic"] = @benchmarkable(
                unsafe_incrementat!($dest, $indices),
                setup = (fill!($dest, 0))
            )
        end
    end
    return suite
end

end # module


text:
module FinEtoolsFlexStructures

__precompile__(true)

using FinEtools

include("RotUtilModule.jl")
include("TransformerModule.jl")
include("AssemblyModule.jl")

include("CrossSectionModule.jl")
include("CompositeLayupModule.jl")

include("FESetCorotBeamModule.jl")
include("FESetShellQ4Module.jl")
include("FESetShellT3Module.jl")

include("MeshFrameMemberModule.jl")

include("FEMMCorotBeamModule.jl")

include("FEMMPointMassModule.jl")

include("FEMMPointGroundedSpringModule.jl")

include("FEMMShellT3FFModule.jl")
include("FEMMShellT3FFCompModule.jl")

end # module

ext: jl


text:
module HashDictionaries

using Dictionaries
using Base: @propagate_inbounds

export HashIndices, HashDictionary

# These can be changed, to trade off better performance for space
const global maxallowedprobe = 16
const global maxprobeshift   = 6

mutable struct HashIndices{T} <: AbstractIndices{T}
    slots::Array{UInt8,1}
    inds::Array{T,1}
    ndel::Int
    count::Int
    idxfloor::Int  # an index <= the indices of all used slots
    maxprobe::Int
end

HashIndices() = HashIndices{Any}()

"""
    HashIndices{I}()

Construct an empty `HashIndices` with indices of type `I`. This container uses hashes for
fast lookup, and is insertable. (See `isinsertable`).
"""
function HashIndices{T}(; sizehint::Int = 16) where {T}
    sz = Base._tablesz(sizehint)
    HashIndices{T}(zeros(UInt8, sz), Vector{T}(undef, sz), 0, 0, 1, 0)
end

## Constructors

"""
    HashIndices(iter)
    HashIndices{I}(iter)

Construct a `HashIndices` with indices from iterable container `iter`.
"""
function HashIndices(iter)
    if Base.IteratorEltype(iter) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        iter = collect(iter)
    end

    return HashIndices{eltype(iter)}(iter)
end

function HashIndices{T}(iter) where {T}
    iter_size = Base.IteratorSize(iter)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        h = HashIndices{T}(; sizehint = length(iter)*2)
    else
        h = HashIndices{T}()
    end

    for i in iter
        insert!(h, i) # should this be `set!` or `insert!`?
    end

    return h
end

function Base.copy(h::HashIndices{T}, ::Type{T}) where {T}
    return HashIndices{T}(copy(h.slots), copy(h.inds), h.ndel, h.count, h.idxfloor, h.maxprobe)
end

## Length
Base.length(h::HashIndices) = h.count

## Token interface

Dictionaries.istokenizable(::HashIndices) = true
Dictionaries.tokentype(::HashIndices) = Int

@propagate_inbounds isslotempty(h::HashIndices, i::Int) = h.slots[i] == 0x0
@propagate_inbounds isslotfilled(h::HashIndices, i::Int) = h.slots[i] == 0x1
@propagate_inbounds isslotdeleted(h::HashIndices, i::Int) = h.slots[i] == 0x2 # deletion marker/tombstone

istokenassigned(h::HashIndices, i::Int) = isslotfilled(h, i)

# iteratetoken

function skip_deleted(h::HashIndices, i)
    L = length(h.slots)
    @inbounds while i <= L && !isslotfilled(h, i)
        i += 1
    end
    return i
end

@propagate_inbounds function Dictionaries.iteratetoken(h::HashIndices{T}) where {T}
    idx = skip_deleted(h, h.idxfloor)
    h.idxfloor = idx # An optimization to skip unnecessary elements when iterating multiple times

    if idx > length(h.inds)
        return nothing
    else
        return (idx, idx + 1)
    end
end

@propagate_inbounds function Dictionaries.iteratetoken(h::HashIndices{T}, idx::Int) where {T}
    idx = skip_deleted(h, idx)

    if idx > length(h.inds)
        return nothing
    else
        return (idx, idx + 1)
    end
end

# gettoken
function hashtoken(key, sz::Int)
    # Given key what is the hash slot? sz is a power of two
    (((hash(key)%Int) & (sz-1)) + 1)::Int
end

function Dictionaries.gettoken(h::HashIndices{T}, key::T) where {T}
    sz = length(h.inds)
    iter = 0
    maxprobe = h.maxprobe
    token = hashtoken(key, sz)
    keys = h.inds

    @inbounds while true
        if isslotempty(h, token)
            break
        end
        if !isslotdeleted(h, token) && (key === keys[token] || isequal(key, keys[token]))
            return (true, token)
        end

        token = (token & (sz-1)) + 1
        iter += 1
        iter > maxprobe && break
    end
    return (false, 0)
end

# gettokenvalue
@propagate_inbounds function Dictionaries.gettokenvalue(h::HashIndices, token::Int)
    return h.inds[token]
end

# insertable interface
Dictionaries.isinsertable(::HashIndices) = true

function Base.empty!(h::HashIndices{T}) where {T}
    fill!(h.slots, 0x0) # It should be OK to reduce this back to some smaller size.
    sz = length(h.slots)
    empty!(h.inds)
    resize!(h.inds, sz)
    h.ndel = 0
    h.count = 0
    h.idxfloor = 1
    return h
end

function rehash!(h::HashIndices, newsz::Int = length(h.inds))
    _rehash!(h, nothing, newsz)
    return h
end

function _rehash!(h::HashIndices{T}, oldv::Union{Nothing, Vector}, newsz::Int) where {T}
    olds = h.slots
    oldk = h.inds
    sz = length(olds)
    newsz = Base._tablesz(newsz)
    h.idxfloor = 1
    if h.count == 0
        resize!(h.slots, newsz)
        fill!(h.slots, 0)
        resize!(h.inds, newsz)
        error()
        oldv === nothing || resize!(oldv, newsz)
        h.ndel = 0
        return oldv
    end

    slots = zeros(UInt8, newsz)
    keys = Vector{T}(undef, newsz)
    vals = oldv === nothing ? nothing : Vector{eltype(oldv)}(undef, newsz)
    count = 0
    maxprobe = h.maxprobe

    for i ∈ 1:sz
        @inbounds if olds[i] == 0x1
            k = oldk[i]
            v = vals === nothing ? nothing : oldv[i]
            index0 = index = hashtoken(k, newsz)
            while slots[index] != 0
                index = (index & (newsz-1)) + 1
            end
            probe = (index - index0) & (newsz-1)
            probe > maxprobe && (maxprobe = probe)
            slots[index] = 0x1
            keys[index] = k
            vals === nothing || (vals[index] = v)
            count += 1
        end
    end

    h.slots = slots
    h.inds = keys
    h.count = count
    h.ndel = 0
    h.maxprobe = maxprobe

    return vals
end

Base.sizehint!(h::HashIndices, newsz::Int) = _sizehint!(h, nothing, newsz)

function _sizehint!(h::HashIndices{T}, values::Union{Nothing, Vector}, newsz::Int) where {T}
    oldsz = length(h.slots)
    if newsz <= oldsz
        # TODO: shrink
        # be careful: rehash!() assumes everything fits. it was only designed
        # for growing.
        return hash
    end
    # grow at least 25%
    newsz = min(max(newsz, (oldsz*5)>>2),
                Base.max_values(T))
    return _rehash!(h, values, newsz)
end

function Dictionaries.gettoken!(h::HashIndices{T}, key::T) where {T}
    (token, _) = _gettoken!(h, nothing, key) # This will make sure a slot is available at `token` (or `-token` if it is new)

    if token < 0
        @inbounds (token, _) = _insert!(h, nothing, key, -token) # This will fill the slot with `key`
        return (false, token)
    else
        return (true, token)
    end
end

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by insert!, set! and get!
function _gettoken!(h::HashIndices{T}, values::Union{Nothing, Vector}, key::T) where {T}
    sz = length(h.inds)
    iter = 0
    maxprobe = h.maxprobe
    token = hashtoken(key, sz)
    avail = 0
    keys = h.inds

    # Search of the key is present or if there is a deleted slot `key` could fill.
    @inbounds while true
        if isslotempty(h, token)
            if avail < 0
                return (avail, values)
            end
            return (-token, values)
        end

        if isslotdeleted(h, token)
            if avail == 0
                # found an available deleted slot, but we need to keep scanning
                # in case `key` already exists in a later collided slot.
                avail = -token
            end
        elseif key === keys[token] || isequal(key, keys[token])
            return (token, values)
        end

        token = (token & (sz-1)) + 1
        iter += 1
        iter > maxprobe && break
    end

    avail < 0 && return (avail, values)

    # The key definitely isn't present, but a slot may become available if we increase
    # `maxprobe` (up to some reasonable global limits).
    maxallowed = max(maxallowedprobe, sz>>maxprobeshift)

    @inbounds while iter < maxallowed
        if !isslotfilled(h,token)
            h.maxprobe = iter
            return (-token, values)
        end
        token = (token & (sz-1)) + 1
        iter += 1
    end

    # If we get here, then all the probable slots are filled, and the only recourse is to
    # increase the size of the hash map and try again
    values = _rehash!(h, values, h.count > 64000 ? sz*2 : sz*4)
    return _gettoken!(h, values, key)
end

@propagate_inbounds function _insert!(h::HashIndices{T}, values::Union{Nothing, Vector}, key::T, token::Int) where {T}
    h.slots[token] = 0x1
    h.inds[token] = key
    h.count += 1
    if token < h.idxfloor
        h.idxfloor = token
    end

    # TODO revisit this...
    #=
    sz = length(h.inds)
    # Rehash now if necessary
    if h.ndel >= ((3*sz)>>2) || h.count*3 > sz*2
        # > 3/4 deleted or > 2/3 full
        values = _rehash!(h, values, h.count > 64000 ? h.count*2 : h.count*4)
        (_, token) = gettoken(h, key)
    end
    =#

    return (token, values)
end

function Dictionaries.deletetoken!(h::HashIndices{T}, token::Int) where {T}
    h.slots[token] = 0x2
    isbitstype(T) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.inds, token-1)

    h.ndel += 1
    h.count -= 1
    return h
end

# Since deleting elements doesn't mess with iteration, we can use `unsafe_filter!``
Base.filter!(pred, h::HashIndices) = Base.unsafe_filter!(pred, h)

# The default insertable indices
Base.empty(d::HashIndices, ::Type{T}) where {T} = HashIndices{T}()

mutable struct HashDictionary{I,T} <: AbstractDictionary{I, T}
    indices::HashIndices{I}
    values::Vector{T}

    HashDictionary{I, T}(indices::HashIndices{I}, values::Vector{T}, ::Nothing) where {I, T} = new(indices, values)
end

"""
    HashDictionary{I, T}()

Construct an empty `HashDictionary` with index type `I` and element type `T`. This type of
dictionary uses hashes for fast lookup and insertion, and is both mutable and insertable.
(See `issettable` and `isinsertable`).
"""
function HashDictionary{I, T}(; sizehint::Int = 16) where {I, T}
    indices = HashIndices{I}(; sizehint=sizehint)
    HashDictionary{I, T}(indices, Vector{T}(undef, length(indices.slots)), nothing)
end
HashDictionary{I}() where {I} = HashDictionary{I, Any}()
HashDictionary() = HashDictionary{Any}()

"""
    HashDictionary{I, T}(indices, undef::UndefInitializer)

Construct a `HashDictionary` with index type `I` and element type `T`. The container is
initialized with `keys` that match the values of `indices`, but the values are unintialized.
"""
function HashDictionary{I, T}(indices, ::UndefInitializer) where {I, T}
    return HashDictionary{I, T}(HashIndices{I}(indices), undef)
end

function HashDictionary{I, T}(h::HashIndices{I}, ::UndefInitializer) where {I, T}
    return HashDictionary{I, T}(h, Vector{T}(undef, length(h.slots)), nothing)
end

function HashDictionary{I, T}(indices::HashIndices{I}, values) where {I, T}
    vals = Vector{T}(undef, length(indices.slots))
    d = HashDictionary{I, T}(indices, vals, nothing)

    @inbounds for (i, v) in zip(tokens(indices), values)
        vals[i] = v
    end

    return d
end

"""
    HashDictionary(indices, values)
    HashDictionary{I}(indices, values)
    HashDictionary{I, T}(indices, values)

Construct a `HashDictionary` with indices from `indices` and values from `values`, matched
in iteration order.
"""
function HashDictionary{I, T}(indices, values) where {I, T}
    iter_size = Base.IteratorSize(indices)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        d = HashDictionary{I, T}(; sizehint = length(indices)*2)
    else
        d = HashDictionary{I, T}()
    end

    for (i, v) in zip(indices, values)
        insert!(d, i, v)
    end

    return d
end
function HashDictionary{I}(indices, values) where {I}
    if Base.IteratorEltype(values) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        values = collect(values)
    end

    return HashDictionary{I, eltype(values)}(indices, values)
end

function HashDictionary(indices, values)
    if Base.IteratorEltype(indices) === Base.EltypeUnknown()
        # TODO: implement automatic widening from iterators of Base.EltypeUnkown
        indices = collect(indices)
    end

    return HashDictionary{eltype(indices)}(indices, values)
end

"""
    HashDictionary(dict::AbstractDictionary)
    HashDictionary{I}(dict::AbstractDictionary)
    HashDictionary{I, T}(dict::AbstractDictionary)

Construct a copy of `dict` with the same keys and values.
(For copying an `AbstractDict` or other iterable of `Pair`s, see `dictionary`).
"""
HashDictionary(dict::AbstractDictionary) = HashDictionary(keys(dict), dict)
HashDictionary{I}(dict::AbstractDictionary) where {I} = HashDictionary{I}(keys(dict), dict)
HashDictionary{I, T}(dict::AbstractDictionary) where {I, T} = HashDictionary{I, T}(keys(dict), dict)

## Implementation

Base.keys(d::HashDictionary) = d.indices
Dictionaries.isinsertable(d::HashDictionary) = true
Dictionaries.issettable(d::HashDictionary) = true

@propagate_inbounds function Dictionaries.gettoken(d::HashDictionary{I}, i::I) where {I}
    return gettoken(keys(d), i)
end

@inline function Dictionaries.gettokenvalue(d::HashDictionary, token)
    return @inbounds d.values[token]
end

function Dictionaries.istokenassigned(d::HashDictionary, token)
    return isassigned(d.values, token)
end

@inline function Dictionaries.settokenvalue!(d::HashDictionary{I, T}, token, value::T) where {I, T}
    @inbounds d.values[token] = value
    return d
end

function Dictionaries.gettoken!(d::HashDictionary{T}, key::T) where {T}
    indices = keys(d)
    (token, values) = _gettoken!(indices, d.values, key)
    if token < 0
        (token, values) = _insert!(indices, values, key, -token)
        d.values = values
        return (false, token)
    else
        d.values = values
        return (true, token)
    end
end

function Base.copy(d::HashDictionary{I, T}, ::Type{I}, ::Type{T}) where {I, T}
    return HashDictionary{I, T}(d.indices, copy(d.values), nothing)
end

Dictionaries.tokenized(d::HashDictionary) = d.values

function Base.empty!(d::HashDictionary)
    empty!(d.indices)
    empty!(d.values)
    resize!(d.values, length(keys(d).slots))
    return d
end

function Dictionaries.deletetoken!(d::HashDictionary{I, T}, token) where {I, T}
    deletetoken!(keys(d), token)
    isbitstype(T) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), d.values, token-1)
    return d
end

function Base.sizehint!(d::HashDictionary, sz::Int)
    d.values = _sizehint!(d.indices, d.values, sz)
    return d
end

function Base.rehash!(d::HashDictionary, newsz::Int = length(d.indices))
    _rehash!(d.indices, d.values, newsz)
    return d
end

Base.filter!(pred, d::HashDictionary) = Base.unsafe_filter!(pred, d)

# For `HashIndices` we don't copy the indices, we allow the `keys` to remain identical (`===`)
function Base.similar(indices::HashIndices{I}, ::Type{T}) where {I, T}
    return HashDictionary{I, T}(indices, undef)
end

function Base.empty(indices::HashIndices, ::Type{I}, ::Type{T}) where {I, T}
    return HashDictionary{I, T}()
end

end # module


text:
using Senuk

X = Matrix{Float64}([12 22; 16 14; 11 19])
y = Vector{Float64}([34; 30; 30])
lasso = Regression.Lasso(1.3)
model = Regression.create(X, y, step_size = 0.0001, max_iter = 500,
                          verbose = true, regularization = lasso)
y_predict = Regression.predict(model, X)
@printf("predicted responses: %s\n", y_predict)
@printf("actual responses: %s\n", "[34; 30; 30]")

max_stars_repo_name: upul/Senuk

max_issues_repo_name: upul/Senuk

max_forks_repo_name: upul/Senuk


text:
import EnergyModels
const EM = EnergyModels

using Test

@testset "EnergyModels" begin
    # @test isempty(detect_ambiguities(EM, Base, Core))

    @testset "WrappedArray" begin
        include("wrappedarray.jl")
    end
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
# ---
# title: 507. Perfect Number
# id: problem507
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/perfect-number/description/>
# hidden: true
# ---
#
# A [**perfect number**](https://en.wikipedia.org/wiki/Perfect_number) is a
# **positive integer** that is equal to the sum of its **positive divisors** ,
# excluding the number itself. A **divisor** of an integer `x` is an integer
# that can divide `x` evenly.
#
# Given an integer `n`, return `true` _if_`n` _is a perfect number, otherwise
# return_`false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: num = 28
#     Output: true
#     Explanation: 28 = 1 + 2 + 4 + 7 + 14
#     1, 2, 4, 7, and 14 are all divisors of 28.
#
#
# **Example 2:**
#
#
#
#     Input: num = 6
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: num = 496
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: num = 8128
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: num = 2
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= num <= 108`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsTeamProfileField(;
        hint=nothing,
        id=nothing,
        is_hidden=nothing,
        label=nothing,
        options=nothing,
        ordering=nothing,
        type=nothing,
    )

    - hint::String
    - id::String
    - is_hidden::Bool
    - label::String
    - options::Any
    - ordering::Float32
    - type::String
"""
mutable struct ObjsTeamProfileField <: SwaggerModel
    hint::Any # spec type: Union{ Nothing, String } # spec name: hint
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    is_hidden::Any # spec type: Union{ Nothing, Bool } # spec name: is_hidden
    label::Any # spec type: Union{ Nothing, String } # spec name: label
    options::Any # spec type: Union{ Nothing, Any } # spec name: options
    ordering::Any # spec type: Union{ Nothing, Float32 } # spec name: ordering
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function ObjsTeamProfileField(;hint=nothing, id=nothing, is_hidden=nothing, label=nothing, options=nothing, ordering=nothing, type=nothing)
        o = new()
        validate_property(ObjsTeamProfileField, Symbol("hint"), hint)
        setfield!(o, Symbol("hint"), hint)
        validate_property(ObjsTeamProfileField, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsTeamProfileField, Symbol("is_hidden"), is_hidden)
        setfield!(o, Symbol("is_hidden"), is_hidden)
        validate_property(ObjsTeamProfileField, Symbol("label"), label)
        setfield!(o, Symbol("label"), label)
        validate_property(ObjsTeamProfileField, Symbol("options"), options)
        setfield!(o, Symbol("options"), options)
        validate_property(ObjsTeamProfileField, Symbol("ordering"), ordering)
        setfield!(o, Symbol("ordering"), ordering)
        validate_property(ObjsTeamProfileField, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type ObjsTeamProfileField

const _property_map_ObjsTeamProfileField = Dict{Symbol,Symbol}(Symbol("hint")=>Symbol("hint"), Symbol("id")=>Symbol("id"), Symbol("is_hidden")=>Symbol("is_hidden"), Symbol("label")=>Symbol("label"), Symbol("options")=>Symbol("options"), Symbol("ordering")=>Symbol("ordering"), Symbol("type")=>Symbol("type"))
const _property_types_ObjsTeamProfileField = Dict{Symbol,String}(Symbol("hint")=>"String", Symbol("id")=>"String", Symbol("is_hidden")=>"Bool", Symbol("label")=>"String", Symbol("options")=>"Any", Symbol("ordering")=>"Float32", Symbol("type")=>"String")
Base.propertynames(::Type{ ObjsTeamProfileField }) = collect(keys(_property_map_ObjsTeamProfileField))
Swagger.property_type(::Type{ ObjsTeamProfileField }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsTeamProfileField[name]))}
Swagger.field_name(::Type{ ObjsTeamProfileField }, property_name::Symbol) =  _property_map_ObjsTeamProfileField[property_name]

const _allowed_ObjsTeamProfileField_type = ["text", "date", "link", "mailto", "options_list", "user"]

function check_required(o::ObjsTeamProfileField)
    (getproperty(o, Symbol("hint")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("label")) === nothing) && (return false)
    (getproperty(o, Symbol("ordering")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsTeamProfileField }, name::Symbol, val)
    if name === Symbol("id")
    end
    if name === Symbol("type")
        Swagger.validate_param(name, "ObjsTeamProfileField", :enum, val, _allowed_ObjsTeamProfileField_type)
    end
end


text:
print("Testing Parameters...")
took_seconds = @elapsed include("ParametersTests.jl")
println("done (took ", took_seconds, " seconds)")

print("Testing Trial/TrialEstimate/TrialRatio/TrialJudgement...")
took_seconds = @elapsed include("TrialsTests.jl")
println("done (took ", took_seconds, " seconds)")

print("Testing BenchmarkGroup...")
took_seconds = @elapsed include("GroupsTests.jl")
println("done (took ", took_seconds, " seconds)")

print("Testing execution...")
took_seconds = @elapsed include("ExecutionTests.jl")
println("done (took ", took_seconds, " seconds)")

print("Testing serialization...")
took_seconds = @elapsed include("SerializationTests.jl")
println("done (took ", took_seconds, " seconds)")


text:
include("filtering.jl");
##
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
s = filter(r->r.Treatment in list &&
    r.Trial < 51 &&
    # r.Trial_Travel_to < 40 &&
    r.MouseID != "pc7",
    streaks)
s[s.Treatment .== "PreVehicle",:Treatment] .= "Control"
s[s.Treatment .== "Saline",:Treatment] = [o ? "Optogenetic" : "Control" for o in s[s.Treatment .== "Saline",:Stim]]
s[s.Treatment .== "SB242084_opt",:Phase] .=  "SB242084_opt"
s[s.Treatment .== "SB242084_opt",:Treatment] = [o ? "SB242084_opt" : "Control" for o in s[s.Treatment .== "SB242084_opt",:Stim]]
# Df = combine([:Treatment,:Stim_Day] => (t,s) -> (treatment = union(t)), groupby(s,:Phase))
######################## All  Manipulations ###################################
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Control, :Drug; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/WilcoxonNumPokes50.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/WilcoxonTravelROI50.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/WilcoxonTravelPOKES50.pdf"))
##
df1 = combine(groupby(s,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/WilcoxonAfterLast50.pdf"))
######################## 5HT Release  Manipulations ###################################
release_man =  filter(r -> r.Phase in ["Citalopram", "Optogenetic"], s)
df1 = combine(groupby(release_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/ReleaseNumPokes50.pdf"))
##
df1 = combine(groupby(release_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/ReleaseTravelROI50.pdf"))
##
df1 = combine(groupby(release_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/ReleaseTravelPOKEES50.pdf"))
##
df1 = combine(groupby(release_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/ReleaseAfterLast50.pdf"))
######################## Global  Manipulations ###################################
glob_man =  filter(r -> r.Phase in ["Citalopram", "Optogenetic", "Methysergide"], s)
df1 = combine(groupby(glob_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/GlobalNumPokes50.pdf"))
##
df1 = combine(groupby(glob_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/GlobalTravelROI50.pdf"))
##
df1 = combine(groupby(glob_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/GlobalTravelPOKES50.pdf"))
##
df1 = combine(groupby(glob_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/GlobalAfterLast50.pdf"))
######################## Selective  Manipulations ###################################
sel_man =  filter(r -> r.Phase in ["SB242084", "Altanserin", "Way_100135"], s)
df1 = combine(groupby(sel_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Num_pokes)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/SelectiveNumPokes50.pdf"))
##
df1 = combine(groupby(sel_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Trial_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/SelectiveTravelROI50.pdf"))
##
df1 = combine(groupby(sel_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:Poking_Travel_to)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/SelectiveTravelPOKES50.pdf"))
##
df1 = combine(groupby(sel_man,:Phase)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    current_drug = Symbol(subdf[1,:Phase])
    rename!(subdf, current_drug => :Drug)
end
df2 = combine(groupby(df1,[:Phase])) do dd
    wilcoxon(dd,:Drug, :Control; f = x -> mean(skipmissing(x)))
end
Drug_colors!(df2)
plot_wilcoxon(df2)
##
savefig(joinpath(figs_loc,"LabMeeting/SelectiveAfterLast50.pdf"))
##
x = collect(1:100)
y = -log.(10,x)
plot(x,y)


text:
using LinearAlgebra

include("byrow.jl")
include("dare.jl")
include("dlyap.jl")
include("simula01.jl")


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((1, 2, true, true),),
    ((1, 2, false, true),),
    ((1, 2, false, false),),
    ((:lotkavolterra, 3, true, true),),
    ((:motzkin_ellipsoid, 3, false, true),),
    ((:schwefel_ball, 2, false, false),),
    ]
insts["fast"] = [
    ((1, 3, true, true),),
    ((1, 30, true, true),),
    ((1, 30, false, true),),
    ((1, 30, false, false),),
    ((2, 8, true, true),),
    ((3, 6, true, true),),
    ((5, 3, true, true),),
    ((10, 1, true, true),),
    ((10, 1, false, true),),
    ((10, 1, false, false),),
    ((4, 4, true, true),),
    ((4, 4, false, true),),
    ((4, 4, false, false),),
    ((:butcher, 2, true, true),),
    ((:caprasse, 4, true, true),),
    ((:goldsteinprice, 7, true, true),),
    ((:goldsteinprice_ball, 6, true, true),),
    ((:goldsteinprice_ellipsoid, 7, true, true),),
    ((:heart, 2, true, true),),
    ((:lotkavolterra, 3, true, true),),
    ((:magnetism7, 2, true, true),),
    ((:magnetism7_ball, 2, true, true),),
    ((:motzkin, 3, true, true),),
    ((:motzkin_ball, 3, true, true),),
    ((:motzkin_ellipsoid, 3, true, true),),
    ((:reactiondiffusion, 4, true, true),),
    ((:robinson, 8, true, true),),
    ((:robinson_ball, 8, true, true),),
    ((:rosenbrock, 5, true, true),),
    ((:rosenbrock_ball, 5, true, true),),
    ((:schwefel, 2, true, true),),
    ((:schwefel_ball, 2, true, true),),
    ((:lotkavolterra, 3, false, true),),
    ((:motzkin, 3, false, true),),
    ((:motzkin_ball, 3, false, true),),
    ((:schwefel, 2, false, true),),
    ((:lotkavolterra, 3, false, false),),
    ((:motzkin, 3, false, false),),
    ((:motzkin_ball, 3, false, false),),
    ]
insts["various"] = [
    ((4, 5, true, true),),
    ((3, 5, true, false),),
    ((4, 5, false, true),),
    ((4, 5, false, false),),
    ((:butcher, 2, true, true),),
    ((:caprasse, 4, true, true),),
    ((:goldsteinprice, 7, true, true),),
    ((:heart, 2, true, true),),
    ((:lotkavolterra, 3, true, true),),
    ((:magnetism7, 2, true, true),),
    ((:motzkin, 3, true, true),),
    ((:reactiondiffusion, 4, true, true),),
    ((:robinson, 8, true, true),),
    ((:rosenbrock, 5, true, true),),
    ((:schwefel, 2, true, true),),
    ]
return (PolyMinJuMP, insts)


text:
using BinaryBuilder

name = "Openlibm"
version = v"0.6.0"
sources = [
    "https://github.com/JuliaMath/openlibm/archive/v0.6.0.tar.gz" =>
    "d45439093d1fd15e2ac3acf69955e462401c7a160d3330256cb4a86c51bdae28",
]

script = raw"""
# Enter the funzone
cd ${WORKSPACE}/srcdir/openlibm-*

# Install into output
flags=("prefix=${prefix}")

# Build ARCH from ${target}
flags+=("ARCH=${target%-*-*}")

# Openlibm build system doesn't recognize our windows cross compilers properly
if [[ ${target} == *mingw* ]]; then
    flags+=("OS=WINNT")
fi

# Add `CC` override, since Openlibm seems to think it knows best:
flags+=("CC=$CC")

# Build the library
make "${flags[@]}" -j${nproc}

# Install the library
make "${flags[@]}" install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms()
platforms = expand_gcc_versions(platforms)

products = prefix -> [
    LibraryProduct(prefix, "libopenlibm", :libopenlibm)
]

dependencies = [
]

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


text:
"""
"""
function DiscreteModel(parent_model::DiscreteModel,cell_to_parent_cell::AbstractVector{<:Integer})
  RestrictedDiscreteModel(parent_model,cell_to_parent_cell)
end

function DiscreteModel(parent_model::DiscreteModel,parent_cell_to_mask::AbstractArray{Bool})
  RestrictedDiscreteModel(parent_model,parent_cell_to_mask)
end

function DiscreteModel(parent_model::DiscreteModel,labels::FaceLabeling;tags)
  parent_cell_to_mask = get_face_mask(labels,tags,num_cell_dims(parent_model))
  DiscreteModel(parent_model,parent_cell_to_mask)
end

function DiscreteModel(parent_model::DiscreteModel;tags)
  labels = get_face_labeling(parent_model)
  DiscreteModel(parent_model,labels,tags=tags)
end

"""
"""
struct RestrictedDiscreteModel{Dc,Dp} <: DiscreteModel{Dc,Dc}
  model::DiscreteModelPortion{Dc,Dp}
end

function RestrictedDiscreteModel(parent_model::DiscreteModel, cell_to_parent_cell::AbstractVector{<:Integer})
  _model = DiscreteModelPortion(parent_model,cell_to_parent_cell)
  RestrictedDiscreteModel(_model)
end

function RestrictedDiscreteModel(parent_model::DiscreteModel, parent_cell_to_mask::AbstractArray{Bool})
  _model = DiscreteModelPortion(parent_model,parent_cell_to_mask)
  RestrictedDiscreteModel(_model)
end

get_grid(model::RestrictedDiscreteModel) = get_grid(model.model)

get_grid_topology(model::RestrictedDiscreteModel) = get_grid_topology(model.model)

get_face_labeling(model::RestrictedDiscreteModel) = get_face_labeling(model.model)

get_face_to_parent_face(model::RestrictedDiscreteModel,d::Integer) = get_face_to_parent_face(model.model,d)

get_cell_to_parent_cell(model::RestrictedDiscreteModel) = get_cell_to_parent_cell(model.model)

get_parent_model(model::RestrictedDiscreteModel) = get_parent_model(model.model)

function Triangulation(model::RestrictedDiscreteModel)
  parent_trian = Triangulation(get_parent_model(model))
  cell_to_parent_cell = get_cell_to_parent_cell(model)
  RestrictedTriangulation(parent_trian,cell_to_parent_cell)
end

function get_triangulation(model::RestrictedDiscreteModel)
  Triangulation(model)
end

function Triangulation(::Type{ReferenceFE{d}},model::RestrictedDiscreteModel) where d
  @notimplemented
end

function Triangulation(model::RestrictedDiscreteModel,cell_to_parent_cell::AbstractVector{<:Integer})
  @notimplemented
end

function Triangulation(model::RestrictedDiscreteModel,cell_to_mask::AbstractArray{Bool})
  @notimplemented
end

function BoundaryTriangulation(
  model::RestrictedDiscreteModel,
  face_to_bgface::AbstractVector{<:Integer},
  bgface_to_lcell::AbstractVector{<:Integer})

  parent_model = get_parent_model(model)
  bgface_to_parent_face = get_face_to_parent_face(model,num_cell_dims(model)-1)
  face_to_parent_face = lazy_map(Reindex(bgface_to_parent_face),face_to_bgface)
  parent_face_to_lcell = fill(Int8(1),num_facets(parent_model))
  parent_face_to_lcell[bgface_to_parent_face] .= bgface_to_lcell

  BoundaryTriangulation(parent_model,face_to_parent_face,parent_face_to_lcell)
end

function BoundaryTriangulation(
  model::RestrictedDiscreteModel,
  bgface_to_mask::AbstractVector{Bool},
  bgface_to_lcell::AbstractVector{<:Integer})

  face_to_bgface = findall(bgface_to_mask)
  BoundaryTriangulation(model,face_to_bgface,bgface_to_lcell)
end

function InterfaceTriangulation(model_in::RestrictedDiscreteModel,model_out::RestrictedDiscreteModel)
  cells_in = get_cell_to_parent_cell(model_in)
  cells_out = get_cell_to_parent_cell(model_out)
  model = get_parent_model(model_in)
  InterfaceTriangulation(model,cells_in,cells_out)
end


text:
function send(c::AbstractConnection, msg)
    error("No send method for connection of type $(typeof(c))")
end

function logmsg(msg, level="info", data=nothing)
    if level == "error" || level == "warn"
        @warn(msg)
    else
        @info(msg)
    end

    Dict("type"=>"log", "message"=>msg, "level"=>level, data=>data)
end

function log(c::AbstractConnection, msg, level="info", data=nothing)
    send(c, logmsg(msg, level, data))
end

function dispatch(conn::AbstractConnection, data)
    # first, check if the message is one of the administrative ones
    cmd = data["command"]
    scopeid = data["scope"]
    if cmd == "_setup_scope"
        if haskey(scopes, scopeid)
            scope = scopes[scopeid]
            addconnection!(scope.pool, conn)
        else
            log(conn, "Client says it has unknown scope $scopeid", "warn")
        end
    elseif cmd == "_acknowledge_message"
        msgid = data["messageId"]
        if haskey(waiting_messages, msgid)
            notify(waiting_messages[msgid])
            delete!(waiting_messages, msgid)
        else
            log(conn,
                "Acknowledgement received for message" *
                "$msgid but not waiting for it",
                "warn")
        end
    else
        if !haskey(scopes, scopeid)
            log(conn, "Message $data received for unknown scope $scopeid", "warn")
            return
        end
        scope = scopes[scopeid]
        dispatch(scope, cmd, data["data"])
    end
end


text:
prelu(x, α) = x > 0 ? x : α*x

"""
    EdgeFtLayer(v_in=>v_out, e_in=>e_out)

Edge Features Layers. This is used in [this article] to compute a state representation of the TSPTW.

[this article]: https://arxiv.org/abs/2006.01610

# Arguments
- `v_in`: the dimension of input features for vertices.
- `v_out`: the dimension of output features for vertices.
- `e_in`: the dimension of input features for edges.
- `e_out`: the dimension of output features for edges.
"""
struct EdgeFtLayer{T <: Real} <: MessagePassing
    W_a::AbstractMatrix{T}
    W_T::AbstractMatrix{T}
    b_T::AbstractVector{T}
    W_e::AbstractMatrix{T}
    W_ee::AbstractMatrix{T}
    prelu_α::T
end

function EdgeFtLayer(;v_dim::Pair{<:Integer,<:Integer}, e_dim::Pair{<:Integer,<:Integer}, init=glorot_uniform, T::DataType=Float32, bias::Bool=true)

    # Used to compute node features
    W_a = T.(init(v_dim[2], 2 * v_dim[1] + e_dim[1]))
    W_T = T.(init(v_dim[2], 2 * v_dim[1] + e_dim[1]))
    b_T = bias ? T.(init(v_dim[2])) : zeros(T, v_dim[2])

    # Used to compute edge features
    W_e = T.(init(e_dim[2], v_dim[1]))
    W_ee = T.(init(e_dim[2], e_dim[1]))

    EdgeFtLayer(W_a, W_T, b_T, W_e, W_ee, init(1)[1])
end

Flux.@functor EdgeFtLayer

function (g::EdgeFtLayer)(fg::FeaturedGraph)
    Zygote.ignore() do
        GraphSignals.check_num_node(graph(fg), node_feature(fg))
    end

    # TODO: when upgrading to the newer version of GeometricFlux, replace `:add` by `+`
    fg_ = GeometricFlux.propagate(g, fg, :add)
end

# TODO: when upgrading to the newer version of GeometricFlux, replace `:add` by `+`
function GeometricFlux.propagate(g::EdgeFtLayer, fg::FeaturedGraph, aggr::Symbol=:add)
    E, X, u = GeometricFlux.propagate(g, adjacency_list(fg), fg.ef, fg.nf, fg.gf, aggr)
    out_channel_v = size(g.W_a, 1)

    # The edge features transmitted are not useful for the node features
    nf = X[1:out_channel_v, :]

    # For the edge features, we only take what is necessary
    ef = E[out_channel_v+1:end, :]

    GeometricFlux.FeaturedGraph(graph(fg), nf=nf, ef=ef)
end

function GeometricFlux.message(l::EdgeFtLayer, x_i::AbstractVector, x_j::AbstractVector, e_ij::AbstractVector)
    x = vcat(x_i, e_ij, x_j)
    attention_logits = prelu.(l.W_a*x, l.prelu_α)
    unattended_node_features = l.W_T*x

    new_edge_features = l.W_e * (x_i + x_j) + l.W_ee * e_ij

    return vcat(attention_logits, unattended_node_features, new_edge_features)
end

function GeometricFlux.update_batch_edge(g::EdgeFtLayer, adj, E::AbstractMatrix, X::AbstractMatrix)
    n = size(adj, 1)
    edge_idx = GeometricFlux.edge_index_table(adj)
    hcat([GeometricFlux.apply_batch_message(g, i, adj[i], edge_idx, E, X) for i in 1:n]...)
end

GeometricFlux.update_batch_edge(g::EdgeFtLayer, adj, E::AbstractMatrix, X::AbstractMatrix, u) = GeometricFlux.update_batch_edge(g, adj, E, X)

function GeometricFlux.apply_batch_message(g::EdgeFtLayer, i, js, edge_idx, E::AbstractMatrix, X::AbstractMatrix)
    # TODO: when upgrading to the newer version of GeometricFlux, replace `get_feature` by `_view`
    mailbox = hcat([GeometricFlux.message(g, GeometricFlux.get_feature(X, i), GeometricFlux.get_feature(X, j), E[:, edge_idx[(i, j)]]) for j = js]...)

    # Get each part of the message separately
    out_channel_v = size(g.W_a, 1)
    attention_logits = mailbox[1:out_channel_v, :]
    unattended_node_features = mailbox[out_channel_v+1:out_channel_v*2, :]
    new_edge_features = mailbox[out_channel_v*2+1:end, :]

    # Apply the attention mechanism
    attention = Flux.softmax(attention_logits, dims=2)
    final_messages = attention .* unattended_node_features

    # Copy the bias over each neighbor
    bias_cloned = hcat([g.b_T/length(js) for j = js]...)

    # The new features are transmitted without any manipulation
    vcat(final_messages + bias_cloned, new_edge_features)
end

function Base.show(io::IO, l::EdgeFtLayer)
    in_channel_v = size(l.W_e, 2)
    out_channel_v = size(l.W_a, 1)
    in_channel_e = size(l.W_ee, 2)
    out_channel_e = size(l.W_ee, 1)
    print(io, "EdgeFtLayer(")
    print(io, "), v_dim=", in_channel_v, "=>", out_channel_v, ", e_dim=", in_channel_e, "=>", out_channel_e)
    print(io, ", PReLU(α=", l.prelu_α)
    print(io, "))")
end


text:
using ClusteredLowRankSolver
using Test

@testset "ClusteredLowRankSolver.jl" begin
    #for these tests we need to register BasesAndSamples first. Or we need to make it a submodule maybe
    @testset "examples" begin
        # these examples test nearly everything
        include("../examples/PolyOpt.jl")
        using  .PolyOpt
        _,sol,_,_ = min_f(2)
        @test sol.dual_objective ≈ -2.113 atol=1e-2

        include("../examples/Delsarte.jl")
        using .Delsarte
        _, sol, _, _ = delsarte(3, 10, 1//2)
        @test sol.dual_objective ≈ 13.158314 atol=1e-5
        _, sol, _, _ = delsarte(8, 10, 1//2)
        @test sol.dual_objective ≈ 240 atol = 1e-10

        include("../examples/SpherePacking.jl")
        using .SpherePacking
        _, sol, _, _ = cohnelkies(8, 25)
        @test sol.dual_objective ≈ BigFloat(pi)^4/384 atol=1e-4 #exact in the limit of d-> ∞, but for this d the error still is relatively large

        include("../examples/ThreePointBound.jl")
        using .ThreePointBound
        _, sol, _, _ = three_point_spherical_cap(3,6, 1//2)
        @test sol.dual_objective ≈ 12.718780 atol=1e-5
    end

    @testset "SampledMPolyElem" begin #this is mostly tested through the examples too
        using AbstractAlgebra
        R, (x,) = PolynomialRing(RealField, ["x"])
        p1 = x^2 + 2
        samples = [[i] for i=0:10]
        p2 = ClusteredLowRankSolver.SampledMPolyElem(p1,samples)
        @testset "addition" begin
            @test (p1+p2)(5) == 2*p1(5)
            @test (p2+p1)(5) == 2*p1(5)
            @test (p2+p2)(5) == 2*p1(5)
        end
        @testset "subtraction" begin
            @test (p1-p2)(4) == 0
            @test (p2-p1)(4) == p2(4)-p1(4)
            @test (p2-p2)(3) == 0
        end
        @testset "multiplication" begin
            @test (p1*p2)(5) == p1(5)^2
            @test (p2*p1)(5) == p1(5)^2
            @test (p2*p2)(5) == p1(5)^2
        end
    end

    @testset "LowRankMat(Pol)" begin
        using AbstractAlgebra
        R, (x,) = PolynomialRing(RealField, ["x"])
        A = LowRankMatPol([x],[[x^2,x^3]]) # the matrix [x^5 x^6; x^6 x^7]
        B = evaluate(A,2)
        @test [B[i,j] for i=1:2,j=1:2] == [2^5 2^6; 2^6 2^7]
        At = transpose(A)
        Bt = evaluate(At,2)
        @test [Bt[i,j] for i=1:2,j=1:2] == [2^5 2^6; 2^6 2^7]
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SetInput <: File
        name::String
        descriptor::Symbol
    end

Read .csv file with specific column names

# Arguments
- `name::String`: input file name
- `descriptor::Symbol`: file descriptor
"""
mutable struct SetInput <: File
    "input file name"
    name::String
    "file descriptor"
    descriptor::Symbol
end

function SetInput(; name, descriptor, )
    SetInput(name, descriptor, )
end

"""Get [`SetInput`](@ref) `name`."""
get_name(value::SetInput) = value.name
"""Get [`SetInput`](@ref) `descriptor`."""
get_descriptor(value::SetInput) = value.descriptor

"""Set [`SetInput`](@ref) `name`."""
set_name!(value::SetInput, val) = value.name = val
"""Set [`SetInput`](@ref) `descriptor`."""
set_descriptor!(value::SetInput, val) = value.descriptor = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using DelimitedFiles
using DifferentialEquations
using ProgressMeter
include("YfindNP.jl")
include("is_it_exit.jl")
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
Hamil(XX,YY,QQ,PP)=( (QQ-XX)^2+(PP-YY)^2 )*( (QQ+XX)^2+(PP+YY)^2 )/((PP^4+2*PP^2*(XX^2-1)+(1+XX^2)^2 )*(QQ^4+2*QQ^2*(YY^2+1)+(YY^2-1)^2 ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
function Eq_of_M(du,u,p,t)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[5]=0
return
end
condition(u,t,integrator)= u[5]==max_hit
function condition2(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
condition3b(u,t,integrator)=  u[5]==2
condition3f(u,t,integrator)=  u[5]==3

affect!(integrator) = terminate!(integrator)
function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end
# @everywhere cb2 = ContinuousCallback(condition2,affect2!, save_positions = (false,true) )
cb2 = ContinuousCallback(condition2, affect2!, rootfind=true, save_positions = (true,true))
cb2f = ContinuousCallback(condition2,  affect2!, rootfind=true, save_positions = (true,true))
cb2b = ContinuousCallback(condition2, affect2!,rootfind=true)

cb1 = DiscreteCallback(condition,affect!)
cb3b = DiscreteCallback(condition3b,affect!)
cb3f = DiscreteCallback(condition3f,affect!)

cb=CallbackSet(cb2,cb1)

cb_backward=CallbackSet(cb2b,cb3b)
cb_forward=CallbackSet(cb2f,cb3f)
max_hit=4
Energy=.25
t_end=1e5
tol_dist=1e-3
H=(2*Energy)^2

Q_start=1e-10
Q_end=1.5
P_start=0
P_end=3

n_iter_P=501
n_iter_Q=501

delta_P=(P_end-P_start)/n_iter_P
delta_Q=(Q_end-Q_start)/n_iter_Q

ArrP= P_start: delta_P: P_end
ArrQ= Q_start: delta_Q: Q_end
#up
Q= 0.4275000000715
P= 1.53
#down
Q_check=zeros(0)
P_check=zeros(0)

figure()
@showprogress 1 "Computing..." for i=1:length(ArrP)
    for j=1:length(ArrQ)
        P=ArrP[i]
        Q=ArrQ[j]
        Y=YfindNP(Q,P,H)
if ~isempty(Y) && Y>0
    u0=zeros(5)
    u0[1]=0 #X=0
    u0[2]=P #P
    u0[3]=Q #Q
    u0[4]=Y #Y
    u0[5]=0 # counter variable
    #are the initial conditions right?
    # prob_test= ODEProblem(Eq_of_M,u0,(0., 100.))
    # sol_test=solve(prob_test, RK4(),maxiters=1e10 )
    if true #sol_test.u[2][1]>0 #check if upcrossing
        #THESE GUYS HIT ONCE AND LEAVE

        prob_forward = ODEProblem(Eq_of_M,u0,(0., t_end))
        sol_forward=solve(prob_forward, Vern9(),maxiters=1e20, reltol=1e-10,abstol=1e-12,callback=cb, save_start=true ,save_end=true,save_everystep=true)
        # plot(sol_forward.t,sol_forward[1,:])
        B=hcat(sol_forward.u[end-10],sol_forward.u[end])'
        if is_it_exit(B,H,tol_dist) && sol_forward.u[end][5]==2
         #THESE GUYS HIT ONCE AND LEAVE
            prob_backward = ODEProblem(Eq_of_M,u0,(0., -t_end))
            sol_backward=solve(prob_backward, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_backward,save_start=true,save_end=true,save_everystep=true)
            uf=zeros(5)
            uf[1]=0 #X
            uf[2]=sol_backward[2,end] #P
            uf[3]=sol_backward[3,end] #Q
            uf[4]=sol_backward[4,end] #Y
            uf[5]=0
            dH=abs(H_test(uf)-H)
            prob_check= ODEProblem(Eq_of_M,uf,(0., t_end))
            sol_check=solve(prob_check, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_forward,save_start=true,save_end=true,save_everystep=true)
            if sol_backward.u[end][5]==2
                push!(Q_check, Q)
                push!(P_check, P)
                plot(sol_check.t,sol_check[1,:])
                # plot(sol_backward.t,sol_backward[1,:])
            end
        end
  end
end
end
end


text:
import Base.Iterators: flatten

function ast_sym_esc!(array_args)
    for i=1:length(array_args)
        if (typeof(array_args[i]) == Symbol)
            array_args[i] = esc(array_args[i])
        elseif (typeof(array_args[i]) == Expr)
            ast_sym_esc!(array_args[i].args)
        end
    end
end

function ast_lookahead_loop_vars(ast_node)
    if (isa(ast_node, Array))
        return union(flatten(map(ast_lookahead_loop_vars, ast_node)))
    elseif (!isa(ast_node, Expr))
        return Vector{Symbol}(undef,0)
    end

    if (ast_node.head == :kw || ast_node.head == :generator || ast_node.head == :(=))
        return [ast_node.args[1]]
    end
    return Vector{Symbol}(undef,0)
end

function ast_sym2_esc!(ast_node, skipsym=Vector{Symbol}(undef,0))
    # TODO(xhub) move thus into a seperate function that would get called first
    if (isa(ast_node, Symbol))
        ast_node = esc(ast_node)
        return ast_node
    elseif (isa(ast_node, Array))
        return map(x -> ast_sym2_esc!(x, skipsym), ast_node)
    elseif (!isa(ast_node, Expr))
        return
    end

    # is we have those instances, do not quote
    if (ast_node.head == :call || ast_node.head == :callmacro || ast_node.head == :kw || ast_node.head == :(=))
        if (ast_node.head == :kw)
            push!(skipsym, ast_node.args[1])
        end
        s=2
    elseif (ast_node.head == :generator)
        # we have faith here that the iterator has a structure like
        # ast_node.args[2] ~ i=1:/N
        push!(skipsym, ast_node.args[2].args[1])
        # now visit ast_node.args[1] and ast_node.args[2].args[2] (i=1:N)
        if (typeof(ast_node.args[1]) != Expr)
            dump(ast_node)
            error("in ast_sym2_esc")
        end
        ast_sym2_esc!(ast_node.args[1], skipsym)
        if (typeof(ast_node.args[2].args[2]) != Expr)
            dump(ast_node.args[2])
            error("in ast_sym2_esc")
        end
        ast_sym2_esc!(ast_node.args[2].args[2], skipsym)
        return ast_node
    else
        s=1
    end

    for i=s:length(ast_node.args)
        if (typeof(ast_node.args[i]) == Symbol && !(ast_node.args[i] in skipsym))
            ast_node.args[i] = esc(ast_node.args[i])
        elseif (typeof(ast_node.args[i]) == Expr)
            ast_sym2_esc!(ast_node.args[i], skipsym)
        end
    end
    return ast_node
end

function ast_args_esc!(array_args)
    for i=1:length(array_args)
        if (typeof(array_args[i]) == Expr)
            array_args[i] = esc(array_args[i])
        end
    end
end

function ast_var_bracket_esc!(expr)
    if expr.head == :kw
        # We want to skip expr.args[1] since it is a loop variable
        ast_sym2_esc!(expr.args[2])
    else
        ast_sym2_esc!(expr)
    end
end

function esc_sym!(ex, start)
    # We want to skip escaping the inner loop var
    skipsym = ast_lookahead_loop_vars(ex.args[start:end])
    for i=start:length(ex.args)
        ex.args[i] = ast_sym2_esc!(ex.args[i], skipsym)
    end
    return skipsym
end

function esc_var!(ex)
    skipsym=Vector{Symbol}(undef,0)
    if (isa(ex, Symbol))
        varname = esc(ex)
    elseif (!isa(ex, Expr))
        dump(ex)
        error("Could not parse variable declaration $(ex); please file a bug report at https://github.com/xhub/EMP.jl/issues/new")
    elseif (ex.head == :ref || ex.head == :typed_vcat)
        # This is the case there we have var[...] or var[...; cond]
        # We can have var[i=1:3] or var[1:3]
        # - case var[1:3]: ex.args[2].head is call; ex.args[2].args is
        #   Array{Any}((3,))
        #     1: Symbol :
        #     2: Int64 1
        #     3: Int64 3
        #
        # case var[i=1:3]: ex.args[2].head is kw

        varname = esc(ex.args[1])
        skipsym = esc_sym!(ex, 2)
    end
    return varname, skipsym
end

function esc_jump_id!(ex)
    if (typeof(ex) == Symbol)
        varname = esc(ex)
        ex = esc(ex)
        skipsym = Vector{Symbol}(undef,0)
    elseif (ex.head == :ref || ex.head == :typed_vcat)
        # This is the case there we have var[...]
        # We can have var[i=1:3] or var[1:3]
        # - case var[1:3]: ex.args[2].head is call; ex.args[2].args is
        #   Array{Any}((3,))
        #     1: Symbol :
        #     2: Int64 1
        #     3: Int64 3
        #
        # case var[i=1:3]: ex.args[2].head is kw

        varname, skipsym = esc_var!(ex)
        # We want to skip escaping the inner loop var
    elseif (ex.head == :call || ex.head == :comparison)
        # This is the case there we have var >= lb or ub >= var[...] >= lb
        if (ex.head == :comparison)
            posvar = 3
        else
            if (length(ex.args) == 3)
                posvar = 2
                # Here we are trying to test for the case 0 >= var or 0 >= var[...]
                if (typeof(ex.args[posvar]) <: Number)
                    posvar = 3
                end
            else
                dump(ex)
                error("Could not parse variable declaration $(ex); please file a bug report at https://github.com/xhub/EMP.jl/issues/new")
            end
        end

        varname, skipsym = esc_var!(ex.args[posvar])
    elseif (ex.head == :vect)
        # This is the case when we have a anonymous variable like
        # @variableMP(m, [1:4])
        skipsym = Vector{Symbol}(undef,0)
        varname = nothing
    elseif (ex.head == :vcat)
        # This is the case when we have a anonymous variable like
        # @variableMP(m, [i=1:4])
        skipsym = esc_sym!(ex, 1)
        varname = nothing
    else
        dump(ex)
        error("Could not parse variable declaration $(ex); please file a bug report at https://github.com/xhub/EMP.jl/issues/new")
    end

    return varname, skipsym
end

"""
Add a variable to a Mathematical Programm. See JuMP `@variable` for examples
"""
macro variableMP(args...)
    mp = esc(args[1])
    extra, kw_args, requestedcontainer = JuMP._extract_kw_args(args[2:end])

    # TODO(xhub) fix this with proper syntax
    # dummyconstr = Expr(:call, @variable, $(mp).emp.m, $(esc(args[1])))
    #then extend the arguments

    if length(extra) == 0
        v = gensym()
        varname = v
        skipsym = Vector{Symbol}(undef,0)
    else
        v = popfirst!(extra)
        varname, skipsym = esc_jump_id!(v)
    end

    if length(extra) > 0
        jump_call = :(JuMP.@variable mmp $v $(extra...))
    else
        jump_call = :(JuMP.@variable mmp $v)
    end

    if length(kw_args) > 0
        ast_sym2_esc!(kw_args, skipsym)
        append!(jump_call.args, kw_args)
    end
    push!(jump_call.args, Expr(:(=), :container, requestedcontainer))

    if (isnothing(varname))
        code = :(tmp = $jump_call; addvar!($mp, tmp); tmp)
    else
        code = :($varname = $jump_call; addvar!($mp, $varname); $varname)
    end

    return quote
        mmp = $mp.emp.model_ds
        $code
    end
end

macro constraintFromExprMP(mp, expr)
    return quote
        cref = @constraint $(esc(mp)).emp.model_ds $expr == 0
        addequ!($(esc(mp)), cref, false)
        cref
    end
end

macro NLconstraintFromExprMP(mp, expr)
    return quote
        cref = @NLconstraint $(esc(mp)).emp.model_ds $expr == 0
        addequ!($(esc(mp)), cref, true)
        cref
    end
end

# This comes from MOI
function _sense(sense)
    if sense == :Min
        expr = MOI.MIN_SENSE
    elseif sense == :Max
        expr = MOI.MAX_SENSE
    else
        expr = esc(sense)
    end
    return expr
end

# TODO(xhub) add check warning for an already defined objective equation
"""
    @NLobjectiveMP(mp, sense, expr)

Add a nonlinear objective to a mathematical programm. Sense is either `Min` or `Max`.
"""
macro NLobjectiveMP(mp, sense, expr)
    ast_sym2_esc!(expr)
    dummyconstr = Expr(:call, :(==), expr, 0)
    mmp = esc(mp)
    sense_sym = _sense(sense)
    return quote
        model = $(mmp).emp.model_ds
        cref = @NLconstraint model $dummyconstr
        addequ!($mmp, cref, true)
        $(mmp).objequ = (cref.index.value, true)
        $(mmp).sense = $sense_sym
        cref
    end
end

"""
    @objectiveMP(mp, sense, expr)

Add a linear objective to a mathematical programm. Sense is either `Min` or `Max`.
"""
macro objectiveMP(mp, sense, expr)
    ast_sym2_esc!(expr)
    dummyconstr = Expr(:call, :(==), expr, 0)
    mmp = esc(mp)
    sense_sym = _sense(sense)
    return quote
        model = $(mmp).emp.model_ds
        cref = @constraint model $dummyconstr
        addequ!($mmp, cref, false)
        $(mmp).objequ = (cref.index.value, false)
        $(mmp).sense = $sense_sym
        cref
    end
end

"""
    @constraintMP(mp, expr)

Add a linear constraint to a mathematical programm
"""
macro constraintMP(mp, expr)
    ast_args_esc!(expr.args)
    return quote
        mmp = $(esc(mp))
        model = mmp.emp.model_ds
        cref = @constraint model $expr
        addequ!(mmp, cref, false)
        cref
    end
end

"""
    @constraintMP(mp, name, expr)

Add a linear constraint (with a identifier `name`) to a mathematical programm
"""
macro constraintMP(mp, name, expr)
    equname, skipsym = esc_jump_id!(name)
    ast_sym2_esc!(expr.args, skipsym)
    return quote
        mmp = $(esc(mp))
        model = mmp.emp.model_ds
        cref = @constraint model $name $expr
        addequ!(mmp, cref, false)
        cref
    end
end

"""
    @NLconstraintMP(mp, expr)

Add a nonlinear constraint to a mathematical programm
"""
macro NLconstraintMP(mp, expr)
    ast_sym2_esc!(expr)
    return quote
        mmp = $(esc(mp))
        model = mmp.emp.model_ds
        cref = @NLconstraint model $expr
        addequ!(mmp, cref, true)
        cref
    end
end

"""
    @NLconstraintMP(mp, name, expr)

Add a nonlinear constraint (with a identifier `name`) to a mathematical programm
"""
macro NLconstraintMP(mp, name, expr)
    equname, skipsym = esc_jump_id!(name)
    ast_sym2_esc!(expr.args, skipsym)
    return quote
        mmp = $(esc(mp))
        model = mmp.emp.model_ds
        cref = @NLconstraint model $name $expr
        addequ!(mmp, cref, true)
        cref
    end
end

"""
    @vipair(mp, expr, var)

Add an affine variational inequality relationship between the variable `var` and the mapping `expr`

If there is no additional constraint of the feasible set of `var` besides lower and upper bounds,
then this is defines a Mixed Complementarity Problem
    expr ⟂ lb ≤ var ≤ ub
Or
    0 ∈ expr + N_[lb, ub] (var)
"""
#function vipair(mp::MathPrgm, expr, var::JuMP.VariableRef)
#    cref, eidx = @constraintFromExprMP(mp, expr)
#    mp.matching[var.index] = eidx
#    cref
#end
#
#function vipair(mp::MathPrgm, expr::Vector, var::Vector{JuMP.VariableRef})
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        vipair(mp, expr[i], var[i])
#    end
#end
#
#function vipair(mp::MathPrgm, expr::Vector, var::JuMP.JuMPArray)
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        vipair(mp, expr[i], var[i])
#    end
#end

macro vipair(mp, expr, var)
    if (expr.head == :vect || expr.head == :generator || expr.head == :comprehension)
        use_broadcast = true
    else
        use_broadcast = false
    end
    ast_sym2_esc!(expr)
    if use_broadcast
        dummyconstr = Expr(:call, :(.==), expr, 0)
    else
        dummyconstr = Expr(:call, :(==), expr, 0)
    end
    mmp = esc(mp)
    mvar = esc(var)
    return quote
        model = $(mmp).emp.model_ds
        cref = @constraint model $dummyconstr
        addequ!($mmp, cref, false)
        if (cref isa Array)
            for i=1:length(cref)
                $(mmp).matching[$(mvar)[i].index.value] = (cref[i].index.value, false)
            end
        else
            $(mmp).matching[$(esc(var)).index.value] = (cref.index.value, false)
        end
        cref
    end
end

#macro vipair(mp, expr, var)
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        vipair(mp, expr[i], var[i])
#    end
#end
#
#macro vipair(mp, expr, var)
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        vipair(mp, expr[i], var[i])
#    end
#end
#
"""
    @NLvipair(mp, expr, var)

Add a nonlinear variational inequality relationship between the variable `var` and the mapping `expr`

If there is no additional constraint of the feasible set of `var` besides lower and upper bounds,
then this is defines a Mixed Complementarity Problem
    expr ⟂ lb ≤ var ≤ ub
Or
    0 ∈ expr + N_[lb, ub] (var)
"""
macro NLvipair(mp, expr, var::JuMP.VariableRef)
    cref = @NLconstraintFromExprMP(mp, expr)
    mp.matching[var.index.value] = (cref.index.value, true)
    cref
end

#function NLvipair(mp::MathPrgm, expr::Vector, var::Vector{JuMP.VariableRef})
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        NLvipair(mp, expr[i], var[i])
#    end
#end


text:
using ImageComponentAnalysis
using Test, Images, JuliaDBMeta, IndexedTables

@testset "ImageComponentAnalysis.jl" begin
    include("test_images.jl")
    include("test_volumes.jl")
    include("test_ellipse_images.jl")
    include("contour_tracing.jl")
    include("measurement.jl")
    include("generic_labelling.jl")
    include("one_component_2d.jl")
    include("one_component_3d.jl")
    include("basic_topology.jl")
    include("region_ellipse.jl")
end


text:
using DICOMUtils
using Test
using TestSetExtensions


text:
## squared_exponential.jl : implementation of squared exponential covariance function

## SquaredExponential ##
struct SquaredExponential{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T

    function SquaredExponential{T}(λ::T, σ::T, p::T) where T
        λ > 0 || throw(DomainError(λ, "correlation length λ of squared exponential covariance cannot be negative or zero"))
        σ > 0 || throw(DomainError(σ, "marginal standard deviation σ of squared exponential covariance cannot be negative or zero"))
        p >= 1 || throw(DomainError(p, "in p-norm, p must be greater than or equal to 1"))
        isinf(p) && throw(DomainError(p, "in p-norm, p cannot be infinity"))

        new{T}(λ, σ, p)
    end
end

"""
    SquaredExponential(λ, [σ = 1], [p = 2])

Squared exponential (Gaussian) covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm, defined as

``C(x, y) = σ \\exp\\left(-\\left(\\displaystyle\\frac{ρ}{λ}\\right)^2\\right)``

with ``ρ = ||x - y||_p``.

# Examples
```jldoctest
julia> SquaredExponential(0.1)
Gaussian (λ=0.1, σ=1.0, p=2.0)

julia> SquaredExponential(1, σ=2.)
Gaussian (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Exponential`](@ref), [`Linear`](@ref), [`Spherical`](@ref), [`Whittle`](@ref), [`Gaussian`](@ref), [`Matern`](@ref)
"""
SquaredExponential(λ::Real; σ::Real=1.0, p::Real=2) =
    SquaredExponential{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ, σ, p)...)

"""
    Gaussian(λ, [σ = 1], [p = 2])

Gaussian (squared exponential) covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm, defined as

``C(x, y) = σ \\exp\\left(-\\left(\\displaystyle\\frac{ρ}{λ}\\right)^2\\right)``

with ``ρ = ||x - y||_p``.

# Examples
```jldoctest
julia> Gaussian(0.1)
Gaussian (λ=0.1, σ=1.0, p=2.0)

julia> Gaussian(1, σ=2.)
Gaussian (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Exponential`](@ref), [`Linear`](@ref), [`Spherical`](@ref), [`Whittle`](@ref), [`SquaredExponential`](@ref), [`Matern`](@ref)
"""
const Gaussian = SquaredExponential

# evaluate Gaussian covariance
apply(s::SquaredExponential, x::Real) = exp(-(x / s.λ)^2)

# short name
shortname(::SquaredExponential) = "Gaussian"


text:
# implementation of MLJ measure interface for LossFunctions.jl

naked(T::Type) = split(string(T), '.') |> last |> Symbol

const MARGIN_LOSSES =
    [:DWDMarginLoss, :ExpLoss, :L1HingeLoss, :L2HingeLoss, :L2MarginLoss,
     :LogitMarginLoss, :ModifiedHuberLoss, :PerceptronLoss, :SigmoidLoss,
     :SmoothedL1HingeLoss, :ZeroOneLoss]
const DISTANCE_LOSSES =
    [:HuberLoss, :L1EpsilonInsLoss, :L2EpsilonInsLoss, :LPDistLoss,
     :LogitDistLoss, :PeriodicLoss, :QuantileLoss]

# Supervised Loss -- measure traits
const LSF = LossFunctions
is_measure_type(::Type{<:SupervisedLoss})          = true
orientation(::Type{<:SupervisedLoss})              = :loss
reports_each_observation(::Type{<:SupervisedLoss}) = true
is_feature_dependent(::Type{<:SupervisedLoss})     = false

MMI.supports_weights(::Type{<:SupervisedLoss}) = true
MMI.docstring(M::Type{<:SupervisedLoss})       = name(M)
instances(M::Type{<:SupervisedLoss}) = [snakecase(string.(naked(M))), ]

## DISTANCE BASED LOSS FUNCTION

MMI.prediction_type(::Type{<:DistanceLoss}) = :deterministic
MMI.target_scitype(::Type{<:DistanceLoss}) = Union{Vec{Continuous},Vec{Count}}

function value(measure::DistanceLoss, yhat, X, y, ::Nothing,
                ::Val{false}, ::Val{true})
    return LSF.value(measure, yhat, y)
end

function value(measure::DistanceLoss, yhat, X, y, w,
                ::Val{false}, ::Val{true})
    return w .* value(measure, yhat, X, y, nothing) ./ (sum(w)/length(y))
end

## MARGIN BASED LOSS FUNCTIONS

MMI.prediction_type(::Type{<:MarginLoss}) = :probabilistic
MMI.target_scitype(::Type{<:MarginLoss})  = AbstractArray{<:Finite{2}}

# rescale [0, 1] -> [-1, 1]:
_scale(p) = 2p - 1

function value(measure::MarginLoss, yhat, X, y, ::Nothing,
                ::Val{false}, ::Val{true})
    check_pools(yhat, y)
    probs_of_observed = broadcast(pdf, yhat, y)
    return (LSF.value).(measure, _scale.(probs_of_observed), 1)
end

function value(measure::MarginLoss, yhat, X, y, w,
                ::Val{false}, ::Val{true})
    return w .* value(measure, yhat, X, y, nothing) ./ (sum(w)/length(y))
end

## KEYWORD CONSTRUCTORS

# distance:
DWDMarginLoss(; q=1.0)         = LossFunctions.DWDMarginLoss(q)
SmoothedL1HingeLoss(; γ=1.0) = LossFunctions.SmoothedL1HingeLoss(γ)

# margin:
HuberLoss(x...; d=1.0)        = LossFunctions.HuberLoss(d)
L1EpsilonInsLoss(; ϵ=1.0) = LossFunctions.L1EpsilonInsLoss(ϵ)
L2EpsilonInsLoss(; ϵ=1.0) = LossFunctions.L2EpsilonInsLoss(ϵ)
LPDistLoss(; P=2)         = LossFunctions.LPDistLoss(P)
QuantileLoss(; τ=0.7)     = LossFunctions.QuantileLoss(τ)

## ADJUSTMENTS

human_name(::Type{<:L1EpsilonInsLoss}) = "l1 ϵ-insensitive loss"
human_name(::Type{<:L2EpsilonInsLoss}) = "l2 ϵ-insensitive loss"
human_name(::Type{<:DWDMarginLoss}) = "distance weighted discrimination loss"

_signature(::Any) = ""
_signature(::Type{<:DWDMarginLoss}) = "`DWDMarginLoss(; q=1.0)`"
_signature(::Type{<:SmoothedL1HingeLoss}) = "`SmoothedL1HingeLoss(; γ=1.0)`"
_signature(::Type{<:L1EpsilonInsLoss}) = "`L1EpsilonInsLoss(; ϵ=1.0)`"
_signature(::Type{<:L2EpsilonInsLoss}) = "`L2EpsilonInsLoss(; ϵ=1.0)`"
_signature(::Type{<:LPDistLoss}) = "`LPDistLoss(; P=2)`"
_signature(::Type{<:QuantileLoss}) = "`QuantileLoss(; τ=0.7)`"

## ALIASES AND DOCSTRINGS

for M_ex in DISTANCE_LOSSES
    eval(quote
         sig = _signature($M_ex)
         isempty(sig) || (sig = "Constructor signature: "*sig)
         @create_aliases $M_ex
         @create_docs($M_ex,
            body="See above for original LossFunctions.jl documentation. ",
            footer=sig)
         end)
end

for M_ex in MARGIN_LOSSES
    eval(quote
         sig = _signature($M_ex)
         isempty(sig) || (sig = "Constructor signature: "*sig)
         @create_aliases $M_ex
         @create_docs($M_ex,
            body="See above for original LossFunctions.jl documentation. ",
            scitype=DOC_FINITE_BINARY,
            footer= sig)
         end)
end

## CALLING BEHAVIOUR

for M in vcat(subtypes(DistanceLoss), subtypes(MarginLoss))
    (m::M)(yhat, y) = MLJBase.value(m, yhat, nothing, y, nothing)
    (m::M)(yhat, y, w) = MLJBase.value(m, yhat, nothing, y, w)
end


text:
"""
    `LineSegment(r1,r2)`

# ARGUMENTS
* `r1::Array{Float64}`        : Starting point of panel
* `r2::Array{Float64}`        : Ending point of panel

# PROPERTIES
* `r1::Array{Float64}`        : Starting point of segment
* `r2::Array{Float64}`        : Ending point of segment
* `L::Float64`                : Length of segment
* `theta::Float64`            : Angle of panel relative to global X-axis (2D ONLY)
* `n_hat::Array{Float64}`     : Normal vector (2D ONLY)
* `t_hat::Array{Float64}`     : Tangent vector
"""
mutable struct LineSegment

    # -- Properties --
    R1::Array{Float64,1}
    R2::Array{Float64,1}
    L::Float64
    theta::Float64
    n_hat::Array{Float64,1}
    t_hat::Array{Float64,1}

    # -- Constructor --
    function LineSegment(r1::Array{Float64,1},r2::Array{Float64,1})

        if length(r1) == 2
            x1, y1 = r1[1], r1[2]
            x2, y2 = r2[1], r2[2]

            L = sqrt((y2 - y1)^2 + (x2 - x1)^2)
            theta = atan((y2 - y1),(x2 - x1))

            t_hat = (r2 - r1)/norm(r2 - r1)
            n_hat = [t_hat[2],-t_hat[1]]

        elseif length(r1) == 3
            x1, y1, z1 = r1[1], r1[2], r1[3]
            x2, y2, z2 = r2[1], r2[2], r2[3]

            L = sqrt((z2 - z1)^2 + (y2 - y1)^2 + (x2 - x1)^2)
            theta = 0.0

            t_hat = (r2 - r1)/norm(r2 - r1)
            n_hat = [0.0, 0.0, 0.0]

        end

        new(r1,r2,L,theta,n_hat,t_hat)
    end
end


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed 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

# Load a CSV file into the given database with the given relation name.

using RAI: Context, HTTPError, load_config, load_csv

include("parseargs.jl")

# Returns the file name without path and extension.
_sansext(fullname) = first(splitext(last(splitdir(fullname))))

function run(database, engine, relation, fullname; profile, kw...)
    isnothing(relation) && (relation = _sansext(fullname))
    data = read(fullname, String)
    cfg = load_config(; profile = profile)
    ctx = Context(cfg)
    rsp = load_csv(ctx, database, engine, relation, data; kw...)
    println(rsp)
end

function main()
    args = parseargs(
        "database", Dict(:help => "database name", :required => true),
        "engine", Dict(:help => "engine name", :required => true),
        "file", Dict(:help => "source file", :required => true),
        "--header-row", Dict(
            :help => "header row number, 0 for no header (default: 1)",
            :arg_type => Int),
        "--delim", Dict(:help => "field delimiter"),
        "--escapechar", Dict(:help => "character used to escape quotes"),
        "--quotechar", Dict(:help => "quoted field character"),
        ["--relation", "-r"], Dict(:help => "relation name (default: file name)"),
        "--profile", Dict(:help => "config profile (default: default)"))
    try
        run(args.database, args.engine, args.relation, args.file;
            delim = args.delim, header_row = args["header-row"],
            escapechar = args.escapechar, quotechar = args.quotechar,
            profile = args.profile)
    catch e
        e isa HTTPError ? show(e) : rethrow()
    end
end

main()

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
## Example from Bradlyn, Science, Eq. (S66)
# SG199
g3 = S"y+1,z,x+1"       # {3₁₁₁⁻|1,0,1}
g1 = S"x-1/2,-y+1/2,-z" # 2₁₀₀|-½,½,0
g2 = S"-x,y+1/2,-z-1/2" # {2₀₁₀|0,½,-½}

G3 = [0.0 0 1; 1 0 0; 0 1 0]
G1 = [-1.0 0 0; 0 -1 0; 0 0 1]
G2 = inv(G3)*G1*(G3)

lg = LittleGroup(199, KVec("0.5,0.5,0.5"), "P", [g1,g2,g3,])#g1, g2, g3]) #
lgir = LGIrrep(
    "P?",
    lg,
    complex.([G1,G2,G3,]),#, G2, G3]), #
    [zeros(3) for _ in 1:3],
    REAL,
    false
    )

H = kdotp(lgir)

## -----
# SG220

g3 = S"z,x,y"
g4 = S"-x+1/2,z+1,-y+1"
m = S"y+1/2,x+1/2,z+1/2"
G3 = -1.0im*[0 1 0; 0 0 -1; 1 0 0]
G4 = [0 0 1; 0 1 0; -1 0 0]*cispi(1/4)
M = [1 0 0; 0 0 1; 0 1 0]*cispi(3/4)

lg = LittleGroup(220, KVec("0.5,0.5,0.5"), "P", [g3, g4, m])
lgir = LGIrrep(
    "P?",
    lg,
    complex.([G3, G4, M]),
    [zeros(3) for _ in 1:3],
    REAL,
    false
    )
kdotp(lgir)


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_grad.jl")
include("calc_rho.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("calc_Energies.jl")

function test_main( Ns )

    const LatVecs = 16.0*diagm( ones(3) )

    pw = PWGrid( Ns, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    #
    # Generate array of distances
    #
    center = sum(LatVecs,2)/2
    dr = gen_dr( r, center )
    #
    # Setup potential
    #
    V_ionic = init_pot_harm_3d( pw, dr )
    #
    const Nstates = 4
    Focc = 2.0*ones(Nstates)

    # Initialize random wave function and orthonormalize in G-space
    srand(1234)
    psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
    psi = ortho_gram_schmidt(psi)

    Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
    rho = calc_rho( pw, Focc, psi )
    Potentials.Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
    Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )

    Energies = calc_Energies( pw, Potentials, Focc, psi, 0.0 )
    E0 = Energies.Total
    g0 = calc_grad( pw, Potentials, Focc, psi )

    println("E0 = ", E0)
    println("sum(g0) = ", sum(g0))

    # random direction
    dW = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)

    for ii = 1:9
        delta = 10.0^(-ii)
        #
        psiw = psi + delta*dW
        psiw = ortho_gram_schmidt(psiw)
        Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
        rho = calc_rho( pw, Focc, psiw )
        Potentials.Hartree = real( G_to_R( Ns, Poisson_solve(pw, rho) ) )
        Potentials.XC = excVWN( rho ) + rho .* excpVWN( rho )
        #
        dE = 2*real(trace(g0'*delta*dW))
        #
        Energies = calc_Energies( pw, Potentials, Focc, psiw, 0.0 )
        E = Energies.Total
        g = calc_grad( pw, Potentials, Focc, psiw )
        #
        println("\ndelta = ", delta)
        # actual change / expected change
        println("Actual / expected = ", (E-E0)/dE)
        # e
        println("Error estimate = ", sqrt(Ngwx)*eps()/abs(dE) )
    end

end

test_main( [30,30,30] )


text:
module TestHeatTransfer

using Modia
@usingModiaPlot
using Test

include("$(Modia.path)/models/HeatTransfer.jl")

SimpleRod = Model(
    fixedT = FixedTemperature | Map(T = 393.15u"K"),
         C = HeatCapacitor    | Map(C = 200u"J/K"),
         G = ThermalConductor | Map(G = 100u"W/K"),
    connect = :[
        (fixedT.port, G.port_a),
        (G.port_b   , C.port)
    ]
)

simpleRod = @instantiateModel(SimpleRod)
simulate!(simpleRod, Tsit5(), stopTime = 15, requiredFinalStates = [393.09468783194814])
plot(simpleRod, ("fixedT.port.T", "C.T"), figure=1)

HeatedRod = Model(
    fixedT     = FixedTemperature | Map(T = 493.15u"K"),
    fixedQflow = FixedHeatFlow,
    rod = InsulatedRod | Map(T = Var(init = fill(293.15,5)u"K")),
    connect = :[
        (fixedT.port, rod.port_a),
        (rod.port_b , fixedQflow.port)
    ]
)

heatedRod = @instantiateModel(HeatedRod)
simulate!(heatedRod, stopTime = 1e5, log=true,
          requiredFinalStates = [492.9629728925529, 492.607421697723, 492.30480548105595, 492.0850627579106, 491.9694736486912])
plot(heatedRod, [("fixedT.port.T", "rod.T"), "der(rod.T)"], figure=2)

end


text:
#=
The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

 1: 1
 3: 1,3
 6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?
=#

function num_divisors(n)
  count,i = 2,2 # accounts for 'n' and '1'
  while (i^2 <= n)
    if n%i==0
      count += 2
    end
    i += 1
  end
  count
end

function triangle_number(n)
  return convert(Int64,(1/2)*n*(n+1))
end

function calc()
  i = 2
  while num_divisors(triangle_number(i)) < 500
    i += 1
  end
  triangle_number(i)
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
file = joinpath(".","pointheatsource.vtu")
vtu = VTUFileHandler.VTUFile(file);
tmp = similar(vtu)
vtu += tmp
true

max_stars_repo_name: nagelt/VTUFileHandler

max_issues_repo_name: nagelt/VTUFileHandler

max_forks_repo_name: nagelt/VTUFileHandler


text:
# Julia wrapper for header: /usr/include/scip/scip_compr.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeCompr(scip, name, desc, priority, minnnodes, comprcopy, comprfree, comprinit, comprexit, comprinitsol, comprexitsol, comprexec, comprdata)
    ccall((:SCIPincludeCompr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_COMPRDATA}), scip, name, desc, priority, minnnodes, comprcopy, comprfree, comprinit, comprexit, comprinitsol, comprexitsol, comprexec, comprdata)
end

function SCIPincludeComprBasic(scip, compr, name, desc, priority, minnnodes, comprexec, comprdata)
    ccall((:SCIPincludeComprBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_COMPR}}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{SCIP_COMPRDATA}), scip, compr, name, desc, priority, minnnodes, comprexec, comprdata)
end

function SCIPsetComprCopy(scip, compr, comprcopy)
    ccall((:SCIPsetComprCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprcopy)
end

function SCIPsetComprFree(scip, compr, comprfree)
    ccall((:SCIPsetComprFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprfree)
end

function SCIPsetComprInit(scip, compr, comprinit)
    ccall((:SCIPsetComprInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprinit)
end

function SCIPsetComprExit(scip, compr, comprexit)
    ccall((:SCIPsetComprExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprexit)
end

function SCIPsetComprInitsol(scip, compr, comprinitsol)
    ccall((:SCIPsetComprInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprinitsol)
end

function SCIPsetComprExitsol(scip, compr, comprexitsol)
    ccall((:SCIPsetComprExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Ptr{Cvoid}), scip, compr, comprexitsol)
end

function SCIPfindCompr(scip, name)
    ccall((:SCIPfindCompr, libscip), Ptr{SCIP_COMPR}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetComprs(scip)
    ccall((:SCIPgetComprs, libscip), Ptr{Ptr{SCIP_COMPR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNCompr(scip)
    ccall((:SCIPgetNCompr, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetComprPriority(scip, compr, priority)
    ccall((:SCIPsetComprPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_COMPR}, Cint), scip, compr, priority)
end