text:
mapleaves(f::F, X...) where {F} = f(X...)
mapleaves(f::F, X::NTuple{N,T}...) where {F,N,T} = ntuple(i -> mapleaves(f, _getindex(X, i)...), Val(N))
mapleaves(f::F, X::AbstractArray...) where {F} = map((x...) -> mapleaves(f, x...), X...)

_getindex(X::NTuple{K,NTuple{N,T}}, i) where {K,N,T} = ntuple(j -> X[j][i], Val(K))

#####
##### Record Type for keeping track of what we're measuring
#####

"""
    Record{name}(data::T) where {T <: Union{Vector, NTuple}}

Create a named `Record` around `data`.
The elements of `data` may themselves be `Record`s, resulting in a hierarchical data structure.
"""
struct Record{name,T <: Union{Vector, NTuple}}
    data::T
end
Record{name}(data::T) where {name,T} = Record{name,T}(data)

name(R::Record{N}) where {N} = N

# Check if there are records below this record.
hassubrecord(R::Record) = eltype(R.data) <: Record

Base.getindex(R::Record, i...) = R.data[i...]
Base.length(R::Record) = length(R.data)

Base.iterate(R::Record) = iterate(R.data)
Base.iterate(R::Record, s) = iterate(R.data, s)

denest(x) = x
denest(x::Tuple{Tuple}) = denest(first(x))

# Recurse down the record stack, but apply the name to whatever is returned.
"""
    mapleaves(f, record::Record) -> Record

Apply `f` to each leaf element of `record`.
This will recursively descend through hierarchies of `Records` and only apply `f` to scalars.

The returned result will have the same hierarchical structure as `record`
"""
function mapleaves(f, R::Record{name}...) where {name}
    return Record{name}(denest(mapleaves(f, _getdata(R...)...)))
end

_getdata(x::Record, y::Record...) = (x.data, _getdata(y...)...)
_getdata() = ()

# Applying a difference to two subsequent records
Base.:-(a::R, b::R) where {R <: Record} = mapleaves(-, a, b)
Base.:+(a::R, b::R) where {R <: Record} = mapleaves(+, a, b)
mapleaves(f) = (x...) -> mapleaves(f, x...)

# Aggregate across records
"""
    aggregate(record::Record)
    aggregate(f, record::Record)

Reduce over all the leaf (terminal) elements of `record`, applying `f` as the reduction function.
If `f` is not supplied, it will defult to `(x,y) -> x .+ y`.
"""
aggregate(f, X::Record) = aggregate(f, X.data)
aggregate(f, X::Vector) = reduce(f, aggregate.(f, X))
aggregate(f, X::NTuple{N,<:Record}) where {N} = reduce(f, aggregate.(f, X))
aggregate(f, X::NTuple{N,<:NTuple}) where {N} = reduce(f, aggregate.(f, X))
aggregate(f, X) = X

# Default
aggregate(X) = aggregate((x,y) -> x .+ y, X)

#####
##### Pretty Printing
#####

titleize(x) = titlecase(replace(String(x), "_" => " "))
function Base.show(io::IO, R::Record, pre = "")
    println(io, pre, titleize(name(R)), " Record with $(length(R)) entries:")
    post = pre * "   "
    if hassubrecord(R)
        show(io, first(R.data), post)
    else
        print(io, post, eltype(R.data))
    end

    # Only show contents from the top level
    if isempty(pre)
        println(io)
        println(io)
        showcontents(io, R)
    end
end

function showcontents(io::IO, R::Record, pre = "")
    post = pre * "   "
    print(io, pre, titleize(name(R)), ": ")
    if hassubrecord(R)
        print(io, "\n")
        showcontents.(Ref(io), R.data, post)
    else
        println(io, R.data)
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Gives a reinterpreted view (of element type T) of the underlying array (of element type S).
If the size of `T` differs from the size of `S`, the array will be compressed/expanded in
the first dimension.
"""
struct ReinterpretArray{T,N,S,A<:AbstractArray{S, N}} <: AbstractArray{T, N}
    parent::A
    readable::Bool
    writable::Bool
    global reinterpret
    function reinterpret(::Type{T}, a::A) where {T,N,S,A<:AbstractArray{S, N}}
        function throwbits(::Type{S}, ::Type{T}, ::Type{U}) where {S,T,U}
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret `$(S)` `$(T)`, type `$(U)` is not a bits type"))
        end
        function throwsize0(::Type{S}, ::Type{T})
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a different size"))
        end
        function thrownonint(::Type{S}, ::Type{T}, dim)
            @_noinline_meta
            throw(ArgumentError("""
                cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
                The resulting array would have non-integral first dimension.
            """))
        end
        function throwaxes1(::Type{S}, ::Type{T}, ax1)
            @_noinline_meta
            throw(ArgumentError("cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first."))
        end
        isbitstype(T) || throwbits(S, T, T)
        isbitstype(S) || throwbits(S, T, S)
        (N != 0 || sizeof(T) == sizeof(S)) || throwsize0(S, T)
        ax1 = axes(a)[1]
        if N != 0 && sizeof(S) != sizeof(T)
            dim = length(ax1)
            rem(dim*sizeof(S),sizeof(T)) == 0 || thrownonint(S, T, dim)
            first(ax1) == 1 || throwaxes1(S, T, ax1)
        end
        readable = array_subpadding(T, S)
        writable = array_subpadding(S, T)
        new{T, N, S, A}(a, readable, writable)
    end
end

function check_readable(a::ReinterpretArray{T, N, S} where N) where {T,S}
    # See comment in check_writable
    if !a.readable && !array_subpadding(T, S)
        throw(PaddingError(T, S))
    end
end

function check_writable(a::ReinterpretArray{T, N, S} where N) where {T,S}
    # `array_subpadding` is relatively expensive (compared to a simple arrayref),
    # so it is cached in the array. However, it is computable at compile time if,
    # inference has the types available. By using this form of the check, we can
    # get the best of both worlds for the success case. If the types were not
    # available to inference, we simply need to check the field (relatively cheap)
    # and if they were we should be able to fold this check away entirely.
    if !a.writable && !array_subpadding(S, T)
        throw(PaddingError(T, S))
    end
end

IndexStyle(a::ReinterpretArray) = IndexStyle(a.parent)

parent(a::ReinterpretArray) = a.parent
dataids(a::ReinterpretArray) = dataids(a.parent)

function size(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
    psize = size(a.parent)
    size1 = div(psize[1]*sizeof(S), sizeof(T))
    tuple(size1, tail(psize)...)
end

function axes(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
    paxs = axes(a.parent)
    f, l = first(paxs[1]), length(paxs[1])
    size1 = div(l*sizeof(S), sizeof(T))
    tuple(oftype(paxs[1], f:f+size1-1), tail(paxs)...)
end

elsize(::Type{<:ReinterpretArray{T}}) where {T} = sizeof(T)
unsafe_convert(::Type{Ptr{T}}, a::ReinterpretArray{T,N,S} where N) where {T,S} = Ptr{T}(unsafe_convert(Ptr{S},a.parent))

@inline @propagate_inbounds getindex(a::ReinterpretArray{T,0}) where {T} = reinterpret(T, a.parent[])
@inline @propagate_inbounds getindex(a::ReinterpretArray) = a[1]

@inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, inds::Vararg{Int, N}) where {T,N,S}
    check_readable(a)
    _getindex_ra(a, inds[1], tail(inds))
end

@inline @propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, i::Int) where {T,N,S}
    check_readable(a)
    if isa(IndexStyle(a), IndexLinear)
        return _getindex_ra(a, i, ())
    end
    # Convert to full indices here, to avoid needing multiple conversions in
    # the loop in _getindex_ra
    inds = _to_subscript_indices(a, i)
    _getindex_ra(a, inds[1], tail(inds))
end

@inline _memcpy!(dst, src, n) = ccall(:memcpy, Cvoid, (Ptr{UInt8}, Ptr{UInt8}, Csize_t), dst, src, n)

@inline @propagate_inbounds function _getindex_ra(a::ReinterpretArray{T,N,S}, i1::Int, tailinds::TT) where {T,N,S,TT}
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return reinterpret(T, a.parent[i1, tailinds...])
    else
        @boundscheck checkbounds(a, i1, tailinds...)
        ind_start, sidx = divrem((i1-1)*sizeof(T), sizeof(S))
        t = Ref{T}()
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            i = 1
            nbytes_copied = 0
            # This is a bit complicated to deal with partial elements
            # at both the start and the end. LLVM will fold as appropriate,
            # once it knows the data layout
            while nbytes_copied < sizeof(T)
                s[] = a.parent[ind_start + i, tailinds...]
                nb = min(sizeof(S) - sidx, sizeof(T)-nbytes_copied)
                _memcpy!(tptr + nbytes_copied, sptr + sidx, nb)
                nbytes_copied += nb
                sidx = 0
                i += 1
            end
        end
        return t[]
    end
end

@inline @propagate_inbounds setindex!(a::ReinterpretArray{T,0,S} where T, v) where {S} = (a.parent[] = reinterpret(S, v))
@inline @propagate_inbounds setindex!(a::ReinterpretArray, v) = (a[1] = v)

@inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, inds::Vararg{Int, N}) where {T,N,S}
    check_writable(a)
    _setindex_ra!(a, v, inds[1], tail(inds))
end

@inline @propagate_inbounds function setindex!(a::ReinterpretArray{T,N,S}, v, i::Int) where {T,N,S}
    check_writable(a)
    if isa(IndexStyle(a), IndexLinear)
        return _setindex_ra!(a, v, i, ())
    end
    inds = _to_subscript_indices(a, i)
    _setindex_ra!(a, v, inds[1], tail(inds))
end

@inline @propagate_inbounds function _setindex_ra!(a::ReinterpretArray{T,N,S}, v, i1::Int, tailinds::TT) where {T,N,S,TT}
    v = convert(T, v)::T
    # Make sure to match the scalar reinterpret if that is applicable
    if sizeof(T) == sizeof(S) && (fieldcount(T) + fieldcount(S)) == 0
        return setindex!(a.parent, reinterpret(S, v), i1, tailinds...)
    else
        @boundscheck checkbounds(a, i1, tailinds...)
        ind_start, sidx = divrem((i1-1)*sizeof(T), sizeof(S))
        t = Ref{T}(v)
        s = Ref{S}()
        GC.@preserve t s begin
            tptr = Ptr{UInt8}(unsafe_convert(Ref{T}, t))
            sptr = Ptr{UInt8}(unsafe_convert(Ref{S}, s))
            nbytes_copied = 0
            i = 1
            # Deal with any partial elements at the start. We'll have to copy in the
            # element from the original array and overwrite the relevant parts
            if sidx != 0
                s[] = a.parent[ind_start + i, tailinds...]
                nb = min(sizeof(S) - sidx, sizeof(T))
                _memcpy!(sptr + sidx, tptr, nb)
                nbytes_copied += nb
                a.parent[ind_start + i, tailinds...] = s[]
                i += 1
                sidx = 0
            end
            # Deal with the main body of elements
            while nbytes_copied < sizeof(T) && (sizeof(T) - nbytes_copied) > sizeof(S)
                nb = min(sizeof(S), sizeof(T) - nbytes_copied)
                _memcpy!(sptr, tptr + nbytes_copied, nb)
                nbytes_copied += nb
                a.parent[ind_start + i, tailinds...] = s[]
                i += 1
            end
            # Deal with trailing partial elements
            if nbytes_copied < sizeof(T)
                s[] = a.parent[ind_start + i, tailinds...]
                nb = min(sizeof(S), sizeof(T) - nbytes_copied)
                _memcpy!(sptr, tptr + nbytes_copied, nb)
                a.parent[ind_start + i, tailinds...] = s[]
            end
        end
    end
    return a
end

# Padding
struct Padding
    offset::Int
    size::Int
end
function intersect(p1::Padding, p2::Padding)
    start = max(p1.offset, p2.offset)
    stop = min(p1.offset + p1.size, p2.offset + p2.size)
    Padding(start, max(0, stop-start))
end

struct PaddingError
    S::Type
    T::Type
end

function showerror(io::IO, p::PaddingError)
    print(io, "Padding of type $(p.S) is not compatible with type $(p.T).")
end

"""
    CyclePadding(padding, total_size)

