text:
### quadgk.jl
#
# Copyright (C) 2019 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file integration with QuadGK.jl.
#
### Code:

# Helper functions to handle result of integration with QuadGK when endpoints are
# Measurement objects.

# The integral itself is a Measurement object, propagate its uncertainty.
quadgk_result(integral::Measurement, derivative::Real, a::Measurement) =
    result(integral.val, (1, derivative.val), (integral, a))
quadgk_result(integral::Measurement, derivatives::Tuple, a::Tuple) =
    result(integral.val, (1, value.(derivatives)...), (integral, a...))
# Only endpoints are are Measurement objects.
quadgk_result(integral::Real, derivative::Real, a::Measurement) =
    result(integral, derivative, a)
quadgk_result(integral::Real, derivatives::Tuple, a::Tuple) =
    result(integral, derivatives, a)

# Upper bound is a Measurement.  The derivative is f(b).
function QuadGK.quadgk(f, a::S, b::Measurement{T}; kws...) where {S,T<:AbstractFloat}
    F = promote_type(S, T)
    bval = b.val
    integral = QuadGK.quadgk(f, convert(F, a), convert(F, bval); kws...)
    return (quadgk_result(integral[1], f(bval), b), integral[2])
end

# Lower bound is a Measurement.  The derivative is -f(a).
function QuadGK.quadgk(f, a::Measurement{T}, b::S; kws...) where {S,T<:AbstractFloat}
    F = promote_type(S, T)
    aval = a.val
    integral = QuadGK.quadgk(f, convert(F, aval), convert(F, b); kws...)
    return (quadgk_result(integral[1], -f(aval), a), integral[2])
end

# Both bounds are Measurement's.  Derivatives are -f(a) and f(b).
function QuadGK.quadgk(f, a::Measurement{T}, b::Measurement{S};
                       kws...) where {T<:AbstractFloat,S<:AbstractFloat}
    F = promote_type(T, S)
    aval = a.val
    bval = b.val
    integral = QuadGK.quadgk(f, convert(F, aval), convert(F, bval); kws...)
    return (quadgk_result(integral[1], (-f(aval), f(bval)), (a, b)), integral[2])
end


text:
using PhomericsTest
using Test

@testset "MaxwellSolver.jl" begin
#    @test Matrix(DX(1.0, 2, 2)) == [0.0 0.5 0.0 0.0; -0.5 0.0 0.0 0.0; 0.0 0.0 0.0 0.5; 0.0 0.0 -0.5 0.0]
end

ext: jl


text:
import lib.src.default.print
import lib.src.default.array

class Owner {
    field name, pets

    def new(name) {
        self.name = name
        self.pets = []
    }

    def adopt(pet)
        self.pets:0 <- pet

}

class Pet {
	field name, say_msg
    def new(name, say_msg) {
        self.name = name
        self.say_msg = say_msg
    }
    def speak print_strln(flatten([self.name,' says ',self.say_msg,'!']))

    def eq(other)
        self.name == other.name & self.say_msg == other.say_msg

    def to_s flatten(['Pet[', name, ']'])

}

class Doge : Pet {
	field name, breed, age

	def new(name, breed, age) : (name, 'wooofff') {
        self.breed = breed
        self.age = age
	}
    def speak {
        super.speak()
        print_strln('BBBBAAAAAARRRRRKKKK!')
    }
}

class Kitty : Pet {
    field size
    def new(name, size) : (name, 'meooowww') self.size = size
}

def doge_info(doge) {
    print_strln(flatten([doge.name,' is a ',doge.breed,' at age ',[doge.age],'.']))
}

class Test {
    field map

    def new(pets) {
        self.map = []
        for i=0, i<|pets|, i=i+1
            self.map:0 <- [pets[i], i]
    }

    def get_map() self.map
}

def main {
    methods = flatten([Doge.methods, Doge.super.methods])
    for method : methods
        println(method)

    owner1 = Owner.new('owner1')
	doge1  = Doge.new('doge1', 'Labrador', 2)
	doge2  = Doge.new('doge2', 'Dachshund', 5)
	kitty1 = Kitty.new('kitty1', 3)

    tmp = owner1

    if owner1 == tmp print_strln('True')
    else print_strln('False')

    doge_info(doge1)

    println(Doge.class.class.class.fields)
    ;       Doge  Pet Object (none)
    println(Doge.super.super.super)

    arr = [1,2,2,4]
    arr[2]++
    println(arr)

    owner1.adopt(doge1)
    owner1.adopt(doge2)
    owner1.adopt(kitty1)

    Test test1 = Test.new([doge1, doge2, kitty1])
    println(test1.get_map())
    test1.get_map()[1][0].speak()

    for pet : owner1.pets
        pet.speak()

}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
function convert_notation(s)
    s == "intersection" && return "\$\\cap\$"
    s == "union" && return "\$\\cup\$"
    s == "subtraction" && return "\$-\$"
    s == "complement" && return "\$\\sim\$"
    s
end

function dict2tex(tree)
    !haskey(tree, "data") && return ""
    isempty(tree["children"]) && return " [" * tree["data"] * " ]"
    mapreduce(dict2tex, *, tree["children"], init = " [" * convert_notation(tree["data"])) * "]"
end

function json2tex(infile, outfile; preamble=false)
    tree = JSON.parse(read(infile, String))
    open(outfile, "w") do f
        println(f, "% ", outfile)
        if preamble
            println(f, """
                       \\documentclass{article}
                       \\usepackage{Forest}
                       \\begin{document}
                       """)
        end
        println(f, "\\begin{forest}")
        println(f, dict2tex(tree))
        println(f, "\\end{forest}")
        if preamble
            println(f, "\\end{document}")
        end
    end
end

function tree2tex(tree, outfile::String; preamble=true)
    js = JSON.parse(JSON.json(tree))
    treestr = tree2tex(tree, preamble)
    open(outfile, "w") do io
        println(io, treestr)
    end
    nothing
end

"""
    tree2tex_orderfix(tree, outfile, pc, extr, convp; preamble=true)

Fix order.
"""
function tree2tex_orderfix(tree::CachedCSGNode, outfile, pc, extr, convp; preamble=true)
    js = JSON.parse(JSON.json(toNiceTree(tree)))
    treestr = tree2tex(js, preamble)
    #order fix
    impls2, _ = ransacresult2implicit(pc, extr, convp)
    impls2o, _ = ransacresult2implicit_ordered(pc, extr, convp)
    _, p12, p23 = pairthem(impls2, impls2o, convp)
    for p in p12
        treestr = replace(treestr, p)
    end
    for p in p23
        treestr = replace(treestr, p)
    end
    open(outfile, "w") do io
        println(io, treestr)
    end
    nothing
end

function tree2tex(tree, preamble::Bool)
    js = JSON.parse(JSON.json(tree))
    io = IOBuffer()
    #println(io, "% ", outfile)
    if preamble
        println(io, """
        \\begin{figure}[h]""")
    end
    println(io, """    \\centering
    \\begin{forest}""")
    println(io, "   ",dict2tex(js))
    println(io, "   \\end{forest}")
    if preamble
        println(io, """
        \\caption{CSG tree}
        \\label{fig:tree_}
        \\end{figure}
        """)
    end
    return String(take!(io))
end


text:
function SeisWindowHeaders(in, out; key=[], minval=[], maxval=[], tmin=0,
                           tmax=99999, ntrace=1000000)

    @compat SeisProcessHeaders(in, out, [WindowHeaders],
                              [Dict(:key=>key,:minval=>minval,:maxval=>maxval)],
                               group="some", key=key, ntrace=ntrace,
                               update_tracenum=false)
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d_out = join([DATAPATH out "@data@"])
    filename_h_out = join([DATAPATH out "@headers@"])
    @compat nhead = length(fieldnames(Header))
    stream_h = open(filename_h_out)
    nx = round(Int,filesize(stream_h)/(nhead*4))
    h = GrabHeader(stream_h,1)
    close(stream_h)
    nt = convert(Int64,round((tmax - tmin)/h.d1)) #+ 1
    if nt > h.n1
	nt = h.n1
    end
    extent = Extent(convert(Int32,nt), convert(Int32,nx), convert(Int32,1),
                    convert(Int32,1), convert(Int32,1), convert(Float32,tmin),
                    convert(Float32,1), convert(Float32,0), convert(Float32,0),
                    convert(Float32,0), convert(Float32,h.d1),
                    convert(Float32,1), convert(Float32,1), convert(Float32,1),
                    convert(Float32,1), "Time", "Trace Number", "", "", "", "s",
                    "index", "", "", "", "")
    WriteTextHeader(out,extent,"native_float",4,filename_d_out,filename_h_out)

end

function WindowHeaders(h_in;key=[],minval=[],maxval=[])

    minval = convert(Array{Float32,1},vec(minval))
    maxval = convert(Array{Float32,1},vec(maxval))
    nx = length(h_in)
    key2 = AbstractString[]
    minval2 = Float32[]
    maxval2 = Float32[]
    for ikey=1:length(key)
	if key[ikey] != "t"
	    key2 = push!(key2,key[ikey])
	    minval2 = push!(minval2,minval[ikey])
	    maxval2 = push!(maxval2,maxval[ikey])
	end
    end
    return RejectHeaders(h_in,key2,minval2,maxval2,length(key2),length(h_in))

end

function RejectHeaders(h_in::Array{Header,1}, key::Array{AbstractString,1},
                       minval::Array{Float32,1}, maxval::Array{Float32,1},
                       nkeys, nx)
    h_out = Header[]
    keep = true
    key_val = 0f0
    for j=1:nx
	keep = true
	for ikey=1:nkeys

	    key_val = convert(Float32,getfield(h_in[j],Symbol(key[ikey])))

	    if (key_val < minval[ikey] || key_val > maxval[ikey])
		keep = false
	    end
	end
	if (keep==true)

	    h_out = push!(h_out,h_in[j])
	end
    end
    return h_out
end


text:
export KDEModel

MLJBase.@mlj_model mutable struct KDEModel <: PDFEstimator end

function MLJBase.fit(model::KDEModel, verbosity, X)
    n_vars = size(X, 2)
    kde_model = n_vars == 1 ? kde(eachcol(X)...) : kde(MLJBase.matrix(X))

    fitresult = (kde_model, n_vars)
    cache = nothing
    report = nothing
    fitresult, cache, report
end

function MLJBase.transform(model::KDEModel, fitresult, Xnew)
    kde_model, n_vars = fitresult

    ynew = n_vars == 1 ? pdf(kde_model, eachcol(Xnew)...) : Diagonal(pdf(kde_model, eachcol(Xnew)...)).diag
    ynew = reshape(ynew, size(ynew, 1), 1)
    ynew = DataFrame(ynew, :auto)
end

function MLJBase.fitted_params(model::KDEModel, fitresult)
    kde_model, n_vars = fitresult

    (
        kde_model=kde_model,
        n_vars=n_vars,
    )
end


text:
using MicroscopyTools
using Test

using Tullio
using Statistics
using NDTools
using IndexFunArrays

include("binning.jl")
include("calculation_tools.jl")
include("utils.jl")
include("grid_generation.jl")

@testset "MicroscopyTools.jl" begin
    # Write your tests here.
end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
insts = OrderedDict()
insts["minimal"] = [
    ((3, 4, true, true, true, true),),
    ((3, 4, false, true, true, true),),
    ((3, 4, true, false, true, true),),
    ((3, 4, false, false, true, true),),
    ((3, 4, true, true, false, true),),
    ((3, 4, false, false, false, true),),
    ((3, 4, true, true, false, false),),
    ((3, 4, false, false, false, false),),
    ]
insts["fast"] = [
    ((12, 24, true, true, true, true),),
    ((12, 24, false, true, true, true),),
    ((12, 24, true, false, true, true),),
    ((12, 24, false, false, true, true),),
    ((12, 24, true, true, false, true),),
    ((12, 24, false, false, false, true),),
    ((12, 24, true, true, false, false),),
    ((12, 24, false, false, false, false),),
    ]
insts["various"] = [
    ((5, 100, true, true, true, true),),
    ((5, 100, true, true, false, true),),
    ((10, 200, true, true, true, true),),
    ((10, 200, true, true, false, true),),
    ]
return (MatrixCompletionNative, insts)


text:
using StaticArrays
using FastBEAST
using Plots
plotlyjs()

points2D = [@SVector rand(2) for i=1:100] + [SVector(1.0, 2.0) for i=1:100]

bbox2D = BoundingBox(points2D)
bboxframe2D = getboxframe(bbox2D)

scatter([points2D[i][1] for i=1:length(points2D)],
        [points2D[i][2] for i=1:length(points2D)])
plot!(bboxframe2D[:,1], bboxframe2D[:,2])

##
points3D = [@SVector rand(3) for i=1:100] + [SVector(0.0, 0.0, 0.0) for i=1:100]

bbox3D = BoundingBox(points3D)
bboxframe3D = getboxframe(bbox3D)

scatter([points3D[i][1] for i=1:length(points3D)],
        [points3D[i][2] for i=1:length(points3D)],
        [points3D[i][3] for i=1:length(points3D)],
        markersize = 2)

plot!(  bboxframe3D[:,1], bboxframe3D[:,2], bboxframe3D[:,3],
        line=(4, :black, :solid),
        label = nothing)

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
"""
$(TYPEDEF)

Base for types defining system bath interactions in open quantum system models.
"""
abstract type AbstractInteraction end

