text:
# This file was generated, do not modify it. # hide
tuning = Grid(resolution=7)
resampling = CV(nfolds=6)

tm = TunedModel(model=rfr, tuning=tuning,
                resampling=resampling, ranges=r, measure=rms)

rfr_tm = machine(tm, X, y);


text:
@testset "Phase and BoundedFluidQueue" begin
    @testset "PhaseSet struct" begin
        @test S[1]==Phase(0.0,1)
        @test_throws DomainError PhaseSet([-1.0],[1])
        @test_throws DomainError PhaseSet([-1.0],[2])
        @test S == [Phase(C[i]) for i in 1:length(C)]
        @test n_phases(S)==3
        @test [rates(S,i) for i in 1:length(C)]==C
        @test rates(S)==C
        @test [membership(S,i) for i in 1:length(m)]==m
        @test membership(S)==m
        @test phases(S)==1:length(C)
        @test N₊(S)==sum(C.>=0.0)
        @test N₋(S)==sum(C.<=0.0)
        @test DiscretisedFluidQueues.checksquare(T)===nothing
        @test_throws DomainError DiscretisedFluidQueues.checksquare(T[1:2,:])
    end

    @testset "BoundedFluidQueue struct" begin
        @test model.T==T
        @test n_phases(model)==length(C)
        @test n_phases(model.S)==length(C)
        @test model.S==S
        @test rates(model)==C
        @test rates(model,1)==C[1]
        @test n_phases(model)==length(C)
        @test phases(model)==1:length(C)

        @test_logs (:warn,"row sums of T should be 0") BoundedFluidQueue(T_warn,S,nodes[end])

        @test_throws DomainError BoundedFluidQueue(T_nz[1:2,:],S,nodes[end])
        @test_throws DomainError BoundedFluidQueue(T,S[1:end-1],nodes[end])
    end

    for f in fieldnames(BoundedFluidQueue)
        @eval @test am.$f==augment_model(model).$f
    end
end

ext: jl


text:
# An implementation of delimited continuations (the shift/reset operators) in
# Julia. Works by transforming all Julia code to continuation passing style.
# The `shift` operator then just has to return the continuation.

# https://en.wikipedia.org/wiki/Delimited_continuation
# https://en.wikipedia.org/wiki/Continuation-passing_style

using IRTools.All

struct Func
  f # Avoid over-specialising on the continuation object.
end

(f::Func)(args...) = f.f(args...)

function captures(ir, vs)
  us = Set()
  for v in vs
    isexpr(ir[v].expr) || continue
    foreach(x -> x isa Variable && push!(us, x), ir[v].expr.args)
  end
  return setdiff(us, vs)
end

rename(env, x) = x
rename(env, x::Variable) = env[x]
rename(env, x::Expr) = Expr(x.head, rename.((env,), x.args)...)
rename(env, x::Statement) = stmt(x, expr = rename(env, x.expr))

excluded = [GlobalRef(Base, :getindex)]

function continuation!(bl, ir, env, vs, ret)
  rename(x) = Main.rename(env, x)
  local v, st
  while true
    isempty(vs) && return return!(bl, rename(Expr(:call, ret, returnvalue(block(ir, 1)))))
    v = popfirst!(vs)
    st = ir[v]
    isexpr(st.expr, :call) && !(st.expr.args[1] ∈ excluded) && break
    isexpr(st.expr, :lambda) &&
      (st = stmt(st, expr = Expr(:lambda, cpslambda(st.expr.args[1]), st.expr.args[2:end]...)))
    env[v] = push!(bl, rename(st))
  end
  cs = [ret, setdiff(captures(ir, vs), [v])...]
  if isempty(vs)
    next = rename(ret)
  else
    next = push!(bl, Expr(:lambda, continuation(ir, vs, cs, v, ret), rename.(cs)...))
    next = xcall(Main, :Func, next)
  end
  ret = push!(bl, stmt(st, expr = xcall(Main, :cps, next, rename(st.expr).args...)))
  return!(bl, ret)
end

function continuation(ir, vs, cs, in, ret)
  bl = empty(ir)
  env = Dict()
  self = argument!(bl)
  env[in] = argument!(bl)
  for (i, c) in enumerate(cs)
    env[c] = pushfirst!(bl, xcall(:getindex, self, i))
  end
  continuation!(bl, ir, env, vs, ret)
end

cpslambda(ir) = cpstransform(ir, true)

function cpstransform(ir, lambda = false)
  lambda || (ir = functional(ir))
  k = argument!(ir, at = lambda ? 2 : 1)
  bl = empty(ir)
  env = Dict()
  for arg in arguments(ir)
    env[arg] = argument!(bl)
  end
  continuation!(bl, ir, env, keys(ir), k)
end

cps(k, f::Core.IntrinsicFunction, args...) = k(f(args...))
cps(k, f::IRTools.Lambda{<:Tuple{typeof(cps),Vararg{Any}}}, args...) = f(k, args...)
cps(k, ::typeof(cond), c, t, f) = c ? cps(k, t) : cps(k, f)
cps(k, ::typeof(cps), args...) = k(cps(args...))

# Speed up compilation
for f in [Broadcast.broadcasted, Broadcast.materialize]
  @eval cps(k, ::typeof($f), args...) = k($f(args...))
end

@dynamo function cps(k, args...)
  ir = IR(args...)
  ir == nothing && return :(args[1](args[2](args[3:end]...)))
  cpstransform(IR(args...))
end

# shift/reset

reset(f) = cps(identity, f)
shift(f) = error("`shift` must be called inside `reset`")
cps(k, ::typeof(shift), f) = f(k)

macro reset(ex)
  :(reset(() -> $(esc(ex))))
end

k = @reset begin
  shift(k -> k)^2
end

k(4) == 16


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module TestLogisticRegression
using ParallelAccelerator

#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc function logistic_regression(iterations::Int64)
    D = 3  # Number of dimensions
    N = 4
    # synthetic but deterministic initial values
    labels = reshape([1.,.9,.4,.8],1,N)
    points = [1. 2. 3. 4.; 4. 3. 2. 1.; 4. 5. 6. 7.]
    w = reshape(2.0.*ones(D)-1.3,1,D)

    for i in 1:iterations
       w -= ((1.0./(1.0.+exp(-labels.*(w*points))).-1.0).*labels)*points'
    end
    w
end

end

using Base.Test
println("Testing logistic regression...")
@test_approx_eq TestLogisticRegression.logistic_regression(10) [0.9461770317849691 0.8827825771001324 1.2035527971160302]
println("Done testing logistic regression...")


text:
# SPDX-License-Identifier: X11
# 2020-11-27
# Sum of gcd of Tuples (Easy) (300pt)

function main()
  K = parse(Int, readline())

  sum = 0
  for i ∈ 1:K
    for j ∈ 1:K
      for k ∈ 1:K
        sum += gcd(i,j,k)
      end
    end
  end

  println(sum)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
"""
$(TYPEDEF)
Proximal operator which implements the soft thresholding operator.

```julia
sign(x) * max(abs(x) - λ, 0)
```
See [by Zheng et. al., 2018](https://ieeexplore.ieee.org/document/8573778).
"""
struct SoftThreshold <: AbstractProximalOperator end;

@inline function (s::SoftThreshold)(x::AbstractArray, λ::T) where T <: Real
    for i in eachindex(x)
        x[i] = sign(x[i]) * max(abs(x[i]) - λ, zero(eltype(x)))
    end
    return
end

@inline function (s::SoftThreshold)(y::AbstractArray, x::AbstractArray, λ::T) where T <: Real
    @assert all(size(y) .== size(x))
    for i in eachindex(x)
        y[i] = sign(x[i]) * max(abs(x[i]) - λ, zero(eltype(x)))
    end
    return
end

"""
$(TYPEDEF)
Proximal operator which implements the hard thresholding operator.

```julia
abs(x) > sqrt(2*λ) ? x : 0
```
See [by Zheng et. al., 2018](https://ieeexplore.ieee.org/document/8573778).
"""
struct HardThreshold <: AbstractProximalOperator end;

@inline function (s::HardThreshold)(x::AbstractArray, λ::T) where T <: Real
    for i in eachindex(x)
        x[i] = abs(x[i]) > sqrt(2*λ) ? x[i] : zero(eltype(x))
    end
    return
end

@inline function (s::HardThreshold)(y::AbstractArray, x::AbstractArray, λ::T) where T <: Real
    @assert all(size(y) .== size(x))
    for i in eachindex(x)
        y[i] = abs(x[i]) > sqrt(2*λ) ? x[i] : zero(eltype(x))
    end
    return
end

"""
$(TYPEDEF)
Proximal operator which implements the (smoothly) clipped absolute deviation operator.

```julia
abs(x) > ρ ? x : sign(x) * max(abs(x) - λ, 0)
```

Where `ρ = 5λ` per default.

#Fields
$(FIELDS)

# Example

```julia
opt = ClippedAbsoluteDeviation()
opt = ClippedAbsoluteDeviation(1e-1)
```

See [by Zheng et. al., 2018](https://ieeexplore.ieee.org/document/8573778).
"""
struct ClippedAbsoluteDeviation{T} <: AbstractProximalOperator where T <: Real
    """Upper threshold"""
    ρ::T
end

ClippedAbsoluteDeviation() = ClippedAbsoluteDeviation(NaN)

function (s::ClippedAbsoluteDeviation)(x::AbstractArray, λ::T) where T <: Real
    ρ = isnan(s.ρ) ? convert(eltype(x), 5)*λ : convert(eltype(x), s.ρ)
    for i in eachindex(x)
        x[i] = abs(x[i]) > ρ ? x[i] : sign(x[i]) * max(abs(x[i]) - λ, 0)
    end
    return
end

function (s::ClippedAbsoluteDeviation)(y::AbstractArray, x::AbstractArray, λ::T) where T <: Real
    @assert all(size(y) .== size(x))
    ρ = isnan(s.ρ) ? convert(eltype(x), 5)*λ : convert(eltype(x), s.ρ)
    for i in eachindex(x)
        y[i] = abs(x[i]) > ρ ? x[i] : sign(x[i]) * max(abs(x[i]) - λ, 0)
    end
    return
end


text:
module Chen1995

using JuMP
#using GLPKMathProgInterface # was used in solve

"""
    build(model, pqr, LWH; only_lxlzwyhz = true)
    build(model, p, q, r, L, W, H; only_lxlzwyhz = true)

Add the variables and constraints of Chen 1995 3D-packing (minimize wasted
container volume, heterogeneous containers, heterogeneous boxes, free rotation
allowed) to the model.

For more information see: Chen, C. S., S. M. Lee, and Q. S. Shen. "An Analytical Model for the Container Loading Problem." European Journal of Operational Research 80, no. 1 (January 5, 1995): 68–76. https://doi.org/10.1016/0377-2217(94)00002-T.

# Arguments
- `model`: JuMP model to which the variables and contraints will be added to.
- `pqr`: three sequences of numbers (with the same length) denoting
  respectively the length/width/height of the boxes/cartons
- `LWH`: three sequences of numbers (with the same length) denoting
  respectively the length/width/height of the containers
- `only_lxlzwyhz`: if the model should use only the orientation variables
  lx, lz, wy, and hz (fast but arcane), or all of them (legible but slow).
- `no_rotation`: if true the boxes will not be allowed to rotate and
  `only_lxlzwyhz` value has no effect; default: false.
- `max_packed_volume`: if true the model will not fail if not all boxes
  cannot be packed, but maximize the box volume that can be packed, however
  there is no guarantee of minimizing wasted container space anymore
  (i.e., if all boxes can be packed in less than all containers, the
  solver can put some boxes in each container anyway, it will not worry
  about trying to pack the boxes in the most dense way).
"""
function build(model, p, q, r, L, W, H;
  only_lxlzwyhz = true,
  no_rotation = false,
  max_packed_volume = false
)
  build(model, (p, q, r), (L, W, H);
    only_lxlzwyhz = only_lxlzwyhz,
    no_rotation = no_rotation,
    max_packed_volume = max_packed_volume
  )