Cylces an iterator of `Padding` structs, restarting the padding at `total_size`.
E.g. if `padding` is all the padding in a struct and `total_size` is the total
aligned size of that array, `CyclePadding` will correspond to the padding in an
infinite vector of such structs.
"""
struct CyclePadding{P}
    padding::P
    total_size::Int
end
eltype(::Type{<:CyclePadding}) = Padding
IteratorSize(::Type{<:CyclePadding}) = IsInfinite()
isempty(cp::CyclePadding) = isempty(cp.padding)
function iterate(cp::CyclePadding)
    y = iterate(cp.padding)
    y === nothing && return nothing
    y[1], (0, y[2])
end
function iterate(cp::CyclePadding, state::Tuple)
    y = iterate(cp.padding, tail(state)...)
    y === nothing && return iterate(cp, (state[1]+cp.total_size,))
    Padding(y[1].offset+state[1], y[1].size), (state[1], tail(y)...)
end

"""
    Compute the location of padding in a type.
"""
function padding(T)
    padding = Padding[]
    last_end::Int = 0
    for i = 1:fieldcount(T)
        offset = fieldoffset(T, i)
        fT = fieldtype(T, i)
        if offset != last_end
            push!(padding, Padding(offset, offset-last_end))
        end
        last_end = offset + sizeof(fT)
    end
    padding
end

function CyclePadding(T::DataType)
    a, s = datatype_alignment(T), sizeof(T)
    as = s + (a - (s % a)) % a
    pad = padding(T)
    s != as && push!(pad, Padding(s, as - s))
    CyclePadding(pad, as)
end

using .Iterators: Stateful
@pure function array_subpadding(S, T)
    checked_size = 0
    lcm_size = lcm(sizeof(S), sizeof(T))
    s, t = Stateful{<:Any, Any}(CyclePadding(S)),
           Stateful{<:Any, Any}(CyclePadding(T))
    isempty(t) && return true
    isempty(s) && return false
    while checked_size < lcm_size
        # Take padding in T
        pad = popfirst!(t)
        # See if there's corresponding padding in S
        while true
            ps = peek(s)
            ps.offset > pad.offset && return false
            intersect(ps, pad) == pad && break
            popfirst!(s)
        end
        checked_size = pad.offset + pad.size
    end
    return true
end

ext: jl

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
struct PAPIError <: Exception
	msg::String
end

function PAPIError(code::Cint)
	str = ccall((:PAPI_strerror, libpapi), Cstring, (Cint,), code)
	msg = if str == C_NULL
		"Unknown error code $code"
	else
		Base.unsafe_string(str)
	end
	PAPIError(msg)
end

macro papichk(expr)
    :((errcode = $(esc(expr))) == PAPI_OK || throw(PAPIError(errcode)))
end


text:
module CPUTime

export
    CPUtime_us,
    CPUtic,
    CPUtoq,
    CPUtoc,
    @CPUtime,
    @CPUelapsed

function CPUtime_us()
    rusage = Libc.malloc(4*sizeof(Clong) + 14*sizeof(UInt64))  # sizeof(uv_rusage_t); this is different from sizeof(rusage)
    ccall(:uv_getrusage, Cint, (Ptr{Nothing},), rusage)
    utime = UInt64(1000000)*unsafe_load(convert(Ptr{Clong}, rusage + 0*sizeof(Clong))) +    # user CPU time
                            unsafe_load(convert(Ptr{Clong}, rusage + 1*sizeof(Clong)))
    stime = UInt64(1000000)*unsafe_load(convert(Ptr{Clong}, rusage + 2*sizeof(Clong))) +    # system CPU time
                            unsafe_load(convert(Ptr{Clong}, rusage + 3*sizeof(Clong)))
    ttime = utime + stime  # total CPU time
    Libc.free(rusage)
    return ttime
end

function CPUtic()
    t0 = CPUtime_us()
    task_local_storage(:CPUTIMERS, (t0, get(task_local_storage(), :CPUTIMERS, ())))
    return t0
end

function CPUtoq()
    t1 = CPUtime_us()
    timers = get(task_local_storage(), :CPUTIMERS, ())
    if timers === ()
        error("CPUtoc() without CPUtic()")
    end
    t0 = timers[1]::UInt64
    task_local_storage(:CPUTIMERS, timers[2])
    (t1-t0)/1e6
end

function CPUtoc()
    t = CPUtoq()
    println("elapsed CPU time: ", t, " seconds")
    return t
end

# print elapsed CPU time, return expression value
macro CPUtime(ex)
    quote
        local t0 = CPUtime_us()
        local val = $(esc(ex))
        local t1 = CPUtime_us()
        println("elapsed CPU time: ", (t1-t0)/1e6, " seconds")
        val
    end
end

# print nothing, return elapsed CPU time
macro CPUelapsed(ex)
    quote
        local t0 = CPUtime_us()
        local val = $(esc(ex))
        (CPUtime_us()-t0)/1e6
    end
end

end # module


text:
# Copyright 2020 Paul Melis (paul.melis@surf.nl)
#
# 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.

using Printf
include("halfedge.jl")

# Some utility functions to handle Julia's 1-based indexing
function get_vertex(V, idx)
    s = 1+3*(idx-1)
    return V[s:s+2]
end

function set_vertex(V, idx, v)
    s = 1+3*(idx-1)
    V[s:s+2] = v
end

function add_edge_face(edge_faces, A, B, left_face_index)
    #println("add_edge_face ", A, " ", B, " ", left_face_index)
    edge_key = (A,B)
    @assert !haskey(edge_faces, edge_key)
    edge_faces[edge_key] = left_face_index
end

const Edge = Tuple{Int32,Int32}
const Face = Int32

function subdivide(vertices::Array, loop_start::Array, loop_total::Array, loops::Array)
    #println(vertices)
    #println(loop_start)
    #println(loop_total)
    #println(loops)

    t0 = time()

    num_vertices, num_faces, num_edges, vertex_start_edges, face_start_edges, polygon_edges = build(vertices, loop_start, loop_total, loops)

    t1 = time()
    @printf("(Julia) Building half edges done in %.3fms\n", 1000*(t1-t0))

    println("(Julia) Input: $(num_vertices) vertices, $(num_faces) polygons, $(num_edges) polygon edges")

    # One new vertex for each input face, one new vertex for each edge
    output_num_vertices = num_vertices + num_faces + num_edges
    # 1 .. NV           Original input vertices (initially, are overwritten later on)
    # NV+1 .. NV+NF     New face points
    # NV+NF+1 .. end    New edge points
    output_vertices = Array{Float32}(undef, 3*output_num_vertices)
    # Copy original input vertex positions, to be modified later on
    output_vertices[1:3*num_vertices] = vertices

    # Output is always all quads, as each input face is split into n quads,
    # where n is the number of vertices in the face
    output_num_quads = sum(loop_total)
    output_loop_start = [1:4:output_num_quads*4;]
    output_loop_total = fill(UInt32(4), output_num_quads)
    output_loops = Array{UInt32}(undef, 4*output_num_quads)

    println("(Julia) Output: $(output_num_vertices) vertices, $(output_num_quads) quads")

    function face_point_index(fi)
        return num_vertices + fi
    end

    function edge_point_index(ei)
        return num_vertices + num_faces + ei
    end

    #
    # Subdivide
    #

    # Add new face points: average of existing original face vertices
    for fi = 1:num_faces

        sum = Vector{Float32}([0, 0, 0])
        n = 0

        he = start = face_start_edges[fi]
        while true
            sum += get_vertex(vertices, he.source)
            n += 1

            he = he.next
            if he == start break end
        end

        set_vertex(output_vertices, face_point_index(fi), sum / n)

    end

    #println("FV ", face_vertices)
    #println("OV ", output_vertices)

    # Add new edge points: average of edge endpoints and neighbouring face points

    for ei = 1:num_edges

        he = polygon_edges[ei]

        edge_point = (
            get_vertex(output_vertices, face_point_index(he.face))
            +
            get_vertex(output_vertices, he.source)
            +
            get_vertex(output_vertices, he.target)
        )
        n = 3

        if he.sibling != nothing
            he2 = he.sibling
            edge_point += get_vertex(output_vertices, face_point_index(he2.face))
            n += 1
        end

        set_vertex(output_vertices, edge_point_index(ei), edge_point / n)
    end

    # Move original input vertices to new positions

    for vi = 1:num_vertices

        if !haskey(vertex_start_edges, vi)
            # Skip unconnected vertices
            continue
        end

        P = get_vertex(output_vertices, vi)
        F_sum = Vector{Float32}([0, 0, 0])
        R_sum = Vector{Float32}([0, 0, 0])
        n = 0

        he = start = vertex_start_edges[vi]
        while true
            F_sum += get_vertex(output_vertices, face_point_index(he.face))
            # XXX could take out P here as use it once in R= below
            R_sum += 0.5*(P + get_vertex(output_vertices, he.target))
            n += 1

            # XXX need to recheck this works correctly
            if he.sibling == nothing break end

            he = he.sibling.next
            if he == start break end
        end

        F = F_sum / n
        R = R_sum / n

        set_vertex(output_vertices, vi, (F + 2*R + (n-3)*P) / n)
    end

    # Create new face loops for the quads for each subdivided original face

    offs = 1
    ofi = 1

    for fi = 1:num_faces

        he = start = face_start_edges[fi]
        while true
            output_loop_start[ofi] = offs
            output_loops[offs] = edge_point_index(he.index)
            output_loops[offs+1] = he.target
            output_loops[offs+2] = edge_point_index(he.next.index)
            output_loops[offs+3] = face_point_index(he.face)
            offs += 4
            ofi += 1

            he = he.next
            if he == start break end
        end
    end

    t2 = time()
    @printf("(Julia) Subdivision done in %.3fms\n", 1000*(t2-t1))

    return output_vertices, output_loop_start, output_loop_total, output_loops
end

max_stars_repo_name: KristofferC/blender-julia-test

max_issues_repo_name: KristofferC/blender-julia-test

max_forks_repo_name: KristofferC/blender-julia-test


text:
#### Constants
const MODEL_FILENAME = "perceptron_model.jld" # jld filename for storing the model
const MODEL_ID       = "perceptron_model"     # if od the model in the filesystem jld data

using Statistics

#### An abstract perceptron model
abstract type PerceptronModel{T} end

#### Linear Perceptron type
mutable struct LinearPerceptron{T<:AbstractFloat} <: PerceptronModel{T}
   α::T
   Θ::Vector{T}
   shuffle_epoch::Bool
   random_state::Integer
   max_epochs::Integer
   last_epoch::Integer
   history::Vector{Integer}
   mx::Matrix{T}          # mean stat after for z-scoring input data (X)
   sx::Matrix{T}          # standard deviation stat after for z-scoring target data (X)
   centralize::Bool
   nfeatures::Integer
end

function LinearPerceptron(X::AbstractArray{T},
                          alpha,
                          shuffle_epoch,
                          random_state,
                          max_epochs,
                          centralize) where T<:AbstractFloat

   return LinearPerceptron(alpha, # I will refactor to a constructor. Cleaner
                           Vector{T}(undef,1),
                           shuffle_epoch,
                           random_state,
                           max_epochs,
                           0,
                           Vector{Integer}(undef,1),
                           mean(X,dims=1),
                           std(X,dims=1),
                           centralize,
                           size(X,2))

end
####################################################################################

#### Linear Perceptron type
mutable struct VotedPerceptron{T<:AbstractFloat} <: PerceptronModel{T}
   α::T
   Θ#::Dict{Integer,Vector{T}}
   c#::Dict{Integer,Integer}
   k::Integer
   shuffle_epoch::Bool
   random_state::Integer
   max_epochs::Integer
   last_epoch::Integer
   history::Vector{Integer}
   mx::Matrix{T}          # mean stat after for z-scoring input data (X)
   sx::Matrix{T}          # standard deviation stat after for z-scoring target data (X)
   centralize::Bool
   nfeatures::Integer
end

function VotedPerceptron(X::AbstractArray{T},
                          alpha,
                          shuffle_epoch,
                          random_state,
                          max_epochs,
                          centralize) where T<:AbstractFloat

   return VotedPerceptron(alpha, # I will refactor to a constructor. Cleaner
                           nothing,
                           nothing,
                           0,
                           shuffle_epoch,
                           random_state,
                           max_epochs,
                           0,
                           Vector{Integer}(undef,1),
                           mean(X,dims=1),
                           std(X,dims=1),
                           centralize,
                           size(X,2))

end

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

#### Kernel Perceptron type
mutable struct KernelPerceptron{T<:AbstractFloat} <: PerceptronModel{T}
   λ::Vector{T} # lagrange vector
   max_epochs::Integer
   last_epoch::Integer
   history::Vector{Integer}
   mx::Matrix{T}          # mean stat after for z-scoring input data (X)
   sx::Matrix{T}          # standard deviation stat after for z-scoring target data (X)
   centralize::Bool
   nfeatures::Integer
   kernel::String
   width::T
   sv_x::Vector{T}
   sv_y::Vector{T}
end

function KernelPerceptron(X::AbstractArray{T},
                          max_epochs,
                          centralize,
                          kernel,
                          width) where T<:AbstractFloat

   return KernelPerceptron(zeros(T,size(X,1)),
                           max_epochs,
                           0,
                           Vector{Integer}(undef,1),
                           mean(X,dims=1),
                           std(X,dims=1),
                           centralize,
                           size(X,2),
                           kernel,
                           width,
                           Vector{T}(undef,1),
                           Vector{T}(undef,1))

end

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

#### Averaged Perceptron type
mutable struct AveragedPerceptron{T<:AbstractFloat} <: PerceptronModel{T}
   α::T
   Θ
   shuffle_epoch::Bool
   random_state::Integer
   max_epochs::Integer
   last_epoch::Integer
   history::Vector{Integer}
   mx::Matrix{T}          # mean stat after for z-scoring input data (X)
   sx::Matrix{T}          # standard deviation stat after for z-scoring target data (X)
   centralize::Bool
   nfeatures::Integer
end

function AveragedPerceptron(X::AbstractArray{T},
                          alpha,
                          shuffle_epoch,
                          random_state,
                          max_epochs,
                          centralize) where T<:AbstractFloat

   return AveragedPerceptron(alpha, # I will refactor to a constructor. Cleaner
                           0,
                           shuffle_epoch,
                           random_state,
                           max_epochs,
                           0,
                           Vector{Integer}(undef,1),
                           mean(X,dims=1),
                           std(X,dims=1),
                           centralize,
                           size(X,2))

end

## choosing types
######################################################################################################
function Model(X,#::AbstractArray{T},
               alpha,
               shuffle_epoch,
               random_state,
               max_epochs,
               centralize,
               kernel,
               width,
               mode) where T<:AbstractFloat

      println("size",size(X))
      if mode == "linear"
         return LinearPerceptron(X,
                                 alpha,
                                 shuffle_epoch,
                                 random_state,
                                 max_epochs,
                                 centralize)
      elseif mode == "kernel"
         return KernelPerceptron(X,
                                 max_epochs,
                                 centralize,
                                 kernel,
                                 width)
      elseif mode == "voted"
      return VotedPerceptron(X,
                           alpha,
                           shuffle_epoch,
                           random_state,
                           max_epochs,
                           centralize)
      elseif mode == "averaged"
      return AveragedPerceptron(X,
                           alpha,
                           shuffle_epoch,
                           random_state,
                           max_epochs,
                           centralize)
      else
         error("Invalid perceptron mode name: $(mode). \n Cadidates are: linear, kernel, voted or averaged")
      end
end

######################################################################################################
## Load and Store models (good for production)
function load(; filename::AbstractString = MODEL_FILENAME, modelname::AbstractString = MODEL_ID)
    local M
    jldopen(filename, "r") do file
        M = read(file, modelname)
    end
    M
end

function save(M::PerceptronModel; filename::AbstractString = MODEL_FILENAME, modelname::AbstractString = MODEL_ID)
    jldopen(filename, "w") do file
        write(file, modelname, M)
    end
end


text:
module TestExamples2DAPE

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Acoustic Perturbation" begin
  @testset "elixir_ape_convergence_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_ape_convergence_test.jl"),
      l2   = [0.0018873176537420558, 0.00228615992541537, 0.003125083031985467],
      linf = [0.007481011083990019, 0.009788636350308355, 0.02522255492947867])
  end

  @testset "elixir_ape_gauss.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_ape_gauss.jl"),
      l2 = [0.08005276517890283, 0.08005276517890268, 0.4187202920734123],
      linf = [0.17261097190220992, 0.17261097190220973, 1.13601894068238])
  end

  @testset "elixir_ape_gaussian_source.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_ape_gaussian_source.jl"),
      l2 = [0.004296394903650806, 0.004241280404758938, 0.006269684906035964],
      linf = [0.03970270697049378, 0.04151096349298151, 0.0640019829058819])
  end
end

end # module


text:
function Base.show(io::IO, c::Cursor)::Void
    @printf(io, "Cursor over a %s", typeof(source(c)))
end


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, @get_numbertype, NUMBERTYPE_NONE, SUPPORTED_PACKAGES, PKG_CUDA
import ParallelStencil.ParallelKernel: @require
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. allocator macros (with default numbertype)" begin
            @require !@is_initialized()
            @init_parallel_kernel($package, Float16)
            @require @is_initialized()
            @testset "mapping to package" begin
                @test @zeros(2,3)         == parentmodule($package).zeros(Float16,2,3)
                @test @zeros(Float32,2,3) == parentmodule($package).zeros(Float32,2,3)
                @test @ones(2,3)          == parentmodule($package).ones(Float16,2,3)
                @test @ones(Float32,2,3)  == parentmodule($package).ones(Float32,2,3)
                @static if $package == $PKG_CUDA
                    @test typeof(@rand(2,3))         == typeof(CUDA.CuArray(rand(Float16,2,3)))
                    @test typeof(@rand(Float64,2,3)) == typeof(CUDA.CuArray(rand(Float64,2,3)))
                else
                    @test typeof(@rand(2,3))         == typeof(parentmodule($package).rand(Float16,2,3))
                    @test typeof(@rand(Float64,2,3)) == typeof(parentmodule($package).rand(Float64,2,3))
                end
            end;
            @reset_parallel_kernel()
        end;
        @testset "2. allocator macros (no default numbertype)" begin
            @require !@is_initialized()
            @init_parallel_kernel(package = $package)
            @require @is_initialized()
            @require @get_numbertype() == NUMBERTYPE_NONE
            @testset "mapping to package" begin
                @test @zeros(Float32,2,3) == parentmodule($package).zeros(Float32,2,3)
                @test @ones(Float32,2,3)  == parentmodule($package).ones(Float32,2,3)
                @static if $package == $PKG_CUDA
                    @test typeof(@rand(Float64,2,3)) == typeof(CUDA.CuArray(rand(Float64,2,3)))
                else
                    @test typeof(@rand(Float64,2,3)) == typeof(parentmodule($package).rand(Float64,2,3))
                end
            end;
            @reset_parallel_kernel()
        end;
        @testset "3. Exceptions (no default numbertype)" begin
            @require !@is_initialized()
            @init_parallel_kernel(package = $package)
            @require @is_initialized()
            @require @get_numbertype() == $NUMBERTYPE_NONE
            @testset "no numbertype" begin
                @test_throws ArgumentError @zeros(2,3)
                @test_throws ArgumentError @ones(2,3)
                @test_throws ArgumentError @rand(2,3)
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;


text:
# Unvalidated code for attempting a stochastic integration of the differential
# equations for the TES. Galen's project.

function stochasticdy!(dy,y,dt,bt::BiasedTES)
   p=bt.p
   T, I = y # tes tempertature, tes current
   r = R(I,T,p) #res resistance
   # generate random variables for use in stocastic integration
   bath_thermal_noise, tes_thermal_noise, rl_johnson_noise, tes_johnson_noise = randn(4)

   johnsontes = tes_johnson_noise*2*sqrt(kb*T*r/dt)
   johnsonrl = rl_johnson_noise*2*sqrt(kb*p.Tbath*p.Rl/dt)

   # k*T^n*dt/kb*T is the average number of phonons going from the TES to the bath
   # assuming all phonons of energy kb*t
   # sqrt(k*T^n*dt/kb*T) is the std deviation of the number of phonons
   # the same reasoning applies in the other direction
   dqdt = -p.k*(T^p.n-p.Tbath^p.n)+I^2*r +
     tes_thermal_noise*sqrt(kb*p.k/dt)*T^(p.n/2 + 0.5) + # variation in number of phonon TES->bath
     bath_thermal_noise*sqrt(kb*p.k/dt)*p.Tbath^(p.n/2+0.5) + # variation in number of phonon bath->TES
     -I*johnsontes # work on the noise voltage source for the TES add/remove energy from the TES

   # the resistors are modeled as noiseless resistors in parallel with
   # voltage sources having white noise
   # the load resistor is assumed to be at Tbath
   ldidt = bt.V-I*(p.Rl+r) + johnsontes + johnsonrl

   if isnan(dqdt)
      @show dy, y, dt, bt
      error()
   end

   dy[1] = dt*dqdt/p.C
   dy[2] = dt*ldidt/p.L
end

"stochastic(nsample::Int, dt::Float64, bt::BiasedTES, E::Number, npresample::Int=0, nrandomizingsample::Int=6000)
Integrate a pulse record with with stocastic noise `E` eV with `nsample` total samples, `npresample` presamples, and use
`nrandomizingsample` samples thrown away to ensure randomized starting `T0` and `I0`. Use `E=0` for noise records."
function stochastic(nsample::Int, dt::Float64, bt, E::Number, npresample::Int=0, nrandomizingsample::Int=6000)
   p = bt.p
   T = Array(Float64, nsample)
   I = Array(Float64, nsample)
   y = [bt.T0, bt.I0]
   dy = Array(Float64,2)
   # manually implement euler's method for integration
   for i = 1:nrandomizingsample # integrate for a while to get well randomized T0 and I0
       stochasticdy!(dy,y,dt,bt)
       y+=dy
   end
   for i = 1:nsample
      if i == npresample+1
         y[1]+=E*J_per_eV/p.C
      end # increase temperature to start pulse
      stochasticdy!(dy,y,dt,bt)
      y+=dy
      T[i] = y[1]
      I[i] = y[2]
   end
   TESRecord(T,I,R(I,T,bt.p),dt)
end


text:
# Let's cover some basics of Makie
#
using CairoMakie
using AlgebraOfGraphics
kwargs = AlgebraOfGraphics.aog_theme()
set_theme!(;resolution=(500, 500),kwargs...)
fig, ax, pltobj = scatterlines(1:10);


text:
using Pkg
Pkg.update("TuringWeb")

using Documenter, TuringWeb
using AdvancedHMC

# cp(joinpath(@__DIR__, "../README.md"), joinpath(@__DIR__, "src/index.md"))

makedocs(
    sitename = "AdvancedHMC",
    format = Documenter.HTML(),
    modules = [AdvancedHMC]
)

deploydocs(
    repo = "github.com/TuringLang/AdvancedHMC.jl.git",
    push_preview = true, # allow PR to deploy
)


text:
# Format strings used in the stream API.

# Complex 64-bit floats (complex double)
const SOAPY_SDR_CF64 = "CF64"

# Complex 32-bit floats (complex float)
const SOAPY_SDR_CF32 = "CF32"

# Complex signed 32-bit integers (complex int32)
const SOAPY_SDR_CS32 = "CS32"

# Complex unsigned 32-bit integers (complex uint32)
const SOAPY_SDR_CU32 = "CU32"

# Complex signed 16-bit integers (complex int16)
const SOAPY_SDR_CS16 = "CS16"

# Complex unsigned 16-bit integers (complex uint16)
const SOAPY_SDR_CU16 = "CU16"

# Complex signed 12-bit integers (3 bytes)
const SOAPY_SDR_CS12 = "CS12"

# Complex unsigned 12-bit integers (3 bytes)
const SOAPY_SDR_CU12 = "CU12"

# Complex signed 8-bit integers (complex int8)
const SOAPY_SDR_CS8 = "CS8"

# Complex unsigned 8-bit integers (complex uint8)
const SOAPY_SDR_CU8 = "CU8"

# Complex signed 4-bit integers (1 byte)
const SOAPY_SDR_CS4 = "CS4"

# Complex unsigned 4-bit integers (1 byte)
const SOAPY_SDR_CU4 = "CU4"

# Real 64-bit floats (double)
const SOAPY_SDR_F64 = "F64"

# Real 32-bit floats (float)
const SOAPY_SDR_F32 = "F32"

# Real signed 32-bit integers (int32)
const SOAPY_SDR_S32 = "S32"

# Real unsigned 32-bit integers (uint32)
const SOAPY_SDR_U32 = "U32"

# Real signed 16-bit integers (int16)
const SOAPY_SDR_S16 = "S16"

# Real unsigned 16-bit integers (uint16)
const SOAPY_SDR_U16 = "U16"

# Real signed 8-bit integers (int8)
const SOAPY_SDR_S8 = "S8"

# Real unsigned 8-bit integers (uint8)
const SOAPY_SDR_U8 = "U8"

# Get the size of a single element in the specified format.
# param format a supported format string
# return the size of an element in bytes
function SoapySDR_formatToSize(format)
    ccall((:SoapySDR_formatToSize, lib), Cint, (Cstring, ), format)
end


text:
using Documenter, ReachabilityModels
using ReachabilityModels: generate_summary

DocMeta.setdocmeta!(ReachabilityModels, :DocTestSetup,
                    :(using ReachabilityModels); recursive=true)

# generate notebooks
include("generate.jl")

# generate bibliography
#include("bibliography.jl")

generate_summary()

makedocs(
    sitename = "ReachabilityModels",
    modules = [ReachabilityModels],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/aligned.css"]),
    strict = false,
    pages = [
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Overview" => "overview.md",
        "Linear ODEs" => Any[
            "Overview" =>  "models/linear_overview.md",
            "Beam" =>  "models/beam.md",
            "Building" =>  "models/building.md",
            "CD Player" =>  "models/cdplayer.md",
            "Crane" =>  "models/crane.md",
            "Elipse" =>  "models/elipse.md",
            "Five Dim System" =>  "models/five_dim_sys.md",
            "Fom" =>  "models/fom.md",
            "Heat" =>  "models/heat.md",
            "Helicopter" =>  "models/helicopter.md",
            "ISS" =>  "models/iss.md",
            "MNA1" =>  "models/mna1.md",
            "MNA5" =>  "models/mna5.md",
            "Motor" =>  "models/motor.md",
            "PDE" =>  "models/pde.md",
            "Projectile" =>  "models/projectile.md"
        ],
        "Nonlinear ODEs" => Any[
            "Overview" =>  "models/nonlinear_overview.md",
            "Biomodel 7d" => "models/biomodel7d.md",
            "Biomodel 9d" => "models/biomodel9d.md",
            "Brusselator" => "models/brusselator.md",
            "Bucking column" => "models/buckling_column.md",
            "Coupled VanDerPol" => "models/coupled_vanderpol.md",
            "Double Gyre" => "models/doublegyre.md",
            "Jet Engine" => "models/jet_engine.md",
            "Henon-Helies" => "models/henon_heiles.md",
            "Lotka-Volterra" => "models/lotka_volterra.md",
            "Lorenz" => "models/lorenz.md",
            "Robot Arm" => "models/robot_arm.md",
            "Roessler" => "models/roessler.md",
            "Spring Pendulum" => "models/spring_pendulum.md",
            "Steam Governor" => "models/steam_governor.md",
            "VanDerPol" => "models/vanderpol.md"
        ],
        "Hybrid ODEs" => Any[
            "Overview" =>  "models/hybrid_overview.md",
            "Bouncing ball" => "models/bouncing_ball.md",
            "Nonlinear Bouncing ball" => "models/bouncing_ball_nonlinear.md",
            "Cardiac Cell" => "models/cardiac_cell.md",
            "Clocked Thermostat" => "models/clocked_thermostat.md",
            "Navigation System" => "models/navigation_system.md",
            "Powertrain Control" => "models/powertrain_control.md",
            "Spiking Neuron" => "models/spiking_neuron.md",
            "Thermostat" => "models/thermostat.md"
        ],
        "References" => "references.md",
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/ReachabilityModels.jl",
    push_preview = true
)

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
const WARN_COLOR = :yellow
const INFO_COLOR = :blue
const ERROR_COLOR = :red

function _print_summary_head(io::IO)
    print(io, "SUMMARY (color code: ")
    printstyled(io, "warning", color = WARN_COLOR)
    printstyled(io, ", info", color = INFO_COLOR)
    printstyled(io, ", error", color = ERROR_COLOR)
    println(io, ")")
end

"""
    Print useful info about the given metnet.