"""
$(TYPEDEF)

An object to hold coupling operator and the corresponding bath object.

$(FIELDS)
"""
struct Interaction <: AbstractInteraction
    "system operator"
    coupling::AbstractCouplings
    "bath coupling to the system operator"
    bath::AbstractBath
end

"""
$(TYPEDEF)

A Lindblad operator, define by a rate ``γ`` and corresponding operator ``L```.

$(FIELDS)
"""
struct Lindblad <: AbstractInteraction
    "Lindblad rate"
    γ::Any
    "Lindblad operator"
    L::Any
    "size"
    size::Tuple
end

Lindblad(γ::Number, L::Matrix) = Lindblad((s) -> γ, (s) -> L, size(L))
Lindblad(γ::Number, L) = Lindblad((s) -> γ, L, size(L(0)))
Lindblad(γ, L::Matrix) = Lindblad(γ, (s) -> L, size(L))

function Lindblad(γ, L)
    if !(typeof(γ(0)) <: Number)
        throw(ArgumentError("γ should return a number."))
    end
    if !(typeof(L(0)) <: Matrix)
        throw(ArgumentError("L should return a matrix."))
    end
    Lindblad(γ, L, size(L(0)))
end

Base.size(lind::Lindblad) = lind.size

"""
$(TYPEDEF)

An container for different system-bath interactions.

$(FIELDS)
"""
struct InteractionSet{T <: Tuple}
    "A tuple of Interaction"
    interactions::T
end

InteractionSet(inters::AbstractInteraction...) = InteractionSet(inters)
Base.length(inters::InteractionSet) = Base.length(inters.interactions)
Base.getindex(inters::InteractionSet, key...) =
    Base.getindex(inters.interactions, key...)
Base.iterate(iters::InteractionSet, state=1) =
    Base.iterate(iters.interactions, state)

# The following functions are used to build different Liouvillians from
# the `InteractionSet`. They are not publicly available in the current
# release.
function redfield_from_interactions(iset::InteractionSet, U, Ta, atol, rtol)
    kernels = [build_redfield_kernel(i) for i in iset if !(typeof(i.bath) <: StochasticBath)]
    [RedfieldLiouvillian(kernels, U, Ta, atol, rtol)]
end

function cg_from_interactions(iset::InteractionSet, U, tf, Ta, atol, rtol)
    if Ta === nothing || ndims(Ta) == 0
        kernels = [build_cg_kernel(i, tf, Ta) for i in iset if !(typeof(i.bath) <: StochasticBath)]
    else
        if length(Ta) != length(iset)
            throw(ArgumentError("Ta should have the same length as the interaction sets."))
        end
        kernels = [build_cg_kernel(i, tf, t) for (i, t) in zip(iset, Ta) if !(typeof(i.bath) <: StochasticBath)]
    end
    [CGLiouvillian(kernels, U, atol, rtol)]
end

function ule_from_interactions(iset::InteractionSet, U, Ta, atol, rtol)
    kernels = [build_ule_kernel(i) for i in iset if !(typeof(i.bath) <: StochasticBath)]
    [ULELiouvillian(kernels, U, Ta, atol, rtol)]
end

function davies_from_interactions(iset::InteractionSet, ω_range, lambshift::Bool, lambshift_S)
    davies_list = []
    for i in iset
        coupling = i.coupling
        bath = i.bath
        if !(typeof(bath) <: StochasticBath)
            γfun = build_spectrum(bath)
            Sfun = build_lambshift(ω_range, lambshift, bath, lambshift_S)
            if typeof(bath) <: CorrelatedBath
                push!(davies_list, CorrelatedDaviesGenerator(coupling, γfun, Sfun, build_inds(bath)))
            else
                push!(davies_list, DaviesGenerator(coupling, γfun, Sfun))
            end
        end
    end
    davies_list
end

function onesided_ame_from_interactions(iset::InteractionSet, ω_range, lambshift::Bool, lambshift_S)
    l_list = []
    for i in iset
        coupling = i.coupling
        bath = i.bath
        if !(typeof(bath) <: StochasticBath)
            γfun = build_spectrum(bath)
            Sfun = build_lambshift(ω_range, lambshift, bath, lambshift_S)
            if typeof(i.bath) <: CorrelatedBath
                inds = build_inds(bath)
            else
                inds = ((i, i) for i in 1:length(coupling))
                γfun = SingleFunctionMatrix(γfun)
                Sfun = SingleFunctionMatrix(Sfun)
            end
            push!(l_list, OneSidedAMELiouvillian(coupling, γfun, Sfun, inds))
        end
    end
    l_list
end

function fluctuator_from_interactions(iset::InteractionSet)
    f_list = []
    for i in iset
        coupling = i.coupling
        bath = i.bath
        if typeof(bath) <: EnsembleFluctuator
            num = length(coupling)
            dist = construct_distribution(bath)
            b0 = [x.b for x in bath.f] .* rand([-1, 1], length(dist), num)
            next_τ, next_idx = findmin(rand(dist, num))
            push!(f_list, FluctuatorLiouvillian(coupling, dist, b0, next_idx, next_τ, sum(b0, dims=1)[:]))
        end
    end
    f_list
end

lindblad_from_interactions(iset::InteractionSet) =
    [LindbladLiouvillian([i for i in iset if typeof(i) <: Lindblad])]

function build_redfield_kernel(i::Interaction)
    coupling = i.coupling
    bath = i.bath
    cfun = build_correlation(bath)
    rinds = typeof(cfun) == SingleFunctionMatrix ?
        ((i, i) for i = 1:length(coupling)) : build_inds(bath)
    # the kernels is current set as a tuple
    (rinds, coupling, cfun)
end

function build_cg_kernel(i::Interaction, tf, Ta)
    coupling = i.coupling
    cfun = build_correlation(i.bath)
    Ta = Ta === nothing ? coarse_grain_timescale(i.bath, tf)[1] : Ta
    rinds = typeof(cfun) == SingleFunctionMatrix ?
        ((i, i) for i = 1:length(coupling)) : build_inds(i.bath)
    # the kernels is current set as a tuple
    (rinds, coupling, cfun, Ta)
end

function build_ule_kernel(i::Interaction)
    coupling = i.coupling
    cfun = build_jump_correlation(i.bath)
    rinds = typeof(cfun) == SingleFunctionMatrix ?
        ((i, i) for i = 1:length(coupling)) : build_inds(i.bath)
    # the kernels is current set as a tuple
    (rinds, coupling, cfun)
end


text:
using TaylorSeries: set_variables
using LazySets.Approximations: box_approximation

import TaylorModels
import IntervalArithmetic
const IA = IntervalArithmetic
const zeroI = IA.Interval(0.0) # [0.0, 0.0]
const oneI = IA.Interval(1.0) # [1.0, 1.0]
const symI = IA.Interval(-1.0, 1.0)
@inline zeroBox(m) = IntervalBox(zeroI, m)
@inline symBox(m) = IntervalBox(symI, m)

function _to_hyperrectangle(tTM, xTM, n)
    N = length(xTM)
    SET_TYPE = Hyperrectangle{Float64, SVector{n, Float64}, SVector{n, Float64}}
    Rsets = Vector{ReachSet{<:SET_TYPE}}(undef, N-1)
    @inbounds for i in 1:N-1
        Hi = convert(Hyperrectangle, xTM[i])
        t0 = tTM[i]
        t1 = tTM[i+1]
        Rsets[i] = ReachSet(Hi, t0, t1)
    end
    return Rsets
end

function _to_intervalbox(tTM, xTM, n)
    N = length(xTM)
    SET_TYPE = IA.IntervalBox{n, Float64}
    Rsets = Vector{ReachSet{<:SET_TYPE}}(undef, N-1)
    @inbounds for i in 1:N-1
        Bi = xTM[i]
        t0 = tTM[i]
        t1 = tTM[i+1]
        Rsets[i] = ReachSet(Bi, t0, t1)
    end
    return Rsets
end

function _to_zonotope(tTM, vTM, n)
    N = length(tTM)
    SET_TYPE = Zonotope{Float64}
    Rsets = Vector{ReachSet{<:SET_TYPE}}(undef, N-1)
    @inbounds for i in 1:N-1 # loop over the reach sets
        # pick the i-th Taylor model
        X = vTM[:, i]

        # pick the time domain of the given TM (same in all dimensions)
        Δt = TaylorModels.domain(X[1])

        # evaluate the Taylor model in time, the coefficents are now intervals
        X_Δt = evaluate(X, Δt)

        # builds the associated taylor model for each coordinate j = 1...n
        X̂ = [TaylorModelN(X_Δt[j], X[j].rem, zeroBox(n), symBox(n)) for j in 1:n]

        # floating point rigorous polynomial approximation
        fX̂ = TaylorModels.fp_rpa.(X̂)

        # LazySets can overapproximate a Taylor model with a Zonotope
        Zi = overapproximate(fX̂, Zonotope)
        t0 = tTM[i]
        t1 = tTM[i+1]

        Rsets[i] = ReachSet(Zi, t0, t1)
    end
    return Rsets
end

function post(𝒜::TMJets,
              𝑃::InitialValueProblem{<:Union{BBCS, CBBCS, CBBCCS}, <:LazySet},
              𝑂_global::Options)

    # ==================================
    # Initialization
    # ==================================

    𝑂 = merge(𝒜.options.defaults, 𝑂_global, 𝒜.options.specified)

    # system of ODEs
    f! = 𝑃.s.f
    n = statedim(𝑃)

    # initial time and time horizon
    t0 = 0.0
    T = 𝑂[:T]

    # maximum allowed number of steps
    max_steps = 𝑂[:max_steps]

    # unrap algorithm-specific options
    abs_tol, orderQ, orderT = 𝑂[:abs_tol], 𝑂[:orderQ], 𝑂[:orderT]

    # initial sets
    box_x0 = box_approximation(𝑃.x0)
    q0 = center(box_x0)
    δq0 = IntervalBox(low(box_x0)-q0, high(box_x0)-q0)

    # fix the working variables and maximum order in the global
    # parameters struct (_params_TaylorN_)
    set_variables("x", numvars=length(q0), order=2*orderQ)

    # define the property
    if 𝑂[:mode] == "check"
        property = 𝑂[:property]
    elseif 𝑂[:mode] == "reach"
        property = (t, x) -> true
    end

    # =====================
    # Flowpipe computation
    # =====================

    info("Reachable States Computation...")
    @timing begin
        tTM, xTM, vTM = validated_integ(f!, q0, δq0, t0, T, orderQ, orderT, abs_tol,
                                        maxsteps=max_steps, check_property=property)
    end

    # convert to hyperrectangle and wrap around the reach solution
    output_type = 𝑂[:output_type]
    if output_type == Hyperrectangle
        Rsets = _to_hyperrectangle(tTM, xTM, n)
    elseif output_type == IntervalBox
        Rsets = _to_intervalbox(tTM, xTM, n)
    elseif output_type == Zonotope
        Rsets = _to_zonotope(tTM, vTM, n)
    end

    Rsol = ReachSolution(Flowpipe(Rsets), 𝑂)

    # ===========
    # Projection
    # ===========

    if 𝑂[:project_reachset]
        info("Projection...")
        Rsol = @timing project(Rsol)
    end

    return Rsol
end


text:
import RecipesBase

"""
    plot(s::SeisSplit.Result) -> ::Plots.Plot

Plot the results of shear wave splitting analysis from the `SeisSplit.splitting` function.
"""
plot

