text:
include("../src/JobSchedulers.jl")

using .JobSchedulers
using Base.Threads
using Dates
using Test

scheduler_start()
scheduler_status()

scheduler_stop()
scheduler_status()

scheduler_start()
sleep(2)
scheduler_status()

job = Job(@task(begin; sleep(2); println("highpriority"); end), name="high_priority", priority = 0)
display(job)
submit!(job)

j2 = job_query_by_id(job.id)
@test j2 === job

job2 = Job(@task(begin; sleep(2); println("lowpriority"); end), name="low_priority", priority = 20)
submit!(job2)
job = Job(@task(begin; sleep(2); println("highpriority"); end), name="high_priority", priority = 0)
submit!(job)
job = Job(@task(begin; sleep(2); println("midpriority"); end), name="mid_priority", priority = 15)
submit!(job)
for i in 1:20
    local job = Job(@task(begin; sleep(2); println(i); end), name="$i", priority = 20)
    submit!(job)
end

jobx = Job(@task(begin; sleep(20); println("run_success"); end), name="to_cancel", priority = 20)
submit!(jobx)
cancel!(jobx)

job2 = Job(@task(begin
    while true
        println(now())
        sleep(1)
    end
end), name="to_cancel", priority = 20)
submit!(job2)
cancel!(job2)

@test_logs (:error,) submit!(job2) # cannot resubmit
@test_throws Exception submit!(job) # cannot resubmit
@test_logs (:error,) submit!(job2)

### dependency
dep1 = Job(@task(begin
    sleep(2)
    println("dep1 ok")
end), name="dep1", priority = 20)

dep2 = Job(@task(begin
    sleep(3)
    println("dep2 ok")
end), name="dep2", priority = 20)

job_with_dep = Job(@task(begin
    println("job with dep1 and dep2 ok")
end), name="job_with_dep", priority = 20,
dependency = [DONE => dep1.id, DONE => dep2.id])

job_with_dep2 = Job(@task(begin
    println("job with dep2 ok")
end), name="job_with_dep", priority = 20,
dependency = DONE => dep2.id)

submit!(dep1)
submit!(dep2)
submit!(job_with_dep)
submit!(job_with_dep2)

while JobSchedulers.DataFrames.nrow(queue()) > 0
    sleep(2)
end

### set backup
rm("/tmp/jl_job_scheduler_backup", force=true)
rm("/tmp/jl_job_scheduler_backup2", force=true)
set_scheduler_backup("/tmp/jl_job_scheduler_backup")

set_scheduler_backup("/tmp/jl_job_scheduler_backup", migrate=true) # do nothing because file not exist

backup()
njobs = JobSchedulers.JOB_QUEUE_OK |> length

deleteat!(JobSchedulers.JOB_QUEUE_OK, 1:3:njobs)

set_scheduler_max_cpu(2)
set_scheduler_backup("/tmp/jl_job_scheduler_backup")
@test njobs == JobSchedulers.JOB_QUEUE_OK |> length
@test JobSchedulers.SCHEDULER_MAX_CPU == (Base.Threads.nthreads() > 1 ? Base.Threads.nthreads()-1 : Sys.CPU_THREADS)

job_queue_backup = deepcopy(JobSchedulers.JOB_QUEUE_OK)

deleteat!(JobSchedulers.JOB_QUEUE_OK, 1:3:njobs)
set_scheduler_backup("/tmp/jl_job_scheduler_backup")

@test njobs == JobSchedulers.JOB_QUEUE_OK |> length

set_scheduler_backup("/tmp/jl_job_scheduler_backup2", migrate=true)
backup()

deleteat!(JobSchedulers.JOB_QUEUE_OK, 2:3:njobs)
set_scheduler_backup("/tmp/jl_job_scheduler_backup")
@test njobs == JobSchedulers.JOB_QUEUE_OK |> length
set_scheduler_backup("/tmp/jl_job_scheduler_backup2", migrate=true, delete_old=true)

@test !isfile("/tmp/jl_job_scheduler_backup")
@test isfile("/tmp/jl_job_scheduler_backup2")

set_scheduler_backup("", delete_old=true)
@test !isfile("/tmp/jl_job_scheduler_backup2")

### Compat Pipelines.jl
using Pipelines
echo = CmdProgram(
    inputs = ["INPUT1", "INPUT2"],
    cmd = `echo INPUT1 INPUT2`
)
inputs = Dict(
    "INPUT1" => "Hello,",
    "INPUT2" => `Pipeline.jl`
)
cmdprog_job = Job(echo, inputs, touch_run_id_file=false)
cmdprog_job2 = Job(echo, inputs=inputs, touch_run_id_file=false)
cmdprog_job3 = Job(echo, touch_run_id_file=false)

submit!(cmdprog_job)
submit!(cmdprog_job2)
submit!(cmdprog_job3)

sleep(8)
@test cmdprog_job3.state == :failed

if Base.Threads.nthreads() > 1
    include("test_thread_id.jl")
end

@info "Test finished."


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
using JuMP
# rosenbrock

let

    m = Model()

    @variable(m, x)
    @variable(m, y)

    @NLobjective(m, Min, (1-x)^2 + 100(y-x^2)^2)

    solve(m)

    println("x = ", getvalue(x), " y = ", getvalue(y))

end

ext: jl


text:
using ExploringOceanData, GLMakie
@analysisplate

local_directory = "/mnt/podaac_drive/Version4/Release4/interp_monthly/"
# Salinity Stuff
salinity_extrema = []
# Variables of intersest
voi = ["SALT"]
variable = voi[1]
# Find available years
years = readdir(local_directory*variable)
year = years[13]
# Find available months
println("Looking at year ", year)
months = readdir(joinpath(local_directory*variable,year))
month = months[5]
tic = time()
println("Looking at month ", month)
full_path = joinpath(local_directory, variable, year, month)
# Load Data Set
ds = Dataset(full_path,"r")
toc = time()
println("Time to load dataset = ", toc-tic , " seconds")

salinity = ds["SALT"]

# Plot variables
latitude = ds["latitude"]
longitude = ds["longitude"]

skip = 1

var = salinity[begin:skip:end,begin:skip:end,1,1]

# lat-lon
lat = latitude[begin:skip:end, 1]
lon = longitude[begin:skip:end, 1]

missing_value = 0.0
bools = (var .== missing_value) .| (var .=== NaN) .| (var .=== missing)

var_max = quantile(abs.(var[.!bools]), 0.99)
clims = quantile.(Ref(var[.!bools]), [0.01,0.99]) # this is what is necessary for statistics

# for plotting make it a NaN
var[bools] .= NaN