end
function build(model, pqr, LWH;
  only_lxlzwyhz = true,
  no_rotation = false,
  max_packed_volume = false
)
  @assert length(pqr) == 3
  @assert (isone ∘ length ∘ unique ∘ map)(length, pqr)
  p, q, r = pqr
  N = length(pqr[1])
  @assert length(LWH) == 3
  @assert (isone ∘ length ∘ unique ∘ map)(length, LWH)
  L, W, H = LWH
  m = length(LWH[1])
  container_volumes = map(prod, zip(LWH...))
  box_volumes = map(prod, zip(pqr...))
  # highest length/width/height of a container
  ML, MW, MH = map((xs)->max(xs...), LWH)
  #ML, MW, MH = (10000, 10000, 10000)

  @variables model begin
    s[1:N, 1:m], Bin
    n[1:m], Bin
    x[1:N] >= 0
    y[1:N] >= 0
    z[1:N] >= 0
    a[1:N, 1:N], Bin
    b[1:N, 1:N], Bin
    c[1:N, 1:N], Bin
    d[1:N, 1:N], Bin
    e[1:N, 1:N], Bin
    f[1:N, 1:N], Bin
  end
  if !no_rotation
    @variables model begin
      lx[1:N], Bin
      lz[1:N], Bin
      wy[1:N], Bin
      hz[1:N], Bin
    end
    if !only_lxlzwyhz
      @variables model begin
        ly[1:N], Bin
        wx[1:N], Bin
        wz[1:N], Bin
        hx[1:N], Bin
        hy[1:N], Bin
      end
    end
  end

  # disable all unused variables (the model only access positions in
  # which i < k, so we can delete all i >= k)
  for var = [:a, :b, :c, :d, :e, :f], i = 1:N, k = 1:i
    setlowerbound(model[var][i, k], 0.0)
    setupperbound(model[var][i, k], 0.0)
  end

  if max_packed_volume
    @objective(model, Max, sum(box_volumes[i] * s[i, j] for i = 1:N, j = 1:m))
  else # minimize wasted space but expects to be possible pack all boxes
    @objective(model, Min, sum(container_volumes[j] * n[j] for j = 1:m))
  end

  for k = 1:N, i = 1:(k-1) # forall i < k ∈ N
    # Constraints (1)-(6) -- no overlap among boxes inside the same container.
    if no_rotation
      @constraints model begin
        x[i] + p[i] <= x[k] + (1 - a[i, k])*ML
        x[k] + p[k] <= x[i] + (1 - b[i, k])*ML
        y[i] + q[i] <= y[k] + (1 - c[i, k])*MW
        y[k] + q[k] <= y[i] + (1 - d[i, k])*MW
        z[i] + r[i] <= z[k] + (1 - e[i, k])*MH
        z[k] + r[k] <= z[i] + (1 - f[i, k])*MH
      end
    elseif only_lxlzwyhz
      @constraints model begin
        x[i] + p[i]*lx[i] + q[i]*(lz[i] - wy[i] + hz[i]) + r[i]*(1 - lx[i] - lz[i] + wy[i] - hz[i]) <= x[k] + (1 - a[i, k])*ML
        # same as above but swapping [i] and [k], and replacing 'a' by 'b'
        x[k] + p[k]*lx[k] + q[k]*(lz[k] - wy[k] + hz[k]) + r[k]*(1 - lx[k] - lz[k] + wy[k] - hz[k]) <= x[i] + (1 - b[i, k])*ML

        y[i] + q[i]*wy[i] + p[i]*(1 - lx[i] - lz[i]) + r[i]*(lx[i] + lz[i] - wy[i]) <= y[k] + (1 - c[i, k])*MW
        # same as above but swapping [i] and [k], and replacing 'c' by 'd'
        y[k] + q[k]*wy[k] + p[k]*(1 - lx[k] - lz[k]) + r[k]*(lx[k] + lz[k] - wy[k]) <= y[i] + (1 - d[i, k])*MW

        z[i] + r[i]*hz[i] + q[i]*(1 - lz[i] - hz[i]) + p[i]*lz[i] <= z[k] + (1 - e[i, k])*MH
        # same as above but swapping [i] and [k], and replacing 'e' by 'f'
        z[k] + r[k]*hz[k] + q[k]*(1 - lz[k] - hz[k]) + p[k]*lz[k] <= z[i] + (1 - f[i, k])*MH
      end
    else
      @constraints model begin
        x[i] + p[i]*lx[i] + q[i]*wx[i] + r[i]*hx[i] <= x[k] + (1 - a[i, k])*ML
        x[k] + p[k]*lx[k] + q[k]*wx[k] + r[k]*hx[k] <= x[i] + (1 - b[i, k])*ML
        y[i] + q[i]*wy[i] + p[i]*ly[i] + r[i]*hy[i] <= y[k] + (1 - c[i, k])*MW
        y[k] + q[k]*wy[k] + p[k]*ly[k] + r[k]*hy[k] <= y[i] + (1 - d[i, k])*MW
        z[i] + r[i]*hz[i] + q[i]*wz[i] + p[i]*lz[i] <= z[k] + (1 - e[i, k])*MH
        z[k] + r[k]*hz[k] + q[k]*wz[k] + p[k]*lz[k] <= z[i] + (1 - f[i, k])*MH
      end
    end

    for j = 1:m # forall i < k ∈ N, j ∈ m
      # Constraint (7) -- boxes in the same container must have a relative
      # position to each other.
      @constraint(model,
        a[i, k] + b[i, k] + c[i, k] + d[i, k] + e[i, k] + f[i, k] >= s[i, j] + s[k, j] - 1
      )
    end
  end # end forall i < k ∈ N

  if !no_rotation && !only_lxlzwyhz
    for i = 1:N
      # Necessary if the model includes all orientation variables.
      @constraints model begin
        lx[i] + ly[i] + lz[i] == 1
        wx[i] + wy[i] + wz[i] == 1
        hx[i] + hy[i] + hz[i] == 1
        lx[i] + wx[i] + hx[i] == 1
        ly[i] + wy[i] + hy[i] == 1
        lz[i] + wz[i] + hz[i] == 1
      end
    end
  end

  # The maximization of the volume of packed boxes make no sense if we already
  # force the model to pack all boxes.
  if !max_packed_volume
    for i = 1:N
      # Constraint (8) -- every box must be inside one container.
      @constraint(model, sum(s[i, j] for j = 1:m) == 1)
    end
  end

  # Constraint (9) -- if there is a box inside a container, the container
  # must be marked as used.
  sorted_bvs = sort!(map((xs)->*(xs...), zip(pqr...)))
  for j = 1:m
    # MJ is the basic volumetric upper bound in the amount of boxes that can be
    # fitted inside the respective container.
    #MJ = 0
    MJ = 10000
    #remaining_cv = container_volumes[j]
    #while MJ < N && remaining_cv >= sorted_bvs[MJ+1]
    #  remaining_cv -= sorted_bvs[MJ+1]
    #  MJ += 1
    #end
    @constraint(model, sum(s[i, j] for i = 1:N) <= MJ*n[j])
  end

  # Constraints (10)-(12) -- the boxes must respect the container boundaries
  # (the back-left-bottom boundary is already enforced by the '>= 0' trivial
  # constraint).
  for i = 1:N, j = 1:m
    if no_rotation
      @constraints model begin
        x[i] + p[i] <= L[j] + (1 - s[i, j])*ML
        y[i] + q[i] <= W[j] + (1 - s[i, j])*MW
        z[i] + r[i] <= H[j] + (1 - s[i, j])*MH
      end
    elseif only_lxlzwyhz
      @constraints model begin
        x[i] + p[i]*lx[i] + q[i]*(lz[i] - wy[i] + hz[i]) + r[i]*(1 - lx[i] - lz[i] + wy[i] - hz[i]) <= L[j] + (1 - s[i, j])*ML
        y[i] + q[i]*wy[i] + p[i]*(1 - lx[i] - lz[i]) + r[i]*(lx[i] + lz[i] - wy[i]) <= W[j] + (1 - s[i, j])*MW
        z[i] + r[i]*hz[i] + q[i]*(1 - lz[i] - hz[i]) + p[i]*lz[i] <= H[j] + (1 - s[i, j])*MH
      end
    else
      @constraints model begin
        x[i] + p[i]*lx[i] + q[i]*wx[i] + r[i]*hx[i] <= L[j] + (1 - s[i, j])*ML
        y[i] + q[i]*wy[i] + p[i]*ly[i] + r[i]*hy[i] <= W[j] + (1 - s[i, j])*MW
        z[i] + r[i]*hz[i] + q[i]*wz[i] + p[i]*lz[i] <= H[j] + (1 - s[i, j])*MH
      end
    end
  end

  model
end

"""
    extract_attributions(model) :: Dict{Int64, Vector{Int64}}

Returns a dictionary with the index of the containers as keys and the index
of the boxes inside that container as values.
"""
function extract_attributions(model) :: Dict{Int64, Vector{Int64}}
  s = getvalue(model[:s])
  N = size(s, 1)
  m = size(s, 2)
  d = Dict{Int64, Vector{Int64}}()
  # builds the dictionary from container index to box indexes
  for i = 1:N, j = 1:m
    if isone(round(s[i, j]))
      if haskey(d, j)
        push!(d[j], i)
      else
        d[j] = [i]
      end
    end
  end

  d
end

"""
    extract_positions(model) :: NTuple{3, Vector{Float64}}

Returns three arrays (x, y, z) of coordinates where the boxes were placed.
This model has the coordinates of he boxes as continuous variables. It is
tempting to truncate the values, as this would keep the validity of the
solution and make it more compatible with other models, but as the model
can be modified to worry about the gravity center, it is best kept as it is.
"""
function extract_positions(model) :: NTuple{3, Vector{Float64}}
  map(s ->getvalue(model[s]), (:x, :y, :z))
end

"""
    extract_oriented_sizes(model, pqr; only_lxlzwyhz = true) :: NTuple{3, Vector{Int64}}

Returns the oriented sizes of the boxes in the solved model. Note that
this function should not be called if the model was built with
`no_rotation = true` as, in this case, the model has no orientation variables
and the oriented size of the boxes is simply the original box sizes.
"""
function extract_oriented_sizes(
  model, pqr; only_lxlzwyhz = true
) :: NTuple{3, Vector{Int64}}
  # computes the oriented pqr values (sizes of the boxes when oriented)
  N = length(model[:x])
  zeroes = repeat([zero(typeof(pqr[1][1]))], N)
  p′, q′, r′ = zeroes, copy(zeroes), copy(zeroes)
  p, q, r = pqr
  lx = getvalue(model[:lx])
  lz = getvalue(model[:lz])
  wy = getvalue(model[:wy])
  hz = getvalue(model[:hz])
  if !only_lxlzwyhz
    ly = getvalue(model[:ly])
    wx = getvalue(model[:wx])
    wz = getvalue(model[:wz])
    hx = getvalue(model[:hx])
    hy = getvalue(model[:hy])
  end
  if only_lxlzwyhz
    for i = 1:N
      p′[i] = p[i]*lx[i] + q[i]*(lz[i] - wy[i] + hz[i]) +
        r[i]*(1 - lx[i] - lz[i] + wy[i] - hz[i])
      q′[i] = p[i]*(1 - lx[i] - lz[i]) + q[i]*wy[i] +
        r[i]*(lx[i] + lz[i] - wy[i])
      r′[i] = p[i]*(lz[i]) + q[i]*(1 - lz[i] - hz[i]) + r[i]*hz[i]
    end
  else
    for i = 1:N
      p′[i] = p[i]*lx[i] + q[i]*wx[i] + r[i]*hx[i]
      q′[i] = p[i]*ly[i] + q[i]*wy[i] + r[i]*hy[i]
      r′[i] = p[i]*lz[i] + q[i]*wz[i] + r[i]*hz[i]
    end
  end

  (p′, q′, r′)
end

"""
    extract_solution(model, pqr; only_lxlzwyhz = true)

Returns the triple of the results of extract\\_{attributions, positions,
oriented\\_sizes} over the solved model.

The value of only_lxlzwyhz needs to be the same passed to Chen1995.build.
"""
function extract_solution(model, pqr; only_lxlzwyhz = true)
  d = extract_attributions(model)
  c = extract_positions(model)
  pqr′ = extract_oriented_sizes(model, pqr; only_lxlzwyhz = only_lxlzwyhz)

  (d, c, pqr′)
end

#= PROBLEM: should we treat the possibility of the model being unfeasible?
"""
    solve(pqr, LWH, solver = GLPKSolverMIP())

Solves the problem defined by pqr and LWH using the Chen 1995 model in the
specified solver, and returns a triple of extract_{attributions, positions,
oriented_sizes} (which is the extracted solution).
"""
function solve(pqr, LWH, solver = GLPKSolverMIP())
  model = Model(solver = solver)
  build(model, pqr, LWH)
  status = solve(model, suppress_warnings = true)
  if
  extract_solution(model, pqr)
end=#

#=function chen1995_1to7(model, pqr, LWH)

end=#

end # module Chen1995


text:
export with_latest

import Base: show

"""

    with_latest(sources...)

Creates a with_latest operator, which combines the source Observable with other
Observables to create an Observable whose values are calculated from the latest values of each, only when the source emits.

# Examples

```jldoctest
using Rocket

source = of(1)
subscribe!(source |> with_latest(from(1:5)), logger())
;

# output

[LogActor] Data: (1, 5)
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref)
"""
with_latest()                                = error("with_latest operator expects at least one inner observable on input")
with_latest(sources...)                      = with_latest(sources)
with_latest(sources::S) where { S <: Tuple } = WithLatestOperator{S}(sources)