RecipesBase.@recipe function f(s::Result)
    grid --> false
    framestyle --> :box
    colorbar := false
    layout := 8
    title := hcat("Input", "Source pol. pre-corr", "Source pol. post-corr",
                  "Fast-slow pre-corr", "Fast-slow corr", "Small eigenvalue",
                  "PM pre-corr", "PM post-corr", "")
    titlefontsize --> 11
    legendfontsize --> 8

    # Input traces in N-E orientation
    n_orig, e_orig = Seis.rotate_through(s.trace1, s.trace2, -s.trace1.sta.azi)

    ## Input traces in original orientation with analysis window highlighted
    # Analysis window in colour behind traces
    # RecipesBase.@series begin
    #     subplot := 1
    #     w1, w2 =s.window_start, s.window_end
    #     []
    # end
    # Traces
    RecipesBase.@series begin
        subplot := 1
        label := [s.trace1.sta.cha s.trace2.sta.cha]
        Seis.times(s.trace1), [Seis.trace(s.trace1) Seis.trace(s.trace2)]
    end
    # Window limits
    RecipesBase.@series begin
        subplot := 1
        seriestype := :vline
        label := ""
        linecolor := :red
        linewidth := 1
        [s.window_start, s.window_end]
    end

    ## Rotated to spol
    RecipesBase.@series begin
        subplot := 2
        spol, spol_90 = Seis.cut.(Seis.rotate_through(n_orig, e_orig, 180-s.spol),
                                 s.window_start, s.window_end)
        label := ""
        linecolor := [:blue :red]
        Seis.times(spol), [Seis.trace(spol), Seis.trace(spol_90)]
    end

    ## Corrected and rotated to spol
    RecipesBase.@series begin
        subplot := 3
        label := ""
        linecolor := [:blue :red]
        n, e = deepcopy.((n_orig, e_orig))
        SeisSplit.apply_split!(Seis.cut!.((n, e), s.window_start, s.window_end)...,
                               s.phi_best, -s.dt_best)
        fast, slow = Seis.rotate_through(n, e, 180-s.spol)
        Seis.times(fast), [Seis.trace(fast), Seis.trace(slow)]
    end

    ## Fast and slow before dt shift
    sfast, sslow = Seis.cut.(Seis.rotate_through(n_orig, e_orig, s.phi_best),
                             s.window_start, s.window_end)
    RecipesBase.@series begin
        subplot := 4
        label := ""
        legend := false
        linecolor := [:blue :red]
        Seis.times(sfast), [Seis.trace(sfast), Seis.trace(sslow)]
    end

    ## Fast and slow after dt shift
    RecipesBase.@series begin
        subplot := 5
        legend := false
        linecolor := [:blue :red]
        xlims := (s.window_start, s.window_end)
        # Fake the shift by changing sample times
        [Seis.times(sfast), Seis.times(sslow) .- s.dt_best], [Seis.trace(sfast), Seis.trace(sslow)]
    end

    ## λ₂ surface
    # 10σ - 100σ contours
    RecipesBase.@series begin
        subplot := 6
        seriestype := :contour
        levels := 10:10:100
        linewidth := 1
        linecolor := :gray
        linealpha := 0.6
        s.dt, s.phi, s.lam2
    end

    # 1σ – 10σ contours
    RecipesBase.@series begin
        subplot := 6
        seriestype := :contour
        levels := 1:10
        linewidth := 1
        linecolor := :black
        s.dt, s.phi, s.lam2
    end
    # 1σ contour
    RecipesBase.@series begin
        subplot := 6
        seriestype := :contour
        levels := [1]
        linewidth := 3
        linecolor := :black
        s.dt, s.phi, s.lam2
    end
    # Best value
    RecipesBase.@series begin
        subplot := 6
        seriestype := :scatter
        xguide --> "dt / s"
        yguide --> "phi / deg"
        xlims := extrema(s.dt)
        ylims := extrema(s.phi)
        label := ""
        [s.dt_best], [s.phi_best]
    end

    ## Particle motion before correction
    amax = 1.2*maximum(x -> sqrt(x[1]^2 + x[2]^2), zip(Seis.trace(e_orig), Seis.trace(n_orig)))
    RecipesBase.@series begin
        subplot := 7
        aspect_ratio := :equal
        label := ""
        linecolor --> :black
        xguide := "East"
        yguide := "North"
        xticks --> nothing
        yticks --> nothing
        xlims := (-amax, amax)
        ylims := (-amax, amax)
        Seis.trace(Seis.cut(e_orig, s.window_start, s.window_end)),
            Seis.trace(Seis.cut(n_orig, s.window_start, s.window_end))
    end

    ## Particle motion after correction
    RecipesBase.@series begin
        subplot := 8
        aspect_ratio := :equal
        label := ""
        linecolor --> :black
        xguide := "East"
        yguide := "North"
        xticks --> nothing
        yticks --> nothing
        xlims := (-amax, amax)
        ylims := (-amax, amax)
        n, e = deepcopy.((n_orig, e_orig))
        SeisSplit.apply_split!(n, e, s.phi_best, -s.dt_best)
        Seis.cut!.((n, e), s.window_start, s.window_end)
        Seis.trace(e), Seis.trace(n)
    end

    ## Empty plot at bottom right
    # RecipesBase.@series begin
    #     subplot := 9
    #     framestyle := :empty
    #     nothing
    # end
end


text:
using PyPlot
using DelimitedFiles

CPU=Sys.cpu_info()[1].model

scalar=readdlm("jlvtriad-scalar.dat",comments=true)'
scalar_shared=readdlm("jlvtriad-scalar-shared.dat",comments=true)'
scalar_avx=readdlm("jlvtriad-scalar-avx.dat",comments=true)'
scalar_shared_avx=readdlm("jlvtriad-scalar-shared-avx.dat",comments=true)'

PyPlot.clf()
PyPlot.semilogx(scalar[1,:],scalar[2,:],"g--",label="Array")
PyPlot.semilogx(scalar_shared[1,:],scalar_shared[2,:],"r--",label="SharedArray")
PyPlot.semilogx(scalar_avx[1,:],scalar_avx[2,:],"g",label="Array + avx")
PyPlot.semilogx(scalar_shared_avx[1,:],scalar_shared_avx[2,:],"r",label="SharedArray+avx")
PyPlot.legend()
PyPlot.grid()
PyPlot.xlabel("Array Size")
PyPlot.ylabel("GFlops/s")
PyPlot.title("SharedArray vs Array (scalar): A[i]=B[i]+C[i]*D[i]\n $(CPU)")
PyPlot.savefig("shared-vs-normal-arrays.pdf")
PyPlot.savefig("shared-vs-normal-arrays.png")

max_stars_repo_name: j-fu/julia-tests

max_issues_repo_name: j-fu/julia-tests

max_forks_repo_name: j-fu/julia-tests


text:
using Test
using ProxyInterfaces

# ProxyInterfaces.@array
# ----------------------

struct MyArray{T}
  value::T
end

ProxyInterfaces.array(::Type{MyArray{T}}) where T = T
ProxyInterfaces.array(myarray::MyArray) = myarray.value
ProxyInterfaces.@array MyArray

myarray = MyArray([11,22,33])
@test myarray.value == [11, 22, 33]
@test 22 in myarray
@test eltype(myarray) == Int
@test length(myarray) == 3
@test size(myarray) == (3,)
@test size(myarray, 1) == 3
@test myarray[1] == 11
@test myarray[end] == 33
@test firstindex(myarray) == 1

@test IndexStyle(typeof(myarray)) == IndexStyle(typeof(myarray.value))

sim1 = similar(myarray, Float32)
sim2 = MyArray(similar(myarray.value, Float32))
@test length(sim1) == length(sim2)
@test eltype(sim1) == eltype(sim2)

@test axes(myarray) == axes(myarray.value)
@test strides(myarray) == strides(myarray.value)

@test isempty(MyArray([]))

# no mutable functionality is defined
@test_throws MethodError myarray[2] = 2
@test_throws MethodError empty!(myarray)
@test_throws MethodError fill!(myarray, 0)

# ProxyInterfaces.@array_mutable
# ------------------------------

struct MyMutableArray{T}
  value::T
end

ProxyInterfaces.array(::Type{MyMutableArray{T}}) where T = T
ProxyInterfaces.array(myarray::MyMutableArray) = myarray.value
ProxyInterfaces.@array_mutable MyMutableArray

myarray_mutable = MyMutableArray([12, 13])
myarray_mutable[2] = 44
@test myarray_mutable[2] == 44

empty!(myarray_mutable)
@test isempty(myarray_mutable)
@test fill!(similar(myarray_mutable, Float32), 0).value == fill!(similar(myarray_mutable.value, Float32), 0)


text:
export KalmanFilter, kalman_filter, filter_likelihood

struct KalmanFilter
    xf
    xp
    Pf
    Pp
    K
    ε
    Σ
    z
end

"""
Apply a Kalman filter specified
by a LinearStateSpaceModel to a time series, `Y`.
"""
function kalman_filter(m::LinearStateSpaceModel,θ,Y,u=fill([0.0],length(Y)))
    T = length(Y)

    μ = initial_mean(m,θ)
    Σ0 = initial_covariance(m,θ)
    F = state_transition_matrix(m,θ)
    G = state_input_matrix(m,θ)
    Q = state_covariance(m,θ)

    Θ = state_noise_transformation_matrix(m,θ)
    Qp = Θ*Q*Θ'

    H = observation_matrix(m,θ)
    R = observation_covariance(m,θ)
    Γ = observation_input_matrix(m,θ)

    d = length(μ)

    dt = eltype(F*μ)
    xf = fill(dt.(μ),T+1)
    xp = fill(zeros(dt,d),T)
    Pf = fill(dt.(Σ0),T+1)
    Pp = fill(zeros(dt,d,d),T)
    Σ = fill(zeros(dt,d,d),T)
    K = fill(zeros(dt,d,d),T)
    ε = fill(zeros(dt,d),T)
    for t in 1:T
        x1,P1,ε1,Σ1,K1,xf1,Pf1 = filter_step(Y[t],u[t],F,G,Qp,H,Γ,R,xf[t],Pf[t])
        xp[t] = x1
        Pp[t] = P1
        ε[t]  = ε1
        Σ[t]  = Σ1
        K[t]  = K1
        xf[t+1]= xf1
        Pf[t+1] = Pf1
    end
    KalmanFilter(xf,xp,Pf,Pp,K,ε,Σ,Y)
end