"""
function summary(io::IO, metnet::MetNet)
    _print_summary_head(io)
    printstyled(io, "model size: $(size(metnet))", "\n", color = INFO_COLOR)
    _summary_bound_state(io, metnet)
end
summary(metnet::MetNet) = summary(stdout, metnet)

function _summary_bound_state(io::IO, metnet::MetNet)
    PRINT_MAX = 50 # TODO make this global?
    M, N = size(metnet)

    all(iszero(metnet.lb)) && all(iszero(metnet.ub)) &&
        (printstyled(io, "lb and ub boths has only zero elements", "\n", color = ERROR_COLOR); return)

    # single checks
    for (name, col) in zip(["lb", "ub"], [metnet.lb, metnet.ub])
        _print_col_summary(io, col, name; expected_l = N, PRINT_MAX = PRINT_MAX)
    end

    # Counple checks
    line_count = 0
    for (i, rxn) in enumerate(metnet.rxns)
        lb = metnet.lb[i]
        ub = metnet.ub[i]
        lb > ub && (printstyled(io, "rxn($i): ($rxn), lb ($lb) > ub ($ub)",
            "\n", color = ERROR_COLOR); line_count += 1)
        lb > 0.0 && (printstyled(io, "rxn($i): ($rxn), lb ($lb) > 0.0",
            "\n", color = WARN_COLOR); line_count += 1)
        ub < 0.0 && (printstyled(io, "rxn($i): ($rxn), ub ($ub) < 0.0",
            "\n", color = WARN_COLOR);  line_count += 1)
        lb == ub && (printstyled(io, "rxn($i): ($rxn), lb ($lb) == ub ($ub)",
            "\n", color = WARN_COLOR); line_count += 1)
        # TODO implement correctly the rev array
        # (isrev(metnet, i) ⊻ metnet.rev[i]) &&
        #         (printstyled(io, "rxn($i): ($rxn), rev and bounds missmatch", "\n",
        #                 color = ERROR_COLOR); line_count += 1)

        if line_count > PRINT_MAX
            printstyled(io, "PRINT_MAX $PRINT_MAX reached!!! ... ", "\n", color = WARN_COLOR)
            break;
        end
        flush(stdout)
    end

    revscount(metnet) > 0 && printstyled(io, "revscount: $(revscount(metnet))", "\n", color = WARN_COLOR)
    fwds_boundedcount(metnet) > 0 && printstyled(io, "fwds_bounded: $(fwds_boundedcount(metnet))", "\n", color = INFO_COLOR)
    bkwds_boundedcount(metnet) > 0 && printstyled(io, "bkwds_bounded: $(bkwds_boundedcount(metnet))", "\n", color = WARN_COLOR)
    blockscount(metnet) > 0 && printstyled(io, "blocks: $(blockscount(metnet))", "\n", color = WARN_COLOR)
    fixxedscount(metnet) > 0 && printstyled(io, "fixxed: $(fixxedscount(metnet))", "\n", color = INFO_COLOR)

    return nothing
end

function _print_col_summary(io::IO, col, name;
        expected_l = length(col),
        PRINT_MAX = 50)

        length(col) != expected_l && (printstyled(io,
            " $name: ($(length(col))) != N ($expected_l), dimention missmatch",
            "\n", color = ERROR_COLOR); return)

        unique_ = sort!(unique(col))
        length(unique_) < PRINT_MAX ?
            printstyled(io, " $name: $(length(unique_)) unique elment(s): ", unique_, "\n", color = INFO_COLOR) :
            printstyled(io, " $name: $(length(unique_)) unique elment(s): min: ",
                first(unique_), " mean: ", mean(col),
                " max: ", last(unique_), "\n", color = INFO_COLOR)

end

function _print_rxn_summary(io::IO, metnet, ider)
    idx = rxnindex(metnet, ider)
    printstyled(io, " rxn[$idx]: ", metnet.rxns[idx], " (", get(metnet.rxnNames, idx, ""), ")\n", color = INFO_COLOR)
    printstyled(io, " lb: ", metnet.lb[idx], ", ub: ", metnet.ub[idx], "\n" , color = INFO_COLOR)
    printstyled(io, " ", rxn_str(metnet, idx), "\n" , color = INFO_COLOR)
end

function _print_met_summary(io::IO, metnet, ider)
    idx = metindex(metnet, ider)
    printstyled(io, " met[$idx]: ", metnet.mets[idx], " (", get(metnet.metNames, idx, ""), ")\n", color = INFO_COLOR)
    printstyled(io, " ", balance_str(metnet, ider), color = INFO_COLOR)
end

function summary(io::IO, model::MetNet, ider::IDER_TYPE)
    _print_summary_head(io)
    try
        _print_rxn_summary(io, model, ider)
    catch end
    try
        _print_met_summary(io, model, ider)
    catch end
end
summary(model::MetNet, ider::IDER_TYPE) = summary(stdout, model, ider)

function summary(io::IO, model::MetNet, state::AbstractMetState;
        PRINT_MAX = 50,
        digits = 4)
    _print_summary_head(io)
    _print_state_head(io, state)
    _print_state_stats(io, model, state, PRINT_MAX)
    _print_state_exchanges_info(io, model, state, digits)
end
summary(model::MetNet, state::AbstractMetState) = summary(stdout, model, state)

function _print_state_stats(io::IO, model::MetNet, state::AbstractMetState, PRINT_MAX)
    println(io, "State stats:")
    av_ = av(state)
    va_ = va(state)
    for (name, col) in zip(["av", "va"], [av_, va_])
        _print_col_summary(io, col, name; expected_l = size(model, 2), PRINT_MAX)
    end
end

function _print_state_head(io::IO, state)
    println(io, " '_print_state_head' not implemented for $(typeof(state))!!!")
end

function _print_state_exchanges_info(io::IO, model::MetNet, state::AbstractMetState, digits = 3)
    println(io, "Exchange info:")
    exchs_ = model.rxns[exchanges(model)] |> sort
    for sense in [1, -1]
        for ider in exchs_
            av_ = av(model, state, ider)
            if sense * av_ > 0
                av_ = round(av_, digits = digits)
                va_ = round(va(model, state, ider), digits = digits)
                rstr = rxn_str(model, ider)
                rbounds = round.(bounds(model, ider), digits = digits)
                printstyled(io, " $ider => av: $av_, va: $va_, bounds: $(rbounds), eq: $(rstr)\n", color = INFO_COLOR)
            end
        end
        println(io)
    end
end


text:
using JuMP
using MosekTools
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

"""
Formulate the problem

      min.  sum_{j=0,...,(n-2)/2} s[j] + t[j] + p[j] + q[j]

      s.t.  (1/2, s[j], x[i]+10x[i+1]) \in Qr,
            (1/2, t[j], 5^{1/2}*(x[i+2]-x[i+3])) \in Qr
            (1/2, r[j], (x[i+1]-2x[i+2])) \in Qr
	    (1/2, u[j], 10^{1/4}*(x[i]-10x[i+3])) \in Qr
	    (1/2, p[j], r[j]) \in Qr
	    (1/2, q[j], u[j]) \in Qr,                       j=0,...,(n-2)/2, i = 2j

	    0.1 <= x[i] <= 1.1,                             i=0,2,...,n-2