##

struct WithLatestOperator{S} <: InferableOperator
    sources :: S
end

function on_call!(::Type{L}, ::Type{R}, operator::WithLatestOperator{S}, source::M) where { L, R, M, S }
    return WithLatestObservable{R, M, S}(source, operator.sources)
end

operator_right(operator::WithLatestOperator, ::Type{L}) where L = Tuple{ L, combined_type(operator.sources).parameters... }

Base.show(io::IO, ::WithLatestOperator) = print(io, "WithLatestOperator()")

##

struct WithLatestInnerActor{L, W, I} <: Actor{L}
    wrapper :: W
end

Base.show(io::IO, inner::WithLatestInnerActor{L, W, I}) where { L, W, I } = print(io, "WithLatestInnerActor($L, $I)")

on_next!(actor::WithLatestInnerActor{L, W, I}, data::L) where { L, W, I } = next_received!(actor.wrapper, data, Val{I}())
on_error!(actor::WithLatestInnerActor{L, W, I}, err)    where { L, W, I } = error_received!(actor.wrapper, err, Val{I}())
on_complete!(actor::WithLatestInnerActor{L, W, I})      where { L, W, I } = complete_received!(actor.wrapper, Val{I}())

##

struct WithLatestActorWrapper{S, A}
    storage :: S
    actor   :: A

    nsize   :: Int
    cstatus :: BitArray{1} # Completion status
    vstatus :: BitArray{1} # Values status

    subscriptions :: Vector{Teardown}

    WithLatestActorWrapper{S, A}(storage::S, actor::A) where { S, A } = begin
        nsize   = length(storage)
        cstatus = falses(nsize)
        vstatus = falses(nsize)
        subscriptions = fill!(Vector{Teardown}(undef, nsize), voidTeardown)
        return new(storage, actor, nsize, cstatus, vstatus, subscriptions)
    end
end

cstatus(wrapper::WithLatestActorWrapper, index) = wrapper.cstatus[index]
vstatus(wrapper::WithLatestActorWrapper, index) = wrapper.vstatus[index]

secondary_cstatus(wrapper::WithLatestActorWrapper) = view(wrapper.cstatus, 2:length(wrapper.cstatus))
secondary_vstatus(wrapper::WithLatestActorWrapper) = view(wrapper.vstatus, 2:length(wrapper.vstatus))

function dispose(wrapper::WithLatestActorWrapper)
    fill!(wrapper.cstatus, true)
    foreach(s -> unsubscribe!(s), wrapper.subscriptions)
end

function next_received!(wrapper::WithLatestActorWrapper, data, index::Val{I}) where I
    setstorage!(wrapper.storage, data, index)
    wrapper.vstatus[I] = true
    if I === 1 && all(wrapper.vstatus) && !all(wrapper.cstatus)
        next!(wrapper.actor, snapshot(wrapper.storage))
    end
end

function error_received!(wrapper::WithLatestActorWrapper, err, index::Val{I}) where I
    if !cstatus(wrapper, I)
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::WithLatestActorWrapper, ::Val{I}) where I
    if !all(wrapper.cstatus)
        wrapper.cstatus[I] = true
        if I === 1 || all(wrapper.cstatus) || wrapper.vstatus[I] === false
            dispose(wrapper)
            complete!(wrapper.actor)
        end
    end
end

##

struct WithLatestObservable{T, M, S} <: Subscribable{T}
    main      :: M
    secondary :: S
end

function on_subscribe!(observable::WithLatestObservable{T, M, S}, actor::A) where { T, M, S, A }
    storage = getmstorage(T)
    wrapper = WithLatestActorWrapper{typeof(storage), A}(storage, actor)

    for (index, source) in enumerate(observable.secondary)
        wrapper.subscriptions[index + 1] = subscribe!(source, WithLatestInnerActor{eltype(source), typeof(wrapper), index + 1}(wrapper))
        if cstatus(wrapper, index + 1) === true && vstatus(wrapper, index + 1) === false
            dispose(wrapper)
            break
        end
    end

    if !all(secondary_cstatus(wrapper)) || all(secondary_vstatus(wrapper))
        source = observable.main
        wrapper.subscriptions[1] = subscribe!(source, WithLatestInnerActor{eltype(source), typeof(wrapper), 1}(wrapper))
    end

    return WithLatestSubscription(wrapper)
end

##

struct WithLatestSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: WithLatestSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::WithLatestSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::WithLatestObservable{D}) where D = print(io, "WithLatestObservable($D)")
Base.show(io::IO, ::WithLatestSubscription)          = print(io, "WithLatestSubscription()")

##


text:
using Random
rng = MersenneTwister(1234)

import Plots
const GR = Plots.GR

GR.__init__()

x = 0:π/100:2π
y = sin.(x)
GR.plot(x, y)

x = LinRange(0, 1, 51)
y = x .- x.^2
GR.scatter(x, y)

sz = LinRange(50, 300, length(x))
c = LinRange(0, 255, length(x))
GR.scatter(x, y, sz, c)

GR.stem(x, y)

GR.histogram(randn(rng, 10000))

GR.plot(randn(rng, 50))

GR.oplot(randn(rng, 50, 3))

x = LinRange(0, 30, 1000)
y = cos.(x) .* x
z = sin.(x) .* x
GR.plot3(x, y, z)

angles = LinRange(0, 2pi, 40)
radii = LinRange(0, 2, 40)
GR.polar(angles, radii)

x = 2 .* rand(rng, 100) .- 1
y = 2 .* rand(rng, 100) .- 1
z = 2 .* rand(rng, 100) .- 1
GR.scatter3(x, y, z)

c = 999 .* rand(rng, 100) .+ 1
GR.scatter3(x, y, z, c)

x = randn(rng, 100000)
y = randn(rng, 100000)
GR.hexbin(x, y)

x = 8 .* rand(rng, 100) .- 4
y = 8 .* rand(rng, 100) .- 4
z = sin.(x) .+ cos.(y)
GR.contour(x, y, z)

x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
GR.contour(x, y, z)

x = 8 .* rand(rng, 100) .- 4
y = 8 .* rand(rng, 100) .- 4
z = sin.(x) .+ cos.(y)
GR.contourf(x, y, z)

x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
GR.contourf(x, y, z)

x = 8 .* rand(rng, 100) .- 4
y = 8 .* rand(rng, 100) .- 4
z = sin.(x) + cos.(y)
GR.tricont(x, y, z)

x = 8 .* rand(rng, 100) .- 4
y = 8 .* rand(rng, 100) .- 4
z = sin.(x) .+ cos.(y)
GR.surface(x, y, z)

x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
GR.surface(x, y, z)

x = 8 .* rand(rng, 100) .- 4
y = 8 .* rand(rng, 100) .- 4
z = sin.(x) .+ cos.(y)
GR.trisurf(x, y, z)

z = GR.peaks()
GR.surface(z)

x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
GR.wireframe(x, y, z)

x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
GR.heatmap(z)

GR.imshow(z)

s = LinRange(-1, 1, 40)
v = 1 .- (s .^ 2 .+ (s .^ 2)' .+ reshape(s,1,1,:) .^ 2) .^ 0.5
GR.isosurface(v, isovalue=0.2)

GR.GR3.terminate()

GR.volume(randn(rng, 50, 50, 50))

N = 1_000_000
x = randn(rng, N)
y = randn(rng, N)
GR.shade(x, y)

GR.setprojectiontype(0)

GR.clearws()
xd = -2 .+ 4 * rand(rng, 100)
yd = -2 .+ 4 * rand(rng, 100)
zd = [xd[i] * exp(-xd[i]^2 - yd[i]^2) for i = 1:100]

GR.setviewport(0.1, 0.95, 0.1, 0.70)
GR.setwindow(-2, 2, -2, 2)
GR.setspace(-0.5, 0.5, 0, 90)
GR.setcolormap(0)
GR.setlinecolorind(1)
GR.setmarkersize(1)
GR.setmarkertype(-1)
GR.setmarkercolorind(1)
GR.setcharheight(0.024)
GR.settextalign(2, 0)
GR.settextfontprec(3, 0)

x, y, z = GR.gridit(xd, yd, zd, 200, 200)
h = -0.5:0.05:0.5
GR.surface(x, y, z, 5)
GR.contour(x, y, h, z, 0)
GR.polymarker(xd, yd)
GR.axes(0.25, 0.25, -2, -2, 2, 2, 0.01)

GR.updatews()

Plots.gr()
x = 1:10; y = rand(rng, 10);
Plots.plot(x,y)

Plots.scatter(x,y)

Plots.histogram(randn(rng, 10000))

z = GR.peaks();
Plots.surface(z)

Plots.contour(z)

x = LinRange(0, 30, 1000); y = cos.(x) .* x; z = sin.(x) .* x
Plots.plot3d(x, y, z)

import StatsPlots
y = rand(rng, 100, 4)
StatsPlots.violin(["Series 1" "Series 2" "Series 3" "Series 4"], y, leg = false)

StatsPlots.boxplot!(["Series 1" "Series 2" "Series 3" "Series 4"], y, leg = false)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ISL"
version = v"0.22.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://isl.gforge.inria.fr/isl-0.22.tar.xz", "6c8bc56c477affecba9c59e2c9f026967ac8bad01b51bdd07916db40a517b9fa"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd isl-0.22/
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/isl-0.14.1-no-undefined.patch
autoreconf
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-gmp-prefix=${prefix}
make -j
make install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libisl", :libisl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "errors" begin

let
    ex = CuError(0)
    @test CUDAdrv.name(ex) == "SUCCESS"
    @test CUDAdrv.description(ex) == "no error"

    io = IOBuffer()
    showerror(io, ex)
    str = String(take!(io))

    @test occursin("0", str)
    @test occursin("no error", str)
end

let
    ex = CuError(0, "foobar")

    io = IOBuffer()
    showerror(io, ex)
    str = String(take!(io))

    @test occursin("foobar", str)
end

end


text:
import Test: @testset, @test, @test_throws
import StaticArrays: SVector
import Statistics: mean

using MultipleScattering
using LinearAlgebra

@testset "Fourier transform" begin

    ω_arr = collect(LinRange(0.0,100.0,200))
    t_arr = ω_to_t(ω_arr)
    N = length(ω_arr) - 1
    t0 = t_arr[N]

    # test invertablity
    @test ω_arr ≈ t_to_ω(ω_to_t(ω_arr))
    fs = rand(-1:0.01:1.,2N+1,1)
    @test fs ≈ frequency_to_time(time_to_frequency(fs,t_arr),ω_arr)

    # Gaussian to Gaussian
    as = [0.5];
    Fs = [exp(-a*ω^2)*exp(im*t0*ω) for ω in ω_arr, a in as]
    fs = frequency_to_time(Fs,ω_arr,t_arr);
    true_fs = (1/(2pi))*[sqrt(pi/a)*exp(-(t-t0)^2/(4*a)) for t in t_arr, a in as]
    @test fs ≈ true_fs
    @test Fs ≈ time_to_frequency(fs,t_arr,ω_arr)

    # sinc to rectangle
    rect(t) = (abs(t)<0.5) ? 1.0 : 0.
    as = [0.7]; #only positive values
    Fs = [ (ω == 0.0im) ? 1.0+0.0im : sin(a*pi*ω)*exp(im*t0*ω)/(a*pi*ω) for ω in ω_arr, a in as]
    fs = frequency_to_time(Fs,ω_arr,t_arr) # only correct for half of time_arr
    true_fs = (1/(2pi))*[(1/abs(a))*rect((t-t0)/(2*pi*a)) for t in t_arr, a in as]
    @test mean(abs.(fs-true_fs)) < 0.02*mean(abs.(true_fs))

    # trapezoidal integration
    ω_arr = LinRange(0.0,130.0,600)
    t_arr = ω_to_t(ω_arr)

    fs = cos.(t_arr).*exp.(-(t_arr .- t_arr[end]/2).^2/(t_arr[end]^2/25))
    fs = reshape(fs,length(fs),1)

    Fs_trap = time_to_frequency(fs, t_arr; method=:trapezoidal)
    Fs = time_to_frequency(fs, t_arr)
    # plot(ω_arr,[real(Fs),real(Fs_trap)])
    @test mean(abs.(Fs-Fs_trap)) < 0.0012*mean(abs.(Fs))

    fs_trap = frequency_to_time(Fs_trap,ω_arr; method=:trapezoidal)
    @test mean(abs.(fs-fs_trap))< 1e-5*mean(abs.(fs))

    ## rectangle to sinc
    # as = [1/(2*(maximum(ω_arr)+ω_arr[2]))]; #only positive values
    # # t0 = t_arr[end]/3
    # Fs = Complex{Float64}[ rect(a*ω)*exp(im*t0*ω) for ω in ω_arr, a in as]
    # fs =  frequency_to_time(Fs, ω_arr, t_arr)
    # true_fs = [ (t==t0) ? 1/(2pi*a) : sin((t-t0)/(2a))/(pi*(t-t0)) for t in t_arr, a in as]
    # @test maximum(abs.(fs-true_fs)) < 1e-3*maximum(abs.(true_fs))

    # pole to wavy heaviside
    # ω_arr = collect(LinRange(0.0,10.0,4000)) # need fine mesh near ω=0
    # dω = (ω_arr[2]-ω_arr[1])
    # N = length(ω_arr)-1;
    # time_arr = LinRange(0,2π/dω,2N+2)[1:(2N+1)]
    # as = [0.1]; #only positive values
    # Fs = Complex{Float64}[1/(a-im*ω) for ω in ω_arr, a in as]
    # fs =  frequency_to_time(Fs,ω_arr,time_arr);
    # true_fs = [exp(-a*t) for t in time_arr, a in as]
    # mean(abs.(fs-true_fs))
    #  1e-3*mean(abs.(true_fs))

end


text:
using SeisIO, SeisNoise, JLD2, Distributed, Dates

@everywhere include("../../src/SeisXcorrelation.jl")
#@everywhere using SeisXcorrelation
include("../../src/pairing.jl")
include("../../src/utils.jl")

# input parameters
InputDict = Dict( "finame"     => "/Users/jared/SCECintern2019/data/seismic/BPnetwork_2003/BPnetwork_2003_combined.jld2",
                  "timeunit"   => 86400,     # unit of time xcorrs are saved in. DL_time_unit if SeisDownload was used
                  "basefoname" => "/Users/jared/SCECintern2019/data/xcorrs/BPnetwork_2003/BPnetwork_2003_xcorrs",
                  "maxReadNum" => 9,         # number of processes to read from JLD2->dict and map over
                  "freqmin"    => 0.1,
                  "freqmax"    => 9.99,
                  "fs"         => 20.0,
                  "cc_len"     => 600,
                  "cc_step"    => 300,
                  "corrtype"   => ["acorr", "xchancorr", "xcorr"],
                  "corrmethod" => "cross-correlation",
                  "maxtimelag" => 100.0,
                  "to_whiten"  => false, # true, false
                  "time_norm"  => false, # false, phase, or one-bit
                  "allstack"   => false) # for now, only stacking over DL_time_unit is supported pre-save

# read data from JLD2
data = jldopen(InputDict["finame"])
# read station and time stamp lists
stlist = data["info/stationlist"]
tstamplist = data["info/DLtimestamplist"]

station_pairs = generate_pairs(stlist) # array
# sort station pairs into autocorr, xcorr, and xchancorr
sorted_pairs = sort_pairs(station_pairs) # dict

# create base output file with network, station, and pairing info
jldopen("$(InputDict["basefoname"]).jld2", "w") do file
    file["info/timestamplist"]   = tstamplist;
    file["info/stationlist"]     = stlist;
    file["info/corrstationlist"] = sorted_pairs;
end

st=time()
# get indices of timesteps to load in each iteration
mapWindows = window_read(length(tstamplist), InputDict["maxReadNum"])
for win in mapWindows
    # load timesteps into array of dicts [Dict(ts/stn=>data for stn in stlist) for ts in tslist]
    dsets = read_JLD22Dict(data, tstamplist[win])
    # map seisxcorrelation over timesteps
    pmap((x,y)->seisxcorrelation(x, y, InputDict), dsets, tstamplist[win])
end
et=time()
close(data)
println("Successfully completed and saved cross-correlations in $(et-st) seconds.")
rmprocs(workers())

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
using Color, Images, FixedPointNumbers

ima = imread("grayscale_image_color.png")
imb = convert(Image{Gray{Ufixed8}}, ima)
imwrite(imb, "grayscale_image_julia.png")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Pkg; Pkg.activate("../")
using WebAssemblyText

function translate(infile, outfile)
    wat = jl2wat(infile)
    open(outfile, "w") do io
        write(io, wat)
    end
end

translate(ARGS[1], ARGS[2])


text:
L = 50
D = UncertainDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])
V = UncertainValueDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])
I = UncertainIndexDataset([UncertainValue(Normal, rand(), rand()) for i = 1:L])