"""
Apply a Kalman filter specified
by a LinearStateSpaceModel to a time series, `Y`.
"""
function kalman_filter(m::TimeVariantLinearStateSpaceModel,Y,u=fill([0.0],length(Y)))
    T = length(Y)

    μ = m.μ
    Σ0 = m.Σ0
    F = m.F
    G = m.G
    Q = m.Q

    Θ = m.Θ
    Qp = [Θ[t]*Q[t]*Θ[t]' for t in eachindex(Θ)]

    H = m.H
    R = m.R
    Γ = m.Γ

    d = length(μ)

    xf = fill(μ,T+1)
    xp = fill(zeros(d),T)
    Pf = fill(Σ0,T+1)
    Pp = fill(zeros(d,d),T)
    K = fill(zeros(d,d),T)
    Σ = fill(zeros(d,d),T)
    ε = fill(zeros(d),T)

    for t in 1:T
        xp[t],Pp[t],ε[t],Σ[t],K[t],xf[t+1],Pf[t+1] = filter_step(Y[t],u[t],F[t],G[t],Qp[t],H[t],Γ[t],R[t],xf[t],Pf[t])
    end
    KalmanFilter(xf,xp,Pf,Pp,K,ε,Σ,Y)
end

function filter_step(Y,u,F,G,Qp,H,Γ,R,xf0,Pf0)
    xp = F*xf0 + G*u
    Pp = F*Pf0*F' .+ Qp

    ε = Y .- H*xp .- Γ*u
    Σ = H*Pp*H'.+R
    K = Pp*H'*inv(Σ)

    xf = xp .+ K * ε
    Pf = (I-K*H)*Pp

    xp,Pp,ε,Σ,K,xf,Pf
end

"""
Calculate the negative log-likelihood
using the innovations representation
"""
function filter_likelihood(kf::KalmanFilter)
    0.5*sum(logdet.(kf.Σ) + [kf.ε[i]'*(kf.Σ[i]\kf.ε[i]) for i in eachindex(kf.ε)] .+ length(kf.ε[1])*log(2pi))
end


text:
# Methods for accessing Cartesian Tuple metadata.

Base.summary(t::CartesianTuple{N,T,U}) where {N,T,U} = "$(N)-dimensional $(T) $(name(U))"
Base.show(io::IO, t::CartesianTuple) = print(io, "$(summary(t)) $(t.coordinates)")

# Iteration, indexing, and field access.

Base.iterate(t::CartesianTuple, state...) = iterate(t.coordinates, state...)
Base.length(t::CartesianTuple{N,T,U}) where {N,T,U} = size(N)
Base.size(::Type{CartesianTuple{N,T,U}}) where {N,T,U} = N
Base.eltype(::Type{CartesianTuple{N,T,U}}) where {N,T,U} = T

Base.getindex(a::CartesianTuple, inds...) = getindex(a.coordinates, inds...)

Base.getproperty(p::CartesianPair, s::Symbol) = @match s begin
        :x => p[1]
        :y => p[2]
        _  => getfield(p, s)
    end

Base.getproperty(t::CartesianTriple, s::Symbol) = @match s begin
        :x => t[1]
        :y => t[2]
        :z => t[3]
        _  => getfield(t, s)
    end

# General operations on Cartesian Tuples.

@inline Base.:(+)(a::CartesianTuple, b::CartesianTuple) = promote_type(typeof(a), typeof(b))(a.coordinates + b.coordinates)

@inline Base.:(-)(a::VectorLike, b::VectorLike) = promote_type(typeof(a), typeof(b))(a.coordinates - b.coordinates)

@inline function Base.:(*)(a::CartesianTuple, b::Number)
                @assert !isnan(b)
                promote_type(typeof(a), typeof(b))(a.coordinates .* b)
        end

@inline Base.:(*)(a::Number, b::CartesianTuple) = b * a

@inline function Base.:(/)(a::CartesianTuple, b::Number)
                @assert b != zero(typeof(b))
                inverse = one(eltype(a)) / b
                a * inverse
        end

@inline Base.:(==)(a::CartesianTuple{N,S,U}, b::CartesianTuple{N,T,U}) where {N,S,T,U} = a.coordinates == b.coordinates

@inline Base.:(==)(a::CartesianTuple{N,S,U}, b::CartesianTuple{N,T,V}) where {N,S,T,U,V} = false

@inline Base.abs(a::CartesianTuple) = typeof(a)(abs.(a.coordinates)...)

@inline Base.:(-)(a::CartesianTuple) = typeof(a)((-).(a.coordinates)...)

# Special behaviors for Point-Vector operations.

@inline Base.:(+)(a::Point{N,S}, b::Vect{N,T}) where {N,S,T} = promote_type(typeof(a), typeof(b))(a.coordinates + b.coordinates)

@inline Base.:(-)(a::Point{N,S}, b::Vect{N,T}) where {N,S,T} = promote_type(typeof(a), typeof(b))(a.coordinates - b.coordinates)


text:
using Test
using MPI
using ClimateMachine.Mesh.BrickMesh

function main()
    MPI.Init()
    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    @assert csize == 3

    (etv, etc, etb, fc) = brickmesh(
        (0:4, 5:9),
        (false, true),
        boundary = ((1, 2), (3, 4), (5, 6)),
        part = crank + 1,
        numparts = csize,
    )
    (etv, etc, etb, fc) = partition(comm, etv, etc, etb, fc)

    if crank == 0
        etv_expect = [
            1 2 7 6 11
            2 3 8 7 12
            6 7 12 11 16
            7 8 13 12 17
        ]
        @test etv == etv_expect
        @test etc[:, :, 1] == [0 1 0 1; 5 5 6 6]
        @test etc[:, :, 2] == [1 2 1 2; 5 5 6 6]
        @test etc[:, :, 3] == [1 2 1 2; 6 6 7 7]
        @test etc[:, :, 4] == [0 1 0 1; 6 6 7 7]
        @test etc[:, :, 5] == [0 1 0 1; 7 7 8 8]
        etb_expect = [
            1 0 0 1 1
            0 0 0 0 0
            0 0 0 0 0
            0 0 0 0 0
        ]
        @test etb == etb_expect
        fc_expect = Array{Int64, 1}[]
        @test fc == fc_expect
    elseif crank == 1
        etv_expect = [
            16 17 12 13 18
            17 18 13 14 19
            21 22 17 18 23
            22 23 18 19 24
        ]
        @test etv == etv_expect
        @test etc[:, :, 1] == [0 1 0 1; 8 8 9 9]
        @test etc[:, :, 2] == [1 2 1 2; 8 8 9 9]
        @test etc[:, :, 3] == [1 2 1 2; 7 7 8 8]
        @test etc[:, :, 4] == [2 3 2 3; 7 7 8 8]
        @test etc[:, :, 5] == [2 3 2 3; 8 8 9 9]
        etb_expect = [
            1 0 0 0 0
            0 0 0 0 0
            0 0 0 0 0
            0 0 0 0 0
        ]
        @test etb == etb_expect
        fc_expect = Array{Int64, 1}[[1, 4, 1, 2], [2, 4, 2, 3], [5, 4, 3, 4]]
        @test fc == fc_expect
    elseif crank == 2
        etv_expect = [
            19 14 9 8 3 4
            20 15 10 9 4 5
            24 19 14 13 8 9
            25 20 15 14 9 10
        ]
        @test etv == etv_expect
        @test etc[:, :, 1] == [3 4 3 4; 8 8 9 9]
        @test etc[:, :, 2] == [3 4 3 4; 7 7 8 8]
        @test etc[:, :, 3] == [3 4 3 4; 6 6 7 7]
        @test etc[:, :, 4] == [2 3 2 3; 6 6 7 7]
        @test etc[:, :, 5] == [2 3 2 3; 5 5 6 6]
        @test etc[:, :, 6] == [3 4 3 4; 5 5 6 6]
        etb_expect = [
            0 0 0 0 0 0
            2 2 2 0 0 2
            0 0 0 0 0 0
            0 0 0 0 0 0
        ]
        @test etb == etb_expect
        fc_expect = Array{Int64, 1}[[1, 4, 4, 5]]
        @test fc == fc_expect
    end
end

main()


text:
# ------------------------------------------------------------------------------------------
# Room shapes
# ------------------------------------------------------------------------------------------
"""
    RoomRectangle(width, height)

Representation of the rectangle room of size (width, height) with origin in (0, 0)
"""
struct RoomRectangle{T<:Real} <: RoomShape
    width::T
    height::T
end

function Base.show(io::IO, r::RoomRectangle)
    print(io, "rectangle (width, height) = $(r.width)m x $(r.height)m")
    return
end

"""
    isvalid(s::RoomRectangle, pos::Point, r = 0)

Checks wheter give point `pos` lies inside of the rectangle room.
"""
function isvalid(s::RoomRectangle, pos::Point, r = 0)
    w, h = s.width, s.height
    x, y = pos
    return r <= x <= w - r && r <= y <= h - r
end

# ------------------------------------------------------------------------------------------
# Simple checkpoints
# ------------------------------------------------------------------------------------------
struct Checkpoint{T<:Real} <: Target
    pos::Point{T}
end

"""
    nearest(c::Checkpoint, pos::Point)

Returns position of the checkpoint `c`.
"""
nearest(c::Checkpoint, ::Point) = c.pos

"""
    distance(p::Pedestrian, c::Checkpoint, pos::Point = p.pos)

Returns euclidian distance of the pedestrian `p` to the checkpoint `c`.
"""
distance(p::Pedestrian, c::Checkpoint, pos::Point = p.pos) = distance(pos, c.pos)

struct HCheckpoint{T<:Real} <: Target
    pos::Point{T}
    width::Float64
end

"""
    nearest(c::HCheckpoint, pos::Point)

Returns the nearest point from the checkpoint `c` to the point `pos`.
"""
function nearest(c::HCheckpoint, pos::Point)
    (xc, yc), w = c.pos, c.width
    x = pos[1]
    return if x < xc
        (xc, yc)
    elseif x > xc + w
        (xc + w, yc)
    else
        (x, yc)
    end
end

"""
    distance(p::Pedestrian, c::HCheckpoint, pos::Point = p.pos)

Returns euclidian distance of the pedestrian `p` to the checkpoint `c`.
"""
function distance(p::Pedestrian, c::HCheckpoint, pos::Point = p.pos)
    return distance(pos, nearest(c, p.pos))
end

# ------------------------------------------------------------------------------------------
# Doors
# ------------------------------------------------------------------------------------------
struct HDoor{T<:Real} <: Door
    pos::Point{T}
    width::Float64
end

"""
    nearest(d::HDoor, pos::Point, r = 0)

Returns the nearest point from the door `d` to the point `pos`. Argument `r` specifies minimum distance from the side of the door.
"""
function nearest(d::HDoor, pos::Point, r = 0)
    (xd, yd), w = d.pos, d.width
    x = pos[1]
    return if x < xd + r
        (xd + r, yd)
    elseif x > xd + w - r
        (xd + w - r, yd)
    else
        (x, yd)
    end
end

"""
    distance(p::Pedestrian, d::HDoor, pos::Point = p.pos)

Returns euclidian distance of the pedestrian `p` to the door `d`.
"""
function distance(p::Pedestrian, d::HDoor, pos::Point = p.pos)
    return distance(pos, nearest(d, p.pos, p.radius))
end

"""
    middle(d::HDoor)

Returns euclidian coordinates of the middle point of the door `d`.
"""
middle(d::HDoor) = (d.pos[1] + d.width/2, d.pos[2])

"""
    isinside(d::HDoor, pos::Point, r = 0)

Checks wheter given position lies inside the door area. Argument `r` specifies minimum distance from the side of the door.
"""
function isinside(d::HDoor, pos::Point, r = 0)
    x, y = pos
    (xd, yd), w = d.pos, d.width

    if xd <= x <= xd + w && yd - r <= y <= yd + r
        if min(distance(d.pos, pos), distance(d.pos .+ (w, 0), pos)) >= r
            return true
        end
    end
    return false
end

"""
    isvalid(d::Door, pos::Point, r = 0)

Checks wheter give position `pos` is valid with respect to pedestrian `p`.
"""
function isvalid(d::Door, pos::Point, r = 0)
    return distance(d.pos, pos) >= r && distance(d.pos .+ (d.width, 0), pos) >= r
end

# ------------------------------------------------------------------------------------------
# Obstacles
# ------------------------------------------------------------------------------------------
struct Rectangle{T<:Real} <: Obstacle
    pos::Point{T}
    width::T
    height::T
end

"""
    isinside(r::Rectangle, pos::Point)

Checks wheter give point `pos` lies inside of the rectangle `r`.
"""
function isinside(r::Rectangle, pos::Point)
    (xr, yr), w, h = r.pos, r.width, r.height
    return xr <= pos[1] <= xr + w && yr <= pos[2] <= yr + h
end

"""
    nearest(r::Rectangle, pos::Point)

Returns the nearest point from the border of the rectangle `r` to the point `pos`.
"""
function nearest(r::Rectangle, pos::Point)
    isinside(r, pos) && throw(ArgumentError("position must be outside of the rectangle"))

    (xr, yr), w, h = r.pos, r.width, r.height
    x, y = pos
    return if x < xr
        if y < yr
            (xr, yr)
        elseif y > yr + h
            (xr, yr + h)
        else
            (xr, y)
        end
    elseif x > xr + w
        if y < yr
            (xr + w, yr)
        elseif y > yr + h
            (xr + w, yr + h)
        else
            (xr + w, y)
        end
    else
        if y < yr
            (x, yr)
        elseif y > yr + h
            (x, yr + h)
        end
    end
end

"""
    isvalid(o::Rectangle, pos::Point, r)

Checks wheter give point `pos` is valid, i.e., if the distance to the nearest point of the rectangle `o` is greater or equal to the minimum distance `r`.
"""
function isvalid(o::Rectangle, pos::Point, r = 0)
    isinside(o, pos) && return false
    return distance(nearest(o, pos), pos) >= r
end


text:
using Test

include("../src/insects.jl")

@testset "insects.jl" begin include("test_insects.jl") end

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
"""Numerically stable symmetric Givens reflection.
Given `a` and `b`, return `(c, s, ρ)` such that

    [ c  s ] [ a ] = [ ρ ]
    [ s -c ] [ b ] = [ 0 ].
"""
function sym_givens(a :: Float64, b :: Float64)
	#
	# Modeled after the corresponding Matlab function by M. A. Saunders and S.-C. Choi.
	# http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf
	# D. Orban, Montreal, May 2015.

  if b == 0.0
    if a == 0.0
      c = 1.0
    else
      c = sign(a)  # In Julia, sign(0) = 0.
    end
    s = 0.0;
    ρ = abs(a);

  elseif a == 0.0
    c = 0.0;
    s = sign(b);
    ρ = abs(b);

  elseif abs(b) > abs(a)
    t = a / b;
    s = sign(b) / sqrt(1.0 + t * t);
    c = s * t;
    ρ = b / s;  # Computationally better than d = a / c since |c| <= |s|.

  else
    t = b / a;
    c = sign(a) / sqrt(1.0 + t * t);
    s = c * t;
    ρ = a / c;  # Computationally better than d = b / s since |s| <= |c|
  end

  return (c, s, ρ)
end

"""Find the real roots of the quadratic

    q(x) = q₂ x² + q₁ x + q₀,

where q₂, q₁ and q₀ are real. Care is taken to avoid numerical
cancellation. Optionally, `nitref` steps of iterative refinement
may be performed to improve accuracy. By default, `nitref=1`.
"""
function roots_quadratic(q₂ :: Float64, q₁ :: Float64, q₀ :: Float64;
                         nitref :: Int=1)
  # Case where q(x) is linear.
  if q₂ == 0.0
    if q₁ == 0.0
      root = [0.0]
      q₀ == 0.0 || (root = Float64[])
    else
      root = [-q₀ / q₁]
    end
    return root
  end

  # Case where q(x) is indeed quadratic.
  rhs = sqrt(eps(Float64)) * q₁ * q₁
  if abs(q₀ * q₂) > rhs
    ρ = q₁ * q₁ - 4.0 * q₂ * q₀
    ρ < 0.0 && return Float64[]
    d = -0.5 * (q₁ + copysign(sqrt(ρ), q₁))
    roots = [d / q₂, q₀ / d]
  else
    # Ill-conditioned quadratic.
    roots = [-q₁ / q₂, 0.0]
  end

  # Perform a few Newton iterations to improve accuracy.
  for k = 1 : 2
    root = roots[k]
    for it = 1 : nitref
      q = (q₂ * root + q₁) * root + q₀
      dq = 2.0 * q₂ * root + q₁
      dq == 0.0 && continue
      root = root - q / dq
    end
    roots[k] = root
  end
  return roots
end

"""Given a trust-region radius `radius`, a vector `x` lying inside the
trust-region and a direction `d`, return `σ1` and `σ2` such that

    ‖x + σi d‖ = radius, i = 1, 2

in the Euclidean norm. If known, ‖x‖² may be supplied in `xNorm2`.

If `flip` is set to `true`, `σ1` and `σ2` are computed such that

    ‖x - σi d‖ = radius, i = 1, 2.
"""
function to_boundary(x :: Vector{Float64}, d :: Vector{Float64},
                     radius :: Float64; flip :: Bool=false, xNorm2 :: Float64=0.0, dNorm2 :: Float64=0.0)
  radius > 0 || error("radius must be positive")

  # ‖d‖² σ² + 2 xᵀd σ + (‖x‖² - radius²).
  xd = dot(x, d)
  flip && (xd = -xd)
  dNorm2 == 0.0 && (dNorm2 = dot(d, d))
  dNorm2 == 0.0 && error("zero direction")
  xNorm2 == 0.0 && (xNorm2 = dot(x, x))
  (xNorm2 <= radius * radius) || error(@sprintf("outside of the trust region: ‖x‖²=%7.1e, Δ²=%7.1e", xNorm2, radius * radius))
  roots = roots_quadratic(dNorm2, 2 * xd, xNorm2 - radius * radius)
  return roots # `σ1` and `σ2`
end

# Call BLAS if possible when using dot, norm, etc.
# Benchmarks indicate that the form BLAS.dot(n, x, 1, y, 1) is substantially faster than BLAS.dot(x, y)

krylov_dot(n :: Int, x :: Vector{T}, dx :: Int, y :: Vector{T}, dy :: Int) where T <: BLAS.BlasReal = BLAS.dot(n, x, dx, y, dy)
krylov_dot(n :: Int, x :: AbstractVector{T}, dx :: Int, y :: AbstractVector{T}, dy :: Int) where T <: Number = dot(x, y)  # ignore dx, dy here

krylov_norm2(n :: Int, x :: Vector{T}, dx :: Int) where T <: BLAS.BlasReal = BLAS.nrm2(n, x, dx)
krylov_norm2(n :: Int, x :: AbstractVector{T}, dx :: Int) where T <: Number = norm(x)  # ignore dx here

krylov_scal!(n :: Int, s :: T, x :: Vector{T}, dx :: Int) where T <: BLAS.BlasReal = BLAS.scal!(n, s, x, dx)
function krylov_scal!(n :: Int, s :: T, x :: AbstractVector{T}, dx :: Int) where T <: Number
  @simd for i = 1:dx:n
    @inbounds x[i] *= s
  end
  return x
end

krylov_axpy!(n :: Int, s :: T, x :: Vector{T}, dx :: Int, y :: Vector{T}, dy :: Int) where T <: BLAS.BlasReal = BLAS.axpy!(n, s, x, dx, y, dy)
function krylov_axpy!(n :: Int, s :: T, x :: AbstractVector{T}, dx :: Int, y :: AbstractVector{T}, dy :: Int) where T <: Number
  # assume dx = dy
  @simd for i = 1:dx:n
    @inbounds y[i] += s * x[i]
  end
  return y
end

krylov_axpby!(n :: Int, s :: T, x :: Vector{T}, dx :: Int, t :: T, y :: Vector{T}, dy :: Int) where T <: BLAS.BlasReal = BLAS.axpby!(n, s, x, dx, t, y, dy)
function krylov_axpby!(n :: Int, s :: T, x :: AbstractVector{T}, dx :: Int, t :: T, y :: AbstractVector{T}, dy :: Int) where T <: Number
  # assume dx = dy
  @simd for i = 1:dx:n
    @inbounds y[i] = s * x[i] + t * y[i]
  end
  return y
end

# the macros are just for readability, so we don't have to write the increments (always equal to 1)

macro kdot(n, x, y)
  return esc(:(krylov_dot($n, $x, 1, $y, 1)))
end

macro knrm2(n, x)
  return esc(:(krylov_norm2($n, $x, 1)))
end

macro kscal!(n, s, x)
  return esc(:(krylov_scal!($n, $s, $x, 1)))
end

macro kaxpy!(n, s, x, y)
  return esc(:(krylov_axpy!($n, $s, $x, 1, $y, 1)))
end

macro kaxpby!(n, s, x, t, y)
  return esc(:(krylov_axpby!($n, $s, $x, 1, $t, $y, 1)))
end


text:
export Image

mutable struct Image
    mHandle::vk.VkImage
    mType::vk.VkImageViewType
    mCreateInfo::vk.VkImageCreateInfo
    mUnowned::Bool
    mDevice::Device
    mMemory::MemoryChunk

    function Image(device::Device, createInfo::vk.VkImageCreateInfo, type::vk.VkImageViewType)
        this = new()
        this.mDevice = device
        this.mType = type
        this.mCreateInfo = createInfo
        this.mUnowned = false

        newImage = Ref{vk.VkImage}()
        infoRef = Ref(createInfo)
        GC.@preserve infoRef begin
            if (vk.vkCreateImage(getLogicalDevice(device), infoRef, C_NULL, newImage) != vk.VK_SUCCESS)
                error("Failed to create image!")
            end
        end
        this.mHandle = newImage[]
        return this
    end
    # Assumes that the Image is managed somewhere else (e.g. in a swapchain)
    # The ImageCreateInfo will not be used to actually create an image, but to
    # derive properties of the image (format, extend, ...)
    function Image(device::Device, createInfo::vk.VkImageCreateInfo, unowned::vk.VkImage, type::vk.VkImageViewType)
        this = new()
        this.mDevice = device
        this.mType = type
        this.mCreateInfo = createInfo
        this.mHandle = unowned
        this.mUnowned = true
        return this
    end
end

@class Image [:handle, :destroy]

function destroy(this::Image)
    if (this.mUnowned)
        return
    end
    vk.vkDestroyImage(getLogicalDevice(this.mDevice), this.mHandle, C_NULL)
    destroy(this.mMemory)
end

function usageToFeatures(usage::vk.VkImageUsageFlags)::vk.VkFormatFeatureFlags
    result::vk.VkFormatFeatureFlags = 0
    if (usage & vk.VK_IMAGE_USAGE_COLOR_ATTACHMENT_BIT != 0)
        result |= vk.VK_FORMAT_FEATURE_COLOR_ATTACHMENT_BIT |
                  vk.VK_FORMAT_FEATURE_COLOR_ATTACHMENT_BLEND_BIT
    end
    if (usage & vk.VK_IMAGE_USAGE_DEPTH_STENCIL_ATTACHMENT_BIT != 0)
        result |= vk.VK_FORMAT_FEATURE_DEPTH_STENCIL_ATTACHMENT_BIT
    end
    if (usage & vk.VK_IMAGE_USAGE_INPUT_ATTACHMENT_BIT != 0)
        result |= vk.VK_FORMAT_FEATURE_COLOR_ATTACHMENT_BIT
    end
    if (usage & vk.VK_IMAGE_USAGE_SAMPLED_BIT != 0)
        result |= vk.VK_FORMAT_FEATURE_SAMPLED_IMAGE_BIT
    end
    if (usage & vk.VK_IMAGE_USAGE_STORAGE_BIT != 0)
        result |= vk.VK_FORMAT_FEATURE_STORAGE_IMAGE_BIT
    end
    if (usage & vk.VK_IMAGE_USAGE_TRANSFER_DST_BIT != 0)
        # not sure if I should use this, because it's hidden behind an
        # extension
        # result |= vk.VK_FORMAT_FEATURE_TRANSFER_DST_BIT_KHR
    end
    if (usage & vk.VK_IMAGE_USAGE_TRANSFER_SRC_BIT != 0)
        # not sure if I should use this, because it's hidden behind an
        # extension
        # result |= vk.VK_FORMAT_FEATURE_TRANSFER_SRC_BIT_KHR
    end
    return result;
end

function createImage(createInfo::ImageCreateInfo, device::Device)::Image
    return Image(device, handleRef(createInfo)[], createInfo.mViewType)
end

function createView(this::Image, range::vk.VkImageSubresourceRange = vk.VkImageSubresourceRange(0, 0, 0, 0, 0))
    return createView(this, this.mType, range)
end

function handle(this::Image)::vk.VkImage
    return this.mHandle
end

function getLogicalDeviceOf(this::Image)::vk.VkDevice
    return getLogicalDevice(this.mDevice)
end

function realizeAttachment(this::Image)
    req = VkExt.getImageMemoryRequirements(getLogicalDevice(this.mDevice), this.mHandle)
    this.mMemory = allocateDedicated(this.mDevice.mSuballocator, req, VRAM)
    bindToImage(this.mMemory, this.mHandle)
end

function getWidth(this::Image)::UInt32
    return this.mCreateInfo.extent.width
end

function getHeight(this::Image)::UInt32
    return this.mCreateInfo.extent.height
end

# TODO
# setDataVRAM
# setDataRAM
# getData
# realizeVRAM
# realizeRAM
# copyFrom
# copyTo
# blitTo
# changeLayout
# changeOwner
# generateMipmaps
# createView
# changeOwnerImpl
# levelPixels
# levelBytes
# prepStagingBuffer


text:
using CutPruners

include("solvers.jl")

isempty(lp_solvers) && warn("Exact Pruning tests not run!")
@testset "Exact Pruning with $solver" for solver in lp_solvers
    optimizer_constructor = optimizer_with_attributes(solver.Optimizer, MOI.Silent() => true)
    for algo in [AvgCutPruningAlgo(20), DecayCutPruningAlgo(20), DeMatosPruningAlgo(20)]
        @testset "Exact pruning" begin
            pruner = CutPruner{2, Int}(algo, :Max)

            # test pruning with one cut
            addcuts!(pruner, [1 0], [0], [true])
            exactpruning!(pruner, optimizer_constructor)
            @test pruner.b == [0]

            # add 10 cuts in a row
            for i in 1:10
                addcuts!(pruner, [1 0], [i], [true])
            end

            exactpruning!(pruner, optimizer_constructor)
            # normally the exact pruning saves only the last cuts
            @test pruner.b == [10]
            # ... and add another set of cuts ...
            for i in 1:10
                addcuts!(pruner, [2 0], [i+10], [true])
            end

            exactpruning!(pruner, optimizer_constructor)
            # perform the same test again
            @test pruner.b == [10, 20]
        end
    end
    for algo in [AvgCutPruningAlgo(20), DecayCutPruningAlgo(20), DeMatosPruningAlgo(20)]
        @testset "Exact pruning" begin
            pruner = CutPruner{2, Int}(algo, :Min)
            # add 10 cuts in a row
            for i in 1:10
                addcuts!(pruner, [1 0], [i], [true])
            end

            exactpruning!(pruner, optimizer_constructor)
            # normally the exact pruning saves only the last cuts
            @test pruner.b == [1]
            # ... and add another set of cuts ...
            for i in 1:10
                addcuts!(pruner, [2 0], [i+10], [true])
            end

            exactpruning!(pruner, optimizer_constructor)
            # perform the same test again
            @test pruner.b == [1, 11]
        end
    end
end


text:
function raindrops(number)
    tst=""
    if (number%3) == 0
        tst="Pling"
    end
    if (number%5) == 0
        tst= tst*"Plang"
    end
    if (number%7) == 0
        tst= tst*"Plong"
    end

    if tst ==""
        tst=string(number)
    end

    return tst
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
# special values

INF(::Type{ArbFloat{P}}) where {P}     =  ArbFloat{P}("Inf")
NAN(::Type{ArbFloat{P}}) where {P}     =  ArbFloat{P}("NaN")
POSINF(::Type{ArbFloat{P}}) where {P}  =  ArbFloat{P}("Inf")
NEGINF(::Type{ArbFloat{P}}) where {P}  =  ArbFloat{P}("-Inf")
ZERO(::Type{ArbFloat{P}}) where {P}    =  ArbFloat{P}(0)
ONE(::Type{ArbFloat{P}}) where {P}     =  ArbFloat{P}(1)
#=
TWO{P}(::Type{ArbFloat{P}})     =  ArbFloat{P}(2)
QRTRPI{P}(::Type{ArbFloat{P}})  =  atan(ArbFloat{P}(1))
PI{P}(::Type{ArbFloat{P}})      =  atan(ArbFloat{P}(1))*4
INVPI{P}(::Type{ArbFloat{P}})   =  ArbFloat{P}( inv(atan(ArbFloat{32+P}(1))*4) )
PHI{P}(::Type{ArbFloat{P}})     =  (sqrt(ArbFloat{P}(5)) + ArbFloat{P}(1)) / ArbFloat{P}(2)
INVPHI{P}(::Type{ArbFloat{P}})  =  (sqrt(ArbFloat{P}(5)) - ArbFloat{P}(1)) / ArbFloat{P}(2)
=#

INF(::Type{T}) where {T <: ArbFloat}     =  INF(ArbFloat{precision(ArbFloat)})
NAN(::Type{T}) where {T <: ArbFloat}     =  NAN(ArbFloat{precision(ArbFloat)})
POSINF(::Type{T}) where {T <: ArbFloat}  =  POSINF(ArbFloat{precision(ArbFloat)})
NEGINF(::Type{T}) where {T <: ArbFloat}  =  NEGINF(ArbFloat{precision(ArbFloat)})
ZERO(::Type{T}) where {T <: ArbFloat}    =  ZERO(ArbFloat{precision(ArbFloat)})
ONE(::Type{T}) where {T <: ArbFloat}     =  ONE(ArbFloat{precision(ArbFloat)})

INF(::Type{ArbFloat}) where {P}     =  ArbFloat{precision(ArbFloat)}("Inf")
NAN(::Type{ArbFloat}) where {P}     =  ArbFloat{precision(ArbFloat)}("NaN")
POSINF(::Type{ArbFloat}) where {P}  =  ArbFloat{precision(ArbFloat)}("Inf")
NEGINF(::Type{ArbFloat}) where {P}  =  ArbFloat{precision(ArbFloat)}("-Inf")
ZERO(::Type{ArbFloat}) where {P}    =  ArbFloat{precision(ArbFloat)}(0)
ONE(::Type{ArbFloat}) where {P}     =  ArbFloat{precision(ArbFloat)}(1)

#=
TWO{T<:ArbFloat}(::Type{T})     =  TWO(ArbFloat{precision(ArbFloat)})
QRTRPI{T<:ArbFloat}(::Type{T})  =  QRTRPI(ArbFloat{precision(ArbFloat)})
PI{T<:ArbFloat}(::Type{T})      =  PI(ArbFloat{precision(ArbFloat)})
PHI{T<:ArbFloat}(::Type{T})     =  PHI(ArbFloat{precision(ArbFloat)})
INVPHI{T<:ArbFloat}(::Type{T})  =  INVPHI(ArbFloat{precision(ArbFloat)})

=#

function PI(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_pi), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
PI(::Type{ArbFloat}) = PI(ArbFloat{precision(ArbFloat)})

function SQRTPI(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_sqrt_pi), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
SQRTPI(::Type{ArbFloat}) = SQRTPI(ArbFloat{precision(ArbFloat)})

function LOG2(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_log2), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
LOG2(::Type{ArbFloat}) = LOG2(ArbFloat{precision(ArbFloat)})

function LOG10(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_log10), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
LOG10(::Type{ArbFloat}) = LOG10(ArbFloat{precision(ArbFloat)})

function EXP1(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_e), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
EXP1(::Type{ArbFloat}) = EXP1(ArbFloat{precision(ArbFloat)})

function EULER(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_euler), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
EULER(::Type{ArbFloat}) = EULER(ArbFloat{precision(ArbFloat)})

function CATALAN(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_catalan), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
CATALAN(::Type{ArbFloat}) = CATALAN(ArbFloat{precision(ArbFloat)})

function KINCHIN(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_kinchin), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
KHINCHIN(::Type{ArbFloat}) = KINCHIN(ArbFloat{precision(ArbFloat)})

function GLAISHER(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_apery), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
GLAISHER(::Type{ArbFloat}) = GLAISHER(ArbFloat{precision(ArbFloat)})

function APERY(::Type{ArbFloat{P}}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_const_apery), Cvoid, (Ref{ArbFloat{P}}, Int), z, P)
    return z
end
APERY(::Type{ArbFloat}) = APERY(ArbFloat{precision(ArbFloat)})


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# int fai_build3(const char *fn, const char *fnfai, const char *fngzi) HTS_RESULT_USED;
@defun fai_build3(fn::Cstring, fnfai::Cstring, fngzi::Cstring)::Cint

# int fai_build(const char *fn) HTS_RESULT_USED;
@defun fai_build(fn::Cstring)::Cint

# void fai_destroy(faidx_t *fai);
@defun fai_destroy(fai::Ptr{faidx_t})::Cvoid

# faidx_t *fai_load3(const char *fn, const char *fnfai, const char *fngzi, int flags);
@defun fai_load3(fn::Cstring, fnfai::Cstring, fngzi::Cstring, flags::Cint)::Ptr{faidx_t}

# faidx_t *fai_load(const char *fn);
@defun fai_load(fn::Cstring)::Ptr{faidx_t}

# faidx_t *fai_load3_format(const char *fn, const char *fnfai, const char *fngzi, int flags, enum fai_format_options format);
@defun fai_load3_format(fn::Cstring, fnfai::Cstring, fngzi::Cstring, flags::Cint, format::fai_format_options)::Ptr{faidx_t}

# faidx_t *fai_load_format(const char *fn, enum fai_format_options format);
@defun fai_load_format(fn::Cstring, format::fai_format_options)::Ptr{faidx_t}

# char *fai_fetch(const faidx_t *fai, const char *reg, int *len);
@defun fai_fetch(fai::Ptr{faidx_t}, reg::Cstring, len::Ptr{Cint})::Ptr{Cchar}

# char *fai_fetch64(const faidx_t *fai, const char *reg, hts_pos_t *len);
@defun fai_fetch64(fai::Ptr{faidx_t}, reg::Cstring, len::Ptr{hts_pos_t})::Ptr{Cchar}

# char *fai_fetchqual(const faidx_t *fai, const char *reg, int *len);
@defun fai_fetchqual(fai::Ptr{faidx_t}, reg::Cstring, len::Ptr{Cint})::Ptr{Cchar}

# char *fai_fetchqual64(const faidx_t *fai, const char *reg, hts_pos_t *len);
@defun fai_fetchqual64(fai::Ptr{faidx_t}, reg::Cstring, len::Ptr{hts_pos_t})::Ptr{Cchar}

# int faidx_fetch_nseq(const faidx_t *fai) HTS_DEPRECATED("Please use faidx_nseq instead");
@defun faidx_fetch_nseq(fai::Ptr{faidx_t})::Cint

# char *faidx_fetch_seq(const faidx_t *fai, const char *c_name, int p_beg_i, int p_end_i, int *len);
@defun faidx_fetch_seq(fai::Ptr{faidx_t}, c_name::Cstring, p_beg_i::Cint, p_end_i::Cint, len::Ptr{Cint})::Ptr{Cchar}

# char *faidx_fetch_seq64(const faidx_t *fai, const char *c_name, hts_pos_t p_beg_i, hts_pos_t p_end_i, hts_pos_t *len);
@defun faidx_fetch_seq64(fai::Ptr{faidx_t}, c_name::Cstring, p_beg_i::hts_pos_t, p_end_i::hts_pos_t, len::Ptr{hts_pos_t})::Ptr{Cchar}

# char *faidx_fetch_qual(const faidx_t *fai, const char *c_name, int p_beg_i, int p_end_i, int *len);
@defun faidx_fetch_qual(fai::Ptr{faidx_t}, c_name::Cstring, p_beg_i::Cint, p_end_i::Cint, len::Ptr{Cint})::Ptr{Cchar}

# char *faidx_fetch_qual64(const faidx_t *fai, const char *c_name, hts_pos_t p_beg_i, hts_pos_t p_end_i, hts_pos_t *len);
@defun faidx_fetch_qual64(fai::Ptr{faidx_t}, c_name::Cstring, p_beg_i::hts_pos_t, p_end_i::hts_pos_t, len::Ptr{hts_pos_t})::Ptr{Cchar}

# int faidx_has_seq(const faidx_t *fai, const char *seq);
@defun faidx_has_seq(fai::Ptr{faidx_t}, seq::Cstring)::Cint

# int faidx_nseq(const faidx_t *fai);
@defun faidx_nseq(fai::Ptr{faidx_t})::Cint

# const char *faidx_iseq(const faidx_t *fai, int i);
@defun faidx_iseq(fai::Ptr{faidx_t}, i::Cint)::Cstring

# int faidx_seq_len(const faidx_t *fai, const char *seq);
@defun faidx_seq_len(fai::Ptr{faidx_t}, seq::Cstring)::Cint

# const char *fai_parse_region(const faidx_t *fai, const char *s, int *tid, hts_pos_t *beg, hts_pos_t *end, int flags);
@defun fai_parse_region(fai::Ptr{faidx_t}, s::Cstring, tid::Ptr{Cint}, beg::Ptr{hts_pos_t}, end_::Ptr{hts_pos_t}, flags::Cint)::Cstring

# void fai_set_cache_size(faidx_t *fai, int cache_size);
@defun fai_set_cache_size(fai::Ptr{faidx_t}, cache_size::Cint)::Cvoid


text:
using CGMethod1D, FastGaussQuadrature
using Test

include("test_interpolation.jl")
include("test_quadrature.jl")


text:
#@testset "test_logger"
using Distributions, InventoryModels
begin
    item1 = Item(Inventory(1, 100), Supplier(100,1), name = "item1", policy = QPolicy())
    item2 = Item(Inventory(1, 50), Supplier(100,1), name = "item2", policy = QPolicy())
    item3 = EndProduct(
        Market(5, CVNormal{0.0}, 5, 0, 10),
        Inventory(50, 15),
        Assembly(100,1, item1 => 1, item2 => 2),
        policy = QPolicy(),
        name = "item3"
        )

    bom = [item1,item2,item3]
    hours_cons = RessourceConstraint(168, item1.sources[1] => 5, item2.sources[1] => 1, item3.sources[1] => 15)
    is = InventorySystem(6, bom, [hours_cons])
    logger = ISLogger(is)
    #=
    is([0,10,0])
    is([0,0,20])
    is([30,0,0])
    is([10,0,10])
    logger(is)
    reset!(is)
    for _ in 1:100
        for it in 1:6
            is(rand(0:20,3))
        end
        logger(is)
        reset!(is)
    end=#
end


text:
# Generic routines to support compilation

##########################################################################
#
#   Types to express evaluation contexts
#
#   A context refers to a specific configuration of the back-end,
#   which can be scalar code, SIMD, CUDA, etc, or even hybrid of them.
#
#   I organize types into a hierarchy, which may simplify later
#   implementation.
#
#   Here, most contexts are empty types. In practice, it is ok to have
#   some information contained in the context
#   (e.g. the capability version of CUDA may be useful for code-gen)
#
##########################################################################

abstract EvalContext

# direct context: directly evaluate the expressions in CPU
abstract DirectContext <: EvalContext

# offshore context: transfer the computation to a specific device (e.g. GPU)
abstract OffshoreContext <: EvalContext

##########################################################################
#
#  Top level compilation skeletons
#
##########################################################################

function compile(ctx::EvalContext, top_expr::Expr)
    # generate codes given an expression

    h = top_expr.head
    if h == :(=) || is_opassign(h) || h == :(block)
        te = texpr(top_expr)
        compile(ctx, te)
    elseif h == :(*=)
        throw(DeError("Devectorize does not support *=, please use .*= for element-wise multiplication."))
    elseif h == :(/=)
        throw(DeError("Devectorize does not support /=, please use ./= for element-wise division."))
    else
        throw(DeError("Top level expression must be either an assignment, op-assignment, or a block."))
    end
end

function add_deps_to_queue(q::Array{TExpr, 1}, ex::TExpr)
    if isa(ex, TFunCall)
        for a in ex.args
            add_deps_to_queue(q, a)
        end
        if ex.deps != nothing && !isempty(ex.deps)
            for d in ex.deps
                add_deps_to_queue(q, d)
                push!(q, d)
            end
        end
    elseif isa(ex, TAssign)
        add_deps_to_queue(q, ex.rhs)
    end
end

function compile(ctx::EvalContext, top_expr::TAssign)
    dep_queue = TExpr[]
    add_deps_to_queue(dep_queue, top_expr)

    if isempty(dep_queue)
        lhs = top_expr.lhs
        rhs = top_expr.rhs

        if is_trivial_assignment(top_expr)
            code_block( assignment(ju_expr(lhs), ju_expr(rhs)) )

        elseif isa(lhs, TGeneralVar)
            # to ensure no alias between left and right hand side
            tmp = gensym("tmp")
            safe_expr = tassign(tvar(tmp), rhs)
            flatten_code_block(
                code_block(compile(ctx, tmode(safe_expr), safe_expr)),
                assignment(ju_expr(lhs), tmp) )

        else
            compile(ctx, tmode(top_expr), top_expr)
        end
    else
        push!(dep_queue, top_expr)
        codes = [compile(ctx, tmode(e), e) for e in dep_queue]
        code_block(codes...)
    end
end

#####################################################################
#
#  Note: the following function:
#
#       compile(ctx, mode, assign_ex)
#
#   should be provided by the back-end
#
#####################################################################

function compile(ctx::EvalContext, ex::TBlock)
    codes = [compile(ctx, s) for s in ex.stmts]
    code_block(codes...)
end


text:
struct BoxFilter <: Filter
    xwidth::Float64
    ywidth::Float64
end

evaluate(f::BoxFilter, x, y) = 1.0

struct TriangleFilter <: Filter
    xwidth::Float64
    ywidth::Float64
end

evaluate(f::TriangleFilter, x, y) = max(0.0, (f.xwidth - abs(x)) * (f.ywidth - abs(y)))

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# This script plots the results from test_flight_interventions.jl

ũ_no_controls = [reshape(u,n,n_s,2) for u in sol_tl_no_controls.u]
cum_I_no_controls =  [sum(u[:,7:8,:]) for u in ũ_no_controls ]
ũ_complete_shutdown = [reshape(u,n,n_s,2) for u in sol_tl_complete_shutdown.u]
cum_I_complete_shutdown =  [sum(u[:,7:8,:]) for u in ũ_complete_shutdown ]
ũ_halved_after_week = [reshape(u,n,n_s,2) for u in sol_tl_halved_after_week.u]
cum_I_halved_after_week =  [sum(u[:,7:8,:]) for u in ũ_halved_after_week ]
ũ_quartered_after_week = [reshape(u,n,n_s,2) for u in sol_tl_quartered_after_week.u]
cum_I_quartered_after_week =  [sum(u[:,7:8,:]) for u in ũ_quartered_after_week ]
ũ_shutdown_after_week = [reshape(u,n,n_s,2) for u in sol_tl_shutdown_after_week.u]
cum_I_shutdown_after_week =  [sum(u[:,7:8,:]) for u in ũ_shutdown_after_week ]

plot(sol_tl_no_controls.t[cum_I_no_controls.>0],cum_I_no_controls[cum_I_no_controls.>0],yaxis=:log,lab="No control measures")
plot!(sol_tl_halved_after_week.t[cum_I_halved_after_week.>0],cum_I_halved_after_week[cum_I_halved_after_week.>0],yaxis=:log,lab="Flights halved after one week")
plot!(sol_tl_quartered_after_week.t[cum_I_quartered_after_week.>0],cum_I_quartered_after_week[cum_I_quartered_after_week.>0],yaxis=:log,lab="Flights quartered after one week")
plot!(sol_tl_shutdown_after_week.t[cum_I_shutdown_after_week.>0],cum_I_shutdown_after_week[cum_I_shutdown_after_week.>0],yaxis=:log,lab="All flights stopped after one week")
plot!(sol_tl_complete_shutdown.t[cum_I_complete_shutdown.>0],cum_I_complete_shutdown[cum_I_complete_shutdown.>0],yaxis=:log,lab="All flights stopped immediately")
xlabel!("Day")
ylabel!("Cumulative incidence")
title!("External prevalence 1e-7")


text:
import InteractiveUtils # stdlib
import Pkg # stdlib
import Test # stdlib

import Random
Random.seed!(999)

@debug(string("Julia depot paths: "), Base.DEPOT_PATH)
@debug(string("Julia load paths: "), Base.LOAD_PATH)

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if !isnothing(logger)
    @debug(string("Julia version info: ",))
    InteractiveUtils.versioninfo(logger.stream; verbose=true)
end

@debug(string("Attempting to import PredictMD...",))
import PredictMD
@debug(string("Successfully imported PredictMD.",))
@debug(string("PredictMD version: "),PredictMD.version(),)
@debug(
    string("PredictMD package directory: "),
    PredictMD.package_directory(),
    )

@debug(string("Attempting to import PredictMDExtra...",))
import PredictMDExtra
@debug(string("Successfully imported PredictMDExtra.",))
@debug(string("PredictMDExtra version: "),PredictMDExtra.version(),)
@debug(
    string("PredictMDExtra package directory: "),
    PredictMDExtra.package_directory(),
    )

@debug(string("Julia depot paths: "), Base.DEPOT_PATH)
@debug(string("Julia load paths: "), Base.LOAD_PATH)

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if !isnothing(logger)
    if ispath(Base.active_project())
        println(logger.stream, "# Location of test environment Project.toml: \"$(Base.active_project())\"")
        println(logger.stream, "# Beginning of test environment Project.toml")
        println(logger.stream, read(Base.active_project(), String))
        println(logger.stream, "# End of test environment Project.toml")
    else
        println(logger.stream, "# File \"$(Base.active_project())\" does not exist")
    end
    if ispath(joinpath(dirname(Base.active_project()), "Manifest.toml"))
        println(logger.stream, "# Location of test environment Manifest.toml: \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\"")
        println(logger.stream, "# Beginning of test environment Manifest.toml")
        println(logger.stream, read(joinpath(dirname(Base.active_project()), "Manifest.toml"),String))
        println(logger.stream, "# End of test environment Manifest.toml")
    else
        println(logger.stream, "# File \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\" does not exist")
    end
end

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if !isnothing(logger)
    if ispath(PredictMD.package_directory("Project.toml"))
        println(logger.stream, "# Location of PredictMD package Project.toml: \"$(PredictMD.package_directory("Project.toml"))\"")
        println(logger.stream, "# Beginning of PredictMD package Project.toml")
        println(logger.stream, read(PredictMD.package_directory("Project.toml"), String))
        println(logger.stream, "# End of PredictMD package Project.toml")
    else
        println(logger.stream, "# File \"$(PredictMD.package_directory("Project.toml"))\" does not exist")
    end
    if ispath(PredictMD.package_directory("Manifest.toml"))
        println(logger.stream, "# Location of PredictMD package Manifest.toml: \"$(PredictMD.package_directory("Manifest.toml"))\"")
        println(logger.stream, "# Beginning of PredictMD package Manifest.toml")
        println(logger.stream, read(PredictMD.package_directory("Manifest.toml"),String))
        println(logger.stream, "# End of PredictMD package Manifest.toml")
    else
        println(logger.stream, "# File \"$(PredictMD.package_directory("Manifest.toml"))\" does not exist")
    end
end

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if !isnothing(logger)
    if ispath(joinpath(dirname(pathof(PredictMDAPI)), "..", "Project.toml"))
        println(logger.stream, "# Location of PredictMDAPI package Project.toml: \"$(joinpath(dirname(pathof(PredictMDAPI)), "..", "Project.toml"))\"")
        println(logger.stream, "# Beginning of PredictMDAPI package Project.toml")
        println(logger.stream, read(joinpath(dirname(pathof(PredictMDAPI)), "..", "Project.toml"), String))
        println(logger.stream, "# End of PredictMDAPI package Project.toml")
    else
        println(logger.stream, "# File \"$(joinpath(dirname(pathof(PredictMDAPI)), "..", "Project.toml"))\" does not exist")
    end
    if ispath(joinpath(dirname(pathof(PredictMDAPI)), "..", "Manifest.toml"))
        println(logger.stream, "# Location of PredictMDAPI package Manifest.toml: \"$(joinpath(dirname(pathof(PredictMDAPI)), "..", "Manifest.toml"))\"")
        println(logger.stream, "# Beginning of PredictMDAPI package Manifest.toml")
        println(logger.stream, read(joinpath(dirname(pathof(PredictMDAPI)), "..", "Manifest.toml"),String))
        println(logger.stream, "# End of PredictMDAPI package Manifest.toml")
    else
        println(logger.stream, "# File \"$(joinpath(dirname(pathof(PredictMDAPI)), "..", "Manifest.toml"))\" does not exist")
    end
end

logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if !isnothing(logger)
    if ispath(PredictMDExtra.package_directory("Project.toml"))
        println(logger.stream, "# Location of PredictMDExtra package Project.toml: \"$(PredictMDExtra.package_directory("Project.toml"))\"")
        println(logger.stream, "# Beginning of PredictMDExtra package Project.toml")
        println(logger.stream, read(PredictMDExtra.package_directory("Project.toml"), String))
        println(logger.stream, "# End of PredictMDExtra package Project.toml")
    else
        println(logger.stream, "# File \"$(PredictMDExtra.package_directory("Project.toml"))\" does not exist")
    end
    if ispath(PredictMDExtra.package_directory("Manifest.toml"))
        println(logger.stream, "# Location of PredictMDExtra package Manifest.toml: \"$(PredictMDExtra.package_directory("Manifest.toml"))\"")
        println(logger.stream, "# Beginning of PredictMDExtra package Manifest.toml")
        println(logger.stream, read(PredictMDExtra.package_directory("Manifest.toml"),String))
        println(logger.stream, "# End of PredictMDExtra package Manifest.toml")
    else
        println(logger.stream, "# File \"$(PredictMDExtra.package_directory("Manifest.toml"))\" does not exist")
    end
end

if group_includes_block(TEST_GROUP, TestBlockUnitTests())
    Test.@testset "Unit tests" begin
        if false
        else
            unit_test_interval_string = "[,)"
        end
        @debug(
            "unit_test_interval_string: ",
            unit_test_interval_string,
            )
        if !is_interval(unit_test_interval_string)
            throw(
                ArgumentError(
                    string(
                        "$(unit_test_interval_string) ",
                        "is not a valid interval",
                        )
                    )
                )
        end
        unit_test_interval = construct_interval(
            unit_test_interval_string
            )
        @debug(
            "unit_test_interval: ",
            unit_test_interval,
            )
        testmodulea_filename = joinpath("TestModuleA","TestModuleA.jl",)
        testmoduleb_filename = joinpath(
            "TestModuleB",
            "directory1",
            "directory2",
            "directory3",
            "directory4",
            "directory5",
            "TestModuleB.jl",
            )
        testmodulec_filename = joinpath(
            PredictMD.maketempdir(),
            "TestModuleC.jl",
            )
        rm(testmodulec_filename; force = true, recursive = true)
        open(testmodulec_filename, "w") do io
            write(io, "module TestModuleC end")
        end
        include(testmodulea_filename)
        include(testmoduleb_filename)
        include(testmodulec_filename)
        test_directory = dirname(@__FILE__)
        unit_test_directory = joinpath(test_directory, "unit")
        include("unit-tests-type-definitions.jl")
        for (root, dirs, files) in walkdir(unit_test_directory)
            Test.@testset "$(root)" begin
                for file in files
                    if endswith(lowercase(strip(file)), ".jl")
                        file_path = joinpath(root, file)
                        if interval_contains_x(unit_test_interval,
                                               strip(file))
                            Test.@testset "$(file_path)" begin
                                @debug("Running $(file_path)")
                                include(file_path)
                            end
                        end
                    end
                end
            end
        end
    end
end

temp_generate_examples_dir = joinpath(
    PredictMD.maketempdir(),
    "generate_examples",
    "PredictMDTEMP",
    "examples",
    )

rm(
    temp_generate_examples_dir;
    force = true,
    recursive = true,
    )

Test.@testset "Integration tests" begin
    Test.@testset "Generate examples      " begin
        Test.@test(
            temp_generate_examples_dir == PredictMD.generate_examples(
                temp_generate_examples_dir;
                scripts = true,
                include_test_statements = true,
                markdown = false,
                notebooks = false,
                execute_notebooks = false,
                )
            )
    end
    Test.@testset "Boston housing regression example (CPU)  " begin
        @debug("Testing Boston housing regression example (CPU)")
        boston_housing_tests = [
            "01_preprocess_data.jl" => TestBlockIntegration1(),
            "02_linear_regression.jl" => TestBlockIntegration1(),
            "03_random_forest_regression.jl" => TestBlockIntegration2(),
            "04_knet_mlp_regression.jl" => TestBlockIntegration2(),
            "05_compare_models.jl" => TestBlockIntegration3(),
            "06_get_model_output.jl" => TestBlockIntegration3(),
            ]
        for test_pair in boston_housing_tests
            test_file = test_pair[1]
            test_block = test_pair[2]
            if group_includes_block(TEST_GROUP, test_block)
                Test.@testset "cpu_examples/bostonhousing/src/$(test_file)" begin
                    include(
                        joinpath(
                            temp_generate_examples_dir,
                            "cpu_examples",
                            "boston_housing",
                            "src",
                            test_file,
                            )
                        )
                end
            end
        end
    end
    Test.@testset "Breast cancer biopsy classification (CPU)" begin
        @debug("Testing breast cancer biopsy classification example (CPU)")
        breast_cancer_tests = [
            "01_preprocess_data.jl" => TestBlockIntegration4(),
            "02_smote.jl" => TestBlockIntegration4(),
            "03_logistic_classifier.jl" => TestBlockIntegration4(),
            "04_random_forest_classifier.jl" => TestBlockIntegration5(),
            "05_c_svc_svm_classifier.jl" => TestBlockIntegration5(),
            "06_nu_svc_svm_classifier.jl" => TestBlockIntegration5(),
            "07_knet_mlp_classifier.jl" => TestBlockIntegration6(),
            "08_compare_models.jl" => TestBlockIntegration7(),
            "09_get_model_output.jl" => TestBlockIntegration7(),
            ]
        for test_pair in breast_cancer_tests
            test_file = test_pair[1]
            test_block = test_pair[2]
            if group_includes_block(TEST_GROUP, test_block)
                Test.@testset "cpu_examples/breastcancer/src/$(test_file)" begin
                    include(
                        joinpath(
                            temp_generate_examples_dir,
                            "cpu_examples",
                            "breast_cancer_biopsy",
                            "src",
                            test_file,
                            )
                        )
                end
            end
        end
    end
end


text:
"""
    solveSOPPeel(num::T,maxWidth::Int;peelSetting::PeelSetting=peellen, rrbSetting::Union{RRBSetting,Nothing}=sop,loggingOn::Bool=true, debugOn::Bool=false, timeLimit::Union{Nothing, T}=nothing, useN::Bool=false, useParallel::Bool=false,fileName::Union{Nothing,String}=nothing) where T<:Int

A method for testing various CVRP solver functions

# Arguments
- `num::T`: which problem from the set to get
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `peelSetting::PeelSetting`:How to pick the peel node
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging statements in other functions
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
- `timeLimit::Union{Nothing, U}`: An optional parameter that adds a time limit
- `useN::Bool`: An optional paramter to use a width of n
- `fileName::Union{Nothing,String}`: The file to write results to
"""
function solveSOPPeel(num::T,maxWidth::Int;peelSetting::PeelSetting=peellen, rrbSetting::Union{RRBSetting,Nothing}=sop,loggingOn::Bool=true, debugOn::Bool=false, timeLimit::Union{Nothing, T}=nothing, useN::Bool=false, useParallel::Bool=false,fileName::Union{Nothing,String}=nothing) where T<:Int
    model = readSOPFile(getSOPFilePath(num))

    if useN
        maxWidth = length(model)
    end
    start = time()
    rrbValues, rrbLastArc = precalculateSOPRRB(model)

    if !isnothing(fileName)
        f = open(fileName,"a")
            write(f, "__________________________________\n\n")
            write(f, string("Peel Setting: ", peelSetting, "\n"))
            write(f, string("Max Width: ", maxWidth, "\n"))
            if useParallel
                write(f, string("Processors: ", Threads.nthreads(), "\n"))
            end
            write(f, "Problem $num: \n")
        close(f)
    end
    if !useParallel
        bestKnownSolution, bestKnownValue, bestBound = peelAndBound(
            model,
            maxWidth,
            orderingType=soporder,
            peelSetting=peelSetting,
            rrbSetting=rrbSetting,
            rrbPacket=(rrbValues, rrbLastArc),
            loggingOn=loggingOn,
            debugOn=debugOn,
            timeElapsed=time()-start,
            timeLimit=timeLimit,
            fileName = fileName
        )
    else
        bestKnownSolution, bestKnownValue, bestBound = parallelPeelAndBound(
            model,
            maxWidth,
            orderingType=soporder,
            peelSetting=peelSetting,
            rrbSetting=rrbSetting,
            rrbPacket=(rrbValues, rrbLastArc),
            loggingOn=loggingOn,
            debugOn=debugOn,
            timeElapsed=time()-start,
            timeLimit=timeLimit,
            fileName = fileName
        )
    end
    if !isnothing(fileName)
        f = open(fileName,"a")
            writeSOPSolution(bestKnownSolution, bestKnownValue, bestBound,fileName)
            write(f, string(bestKnownSolution,",",bestKnownValue,",",bestBound, "\n"))
        close(f)
    end
    return bestKnownSolution, bestKnownValue, bestBound
end

"""
    solveSOPBnB(num::T,maxWidth::Int = 32; rrbSetting::Union{RRBSetting,Nothing}=sop,loggingOn::Bool=true, debugOn::Bool=false, timeLimit::Union{Nothing, T}=nothing, useN::Bool=false,fileName::Union{Nothing,String}=nothing) where T<:Int

A method for testing various CVRP solver functions

# Arguments
- `num::T`: which problem from the set to get
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging statements in other functions
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
- `timeLimit::Union{Nothing, U}`: An optional parameter that adds a time limit
- `useN::Bool`: An optional paramter to use a width of n
- `fileName::Union{Nothing,String}`: The file to write results to
"""
function solveSOPBnB(num::T,maxWidth::Int = 32; rrbSetting::Union{RRBSetting,Nothing}=sop,loggingOn::Bool=true, debugOn::Bool=false, timeLimit::Union{Nothing, T}=nothing, useN::Bool=false,fileName::Union{Nothing,String}=nothing) where T<:Int
    model = readSOPFile(getSOPFilePath(num))

    if useN
        maxWidth = length(model)
    end
    start = time()
    rrbValues, rrbLastArc = precalculateSOPRRB(model)

    if !isnothing(fileName)
        f = open(fileName,"a")
        write(f, "__________________________________\n\n")
            write(f, "Problem $num: \n")
            write(f, string("Max Width: ", maxWidth, "\n"))
        close(f)
    end
    bestKnownSolution, bestKnownValue, bestBound = branchAndBound(
        model,
        maxWidth,
        orderingType=soporder,
        rrbSetting=rrbSetting,
        rrbPacket=(rrbValues, rrbLastArc),
        loggingOn=loggingOn,
        debugOn=debugOn,
        timeElapsed=time()-start,
        timeLimit=timeLimit,
        fileName=fileName
    )
    if !isnothing(fileName)
        f = open(fileName,"a")
            writeSOPSolution(bestKnownSolution, bestKnownValue, bestBound,fileName)
            write(f, string(bestKnownSolution,",",bestKnownValue,",",bestBound, "\n"))
        close(f)
    end
    return bestKnownSolution, bestKnownValue, bestBound
end

function benchmarkSOPs(width::Int, fileName::String, timeLimit::Int;numStart=1,numEnd=41,usePeel::Bool=true, peelSetting::PeelSetting=peellen, useN::Bool=false, useParallel::Bool=false)
    if usePeel
        solveSOPPeel(8,64, loggingOn=true, debugOn=false)
    else
        solveSOPBnB(8,324, loggingOn=true, debugOn=false)
    end
    for i in numStart:numEnd
        if usePeel
            block = @timed solveSOPPeel(i,width,peelSetting=peelSetting, loggingOn=true, debugOn=false, timeLimit=timeLimit, useN=useN, useParallel=useParallel, fileName=fileName)
        else

            block = @timed solveSOPBnB(i,width, loggingOn=true, debugOn=false, timeLimit=timeLimit, fileName=fileName)
        end
        f = open(fileName,"a")
            write(f, string(block[2], " seconds \n"))
        close(f)
    end
end

"""
    getSOPFilePath(num::T;getSet::Bool=false) where T<:Int

The locations of the benchmark problems for convenience

# Arguments
- `num::T`: which problem from the selected set to get
- `getSet::Bool`: An optional parameter that can be used to get an entire set instead of just one
"""
function getSOPFilePath(num::T;getSet::Bool=false) where T<:Int
    #SOPLIB: http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/
    paths = String[
        "./examples/sop/ESC07.sop",#2125 (1)
        "./examples/sop/ESC11.sop",#2075 (2)
        "./examples/sop/ESC12.sop",#1675 (3)
        "./examples/sop/ESC25.sop",#1681 (4)
        "./examples/sop/ESC47.sop",#1288 (5)
        "./examples/sop/ESC63.sop",#62 (6)
        "./examples/sop/ESC78.sop",#18230 (7)
        "./examples/sop/br17.10.sop",#55 (8)
        "./examples/sop/br17.12.sop",#55 (9)
        "./examples/sop/ft53.1.sop",#7531 (10)
        "./examples/sop/ft53.2.sop",#8026 (11)
        "./examples/sop/ft53.3.sop",#10262 (12)
        "./examples/sop/ft53.4.sop",#14425 (13)
        "./examples/sop/ft70.1.sop",#39313 (14)
        "./examples/sop/ft70.2.sop",#[40101,40419]  (15)
        "./examples/sop/ft70.3.sop",#42535  (16)
        "./examples/sop/ft70.4.sop",#53530  (17)
        "./examples/sop/kro124p.1.sop",#[38762,39420]  (18)
        "./examples/sop/kro124p.2.sop",#[39841,41336]  (19)
        "./examples/sop/kro124p.3.sop",#[43904,49499]  (20)
        "./examples/sop/kro124p.4.sop",#[73021,76103]  (21)
        "./examples/sop/p43.1.sop",#28140  (22)
        "./examples/sop/p43.2.sop",#28480  (23)
        "./examples/sop/p43.3.sop",#28835 (24)
        "./examples/sop/p43.4.sop",#83005  (25)
        "./examples/sop/prob.42.sop",#243 (26)
        "./examples/sop/prob.100.sop",#[1045,1163]  (27)
        "./examples/sop/rbg048a.sop",#351 (28)
        "./examples/sop/rbg050c.sop",#467 (29)
        "./examples/sop/rbg109a.sop",#1038 (30)
        "./examples/sop/rbg150a.sop",#1750 (31)
        "./examples/sop/rbg174a.sop",#2033 (32)
        "./examples/sop/rbg253a.sop",#2950 (33)
        "./examples/sop/rbg323a.sop",#3140 (34)
        "./examples/sop/rbg341a.sop",#2568 (35)
        "./examples/sop/rbg358a.sop",#2545 (36)
        "./examples/sop/rbg378a.sop",#[2809,2816]  (37)
        "./examples/sop/ry48p.1.sop",#15805 (38)
        "./examples/sop/ry48p.2.sop",#[16074,16666]  (39)
        "./examples/sop/ry48p.3.sop",#[19490,19894]  (40)
        "./examples/sop/ry48p.4.sop"#31446  (41)
    ]
    if getSet
        return paths
    else
        return paths[num]
    end
end


text:
abstract type AbstractSection end

update!(sh::AbstractSection) = error("Not yet implemented. You can use the dummy function `sh->true`.")

####################################################################################################
@views function sectionShooting(x::AbstractVector, T, normal::AbstractVector, center::AbstractVector)
	N = length(center)
	# we only constrain the first point to lie on a specific hyperplane
	# this avoids the temporary xc - centers
	return (dot(x, normal) - dot(center, normal)) * T
end

# section for Standard Shooting
"""
$(TYPEDEF)

This composite type (named for Section Standard Shooting) encodes a type of section implemented by a hyperplane. It can be used in conjunction with [`ShootingProblem`](@ref). The hyperplane is defined by a point `center` and a `normal`.

$(TYPEDFIELDS)

"""
struct SectionSS{Tn, Tc}  <: AbstractSection
	"Normal to define hyperplane"
	normal::Tn

	"Representative point on hyperplane"
	center::Tc
end

(sect::SectionSS)(u, T) = sectionShooting(u, T, sect.normal, sect.center)

isEmpty(sect::SectionSS{Tn, Tc}) where {Tn, Tc} = (Tn == Nothing) || (Tn == Nothing)

# we update the field of Section, useful during continuation procedure for updating the section
function update!(sect::SectionSS, normal, center)
	copyto!(sect.normal, normal)
	copyto!(sect.center, center)
	sect
end

####################################################################################################
# Poincare shooting based on Sánchez, J., M. Net, B. Garcı́a-Archilla, and C. Simó. “Newton–Krylov Continuation of Periodic Orbits for Navier–Stokes Flows.” Journal of Computational Physics 201, no. 1 (November 20, 2004): 13–33. https://doi.org/10.1016/j.jcp.2004.04.018.

function sectionHyp!(out, x, normals, centers)
	for ii = 1:length(normals)
		out[ii] = dot(normals[ii], x) - dot(normals[ii], centers[ii])
	end
	out
end

"""
$(TYPEDEF)

This composite type (named for SectionPoincaréShooting) encodes a type of Poincaré sections implemented by hyperplanes. It can be used in conjunction with [`PoincareShootingProblem`](@ref). Each hyperplane is defined par a point (one example in `centers`) and a normal (one example in `normals`).

$(TYPEDFIELDS)

# Constructor(s)
    SectionPS(normals::Vector{Tv}, centers::Vector{Tv})

"""
struct SectionPS{Tn, Tc, Tnb, Tcb} <: AbstractSection
	M::Int64				# number of hyperplanes
	normals::Tn 			# normals to define hyperplanes
	centers::Tc 			# representative point on each hyperplane
	indices::Vector{Int64} 	# indices to be removed in the operator Ek

	normals_bar::Tnb
	centers_bar::Tcb

	function SectionPS(normals, centers)
		@assert length(normals) == length(centers)
		M = length(normals)
		indices = zeros(Int64, M)
		for ii in 1:M
			indices[ii] = argmax(abs.(normals[ii]))
		end
		nbar = [R(normals[ii], indices[ii]) for ii=1:M]
		cbar = [R(centers[ii], indices[ii]) for ii=1:M]

		return new{typeof(normals), typeof(centers), typeof(nbar), typeof(cbar)}(M, normals, centers, indices, nbar, cbar)
	end

	SectionPS(M = 0) = new{Nothing, Nothing, Nothing, Nothing}(M, nothing, nothing, Int64[], nothing, nothing)
end

(hyp::SectionPS)(out, u) = sectionHyp!(out, u, hyp.normals, hyp.centers)
isEmpty(sect::SectionPS{Tn, Tc, Tnb, Tcb}) where {Tn, Tc, Tnb, Tcb} = (Tn == Nothing) || (Tn == Nothing)

"""
	update!(hyp::SectionPS, normals, centers)

Update the hyperplanes saved in `hyp`.
"""
function update!(hyp::SectionPS, normals, centers)
	M = hyp.M
	@assert length(normals) == M "Wrong number of normals"
	@assert length(centers) == M "Wrong number of centers"
	for ii in 1:M
		hyp.normals[ii] .= normals[ii]
		hyp.centers[ii] .= centers[ii]
		hyp.indices[ii] = argmax(abs.(normals[ii]))
		R!(hyp.normals_bar[ii], normals[ii], hyp.indices[ii])
		R!(hyp.centers_bar[ii], centers[ii], hyp.indices[ii])
	end
	return hyp
end

# Operateur Rk from the paper above
@views function R!(out, x::AbstractVector, k::Int)
	out[1:k-1] .= x[1:k-1]
	out[k:end] .= x[k+1:end]
	return out
end

R!(hyp::SectionPS, out, x::AbstractVector, k::Int) = R!(out, x, hyp.indices[k])
R(x::AbstractVector, k::Int) = R!(similar(x, length(x) - 1), x, k)
R(hyp::SectionPS, x::AbstractVector, k::Int) = R!(hyp, similar(x, length(x) - 1), x, k)

# differential of R
dR!(hyp::SectionPS, out, dx::AbstractVector, k::Int) = R!(hyp, out, dx, k)

# Operateur Ek from the paper above
function E!(hyp::SectionPS, out, xbar::AbstractVector, ii::Int)
	@assert length(xbar) == length(hyp.normals[1]) - 1 "Wrong size for the projector / expansion operators, length(xbar) = $(length(xbar)) and length(normal) = $(length(hyp.normals[1]))"
	k = hyp.indices[ii]
	nbar = hyp.normals_bar[ii]
	xcbar = hyp.centers_bar[ii]
	coord_k = hyp.centers[ii][k] - (dot(nbar, xbar) - dot(nbar, xcbar)) / hyp.normals[ii][k]

	@views out[1:k-1] .= xbar[1:k-1]
	@views out[k+1:end] .= xbar[k:end]
	out[k] = coord_k
	return out
end

function E(hyp::SectionPS, xbar::AbstractVector, ii::Int)
	out = similar(xbar, length(xbar) + 1)
	E!(hyp, out, xbar, ii)
end

# differential of E!
function dE!(hyp::SectionPS, out, dxbar::AbstractVector, ii::Int)
	k = hyp.indices[ii]
	nbar = hyp.normals_bar[ii]
	xcbar = hyp.centers_bar[ii]
	coord_k = - dot(nbar, dxbar) / hyp.normals[ii][k]

	@views out[1:k-1]   .= dxbar[1:k-1]
	@views out[k+1:end] .= dxbar[k:end]
	out[k] = coord_k
	return out
end

function dE(hyp::SectionPS, dxbar::AbstractVector, ii::Int)
	out = similar(dxbar, length(dxbar) + 1)
	dE!(hyp, out, dxbar, ii)
end