# plt = heatmap(lat, lon, var, colormap = :balance, colorrange = clims)
plt = heatmap(var, colormap = :balance, colorrange = clims)
plt.axis.xlabel = "longitude"
plt.axis.ylabel = "latitude"
Colorbar(plt.figure[1, 2], plt.plot)
println("The quantiles are ", clims)

display(plt)


text:
function draw_neuron_panel(gui_data::GuiData, app_data::Model.AppData, sock::Comm.SocClient)
    if CImGui.TreeNode("Neuron")
        CImGui.PushItemWidth(80)

        # Row 1 *****************************************************
        gui_data.buffer = Model.prep_field(Model.refractory_period(app_data.model), 16)
        returned = CImGui.InputText("RefractoryPeriod", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_refractory_period!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(250)

        gui_data.buffer = Model.prep_field(Model.ap_max(app_data.model), 16)
        returned = CImGui.InputText("APMax", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_ap_max!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(450)

        gui_data.buffer = Model.prep_field(Model.threshold(app_data.model), 16)
        returned = CImGui.InputText("Threshold", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_threshold!(app_data.model, gui_data.buffer)
        end

        # Row 2 *****************************************************
        gui_data.buffer = Model.prep_field(Model.fast_surge(app_data.model), 16)
        returned = CImGui.InputText("Fast Surge", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_fast_surge!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(250)

        gui_data.buffer = Model.prep_field(Model.slow_surge(app_data.model), 16)
        returned = CImGui.InputText("Slow Surge", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_slow_surge!(app_data.model, gui_data.buffer)
        end

        # Row 3 *****************************************************
        gui_data.buffer = Model.prep_field(Model.tao(app_data.model), 16)
        returned = CImGui.InputText("Tao", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_tao!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(250)

        gui_data.buffer = Model.prep_field(Model.tao_j(app_data.model), 16)
        returned = CImGui.InputText("Tao J", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_tao_j!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(450)

        gui_data.buffer = Model.prep_field(Model.tao_s(app_data.model), 16)
        returned = CImGui.InputText("Tao S", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_tao_s!(app_data.model, gui_data.buffer)
        end

        # Row 4 *****************************************************
        gui_data.buffer = Model.prep_field(Model.weight_divisor(app_data.model), 16)
        returned = CImGui.InputText("Weight divisor", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_weight_divisor!(app_data.model, gui_data.buffer)
        end
        CImGui.SameLine(250)

        gui_data.buffer = Model.prep_field(Model.weight_max(app_data.model), 16)
        returned = CImGui.InputText("Weight max", gui_data.buffer, 16, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_weight_max!(app_data.model, gui_data.buffer)
        end

        gui_data.buffer = Model.prep_field(Model.active_synapse(app_data.model), 8)
        returned = CImGui.InputText("Active Synapse", gui_data.buffer, 8, CImGui.ImGuiInputTextFlags_CharsDecimal | CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
        if returned
            Model.set_active_synapse!(app_data.model, gui_data.buffer)
        end

        CImGui.PopItemWidth()

        CImGui.TreePop()
    end
end


text:
using Reexport

@reexport module UncertainDataPlotRecipes
    import ..Resampling.resample
    using RecipesBase

    # Recipes for uncertain values
    include("recipes_certainvalues.jl")
    include("recipes_populations.jl")
    include("recipes_uncertainvalues_theoretical.jl")
    include("recipes_uncertainvalues_kde.jl")
    include("recipes_uncertaindatasets.jl")

    # Uncertainties for values with uncertainties both in index and value
    include("recipes_uncertain_index_and_value.jl")
    include("recipes_uncertainindexvaluedataset.jl")
    include("recipes_vectors_of_uncertainvalues.jl")

    # Recipes for resampled statistics
    #include("recipes_resampledstatistics.jl")

end # module

"""
    UncertainDataPlotRecipes

Plot recipes for uncertain values and uncertain datasets.
"""
UncertainDataPlotRecipes


text:
import Base.complex

export Pauli, Id, X, Y, Z, allpaulis, paulieye, weight, complex

# Paulis's are represented by a vector of numbers (0-3) corresponding to
# single-qubit Paulis, along with a phase parameter.
immutable Pauli
    v::Vector{UInt8} # 0 = I, 1 = X, 2 = Z, 3 = Y
    s::UInt8 # 0 = +1, 1 = +i, 2 = -1, 3 = -i (im^s)
    function Pauli(v::Vector{UInt8}, s::UInt8)
        new(v,mod(s,4))
    end
end

Pauli(v::Vector, s = 0) = Pauli(convert(Vector{UInt8}, v), convert(UInt8, s))
Pauli(v::Integer, s = 0) = Pauli([v], s)
Pauli(m::Matrix) = convert(Pauli, m)

weight(p::Pauli) = sum( p.v .> 0 )

show(io::IO, p::Pauli) = print(io,convert(String,p))

==(a::Pauli, b::Pauli) = (a.v == b.v && a.s == b.s)
isequal(a::Pauli, b::Pauli) = (a == b)
hash(a::Pauli, h::UInt) = hash(a.v, hash(a.s, h))
isid(a::Pauli) = isempty(find(a.v))

function isless(a::Pauli, b::Pauli)
    # canonical total order defined by weight and then "lexicographic":
    # Id < X < Y < Z
    if weight(a) != weight(b)
        return weight(a) < weight(b)
    else
        return lex_tuple(a) < lex_tuple(b)
    end
end

function lex_tuple(p::Pauli)
    const pauli_lex = (1, 2, 4, 3)
    tuple([pauli_lex[x+1] for x in p.v]...)
end

function convert(::Type{String}, p::Pauli)
    phases = ["+","i","-","-i"]
    paulis = "IXZY"
    phases[p.s+1] * join([paulis[i+1] for i in p.v])
end

function convert{T}(::Type{Matrix{Complex{T}}}, p::Pauli)
    const mats = @compat Dict(
        0x00 => eye(Complex{T},2),
        0x01 => Complex{T}[0 1; 1 0],
        0x02 => Complex{T}[1 0; 0 -1],
        0x03 => Complex{T}[0 -im; im 0])

    return phase(p)*reduce(kron,[mats[x] for x in p.v])
end

complex(p::Pauli) = convert(Matrix{Complex128},p)

function convert(::Type{Pauli}, m::Matrix)
    d = size(m,1)
    n = log(2,size(m,1))
    for p in allpaulis(n)
        overlap = trace(m*convert(typeof(complex(m)),p)) / d
        if isapprox(abs(overlap),1,atol=d*eps(Float64))
            return (round(real(overlap))+im*round(imag(overlap)))*p
        elseif !isapprox(abs(overlap),0,atol=d*eps(Float64))
            println(m, overlap,isapprox(abs(overlap),0))
            error("Trying to convert non-Pauli matrix to a Pauli object")
        end
    end
end

promote_rule{T<:Real}(::Type{Pauli}, ::Type{Matrix{T}}) = Matrix{Complex{T}}
promote_rule{T<:Complex}(::Type{Pauli}, ::Type{Matrix{T}}) = Matrix{T}

function levicivita(a::UInt8, b::UInt8)
    # an unusual Levi-Civita pseudo-tensor for the (1,2,3) = (X,Z,Y) convention
    if (a,b) == (1,3) || (a,b) == (3,2) || (a,b) == (2,1)
        0x01
    elseif (a,b) == (1,2) || (a,b) == (2,3) || (a,b) == (3,1)
        0x03
    else
        0x00
    end
end
levicivita(x::@compat Tuple{UInt8,UInt8}) = levicivita(x...)
levicivita(a::Vector{UInt8}, b::Vector{UInt8}) = mapreduce(levicivita, +, zip(a,b))

# with our Pauli representation, multiplication is the sum (mod 4), or equivalently, the
# XOR of the bits
*(a::Pauli, b::Pauli) = Pauli(a.v $ b.v, a.s + b.s + levicivita(a.v, b.v))

const phases_ = [1, im, -1, -im]

function *(n::Number, p::Pauli)
    ns = findfirst(n .== phases_) - 1
    @assert(ns >= 0, "Multiplication only supported for +/- 1, +/- im")
    Pauli(p.v, p.s + ns)
end
*(p::Pauli, n::Number) = n * p
*(p::Pauli, u::Matrix) = *(promote(p, u)...)
*(u::Matrix, p::Pauli) = *(promote(u, p)...)

+(p::Pauli) = p
-(p::Pauli) = Pauli(p.v, p.s + 2)

abs(p::Pauli) = Pauli(p.v, 0)
phase(p::Pauli) = phases_[p.s+1]

length(p::Pauli) = length(p.v)
vec(p::Pauli) = vec(convert(Matrix{Complex{Int}}, p))

kron(a::Pauli, b::Pauli) = Pauli([a.v; b.v], a.s + b.s)

function expand(a::Pauli, subIndices, n)
    v = zeros(n)
    for (ct, i) in enumerate(subIndices)
        v[i] = a.v[ct]
    end
    Pauli(v, a.s)
end

function generators(a::Pauli)
    G = Pauli[]
    all(a.v .== 0) && return abs(a)
    s = phase(a)
    for (idx, p) in enumerate(a.v)
        if p == 0 # I, skip it
            continue
        elseif p == 1 || p == 2 # X or Z
            push!(G, expand(Pauli(p), [idx], length(a.v)))
        else # Y
            s *= im
            push!(G, expand(X, [idx], length(a.v)))
            push!(G, expand(Z, [idx], length(a.v)))
        end
    end
    G[1] *= s # give phase to first generator
    return G
end

function paulieye(n)
    expand(Id, [1], n)
end

# 1-qubit Paulis
const Id = Pauli(0)
const X = Pauli(1)
const Y = Pauli(3)
const Z = Pauli(2)

# 2-qubit Paulis
labelOpPairs = [("I", Id), ("X", X), ("Y", Y), ("Z", Z)]
for (a, ao) in labelOpPairs, (b, bo) in labelOpPairs
    @eval const $(Symbol(a*b)) = kron($ao, $bo)
end

function allpaulis(n)
    if n <= 0
        error("Need at least 1 qubit")
    elseif n==1
        return [Id,X,Y,Z]
    else
        return map(x->kron(x[1],x[2]),product([Id,X,Y,Z],allpaulis(n-1)))
    end
end


text:
using ReachabilityAnalysis
using LazySets.Arrays

# functions for Kronecker powers and sums
using ReachabilityAnalysis: kron_pow, kron_pow_stack

using LaTeXStrings

__toL(x, digits) = L"%$(round(x, digits=digits))"
xticklatex(vec, digits) = (vec, __toL.(vec, Ref(digits)));

# functions for error bounds
using CarlemanLinearization
using CarlemanLinearization: kron_sum,
                             error_bound_pseries,
                             error_bound_specabs,
                             convergence_radius_specabs,
                             _error_bound_specabs_R,
                             build_matrix

using ReachabilityAnalysis: ReachSolution

# TODO remove after ReachabilityAnalysis#551
LazySets.Approximations.box_approximation(X::Vector{IntervalArithmetic.Interval{Float64}}) = box_approximation(IntervalBox(X))

function _template(; n, N)
    dirs = Vector{SingleEntryVector{Float64}}()
    d = sum(n^i for i in 1:N)

    for i in 1:n
        x = SingleEntryVector(i, d, 1.0)
        push!(dirs, x)
    end
    for i in 1:n
        x = SingleEntryVector(i, d, -1.0)
        push!(dirs, x)
    end
    return CustomDirections(dirs)
end

function _project(sol, vars)
    πsol_1n = Flowpipe([ReachSet(set(project(R, vars)), tspan(R)) for R in sol])
end

# general method given a reachability algorithm for the linear system
function _solve_CARLIN_alg(X0, F1, F2, alg;
                           resets=0, N, T, Δt0=interval(0),
                           error_bound_func=error_bound_specabs)
    if resets == 0
        _solve_CARLIN(X0, F1, F2; alg=alg, N=N, T=T, Δt0=interval(0), error_bound_func=error_bound_func)
    else
        _solve_CARLIN_resets(X0, F1, F2; resets=resets, alg=alg, N=N, T=T, Δt0=interval(0), error_bound_func=error_bound_func)
    end
end

# method using the LGG09 reachability algorithm
function _solve_CARLIN_LGG09(X0, F1, F2;
                             resets=0, N, T, δ, Δt0=interval(0),
                             error_bound_func=error_bound_specabs)
    n = dim(X0)
    dirs = _template(n=n, N=N)
    alg = LGG09(δ=δ, template=dirs)
    return _solve_CARLIN_alg(X0, F1, F2, alg;
                             resets=resets, N=N, T=T, Δt0=Δt0,
                             error_bound_func=error_bound_func)
end

function _solve_CARLIN(X0, F1, F2; alg, N, T, Δt0=interval(0), A=nothing,
                                   error_bound_func=error_bound_specabs,
                                   bloat=true, kron_pow_algorithm="explicit")

    # lift initial states
    n = dim(X0)
    ŷ0 = kron_pow_stack(X0, N) |> box_approximation

    # solve continuous ODE
    if isnothing(A)
        A = build_matrix(F1, F2, N)
    end
    prob = @ivp(ŷ' = Aŷ, ŷ(0) ∈ ŷ0)
    sol = solve(prob, T=T, alg=alg, Δt0=Δt0)

    # projection onto the first n variables
    πsol_1n = _project(sol, 1:n)

    if !bloat
        return πsol_1n
    end

    # compute errors
    errfunc = error_bound_func(X0, Matrix(F1), Matrix(F2), N=N)

    # evaluate error bounds for each reach-set in the solution
    E = [errfunc.(tspan(R)) for R in sol]

    # if the interval is always > 0 then we can just take max(Ei)

    # symmetrize intervals
    E_rad = [symmetric_interval_hull(Interval(ei)) for ei in E]
    E_ball = [BallInf(zeros(n), max(ei)) for ei in E_rad]

    # sum the solution with the error
    fp_bloated = [ReachSet(set(Ri) ⊕ Ei, tspan(Ri)) for (Ri, Ei) in zip(πsol_1n, E_ball)] |> Flowpipe

    return ReachSolution(fp_bloated, alg)
end

function _compute_resets(resets::Int, T)
    return mince(0 .. T, resets+1)
end

function _compute_resets(resets::Vector{Float64}, T)
    # assumes initial time is 0
    aux = vcat(0.0, resets, T)
    return [interval(aux[i], aux[i+1]) for i in 1:length(aux)-1]
end

function _solve_CARLIN_resets(X0, F1, F2; resets, alg, N, T, Δt0=interval(0),
                                          error_bound_func=error_bound_specabs,
                                          bloat=true)

    # build state matrix (remains unchanged upon resets)
    A = build_matrix(F1, F2, N)

    # time intervals to compute
    time_intervals = _compute_resets(resets, T)

    # compute until first chunk
    T1 = sup(first(time_intervals))
    sol_1 = _solve_CARLIN(X0, F1, F2; alg=alg, N=N, T=T1, Δt0=interval(0), A=A, error_bound_func=error_bound_func, bloat=bloat)

    # preallocate output flowpipe
    fp_1 = flowpipe(sol_1)
    out = Vector{typeof(fp_1)}()
    push!(out, fp_1)

    # approximate final reach-set
    Rlast = sol_1[end]
    X0 = box_approximation(set(Rlast))

    # compute remaining chunks
    for i in 2:length(time_intervals)
        T0 = T1
        Ti = sup(time_intervals[i])
        sol_i = _solve_CARLIN(X0, F1, F2; alg=alg, N=N, T=Ti-T0, Δt0=interval(T0), A=A, error_bound_func=error_bound_func, bloat=bloat)
        push!(out, flowpipe(sol_i))
        X0 = box_approximation(set(sol_i[end]))
        T1 = Ti
    end

    return ReachSolution(MixedFlowpipe(out), alg)
end


text:
# ENV["JULIA_DEBUG"] = "all"

Base.push!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
import VSCodeDebugger
pop!(LOAD_PATH)

VSCodeDebugger.startdebugger()

max_stars_repo_name: fredrikekre/julia-vscode

max_issues_repo_name: fredrikekre/julia-vscode

max_forks_repo_name: fredrikekre/julia-vscode


text:
@testset "Testing Elements" begin
    elements = [
        Bool,
        #Float16,
        Float32,
        Float64,
        Int8,
        Int16,
        Int32,
        Int64,
        UInt8,
        UInt16,
        UInt32,
        UInt64,
    ]

    # Does this survive the round trip?
    for element in elements
        @test nGraph.back(nGraph.Element(element)) == element
    end
end

@testset "Testing Nodes" begin
    # Create a dummy parameter.
    param = nGraph.parameter(Float32, (10, 10))

    @test isa(param, nGraph.Node)

    # Test printing
    println(devnull, param)

    @test ndims(param) == 2
    @test size(param) == (10,10)
    @test size(param, 1) == 10
    @test size(param, 2) == 10
    @test eltype(param) == Float32
    @test eltype(param.obj) == Float32
    @test nGraph.description(param) == "Parameter"
    println(devnull, nGraph.name(param))

    # Go into higher dimensions.
    x = nGraph.parameter(Int32, (10, 20, 30))
    @test size(x) == (10, 20, 30)
    @test size(x, 3) == 30
end

@testset "Testing Function" begin
    # For this test, we need the `+` funcion defined in `ops`.
    a = nGraph.parameter(Float32, (10,10))
    b = nGraph.parameter(Float32, (10,10))

    z = a + b
    @test isa(z, nGraph.Node)
    @test eltype(z) == Float32
    @test size(z) == (10,10)
    @test nGraph.description(z) == "Add"

    fn = nGraph.NGFunction([a,b], [z])

    println(devnull, nGraph.name(fn))

    # At this point, we haven't compiled the function, so the poolsize should just default
    # to zero.
    @test iszero(nGraph.poolsize(fn))
end

@testset "Testing Backend" begin
    backend = nGraph.Backend{:CPU}()
    @test nGraph.version(backend) == "0.0.0"
end

@testset "Testing Tensor View" begin
    backend = nGraph.Backend{:CPU}()
    x = Array{Float32}(undef, 10, 10)

    tv = nGraph.Tensor(backend, x)
    @test parent(tv) === x
    @test eltype(tv) == eltype(parent(tv))
    @test sizeof(tv) == sizeof(parent(tv))
    @test size(tv) == size(parent(tv))
    println(devnull, tv)

    # Construct a Tensor from a scalar.
    tv = nGraph.Tensor(backend, 1.0)
    @test size(tv) == ()
    @test eltype(tv) == Float64

    # Construct a Tensor from a Node
    node = nGraph.parameter(Int32, (20, 20))
    tv = nGraph.Tensor(backend, node)
    @test size(node) == size(parent(node)) == (20, 20)
    @test eltype(node) == eltype(parent(node)) == Int32
end

@testset "Testing Executable" begin
    backend = nGraph.Backend{:CPU}()
    a = rand(Float32, 10, 10)
    b = rand(Float32, 10, 10)

    na = nGraph.Node(a)
    nb = nGraph.Node(b)

    nz = na + nb

    ex = nGraph.compile(backend, [na, nb], [nz])

    # At this point, we should have an executable we can call.
    # Create Tensor of the inputs and outputs and verify that the addition operation
    # is the same in both Julia land and ngraph land.
    tva = nGraph.Tensor(backend, a)
    tvb = nGraph.Tensor(backend, b)
    tvz = nGraph.Tensor(backend, nz)

    ex([tva, tvb], [tvz])

    @test isapprox(parent(tvz), a + b)
end


text:
module SimpleSVMs

import StatsBase
using JuMP
using LinearAlgebra: ⋅

export SVM, build_svm, L1Penalty, L2Penalty

struct SVM{W, B}
    w::W
    b::B
    SVM(w::W, b::B) where {W, B} = new{W,B}(w,b)
end

function StatsBase.fit(::Type{<:SVM}, penalty_func, X, y, optimizer)
    (m, w, b) = build_svm(penalty_func, X, y, optimizer)
    optimize!(m)
    @assert termination_status(m) == JuMP.MOI.OPTIMAL
    return SVM(JuMP.value.(w), JuP.value.(b))
end

function build_svm(penalty, X, y, optimizer)
    (nobs, nfeat) = size(X)
    length(y) == nobs || throw(DimensionMismatch("X and y need the same number of observations"))
    m = JuMP.Model(optimizer)
    @variable(m, w[1:nfeat])
    @variable(m, b)
    @variable(m, l[1:nobs] ≥ 0)
    @constraint(m, hinge_loss[i=1:nobs], l[i] ≥ 1 - y[i] * (X[i,:]⋅w + b))
    pen_cons = restrict_weights(m, penalty, w, b)
    @objective(m, Min, sum(l))
    return (m, w, b, pen_cons)
end

struct L1Penalty{R}
    rhs::R
    L1Penalty(rhs::R) where {R} = new{R}(rhs)
end

function restrict_weights(m::JuMP.AbstractModel, l1::L1Penalty, w, b)
    nfeat = length(w)
    @variable(m, wm[1:nfeat] ≥ 0)
    @variable(m, wp[1:nfeat] ≥  0)
    @constraint(m, weq[i=1:nfeat], wp[i] - wm[i] == w[i])
    @variable(m, bm ≥ 0)
    @variable(m, bp ≥ 0)
    @constraint(m, beq, bp - bm == b)
    @constraint(m, sum(wm) + sum(wp) + bp + bm ≤ l1.rhs)
end

struct L2Penalty{R}
    rhs::R
    L2Penalty(rhs::R) where {R} = new{R}(rhs)
end

function restrict_weights(m::JuMP.AbstractModel, l2::L2Penalty, w, b)
    @constraint(m, w ⋅ w  + b^2 ≤ l2.rhs)
end

end # module

max_stars_repo_name: matbesancon/SimpleSVMs


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Getting started with JuMP

# This tutorial is aimed at providing a quick introduction to writing JuMP code.
# If you're new to Julia, you should start with [Getting started with Julia](@ref).

# ## What is JuMP?

# JuMP ("Julia for Mathematical Programming") is an open-source modeling
# language that is embedded in Julia. It allows users to users formulate various
# classes of optimization problems (linear, mixed-integer, quadratic, conic
# quadratic, semidefinite, and nonlinear) with easy-to-read code. These problems
# can then be solved using state-of-the-art open-source and commercial solvers.

# JuMP also makes advanced optimization techniques easily accessible from a
# high-level language.

# ## Installation

# JuMP is a package for Julia. From Julia, JuMP is installed by using the
# built-in package manager.

# ```julia
# import Pkg
# Pkg.add("JuMP")
# ```

# You also need to include a Julia package which provides an appropriate solver.
# One such solver is `GLPK.Optimizer`, which is provided by the
# [GLPK.jl package](https://github.com/JuliaOpt/GLPK.jl).
# ```julia
# import Pkg
# Pkg.add("GLPK")
# ```
# See [Installation Guide](@ref) for a list of other solvers you can use.

# ## An example

# Let's try to solve the following linear programming problem by using JuMP and
# GLPK. We will first look at the complete code to solve the problem and then go
# through it step by step.

# ```math
# \begin{aligned}
# & \min & 12x + 20y \\
# & \;\;\text{s.t.} & 6x + 8y \geq 100 \\
# & & 7x + 12y \geq 120 \\
# & & x \geq 0 \\
# & & y \in [0, 3] \\
# \end{aligned}
# ```

using JuMP
using GLPK
model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, 0 <= y <= 3)
@objective(model, Min, 12x + 20y)
@constraint(model, c1, 6x + 8y >= 100)
@constraint(model, c2, 7x + 12y >= 120)
print(model)
optimize!(model)
@show termination_status(model)
@show primal_status(model)
@show dual_status(model)
@show objective_value(model)
@show value(x)
@show value(y)
@show shadow_price(c1)
@show shadow_price(c2)
nothing #hide

# ## Step-by-step

# Once JuMP is installed, to use JuMP in your programs, we just need to write:

using JuMP

# We also need to include a Julia package which provides an appropriate solver.
# We want to use `GLPK.Optimizer` here which is provided by the `GLPK.jl`
# package.

using GLPK

# A model object is a container for variables, constraints, solver options, etc.
# Models are created with the [`Model`](@ref) function. The model can be created
# with an optimizer attached with default arguments by calling the constructor
# with the optimizer type, as follows:

model = Model(GLPK.Optimizer)

# Variables are modelled using [`@variable`](@ref):

@variable(model, x >= 0)

# They can have lower and upper bounds.

@variable(model, 0 <= y <= 30)

# The objective is set using [`@objective`](@ref):

@objective(model, Min, 12x + 20y)

# Constraints are modelled using [`@constraint`](@ref). Here `c1` and `c2` are
# the names of our constraint.

@constraint(model, c1, 6x + 8y >= 100)

#-

@constraint(model, c2, 7x + 12y >= 120)

#- Call `print` to display the model:

print(model)

# To solve the optimization problem, call the [`optimize!`] function.

optimize!(model)

# !!! info
#     The `!` after optimize is just part of the name. It's nothing special.
#     Julia has a convention that functions which mutate their arguments should
#     end in `!`. A common example is `push!`.

# Now let's see what information we can query about the solution.

# [`termination_status`](@ref) tells us why the solver stopped:

termination_status(model)

# In this case, the solver found an optimal solution. We should also check
# [`primal_status`](@ref) to see if the solver found a primal feasible point:

primal_status(model)

# and [`dual_status`](@ref) to see if the solver found a dual feasible point:

dual_status(model)

# Now we know that our solver found an optimal solution, and has a primal and a
# dual solution to query.

# Query the objective value using [`objective_value`](@ref):

objective_value(model)

# The primal solution using [`value`](@ref):

value(x)

#-

value(y)

# and the dual solution using [`shadow_price`](@ref):

shadow_price(c1)

#-

shadow_price(c2)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

Base.showerror(io::IO, e::BoundaryConditionError) = print(
    io, "BoundaryConditionError: ", e.msg
)


text:
const R = MLJLinearModels
const CI = get(ENV, "CI", "false") == "true"

DO_COMPARISONS = DO_COMPARISONS && !CI
DO_COMPARISONS && (using PyCall; using RCall)
SKLEARN_LM = nothing
PY_RND     = nothing
if DO_COMPARISONS
    SKLEARN_LM = pyimport("sklearn.linear_model")
    PY_RND     = pyimport("random")
    QUANTREG   = rimport("quantreg")
end

m(s, p=true) = println("\n== $s ==" * ifelse(p, "\n", ""))
mm(s) = println("\n > $s < \n")

nnz(θ) = sum(abs.(θ) .> 0)

"""Make portion s of vector 0."""
sparsify!(θ, s) = (θ .*= (rand(length(θ)) .< s))

"""Add outliers to portion s of vector."""
outlify(y, s) = (n = length(y); y .+ 20 * randn(n) .* (rand(n) .< s))

"""Generate continuous (X, y) with and without intercept."""
function generate_continuous(n, p; seed=0, sparse=1)
    Random.seed!(seed)
    X  = randn(n, p)
    X1 = R.augment_X(X, true)
    θ  = randn(p)
    θ1 = randn(p+1)
    sparse < 1 && begin
        sparsify!(θ, sparse)
        sparsify!(θ1, sparse)
    end
    y  = X*θ + 0.1 * randn(n)
    y1 = X1*θ1 + 0.1 * randn(n)
    return ((X, y, θ), (X1, y1, θ1))
end

"""Generate continuous X and binary y with and without intercept."""
function generate_binary(n, p; seed=0, sparse=1)
    Random.seed!(seed)
    X  = randn(n, p)
    X1 = R.augment_X(X, true)
    θ  = randn(p)
    θ1 = randn(p+1)
    sparse < 1 && begin
        sparsify!(θ, sparse)
        sparsify!(θ1, sparse)
    end
    y  = rand(n) .< R.σ.(X*θ)
    y  = y .* ones(Int, n) .- .!y .* ones(Int, n)
    y1 = rand(n) .< R.σ.(X1*θ1)
    y1 = y1 .* ones(Int, n) .- .!y1 .* ones(Int, n)
    return ((X, y, θ), (X1, y1, θ1))
end

"""Simple function to sample from a multinomial."""
function multi_rand(Mp)
    # Mp[i, :] sums to 1
    n, c = size(Mp)
    be   = reshape(rand(length(Mp)), n, c)
    y    = zeros(Int, n)
    @inbounds for i in eachindex(y)
        rp = 1.0
        for k in 1:c-1
            if (be[i, k] < Mp[i, k] / rp)
                y[i] = k
                break
            end
            rp -= Mp[i, k]
        end
    end
    y[y .== 0] .= c
    return y
end

"""Generate continuous X and multiclass y with and without intercept."""
function generate_multiclass(n, p, c; seed=0, sparse=1)
    Random.seed!(seed)
    X   = randn(n, p)
    X1  = R.augment_X(X, true)
    θ   = randn(p * c)
    θ1  = randn((p+1) * c)
    sparse < 1 && begin
        sparsify!(θ, sparse)
        sparsify!(θ1, sparse)
    end
    y   = zeros(Int, n)
    y1  = zeros(Int, n)

    P   = R.apply_X(X, θ, c)
    M   = exp.(P)
    Mn  = M ./ sum(M, dims=2)
    P1  = R.apply_X(X, θ1, c)
    M1  = exp.(P1)
    Mn1 = M1 ./ sum(M1, dims=2)

    y  = multi_rand(Mn)
    y1 = multi_rand(Mn1)
    return ((X, y, θ), (X1, y1, θ1))
end


text:
@testset "Problems" begin
  grid2D = RegularGrid{Float64}(100,100)
  grid3D = RegularGrid{Float64}(100,100,100)

  @testset "Estimation" begin
    problem = EstimationProblem(data3D, grid3D, :value)
    @test data(problem) == data3D
    @test domain(problem) == grid3D
    @test variables(problem) == Dict(:value => Float64)
    @test_throws AssertionError EstimationProblem(data3D, grid2D, :value)

    # show methods
    grid2D = RegularGrid{Float64}(100,200)
    problem = EstimationProblem(data2D, grid2D, :value)
    @test sprint(show, problem) == "2D EstimationProblem"
    @test sprint(show, MIME"text/plain"(), problem) == "2D EstimationProblem\n  data:      3×3 GeoDataFrame (x and y)\n  domain:    100×200 RegularGrid{Float64,2}\n  variables: value (Float64)"
  end

  @testset "Simulation" begin
    problem = SimulationProblem(data3D, grid3D, :value, 100)
    @test data(problem) == data3D
    @test domain(problem) == grid3D
    @test variables(problem) == Dict(:value => Float64)
    @test hasdata(problem)
    @test nreals(problem) == 100
    @test_throws AssertionError SimulationProblem(data3D, grid2D, :value, 100)

    # show methods
    grid2D = RegularGrid{Float64}(100,200)
    problem = SimulationProblem(data2D, grid2D, :value, 100)
    @test sprint(show, problem) == "2D SimulationProblem (conditional)"
    @test sprint(show, MIME"text/plain"(), problem) == "2D SimulationProblem (conditional)\n  data:      3×3 GeoDataFrame (x and y)\n  domain:    100×200 RegularGrid{Float64,2}\n  variables: value (Float64)\n  N° reals:  100"
  end
end


text:
using Postprocessing
using Test
using LinearAlgebra, DelimitedFiles

@testset "Fourier Padé reconstruction" begin
  println(devnull, reference(fourier_pade))

  for T in (Float32, Float64)
    data = readdlm(joinpath(@__DIR__, "test_u_piecewise_constant.txt"), comments=true)
    x = T.(data[:, 1])
    u = T.(data[:, 2])
    u_true = @. ifelse(abs(x - 0.6) < 0.2, one(T), zero(T))

    u_reconstructed = fourier_pade(u, 50, 50)
    @test norm(u_reconstructed - u_true) < norm(u - u_true)
    @test total_variation(u_reconstructed) < total_variation(u)

    u_reconstructed = fourier_pade(u, 50, 45)
    @test norm(u_reconstructed - u_true) < norm(u - u_true)
    @test total_variation(u_reconstructed) < total_variation(u)

    # this choice of the degrees leads to increased oscillations
    # near the discontinuity but should still be computable
    u_reconstructed = fourier_pade(u, 50, 52)

    N = length(u)
    @test_throws ArgumentError fourier_pade(u[1:end-1], N-(N÷2)+1, N÷2)
  end
end


text:
module SossGen

using Reexport

@reexport using Soss
@reexport using Gen
const logpdf = Soss.logpdf
const ifelse = Base.ifelse

import Gen

include("utils.jl")

# Following https://probcomp.github.io/Gen/dev/ref/extending/#Custom-generative-functions-1
# Define a trace data type

struct SossTrace{F, NT} <: Gen.Trace
    gen_fn :: F
    choices :: NT
    logprob :: Float64
    args
end

# Define a data type for your generative function

export SossGenerativeFunction
struct SossGenerativeFunction{M,F,V}  <: Gen.GenerativeFunction{V,SossTrace{SossGenerativeFunction}}
    model  :: M
    logpdf :: F
end

# Decide what the arguments to a generative function should be

function SossGenerativeFunction(m::Model{A,B,M}, ::typeof(codegen)) where {A,B,M}
    ℓ(args, data) = logpdf(m(args...), data, codegen)::Float64
    SossGenerativeFunction{Model{A,B,M},typeof(ℓ), typeof(variables(m))}(m, ℓ)
end

function SossGenerativeFunction(m::Model{A,B,M}) where {A,B,M}
    ℓ(args, data) = logpdf(m(args...), data)::Float64
    SossGenerativeFunction{Model{A,B,M},typeof(ℓ), typeof(variables(m))}(m, ℓ)
end

# Implement methods of the Generative Function Interface

# simulate

function Gen.simulate(gen_fn::SossGenerativeFunction{M,F,V}, args::Tuple) where {M,F,V}
    choices = rand(gen_fn.model(args...))
    score = gen_fn.logpdf(args, choices)
    SossTrace{SossGenerativeFunction,typeof(choices)}(gen_fn, choices, score, args)
end

# has_argument_grads

# accepts_output_grad

# get_args

function Gen.get_args(t::SossTrace{SossGenerativeFunction,NT}) where {NT}
    t.args
end

# get_retval

function Gen.get_retval(t::SossTrace{SossGenerativeFunction,NT}) where {NT}
    t.choices
end

# get_choices

function Gen.get_choices(t::SossTrace{SossGenerativeFunction,NT}) where {NT}
    c = Gen.choicemap()
    Ts = NTtypes(NT)
    for (k::Symbol, v) in pairs(t.choices)
        c[k] = v :: Ts[k]
    end
    c
end

# get_score

function Gen.get_score(t::SossTrace{SossGenerativeFunction,NT}) where {NT}
    t.logprob
end

# get_gen_fn

# project

Gen.project(t::SossTrace{SossGenerativeFunction,NT}, ::Gen.EmptySelection) where {NT} = 0.0

# generate

function Gen.generate(gen_fn::SossGenerativeFunction{M,F,V}, args::Tuple, constraints::Gen.ChoiceMap) where {M,F,V}
    kvs = Gen.get_values_shallow(constraints)
    data = namedtuple(Dict{Symbol, Any}(kvs))
    weight, choices = weightedSample(gen_fn.model(args...), data)
    logprob = gen_fn.logpdf(args, choices)
    SossTrace{SossGenerativeFunction,typeof(choices)}(gen_fn, choices, logprob, args), weight
end

# update

function Gen.update(t::SossTrace{SossGenerativeFunction,NT}, new_args::Tuple, argdiffs::Tuple, constraints::Gen.ChoiceMap) where {NT}
    retdiff = ifelse(isempty(constraints), Gen.NoChange(), Gen.UnknownChange())
    new_choices = NT(merge(t.choices, namedtuple(Dict{Symbol,Any}(Gen.get_values_shallow(constraints))))) # :: NT
    new_logprob = t.gen_fn.logpdf(new_args, new_choices) :: Float64
    new_trace = SossTrace{SossGenerativeFunction,NT}(t.gen_fn, new_choices, new_logprob, new_args)
    weight = new_logprob - t.logprob
    discard = Gen.choicemap()
    for (k::Symbol, v) in Gen.get_values_shallow(constraints)
        discard[k] = t.choices[k]
    end
    return (new_trace, weight, retdiff, discard)
end

# regenerate

# choice_gradients

# accumulate_param_gradients!

end

max_stars_repo_name: cscherrer/SossGen

max_issues_repo_name: cscherrer/SossGen

max_forks_repo_name: cscherrer/SossGen


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

mutable struct RecoveryWalkResponse <: SwaggerModel
    walkPerformed::Any # spec type: Union{ Nothing, Bool } # spec name: walkPerformed
    nextPlatformUpdateDomain::Any # spec type: Union{ Nothing, Int32 } # spec name: nextPlatformUpdateDomain

    function RecoveryWalkResponse(;walkPerformed=nothing, nextPlatformUpdateDomain=nothing)
        o = new()
        validate_property(RecoveryWalkResponse, Symbol("walkPerformed"), walkPerformed)
        setfield!(o, Symbol("walkPerformed"), walkPerformed)
        validate_property(RecoveryWalkResponse, Symbol("nextPlatformUpdateDomain"), nextPlatformUpdateDomain)
        setfield!(o, Symbol("nextPlatformUpdateDomain"), nextPlatformUpdateDomain)
        o
    end
end # type RecoveryWalkResponse

const _property_map_RecoveryWalkResponse = Dict{Symbol,Symbol}(Symbol("walkPerformed")=>Symbol("walkPerformed"), Symbol("nextPlatformUpdateDomain")=>Symbol("nextPlatformUpdateDomain"))
const _property_types_RecoveryWalkResponse = Dict{Symbol,String}(Symbol("walkPerformed")=>"Bool", Symbol("nextPlatformUpdateDomain")=>"Int32")
Base.propertynames(::Type{ RecoveryWalkResponse }) = collect(keys(_property_map_RecoveryWalkResponse))
Swagger.property_type(::Type{ RecoveryWalkResponse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RecoveryWalkResponse[name]))}
Swagger.field_name(::Type{ RecoveryWalkResponse }, property_name::Symbol) =  _property_map_RecoveryWalkResponse[property_name]

function check_required(o::RecoveryWalkResponse)
    true
end

function validate_property(::Type{ RecoveryWalkResponse }, name::Symbol, val)
end


text:
struct ClopperPearsonInterval{TO<:Ordering,TD<:AbstractDesign,TI<:Integer,TR<:Real} <: IntervalEstimator
    ordering::TO
    design::TD
    XX::Array{TI,2}
    bounds::Array{TR,2}
    alpha::TR
end

string(ci::ClopperPearsonInterval) = @sprintf "ClopperPearsonInterval<%s>" string(ci.ordering)

function ClopperPearsonInterval(ordering::Ordering, design::AbstractDesign, α::Real; epsilon = 1e-6)

    ϵ = epsilon
    XX = sample_space(design)
    nn = size(XX, 1)

    pval_sup(p0, x1, x2) = p_value(x1, x2, p0, ordering, design; orientation = :superiority)
    pval_inf(p0, x1, x2) = p_value(x1, x2, p0, ordering, design; orientation = :inferiority)

    bounds = zeros(Float64, nn, 2)
    for i in 1:nn
        x1, x2 = XX[i,:]
        l = try
                Roots.find_zero(p -> pval_sup(p, x1, x2) - (α + ϵ), (ϵ, 1 - ϵ))
            catch e
                isa(e, Roots.ConvergenceFailed) | isa(e, Roots.ArgumentError) ? 0.0 : e
        end
        u = try
                Roots.find_zero(p -> pval_inf(p, x1, x2) - (α + ϵ), (ϵ, 1 - ϵ))
            catch e
                isa(e, Roots.ConvergenceFailed) | isa(e, Roots.ArgumentError) ? 1.0 : e
        end
        u <= l ? error("confidence interval collapsed") : nothing
        bounds[i,:] .= [l, u]
    end
    return ClopperPearsonInterval{typeof(ordering),typeof(design),eltype(XX),Float64}(ordering, design, XX, bounds, α)
end

function ClopperPearsonInterval(estimator::Estimator, design::AbstractDesign, α::Real; ϵ = 1e-6)

    return ClopperPearsonInterval(EstimatorOrdering(estimator), design, α; ϵ = ϵ)
end


text:
using LinearAlgebra
using SparseArrays
using ProximalOperators
using Test

@testset "LeastSquares" begin

@testset "$(T), $(s), $(matrix_type), $(mode)" for (T, s, matrix_type, mode) in Iterators.product(
    [Float64, ComplexF64],
    [(10, 29), (29, 10), (10, 29, 3), (29, 10, 3)],
    [:dense, :sparse],
    [:direct, :iterative],
)

if mode == :iterative && length(s) == 3
    # FIXME this case is currently not supported due to cg! in IterativeSolvers
    # See https://github.com/JuliaMath/IterativeSolvers.jl/issues/248
    # The fix is simple, but affects a lot of solvers in IterativeSolvers
    # Maybe we can use our own CG here and drop the dependency
    continue
end

R = real(T)

shape_A = s[1:2]
shape_b = if length(s) == 2 s[1] else s[[1, 3]] end
shape_x = if length(s) == 2 s[2] else s[2:3] end

A = if matrix_type == :sparse sparse(randn(T, shape_A...)) else randn(T, shape_A...) end
b = randn(T, shape_b...)
x = randn(T, shape_x...)

f = LeastSquares(A, b, iterative=(mode == :iterative))
predicates_test(f)

@test ProximalOperators.is_smooth(f) == true
@test ProximalOperators.is_quadratic(f) == true
@test ProximalOperators.is_generalized_quadratic(f) == true
@test ProximalOperators.is_set(f) == false

grad_fx, fx = gradient(f, x)
lsres = A*x - b
@test fx ≈ 0.5*norm(lsres)^2
@test all(grad_fx .≈ (A'*lsres))

call_test(f, x)
prox_test(f, x)
prox_test(f, x, R(1.5))

lam = R(0.1) + rand(R)
f = LeastSquares(A, b, lam, iterative=(mode == :iterative))
predicates_test(f)

grad_fx, fx = gradient(f, x)
@test fx ≈ (lam/2)*norm(lsres)^2
@test all(grad_fx .≈ lam*(A'*lsres))

call_test(f, x)
prox_test(f, x)
prox_test(f, x, R(2.1))

end

end


text:
### RLEVectors Type

"""
# RLEVectors
`RLEVectors` is an alternate implementation of the Rle type from Bioconductor's
IRanges package by H. Pages, P. Aboyoun and M. Lawrence. RLEVectors represent a
vector with repeated values as the ordered set of values and repeat extents. In
the field of genomics, data of various types measured across the ~3 billion
letters in the human genome can often be represented in a few thousand runs.
It is useful to know the bounds of genome regions covered by these runs, the
values associated with these runs, and to be able to perform various
mathematical operations on these values.

`RLEVectors` can be created from a single vector or a vector of values and a
vector of run ends. In either case runs of values or zero length runs will
be compressed out. RLEVectors can be expanded to a full vector like a
`Range` with `collect`.

### Examples
 * `x = RLEVector([1,1,2,2,3,3,4,4,4])`
 * `x = RLEVector([4,5,6],[3,6,9])`

"""
## Types and constructors

type RLEVector{T1,T2 <: Integer} <: AbstractArray{T1, 1}
  runvalues::Vector{T1}
  runends::Vector{T2}
  function RLEVector(runvalues, runends)
    rle = new(runvalues,runends)
    return(rle)
  end
end

function RLEVector{T1,T2 <: Integer}(runvalues::Vector{T1}, runends::Vector{T2})
  nrun = numruns(runvalues,runends)
  if nrun != length(runends)
    runvalues, runends = ree(runvalues,runends,nrun)
  end
  RLEVector{T1,T2}(runvalues, runends)
end

function RLEVector{T2 <: Integer}(runvalues::BitVector, runends::Vector{T2})
  nrun = numruns(runvalues,runends)
  if nrun != length(runends)
    runvalues, runends = ree(runvalues,runends,nrun)
  end
  RLEVector{Bool,T2}(runvalues, runends)
end

function RLEVector(vec::Vector)
  runvalues, runends = ree(vec)
  RLEVector(runvalues, runends)
end

#  Having specific types of Rle would be useful for lists of the same type, but Julia does a good job noticing that
#  Could also be useful for method definitions
typealias FloatRle RLEVector{Float64,UInt32}
typealias IntegerRle RLEVector{Int64,UInt32}
typealias BoolRle RLEVector{Bool,UInt32}
typealias StringRle RLEVector{String,UInt32}

# similar
function similar(x::RLEVector, element_type::Type, dims::Dims)
    len = dims[1]
    if len == 0
        return( RLEVector(element_type[], eltype(x.runends)[]) )
    else
        return( RLEVector(zeros(element_type, 1), eltype(x.runends)[len]) )
    end
end

# show
Base.show(io::IO, ::MIME"text/plain",  a::RLEVector) = show(io, a)
function show(io::IO, x::RLEVector)
    t = typeof(x)::DataType
    show(io, t)
    print("\n")
    n = nrun(x)
    if n > 10
        rv = x.runvalues
        re = x.runends
        print("run values: [$(rv[1]),$(rv[2]),$(rv[5]) \u2026 $(rv[n-4]),$(rv[n-1]),$(rv[n])]\n")
        print("run ends:   [$(re[1]),$(re[2]),$(re[5]) \u2026 $(re[n-4]),$(re[n-1]),$(re[n])]")
    else
        println("run values: ", x.runvalues)
        println("run ends:   ", x.runends)
    end
end

# conversions
convert(::Type{Vector}, x::RLEVector) = collect(x)
convert(::Type{Set}, x::RLEVector) = Set(rvalue(x))
promote_rule(::Type{Set}, ::Type{RLEVector}) = Set

# the basics
function collect(x::RLEVector)
  inverse_ree(x.runvalues,x.runends)
end

function isequal(x::RLEVector, y::RLEVector)
  isequal(x.runends,y.runends) && isequal(x.runvalues, y.runvalues)
end

Base.hash(a::RLEVector) = hash(a.runvalues, hash(a.runlengths, hash(:RLEVector)))
==(a::RLEVector, b::RLEVector) = isequal(a.runvalues, b.runvalues) && isequal(a.runends, b.runends) && true