@test randomshuffle(D) isa Vector{Float64}
@test randomshuffle(V) isa Vector{Float64}
@test randomshuffle(I) isa Vector{Float64}

@test randomphases(D) isa Vector{Float64}
@test randomphases(V) isa Vector{Float64}
@test randomphases(I) isa Vector{Float64}

@test randomamplitudes(D) isa Vector{Float64}
@test randomamplitudes(V) isa Vector{Float64}
@test randomamplitudes(I) isa Vector{Float64}

@test aaft(D) isa Vector{Float64}
@test aaft(V) isa Vector{Float64}
@test aaft(I) isa Vector{Float64}

@test iaaft(D) isa Vector{Float64}
@test iaaft(V) isa Vector{Float64}
@test iaaft(I) isa Vector{Float64}

method_AAFT = AAFT()
method_IAAFT = IAAFT()
method_RandomShuffle = RandomShuffle()
method_RandomFourier = RandomFourier()

@test surrogate(D, method_AAFT) isa Vector{Float64}
@test surrogate(D, method_IAAFT) isa Vector{Float64}
@test surrogate(D, method_RandomShuffle) isa Vector{Float64}
@test surrogate(D, method_RandomFourier) isa Vector{Float64}

@test surrogate(V, method_AAFT) isa Vector{Float64}
@test surrogate(V, method_IAAFT) isa Vector{Float64}
@test surrogate(V, method_RandomShuffle) isa Vector{Float64}
@test surrogate(V, method_RandomFourier) isa Vector{Float64}

@test surrogate(I, method_AAFT) isa Vector{Float64}
@test surrogate(I, method_IAAFT) isa Vector{Float64}
@test surrogate(I, method_RandomShuffle) isa Vector{Float64}
@test surrogate(I, method_RandomFourier) isa Vector{Float64}


text:
"""
    NCBIInfluenzaData

Package for parsing and filtering NCBI's bulk influenza data. This can be used to
extract and work with tens of thousands of sequences in bulk.

It has been confirmed to work with the NCBI data of 2020-10-13
"""
module NCBIInfluenzaData

using BioSequences: BioSequences, BioSequence, LongDNASeq, @dna_str,
    AA_Term, reverse_complement
using CodecZlib: GzipCompressorStream, GzipDecompressorStream
using ErrorTypes: Option, none, some, unwrap, unwrap_or, @unwrap_or, is_error, @?
using FASTX: FASTA
using Influenza: Influenza, Protein, Segment, SeroType, Segments, Proteins, source
import Downloads

import Influenza: ReferenceProtein

imap(f) = x -> Iterators.map(f, x)
ifilter(f) = x -> Iterators.filter(f, x)

function try_split!(
    v::Vector{SubString{String}},
    s::Union{String, SubString{String}},
    sep::UInt8
)::Option{UInt}
    n = UInt(0)
    start = 1
    isempty(v) && return none
    @inbounds for i in 1:ncodeunits(s)
        if codeunit(s, i) == sep
            n += 1
            n >= length(v) && return none
            substr = SubString(s, start, i-1)
            @inbounds v[n] = substr
            start = i + 1
        end
    end
    @inbounds v[n+1] = SubString(s, start, lastindex(s))
    some(n + 1)
end

function try_split!(
    v::Vector{SubString{String}},
    s::Union{String, SubString{String}},
    sep::Char
)::Option{UInt}
    try_split!(v, s, UInt8(sep))
end

# This is pretty primitive, but fuck it
function maybe_gzip(f, path::AbstractString)
    if endswith(path, ".gz")
        io = GzipDecompressorStream(open(path))
        try
            f(io)
        finally
            close(io)
        end
    else
        open(f, path)
    end
end

"""
    download_influenza_data(dst::AbstractString, force::Bool=false)

Download data from the NCBI FTP server to an existing directory `dst`.
If `force`, re-download, even if the files already exist in the directory.
"""
function download_influenza_data(dst_dir::AbstractString; force::Bool=false)
    isdir(dst_dir) || error("Path is not a directory: \"$dst_dir\"")
    ftp_address = "https://ftp.ncbi.nih.gov/genomes/INFLUENZA/"
    for filename in [
        "genomeset.dat.gz",
        "influenza.dat.gz",
        "influenza.fna.gz",
        "influenza_aa.dat.gz",
        "ANNOTATION/blastDB.fasta"
        ]
        dstpath = joinpath(dst_dir, basename(filename))
        if !isfile(dstpath) || force
            Downloads.download(joinpath(ftp_address, filename), dstpath)
        end
    end
end

baremodule Hosts
import Base: @enum
@enum Host::UInt8 human swine avian other
export Host
end # baremodule
using .Hosts

function Base.parse(::Type{Host}, s::Union{String, SubString})
    lcase = lowercase(s)
    lcase == "human" && return Hosts.human
    lcase == "swine" && return Hosts.swine
    lcase == "avian" && return Hosts.avian
    return Hosts.other
end

"""
    IncompleteSegmentData

Placeholder struct used to incrementally build `SegmentData`. You can construct
`SegmentData` using `SegmentData(::IncompleteSegmentData)`, after which the created
struct takes ownership of all fields of the incomplete data.
"""
mutable struct IncompleteSegmentData
    id::String
    host::Host
    segment::Segment
    serotype::Option{SeroType}
    year::Int16
    isolate::String
    proteins::Vector{ReferenceProtein}
    seq::Option{LongDNASeq}
end

"""
    SegmentData

Data structure representing an influenza segment from NCBI, and relevant information
that can be extracted from the files provided by NCBI.

Some fields in some records of input files may be missing. If so, the record
will be skipped, and not parsed to a `SegmentData`.

Notes:
* The SeroType is `none(SeroType)` if the segment is not Influenza A.
* The isolate is the isolate name, e.g. "A/Denmark/1/2021".
"""
struct SegmentData
    id::String
    host::Host
    segment::Segment
    serotype::Option{SeroType}
    year::Int16
    isolate::String
    proteins::Vector{ReferenceProtein}
    seq::LongDNASeq
end

# What we need to check here is the presence of proteins and sequence, since these
# are added separately to the incomplete records. The rest of the fields are checked
# when instantiating the IncompleteSegmentData. The biological plausibility of the
# data is not addressed at instantiation, only whether it's parsable.
function SegmentData(x::IncompleteSegmentData)
    isempty(x.proteins) && error("Cannot construct SegmentData from empty proteins")
    SegmentData(x.id, x.host, x.segment, x.serotype, x.year, x.isolate, x.proteins, unwrap(x.seq))
end

function Influenza.Reference(x::SegmentData)
    Influenza.Reference(x.id, x.segment, copy(x.seq), copy(x.proteins))
end

include("parse_genomeset.jl")
include("parse_fasta.jl")
include("parse_orfs.jl")
include("filter.jl")

"""
    parse_ncbi_records(genomeset, fasta, influenza_aa_dat, influenza_dat, filter_termini)

Parse NCBI records from the paths of four files: The cleaned "genomeset.dat.gz"
output of `clean_genomeset`, "influenza.fna.gz", "influenza_aa.dat.gz" and
"influenza.dat.gz".
Return a `Dict{String, SegmentData}` with keys being the accession number.

If `filter_termini`, remove all seqs without proper termini, e.g. if only the coding
sequence is provided. This is set to `false` by default, as most submitted sequences
discard the part of the sequence covered by the primers.

Note that minimal filtering is done on the records. You may want to filter for
`isok_all(::Record)` to keep only those that pass all checks.
"""
function parse_ncbi_records(
    genomeset::String,
    fasta::String,
    influenza_aa_dat::String,
    influenza_dat::String,
    filter_termini::Bool,
)
    dict = open(parse_cleaned_genomeset, GzipDecompressorStream, genomeset)
    add_sequences!(dict, fasta)
    add_orfs!(dict, influenza_aa_dat, influenza_dat)

    # Now filter away any segments without seqs or orfs
    filter!(dict) do (_, v)
        !is_error(v.seq) && !isempty(v.proteins)
    end

    # This is a fairly slow operation, but necessary
    strip_false_termini!(dict, filter_termini)

    return Dict(k => SegmentData(v) for (k, v) in dict)
end