"""
function chainsing1(M :: Model,n :: Int)
    m = (n-2) >> 1

    @variable(M,x[1:n])
    @variable(M,p[1:m])
    @variable(M,q[1:m])
    @variable(M,r[1:m])
    @variable(M,s[1:m])
    @variable(M,t[1:m])
    @variable(M,u[1:m])

    for j in 1:m
        i = ((j - 1) << 1) + 1

        # s[j] >= (x[i] + 10*x[i+1])^2
        @constraint(M, [ 0.5, s[j], x[i]+10*x[i+1] ]            in MOI.RotatedSecondOrderCone(3))
        # t[j] >= 5*(x[i+2] - x[i+3])^2
        @constraint(M, [ 0.5, t[j], sqrt(5.0)*(x[i+2] - x[i+3]) ] in MOI.RotatedSecondOrderCone(3))
        # r[j] >= (x[i+1] - 2*x[i+2])^2
        @constraint(M, [ 0.5, r[j], x[i+1] - 2.0*x[i+2] ]         in MOI.RotatedSecondOrderCone(3))
        # 1/10 * u[j] >= (x[i] - 10*x[i+3])^2
        @constraint(M, [ 0.5 * 10^(-1.0), u[j], x[i] - 10.0*x[i+3] ] in MOI.RotatedSecondOrderCone(3))
        # p[j] >= r[j]^2
        @constraint(M, [ 0.5, p[j], r[j] ]                      in MOI.RotatedSecondOrderCone(3))
        # q[j] >= u[j]^2
        @constraint(M, [ 0.5, q[j], u[j] ]                      in MOI.RotatedSecondOrderCone(3))
    end

    #0.1 <= x[i] <= 1.1, i=0,2,...,n-2
    for i in 1:2:n-1
        @constraint(M,x[i] <= 1.1)
        @constraint(M,x[i] >= 0.1)
    end

    @objective(M, Min, sum(s) + sum(t) + sum(p) + sum(q))
end

function chainsing2(M :: Model,n :: Int)
    m = (n-2) >> 1

    @variable(M,x[1:n])
    @variable(M,p[1:m])
    @variable(M,q[1:m])
    @variable(M,r[1:m])
    @variable(M,s)
    @variable(M,u[1:m])

    se = Array{Any}(2 * m + 2)

    #@constraint(M, [ 0.5, s,
    #                 @expression(M,    x[i]   + 10*x[i+1]  for i in 1:2:n-3 )...,
    #                 @expression(M, 5*(x[i+2] -    x[i+3]) for i in 1:2:n-3 )... ] in MOI.RotatedSecondOrderCone(2*m+2))
    @constraint(M, @expression(M,  [ 0.5; s; [ x[i]   + 10*x[i+1]  for i in 1:2:n-3 ] ; [ 5*(x[i+2] -    x[i+3]) for i in 1:2:n-3 ] ] ) in MOI.RotatedSecondOrderCone(2*m+2))
    for j in 1:m
        i = ((j - 1) << 1) + 1

        # r[j] >= (x[i+1] - 2*x[i+2])^2
        @constraint(M, [ 0.5, r[j], x[i+1] - 2*x[i+2] ]         in MOI.RotatedSecondOrderCone(3))
        # u[j] >= sqrt(10)*(x[i] - 10*x[i+3])^2
        @constraint(M, [ 0.5 * 10.0^(-1.0), u[j], x[i] - 10*x[i+3] ] in MOI.RotatedSecondOrderCone(3))
        # p[j] >= r[j]^2
        @constraint(M, [ 0.5, p[j], r[j] ]                      in MOI.RotatedSecondOrderCone(3))
        # q[j] >= u[j]^2
        @constraint(M, [ 0.5, q[j], u[j] ]                      in MOI.RotatedSecondOrderCone(3))
    end

    #0.1 <= x[i] <= 1.1, i=0,2,...,n-2
    for i in 1:2:n-1
        @constraint(M,x[i] <= 1.1)
        @constraint(M,x[i] >= 0.1)
    end

    @objective(M, Min, sum(s) + sum(p) + sum(q))
end

function chainsing3(M :: Model,n :: Int)
    m = (n-2) >> 1

    @variable(M,x[1:n])
    @variable(M,p[1:m])
    @variable(M,q[1:m])
    @variable(M,r[1:m])
    @variable(M,s)
    @variable(M,u[1:m])

    se = Array{Any}(2 * m + 2)

    @constraint(M, @expression(M, [ 0.5; s;
                                    [ x[i]   + 10*x[i+1]  for i in 1:2:n-3 ] ;
                                    [ 5*(x[i+2] -    x[i+3]) for i in 1:2:n-3 ] ;
                                    p ; q ]) in MOI.RotatedSecondOrderCone(4*m+2))
    for j in 1:m
        i = ((j - 1) << 1) + 1

        # r[j] >= (x[i+1] - 2*x[i+2])^2
        @constraint(M, [ 0.5, r[j], x[i+1] - 2*x[i+2] ]         in MOI.RotatedSecondOrderCone(3))
        # u[j] >= sqrt(10)*(x[i] - 10*x[i+3])^2
        @constraint(M, [ 0.5 * 10.0^(-1.0), u[j], x[i] - 10*x[i+3] ] in MOI.RotatedSecondOrderCone(3))
    end

    #0.1 <= x[i] <= 1.1, i=0,2,...,n-2
    for i in 1:2:n-1
        @constraint(M,x[i] <= 1.1)
        @constraint(M,x[i] >= 0.1)
    end

    @objective(M, Min, s)
end

function main(argv :: Vector{String})
    n = 8192
    backend = :generic
    method = 1
    timefile = nothing
    taskfile = nothing
    opt = false

    for a in argv
        m = match(r"--([^=]+)(?:=(.*))?",a)
        if m == nothing
            n = parse(Int,argv[1])
        else
            key = m.captures[1]
            val = m.captures[2]

            if key == "optimize"
                opt = true
            elseif key == "method"
                method = parse(Int,val)
            elseif key == "timefile"
                timefile = val
            elseif key == "out"
                taskfile = val
            elseif key == "backend"
                if val     == "generic" val == nothing
                    backend = :generic
                elseif val == "mock"
                    backend = :mock
                elseif val == "mosek"
                    backend = :mosek
                end
            elseif key == "help"
                println("Usage: chainsing n [options]")
                println("  n The problem scale, the default being 8192. Must be even.")
                println("Options:")
                println("\t--method=[1|2|3] Which formulation to use")
                println("\t--timefile=filename Which formulation to use")
                println("\t--backend=[generic|mosek|mock] Which backend to use.")
                println("\t\tgeneric Use generic backend, then copy to MOSEK optimizer")
                println("\t\tmosek   Use mosek as backend")
                println("\t\tmock    Use generic backend, then copy to JuMP mock optimizer")
                println("\t--out=filename Write task to this file")
                println("\t--optimize Call optimizer")
                return
            end
        end
    end

    #println("Chainsing")
    #println("\tn        : $n")
    #println("\tmethod   : $method")
    #println("\ttimefile : $timefile")
    #println("\ttaskfile : $taskfile")
    #println("\toptimize : $opt")
    #println("\tbackend  : $backend")

    T = 0.0

    m =
        if  backend == :generic
            solver = Mosek.Optimizer()
            Model( optimizer = solver)
        elseif backend == :mosek
            solver = Mosek.Optimizer()
            Model( mode = JuMP.Direct, backend = solver)
        elseif backend == :mock
            solver = MOIU.MockOptimizer(JuMP.JuMPMOIModel{Float64}())
            Model( optimizer = solver)
        end

    if     backend == :generic
        println("Use: Generic backend, then copy to Mosek")
    elseif backend == :mosek
        println("Use: Mosek directly as backend")
    elseif backend == :mock
        println("Use: Generic backend, then copy to mock solver")
    end

    T0 = time()
    if method == 1
        chainsing1(m,n)
    elseif method == 2
        chainsing2(m,n)
    elseif method == 3
        chainsing3(m,n)
    end

    if backend != :mosek
        MOIU.attachoptimizer!(m)
    end
    T1 = time()
    T = T1-T0

    @printf("Model build time: %.2f secs\n",T)
    if timefile != nothing
        open(timefile, "w") do f
            @printf(f,"%.2f", T)
        end
    end

end

main(ARGS)


text:
using MPI
using StaticArrays
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGMethods

using ClimateMachine.BalanceLaws

import ClimateMachine.BalanceLaws: vars_state, nodal_init_state_auxiliary!

using ClimateMachine.Mesh.Geometry: LocalGeometry

struct GradTestModel{dim, dir} <: BalanceLaw end

vars_state(m::GradTestModel, ::Auxiliary, T) = @vars begin
    a::T
    ∇a::SVector{3, T}
    ∇a_exact::SVector{3, T}
end
vars_state(::GradTestModel, ::Prognostic, T) = @vars()

function nodal_init_state_auxiliary!(
    ::GradTestModel{dim, dir},
    aux::Vars,
    tmp::Vars,
    g::LocalGeometry,
) where {dim, dir}
    x, y, z = g.coord
    if dim == 2
        aux.a = x^2 + y^3 - x * y
        if dir isa EveryDirection
            aux.∇a_exact = SVector(2 * x - y, 3 * y^2 - x, 0)
        elseif dir isa HorizontalDirection
            aux.∇a_exact = SVector(2 * x - y, 0, 0)
        elseif dir isa VerticalDirection
            aux.∇a_exact = SVector(0, 3 * y^2 - x, 0)
        end
    else
        aux.a = x^2 + y^3 + z^2 * y^2 - x * y * z
        if dir isa EveryDirection
            aux.∇a_exact = SVector(
                2 * x - y * z,
                3 * y^2 + 2 * z^2 * y - x * z,
                2 * z * y^2 - x * y,
            )
        elseif dir isa HorizontalDirection
            aux.∇a_exact =
                SVector(2 * x - y * z, 3 * y^2 + 2 * z^2 * y - x * z, 0)
        elseif dir isa VerticalDirection
            aux.∇a_exact = SVector(0, 0, 2 * z * y^2 - x * y)
        end
    end
end

using Test
function test_run(mpicomm, dim, direction, Ne, N, FT, ArrayType)

    brickrange = ntuple(j -> range(FT(0); length = Ne[j] + 1, stop = 3), dim)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, dim),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    model = GradTestModel{dim, direction}()
    dg = DGModel(
        model,
        grid,
        nothing,
        nothing,
        nothing;
        state_gradient_flux = nothing,
    )

    continuous_field_gradient!(
        model,
        dg.state_auxiliary,
        (:∇a,),
        dg.state_auxiliary,
        (:a,),
        grid,
        direction,
    )

    # Wrapping in Array ensure both GPU and CPU code use same approx
    @test Array(dg.state_auxiliary.∇a) ≈ Array(dg.state_auxiliary.∇a_exact)
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 5, 1)
    lvls = 1
    polynomialorder = 4

    @testset for FT in (Float64, Float32)
        @testset for dim in 2:3
            @testset for direction in (
                EveryDirection(),
                HorizontalDirection(),
                VerticalDirection(),
            )
                for l in 1:lvls
                    test_run(
                        mpicomm,
                        dim,
                        direction,
                        ntuple(j -> 2^(l - 1) * numelem[j], dim),
                        polynomialorder,
                        FT,
                        ArrayType,
                    )
                end
            end
        end
    end
end

nothing

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module HermiteGF

  export Hermite, Radial
  export Chebyshev, Uniform
  export interpolate

  " Interpolation type (Hermite or Radial) "
  abstract type InterpolationType end

  " Node positions (Uniform or Chebyshev) "
  abstract type NodesType end

  include("chebyshev.jl")
  include("uniform.jl")
  include("hermite.jl")
  include("radial.jl")
  include("interpolation.jl")

end


text:
module HexaLatticeModel

using StaticArrays

mutable struct ParamDict
    nrows::Int
    ncols::Int
    N::Int
end

function genFunctor(d)
    d = d
    # cell size
    w = 0.5
    h = 2w/sqrt(3)

    function roundN(n, N)
        m = (n - 1 + N) % N
        m + 1
    end

    function getXY(k, d=d)
        m = roundN(k, d.ncols)
        n = roundN(div(k-1, d.ncols)+1, d.nrows)
        if n % 2 == 1
            x = m + 0.5
            y = n*(1.5h)
        else
            x = m
            y = n*(1.5h)
        end
        (x, y)
    end

    function getK( x, y, d=d)
        # must round each x, y
        n = (round(Int,y/(1.5h)) - 1 + d.nrows) % d.nrows + 1
        #
        if n % 2 == 1
            xn = Int(x - 0.5)
        else
            xn = Int(x)
        end
        m = (xn - 1 + d.ncols) % d.ncols + 1
        k = (n - 1) * d.ncols + m
    end

    function getNeighbors( k, d=d)
        (x, y) = getXY(k, d)
        @SVector [
         getK(x-2w,      y, d), getK(x+2w,      y, d),
         getK(x- w, y-1.5h, d), getK(x+ w, y-1.5h, d),
         getK(x- w, y+1.5h, d), getK(x+ w, y+1.5h, d)]
    end

    # polygon for each cell
    function genPoly(k, d=d)
        (x, y) = getXY(k, d)
        (
         @SVector [   x-w,   x,   x+w,    x+w,  x,   x-w,  x-w]
            ,
         @SVector [ y+h/2, y+h, y+h/2, y-h/2, y-h, y-h/2, y+h/2]
        )
    end

    return (getNeighbors, genPoly)
end # of functor

end # of module


text:
(eos::BulkModulusEquation{<:Murnaghan1st})(v) = getparam(eos).b0 + PressureEquation(eos)(v)
function (eos::BulkModulusEquation{<:BirchMurnaghan2nd})(v)
    @unpack v0, b0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return b0 * (7f + 1) * (2f + 1)^_2½
end
# The formula in the Gibbs2 paper is wrong! See the EosFit7c paper!
function (eos::BulkModulusEquation{<:BirchMurnaghan3rd})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return b0 * (2f + 1)^_2½ * (27f^2 / 2 * (b′0 - 4) + (3b′0 - 5) * f + 1)
end
function (eos::BulkModulusEquation{<:BirchMurnaghan4th})(v)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return b0 / 6 *
           (2f + 1)^_2½ *
           ((99h - 693b′0 + 1573) * f^3 + (27h - 108b′0 + 105) * f^2 + 6f * (3b′0 - 5) + 6)
end
function (eos::BulkModulusEquation{<:PoirierTarantola2nd})(v)
    @unpack v0, b0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return b0 * (1 - 3f) * exp(-3f)
end
function (eos::BulkModulusEquation{<:PoirierTarantola3rd})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return -b0 / 2 * exp(-3f) * (9f^2 * (b′0 - 2) - 6f * (b′0 + 1) - 2)
end
function (eos::BulkModulusEquation{<:PoirierTarantola4th})(v)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return -b0 / 2 *
           exp(-3f) *
           (9f^3 * (h + 3b′0 + 3) - 9f^2 * (h + 2b′0 + 1) - 6f * (b′0 + 1) - 2)
end
function (eos::BulkModulusEquation{<:Vinet})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    x, ξ = (v / v0)^_⅓, 3 / 2 * (b′0 - 1)
    return -b0 / (2 * x^2) * (3x * (x - 1) * (b′0 - 1) + 2 * (x - 2)) * exp(-ξ * (x - 1))
end
function (eos::BulkModulusEquation{<:AntonSchmidt})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    x, n = v / v0, -b′0 / 2
    return b0 * x^n * (1 + n * log(x))
end


text:
## MEtropolis
export Metropolis
struct Metropolis <: MCMCRule
    num_switches::Int

    # Inner constructor to assert that the number of switches is always odd.
    Metropolis(i) = begin
        if iseven(i)
            @warn """
                Metropolis sampling with an even number of switches per step
                is known to be non-ergodic, because it can potentially never
                change the configuration.
                Setting the number of switches to $(i+1).
                """
            i += 1
        end
        return new(i)
    end
end

"""
    Metropolis([n_switches=1])

