text:
@muladd function perform_step!(integrator,
                               cache::Union{ImplicitEMConstantCache,
                                            ImplicitEulerHeunConstantCache,
                                            ImplicitRKMilConstantCache},
                                            f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack uf = cache
  nlsolve! = cache.nlsolve; nlcache = nlsolve!.cache
  alg = unwrap_alg(integrator, true)
  theta = alg.theta
  alg.symplectic ? a = dt/2 : a = dt
  γdt = dt*theta
  repeat_step = false
  if nlsolve! isa NLNewton
    uf.t = t
    J, nlcache.W = calc_W!(integrator, cache, γdt, repeat_step)
  end

  # TODO: Stochastic extrapolants?
  u = uprev

  L = integrator.g(uprev,p,t)
  ftmp = integrator.f(uprev,p,t)
  gtmp = L.*integrator.W.dW

  if typeof(cache) <: ImplicitEulerHeunConstantCache
    utilde = uprev + gtmp
    gtmp = ((integrator.g(utilde,p,t) + L)/2)*integrator.W.dW
  end

  if typeof(cache) <: ImplicitRKMilConstantCache || integrator.opts.adaptive == true
    if alg_interpretation(alg) == :Ito ||
       typeof(cache) <: ImplicitEMConstantCache
      K = @. uprev + dt * ftmp
      utilde =  K + L*integrator.sqdt
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      mil_correction = ggprime .* (integrator.W.dW.^2 .- dt)./2
      gtmp += mil_correction
    elseif alg_interpretation(alg) == :Stratonovich ||
           typeof(cache) <: ImplicitEulerHeunConstantCache
      utilde = uprev + L*integrator.sqdt
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      mil_correction = ggprime.*(integrator.W.dW.^2)./2
      gtmp += mil_correction
    end
  end

  if alg.symplectic
    z = zero(u) # constant extrapolation, justified by ODE IM
  else
    z = dt*ftmp # linear extrapolation
  end
  nlcache.z = z

  nlcache.c = a
  if alg.symplectic
    # u = uprev + z then  u = (uprev+u)/2 = (uprev+uprev+z)/2 = uprev + z/2
    #u = uprev + z/2 + gtmp/2
    tmp = uprev + gtmp/2
  else
    #u = uprev + dt*(1-theta)*ftmp + theta*z + gtmp
    tmp = uprev + dt*(1-theta)*ftmp + gtmp
  end
  nlcache.tmp = tmp
  (z, η, iter, fail_convergence) = nlsolve!(integrator); fail_convergence && return nothing

  if alg.symplectic
    u = uprev + z + gtmp
  else
    #u = uprev + dt*(1-theta)*ftmp + theta*z + gtmp
    u = tmp + theta*z
  end

  nlcache.ηold = η
  nlcache.nl_iters = iter

  if integrator.opts.adaptive

    Ed = _reshape(dt*(J*_vec(ftmp))/2, axes(ftmp))
    if typeof(cache) <: Union{ImplicitEMConstantCache,ImplicitEulerHeunConstantCache}
        En = mil_correction
    else
        En = integrator.opts.internalnorm.(integrator.W.dW.^3) .*
             integrator.opts.internalnorm.(ggprime).^2 ./ 6
    end

    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(resids)
  end

  integrator.u = u
  return nothing
end

@muladd function perform_step!(integrator,
                               cache::Union{ImplicitEMCache,
                                            ImplicitEulerHeunCache,
                                            ImplicitRKMilCache},
                               f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack du1,dz,z,k,J,W,jac_config,gtmp,gtmp2,tmp = cache
  nlsolve! = cache.nlsolve; nlcache = nlsolve!.cache
  alg = unwrap_alg(integrator, true)
  alg.symplectic ? a = dt/2 : a = dt
  dW = integrator.W.dW
  mass_matrix = integrator.f.mass_matrix
  theta = alg.theta

  repeat_step = false

  if integrator.success_iter > 0 && !integrator.u_modified && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
  elseif alg.extrapolant == :linear
    @. u = uprev + integrator.fsalfirst*dt
  else # :constant
    copyto!(u,uprev)
  end

  nlsolve! isa NLNewton && calc_W!(integrator, cache, dt, repeat_step)

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

  # Handle noise computations

  integrator.g(gtmp,uprev,p,t)
  integrator.f(tmp,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @. gtmp2 = gtmp*dW
  else
    mul!(gtmp2,gtmp,dW)
  end

  if typeof(cache) <: ImplicitEulerHeunCache
    gtmp3 = cache.gtmp3
    @. z = uprev + gtmp2
    integrator.g(gtmp3,z,p,t)
    @. gtmp = (gtmp3 + gtmp)/2
    if is_diagonal_noise(integrator.sol.prob)
      @. gtmp2 = gtmp*dW
    else
      mul!(gtmp2,gtmp,dW)
    end
  end

  if typeof(cache) <: ImplicitRKMilCache
    gtmp3 = cache.gtmp3
    if alg_interpretation(alg) == :Ito
      @. z = uprev + dt * tmp + integrator.sqdt * gtmp
      integrator.g(gtmp3,z,p,t)
      @. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation
      @. gtmp2 += gtmp3*(dW.^2 - dt)/2
    elseif alg_interpretation(alg) == :Stratonovich
      @. z = uprev + integrator.sqdt * gtmp
      integrator.g(gtmp3,z,p,t)
      @. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation
      @. gtmp2 += gtmp3*(dW.^2)/2
    end
  end

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

  if alg.symplectic
    @. z = zero(eltype(u)) # Justified by ODE solvers, constrant extrapolation when IM
  else
    @. z = dt*tmp # linear extrapolation
  end

  nlcache.c = a
  if alg.symplectic
    #@. u = uprev + z/2 + gtmp2/2
    @. tmp = uprev + gtmp2/2
  else
    #@. u = uprev + dt*(1-theta)*tmp + theta*z + gtmp2
    @. tmp = uprev + dt*(1-theta)*tmp + gtmp2
  end
  z,η,iter,fail_convergence = nlsolve!(integrator)

  if alg.symplectic
    @. u = uprev + z + gtmp2
  else
    @. u = tmp + theta*z
  end

  nlcache.ηold = η
  nlcache.nl_iters = iter

  if integrator.opts.adaptive

    if has_invW(f)
      # This means the Jacobian was never computed!
      f.jac(J,uprev,p,t)
    end

    mul!(vec(z),J,vec(tmp))
    @. k = dt*dt*z/2

    # k is Ed
    # dz is En
    if typeof(cache) <: Union{ImplicitEMCache,ImplicitEulerHeunCache}
      dW_cache = cache.dW_cache
      if !is_diagonal_noise(integrator.sol.prob)
        g_sized = norm(gtmp,2)
      else
        g_sized = gtmp
      end

      if typeof(cache) <: ImplicitEMCache
        @. z = uprev + dt * tmp + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @. dW_cache = dW.^2 - dt
          diff_tmp = integrator.opts.internalnorm(dW_cache)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2 - dt)
        end

      elseif typeof(cache) <: ImplicitEulerHeunCache
        @. z = uprev + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @. dW_cache = dW.^2
          diff_tmp = integrator.opts.internalnorm(dW_cache)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2)
        end

      end

    elseif typeof(cache) <: ImplicitRKMilCache
      # gtmp3 is ggprime
      @. dz = integrator.opts.internalnorm(dW^3)*integrator.opts.internalnorm(gtmp3)^2 / 6
    end

    calculate_residuals!(tmp, k, dz, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(tmp)

  end
end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

"""
    register!(lks::Vector{Link}, lk::Link)

Register a link `lk` to a vector of links `lks`.
"""
register!(lks::Vector{Link}, lk::Link) = push!(lks, lk)
register!(lks::Vector{Channel}, lk::Link) = push!(lks, lk)

"""
	istaskfailed(lk::Link)

Returns true if a task associated with `lk` has failed.
"""
Base.istaskfailed(lk::Link) = !isnothing(lk.chn.excp)

"""
    istaskfailed(lks::Vector{Link})

Returns true if any task associated with a vector `lks` of links
has failed.
"""
Base.istaskfailed(lks::Vector{Link}) = any(lk->istaskfailed(lk), lks)

"""
	info(lk::Link)

Return the state (eventually the stacktrace) of a task associated
with `lk`.
"""
function info(lk::Link)
	if istaskfailed(lk)
		return hasfield(typeof(lk.chn.excp), :task) ? lk.chn.excp.task : lk.chn.excp
	else
		return lk.chn.cond_take.waitq.head.donenotify.waitq.head.code.task.state
	end
end

"""
```
diag!(lk::Link)
diag!(name::Symbol)
```
Return the internal `_ACT` variable of the `lk` actor.
This is only for diagnosis and testing.
"""
diag!(lk::Link) = request!(lk, Diag)
diag!(name::Symbol) = request!(name, Diag)

ext: jl


text:
using JuliaDB, IndexedTables, OnlineStats
using Dates

path =  joinpath(homedir(), "PacktPub","Chp09","Files","Stocks")
tickdata = loadndsparse(path, indexcols=["date", "ticker"])
tickdata[Date("2010-06-01"), "GOOGL"]
tickdata[Date("2012-01"):Dates.Month(1):Date("2014-12"), ["GOOGL", "KO"]]

filter(:date=>Dates.ismonday, tickdata)
filter((1=>Dates.ismonday, 2=>x->startswith(x, "G")), tickdata)
filter(x -> x.low > 10.0, tickdata)

diffs = map(x -> x.high - x.low, tickdata)
volumes = map(x -> x.volume, tickdata)

BigInt(reduce(+, map(x -> x.volume, tickdata)))

using MySQL, IndexedTables
using Tables

conn = MySQL.connect("localhost", "root", "LQaLxxy6",db="Chinook",
            unix_socket="/Applications/MAMP/tmp/mysql/mysql.sock");

sql = """select a.FirstName,a.LastName, sum(b.Total) as TotalSpent
 from Customer a
 join Invoice b on a.CustomerId = b.CustomerId
 group by a.FirstName, a.LastName""";

tbl = MySQL.query(conn,sql) |> table

rows = Tables.rows(tbl)
sch = Tables.schema(rows)
names = sch.names

for row in rows
    Tables.eachcolumn(sch, row) do val, col, name
        println(col," :: ",name)
    end
end


text:
function error(Mdim::Int, Ndim::Int, Pdim::Int, C::Array{T}) where T
    cval  = Float32(Pdim * AVAL * BVAL)
    errsq = 0f0
    for i in 1:Ndim
        for j in 1:Mdim
            err = C[(i-1)*Ndim+j] - cval
            if isnan(err)
                println((i,j))
            end
            errsq += err^2
        end
    end
    return errsq
end

function results(Mdim::Int, Ndim::Int, Pdim::Int, C::Array{T}, run_time) where T
    mflops = 2.0 * Mdim * Ndim * Pdim/(1000000.0* run_time)
    println("$run_time seconds at $mflops MFLOPS")
    errsq = error(Mdim, Ndim, Pdim, C)
    if isnan(errsq) || errsq > TOL
        println("Errors in multiplication: $errsq")
    end
end


text:
const categoricalplottypes = [BarPlot, Heatmap, Volume]

to_weights(v) = weights(v)
to_weights(v::AbstractWeights) = v

function compute_edges(extrema, bins::Tuple{Vararg{Integer}}, closed)
    return map(extrema, bins) do (min, max), n
        histrange(float(min), float(max), n, closed)
    end
end
compute_edges(extrema, bins::Tuple{Vararg{AbstractArray}}, closed) = bins

function midpoints(edges::AbstractRange)
    min, s, l = minimum(edges), step(edges), length(edges)
    return range(min + s / 2, step=s, length=l - 1)
end

function _histogram(data...; bins=sturges(length(data[1])), weights=automatic,
    normalization::Symbol, datalimits::Tuple, closed::Symbol)

    bins_tuple = bins isa Tuple ? bins : map(_ -> bins, data)
    edges = compute_edges(datalimits, bins_tuple, closed)
    w = weights === automatic ? () : (to_weights(weights),)
    h = fit(Histogram, data, w..., edges)
    return normalize(h, mode=normalization)
end

Base.@kwdef struct HistogramAnalysis{D, B}
    datalimits::D=automatic
    bins::B=automatic
    closed::Symbol=:left
    normalization::Symbol=:none
end

function (h::HistogramAnalysis)(input::ProcessedLayer)
    datalimits = compute_datalimits(input.positional, h.datalimits)
    options = valid_options(; datalimits, h.bins, h.closed, h.normalization)

    output = map(input) do p, n
        hist = _histogram(p...; pairs(n)..., pairs(options)...)
        return (map(midpoints, hist.edges)..., hist.weights), (;)
    end

    N = length(input.positional)
    label = h.normalization == :none ? "count" : string(h.normalization)
    labels = set(output.labels, N+1 => label)
    attributes = if N == 1
        set(output.attributes, :dodge_gap => 0, :x_gap => 0)
    else
        output.attributes
    end
    default_plottype = categoricalplottypes[N]
    plottype = Makie.plottype(output.plottype, default_plottype)
    return ProcessedLayer(output; plottype, labels, attributes)
end

"""
    histogram(; bins=automatic, datalimits=automatic, closed=:left, normalization=:none)

Compute a histogram. `bins` can be an `Int` to create that
number of equal-width bins over the range of `values`. In that case, the range covered
by the `bins` is defined by `datalimits` (defaults to the extrema of the data).
Alternatively, `bins` can be a sorted iterable of bin edges.
`closed` determines whether the the intervals are closed to the left or to the right.
The histogram can be normalized by setting `normalization`. Possible values are:
*  `:pdf`: Normalize by sum of weights and bin sizes. Resulting histogram
   has norm 1 and represents a PDF.
* `:density`: Normalize by bin sizes only. Resulting histogram represents
   count density of input and does not have norm 1.
* `:probability`: Normalize by sum of weights only. Resulting histogram
   represents the fraction of probability mass for each bin and does not have
   norm 1.
*  `:none`: Do not normalize.
Weighted data is supported via the keyword `weights`.

!!! note

    Normalizations are computed withing groups. For example, in the case of
    `normalization=:pdf`, sum of weights *within each group* will be equal to `1`.
"""
histogram(; options...) = transformation(HistogramAnalysis(; options...))


text:
@testset "radial_profile" begin

    X = [0 1 1 1 0
         1 2 2 2 1
         1 2 3 2 1
         1 2 2 2 1
         0 1 1 1 0];

    r, prof = @inferred radial_profile(X);
    @test r == [0, 1, 2, 3]
    @test prof ≈ [3, 2, 1, 0]

    X = [1 2
         2 2]
    r, prof = @inferred radial_profile(X, center=(1, 1))
    @test r == [0, 1]
    @test prof ≈ [1, 2]

    img = randn(rng, 1001, 1001)
    r, prof = @inferred radial_profile(img)
    @test mean(prof) ≈ 0 atol=1e-2
end


text:
"""
NsolPDE(n; sham=1, resdec=.5, rtol=1.e-7, atol=1.e-10)

Solve the Elliptic PDE using nsol.jl on an n x n grid. You give me
n and (optionally) sham and resdec and I return the output of nsol.
"""
function NsolPDE(n; sham = 1, resdec = 0.5, rtol = 1.e-7, atol = 1.e-10)
    # Get some room for the residual
    u0 = zeros(n * n)
    FV = copy(u0)
    # Get the precomputed data from pdeinit
    pdata = pdeinit(n)
    # Storage for the Jacobian
    J = copy(pdata.D2)
    # Call the solver
    hout = nsol(
        pdeF!,
        u0,
        FV,
        J,
        pdeJ!;
        rtol = rtol,
        atol = atol,
        pdata = pdata,
        sham = sham,
        resdec = resdec,
    )
    return hout
end


text:
mutable struct B
    id::Int
    name::String
end

import Base.==
==(a::B, b::B) = a.id == b.id && a.name == b.name

mutable struct A
    int8::Int8
    int::Int
    float::Float64
    str::String
    nullint::Nullable{Int}
    nullnullint::Nullable{Int}
    nullstr::Nullable{String}
    nullnullstr::Nullable{String}
    void::Void
    truebool::Bool
    falsebool::Bool
    b::B

    ints::Vector{Int}
    emptyarray::Vector{Int}
end
mutable struct E
    bs::Vector{B}
    dict::Dict{String,Int}
    emptydict::Dict{String,Int}
end

b1 = B(1, "harry")
b2 = B(2, "hermione")
b3 = B(3, "ron")

a = A(0, -1, 3.14, "string \\\" w/ escaped double quote", Nullable(4), Nullable{Int}(),
        Nullable("null string"), Nullable{String}(), nothing, true, false, b1, [1,2,3], Int[])#, [b2, b3],
        #Dict("1"=>1, "2"=>2), Dict{String,Int}())

io = IOBuffer()
Endpoints.JSON2.write(io, a)
seekstart(io)
a2 = Endpoints.JSON2.read(A, io)

@test a.int8 == a2.int8
@test a.int == a2.int
@test a.float == a2.float
@test a.str == a2.str
@test get(a.nullint) == get(a2.nullint)
@test isnull(a.nullnullint) && isnull(a2.nullnullint)
@test get(a.nullstr) == get(a2.nullstr)
@test isnull(a.nullnullstr) && isnull(a2.nullnullstr)
@test a.void == a2.void
@test a.truebool == a2.truebool
@test a.falsebool == a2.falsebool
@test a.b == a2.b
@test a.ints == a2.ints
@test a.emptyarray == a2.emptyarray
# @test a.bs == a2.bs
# @test a.dict == a2.dict
# @test a.emptydict == a2.emptydict


text:
using WignerSymbols
using JLD2

Nmax = 250;
A_lnm = zeros(Float32,2Nmax+1,Nmax,2Nmax+1);
B_lnm = zeros(Float32,2Nmax+1,Nmax,2Nmax+1);

#arr_to_fill = Array{Any}(nothing, l+n+1, n+1);

function f(Nmax)
    cN = 0;
    cM = 0;
    for l = 0:2Nmax
        for n = 1:Nmax
            mᵒ    = max(l-n,n+1)
            m_max = min(l+n,Nmax)
            #println("l,n: ", l, " ", n, " " , m_max-mᵒ, " " , (n+l)-mᵒ)
            cN += max(0,m_max-mᵒ)
            cM += max(0,(n+l)-mᵒ)
            #= for m = mᵒ:n+l
                println()
                fac = sqrt(2(2l+1)*(2n+1)*(2m+1));
                #println(m," ",l, " ", n)
                @inbounds A_lnm[l+1,n,m-n+1] = fac*wigner3j(m,n,l,-1,1,0);
                if l>1
                    @inbounds B_lnm[l+1,n,m-n+1] = fac*wigner3j(m,n,l,-1,-1,2);
                end
            end
            println(l," ",n) =#
        end
    end
    return cN,cM
end

global
global

@save "wignerTable.jld2" A_lnm B_lnm


text:
function tellYourWorld()
    print("Halo World!\n")
    print("Everyone,Creator!\n")
end

tellYourWorld()

max_stars_repo_name: mageGf/Hello-world

max_issues_repo_name: mageGf/Hello-world

max_forks_repo_name: mageGf/Hello-world


text:
module StressTest

Base.Experimental.@optlevel 0

"""
    busywait(seconds)

Like Base.sleep() except maxes out the thread for a specified number of seconds. The minimum busywait time is 1
millisecond or input of `0.001`.
"""
function busywait(sec::Real; yield_loop = false)
    sec ≥ 0 || throw(ArgumentError("cannot busywait for $sec seconds"))
    ns = sec * 1e9
    t = Base.time_ns()
    while Base.time_ns() - t < ns
        yield_loop && yield()
    end
    nothing
end

const TEST_MAT = rand(1000,1000)
testfunc() = TEST_MAT .* TEST_MAT

function check_parallel_threading(nthreads = Threads.nthreads())
    times = Vector{Float64}(undef, nthreads)
    t_loop = @elapsed Threads.@threads for i in 1:nthreads
        times[i] = @elapsed testfunc()
    end
    t_max = maximum(times)
    concurrency = t_loop <= (t_max * ((nthreads + 1) / nthreads))
    @show t_loop, t_max, concurrency
    return concurrency
end

export busywait, check_parallel_threading

function _precompile()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(busywait, (Int,))
    precompile(testfunc, ())
    precompile(check_parallel_threading, ())
end
_precompile()

end # module


text:
using AutoAccessorMethods
using AutoAccessorMethods: @noop
using Test

@testset "AutoAccessorMethods.jl" begin
    @testset "Accessor methods" begin
        include("test_accessor_methods.jl")
    end

    @testset "Exports" begin
        @testset "Don't export" begin
            include("test_do_not_export_methods.jl")
        end
        @testset "Do export" begin
            include("test_export_methods.jl")
        end
    end
end


text:
using SparseArrays, Test, JuMP, Polyhedra

function _hypercube(n, base_name)
    model = Model()
    @variable(model, [1:n], lower_bound = 0, upper_bound = 1, base_name = base_name)
    h = hrep(model)
    @test dimension_names(h) == ["$base_name[$i]" for i in 1:n]
    return h
end
function hypercube_name()
    h = _hypercube(1, "x") * _hypercube(2, "y")
    @test dimension_names(h) == ["x[1]", "y[1]", "y[2]"]
    @test dimension_names(spzeros(3, 4) \ h) === nothing
end
@testset "Hypercube name" begin
    hypercube_name()
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
struct Kron{T,N,I} <: AbstractArray{T,N}
    arrays::I
    function Kron{T,N,I}(A::I) where {T,N,I}
        isempty(A) && throw(ArgumentError("Cannot take kronecker product of empty vectors"))
        new{T,N,I}(A)
    end
end

_kron_dims() = 0
_kron_dims(A, B...) = max(ndims(A), _kron_dims(B...))

Kron{T,N}(A::AbstractArray...) where {T,N} = Kron{T,N,typeof(A)}(A)
Kron{T}(A::AbstractArray...) where {T} = Kron{T,_kron_dims(A...)}(A...)
Kron(A...) = Kron{promote_type(eltype.(A)...)}(A...)

size(K::Kron, j::Int) = prod(size.(K.arrays, j))
size(a::Kron{<:Any,1}) = (size(a,1),)
size(a::Kron{<:Any,2}) = (size(a,1), size(a,2))
size(a::Kron{<:Any,N}) where {N} = (@_inline_meta; ntuple(M -> size(a, M), Val(N)))

getindex(K::Kron{T,1,<:Tuple{<:AbstractVector}}, k::Int) where T =
    first(K.arrays)[k]

function getindex(K::Kron{T,1,<:Tuple{<:AbstractVector,<:AbstractVector}}, k::Int) where T
    A,B = K.arrays
    K,κ = divrem(k-1, length(B))
    A[K+1]*B[κ+1]
end

getindex(K::Kron{T,2,<:Tuple{<:AbstractMatrix}}, k::Int, j::Int) where T =
    first(K.arrays)[k,j]

function getindex(K::Kron{T,2,<:Tuple{<:AbstractArray,<:AbstractArray}}, k::Int, j::Int) where T
    A,B = K.arrays
    K,κ = divrem(k-1, size(B,1))
    J,ξ = divrem(j-1, size(B,2))
    A[K+1,J+1]*B[κ+1,ξ+1]
end

## Adapted from julia/stdlib/LinearAlgebra/src/dense.jl kron definition
function _kron2!(R, K)
    size(R) == size(K) || throw(DimensionMismatch("Matrices have wrong dimensions"))
    a,b = K.arrays
    @assert !has_offset_axes(a, b)
    m = 1
    @inbounds for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
        aij = a[i,j]
        for k = 1:size(b,1)
            R[m] = aij*b[k,l]
            m += 1
        end
    end
    R
end

copyto!(R::AbstractMatrix, K::Kron{<:Any,2,<:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) =
    _kron2!(R, K)
copyto!(R::AbstractVector, K::Kron{<:Any,1,<:Tuple{<:AbstractVector,<:AbstractVector}}) =
    _kron2!(R, K)

struct Diff{T, N, Arr} <: AbstractArray{T, N}
    v::Arr
    dims::Int
end

Diff(v::AbstractVector{T}) where T = Diff{T,1,typeof(v)}(v, 1)

function Diff(A::AbstractMatrix{T}; dims::Integer) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Diff{T,2,typeof(A)}(A, dims)
end

IndexStyle(::Type{<:Diff{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Diff{<:Any,2}}) = IndexCartesian()

size(D::Diff{<:Any,1}) = (length(D.v)-1,)
function size(D::Diff{<:Any,2})
    if D.dims == 1
        (size(D.v,1)-1,size(D.v,2))
    else #dims == 2
        (size(D.v,1),size(D.v,2)-1)
    end
end

getindex(D::Diff{<:Any, 1}, k::Integer) = D.v[k+1] - D.v[k]
function getindex(D::Diff, k::Integer, j::Integer)
    if D.dims == 1
        D.v[k+1,j] - D.v[k,j]
    else # dims == 2
        D.v[k,j+1] - D.v[k,j]
    end
end

struct Cumsum{T, N, Arr} <: AbstractArray{T, N}
    v::Arr
    dims::Int
end

Cumsum(v::AbstractVector{T}) where T = Cumsum{T,1,typeof(v)}(v, 1)

function Cumsum(A::AbstractMatrix{T}; dims::Integer) where T
    dims == 1 || dims == 2 || throw(ArgumentError("dimension must be 1 or 2, got $dims"))
    Cumsum{T,2,typeof(A)}(A, dims)
end

IndexStyle(::Type{<:Cumsum{<:Any,1}}) = IndexLinear()
IndexStyle(::Type{<:Cumsum{<:Any,2}}) = IndexCartesian()

size(Q::Cumsum) = size(Q.v)

getindex(Q::Cumsum{<:Any, 1}, k::Integer) = k == 1 ? Q.v[1] : Q.v[k] + Q[k-1]
function getindex(Q::Cumsum, k::Integer, j::Integer)
    if Q.dims == 1
        k == 1 ? Q.v[1,j] : Q.v[k,j] + Q[k-1,j]
    else # dims == 2
        j == 1 ? Q.v[k,1] : Q.v[k,j] + Q[k,j-1]
    end
end


text:
# TODO
# Add browser specific options
"""
	Capabilities(browserName::AbstractString;
                 browserVersion::AbstractString = "",
                 platformName::AbstractString = "",
                 acceptInsecureCerts::Union{Nothing, Bool} = nothing,
                 pageLoadStrategy::AbstractString = "",
                 proxy::Union{Nothing, Proxy} = nothing,
                 setWindowRect::Union{Nothing, Bool} = nothing,
                 timeouts::Timeouts = Timeouts(),
                 strictFileInteractability::Union{Nothing, Bool} = nothing,
                 unhandledPromptBehavior::AbstractString = "dismiss and notify"
                 )::Capabilities

According to [W3C](https://www.w3.org/TR/webdriver/#dfn-extension-capability).

# Example

```jldoctest
julia> capabilities = Capabilities("chrome")
Remote WebDriver Capabilities
browserName: chrome
```
"""
struct Capabilities{
    aIC<:Union{Nothing,Bool},
    P<:Union{Nothing,Proxy},
    sWR<:Union{Nothing,Bool},
    sFI<:Union{Nothing,Bool},
}
    browserName::String
    browserVersion::String
    platformName::String
    acceptInsecureCerts::aIC
    pageLoadStrategy::String
    proxy::P
    setWindowRect::sWR
    timeouts::Timeouts
    strictFileInteractability::sFI
    unhandledPromptBehavior::String
    function Capabilities(
        browserName::AbstractString;
        browserVersion::AbstractString = "",
        platformName::AbstractString = "",
        acceptInsecureCerts::Union{Nothing,Bool} = nothing,
        pageLoadStrategy::AbstractString = "",
        proxy::Union{Nothing,Proxy} = nothing,
        setWindowRect::Union{Nothing,Bool} = nothing,
        timeouts::Timeouts = Timeouts(),
        strictFileInteractability::Union{Nothing,Bool} = nothing,
        unhandledPromptBehavior::AbstractString = "dismiss and notify",
    )
        new{
            typeof(acceptInsecureCerts),
            typeof(proxy),
            typeof(setWindowRect),
            typeof(strictFileInteractability),
        }(
            browserName,
            browserVersion,
            platformName,
            acceptInsecureCerts,
            pageLoadStrategy,
            proxy,
            setWindowRect,
            timeouts,
            strictFileInteractability,
            unhandledPromptBehavior,
        )
    end
end
summary(io::IO, obj::Capabilities) = println(io, "Remote WebDriver Capabilities")
function show(io::IO, obj::Capabilities)
    print(io, summary(obj))
    println(io, "browserName: $(obj.browserName)")
end
JSON3.StructType(::Type{<:Capabilities}) = JSON3.Struct()
omitempties(::Type{<:Capabilities}) = fieldnames(Capabilities)


text:
@testset "multi_process" begin

end


text:
#InspectDR: Heatmap tools
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

#==Main types
===============================================================================#
struct ColorScale
	e::Vector{ARGB32} #Elements of color scale
end
ColorScale() = ColorScale([ARGB32(0, 0, 0, 0)])

#Transform used to map a 1D "z" value to a colour.
struct Transform1DToARGB
	vmin::DReal; span::DReal #Min value; span of colour scale
	s::DReal #scaling factor
	colorscale::ColorScale
end
function Transform1DToARGB(colorscale::ColorScale; vmin::DReal=0, vmax::DReal=1)
	vmin, vmax = minmax(vmin, vmax)
	span = vmax - vmin
	N = length(colorscale.e)
	s = (N-1) / span
	return Transform1DToARGB(vmin, span, s, colorscale)
end

#==Constructor-like functions
===============================================================================#

#Construct ColorScale from vector of color values
#ex: Colors.colormap("blues") returns a vector of RGB{Float64} values
function ColorScale(v::Vector{T}) where T<:Colorant
	return ColorScale([ARGB32(c) for c in v])
end

#==Main functions
===============================================================================#
function apply(xf::Transform1DToARGB, d::Array{T}) where T<:DReal
	result = Array{ARGB32}(undef,size(d))
	for i in 1:length(d)
		v = max(min(d[i] - xf.vmin, xf.span), 0)
		idx = Int(round(v*xf.s)) + 1
		result[i] = xf.colorscale.e[idx]
	end
	return result
end

#Last line


text:
"""
    p(sem_fit::SemFit)

Return the p value computed from the χ² test statistic.
"""
p_value(sem_fit::SemFit) = 1 - cdf(Chisq(df(sem_fit)), χ²(sem_fit))


text:
using Elly
using Base.Test

function test_yarn(host="localhost", port=8032)
    yarnclnt = YarnClient(host, port)
    nnodes = nodecount(yarnclnt)
    @test nnodes > 0
    println("number of yarn nodes: $nnodes")

    nlist = nodes(yarnclnt)
    show(STDOUT, nlist)
end


text: Dict{String, Any}("clients" => Dict{String, Any}("hosts" => Any[Dict{String, Any}("value" => "alpha", "type" => "string"), Dict{String, Any}("value" => "omega", "type" => "string")], "data" => Any[Any[Dict{String, Any}("value" => "gamma", "type" => "string"), Dict{String, Any}("value" => "delta", "type" => "string")], Any[Dict{String, Any}("value" => "1", "type" => "integer"), Dict{String, Any}("value" => "2", "type" => "integer")]]), "servers" => Dict{String, Any}("alpha" => Dict{String, Any}("dc" => Dict{String, Any}("value" => "eqdc10", "type" => "string"), "ip" => Dict{String, Any}("value" => "10.0.0.1", "type" => "string")), "beta" => Dict{String, Any}("dc" => Dict{String, Any}("value" => "eqdc10", "type" => "string"), "ip" => Dict{String, Any}("value" => "10.0.0.2", "type" => "string"))), "owner" => Dict{String, Any}("name" => Dict{String, Any}("value" => "Lance Uppercut", "type" => "string"), "dob" => Dict{String, Any}("value" => "1979-05-27T07:32:00Z", "type" => "datetime")), "title" => Dict{String, Any}("value" => "TOML Example", "type" => "string"), "database" => Dict{String, Any}("server" => Dict{String, Any}("value" => "192.168.1.1", "type" => "string"), "connection_max" => Dict{String, Any}("value" => "5000", "type" => "integer"), "enabled" => Dict{String, Any}("value" => "true", "type" => "bool"), "ports" => Any[Dict{String, Any}("value" => "8001", "type" => "integer"), Dict{String, Any}("value" => "8001", "type" => "integer"), Dict{String, Any}("value" => "8002", "type" => "integer")]))


text:
using myJuliaUtils
using Test
using LinearAlgebra
using DelimitedFiles
using Statistics

path = dirname(@__FILE__)
# path = "/Users/zymantas/TVPVARPkg/test/testAcf"

y = readdlm("$path/y.txt", ',');
acf2match = readdlm("$path/acf.txt", ',');

@test acf(y[:,1], 0:40) ≈ acf2match[:,1]

@test acf(y, 0:40) ≈ acf2match

for i=1:40
    @test acf(y, i) ≈ acf2match[i+1,:]'
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
import .ExpressionExplorer: SymbolsState, FunctionName, FunctionNameSignaturePair, try_compute_symbolreferences

"Every cell is a node in the reactive graph. The nodes/point/vertices are the _cells_, and the edges/lines/arrows are the _dependencies between cells_. In a reactive notebook, these dependencies are the **global variable references and definitions**. (For the mathies: a reactive notebook is represented by a _directed multigraph_. A notebook without reactivity errors is an _acyclic directed multigraph_.) This struct contains the back edges (`references`) and forward edges (`definitions`, `funcdefs_with_signatures`, `funcdefs_without_signatures`) of a single node.

Before 0.12.0, we could have written this struct with just two fields: `references` and `definitions` (both of type `Set{Symbol}`) because we used variable names to form the reactive links. However, to support defining _multiple methods of the same function in different cells_ (https://github.com/fonsp/Pluto.jl/issues/177), we needed to change this. You might want to think about this old behavior first (try it on paper) before reading on.

The essential idea is that edges are still formed by variable names. Simple global variables (`x = 1`) are registered by their name as `Symbol`, but _function definitions_ `f(x::Int) = 5` are sometimes stored in two ways:
- by their name (`f`) as `Symbol`, in `funcdefs_without_signatures`, and
- by their name with its method signature as `FunctionNameSignaturePair`, in `funcdefs_with_signatures`.

The name _without_ signature is most important: it is used to find the reactive dependencies between cells. The name _with_ signature is needed to detect multiple cells that define methods with the _same_ signature (`f(x) = 1` and `f(x) = 2`) - this is illegal. This is why we do not collect `definitions`, `funcdefs_with_signatures` and `funcdefs_without_signatures` onto a single pile: we need them separately for different searches.
"
Base.@kwdef struct ReactiveNode
    references::Set{Symbol} = Set{Symbol}()
    definitions::Set{Symbol} = Set{Symbol}()
    funcdefs_with_signatures::Set{FunctionNameSignaturePair} = Set{FunctionNameSignaturePair}()
    funcdefs_without_signatures::Set{Symbol} = Set{Symbol}()
    macrocalls::Set{Symbol} = Set{Symbol}()
end

function Base.union!(a::ReactiveNode, bs::ReactiveNode...)
	union!(a.references, (b.references for b in bs)...)
	union!(a.definitions, (b.definitions for b in bs)...)
	union!(a.funcdefs_with_signatures, (b.funcdefs_with_signatures for b in bs)...)
	union!(a.funcdefs_without_signatures, (b.funcdefs_without_signatures for b in bs)...)
	union!(a.macrocalls, (b.macrocalls for b in bs)...)
	return a
end

"Turn a `SymbolsState` into a `ReactiveNode`. The main differences are:
- A `SymbolsState` is a nested structure of function definitions inside function definitions inside... This conversion flattens this structure by merging `SymbolsState`s from defined functions.
- `ReactiveNode` functions as a cache to improve efficienty, by turning the nested structures into multiple `Set{Symbol}`s with fast lookups."
function ReactiveNode(symstate::SymbolsState)
	macrocalls = join_funcname_parts.(symstate.macrocalls) |> Set{Symbol}
	result = ReactiveNode(
		references=Set{Symbol}(symstate.references),
		definitions=Set{Symbol}(symstate.assignments),
		macrocalls=macrocalls,
		)

	# defined functions are 'exploded' into the cell's reactive node
	union!(result, (ReactiveNode(body_symstate) for (_, body_symstate) in symstate.funcdefs)...)

	# now we will add the function names to our edges:
	union!(result.references, symstate.funccalls .|> join_funcname_parts)
	union!(result.references, macrocalls)

	for (namesig, body_symstate) in symstate.funcdefs
		push!(result.funcdefs_with_signatures, namesig)
		push!(result.funcdefs_without_signatures, join_funcname_parts(namesig.name))
	end

	return result
end

# Convenience functions
ReactiveNode(code::String) = ReactiveNode(try_compute_symbolreferences(Meta.parse(code)))


text:
using JSON
using DataStructures
using NearestNeighbors
using Serialization
include("SpatialRegionTools.jl")

param  = JSON.parsefile("../hyper-parameters.json")
regionps = param[ARGS[2]]
cityname = regionps["cityname"]
cellsize = parse(Float64, ARGS[1])
file_path = regionps["filepath"]
if !isfile(file_path)
    println("Please provide the correct hdf5 file $file_path")
    exit(1)
end

region = SpatialRegion(cityname,
                       regionps["minlon"], regionps["minlat"],
                       regionps["maxlon"], regionps["maxlat"],
                       cellsize, cellsize,
                       regionps["minfreq"], # minfreq
                       40_000, # maxvocab_size
                       10, # k
                       4) # vocab_start

println("Building spatial region with:
        cityname=$(region.name),
        minlon=$(region.minlon),
        minlat=$(region.minlat),
        maxlon=$(region.maxlon),
        maxlat=$(region.maxlat),
        xstep=$(region.xstep),
        ystep=$(region.ystep),
        minfreq=$(region.minfreq)")

paramfile = "../data/$(region.name)-param-cell$(Int(cellsize))"
if isfile(paramfile)
    println("Reading parameter file from $paramfile")
    #loadregion!(region, paramfile)
    region = deserialize(open(paramfile, "r"))
else
    println("Creating paramter file $paramfile")
    num_out_region = makeVocab!(region, file_path)
    #saveregion(region, paramfile)
    serialize(open(paramfile, "w"), region)
end

#println("Making Vocabulary...")
#num_out_region = makeVocab!(region, "../data/$cityname.h5")

println("Vocabulary size $(region.vocab_size) with cell size $cellsize (meters)")
# createTrainVal(region, file_path, downsamplingDistort, 1_000_000, 10_000)
createTrainVal(region, file_path, downsamplingDistort, parse(Int, ARGS[3]), parse(Int, ARGS[4]), max_length=1000)
saveKNearestVocabs(region)

#region.cellcount |> keys |> maximum
#region.cellcount |> keys |> minimum
#region.cellcount |> length
#region.cellcount |> values |> sum
#region.cellcount |> values |> minimum
#region.cellcount |> values |> maximum
#region.vocab2hotcell[1638]
#region.vocab2hotcell |> keys |> maximum
#region.vocab_size
## test cell2gps(), gps2cell()
#for cell in 1:region.numx*region.numy-1
#    gps = cell2gps(region, cell)
#    @assert cell == gps2cell(region, gps...) "cell $cell bad"
#end

#cell = region.hotcell[10]
#kcells, dists = knearestHotcells(region, cell, 10)
#kvocabs = map(x->region.hotcell2vocab[x], kcells)
#[(v,d) for (v, d) in zip(kvocabs, dists)]

#saveKNearestVocabs(region)
#V, D = h5open("$cityname-vocab-dist.h5", "r") do f
#           f["V"]|>read, f["D"]|>read
#end

#seq = " 793 14350  8606 23750 13308" |> split |> x->map(parse, x)
#seq = [region.hotcell2vocab[c] for c in region.hotcell]
#trip = seq2trip(region, seq)
#h5open("/tmp/trip.h5", "w") do f f["trip"] = trip end

#a = sort(collect(region.cellcount), by=last, rev=true)[1:40000]
#region.hotcell
#a[end] |> last

max_stars_repo_name: lzzppp/DERT

max_issues_repo_name: lzzppp/DERT

max_forks_repo_name: lzzppp/DERT


text:
"""
Abstract problem.
"""
abstract type AbstractProblem{T} end

"""
    GeneralBTPDE(; model, matrices)

Bloch-Torrey PDE problem.
"""
Base.@kwdef struct BTPDE{T} <: AbstractProblem{T}
    model::Model{T}
    matrices::NamedTuple
end

"""
    HADC(; model, matrices)

HADC problem.
"""
Base.@kwdef struct HADC{T} <: AbstractProblem{T}
    model::Model{T}
    matrices::NamedTuple
end

"""
    Karger(; model, difftensors)

Karger problem.
"""
Base.@kwdef struct Karger{T} <: AbstractProblem{T}
    model::Model{T}
    difftensors::Vector{SMatrix{3,3,T,9}}
end

"""
    Laplace(; model, matrices, neig_max)

Laplace eigenvalue problem.
"""
Base.@kwdef struct Laplace{T} <: AbstractProblem{T}
    model::Model{T}
    matrices::NamedTuple
    neig_max::Int
end

"""
    MatrixFormalism(; model, matrices, lap_eig)

Matrix formalism problem. Given a Laplace eigendecomposition `lap_eig`, this problem
consists of computing the MF magnetization.
"""
Base.@kwdef struct MatrixFormalism{T} <: AbstractProblem{T}
    model::Model{T}
    matrices::NamedTuple
    lap_eig::NamedTuple
end

"""
    AnalyticalLaplace(; ρ, r, D, W, T₂, γ, dim, eiglim, eigstep)

Analytical radial Laplace eigenvalue problem.
"""
Base.@kwdef struct AnalyticalLaplace{T} <: AbstractProblem{T}
    ρ::Vector{Complex{T}}
    r::Vector{T}
    D::Vector{T}
    W::Vector{T}
    T₂::Vector{T}
    γ::T
    dim::Int
    eiglim::T
    eigstep::T
end

"""
    AnalyticalMatrixFormalism(; analytical_laplace, lap_mat, volumes)

Analytical Matrix formalism problem. Given a radial Laplace eigendecomposition
`analytical_laplace`, this problem consists of computing the MF compartment signals.
"""
Base.@kwdef struct AnalyticalMatrixFormalism{T} <: AbstractProblem{T}
    analytical_laplace::AnalyticalLaplace{T}
    lap_mat::NamedTuple
    volumes::Vector{T}
end


text:
# canonical data version: 1.2.0
# This file was generated by a script.

using Test

include("circular-buffer.jl")

@testset "reading empty buffer should fail" begin
    cb = CircularBuffer{Int}(1)
    @test_throws BoundsError popfirst!(cb)
end

@testset "can read an item just written" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test popfirst!(cb) == 1
end

@testset "each item may only be read once" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test popfirst!(cb) == 1
    @test_throws BoundsError popfirst!(cb)
end

@testset "items are read in the order they are written" begin
    cb = CircularBuffer{Int}(2)
    @test push!(cb, 1) === cb
    @test push!(cb, 2) === cb
    @test popfirst!(cb) == 1
    @test popfirst!(cb) == 2
end

@testset "full buffer can't be written to" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test_throws BoundsError push!(cb, 2)
end

@testset "a read frees up capacity for another write" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test popfirst!(cb) == 1
    @test push!(cb, 2) === cb
    @test popfirst!(cb) == 2
end

@testset "read position is maintained even across multiple writes" begin
    cb = CircularBuffer{Int}(3)
    @test push!(cb, 1) === cb
    @test push!(cb, 2) === cb
    @test popfirst!(cb) == 1
    @test push!(cb, 3) === cb
    @test popfirst!(cb) == 2
    @test popfirst!(cb) == 3
end

@testset "items cleared out of buffer can't be read" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test empty!(cb) === cb
    @test_throws BoundsError popfirst!(cb)
end

@testset "clear frees up capacity for another write" begin
    cb = CircularBuffer{Int}(1)
    @test push!(cb, 1) === cb
    @test empty!(cb) === cb
    @test push!(cb, 2) === cb
    @test popfirst!(cb) == 2
end

@testset "clear does nothing on empty buffer" begin
    cb = CircularBuffer{Int}(1)
    @test empty!(cb) === cb
    @test push!(cb, 1) === cb
    @test popfirst!(cb) == 1
end

@testset "overwrite acts like write on non-full buffer" begin
    cb = CircularBuffer{Int}(2)
    @test push!(cb, 1) === cb
    @test push!(cb, 2; overwrite=true) === cb
    @test popfirst!(cb) == 1
    @test popfirst!(cb) == 2
end

@testset "overwrite replaces the oldest item on full buffer" begin
    cb = CircularBuffer{Int}(2)
    @test push!(cb, 1) === cb
    @test push!(cb, 2) === cb
    @test push!(cb, 3; overwrite=true) === cb
    @test popfirst!(cb) == 2
    @test popfirst!(cb) == 3
end

@testset "overwrite replaces the oldest item remaining in buffer following a read" begin
    cb = CircularBuffer{Int}(3)
    @test push!(cb, 1) === cb
    @test push!(cb, 2) === cb
    @test push!(cb, 3) === cb
    @test popfirst!(cb) == 1
    @test push!(cb, 4) === cb
    @test push!(cb, 5; overwrite=true) === cb
    @test popfirst!(cb) == 3
    @test popfirst!(cb) == 4
    @test popfirst!(cb) == 5
end

@testset "initial clear does not affect wrapping around" begin
    cb = CircularBuffer{Int}(2)
    @test empty!(cb) === cb
    @test push!(cb, 1) === cb
    @test push!(cb, 2) === cb
    @test push!(cb, 3; overwrite=true) === cb
    @test push!(cb, 4; overwrite=true) === cb
    @test popfirst!(cb) == 3
    @test popfirst!(cb) == 4
    @test_throws BoundsError popfirst!(cb)
end

@testset "overwrite=true and overwrite=false" begin
    cb = CircularBuffer{Int}(2)
    @test empty!(cb) === cb
    @test push!(cb, 1; overwrite=true) === cb
    @test push!(cb, 2; overwrite=true) === cb
    @test_throws BoundsError push!(cb, 3; overwrite=false)
    @test push!(cb, 4; overwrite=true) === cb
    @test popfirst!(cb) == 2
    @test popfirst!(cb) == 4
    @test_throws BoundsError popfirst!(cb)
end

@testset "type parameter" begin
    cb = CircularBuffer{Int}(1)
    @test_throws Exception push!(cb, "0")
    cb = CircularBuffer{String}(1)
    @test_throws Exception push!(cb, 0)
    @test push!(cb, "0") === cb
    @test popfirst!(cb) == "0"
end

# Uncomment the following line to enable bonus tests.
# enable_bonus_tests = true

if @isdefined(enable_bonus_tests) && enable_bonus_tests
    println("\nBonus tests enabled.\n")

    @testset "is subtype of AbstractVector" begin
        @test CircularBuffer <: AbstractVector
    end

    # Copied from DataStructures.jl and slightly modified.
    @testset "Bonus test set taken from DataStructures.jl (CircularBuffer)" begin
        @testset "Core Functionality" begin
            cb = CircularBuffer{Int}(5)
            @testset "When empty" begin
                @test length(cb) == 0
                @test capacity(cb) == 5
                @test_throws BoundsError first(cb)
                @test_throws BoundsError last(cb)
                @test isempty(cb) == true
                @test isfull(cb) == false
                @test eltype(cb) == Int
                @test eltype(typeof(cb)) == Int
            end

            @testset "With 1 element" begin
                push!(cb, 1)
                @test length(cb) == 1
                @test capacity(cb) == 5
                @test isfull(cb) == false
                @test first(cb) == last(cb)
            end

            @testset "Appending many elements" begin
                append!(cb, 2:8; overwrite=true)
                @test length(cb) == capacity(cb)
                @test size(cb) == (length(cb),)
                @test isempty(cb) == false
                @test isfull(cb) == true
                @test convert(Array, cb) == Int[4,5,6,7,8]
            end

            @testset "getindex" begin
                @test cb[1] == 4
                @test cb[2] == 5
                @test cb[3] == 6
                @test cb[4] == 7
                @test cb[5] == 8
                @test_throws BoundsError cb[6]
                @test_throws BoundsError cb[3:6]
                @test cb[3:4] == Int[6,7]
                @test cb[[1,5]] == Int[4,8]
                @test first(cb) == 4
                @test last(cb) == 8
            end

            @testset "setindex" begin
                cb[3] = 999
                @test convert(Array, cb) == Int[4,5,999,7,8]
            end
        end

        @testset "pushfirst" begin
            cb = CircularBuffer{Int}(5)  # New, empty one for full test coverage
            for i in -5:5
                pushfirst!(cb, i; overwrite=true)
            end
            arr = convert(Array, cb)
            @test arr == Int[5, 4, 3, 2, 1]
            for (idx, n) in enumerate(5:1)
                @test arr[idx] == n
            end
        end

        @testset "Issue 429" begin
            cb = CircularBuffer{Int}(5)
            map(x -> pushfirst!(cb, x; overwrite=true), 1:8)
            pop!(cb)
            pushfirst!(cb, 9)
            arr = convert(Array, cb)
            @test arr == Int[9, 8, 7, 6, 5]
        end

        @testset "Issue 379" begin
            cb = CircularBuffer{Int}(5)
            pushfirst!(cb, 1)
            @test cb == [1]
            pushfirst!(cb, 2)
            @test cb == [2, 1]
        end

        @testset "empty!" begin
            cb = CircularBuffer{Int}(5)
            push!(cb, 13)
            empty!(cb)
            @test length(cb) == 0
        end

        @testset "pop!" begin
            cb = CircularBuffer{Int}(5)
            for i in 0:5    # one extra to force wraparound
                push!(cb, i; overwrite=true)
            end
            for j in 5:-1:1
                @test pop!(cb) == j
                @test convert(Array, cb) == collect(1:j-1)
            end
            @test isempty(cb)
            @test_throws BoundsError pop!(cb)
        end

        @testset "popfirst!" begin
            cb = CircularBuffer{Int}(5)
            for i in 0:5    # one extra to force wraparound
                push!(cb, i; overwrite=true)
            end
            for j in 1:5
                @test popfirst!(cb) == j
                @test convert(Array, cb) == collect(j+1:5)
            end
            @test isempty(cb)
            @test_throws BoundsError popfirst!(cb)
        end
    end
end


text:
using DelimitedFiles
using Plots
using LaTeXStrings
pyplot()

wavesInFile = readdlm("src/wavelist_formated.txt")

path_to_thresholds = "src/thresholds_formated.txt"

thresholds = fill(0.0,size(wavesInFile,1))
let path = path_to_thresholds
    if isfile(path) && size(wavesInFile,1)==88
        v = readdlm(path)
        for i in 1:size(v,1)
            global thresholds[Int64(v[i,1])] = v[i,2]
        end
    else
        warn("Do not consider thresholds!")
    end
end

# gr()
# default(fmt = :png)
let mv = 0.5:0.02:2.5
    bar(mv, [sum((thresholds[i]<m) for i in 1:88) for m in mv],
        xlab = L"M_{3π}\,\,\mathrm{bins}", ylab="Number of waves", lab="",
        c=:lightblue, lw=0, yticks=[0,40,61,88])
end
savefig("plots/thresholds.pdf")

# plot(rand(10))

# [sum((thresholds[i]<m) for i in 1:88) for m in 0.5:0.02:2.5]
get_color_palette


text:
using FMMTrees
using StaticArrays
using Test

const P = SVector{3,Float64}

using DelimitedFiles
Q = readdlm(joinpath(@__DIR__,"points.dlm"), Float64)
points = [SVector{3,Float64}(Q[i,:]) for i in axes(Q,1)]

root_center = P(0.5,0.5,0.5)
root_size = 0.5
root_node = FMMTrees.LevelledTrees.HNode(FMMTrees.PointerBasedTrees.Node(FMMTrees.LevelledTrees.Data(0,Int[]), 0, 0, 0, 0), 0)
tree = FMMTrees.LevelledTrees.LevelledTree([root_node], 1, root_center, root_size, Int[1])

smallest_box_size = 0.1
root_sector = 0
root_sfc_state = 1
root_depth = 1
for i in 1:length(points)
    router = FMMTrees.LevelledTrees.Router(smallest_box_size, points[i])
    root_state = root(tree), root_center, root_size, root_sfc_state, root_depth
    FMMTrees.update!(tree, root_state, i, router)  do tree, node, data
        push!(FMMTrees.data(tree,node).values, data)
    end
end

# FMMTrees.print_tree(tree)

num_bins = 3
bins = [P[] for i in 1:num_bins]

num_printed = 0
num_points = 0
num_nodes = length(tree.nodes)
for (i,node) in enumerate(FMMTrees.DepthFirstIterator(tree, root(tree)))
    # println(node, ": ", tree.nodes[node])
    b = div((i-1)*num_bins, num_nodes) + 1
    append!(bins[b], points[FMMTrees.data(tree,node).values])
    global num_printed += 1
    global num_points += length(data(tree,node).values)
end

ordered_points = reduce(append!, bins)
ordered_points = [p[i] for p in ordered_points, i = 1:3]
@test num_printed == length(tree.nodes)
@test num_points == length(points)

num_points = 0
ordered_points = P[]
for b in FMMTrees.LevelledTrees.LevelIterator(tree, length(tree.levels))
    global num_points += length(FMMTrees.data(tree,b).values)
    append!(ordered_points, points[FMMTrees.data(tree,b).values])
end
@test num_points == length(points)

@test length(tree.levels) == 4
@test FMMTrees.LevelledTrees.numlevels(tree) == 4
lvl1 = collect(FMMTrees.LevelledTrees.LevelIterator(tree,1))
@test lvl1 == [1]
lvl2 = collect(FMMTrees.LevelledTrees.LevelIterator(tree,2))
@test lvl2 == [10,13,2,5,23]
lvl3 = collect(FMMTrees.LevelledTrees.LevelIterator(tree,3))
@test lvl3 == [19,11,21,16,14,3,8,6,24]
lvl4 = collect(FMMTrees.LevelledTrees.LevelIterator(tree,4))
@test lvl4 == [20,18,12,22,17,15,4,9,7,25]


text:
#!/usr/bin/env julia

using JSON, CSV, DataFrames, DataStructures, ArgParse

function make_course_json(tees_file, pins_file, id, name)

    println("Note: par is initialized to 3 for all hole combinations. Edit resulting json manually.")

    holes_dict = OrderedDict()

    tees_df = CSV.read(tees_file, DataFrame)
    pins_df = CSV.read(pins_file, DataFrame)

    hole_names = unique(tees_df.hole)

    for h in hole_names
        hole_dict = OrderedDict()

        td = OrderedDict()
        for r in eachrow(tees_df)
            if r.hole == h
                ismissing(r.variation) ? key = "" : key = r.variation
                td[key] = [r.lat, r.lon]
            end
        end
        hole_dict[:tees] = td

        pd = OrderedDict()
        for r in eachrow(pins_df)
            if r.hole == h
                ismissing(r.variation) ? key = "" : key = r.variation
                pd[key] = [r.lat, r.lon]
            end
        end
        hole_dict[:pins] = pd

        sd = OrderedDict()
        for r in eachrow(tees_df)
            if r.hole == h
                ismissing(r.variation) ? key = "" : key = r.variation

                sdd = OrderedDict()
                for rr in eachrow(pins_df)
                    if rr.hole == h
                        ismissing(rr.variation) ? key2 = "" : key2 = rr.variation
                        sdd[key2] = 3
                    end
                end

                sd[key] = sdd
            end
        end
        hole_dict[:pars] = sd

        holes_dict[h] = hole_dict
    end

    D = OrderedDict()
    D[:id] = id
    D[:name] = name
    D[:loc] = [tees_df.lat[1], tees_df.lon[1]]
    D[:holes] = holes_dict

    open("$id.json", "w") do f
        JSON.print(f, D, 2)
    end

end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "tees"
            help = "Tees file"
            required = true
        "pins"
            help = "Pins file"
            required = true
        "id"
            help = "Course ID"
            required = true
        "name"
            help = "Course Name"
            required = true

    end
    return parse_args(s)
end

args = parse_commandline()

make_course_json(args["tees"], args["pins"], args["id"], args["name"])


text:
using Compat

mutable struct TreeSet{LABEL, RT, NL, N, B, TREE <:
                       AbstractTree{OneTree, RT, NL, N, B}} <:
    AbstractTree{ManyTrees, RT, NL, N, B}
    trees::Dict{LABEL, TREE}
    treeinfo::Dict{String, Dict{String, Any}}
end

import Compat.IteratorSize
function IteratorSize(::Type{TreeSet})
    return HasLength()
end

import Compat.IteratorEltype
function IteratorEltype(::Type{TreeSet})
    return HasEltype()
end

import Phylo.API: _ntrees
_ntrees(ts::TreeSet) = length(ts.trees)

import Phylo.API: _gettrees
_gettrees(ts::TreeSet) = values(ts.trees)

import Phylo.API: _gettreenames
_gettreenames(ts::TreeSet) = keys(ts.trees)

gettreeinfo(ts::TreeSet) = ts.treeinfo
gettreeinfo(ts::TreeSet, name) = ts.treeinfo[name]

import Phylo.API: _gettree
_gettree(ts::TreeSet, name) = ts.trees[name]

import Base.eltype
eltype(::TreeSet{LABEL, RT, NL, N, B, TREE}) where
    {LABEL, RT, NL, N, B,
     TREE <: AbstractTree{OneTree, RT, NL, N, B}} = TREE

import Base.length
length(ts::TreeSet) = length(ts.trees)

import Base.getindex
getindex(ts::TreeSet, idx) = ts.trees[idx]

import Phylo.API: _getleafnames
function _getleafnames(ts::TREESET) where {LABEL, RT, NL, N, B,
                    TREE <: AbstractTree{OneTree, RT, NL, N, B},
                    TREESET <: TreeSet{LABEL, RT, NL, N, B, TREE}}
    lns = unique(_getleafnames(t) for t in values(ts.trees))
    if length(lns) > 1
        error("Inconsistent leaf names in TreeSet")
    elseif isempty(lns)
        return NL[]
    end
    return first(lns)
end

import Phylo.API: _nleaves
function _nleaves(ts::TREESET) where {LABEL, RT, NL, N, B,
                    TREE <: AbstractTree{OneTree, RT, NL, N, B},
                    TREESET <: TreeSet{LABEL, RT, NL, N, B, TREE}}
    nls = unique(_nleaves(t) for t in values(ts.trees))
    if length(nls) > 1
        error("Inconsistent leafinfo in TreeSet")
    elseif isempty(nls)
        return 0
    end
    return first(nls)
end

import Base.keys
keys(ts::TreeSet) = keys(ts.trees)


text:
# module Angles

export
    wrap_to_2pi,
    wrap_to_pi,
    angular_difference

wrap_to_2pi(θ) = mod(θ,2π)

"""
    `wrap_to_pi(θ₀)`

    wraps the angle θ₀ to a value in (-π,π]
"""
function wrap_to_pi(θ₀)
    # θ = θ₀
    θ = wrap_to_2pi(θ₀)
    # while θ > π
    if θ > π
        θ = θ - 2π
    # end
    # while θ <= -π
    elseif θ <= -π
        θ += 2π
    end
    return θ
end

function angular_difference(θ1,θ2)
    t1 = wrap_to_2pi(θ1)
    t2 = wrap_to_2pi(θ2)
    if t2 < t1
        return -angular_difference(t2,t1)
    end
    wrap_to_pi(t2 - t1)
end

"""
    intersections_between_circles(v1,v2,r1,r2)
"""
function intersections_between_circles(v1,v2,r1,r2)
    d = norm(v1-v2)
    if d > r1+r2
        return nothing
    elseif d < 1e-4 # almost collocated
        return nothing
    elseif d + r1 < r2 || d + r2 < r1 # completely contained
        return nothing
    else
        x1,y1 = v1[1], v1[2]
        x2,y2 = v2[1], v2[2]
        l = (r1^2 - r2^2 + d^2)/(2*d)
        h = sqrt(r1^2-l^2)
        v0 = (l/d) * [(x2-x1), (y2-y1)] .+ v1
        dv = (h/d) * [(y2-y1), -(x2-x1)]

        return v0 .+ dv, v0 .- dv
    end
end

function nearest_points_between_circles(v1,v2,r1,r2)
    d = norm(v1-v2)
    vec = v2 .- v1
    if d > r1+r2
        return v1 + normalize(vec) * r1,  v2 - normalize(vec) * r2
    elseif d + r1 < r2 # circle 1 contained in circle 2
        return v1 - normalize(vec) * r1,  v2 - normalize(vec) * r2
    elseif d + r2 < r1 # circle 2 contained in circle 1
        return v1 + normalize(vec) * r1,  v2 + normalize(vec) * r2
    elseif d < 1e-4 # almost collocated
        return nothing
    else
        return intersections_between_circles(v1,v2,r1,r2)
    end
end


text:
using DataDeps

register(DataDep(
  "blood-transfusion-service-center",
  "http://archive.ics.uci.edu/ml/datasets/Blood+Transfusion+Service+Center",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/blood-transfusion/transfusion.data",
  "96c8e1091b9c037bcaf25a19b24b49d07771cd88689ffd273e056e9e8845ffe7",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=5,
        feature_indices=1:4,
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
### A Pluto.jl notebook ###
# v0.17.1

using Markdown
using InteractiveUtils

# ╔═╡ dc3c2030-144f-11eb-08ab-0b68dc2698ac
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using BenchmarkTools
end

# ╔═╡ 615d2b78-9a47-41bf-bae4-f0f7fb875956
# Make cells wider
html"""<style>
/*              screen size more than:                     and  less than:                     */
@media screen and (max-width: 699px) { /* Tablet */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/

}

@media screen and (min-width: 700px) and (max-width: 1199px) { /* Laptop*/
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1000px) and (max-width: 1920px) { /* Desktop */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1000px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}

@media screen and (min-width:1921px) { /* Stadium */
  /* Nest everything into here */
    main { /* Same as before */
        max-width: 1200px !important; /* Same as before */
        margin-right: 100px !important; /* Same as before */
    } /* Same as before*/
}
</style>
"""

# ╔═╡ 1c49d8ae-144e-11eb-024a-137460f78b15
md"""
# Performance
When building more complex programs, performance can be an issue. The [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/) website already contains a lot tips. We will illustrate potential performance gains for a function by providing one or more alternatives. In order to make sure the the functions produce the same result, the results are verified to be identical before starting the benchmarks.

The actual performance is measured using [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl). More advanced methods (that we cannot cover due to time constraints) included profiling your code. In Visual Studio Code you can get a detailed overview with the [`@profview`](https://discourse.julialang.org/t/julia-vs-code-extension-version-v0-17-released/42865) macro. For a REPL based method you could use [Traceur](https://github.com/JunoLab/Traceur.jl).

## Tools
Below you can find a set of very generic functions that we will be using to evaluate the functions that we will be benchmarking.

### Notes
At some point we try to make use of multithreading. You can see the number of threads that are currently being used with
```julia
Threads.nthreads()
```

Should this be only 1, you can set the number of threads to use by setting the environment variable `JULIA_NUM_THREADS` to the required before starting julia.

In a windows terminal:
```
$env:JULIA_NUM_THREADS = 4
julia
```
In a Linux/Mac OS terminal:
```
export JULIA_NUM_THREADS=4
julia
```

both examples above suppose "julia" is known in your terminal. If this is not the case, simply replace it by the path to the julia executable.

*Remark*: starting from julia 1.5 you can use command line arguments to specify the number of threads e.g. `julia --threads 4`
"""

# ╔═╡ a4ed00b8-144f-11eb-0c66-932242a1767c
begin
	"""
    testfuns(F, args...; kwargs...)

Benchmark the functions in F using `args` and `kwargs`.
"""
function testfuns(F::Union{Array{Function,1}, Array{Type,1}}, args...; kwargs...)
    # set up benchmark suite
    suite = BenchmarkGroup()
    for f in F
        suite[string(f)] = @benchmarkable $(f)($(args)...; $(kwargs)...)
    end
    # tune and run suite
    tune!(suite)
    results = run(suite, verbose=get(kwargs, :verbose, false))
    # show results
    for f in F
        println("\nResult for $(f):\n$("-"^length("Result for $(f):"))")
        display(results[string(f)])
    end
    medians = sort!([(fun, median(trial.times)) for (fun, trial) in results], by=x->x[2])
    speedup = round(Int,medians[end][2] / medians[1][2])
    println("\n$(speedup)x speedup (fastest: $(medians[1][1]), slowest: $(medians[end][1]))\n")
    return

end
	nothing
end

# ╔═╡ e32569b0-144f-11eb-33e9-8fd8755c208f
begin
	"""
		equality(F::Array{Function,1}, args...; kwargs...)

	Verify the result of an array of functions F is the same using `args` and `kwargs`.

	### Notes
	The first function in the array is used as the reference
	"""
	function equality(F::Array{Function,1}, args...; kwargs...)
		return all([equality(F[1], F[i], args...; kwargs...) for i in 2:length(F)])
	end

	"""
		equality(f::Function, g::Function, args...; kwargs...)

	Verify the result of function f and g is the same (element-wise) using `args` and `kwargs`.
	"""
	function equality(f::Function, g::Function, args...; kwargs...)
		return all(isequal.(map(f->f(args...; kwargs...), [f, g])...))
	end

	nothing
end

# ╔═╡ ff480c88-144f-11eb-2839-598b8452a945
md"""
## Hands-on
"""

# ╔═╡ 1e0ef46a-1450-11eb-2676-d94d5cc1f073
begin

	"""
		renew(char="", n::Int=0)

	Slow implementation of string joining
	"""
	function renew(char="", n::Int=0; kwargs...)
		res = ""
		for i in 1:n
			res *= char
		end

		return res
	end

	"""
		buffer(char="", n::Int=0)

	Fast implementation of string joining
	"""
	function buffer(char="", n::Int=0; kwargs...)
		res = IOBuffer()
		for i in 1:n
			print(res, char)
		end

		return String(take!(res))
	end

	"""
		demo1

	- Illustration of implementation issues

	speedup ∝ length
	"""
	function demo1()
		F = [renew, buffer];
		args=("a", 10000)
		# test validity
		@info "Testing demo 1"
		@assert equality(F, args...)
		# Run benchmark
		@info "Benching demo 1"
		testfuns(F, args...)
	end

	demo1()

end

# ╔═╡ 41ef249a-1450-11eb-2c52-0930531fa46e
begin
	function f_alloc(x::Vector{T}, y::Array{T,2}; props::Vector{T}, kwargs...) where T
		# read properties
		prop_A = props[1]
		prop_B = props[2]
		prop_C = props[3]
		prop_D = props[4]
		# initialise
		result = zeros(T, size(x))
		result .+= x
		# update result
		result = result .+ prop_A .* y[:,1]
		result = result .+ prop_B .* y[:,2]
		result = result .+ prop_C .* y[:,3]
		result = result .+ prop_D .* y[:,4]
		return result
	end

	struct Props{T<:AbstractFloat}
		prop_A::T
		prop_B::T
		prop_C::T
		prop_D::T
	end

	@views function f_noalloc(x::Vector{T}, y::Array{T,2}; p::Props, kwargs...) where T
		return x .+ p.prop_A .* y[:,1] .+ p.prop_B .* y[:,2] + p.prop_C .* y[:,3] + p.prop_D .* y[:,4]
	end

	"""
		demo2()

	- Avoiding unnecessary allocations
	- Using views
	- Avoiding fields with abstract types
	"""
	function demo2()
		F = [f_alloc, f_noalloc]
		N = 10000
		x = rand(N)
		y = rand(N,4)
		args = [x, y]
		props = Float64[1,2,3,4]
		p = Props(props...)
		kwargs = Dict(:props => props, :p => p)
		# test validity
		@info "Testing demo 2"
		@assert equality(F, args...; kwargs...)
		# run benchmark
		@info "Benching demo 2"
		testfuns(F, args...; kwargs...)
	end

	demo2()

end

# ╔═╡ 6332ef7e-1450-11eb-0a36-41818207b78f
begin
	# avoid changing type of a variable
	function f_typechange(x)
		res = 0
		for val in x
			res += val
		end
		return res
	end

	function f_notypechange(x::Array{T,1}) where T <: Number
		res = zero(T)
		for val in x
			res += val
		end
		return res
	end

	"""
		demo3()

	- avoid changing type of a variable
	"""
	function demo3()
		F = [f_typechange, f_notypechange]
		N = 2
		args = [rand(N)]
		# test validity
		@info "Testing demo 3"
		@assert equality(F, args...)
		# run benchmark
		@info "Benching demo 3"
		testfuns(F, args...)
	end

	demo3()
end

# ╔═╡ 75948388-1450-11eb-0a0e-431ad5929504
begin
	f_norm(x) = 3x.^2 + 4x + 7x.^3;
	f_vec(x) = @. 3x.^2 + 4x + 7x.^3;

	"""
		demo 4

	- use vectorisation where possible
	"""
	function demo4()
		F = [f_norm, f_vec]
		N::Int = 1e6
		args = [rand(N)]
		# test validity
		@info "Testing demo 4"
		@assert equality(F, args...)
		# run benchmark
		@info "Benching demo 4"
		testfuns(F, args...)
	end

	demo4()

end

# ╔═╡ 98e356f4-1450-11eb-0229-eb57f262ea4d
begin
	function f_nopre(x::Vector; kwargs...)
		y = similar(x)
		for i in eachindex(x)
			y[i] = sum(x[1:i])
		end
		return y
	end

	function f_pre!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
		for i in eachindex(x)
			y[i] = sum(x[1:i])
		end
		return y
	end

	function f_preboost!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
		Threads.@threads for i in eachindex(x)
			@inbounds y[i] = sum(x[1:i])
		end
		return y
	end

	@views function f_preturboboost!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
		Threads.@threads for i in eachindex(x)
			@inbounds y[i] = sum(x[1:i])
		end
		return y
	end

	"""
		demo5()

	- preallocate outputs if possible
	- make use of threads if possible
	- make use of views if possible
	- skip index check when sure of dimensions
	"""
	function demo5()
		F = [f_nopre, f_pre!, f_preboost!, f_preturboboost!]
		@info "Currently running on $(Threads.nthreads()) threads"
		N::Int = 1e4
		x::Vector = rand(N)
		y_pre = similar(x)
		args = [x]
		kwargs = Dict( :y => y_pre)
		# test validity
		@info "Testing demo 5"
		@assert equality(F, args...; kwargs...)
		# run benchmark
		@info "Benching demo 5"
		testfuns(F, args...; kwargs...)
	end

	demo5()
end

# ╔═╡ b341b4d2-1450-11eb-0549-7db5bd6756ba
begin
	function grow(N::Int; kwargs...)
		res = []
		for _ in 1:N
			push!(res, 1)
		end
		return res
	end

	function nogrow(N::Int; x::Array{T,1}) where{T}
		res = similar(x)
		for i in 1:N
			res[i] = 1
		end
		return res
	end

	"""
		demo6()

	- growing vs. allocating for known dimensions
	"""
	function demo6()
		N = 10000
		x = rand(N)
		F = [grow, nogrow]
		args = [N]
		kwargs = Dict(:x => x)
		# test validity
		@info "Testing demo 6"
		@assert equality(F, args...; kwargs...)
		# run benchmark
		@info "Benching demo 6"
		testfuns(F, args...; kwargs...)
	end

	demo6()
end

# ╔═╡ c6702048-1450-11eb-1f1d-a54089545949
begin
	struct mystruct
		a
		b
		c
	end

	struct mytypedstruct{T}
		a::T
		b::T
		c::T
	end

	mutable struct mymutablestruct
		a
		b
		c
	end

	"""
		demo7

	- creation time of mutable vs unmutable structs and the parameter version
	"""
	function demo7()
		structs = [mystruct, mytypedstruct, mymutablestruct]
		@info "Benching demo 7"
		testfuns(structs, [1;2;3]...)
	end

	demo7()

end

# ╔═╡ df4563da-1450-11eb-08a7-25ed365f578c
begin
	function reader(s::mystruct)
		return (s.a, s.b, s.c)
	end

	function reader(s::mytypedstruct)
		return (s.a, s.b, s.c)
	end

	function reader(sm::mymutablestruct)
		return (sm.a, sm.b, sm.c)
	end

	"""
		demo8

	- read times of fields of mutable vs unmutable structs and the parameter version
	"""
	function demo8()
		vals = [1;2;3]
		@info "Benching demo8"
		for arg in [mystruct, mytypedstruct, mymutablestruct]
			println("\nResult for $(arg):\n$("-"^length("Result for $(arg):"))")
			x = arg(vals...)
			display(@benchmark reader($(x)))
		end
	end

	demo8()
end

# ╔═╡ Cell order:
# ╟─615d2b78-9a47-41bf-bae4-f0f7fb875956
# ╟─1c49d8ae-144e-11eb-024a-137460f78b15
# ╠═dc3c2030-144f-11eb-08ab-0b68dc2698ac
# ╠═a4ed00b8-144f-11eb-0c66-932242a1767c
# ╠═e32569b0-144f-11eb-33e9-8fd8755c208f
# ╟─ff480c88-144f-11eb-2839-598b8452a945
# ╠═1e0ef46a-1450-11eb-2676-d94d5cc1f073
# ╠═41ef249a-1450-11eb-2c52-0930531fa46e
# ╠═6332ef7e-1450-11eb-0a36-41818207b78f
# ╠═75948388-1450-11eb-0a0e-431ad5929504
# ╠═98e356f4-1450-11eb-0229-eb57f262ea4d
# ╠═b341b4d2-1450-11eb-0549-7db5bd6756ba
# ╠═c6702048-1450-11eb-1f1d-a54089545949
# ╠═df4563da-1450-11eb-08a7-25ed365f578c

max_stars_repo_path: Exercises/PS10 - Performance.jl

max_issues_repo_path: Exercises/PS10 - Performance.jl

max_forks_repo_path: Exercises/PS10 - Performance.jl


text:
@scenario EmptyReservoirsScenario = begin
    ρ::PriceCurve{Float64}

    @zero begin
        return EmptyReservoirsScenario(PriceCurve(zeros(24)), probability = 1.0)
    end

    @expectation begin
        π = 1.0
        ρ = [mean([s.ρ[i] for s in scenarios]) for i in 1:length(scenarios[1].ρ)]
        return EmptyReservoirsScenario(PriceCurve(ρ), probability = 1.0)
    end
end

@sampler RecurrentEmptyReservoirsSampler = begin
    date::Date
    horizon::Horizon
    price_forecaster::Forecaster{:price}

    @sample EmptyReservoirsScenario begin
        price_curve = Vector{Float64}()
        for d = 1:nhours(sampler.horizon)
            daily_curve = forecast(sampler.price_forecaster, month(sampler.date + Dates.Day(d)))
            while !all(daily_curve .>= 0)
                daily_curve = forecast(sampler.price_forecaster, month(sampler.date + Dates.Day(d)))
            end
            append!(price_curve, daily_curve)
        end
        return EmptyReservoirsScenario(PriceCurve(price_curve))
    end
end


text:
using ImageFiltering, OffsetArrays, ImageCore, Random
import OffsetArrays.IdentityUnitRange

using Test

@testset "Border" begin
    borders = []
    push!(borders, Inner())
    push!(borders, Fill(0, (1,2), (3,4)))
    for _ in 1:10
        s = rand([:replicate, :circular, :symmetric, :reflect])
        pad = Pad(s, (rand(0:4),rand(0:4)), (rand(0:4), rand(0:4)))
        push!(borders, pad)
    end
    @testset "BorderArray consistent with padarray $b" for b in borders
        A = randn(5,5)
        B = BorderArray(A, b)
        @test padarray(A, b) == B
        @inferred B[1,1]
        for I in CartesianIndices(A)
            @test A[I] == B[I]
        end
    end
    @testset "BorderArray" begin
        @test_throws ArgumentError BorderArray(randn(3), Fill(0.))
        @test_throws ArgumentError BorderArray(randn(3), Fill(nothing, (1,)))
        @inferred BorderArray(randn(3), Fill(0., (1,)))
        @inferred BorderArray(randn(3), Fill(0, (1,)))
        ba1 = BorderArray(randn(3), Fill(zero(Int    ), (1,)))
        ba2 = BorderArray(randn(3), Fill(zero(Float64), (1,)))
        @test typeof(ba1) === typeof(ba2)
        @inferred ba1[1]
    end

    @testset "padarray" begin
        A = reshape(1:25, 5, 5)
        @test @inferred(padarray(A, Fill(0,(2,2),(2,2)))) == OffsetArray(
            [0  0  0   0   0   0   0  0  0;
             0  0  0   0   0   0   0  0  0;
             0  0  1   6  11  16  21  0  0;
             0  0  2   7  12  17  22  0  0;
             0  0  3   8  13  18  23  0  0;
             0  0  4   9  14  19  24  0  0;
             0  0  5  10  15  20  25  0  0;
             0  0  0   0   0   0   0  0  0;
             0  0  0   0   0   0   0  0  0], (-2,-2))
        @test @inferred(padarray(A, Pad(:replicate,(2,2),(2,2)))) == OffsetArray(
            [ 1  1  1   6  11  16  21  21  21;
              1  1  1   6  11  16  21  21  21;
              1  1  1   6  11  16  21  21  21;
              2  2  2   7  12  17  22  22  22;
              3  3  3   8  13  18  23  23  23;
              4  4  4   9  14  19  24  24  24;
              5  5  5  10  15  20  25  25  25;
              5  5  5  10  15  20  25  25  25;
              5  5  5  10  15  20  25  25  25], (-2,-2))
        @test @inferred(padarray(A, Pad(:circular,(2,2),(2,2)))) == OffsetArray(
            [19  24  4   9  14  19  24  4   9;
             20  25  5  10  15  20  25  5  10;
             16  21  1   6  11  16  21  1   6;
             17  22  2   7  12  17  22  2   7;
             18  23  3   8  13  18  23  3   8;
             19  24  4   9  14  19  24  4   9;
             20  25  5  10  15  20  25  5  10;
             16  21  1   6  11  16  21  1   6;
             17  22  2   7  12  17  22  2   7], (-2,-2))
        @test @inferred(padarray(A, Pad(:symmetric,(2,2),(2,2)))) == OffsetArray(
            [ 7  2  2   7  12  17  22  22  17;
              6  1  1   6  11  16  21  21  16;
              6  1  1   6  11  16  21  21  16;
              7  2  2   7  12  17  22  22  17;
              8  3  3   8  13  18  23  23  18;
              9  4  4   9  14  19  24  24  19;
             10  5  5  10  15  20  25  25  20;
             10  5  5  10  15  20  25  25  20;
              9  4  4   9  14  19  24  24  19], (-2,-2))
        @test @inferred(padarray(A, Pad(:reflect,(2,2),(2,2)))) == OffsetArray(
            [13   8  3   8  13  18  23  18  13;
             12   7  2   7  12  17  22  17  12;
             11   6  1   6  11  16  21  16  11;
             12   7  2   7  12  17  22  17  12;
             13   8  3   8  13  18  23  18  13;
             14   9  4   9  14  19  24  19  14;
             15  10  5  10  15  20  25  20  15;
             14   9  4   9  14  19  24  19  14;
             13   8  3   8  13  18  23  18  13], (-2,-2))
        ret = @test_throws ArgumentError padarray(A, Fill(0))
        @test occursin("lacks the proper padding", ret.value.msg)
        for Style in (:replicate, :circular, :symmetric, :reflect)
            ret = @test_throws ArgumentError padarray(A, Pad(Style,(1,1,1),(1,1,1)))
            @test occursin("lacks the proper padding", ret.value.msg)
            ret = @test_throws ArgumentError padarray(A, Pad(Style))
            @test occursin("lacks the proper padding", ret.value.msg)
        end
        # same as above, but input arrays are of type OffsetArray
        Ao =  OffsetArray(reshape(1:25, 5, 5),-1,-1)
        @test @inferred(padarray(Ao, Fill(0,(2,2),(2,2)))) == OffsetArray(
            [0  0  0   0   0   0   0  0  0;
             0  0  0   0   0   0   0  0  0;
             0  0  1   6  11  16  21  0  0;
             0  0  2   7  12  17  22  0  0;
             0  0  3   8  13  18  23  0  0;
             0  0  4   9  14  19  24  0  0;
             0  0  5  10  15  20  25  0  0;
             0  0  0   0   0   0   0  0  0;
             0  0  0   0   0   0   0  0  0], (-3,-3))
        @test @inferred(padarray(Ao, Pad(:replicate,(2,2),(2,2)))) == OffsetArray(
            [ 1  1  1   6  11  16  21  21  21;
              1  1  1   6  11  16  21  21  21;
              1  1  1   6  11  16  21  21  21;
              2  2  2   7  12  17  22  22  22;
              3  3  3   8  13  18  23  23  23;
              4  4  4   9  14  19  24  24  24;
              5  5  5  10  15  20  25  25  25;
              5  5  5  10  15  20  25  25  25;
              5  5  5  10  15  20  25  25  25], (-3,-3))
        @test @inferred(padarray(Ao, Pad(:circular,(2,2),(2,2)))) == OffsetArray(
            [19  24  4   9  14  19  24  4   9;
             20  25  5  10  15  20  25  5  10;
             16  21  1   6  11  16  21  1   6;
             17  22  2   7  12  17  22  2   7;
             18  23  3   8  13  18  23  3   8;
             19  24  4   9  14  19  24  4   9;
             20  25  5  10  15  20  25  5  10;
             16  21  1   6  11  16  21  1   6;
             17  22  2   7  12  17  22  2   7], (-3,-3))
        @test @inferred(padarray(Ao, Pad(:symmetric,(2,2),(2,2)))) == OffsetArray(
            [ 7  2  2   7  12  17  22  22  17;
              6  1  1   6  11  16  21  21  16;
              6  1  1   6  11  16  21  21  16;
              7  2  2   7  12  17  22  22  17;
              8  3  3   8  13  18  23  23  18;
              9  4  4   9  14  19  24  24  19;
             10  5  5  10  15  20  25  25  20;
             10  5  5  10  15  20  25  25  20;
              9  4  4   9  14  19  24  24  19], (-3,-3))
        @test @inferred(padarray(Ao, Pad(:reflect,(2,2),(2,2)))) == OffsetArray(
            [13   8  3   8  13  18  23  18  13;
             12   7  2   7  12  17  22  17  12;
             11   6  1   6  11  16  21  16  11;
             12   7  2   7  12  17  22  17  12;
             13   8  3   8  13  18  23  18  13;
             14   9  4   9  14  19  24  19  14;
             15  10  5  10  15  20  25  20  15;
             14   9  4   9  14  19  24  19  14;
             13   8  3   8  13  18  23  18  13], (-3,-3))
        ret = @test_throws ArgumentError padarray(Ao, Fill(0))
        @test occursin("lacks the proper padding", ret.value.msg)
        for Style in (:replicate, :circular, :symmetric, :reflect)
            ret = @test_throws ArgumentError padarray(Ao, Pad(Style,(1,1,1),(1,1,1)))
            @test occursin("lacks the proper padding", ret.value.msg)
            ret = @test_throws ArgumentError padarray(Ao, Pad(Style))
            @test occursin("lacks the proper padding", ret.value.msg)
        end
        # arrays smaller than the padding
        A = [1 2; 3 4]
        @test @inferred(padarray(A, Pad(:replicate,(3,3),(3,3)))) == OffsetArray(
            [ 1  1  1  1  2  2  2  2;
              1  1  1  1  2  2  2  2;
              1  1  1  1  2  2  2  2;
              1  1  1  1  2  2  2  2;
              3  3  3  3  4  4  4  4;
              3  3  3  3  4  4  4  4;
              3  3  3  3  4  4  4  4;
              3  3  3  3  4  4  4  4], (-3,-3))
        @test @inferred(padarray(A, Pad(:circular,(3,3),(3,3)))) == OffsetArray(
            [ 4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1], (-3,-3))
        @test @inferred(padarray(A, Pad(:symmetric,(3,3),(3,3)))) == OffsetArray(
            [ 4  4  3  3  4  4  3  3;
              4  4  3  3  4  4  3  3;
              2  2  1  1  2  2  1  1;
              2  2  1  1  2  2  1  1;
              4  4  3  3  4  4  3  3;
              4  4  3  3  4  4  3  3;
              2  2  1  1  2  2  1  1;
              2  2  1  1  2  2  1  1], (-3,-3))
        @test @inferred(padarray(A, Pad(:reflect,(3,3),(3,3)))) == OffsetArray(
            [ 4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1;
              4  3  4  3  4  3  4  3;
              2  1  2  1  2  1  2  1], (-3,-3))
        for Style in (:replicate, :circular, :symmetric, :reflect)
            @test @inferred(padarray(A, Pad(Style,(0,0), (0,0)))) == A
        end
        @test @inferred(padarray(A, Fill(0, (0,0), (0,0)))) == A
        @test @inferred(padarray(A, Pad(:replicate,(1,2), (2,0)))) == OffsetArray([1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4], (-1,-2))
        @test @inferred(padarray(A, Pad(:circular,(2,1), (0,2)))) == OffsetArray([2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4], (-2,-1))
        @test @inferred(padarray(A, Pad(:symmetric,(1,2), (2,0)))) == OffsetArray([2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2], (-1,-2))
        @test @inferred(padarray(A, Fill(-1,(1,2), (2,0)))) == OffsetArray([-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1], (-1,-2))
        A = [1 2 3; 4 5 6]
        @test @inferred(padarray(A, Pad(:reflect,(1,2), (2,0)))) == OffsetArray([6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6], (-1,-2))
        A = [1 2; 3 4]
        @test @inferred(padarray(A, Pad(1,1))) == OffsetArray([1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4], (-1,-1))
        @test @inferred(padarray(A, Pad(:circular,(1,1), ()))) == OffsetArray([4 3 4; 2 1 2; 4 3 4], (-1,-1))
        @test @inferred(padarray(A, Pad(:symmetric,(), (1,1)))) == [1 2 2; 3 4 4; 3 4 4]
        A = ["a" "b"; "c" "d"]
        @test @inferred(padarray(A, Pad(1,1))) == OffsetArray(["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"], (-1,-1))
        @test_throws ErrorException padarray(A, Pad(:unknown,1,1))
        A = trues(3,3)
        @test isa(parent(padarray(A, Pad((1,2), (2,1)))), BitArray{2})
        A = falses(10,10,10)
        B = view(A,1:8,1:8,1:8)
        @test isa(parent(padarray(A, Pad((1,1,1)))), BitArray{3})
        @test isa(parent(padarray(B, Pad((1,1,1)))), BitArray{3})
        A = reshape(1:15, 3, 5)
        B = @inferred(padarray(A, Inner((1,1))))
        @test B == A
        # test that it's a copy
        B[1,1] = 0
        @test B != A
        A = rand(RGB{N0f8}, 3, 5)
        ret = @test_throws ArgumentError padarray(A, Fill(0, (0,0), (0,0)))
        @test occursin("RGB", ret.value.msg)
        @test occursin("convert", ret.value.msg)
        A = bitrand(3, 5)
        ret = @test_throws ArgumentError padarray(A, Fill(7, (0,0), (0,0)))
        @test occursin("Bool", ret.value.msg)
        @test occursin("convert", ret.value.msg)
        @test isa(parent(padarray(A, Fill(false, (1,1), (1,1)))), BitArray)
    end

    @testset "Pad" begin
        @test Pad(:replicate,[1,2], [5,3]) == Pad(:replicate,(1,2), (5,3))
        @test @inferred(Pad{2}(:replicate, [1,2], [5,3])) == Pad(:replicate,(1,2), (5,3))
        @eval @test_throws MethodError Pad{3}(:replicate, [1,2], [5,3])
        @test @inferred(Pad(:circular)(rand(3,5))) == Pad(:circular, (0,0),(3,5))
        @test @inferred(Pad(:circular)(centered(rand(3,5)))) == Pad{2}(:circular, (1,2),(1,2))
        @test @inferred(Pad(:symmetric)(Kernel.Laplacian())) == Pad{2}(:symmetric, (1,1),(1,1))
        @test @inferred(Pad(:symmetric)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Pad(:symmetric,(1,1),(1,1))

        a = reshape(1:15, 3, 5)
        targetinds = (OffsetArray([1,1,2,3,3], 0:4), OffsetArray([1:5;], 1:5))
        ret = @test_throws ArgumentError ImageFiltering.padindices(rand(3,5), Pad(:replicate,(1,)))
        @test occursin("lacks the proper padding sizes", ret.value.msg)
    end

    @testset "Inner" begin
        @test @inferred(Inner((1,2))) == Inner((1,2), (1,2))
        @test @inferred(Inner((1,2), ())) == Inner((1,2), (0,0))
        @test @inferred(Inner((), (1,2))) == Inner((0,0), (1,2))
        @test @inferred(Inner{2}([1,2],[5,6])) == Inner((1,2), (5,6))
        @test Inner([1,2],[5,6]) == Inner((1,2), (5,6))
        @test @inferred(Inner()(rand(3,5))) == Inner((0,0),(3,5))
        @test @inferred(Inner()(centered(rand(3,5)))) == Inner((1,2),(1,2))
    end

    @testset "Fill" begin
        @test Fill(1,[1,2],[5,6]) == Fill(1, (1,2), (5,6))
        @test Fill(1,(3,2)) == Fill(1, (3,2), (3,2))
        @test @inferred(Fill(-1)(rand(3,5))) == Fill(-1, (0,0), (3,5))
        @test @inferred(Fill(-1)(centered(rand(3,5)))) == Fill(-1, (1,2), (1,2))
        @test @inferred(Fill(-1, rand(3,5))) == Fill(-1, (0,0), (3,5))
        @test @inferred(Fill(-1, centered(rand(3,5)))) == Fill(-1, (1,2), (1,2))
        @test @inferred(Fill(-1)(rand(3,5))) == Fill(-1, (0,0), (3,5))
        @test @inferred(Fill(-1)(centered(rand(3,5)))) == Fill(-1, (1,2), (1,2))
        @test @inferred(Fill(2)(Kernel.Laplacian(), rand(5,5), Algorithm.FIR())) == Fill(2, (1,1),(1,1))
    end

    @testset "misc" begin
        a0 = reshape([1])  # 0-dimensional
        @test ImageFiltering.accumulate_padding((), a0) == ()
        @test ImageFiltering.accumulate_padding((0:1, -1:1), a0) == (0:1, -1:1)
        # ImageFiltering #137 and #139
        b_front = BorderArray(zeros(10), Fill(1.0, (3,), (0,)))
        b_both = BorderArray(b_front, Fill(1.0, (0,), (5,)))
        b_sim = similar(b_both)
        @test axes(b_sim) == axes(b_both)
	@test axes(b_sim) === (IdentityUnitRange(-2:15),)
	@test axes(similar(b_sim, Int32, IdentityUnitRange(-1:3))) === (IdentityUnitRange(-1:3),)
        @test all(LinearIndices(b_front) .== LinearIndices(OffsetArray(zeros(13), -3)))
    end
end

nothing


text:
export AddStoppingCondition, SetIterationLimit, BumpIterationLimit, SetFunctionEvaluationLimit,
       BumpFunctionEvaluationLimit, SetMinimumMeshSize, SetMinimumPollSize, RuntimeStoppingCondition,
       SetButtonLimit,SetHypervolumeLimit,HypervolumeStoppingCondition, checkKeyInterrupt

"""
    AddStoppingCondition(p::DSProblem, c::T) where T <: AbstractStoppingCondition

Add the stopping condition `c` to the problem `p`.
"""
function AddStoppingCondition(p::DSProblem, c::T) where T <: AbstractStoppingCondition
    push!(p.stoppingconditions, c)
end

_check_stoppingconditions(p::DSProblem) = _check_stoppingconditions(p, p.stoppingconditions)
function _check_stoppingconditions(p::DSProblem, c::Vector{T}) where T <: AbstractStoppingCondition
    for condition in c
        #Stopping conditions return true if optimisation should continue,
        #false if the condition is met and optimisation should stop
        if CheckStoppingCondition(p, condition) == false
            setstatus(p, condition)
            return false
        end
    end
        !checkKeyInterrupt(p) && return false

    return true
end

function setstatus(p, s::T) where T <: AbstractStoppingCondition
    p.status.optimization_status_string = StoppingConditionStatus(s)

    if s isa IterationStoppingCondition
        p.status.optimization_status = IterationLimit
    elseif s isa MeshPrecisionStoppingCondition
        p.status.optimization_status = MeshPrecisionLimit
    elseif s isa PollPrecisionStoppingCondition
        p.status.optimization_status = PollPrecisionLimit
    elseif s isa FunctionEvaluationStoppingCondition
        p.status.optimization_status = FunctionEvaluationLimit
    elseif s isa RuntimeStoppingCondition
        p.status.optimization_status = RuntimeLimit
    else
        p.status.optimization_status = OtherStoppingCondition
    end
end

StoppingConditionStatus(::T) where T <: AbstractStoppingCondition = "Unknown stopping condition status"

function _init_stoppingconditions(p::DSProblem)
    for c in p.stoppingconditions
        init_stoppingcondition(p, c)
    end
end

init_stoppingcondition(p::DSProblem, ::AbstractStoppingCondition) = nothing

_get_conditionindexes(p::DSProblem, target::Type) = _get_conditionindexes(p.stoppingconditions, target)
_get_conditionindexes(s::Vector{AbstractStoppingCondition}, target::Type) =
    [i for (i,v) in enumerate(s) if typeof(v) == target]

#===== Built-in stopping conditions =====#
# Key Interrupt

# for MADS
function checkKeyInterrupt(p::DSProblem)::Bool
    bb = bytesavailable(stdin)
    data = read(stdin, bb)
    if bb>0 && data==UInt8[0x71] #q for quit
        println("quit")
        p.status.optimization_status=KeyInterrupt
        p.status.optimization_status_string = "Key Interrupt"
        return false
    end
    return true
end

#Iteration limit
mutable struct IterationStoppingCondition <: AbstractStoppingCondition
    limit::Int64
end

StoppingConditionStatus(::IterationStoppingCondition) = "Iteration limit"

CheckStoppingCondition(p::DSProblem, s::IterationStoppingCondition) = p.status.iteration < s.limit

function init_stoppingcondition(p::DSProblem, s::IterationStoppingCondition)
    if s.limit == -1
        error("Please set a maximum number of iterations")
    end
end

"""
    SetIterationLimit(p::DSProblem, i::Int)

Set the maximum number of iterations to `i`.
"""
function SetIterationLimit(p::DSProblem, i::Int)
    if i < p.status.iteration
        error("Cannot set iteration limit to lower than the number of iterations that have run")
    else
        iteration_indexes = _get_conditionindexes(p, IterationStoppingCondition)
        for index in iteration_indexes
            p.stoppingconditions[index].limit = i
        end
    end
end

"""
    BumpIterationLimit(p::DSProblem, i::Int)

Increase the iteration limit by `i`.
"""
function BumpIterationLimit(p::DSProblem, i::Int)
    iteration_indexes = _get_conditionindexes(p, IterationStoppingCondition)
    for index in iteration_indexes
        p.stoppingconditions[index].limit += i
    end
end

#Mesh size precision limit
mutable struct MeshPrecisionStoppingCondition{T} <: AbstractStoppingCondition
    cont_min_mesh_size::T
    min_mesh_sizes::Vector{T}

    function MeshPrecisionStoppingCondition{T}(min_mesh_size::Union{T, Nothing}=nothing) where T
        c = new()

        if min_mesh_size != nothing
            c.cont_min_mesh_size = min_mesh_size
        else
            c.cont_min_mesh_size = get_min_mesh_size(T)
        end

        return c
    end
end

function init_stoppingcondition(p::DSProblem, s::MeshPrecisionStoppingCondition)
    s.min_mesh_sizes = map(δ_min -> δ_min > 0 ? δ_min : s.cont_min_mesh_size, p.config.mesh.δ_min)
end

StoppingConditionStatus(::MeshPrecisionStoppingCondition) = "Mesh Precision limit"

function CheckStoppingCondition(p::DSProblem, s::MeshPrecisionStoppingCondition)
    if p.config.mesh.only_granular
        p.config.mesh.l > -50
    else
        any(GetMeshSizeVector(p) .> s.min_mesh_sizes)
    end
end

(GetMeshSizeVector(p::DSProblem{T})::Vector{T}) where T = p.config.mesh.δ

get_min_mesh_size(::Type{Float64}) = 1.1102230246251565e-16
get_min_mesh_size(::Type{T}) where T = eps(T)/2

"""
    SetMinimumMeshSize(p::DSProblem{T}, i::T) where T

Set the minimum mesh size for continuous variables.
"""
function SetMinimumMeshSize(p::DSProblem{T}, i::T) where T
    if i <= 0
        error("Minimum mesh size must be positive.")
    else
        mesh_precision_indexes = _get_conditionindexes(p, MeshPrecisionStoppingCondition)
        for index in mesh_precision_indexes
            p.stoppingconditions[index].cont_min_mesh_size = i
        end
    end
end

#Poll size precision limit
mutable struct PollPrecisionStoppingCondition{T} <: AbstractStoppingCondition
    cont_min_poll_size::T
    min_poll_sizes::Vector{T}

    function PollPrecisionStoppingCondition{T}(min_poll_size::Union{T, Nothing}=nothing) where T
        c = new()

        if min_poll_size != nothing
            c.cont_min_poll_size = min_poll_size
        else
            c.cont_min_poll_size = get_min_poll_size(T)
        end

        return c
    end
end

function init_stoppingcondition(p::DSProblem, s::PollPrecisionStoppingCondition)
    s.min_poll_sizes = map(δ_min -> δ_min > 0 ? δ_min : s.cont_min_poll_size, p.config.mesh.δ_min)
end

StoppingConditionStatus(::PollPrecisionStoppingCondition) = "Poll Precision limit"

function CheckStoppingCondition(p::DSProblem, s::PollPrecisionStoppingCondition)
    if p.config.mesh.only_granular
        p.config.mesh.l > -50
    else
        any(GetPollSizeVector(p) .> s.min_poll_sizes)
    end
end

(GetPollSizeVector(p::DSProblem{T})::Vector{T}) where T = p.config.mesh.Δ

get_min_poll_size(::Type{Float64}) = 1.1102230246251565e-16
get_min_poll_size(::Type{T}) where T = eps(T)/2

"""
    SetMinimumPollSize(p::DSProblem{T}, i::T) where T

Set the minimum poll size for continuous variables.
"""
function SetMinimumPollSize(p::DSProblem{T}, i::T) where T
    if i <= 0
        error("Minimum poll size must be positive.")
    else
        poll_precision_indexes = _get_conditionindexes(p, PollPrecisionStoppingCondition)
        for index in poll_precision_indexes
            p.stoppingconditions[index].cont_min_poll_size = i
        end
    end
end

#Function evaluation limit
mutable struct FunctionEvaluationStoppingCondition <: AbstractStoppingCondition
    limit::Int64
end

StoppingConditionStatus(::FunctionEvaluationStoppingCondition) = "Function evaluation limit"

CheckStoppingCondition(p::DSProblem, s::FunctionEvaluationStoppingCondition) = p.status.function_evaluations < s.limit

function init_stoppingcondition(::DSProblem, s::FunctionEvaluationStoppingCondition)
    if s.limit == -1
        error("Please set a maximum number of function evaluations")
    end
end

"""
    SetFunctionEvaluationLimit(p::DSProblem, i::Int)

Set the maximum number of function evaluations to `i`.
"""
function SetFunctionEvaluationLimit(p::DSProblem, i::Int)
    if i < p.status.function_evaluations
        error("Cannot set function evaluation limit to lower than the number of function evaluations that have run")
    else
        function_evaluation_indexes = _get_conditionindexes(p, FunctionEvaluationStoppingCondition)
        for index in function_evaluation_indexes
            p.stoppingconditions[index].limit = i
        end
    end
end

"""
    BumpFunctionEvaluationLimit(p::DSProblem, i::Int)

Increase the function evaluation limit by `i`.
"""
function BumpFunctionEvaluationLimit(p::DSProblem, i::Int)
    function_evaluation_indexes = _get_conditionindexes(p, FunctionEvaluationStoppingCondition)
    for index in function_evaluation_indexes
        p.stoppingconditions[index].limit += i
    end
end

#===== Optional stopping conditions =====#

#Runtime limit
mutable struct RuntimeStoppingCondition <: AbstractStoppingCondition
    limit::Float64
end

StoppingConditionStatus(::RuntimeStoppingCondition) = "Runtime limit"

CheckStoppingCondition(p::DSProblem, s::RuntimeStoppingCondition) = (time() - p.status.start_time) < s.limit

function init_stoppingcondition(::DSProblem, s::RuntimeStoppingCondition)
    if s.limit <= 0
        error("Runtime limit must be positive.")
    end
end

"""
    SetRuntimeLimit(p::DSProblem, i::Float64)

Set the runtime limit to `i`.
"""
function SetRuntimeLimit(p::DSProblem, i::Float64)
    if i < p.status.runtime_total
        error("Cannot set runtime limit to lower than the runtime of the previous run")
    elseif i <= 0
        error("Runtime limit has to be positive.")
    else
        runtime_indexes = _get_conditionindexes(p, RuntimeStoppingCondition)
        for index in runtime_indexes
            p.stoppingconditions[index].limit = i
        end
    end
end

#Pareto front hyper-volume limit (Bimads only)
mutable struct HypervolumeStoppingCondition <: AbstractStoppingCondition
    limit::Float64
end

StoppingConditionStatus(::HypervolumeStoppingCondition) = "Hyper-Volume limit"
CheckStoppingCondition(p::DSProblem, s::HypervolumeStoppingCondition) = true

function init_stoppingcondition(::DSProblem, s::HypervolumeStoppingCondition)
    if s.limit <= 0
        error("Hypervolume limit must be positive.")
    end
end

"""
    SetHypervolumeLimit(p::DSProblem, i::Float64)

Set the Hyper-Volume limit to `i`.
"""
function SetHypervolumeLimit(p::DSProblem, i::Float64)
    if i <= 0
        error("Hyper-Volume limit has to be positive.")
    else
        runtime_indexes = _get_conditionindexes(p, HypervolumeStoppingCondition)
        for index in runtime_indexes
            p.stoppingconditions[index].limit = i
        end
    end
end