"""
    run_all(dir::AbstractString filter_termini=false)

Convenience function: Downloads influenza data to `dir` if not already present,
then cleans the genomeset, then parses the records, filters them.

Returns a dict of `identifier => SegmentData`
"""
function run_all(dir::AbstractString, filter_termini::Bool=false)
    println("Downloading...")
    @time download_influenza_data(dir)

    println("Cleaning genomeset.dat.gz...")
    cleaned_path = "$dir/genomeset.clean.dat.gz"
    @time if !isfile(cleaned_path)
        open(GzipDecompressorStream, "$dir/genomeset.dat.gz") do inio
            open(GzipCompressorStream, cleaned_path, "w") do outio
                clean_genomeset(outio, inio)
            end
        end
    end

    println("Parsing records...")
    @time data = parse_ncbi_records(
        "$dir/genomeset.clean.dat.gz",
        "$dir/influenza.fna.gz",
        "$dir/influenza_aa.dat.gz",
        "$dir/influenza.dat.gz",
        filter_termini
    )

    println("Filtering records...")
    @time filter!(data) do (k, v)
        isok_all(v)
    end

    println("Done!")
    return data
end

end # module


text:
module NonconvexSearch

export MTS, MTSAlg, LS1Alg, MTSOptions, LS1Options

using Reexport, Parameters
using Random: randperm
@reexport using NonconvexCore
using NonconvexCore: @params, AbstractOptimizer, AbstractResult
using NonconvexCore: AbstractModel, VecModel, debugging, getdim
using NonconvexCore: clamp_and_evaluate!
import NonconvexCore: optimize!, Workspace

include("mts.jl")

end


text:
########################
## plotting utilities ##
########################

## process dates to strings
##-------------------------

function dat2str(dats::Array{Date, 1})
    nObs = size(dats, 1)
    datsAsStr = Array(ASCIIString, nObs)
    for ii=1:nObs
        datsAsStr[ii] = string(dats[ii])
    end
    return datsAsStr
end

function dat2str(dats::Array{DateTime, 1})
    nObs = size(dats, 1)
    datsAsStr = Array(ASCIIString, nObs)
    for ii=1:nObs
        datsAsStr[ii] = string(dats[ii])
    end
    return datsAsStr
end

function dat2str(dats::Array{Integer, 1})
    nObs = size(dats, 1)
    datsAsStr = Array(ASCIIString, nObs)
    for ii=1:nObs
        datsAsStr[ii] = string(dats[ii])
    end
    return datsAsStr
end

function dat2str(dats::Array{Float64, 1})
    nObs = size(dats, 1)
    datsAsStr = Array(ASCIIString, nObs)
    for ii=1:nObs
        datsAsStr[ii] = string(dats[ii])
    end
    return datsAsStr
end

function dat2str(tm::AbstractTimedata)
    return dat2str(idx(tm))
end

## process dates to floating point numbers
##----------------------------------------

function dat2num(dats::Array{Date, 1})
    nObs = size(dats, 1)
    datsAsFloat = Array(Float64, nObs)
    for ii=1:nObs
        datsAsFloat[ii] = Dates.dayofyear(dats[ii])./366 .+ Dates.year(dats[ii])
    end
    return datsAsFloat
end

function dat2num(dats::Array{DateTime, 1})
    nObs = size(dats, 1)
    datsAsFloat = Array(Float64, nObs)
    for ii=1:nObs
        datsAsFloat[ii] = Dates.year(dats[ii]) .+
        Dates.dayofyear(dats[ii])./366 .+
        Dates.hour(dats[ii])./(366*24) .+
        Dates.minute(dats[ii])./(366*24*60) .+
        Dates.second(dats[ii])./(366*24*60*60)
    end
    return datsAsFloat
end

function dat2num(dats::Array{Float64, 1})
    return dats
end

function dat2num(dats::Array{Integer, 1})
    return dats
end

function dat2num(tm::AbstractTimedata)
    return dat2num(idx(tm))
end

function dat2num(da::DataArray)
    datsAsFloat = Array(Float64, length(da))
    datsAsFloat[!da.na] = dat2num(da.data[!da.na])
    datsAsFloat[da.na] = NaN
    return datsAsFloat
end

## remove rows with NAs only
##--------------------------

function narowrm(tm::Timematr)
    return tm # nothing to do
end

function narowrm(tn::AbstractTimedata)
    nObs, nVars = size(tn)
    onlyNAs = trues(nObs)
    for ii=1:nVars
        if isa(tn.vals.columns[ii], DataArray)
            onlyNAs = onlyNAs & tn.vals.columns[ii].na
        else
            onlyNAs = falses(nObs)
            break
        end
    end
    return tn[!onlyNAs, :]
end

#########################
## unit test utilities ##
#########################

function testcase(timedataType, number::Int)
    ## easily create TimeData objects for testing

    if number == 1
        dats = [Date(2010, 1, 1):Date(2010, 1, 4)]
        df = DataFrame()
        df[:prices1] = @data([100, 120, 110, 170])
        df[:prices2] = @data([110, 120, 100, 130])
        testInst = timedataType(df, dats)

    elseif number == 2
        dats = [Date(2010, 1, 1):Date(2010, 1, 5)]
        df = DataFrame()
        df[:prices1] = @data([NA, 120, 140, 170, 200])
        df[:prices2] = @data([110, 120, NA, 130, NA])
        testInst = timedataType(df, dats)

    elseif number == 3
        dats = [Date(2010, 1, 1):Date(2010, 1, 3)]
        df = DataFrame()
        df[:color] = @data([NA, "green", "red"])
        df[:value] = @data([110, 120, NA])
        testInst = timedataType(df, dats)
    elseif number == 4
        dats = [Date(2010, 1, 1):Date(2010, 1, 5)]
        df = DataFrame()
        df[:prices1] = @data([100, 120, 140, 170, 200])
        df[:prices2] = @data([110, 120, NA, 130, 150])
        testInst = timedataType(df, dats)
    elseif number == 5
        dats = [Date(2010, 1, 1):Date(2010, 1, 5)]
        df = DataFrame()
        df[:prices1] = @data([100, 120, NA, 170, 200])
        df[:prices2] = @data([110, 120, NA, 130, 150])
        testInst = timedataType(df, dats)
    else
        error("testInst number not implemented yet")
    end
    return deepcopy(testInst)
end

#########################
## DataArray utilities ##
#########################

function findstub_vector(x::Vector)
    ## return first value that is not NA
    if !isna(x[1])
        return x[1]
    else
        for ii=2:length(x)
            if !isna(x[ii])
                return x[ii]
            end
        end
    end
    return NA
end

function anyToDa(x::Vector)
    nObs = size(x, 1)
    stubVal = findstub_vector(x)
    if isna(stubVal) # only NAs
        return DataArray(NAtype, nObs)
    else
        # split up in values and NAs
        naIndicator = falses(nObs)
        data = Array(Any, nObs)
        for ii=1:nObs
            if isna(x[ii])
                naIndicator[ii] = true
                data[ii] = stubVal
            else
                data[ii] = x[ii]
            end
        end
    end
    return DataArray([data...], naIndicator)
end

function anyToDf(x::Matrix, nams::Array{Symbol, 1})
    nObs, nCols = size(x)
    if length(nams) != nCols
        error("number of names and columns must match")
    end
    df = DataFrame()
    for ii=1:nCols
        df[nams[ii]] = anyToDa(x[:, ii])
    end
    return df
end

## applying array function to DataArray
##-------------------------------------

function arrFunc2Da(f::Function, da::DataArray)
    ## apply a function to all no-NA elements of DataArray and return
    ## result as DataArray
    ## Note: function must preserve dimensions: n values need to be
    ## mapped to n values again!
    res = f(da.data[!da.na])
    resDa = DataArray(eltype(res), length(da))
    resDa[!da.na] = res
    return resDa
end

function arrFunc2Da(f::Function, da::DataArray, replaceNA = NA)
    ## apply a function to all no-NA elements of DataArray and return
    ## result as array
    ## Note: one must explicitly deal with NAs!
    ## Note: function must preserve dimensions: n values need to be
    ## mapped to n values again!
    res = f(da.data[!da.na])
    resArr = Array(eltype(res), length(da))
    resArr[!da.na] = replaceNA
    return resArr
end