Metropolis rule for Markov Chain Monte Carlo sampling where at every step
`n_switches` sites are switched.

The number of switches must be odd, otherwise there might
"""
Metropolis() = Metropolis(1)

function markov_chain_step!(σ, s::MCMCSampler{Metropolis}, net::Union{MatrixNet,KetNet}, c)
    # the denominator
    logψ_σ = 2*real(net(σ))

    for i=1:s.rule.num_switches
        # 1- Flip a spin
        sites_to_flip = 1:nsites(σ)
        flipat = rand(c.rng, sites_to_flip)
        old_val, new_val =flipat!(c.rng, σ, flipat)

        # 2- Compute the probability to be in that state
        logψ_σp = 2*real(net(σ)) #log_prob_ψ(net, σ)
        log_ratio = logψ_σp - logψ_σ
        prob = exp(log_ratio)
        rv = rand(c.rng)

        # If the move is rejected, revert it
        if !(rv < prob)
            setat!(σ, flipat, old_val)
        else # else accept it
            logψ_σ = logψ_σp
            c.steps_accepted += 1
        end
    end
    c.steps_done += 1

    return true
end

function markov_chain_step!(σ::LUState, s::MCMCSampler{Metropolis}, net::Union{MatrixNet,KetNet}, c)
    # Switch num_switches times
    for i=1:s.rule.num_switches
        # 1- Flip a spin
        sites_to_flip = 1:nsites(σ)
        flipat = rand(c.rng, sites_to_flip) #1:nsites(σ))
        old_val, new_val = flipat_fast!(c.rng, σ, flipat)

        # 2- Compute the probability to be in that state
        log_ratio = 2*real(logψ_Δ(net, σ))
        prob = exp(log_ratio)
        rv = rand(c.rng)
        # If the move is rejected, revert it
        if !(rv < prob)
            clear_changes!(σ)
        else
            update_lut!(σ, net)
            c.steps_accepted += 1
        end
    end
    c.steps_done += 1

    return true
end


text:
# set default values for most options
# these values should never really be used because PDESolver should set them

"""
  This function supplies default values for the options dictionary.
  Options already defined are not overwritten.

  The keys it does not supply default values for are:

   * smb_name
   * dmg_name
   * order
   * coloring_distance
   * numBC
   * other boundary condition related keys

  **Inputs**

   * opts: a dictionary