text:
function decompose_hadamard()

    println(">>>>> Hadamard Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 3,
    "elementary_gates" => ["U3_1", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.get_full_sized_gate("H_1", 2),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => [0, π/2],
    "U3_ϕ_discretization" => [0, π/2],
    "U3_λ_discretization" => [0, π],
    )

end

function decompose_controlled_Z()

    println(">>>>> Controlled-Z Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 4,
    "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.CZGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => -π:π/2:π,
    "U3_ϕ_discretization" => -π:π/2:π,
    "U3_λ_discretization" => -π:π/2:π,
    )

end

function decompose_controlled_V()

    println(">>>>> Controlled-V Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 7,
    "elementary_gates" => ["H_1", "H_2", "T_1", "T_2", "Tdagger_1", "Tdagger_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCOpt.CVGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",
    )

end

function decompose_controlled_H()

    println(">>>>> Controlled-H Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "CNot_2_1", "Identity"],
    "target_gate" => QCOpt.CHGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "U3_θ_discretization" => -2*π:π/4:2*π,
    "U3_ϕ_discretization" => [0],
    "U3_λ_discretization" => [0],
    )

end

function decompose_controlled_H_with_R()

    println(">>>>> Controlled-H with R Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["RY_1", "RY_2", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.CHGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",

    "RY_discretization" => [-π/4, π/4, π/2, -π/2],
    "set_cnot_lower_bound" => 1,
    "set_cnot_upper_bound" => 2,
    )

end

function decompose_magic()

    println(">>>>> M Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 4,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π/2:π,
        "U3_λ_discretization" => -π:π/2:π,
        )

end

function decompose_magic_using_CNOT_1_2()

    println(">>>>> M Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π:π,
        "U3_λ_discretization" => -π:π/2:π,
        )

end

function decompose_magic_using_SHCnot()

    println(">>>>> M gate using S, H and CNOT Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["S_1", "S_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.MGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_S()

    println(">>>>> S Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 3,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.get_full_sized_gate("S_1", 2),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => [0, π/4, π/2],
        "U3_ϕ_discretization" => [-π/2, 0, π/2],
        "U3_λ_discretization" => [0, π],
        )
end

function decompose_revcnot()

    println(">>>>> Reverse CNOT Gate <<<<<")

    return Dict{String, Any}(

    "num_qubits" => 2,
    "maximum_depth" => 5,
    "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],
    "target_gate" => QCOpt.CNotRevGate(),
    "objective" => "minimize_depth",
    "decomposition_type" => "exact",
    )

end

function decompose_revcnot_with_U()

    println(">>>>> Reverse CNOT using U3 and CNOT Gates <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["U3_1", "U3_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.CNotRevGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/4:π,
        "U3_ϕ_discretization" => [0],
        "U3_λ_discretization" => [0],
        )

end

function decompose_swap()

    println(">>>>> SWAP Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.SwapGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_W()

    println(">>>>> W Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 5,
        "elementary_gates" => ["U3_2", "CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.WGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => [-π/4, 0, π/4],
        "U3_ϕ_discretization" => [0, π/4],
        "U3_λ_discretization" => [0, π/2],
        )

end

function decompose_W_using_HCnot()

    println(">>>>> W Gate using H and CNOT gates <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 6,
        "elementary_gates" => ["CH_1_2", "CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.WGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_HCoinGate()

    println(">>>>> Hadamard Coin gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["Y_1", "Y_2", "Z_1", "Z_2", "T_2", "Tdagger_1", "Sdagger_1", "SX_1", "SXdagger_2", "CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.HCoinGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_GroverDiffusion_using_HX()

    println(">>>>> Grover's Diffusion Operator <<<<<")

    return Dict{String, Any}(

        #Reference: https://arxiv.org/pdf/1804.03719.pdf (Fig 6)

        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["X_1", "X_1xX_2", "H_1xH_2", "X_2", "H_1", "H_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.GroverDiffusionGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_GroverDiffusion_using_Clifford()

    println(">>>>> Grover's Diffusion Operator <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 6,
        "elementary_gates" => ["X_1", "X_2", "H_1", "H_2", "S_1", "S_2", "T_1", "T_2", "Y_1", "Y_2", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.GroverDiffusionGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_GroverDiffusion_using_U3()

    println(">>>>> Grover's Diffusion Operator using U3 gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["U3_1", "CNot_1_2", "Identity"],
        "target_gate" => QCOpt.GroverDiffusionGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",

        "U3_θ_discretization" => -π:π/2:π,
        "U3_ϕ_discretization" => -π:π/2:π,
        "U3_λ_discretization" => [0],
        )

end

function decompose_iSwap()

    println(">>>>> iSwap Gate <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["T_1", "T_2", "Tdagger_1", "Tdagger_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
        # "elementary_gates" => ["S_1", "S_2", "Sdagger_1", "Sdagger_2", "H_1", "H_2", "CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.iSwapGate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end

function decompose_qft2_using_R()
    println(">>>>> QFT2 Gate using R and CNOT gates <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 10,
        "elementary_gates" => ["RX_1", "RZ_1", "RZ_2", "CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.QFT2Gate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "approximate",

        "RX_discretization" => [π/2],
        "RZ_discretization" => [-π/4, π/2, 3*π/4, 7*π/4],
        )

end

function decompose_qft2_using_HT()
    println(">>>>> QFT2 Gate using H, T, CNOT gates <<<<<")

    return Dict{String, Any}(

        "num_qubits" => 2,
        "maximum_depth" => 8,
        "elementary_gates" => ["H_1", "H_2", "T_1", "T_2", "Tdagger_1", "Tdagger_2", "CNot_1_2", "CNot_2_1", "Identity"],
        "target_gate" => QCOpt.QFT2Gate(),
        "objective" => "minimize_depth",
        "decomposition_type" => "exact",
        )

end


text:
"""
optimize(mins::Array, maxs::Array, nutamounts::Array, calories::Array)

This is also called the Stigler's Diet Problem

# Parameters

* mins: minimum allowed intake of nutrients for a given age and sex per day
* maxs: maximum allowed intake of nutrients for a given age and sex per day
* nutamounts: A nested list where each food (the outer list) has a combination of different nutrients (the inner list)
* calories: A list of amount of energy (kcal) in each food
"""
function optimize(mins, maxs, nutamounts, calories, dri_ids, nfoods, exclude_indices)

	optimizer = Juniper.Optimizer
	params = Dict{Symbol,Any}()
	params[:nl_solver] = with_optimizer(Ipopt.Optimizer, print_level=0)
	# params[:mip_solver] = with_optimizer(Cbc.Optimizer, logLevel=0)
	m = Model(with_optimizer(optimizer, params))

	nnutrients = length(dri_ids)

	@variable(m, a[i=1:nfoods,j=1:nnutrients] == nutamounts[i][j]);  # amount of nutrient Nj in one unit of food Fi.
	@variable(m, y[1:nfoods]) # amount of food i per day. the amounts are per 100 grams per day
	@variable(m, cc[c=1:nfoods] == calories[c]); # calorie of each food

	@constraint(m, y[[i for i in 1:nfoods if !in(i, exclude_indices)]] .>= 0)  # all food amounts except those in the exclude_indices can have a value more than 0.
	@constraint(m, y[exclude_indices] .== 0.0)  # set those foods that are excluded to 0
	for nn in 1:nnutrients
		@constraint(m, mins[nn] <= sum(a[i,nn]*y[i] for i=1:nfoods) <= maxs[nn])
	end

	@objective(m, Min, sum(cc[i]*y[i] for i=1:nfoods));

	optimize!(m)
	status = termination_status(m)
	primal_status(m)
	dual_status(m)
	objval = JuMP.objective_value(m)
	# JuMP.value.(a)
	yval = JuMP.value.(y);

	return m, status, objval, yval
end


text:
# This file was generated, do not modify it. # hide
ẑ[:lambda] = 5.0;


text:
using XRayQuadrics

struct Particle
    r0::Vector{Float64}     # initial position
    v::Vector{Float64}      # propagation direction (unit)
    E::Float64              # energy
    id::Int64               # id number for tracking
    Particle(r0,v,E) = begin
        if length(r0) != 3 || length(v) != 3
            error("3-vectors please")
        else
            new(r0,v/norm(v),E)
        end
    end
end

struct Attenuator   # UNUSED
    density
    toppoint
    bottompoint
    normal
    massattenuation
end

struct PixelatedAttenuator
    holes
    largeradius # largest hole radius in attenuator
    density     # in kg/m^3
    toppoint    # a point (x,y,z) on the one side of the attenuator (in m)
    bottompoint # a point on the other side of the attenuator
    normal      # normal vector to attenuator surface
    massattenuation # table of mass attenuation coefficients per energy
end

struct Mirror
    tq::TruncatedQuadric
    reflection_probability::Matrix{Float64}
    matl::String
end

struct WolterIShell
    hyp::Mirror
    par::Mirror
end

struct Optic
    shells::Vector{WolterIShell}
    blocks::Vector{Attenuator}
end

# function Particle(r0::Vector{Real}, v::Vector{Real}, E::Real)
#     if length(r0) != 3 || length(v) != 3
#         error("3-vectors, please")
#     end

#     return Particle(r0, v/norm(v), E)
# end


text:
module FilePathsBase

using Dates
using Mmap
using Printf
using UUIDs

import Base: ==
export
    # Types
    AbstractPath,
    Path,
    SystemPath,
    PosixPath,
    WindowsPath,
    Mode,
    Status,
    FileBuffer,

    # Methods
    cwd,
    home,
    hasparent,
    parents,
    isascendant,
    isdescendant,
    filename,
    extension,
    extensions,
    exists,
    absolute,
    isabsolute,
    mode,
    created,
    modified,
    normalize,
    canonicalize,
    relative,
    isrelative,
    ismount,
    islink,
    cp,
    mv,
    sync,
    tmpname,
    tmpdir,
    mktmp,
    mktmpdir,
    chown,
    executable,
    readable,
    writable,
    raw,
    readpath,
    walkpath,
    diskusage,

    # Macros
    @p_str,
    @__PATH__,
    @__FILEPATH__,

    # Constants
    READ,
    WRITE,
    EXEC

export isexecutable

const PATH_TYPES = Type[]

function __init__()
    # Register the default fallback path type based on the os.
    register(Sys.iswindows() ? WindowsPath : PosixPath)
end

"""
    AbstractPath

Defines an abstract filesystem path.

# Properties

- `segments::Tuple{Vararg{String}}` - path segments (required)
- `root::String` - path root (defaults to "/")
- `drive::String` - path drive (defaults to "")
- `separator::String` - path separator (defaults to "/")

# Required Methods
- `tryparse(::Type{T}, str::String)` - For parsing string representations of your path
- `read(path::T)`
- `write(path::T, data)`
- `exists(path::T` - whether the path exists
- `stat(path::T)` - File status describing permissions, size and creation/modified times
- `mkdir(path::T; kwargs...)` - Create a new directory
- `rm(path::T; kwags...)` - Remove a file or directory
- `readdir(path::T)` - Scan all files and directories at a specific path level
"""
abstract type AbstractPath end  # Define the AbstractPath here to avoid circular include dependencies

"""
    register(::Type{<:AbstractPath})

Registers a new path type to support using `Path("...")` constructor and p"..." string
macro.
"""
function register(T::Type{<:AbstractPath})
    # We add the type to the beginning of our PATH_TYPES,
    # so that they can take precedence over the Posix and
    # Windows paths.
    pushfirst!(PATH_TYPES, T)
end

"""
    ispathtype(::Type{T}, x::AbstractString) where T <: AbstractPath

Return a boolean as to whether the string `x` fits the specified the path type.
"""
function ispathtype end

include("constants.jl")
include("utils.jl")
include("libc.jl")
include("mode.jl")
include("status.jl")
include("buffer.jl")
include("path.jl")
include("aliases.jl")
include("system.jl")
include("posix.jl")
include("windows.jl")
include("test.jl")
include("deprecates.jl")

end # end of module


text:
using DataFrames
using Distributions
using Gillespie
using BenchmarkTools

function sir(beta,gamma,N,S0,I0,R0,tf)
    t = 0
    S = S0
    I = I0
    R = R0
    ta=Vector{Float64}(undef,0)
    Sa=Vector{Float64}(undef,0)
    Ia=Vector{Float64}(undef,0)
    Ra=Vector{Float64}(undef,0)
    while t < tf
        push!(ta,t)
        push!(Sa,S)
        push!(Ia,I)
        push!(Ra,R)
        pf1 = beta*S*I
        pf2 = gamma*I
        pf = pf1+pf2
        dt = rand(Exponential(1/pf))
        t = t+dt
        if t>tf
            break
        end
        ru = rand()
        if ru<(pf1/pf)
            S=S-1
            I=I+1
        else
            I=I-1
            R=R+1
        end
    end
    results = DataFrame()
    results[!,:t] = ta
    results[!,:S] = Sa
    results[!,:I] = Ia
    results[!,:R] = Ra
    return(results)
end

function F(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = [9999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/10000.0,0.05]
tf = 1000.0

ssa(x0,F,nu,parms,tf) # compile
Random.seed!(1234)
@benchmark ssa($x0,$F,$nu,$parms,$tf) samples=1000 seconds=100

function F2(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = SA[9999,1,0]
nu = SA[-1 1 0
         0 -1 1]
parms = SA[0.1/10000.0,0.05]
tf = 1000.0

ssa(x0,F2,nu,parms,tf) # compile
Random.seed!(1234)
@benchmark ssa($x0,$F2,$nu,$parms,$tf) samples=1000 seconds=100

sir(0.1/10000,0.05,10000,9999,1,0,1000) # compile
Random.seed!(1234)
@benchmark sir(0.1/10000,0.05,10000,9999,1,0,1000) samples=1000 seconds=100

using DiffEqBiological

sir_model = @reaction_network rn begin
  0.1/10000.0, s + i --> 2i
  0.05, i --> r
end
sir_prob = DiscreteProblem([9999,1,0],(0.0,tf))
sir_jump_prob = JumpProblem(sir_prob,Direct(),sir_model)
sir_sol = solve(sir_jump_prob,SSAStepper()) # compile
@benchmark solve(sir_jump_prob,SSAStepper()) samples=1000 seconds=100

sir_prob = DiscreteProblem(SA[9999,1,0],(0.0,tf))
sir_jump_prob = JumpProblem(sir_prob,Direct(),sir_model)
sir_sol = solve(sir_jump_prob,SSAStepper()) # compile
@benchmark solve(sir_jump_prob,SSAStepper()) samples=1000 seconds=100


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
@recipe function f(
    marg::MarginalDist,
    vsel::NTuple{2,Union{Symbol, Expr, Integer}};
    intervals = default_credibilities,
    colors = default_colors,
    smoothing = 0,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict(),
    interval_labels = [],
    normalize = true
)
    _plots_module() != nothing || throw(ErrorException("Package Plots not available, but required for this operation"))
    hist = convert(Histogram, marg.dist)
    seriestype = get(plotattributes, :seriestype, :histogram2d)

    xlabel = get(plotattributes, :xguide, "x$(vsel[1])")
    ylabel = get(plotattributes, :yguide, "x$(vsel[2])")

    # histogram / heatmap
    if seriestype == :histogram2d || seriestype == :histogram || seriestype == :hist
        @series begin
            seriestype := :bins2d
            xguide --> xlabel
            yguide --> ylabel
            colorbar --> true

            hist.edges[1], hist.edges[2], _plots_module().Surface(hist.weights)
        end

    # smallest interval contours
    elseif seriestype == :smallest_intervals_contour || seriestype == :smallest_intervals_contourf

        colors = colors[sortperm(intervals, rev=true)]

        if seriestype == :smallest_intervals_contour
            plotstyle = :contour
        else
            plotstyle = :contourf
        end

        lev = calculate_levels(hist, intervals)
        x, y = get_bin_centers(marg)
        m = hist.weights

        if smoothing != 0
            ker = gaussian_kernel(smoothing)
            m = convolution(m, ker, padding=:same)
        end

        # quick fix: needed when plotting contour on top of histogram
        # otherwise scaling of histogram colorbar would change scaling
        lev = lev/10000
        m = m/10000

        colorbar --> false
        xguide --> xlabel
        yguide --> ylabel

        if _plots_module().backend() == _plots_module().PyPlotBackend()
            @series begin
                seriestype := plotstyle
                levels --> lev
                linewidth --> 2
                seriescolor --> colors # currently only works with pyplot
                (x, y, m')
            end
        else
             @series begin
                seriestype := plotstyle
                levels --> lev
                linewidth --> 2
                (x, y, m')
            end
        end

    # smallest intervals heatmap
    elseif seriestype == :smallest_intervals
        colors = colors[sortperm(intervals, rev=true)]

        hists, realintervals = get_smallest_intervals(hist, intervals)

        for (i, int) in enumerate(realintervals)
            @series begin
                seriestype := :bins2d
                seriescolor --> _plots_module().cgrad([colors[i], colors[i]])
                xguide --> xlabel
                yguide --> ylabel

                hists[i].edges[1], hists[i].edges[2], _plots_module().Surface(hists[i].weights)
            end

            # fake a legend
            interval_label = isempty(interval_labels) ? "smallest $(@sprintf("%.2f", realintervals[i]*100))% interval(s)" : interval_labels[i]

            @series begin
                seriestype := :shape
                fillcolor --> colors[i]
                linewidth --> 0
                label --> interval_label
                colorbar --> false
                [hists[i].edges[1][1], hists[i].edges[1][1]], [hists[i].edges[2][1], hists[i].edges[2][1]]
            end
        end

    # marginal histograms
    elseif seriestype == :marginal

        layout --> _plots_module().grid(2,2, widths=(0.8, 0.2), heights=(0.2, 0.8))
        link --> :both

        if get(diagonal, "seriestype", :histogram) != :histogram
            colorbar --> false
        end

        @series begin
            subplot := 1
            xguide := xlabel
            yguide := "p("*xlabel*")"
            seriestype := get(upper, "seriestype", :histogram)
            bins --> get(upper, "nbins", 200)
            normalize --> get(upper, "normalize", true)
            colors --> get(upper, "colors", default_colors)
            intervals --> get(upper, "intervals", default_credibilities)
            legend --> get(upper, "legend", true)

            marg, vsel[1]
        end

        # empty plot (needed since @layout macro not available)
        @series begin
            seriestype := :scatter
            subplot := 2
            grid := false
            xaxis := false
            yaxis := false
            markersize := 0.001
            markerstrokewidth := 0
            markeralpha := 1
            markerstrokealpha := 1
            legend := false
            label := ""
            xguide := ""
            yguide := ""
            [(0,0)]
        end

        @series begin
            subplot := 3
            seriestype := get(diagonal, "seriestype", :histogram)
            xguide --> xlabel
            yguide --> ylabel
            normalize --> get(diagonal, "normalize", true)
            bins --> get(diagonal, "nbins", 200)
            colors --> get(diagonal, "colors", default_colors)
            intervals --> get(diagonal, "intervals", default_credibilities)
            legend --> get(diagonal, "legend", false)

            marg, (vsel[1], vsel[2])
        end

        @series begin
            subplot := 4
            seriestype := get(right, "seriestype", :histogram)
            orientation := :horizontal
            xguide := ylabel
            yguide := "p("*ylabel*")"
            normalize --> get(right, "normalize", true)
            bins --> get(right, "nbins", 200)
            colors --> get(right, "colors", default_colors)
            intervals --> get(right, "intervals", default_credibilities)
            legend --> get(right, "legend", true)

            marg, vsel[2]
        end

    else
        error("seriestype $seriestype not supported")
    end

end

# rectangle bounds
@recipe function f(bounds::HyperRectBounds, vsel::NTuple{2,Integer})
    pi_x, pi_y = vsel

    vol = spatialvolume(bounds)
    vhi = vol.hi[[pi_x, pi_y]]; vlo = vol.lo[[pi_x, pi_y]]
    rect_xy = rectangle_path(vlo, vhi)
    bext = 0.1 * (vhi - vlo)
    xlims = (vlo[1] - bext[1], vhi[1] + bext[1])
    ylims = (vlo[2] - bext[2], vhi[2] + bext[2])

    @series begin
        seriestype := :path
        label --> "bounds"
        seriescolor --> :darkred
        seriesalpha --> 0.75
        linewidth --> 2
        xlims --> xlims
        ylims --> ylims
        (rect_xy[:,1], rect_xy[:,2])
    end

    nothing
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4109a39c-00b7-4954-bb73-35407f82a6f7
using Pkg; Pkg.activate("..")

# ╔═╡ 24b76a7c-63dd-11eb-1b78-d5a20557e5cd
# edit the code below to set your name and UGent username

student = (name = "Jim Janssen", email = "Jim.Janssen@UGent.be");

# press the ▶ button in the bottom right of this cell to run your edits
# or use Shift+Enter

# you might need to wait until all other cells in this notebook have completed running.
# scroll down the page to see what's up

# ╔═╡ 446c3e4e-63dd-11eb-3300-fb006008ff1f
begin
	using DSJulia;
	using PlutoUI;
	tracker = ProgressTracker(student.name, student.email);
	md"""

	Submission by: **_$(student.name)_**
	"""
end

# ╔═╡ 5368424c-63dd-11eb-3a1b-05a61e266d2b
md"""
# Additional exercises for basics

First of all, welcome to these additional exercises, **$(student[:name])**!

This exercise is based on [the advent of code](https://adventofcode.com/), a yearly programming challenge during the advent period.
"""

# ╔═╡ 6c0f2d10-63dd-11eb-2c17-fddf9cd51bfe
md"""## 1. Trajecting trajectories

Suppose you find yourself on top of a snowy mountain, sparsely covered with trees. You want to reach the bottom of the mountain using your sled, but the situation you find yourself in is certainly not safe: there's very minimal steering and the area is covered in trees. You'll need to see which angles will take you near the fewest trees.

Due to the local geology, trees in this area only grow on exact integer coordinates in a grid. You make a map (this will serve as the input later) of the open squares (.) and trees (#) you can see. For example:

```
..##.......
#...#...#..
.#....#..#.
..#.#...#.#
.#...##..#.
..#.##.....
.#.#.#....#
.#........#
#.##...#...
#...##....#
.#..#...#.#
```

These aren't the only trees, though; due to something you read about once involving arboreal genetics and biome stability, the same pattern repeats to the right many times 🤯:

```
..##.........##.........##.........##.........##.........##.......  --->
#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..
.#....#..#..#....#..#..#....#..#..#....#..#..#....#..#..#....#..#.
..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#
.#...##..#..#...##..#..#...##..#..#...##..#..#...##..#..#...##..#.
..#.##.......#.##.......#.##.......#.##.......#.##.......#.##.....  --->
.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#
.#........#.#........#.#........#.#........#.#........#.#........#
#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...#.##...#...
#...##....##...##....##...##....##...##....##...##....##...##....#
.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#.#..#...#.#  --->
```

You start on the open square `(.)` in the top-left corner and need to reach the bottom (below the bottom-most row on your map).

Your sled can only follow a few specific slopes (you opted for a cheaper model that prefers rational numbers); start by counting all the trees you would encounter for the slope right 3, down 1:

From your starting position at the top-left, check the position that is right 3 and down 1. Then, check the position that is right 3 and down 1 from there, and so on until you go past the bottom of the map.

The locations you would check in the above example are marked here with `O` where there was an open square and `X` where there was a tree:

```
..##.........##.........##.........##.........##.........##.......  --->
#..O#...#..#...#...#..#...#...#..#...#...#..#...#...#..#...#...#..
.#....X..#..#....#..#..#....#..#..#....#..#..#....#..#..#....#..#.
..#.#...#O#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#..#.#...#.#
.#...##..#..X...##..#..#...##..#..#...##..#..#...##..#..#...##..#.
..#.##.......#.X#.......#.##.......#.##.......#.##.......#.##.....  --->
.#.#.#....#.#.#.#.O..#.#.#.#....#.#.#.#....#.#.#.#....#.#.#.#....#
.#........#.#........X.#........#.#........#.#........#.#........#
#.##...#...#.##...#...#.X#...#...#.##...#...#.##...#...#.##...#...
#...##....##...##....##...#X....##...##....##...##....##...##....#
.#..#...#.#.#..#...#.#.#..#...X.#.#..#...#.#.#..#...#.#.#..#...#.#  --->
```

In this example, traversing the map using this slope would cause you to encounter 7 trees.

"""

# ╔═╡ e1259736-63df-11eb-0010-3f2fc5719596
md"""
Let us first give you a hand with the data input:
"""

# ╔═╡ 36e51b02-63df-11eb-34a5-39e6e0650541
begin
	download("https://gist.githubusercontent.com/DaanVanHauwermeiren/2b02f275eef08643acfbf8fa9e9e4097/raw/e4e68d69d9c54a9b33a319d7b56dc8daf83fea5b/data.txt", "map.txt")

	function parse_input(input)
		lines = split(input, "\n")
		m = length(first(lines))
		#return [l[j]=='#' for l in lines, j in 1:m]
		return [l[j] for l in lines, j in 1:m]
	end

	map = read("map.txt", String) |> rstrip |> parse_input;
end

# ╔═╡ e677ca7e-63df-11eb-186b-8d20826ec916
map

# ╔═╡ c34af248-63e3-11eb-2bfd-898ef2cfbea5
md"## Questions"

# ╔═╡ 4b4bbeb8-63e2-11eb-0ea2-6de494f61dc1
begin
	# this is where you start
	x, y = 1, 1
	# how big is the mountain?
	n, m = size(map)
	# you start your tree hitting counter at zero of course!
  	n_trees = 0

	# go over the mountain and count the trees!

	n_trees = missing
end

# ╔═╡ dcac803a-648d-11eb-21f2-d72e56528e92

# ╔═╡ fb39928c-63e2-11eb-2b85-15044c04fbe6
product_all_trees = missing

# ╔═╡ ede6342c-648d-11eb-0712-836637408e53

# ╔═╡ b9e97832-63e3-11eb-0e39-b5c9129eabbe
md"## Overly engineered solution to these questions showcasing the 💪 of Julia"

# ╔═╡ 9c015578-648e-11eb-0768-d1e5edb4f65f
md"This is already some advanced use of Julia, so do not worry if this is not yet fully clear. You will learn about this later in this course."

# ╔═╡ f0f0da0a-63e3-11eb-2065-838c77fc57a6
begin
	qalt = Question(;
			description=md"""
		The alternative solutions makes use of the fact that the boundaries are circular. As an example, this would be valid:
		```julia
		A = ['a', 'b', 'c']
		A[1] == A[4] == A[7]
		```
		So we need some datastructure that behaves like a normal array, but with wonky out-of-bounds behaviour!

		In julia, a custom datatype can be defined as:
		```
		struct CircularArray
			data
		end
		```
		Now we have our data type `CircularArray`, which has a field `data`.
		Using our knowledge from the simple example above, let us define a function which takes an index and the length of the array and returns the correct index using circular bounds. We will also create a function that takes an array of indices and returns an array of correct indices. Take a moment to think about how you would do this, and afterwards check a possible implementation in hint 1.

		So now we know how the compute the **normal** indices from the circular ones.
		We now need to define some things for our new datatype `CircularArray`.
		1. We need the `size` of `CircularArray` and
		2. we need to define how `A[1]` works on our datatype (this is the `getindex` function in julia).
		For the first problem, note that the size of a `CircularArray` is actually the size of the underlying array! So how do you think the function `Base.size(A::CircularArray) = ...` looks like?
		For the second problem we need a definition of `Base.getindex(A::CircularArray, i::Int)`. How can we define this using the `circindex` function?
		If you have given this some thought, check out hint no 2 for the answer!.
				""")

   QuestionBlock(;
	title=md"**Alternative answer: defining a new data type**",
	questions = [qalt],
	hints=[hint(md"
```julia
circindex(i::Int, N::Int) = 1 + mod(i-1,N)
circindex(I, N::Int) = [circindex(i,N)::Int for i in I]
```"),
			hint(md"
Note that I use `length`. This is because I already know that `data` will be a `string`, as you will see below.
```julia
Base.size(A::CircularArray) = length(A.data)

Base.getindex(A::CircularArray, i::Int) = getindex(A.data, circindex(i, size(A)))
Base.getindex(A::CircularArray, I) = getindex(A.data, circindex(I, size(A)))
```")
			],
	)
end

# ╔═╡ ec0aa2c6-63ea-11eb-062a-cdd3e81a743a
md"""
Let us defined what we talked about before:
"""

# ╔═╡ de9cf7b0-63ea-11eb-12fb-773d79c2463b
begin
	struct CircularArray
		data
	end

	circindex(i::Int,N::Int) = 1 + mod(i-1,N)
	circindex(I,N::Int) = [circindex(i,N)::Int for i in I]

	Base.size(A::CircularArray) = length(A.data)

	Base.getindex(A::CircularArray, i::Int) = getindex(A.data, circindex(i, size(A)))
	Base.getindex(A::CircularArray, I) = getindex(A.data, circindex(I, size(A)))
end

# ╔═╡ 424c1470-63e0-11eb-323f-454ac8ad273e
begin
	q1 = Question(validators = @safe[n_trees == Solutions.n_trees],
			description=md"""
			Starting at the top-left corner of your map and following a slope of right 3 and down 1, how many trees would you encounter?
			""")

   qb1 = QuestionBlock(;
	title=md"**Question 1: counting trees**",
	questions = [q1],
	hints=[hint(md"This question can be answered without defining your own datatypes and using all that is available from the `Julia.Base` module. This can also be completed without knowing a lot about collections. Since this is only explained in the next part of the course. "),
			hint(md"An alterative solution showcasing the power of Julia will be shown below, it includes some teasers of what is yet to come in this course")
			],
	)
	validate(qb1, tracker)
end

# ╔═╡ c1206434-63e2-11eb-16c1-d7fafb48bbbe
begin
	q2 = Question(validators = @safe[product_all_trees == Solutions.product_all_trees],
			description=md"""
Time to check the rest of the slopes - you need to minimize the probability of a sudden arboreal stop, after all.

Determine the number of trees you would encounter if, for each of the following slopes, you start at the top-left corner and traverse the map all the way to the bottom:

    Right 1, down 1.
    Right 3, down 1. (This is the slope you already checked.)
    Right 5, down 1.
    Right 7, down 1.
    Right 1, down 2.

In the above example, these slopes would find 2, 7, 3, 4, and 2 tree(s) respectively; multiplied together, these produce the answer 336.

What do you get if you multiply together the number of trees encountered on each of the listed slopes?""")

   qb2 = QuestionBlock(;
	title=md"**Question 2: counting trees**",
	questions = [q2],
	)
	validate(qb2, tracker)
end

# ╔═╡ f752b934-63ea-11eb-09cc-99172ed4870d
md"""
Now, let us read in the data, but this time each horizontal line is a `CircularArray`! (Try to understand what each line does in the code below)
"""

# ╔═╡ 5c2438ce-63eb-11eb-0892-272fb362a2e8
begin
	# init 1D array of data type CircularArrays with len=number of lines
	# in the file and leave it undefined for now
	# Note that the number of lines is hardcoded (323).
	# This could be done more elegantly in an other way, but we will leave that
	# exercise for the user.
	slope = Array{CircularArray, 1}(undef, 323)

	import Base.Iterators.enumerate

	# open the file
	open("map.txt") do file
		# if you have some programming experience, you can probably guess
		# what enumerate does.
		# the eachline function speaks for itself: iterate over each line in the file!
		for (idx, ln) in enumerate(eachline(file))
			slope[idx] = CircularArray(ln)
		end
	end
end

# ╔═╡ 9cafdef2-63eb-11eb-326a-5d51766d17d5
md"This is how the data now looks like: (Click to expand and have a nicer view of the circular mountain)"

# ╔═╡ 9371b694-63eb-11eb-30a5-d3542ca8f4e5
slope

# ╔═╡ 18054676-63f0-11eb-1a53-5b65f8477882
slope[1]

# ╔═╡ 73f03a14-63f0-11eb-1a14-17e4e69745af
md"Let us look at one horizontal line of the mountain. Check the output below of data at index 1, 32 and 63. Does this match your intuition?"

# ╔═╡ 2164576e-63f0-11eb-3e57-4fad123c494b
firstline = slope[1]

# ╔═╡ 8eda596a-63f0-11eb-03ba-abd147192428
size(firstline)

# ╔═╡ 62dfc3ae-63f0-11eb-350c-b3040f8a6d9b
firstline[1], firstline[32], firstline[63]

# ╔═╡ 224c26ce-63ec-11eb-0744-e97f5196972a
md"""
If we want to go down the slope using our sled in the configuration of 1 down + 3 right, we can write down the *mountain coordinates* like this:

```julia
(1, 1)
(2, 4)
(3, 7)
(4, 10)
(5, 13)
...
```
We have not talked about collections and ranges, yet so this will be a sneak peak. You can define a range in julia like this: `1:5`. This will give you the values `[1,2,3,4,5]`. On the other hand, `1:2:5` will give you `[1,3,5]`. Do you see the difference?
change some values in the two cells below to check if you understand the basics on this concept. Then, think about how you would write up the list of x-coordinates and y-coordinates if you go down the slope 1 down, 3 right.
"""

# ╔═╡ 9ec3845e-63ec-11eb-1a81-c141b4ff37d8
collect(1:5)

# ╔═╡ b590457a-63ec-11eb-3b3b-5d5b0de1f98c
collect(1:2:5)

# ╔═╡ a0f50fee-63ed-11eb-0208-0b5f7036a7b1
md"Moving on, let us define our function to go down the slope. Does this match your intuition of before?"

# ╔═╡ ae58f724-63eb-11eb-0541-3584d99d9fc5
function goingdowntheslope(right, down, slope)
    xvals = 1:right:right*323
    yvals = 1:down:323
    thingsyoufind = [slope[yval][xval] for (xval, yval) in zip(xvals, yvals)]
    return sum(thingsyoufind .== '#')
end

# ╔═╡ 44885292-63ee-11eb-048b-2b7481a83ae8
md"""Now that you possess all this new knowledge, can find the answers to the two slope-questions in this notebook?"""

# ╔═╡ 10e183ca-63ec-11eb-3d4b-57ed0543cda1
goingdowntheslope(3, 1, slope)

# ╔═╡ dd716abc-63ed-11eb-36e9-fdbd1187146e
slopestyles = (
    (1, 1),
    (3, 1),
    (5, 1),
    (7, 1),
    (1, 2),
)

# ╔═╡ e0569568-63ed-11eb-33b9-51a8b179292c
prod([goingdowntheslope(direction..., slope) for direction in slopestyles])

# ╔═╡ Cell order:
# ╠═24b76a7c-63dd-11eb-1b78-d5a20557e5cd
# ╟─4109a39c-00b7-4954-bb73-35407f82a6f7
# ╠═446c3e4e-63dd-11eb-3300-fb006008ff1f
# ╟─5368424c-63dd-11eb-3a1b-05a61e266d2b
# ╟─6c0f2d10-63dd-11eb-2c17-fddf9cd51bfe
# ╟─e1259736-63df-11eb-0010-3f2fc5719596
# ╠═36e51b02-63df-11eb-34a5-39e6e0650541
# ╠═e677ca7e-63df-11eb-186b-8d20826ec916
# ╟─c34af248-63e3-11eb-2bfd-898ef2cfbea5
# ╟─424c1470-63e0-11eb-323f-454ac8ad273e
# ╠═4b4bbeb8-63e2-11eb-0ea2-6de494f61dc1
# ╟─dcac803a-648d-11eb-21f2-d72e56528e92
# ╟─c1206434-63e2-11eb-16c1-d7fafb48bbbe
# ╠═fb39928c-63e2-11eb-2b85-15044c04fbe6
# ╟─ede6342c-648d-11eb-0712-836637408e53
# ╟─b9e97832-63e3-11eb-0e39-b5c9129eabbe
# ╟─9c015578-648e-11eb-0768-d1e5edb4f65f
# ╟─f0f0da0a-63e3-11eb-2065-838c77fc57a6
# ╟─ec0aa2c6-63ea-11eb-062a-cdd3e81a743a
# ╠═de9cf7b0-63ea-11eb-12fb-773d79c2463b
# ╟─f752b934-63ea-11eb-09cc-99172ed4870d
# ╠═5c2438ce-63eb-11eb-0892-272fb362a2e8
# ╟─9cafdef2-63eb-11eb-326a-5d51766d17d5
# ╠═9371b694-63eb-11eb-30a5-d3542ca8f4e5
# ╠═18054676-63f0-11eb-1a53-5b65f8477882
# ╟─73f03a14-63f0-11eb-1a14-17e4e69745af
# ╠═2164576e-63f0-11eb-3e57-4fad123c494b
# ╠═8eda596a-63f0-11eb-03ba-abd147192428
# ╠═62dfc3ae-63f0-11eb-350c-b3040f8a6d9b
# ╟─224c26ce-63ec-11eb-0744-e97f5196972a
# ╠═9ec3845e-63ec-11eb-1a81-c141b4ff37d8
# ╠═b590457a-63ec-11eb-3b3b-5d5b0de1f98c
# ╟─a0f50fee-63ed-11eb-0208-0b5f7036a7b1
# ╠═ae58f724-63eb-11eb-0541-3584d99d9fc5
# ╟─44885292-63ee-11eb-048b-2b7481a83ae8
# ╠═10e183ca-63ec-11eb-3d4b-57ed0543cda1
# ╠═dd716abc-63ed-11eb-36e9-fdbd1187146e
# ╠═e0569568-63ed-11eb-33b9-51a8b179292c


text:
using Documenter, MetaGraphs

# index is equal to the README for the time being
cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/index.md"); force=true)

# same for license
cp(normpath(@__FILE__, "../../LICENSE.md"), normpath(@__FILE__, "../src/license.md"); force=true)

makedocs(
    modules = [MetaGraphs],
    sitename = "MetaGraphs",
    format = Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    pages = Any[
        "Overview"             => "index.md",
        "MetaGraphs Functions" => "metagraphs.md",
        "License Information"  => "license.md",
    ],
)

deploydocs(repo="github.com/JuliaGraphs/MetaGraphs.jl.git", push_preview=true)

rm(normpath(@__FILE__, "../src/index.md"))
rm(normpath(@__FILE__, "../src/license.md"))


text:
using PyCall
using DataFrames
using CSV
using BenchmarkTools
using BenchmarkTools: Trial, TrialEstimate, median, mean

py"""
import torch
import torch_scatter as sc
torch.set_num_threads(12)
cuda = torch.device("cuda:0")
d = 50
nbins = 20
"""

d = 50
nbins = 20
getinfo(te::TrialEstimate) = te.time, te.gctime, te.memory
getstats(t::Trial) = [getinfo(minimum(t)), getinfo(median(t)), getinfo(mean(t)),
                      getinfo(maximum(t))]

metadata = DataFrame(device=String[], dim=Int[], sample=Int[], bins=Int[])
mintime = DataFrame(min_time=Float64[], min_gc=Float64[], min_mem=Int[])
medtime = DataFrame(med_time=Float64[], med_gc=Float64[], med_mem=Int[])
meantime = DataFrame(mean_time=Float64[], mean_gc=Float64[], mean_mem=Int[])
maxtime = DataFrame(max_time=Float64[], max_gc=Float64[], max_mem=Int[])

for l = [2^5, 2^10, 2^15, 2^20]
    py"""
    hist = torch.zeros([d, nbins], dtype=torch.float32)
    delta = torch.rand([d, $(l)], dtype=torch.float32)
    idx = torch.randint(0, nbins, size=($(l),))

    hist_gpu = torch.zeros([d, nbins], dtype=torch.float32, device=cuda)
    delta_gpu = torch.rand([d, $(l)], dtype=torch.float32, device=cuda)
    idx_gpu = torch.randint(0, nbins, size=($(l),), device=cuda)

    sc.scatter_add(delta, idx, out=hist)
    sc.scatter_add(delta_gpu, idx_gpu, out=hist_gpu)
    """

    b_cpu = @benchmark py"sc.scatter_add(delta, idx, out=hist)";
    b_gpu = @benchmark py"sc.scatter_add(delta_gpu, idx_gpu, out=hist_gpu)";
    s_cpu = getstats(b_cpu)
    s_gpu = getstats(b_gpu)

    push!(metadata, ("cpu", d, l, nbins))
    push!(mintime, s_cpu[1])
    push!(medtime, s_cpu[2])
    push!(meantime, s_cpu[3])
    push!(maxtime, s_cpu[4])

    push!(metadata, ("gpu", d, l, nbins))
    push!(mintime, s_gpu[1])
    push!(medtime, s_gpu[2])
    push!(meantime, s_gpu[3])
    push!(maxtime, s_gpu[4])
end

data = hcat(metadata, mintime, medtime, meantime, maxtime)
CSV.write("benchmark_scatter_pytorch.tsv", data; delim="\t")


text:
@testset "Downloader selection" begin
    @test LightOSM.osm_network_downloader(:place_name) == LightOSM.osm_network_from_place_name
    @test LightOSM.osm_network_downloader(:bbox) == LightOSM.osm_network_from_bbox
    @test LightOSM.osm_network_downloader(:point) == LightOSM.osm_network_from_point
    @test LightOSM.osm_network_downloader(:polygon) == LightOSM.osm_network_from_polygon
    @test_throws ArgumentError LightOSM.osm_network_downloader(:doesnt_exist)
end

function wait_for_overpass()
    count = 0
    while !LightOSM.is_overpass_server_availabile()
        count == 7 && break # Can't wait indefinitely, and the failure will be caught in the tests
        count += 1
        @info "Waiting for overpass server..."
        sleep(5 * count)
    end
end

@testset "Downloads" begin
    filenames = ["map.osm", "map.json"]
    formats = [:osm, :json]
    for (filename, format) in zip(filenames, formats)
        try
            wait_for_overpass()
            data = download_osm_network(:point,
                                        radius=0.5,
                                        point=GeoLocation(-37.8136, 144.9631),
                                        network_type=:drive,
                                        download_format=format,
                                        save_to_file_location=filename);
            @test isfile(filename)
            g = graph_from_file(filename) # Check it doesn't error
        catch err
            # Sometimes gets HTTP.ExceptionRequest.StatusError in tests due to connection to overpass
            !isa(err, HTTP.ExceptionRequest.StatusError) && rethrow()
            @error "Test failed due to connection issue" exception=(err, catch_backtrace())
        end

        try
            rm(filename)
        catch
        end
    end
end