"""
function set_defaults(opts)

  get!(opts, "run_type", 4)
  get!(opts, "verify_coloring", true)
  get!(opts, "use_edge_res", false)
  get!(opts, "write_edge_vertnums", false)
  get!(opts, "write_boundarynums", false)
  get!(opts, "write_dxidx", false)
  get!(opts, "write_jac", false)
  get!(opts, "write_sparsity", false)
  get!(opts, "write_offsets", false)
  get!(opts, "write_counts", false)
  get!(opts, "write_sparsity_nodebnds", false)
  get!(opts, "write_dofs", false)
  get!(opts, "write_interfaces", false)
  get!(opts, "write_boundaries", false)
  get!(opts, "write_sharedboundaries", false)
  get!(opts, "write_face_vertnums", false)
  get!(opts, "write_coords", false)
  get!(opts, "reordering_algorithm", "default")
  get!(opts, "write_element_vertnums", false)
  get!(opts, "exact_visualization", false)
  get!(opts, "error_undefined_bc", true)

  # force the use of the metric calculation method for linear elements
  # typically the curvilinear method is used even for linear (straight-sided)
  # elements.
  get!(opts, "use_linear_metrics", false)
  get!(opts, "reordering_start_coords", [0.0, 0.0, 0.0])

  return nothing
end

"""
  Create an options dictionary with the default options

  **Outputs**

   * opts: the dictionary
"""
function get_defaults()
  opts = Dict{Any, Any}()
  set_defaults(opts)
  return opts
end


text:
abstract AbstractJobSegment

### Factory function
segment(p::AbstractPolicy,args...;keyargs...) = _segment(p,args...;keyargs...)

numproposals(seg::AbstractJobSegment) = numsamples(proposals(seg.samplerstate))

### Interface for AbstractJobSegments
iterate!(seg::AbstractJobSegment,indicatorstate::AbstractSamplerState) = throw(MethodError(iterate!, (seg,indicatorstate)))
prepare!(seg::AbstractJobSegment,indicatorstate::AbstractSamplerState,i::Int) = throw(MethodError(prepare!, (seg,indicatorstate,i)))
tune!(seg::AbstractJobSegment,tvals...) = throw(MethodError(tune!, (seg,tvals...)))
getsamples(seg::AbstractJobSegment,sampleindex) = throw(MethodError(getsamples, (seg,sampleindex)))
getsamplerstatevars(seg::AbstractJobSegment) = throw(MethodError(getsamplerstatevars, (seg)))

### Interface which maps from Future remote references to AbstractJobSegments
@compat function iterate!(remote::Future,indicatorstate::AbstractSamplerState)
    r = fetch(remote)
    if ~isa(r,RemoteException)
        iterate!(r,indicatorstate)
    else
        throw(r)
    end
end

@compat function prepare!(remote::Future,indicatorstate::AbstractSamplerState,i::Int)
    r = fetch(remote)
    if ~isa(r,RemoteException)
        prepare!(r,indicatorstate,i)
    else
        throw(r)
    end
end

@compat function tune!(remote::Future,tvals...)
    r = fetch(remote)
    if ~isa(r,RemoteException)
        tune!(r,tvals...)
    else
        throw(r)
    end
end

@compat function getsamples(remote::Future,sampleindex)
    r = fetch(remote)
    if ~isa(r,RemoteException)
        getsamples(r,sampleindex)
    else
        throw(r)
    end
end

@compat function getsamplerstatevars(remote::Future)
    r = fetch(remote)
    if ~isa(r,RemoteException)
        getsamplerstatevars(r)
    else
        throw(r)
    end
end


text:
# MultithreadEvaluator worker
# that calculates fitnesses in it's own separate thread
mutable struct MTEvaluatorWorker
    task::Task
    num_evals::Int
    is_available::Bool

    MTEvaluatorWorker(task::Task) = new(task, 0, false)
end

# check whether the worker task is running
isactive(worker::MTEvaluatorWorker) =
    istaskstarted(worker.task) && !(istaskdone(worker.task) || istaskfailed(worker.task))
# check whether the worker is waiting for the job to be assigned
isavailable(worker::MTEvaluatorWorker) = worker.is_available

# Job object returned by async_update_fitness() method of MultithreadEvaluator
mutable struct MTFitnessEvalJob{FA, NFA, I, S} <: AbstractFitnessEvaluationJob{FA}
    candidates::I   # candidates iterator
    next_candidate::Union{Nothing, Tuple{Candidate{FA}, S}} # next candidate and iter state
    nafitness::NFA

    npending::Threads.Atomic{Int}   # how many candidates are currently being evaluated
    results_lock::Threads.SpinLock
    results::Vector{Candidate{FA}}
    discarded::Bool
end

# all candidates have been iterated over by workers
has_next_candidate(job::MTFitnessEvalJob) = !isnothing(job.next_candidate)
# fitness calculation of all candidates is done
iscomplete(job::MTFitnessEvalJob) = isnothing(job.next_candidate) && (job.npending[] == 0)
Base.isready(job::MTFitnessEvalJob) = iscomplete(job)

# all calculated fitness candidates are claimed by the caller
isclaimed(job::MTFitnessEvalJob) = iscomplete(job) && isempty(job.results)

# wrapper around iterate()
iterate_candidates(candidates, state, nafitness::Nothing) =
    isnothing(state) ? iterate(candidates) : iterate(candidates, state)

# wrapper around iterate() that skips candidate, if their fitness
# doens't match nafitness
function iterate_candidates(candidates, state, nafitness)
    next = isnothing(state) ? iterate(candidates) : iterate(candidates, state)
    while !isnothing(next) && !isequal(next[1].fitness, nafitness)
        next = iterate(candidates, next[2])
    end
    return next
end

# pick the next candidate from the iterator, requires Evaluator's jobs_queue_lock
function take_candidate!(job::MTFitnessEvalJob)
    isnothing(job.next_candidate) && error("No candidates to take")
    candi = job.next_candidate[1]
    # increment before iterating to avoid job being complete for a short while
    Threads.atomic_add!(job.npending, 1)
    # prepare next candidate
    job.next_candidate = iterate_candidates(job.candidates, job.next_candidate[2], job.nafitness)
    return candi
end

# put the candidate to the results, requires job's results_lock
function put_candidate!(job::MTFitnessEvalJob, candi::Candidate)
    job.npending[] > 0 || error("Job doesn't expect candidates to be put back")
    push!(job.results, candi)
    # decrement after adding to avoid job being subject for discarding
    Threads.atomic_sub!(job.npending, 1)
    return job
end

"""
Fitness evaluator that asynchronously distributes calculation
among several worker threads.
"""
mutable struct MultithreadEvaluator{F, FA, FS, P<:OptimizationProblem, A<:Archive} <: AbstractAsyncEvaluator{P}
    problem::P          # optimization problem
    archive::A          # archive where good candidates are automatically stored
    num_evals::Int      # fitness evaluations counter
    num_jobs::Int       # counter of completed claimed jobs
    last_fitness::F     # last fitness
    arch_nafitness::FA  # NA fitness

    is_stopping::Bool       # whether the evaluator is in the shutdown sequence
    workers::Vector{MTEvaluatorWorker}
    last_woken_worker::Int  # index of the last worker, which thread was waken by the master

    jobs_queue::Vector{MTFitnessEvalJob} # jobs for workers to pick up
    jobs_queue_lock::Threads.SpinLock
    iterated_jobs::Vector{MTFitnessEvalJob} # jobs with all candidates iterated (some calculations might be pending), but not yet fully claimed
    iterated_jobs_lock::Threads.SpinLock

    resources_lock::Threads.SpinLock
    atomics_pool::Vector{Threads.Atomic{Int}}
    locks_pool::Vector{Threads.SpinLock}
    results_pool::Vector{Vector{Candidate{FA}}}

    function MultithreadEvaluator(
        problem::P, archive::A;
        nworkers::Integer = Threads.nthreads() - 1
    ) where {P<:OptimizationProblem, A<:Archive}
        nworkers > 0 || throw(ArgumentError("nworkers must be positive"))
        nworkers < Threads.nthreads() ||
            throw(ArgumentError("nworkers must be less than threads count ($(Threads.nthreads()))"))
        fs = fitness_scheme(problem)
        F = fitness_type(fs)
        FA = fitness_type(archive)

        eval = new{F, FA, typeof(fs), P, A}(
            problem, archive,
            0, 0, nafitness(fs), nafitness(FA), false,
            Vector{MTEvaluatorWorker}(), 0,
            Vector{MTFitnessEvalJob}(), Threads.SpinLock(),
            Vector{MTFitnessEvalJob}(), Threads.SpinLock(),
            Threads.SpinLock(),
            Vector{Threads.Atomic{Int}}(),
            Vector{Threads.SpinLock}(),
            Vector{Vector{Candidate{FA}}}(),
        )
        create_workers!(eval, nworkers)

        finalizer(_shutdown!, eval)
        return eval
    end
end

MultithreadEvaluator(
    problem::OptimizationProblem;
    nworkers::Integer = Threads.nthreads() - 1,
    archiveCapacity::Integer = 10) =
    MultithreadEvaluator(problem, TopListArchive(fitness_scheme(problem), numdims(problem), archiveCapacity),
                         nworkers=nworkers)

num_evals(eval::MultithreadEvaluator) = eval.num_evals

# FIXME move these method to abstract Evaluator (it would need to support A and FA)
archfitness_type(::Type{<:MultithreadEvaluator{F,FA}}) where {F, FA} = FA
archfitness_type(eval::MultithreadEvaluator) = archfitness_type(typeof(eval))
candidate_type(::Type{T}) where T<:MultithreadEvaluator = Candidate{archfitness_type(T)}
candidate_type(eval::MultithreadEvaluator) = candidate_type(typeof(eval))

nworkers(eval::MultithreadEvaluator) = length(eval.workers)

is_stopping(eval::MultithreadEvaluator) = eval.is_stopping

# what the worker thread should do while waiting for a lock/condition
pause(worker::MTEvaluatorWorker; longer::Bool=false) =
    longer ? yield() : ccall(:jl_cpu_pause, Cvoid, ())

# what the master evaluator thread should do while waiting for a lock/condition
pause(eval::MultithreadEvaluator; longer::Bool=false) =
    longer ? yield() : ccall(:jl_cpu_pause, Cvoid, ());

# Locks SpinLock in a "smart" way: allows yielding to other task from time to
# time, cancels waiting for a lock if the evaluator is being shut down.
# Returns true if the lock was acquired
function smartlock(spinlock::Threads.SpinLock, eval::MultithreadEvaluator,
                   worker::Union{MTEvaluatorWorker, Nothing} = nothing;
                   adaptive::Bool=true)
    locked = false
    threadowner = isnothing(worker) ? eval : worker
    if adaptive
        i = 0
        while !is_stopping(eval) && !(locked = trylock(spinlock))
            pause(threadowner, longer=((i += 1) & 0xFFFF) == 0)
        end
    else
        while !is_stopping(eval) && !(locked = trylock(spinlock))
            pause(threadowner, longer=false)
        end
    end
    return locked
end

# constructs Job object given the candidates iterator.
# Skips candidates with non-natifness, if nafitness !== nothing.
function MTFitnessEvalJob(eval::MultithreadEvaluator, candidates::Any, nafitness::Any)
    next_candidate = iterate_candidates(candidates, nothing, nafitness)
    isnothing(next_candidate) && return nothing

    if smartlock(eval.resources_lock, eval)
        res = MTFitnessEvalJob{archfitness_type(eval), typeof(nafitness),
                               typeof(candidates), typeof(last(next_candidate))}(
            candidates, next_candidate, nafitness,
            isempty(eval.atomics_pool) ? Threads.Atomic{Int}(0) : pop!(eval.atomics_pool),
            isempty(eval.locks_pool) ? Threads.SpinLock() : pop!(eval.locks_pool),
            isempty(eval.results_pool) ? eltype(eval.results_pool)() : pop!(eval.results_pool),
            false)
        unlock(eval.resources_lock)
        return res
    else
        error("Failed to lock MultithreadEvaluator resources")
    end
end

# create MTEvaluatorWorker and assign it to a given tread
# HACK (ab)use julia internals to make sure that workers are spawned on different threads
# HACK "inspired" by enq_work() (base/task.jl) and Channel ctor (base/channels.jl)
function MTEvaluatorWorker(eval::MultithreadEvaluator, workerix::Integer, tid::Integer,
                           readysteadygo::Threads.Atomic{Int})
    task = Task(() -> run_mteval_worker(eval, workerix, readysteadygo))
    task.sticky = true
    ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, tid-1)
    push!(Base.Workqueues[tid], task)
    ccall(:jl_wakeup_thread, Cvoid, (Int16,), (tid - 1) % Int16)
    return MTEvaluatorWorker(task)
end

# creates workers and assigns them to different threads
function create_workers!(eval::MultithreadEvaluator, nworkers::Integer)
    @debug "Initializing $nworkers multithread workers..."
    readysteadygo = Threads.Atomic{Int}(-nworkers)
    eval.workers = [MTEvaluatorWorker(eval, i, i >= Threads.threadid() ? i+1 : i, readysteadygo)
                    for i in 1:nworkers]
    @debug "Waiting for workers initialization..."
    while !any(w -> istaskfailed(w.task) || istaskdone(w.task), eval.workers) &&
        (readysteadygo[] < 0)
        pause(eval, longer=true)
    end
    if !all(isactive, eval.workers) || readysteadygo[] < 0
        @info "MultithreadEvaluator: workers initialization failed, shutting down"
        shutdown!(eval)
    else
        @info "MultithreadEvaluator: $nworkers workers ready, starting"
        Threads.atomic_add!(readysteadygo, 1)
    end
    return eval.workers
end

# Runs the evaluator worker until the evaluator is shut down.
# The workers look and pick up the first job in the jobs_queue and send it
# to calculate_fitness()
function run_mteval_worker(
    eval::MultithreadEvaluator,
    workerix::Int,
    readysteadygo::Threads.Atomic{Int}
)
    @debug "Initializing MultithreadEvaluator worker #$workerix at thread=#$(Threads.threadid())"
    Threads.atomic_add!(readysteadygo, 1)
    while !is_stopping(eval) && readysteadygo[] <= 0
        yield()
    end
    @debug "MultithreadEvaluator worker #$workerix started"
    try
        worker = eval.workers[workerix]
        worker.is_available = true
        i = 0
        nlock_failed = 0
        while !is_stopping(eval)
            queuelock = false
            if isempty(eval.jobs_queue) || islocked(eval.jobs_queue_lock)
                nlock_failed = 0
            else
                if !(queuelock = trylock(eval.jobs_queue_lock))
                    if (nlock_failed += 1) >= 1000000
                        nlock_failed = 0
                        @warn("worker #$workerix: waiting for job for too long (jobs_queue=$(length(eval.jobs_queue)), iterated_jobs=$(length(eval.iterated_jobs)), worker_num_evals=$(worker.num_evals), total_num_evals=$(num_evals(eval)), num_jobs=$(eval.num_jobs))")
                    end
                end
            end
            if queuelock
                nlock_failed = 0
                if !is_stopping(eval) && !isempty(eval.jobs_queue)
                    #@debug "worker #$workerix: got job, calculating"
                    worker.is_available = false
                    calculate_fitness(first(eval.jobs_queue), eval, workerix)
                    worker.is_available = true
                    # jobs_queue unlocked by calculate_fitness()
                else
                    @debug "worker #$workerix: locked empty jobs queue ($(length(eval.jobs_queue))) or stopping evaluator ($(is_stopping(eval))), worker_num_evals=$(worker.num_evals), total_num_evals=$(num_evals(eval)), num_jobs=$(eval.num_jobs)"
                    unlock(eval.jobs_queue_lock)
                end
            else
                pause(worker, longer=((i += 1) & 0xFFFF) == 0)
            end
        end
    catch ex
        @warn "Exception at MultithreadEvaluator worker #$workerix" exception=ex
        worker.is_available = false
        showerror(stderr, ex, catch_backtrace())
        rethrow(ex)
    end
    @debug "MultithreadEvaluator worker #$workerix finished"
    nothing
end

# Calculates the fitness of the next candidate in the job and puts the
# candidate with the evaluated fitness to job.results.
# Moves the job from eval.jobs_queue to eval.iterated_jobs if it has no
# other candidates.
function calculate_fitness(
    job::MTFitnessEvalJob,
    eval::MultithreadEvaluator,
    workerix::Int
)
    worker = eval.workers[workerix]
    #@debug "worker #$workerix calculate(): taking candidate"
    candi = take_candidate!(job)
    @assert !isnothing(candi)
    no_next = !has_next_candidate(job)
    if no_next
        #@debug "worker #$workerix calculate(): removing iterated job from the jobs_queue"
        first_job = popfirst!(eval.jobs_queue)
        @assert first_job === job # the job should be the first in the running list
    end
    unlock(eval.jobs_queue_lock) # release the jobs queue lock as soon as possible

    if no_next
        # move the job with no further candidates to iterated
        #@debug "worker #$workerix calculate(): move the job to iterated"
        if smartlock(eval.iterated_jobs_lock, eval, worker)
            push!(eval.iterated_jobs, job)
            if length(eval.iterated_jobs) > 100
                @warn "$(length(eval.iterated_jobs)) unclaimed iterated job(s)"
            end
            unlock(eval.iterated_jobs_lock)
        elseif !is_stopping(eval)
            error("Failed to lock iterated_jobs")
        end
    end
    #@debug "worker #$workerix calculate(): calculating fitness..."
    eval.last_fitness = candi_fitness = fitness(params(candi), eval.problem)
    candi.fitness = archived_fitness(candi_fitness, eval.archive)
    worker.num_evals += 1

    #@debug "worker #$workerix calculate(): putting the result back to the job..."
    if smartlock(job.results_lock, eval, worker)
        put_candidate!(job, candi)
        unlock(job.results_lock)
    elseif !is_stopping(eval)
        error("Failed to lock job results")
    end
    return nothing
end

# silent (no IO) shutdown of the evaluator, to please the finalizer()
# set the is_stopping flag so that the workers can quit from their locks
function _shutdown!(eval::MultithreadEvaluator)
    eval.is_stopping = true
end

# shutdown the evaluator, automatically called when the error occurs
function shutdown!(eval::MultithreadEvaluator)
    @debug "shutdown!(MultithreadEvaluator)"
    eval.is_stopping && error("Cannot shutdown!(MultithreadEvaluator) twice")
    eval.is_stopping = true
    # make sure all the workers tasks are done
    for (i, worker) in enumerate(eval.workers)
        if isactive(worker)
            @debug "shutdown!(MultithreadEvaluator): worker #$i is still running, waiting..."
            wait(worker.task)
            @debug "shutdown!(MultithreadEvaluator): worker #$i finished"
        end
    end
    @info "shutdown!(MultithreadEvaluator): all $(nworkers(eval)) workers stopped"
    @info "shutdown!(MultithreadEvaluator): num_evals=$(num_evals(eval)), num_jobs=$(eval.num_jobs)"
    @info "shutdown!(MultithreadEvaluator): function evals per worker: $(join([worker.num_evals for worker in eval.workers], ", "))"
end

# Discards the job object releasing the reusable resources back to the
# evaluator's resources pool.
# Called from sync_update_fitness()
function discard!(eval::MultithreadEvaluator, job::MTFitnessEvalJob)
    #@debug "discard!(): discarding job..."
    job.discarded && error("Job is already discarded")
    @assert isclaimed(job)
    job.discarded = true

    #@debug "discard!(): removing the job from iterated_jobs"
    if smartlock(eval.iterated_jobs_lock, eval)
        ix = findfirst(j -> j===job, eval.iterated_jobs)
        @assert !isnothing(ix)
        deleteat!(eval.iterated_jobs, ix)
        unlock(eval.iterated_jobs_lock)
    elseif !is_stopping(eval)
        error("Failed to lock iterated_jobs")
    end

    #@debug "discard!(): return the job resources to the pool"
    if smartlock(eval.resources_lock, eval)
        push!(eval.results_pool, empty!(job.results))
        @assert !islocked(job.results_lock)
        @assert job.npending[] == 0
        push!(eval.locks_pool, job.results_lock)
        push!(eval.atomics_pool, job.npending)
        unlock(eval.resources_lock)
    elseif !is_stopping(eval)
        error("Failed to lock resources")
    end
    #@debug "discard!(): done"
    return nothing
end

# Creates MTFitnessEvalJob object and puts it to the jobs_queue
function async_update_fitness!(eval::MultithreadEvaluator, candidates::Any; force::Bool = false)
    #@debug "async_update_fitness!(): creating job..."
    job = MTFitnessEvalJob(eval, candidates, force ? nothing : eval.arch_nafitness)
    isnothing(job) && return job # empty job, do not process further
    #@debug "async_update_fitness!(): adding the job to the queue"
    if smartlock(eval.jobs_queue_lock, eval)
        push!(eval.jobs_queue, job)
        # hint the scheduler on which thread to run the job
        wake_worker = findnext(isavailable, eval.workers, eval.last_woken_worker+1)
        # if not found, search from the beginning
        isnothing(wake_worker) && (wake_worker = findfirst(isavailable, eval.workers))
        if wake_worker !== nothing
            wake_task = eval.workers[wake_worker].task
            ccall(:jl_wakeup_thread, Cvoid, (Int16,), (Threads.threadid(wake_task) - 1) % Int16)
            eval.last_woken_worker = wake_worker
        end
        unlock(eval.jobs_queue_lock)
    elseif !is_stopping(eval)
        error("Failed to lock jobs_queue")
    end
    #@debug "async_update_fitness!(): done"
    return job
end

sync_update_fitness(f::Any, job::Nothing, eval::MultithreadEvaluator) = nothing

# Waits until MTFitnessEvalJob is fully processed.
# Constantly monitors job.results for the new candidates and adds them
# to the fitnesses archive.
function sync_update_fitness(f::Any, job::MTFitnessEvalJob, eval::MultithreadEvaluator)
    while !is_stopping(eval) && !isclaimed(job) # continue until all candidates processed
        #@debug "sync_update_fitness(): wait until new results: hasnext=$(has_next_candidate(job)), npending=$(job.npending[]), results=$(length(job.results))"
        i = 0
        while !is_stopping(eval) && isempty(job.results)
            if ((i+=1) & 0xFFFF) == 0
                if all(isactive, eval.workers)
                    pause(eval, longer=true)
                else
                    @warn "MultithreadEvaluator: some workers have stopped, shutting down"
                    shutdown!(eval)
                end
            else
                pause(eval, longer=false)
            end
        end
        is_stopping(eval) && continue
        #@debug "sync_update_fitness(): locking results and resources"
        reslock = false
        if (reslock = smartlock(job.results_lock, eval)) && smartlock(eval.resources_lock, eval)
            #@debug "sync_update_fitness(): swapping results with an empty vector"
            @assert !isempty(job.results)
            results = job.results
            job.results = isempty(eval.results_pool) ?
                similar(results, 0) : pop!(eval.results_pool)
            @assert isempty(job.results)
            unlock(job.results_lock) # the workers can put their results
            reslock = false
            unlock(eval.resources_lock)
            #@debug "sync_update_fitness(): adding $(length(results)) candidate(s) to archive"
            for candi in results
                add_candidate!(eval.archive, candi.fitness, candi.params, candi.tag, eval.num_evals += 1)
                f !== nothing && f(candi)
            end
            #@debug "sync_update_fitness(): returning results vector to the pool"
            if smartlock(eval.resources_lock, eval)
                push!(eval.results_pool, empty!(results))
                unlock(eval.resources_lock)
            elseif !is_stopping(eval)
                error("Failed to lock resources")
            end
        elseif !is_stopping(eval)
            error("Failed to lock resources and job results")
        end
        reslock && unlock(job.results_lock)
    end
    #@debug "sync_update_fitness(): job is done, discarding"
    discard!(eval, job)
    eval.num_jobs += 1

    return job.candidates
end

update_fitness!(f::Any, eval::MultithreadEvaluator, candidates::Any; force::Bool=false) =
    sync_update_fitness(f, async_update_fitness!(eval, candidates, force=force), eval)

update_fitness!(f::Any, eval::MultithreadEvaluator, candidate::Candidate; force::Bool=false) =
    update_fitness!(f, eval, (candidate,), force=force)[1]

# WARNING it's not efficient to synchronously calculate single fitness using
# asynchronous `MultithreadEvaluator`
function fitness(params::Individual, eval::MultithreadEvaluator, tag::Int=0)
    candi = candidate_type(eval)(params, -1, eval.arch_nafitness, nothing, tag)
    update_fitness!(eval, (candi,))
    return fitness(candi)
end


text:
"""
    plotSRF(X::AbstractVector, Y::AbstractVector, Z::AbstractVector)

Very simple plotting of surface results
"""
function plotSRF(X::AbstractVector, Y::AbstractVector, Z::AbstractVector)

    # TODO: remove Plots. once all other plotting code is gone

    a = Plots.scatter(X, -Y)
    b = Plots.scatter(Z, -Y)
    c = Plots.scatter(Z, X)

    Plots.plot(a, b, c, layout = @layout([a b c]))

end


text:
function digitalmultroot{S<:Integer,T<:Integer}(n::S, bs::T=10)
    -1 < n && 1 < bs || throw(DomainError())
    ds = n
    pers = 0
    while bs <= ds
        ds = prod(digits(ds, bs))
        pers += 1
    end
    return (pers, ds)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module UAF2

using NLopt, PhageData, PhageFields, Statistics, Random, StatsBase, PyPlot

include("log_approx.jl")
include("import_datasets.jl")
include("model.jl")
include("single_seq.jl")
include("likelihood.jl")
include("preprocess_utils.jl")
include("setters.jl")
include("learn.jl")
include("validation_utils.jl")
include("likelihood_rare.jl")
include("learn_rare.jl")

    module BioseqUtils_dev
        include("BioseqUtils_dev.jl")
    end

end