text:
#///////////////////////////////////////
#// File Name: cross_entropy_bilevel_optimization.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/10/28
#// Description: Cross Entropy Bilevel Optimization (RATiLQR) algorithm
#///////////////////////////////////////

using Random
using Distributions
using LinearAlgebra
using Distributed

"""
    CrossEntropyBilevelOptimizationSolver(kwargs...)

RAT iLQR (i.e. Cross Entropy Method + iLEQG) Solver.

# Optional Keyword Arguments
## iLEQG Solver Parameters
- `μ_min_ileqg::Float64` -- minimum value for Hessian regularization parameter `μ` (> 0).
  Default: `1e-6`.
- `Δ_0_ileqg::Float64` -- minimum multiplicative modification factor (> 0) for `μ`.
  Default: `2.0`.
- `λ_ileqg::Float64` -- multiplicative modification factor in (0, 1) for line search
  step size `ϵ`. Default: `0.5`.
- `d_ileqg::Float64` -- convergence error norm threshold (> 0). If the maximum l2
  norm of the change in nominal control over the horizon is less than `d`, the
  solver is considered to be converged. Default: `1e-2`.
- `iter_max_ileqg::Int64` -- maximum iteration number. Default: 100
- `ϵ_init_ileqg::Float64` -- initial step size in (`ϵ_min`, 1] to start the
  backtracking line search with. If `adaptive_ϵ_init` is `true`, then this
  value is overridden by the solver's adaptive initialization functionality
  after the first iLEQG iteration. If `adaptive_ϵ_init` is `false`, the
  specified value of `ϵ_init` is used across all the iterations as the initial
  step size. Default:`1.0`.
- `adaptive_ϵ_init_ileqg::Bool` -- if `true`, `ϵ_init` is adaptively changed based on
  the last step size `ϵ` of the previous iLEQG iteration. Default: `false`.
   - If the first line search iterate `ϵ_init_prev` in the previous iLEQG
     iteration is successful, then `ϵ_init` for the next iLEQG iteration is set
     to `ϵ_init = ϵ_init_prev / λ` so that the initial line search step increases.
   - Otherwise `ϵ_init = ϵ_last` where `ϵ_last` is the line search step accepted
     in the previous iLEQG iteration.
- `ϵ_min_ileqg::Float64` -- minimum value of step size `ϵ` to terminate the line
  search. When `ϵ_min` is reached, the last candidate nominal trajectory is accepted
  regardless of the Armijo condition and the current iLEQG iteration is
  finished. Default: `1e-6`.
- `f_returns_jacobian::Bool` -- if `true`, Jacobian matrices of the dynamics function
  are user-provided. This can reduce computation time since automatic
  differentiation is not used. Default: `false`.

## Cross Entropy Solver Parameters
- `μ_init::Float64` -- initial value of the mean parameter `μ` used in the first
  Cross Entropy iteration. Default: `1.0`.
- `σ_init::Float64` -- initial value of the standard deviation parameter `σ`
  used in the first Cross Entropy iteration. Default: `2.0`.
- `num_samples::Int64` -- number of Monte Carlo samples for the risk-sensitivity
  parameter `θ`. Default: `10`.
- `num_elite::Int64` -- number of elite samples. Default: `3`.
- `iter_max::Int64` -- maximum iteration number. Default: `5`.
- `λ::Float64` -- multiplicative modification factor in (0, 1) for ``μ_init` and
  `σ_init`. Default: `0.5`.
- `use_θ_max::Bool` -- if `true`, the maximum feasible `θ` found is used to
  perform the final iLEQG optimization instead of the optimal one. Default: `false`.

# Notes
- The values of `μ_init` and `σ_init`, which may be modified during optimization,
  are stored internally in the solver and　carried over to the next call to `solve!`.
"""
mutable struct CrossEntropyBilevelOptimizationSolver
    # ileqg solver parameters
    μ_min_ileqg::Float64           # Minimum damping parameter for regularization
    Δ_0_ileqg::Float64             # Minimum Modification factor for μ
    λ_ileqg::Float64               # Multiplicative factor for line search step parameter in (0, 1)
    d_ileqg::Float64               # Convergence error norm thresholds
    iter_max_ileqg::Int64          # Maximum iteration
    # β_ileqg::Float64               # Armijo search condition
    ϵ_init_auto_ileqg::Bool        # Automatic initialization of ϵ_init from the previous iLEQG iteration.
    ϵ_init_ileqg::Float64          # Initial step size for backtracking line search
    ϵ_min_ileqg::Float64           # Minimum step size for backtracking line search
    f_returns_jacobian::Bool       # Whether the dynamics funcion f also returns jacobians or not

    # CE solver parameters
    num_samples::Int64             # Number of CE samples
    num_elite::Int64               # Number of elite samples
    iter_max::Int64                # Maximum iteration
    λ::Float64                     # Multiplicative factor for μ and σ in (0, 1)
    use_θ_max::Bool                # Use maximum feasible θ instead of optimal one

    # CE solver mutable parameters
    μ_init::Float64                # Initial mean parameter
    σ_init::Float64                # Initial standard deviation parameter
    μ::Float64                     # Current mean parameter
    σ::Float64                     # Current standard deviation parameter
    θ_max::Float64                 # Maximum valid θ encountered so far
    θ_min::Float64                 # Minimum valid θ encountered so far
    iter_current::Int64            # Current CE iteration
end

function CrossEntropyBilevelOptimizationSolver(;μ_min_ileqg=1e-6,
                                               Δ_0_ileqg=2.0,
                                               λ_ileqg=0.5,
                                               d_ileqg=1e-2,
                                               iter_max_ileqg=100,
                                               # β_ileqg=1e-4,
                                               adaptive_ϵ_init_ileqg=false,
                                               ϵ_init_ileqg=1.0,
                                               ϵ_min_ileqg=1e-6,
                                               μ_init=1.0,
                                               σ_init=2.0,
                                               num_samples=10,
                                               num_elite=3,
                                               iter_max=5,
                                               λ=0.5,
                                               f_returns_jacobian=false,
                                               use_θ_max=false)

    μ, σ = μ_init, σ_init
    θ_max, θ_min = 0.0, Inf
    iter_current = 0

    return CrossEntropyBilevelOptimizationSolver(μ_min_ileqg, Δ_0_ileqg, λ_ileqg, d_ileqg,
                                                 iter_max_ileqg, adaptive_ϵ_init_ileqg,
                                                 ϵ_init_ileqg, ϵ_min_ileqg, f_returns_jacobian,
                                                 num_samples, num_elite, iter_max, λ, use_θ_max,
                                                 μ_init, σ_init, μ, σ, θ_max, θ_min, iter_current)
end

"""
Initialize RATiLQR Solver
"""
function initialize!(ce_solver::CrossEntropyBilevelOptimizationSolver)
    ce_solver.iter_current = 0;
    ce_solver.μ, ce_solver.σ = ce_solver.μ_init, ce_solver.σ_init
    ce_solver.θ_max = 0.0;
    ce_solver.θ_min = Inf;
end

"""
Compute iLEQG value on a worker process
"""
function compute_value_worker(ce_solver::CrossEntropyBilevelOptimizationSolver,
                              problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                              x::Vector{Float64}, u_array::Vector{Vector{Float64}},
                              θ::Float64)
    ileqg = ILEQGSolver(problem,
                        μ_min=ce_solver.μ_min_ileqg,
                        Δ_0=ce_solver.Δ_0_ileqg,
                        λ=ce_solver.λ_ileqg,
                        d=ce_solver.d_ileqg,
                        iter_max=ce_solver.iter_max_ileqg,
                        # β=ce_solver.β_ileqg,
                        adaptive_ϵ_init=ce_solver.ϵ_init_auto_ileqg,
                        ϵ_init=ce_solver.ϵ_init_ileqg,
                        ϵ_min=ce_solver.ϵ_min_ileqg,
                        f_returns_jacobian=ce_solver.f_returns_jacobian)
    #initialize!(ileqg, problem, x, u_array, θ)
    value = 0.0
    try
        value = solve!(ileqg, problem, x, u_array, θ=θ, verbose=false)[4];
    catch
        value = Inf
    end
    return value
end

"""
Compute iLEQG values in parallel on multiple worker processes
"""
function compute_cost(ce_solver::CrossEntropyBilevelOptimizationSolver,
                      problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                      x::Vector{Float64},
                      u_array::Vector{Vector{Float64}},
                      θ_array::Vector{Float64},
                      kl_bound::Float64)
    num_samples = length(θ_array)
    if nprocs() > 1
        proc_id_array = 2 .+ [mod(ii, nprocs() - 1) for ii = 0 : num_samples - 1]
    else
        proc_id_array = [1 for ii = 0 : num_samples - 1]
    end
    value_array = Vector{Float64}(undef, num_samples)
    @sync begin
        for ii = 1 : num_samples
            @inbounds @async value_array[ii] =
                remotecall_fetch(compute_value_worker, proc_id_array[ii],
                                ce_solver, problem, x, u_array, θ_array[ii])
        end
    end
    cost_array = value_array .+ kl_bound./θ_array
    return cost_array
end

# For debugging only. It should work in the same way as compute_cost
function compute_cost_serial(ce_solver::CrossEntropyBilevelOptimizationSolver,
                             problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                             x::Vector{Float64},
                             u_array::Vector{Vector{Float64}},
                             θ_array::Vector{Float64},
                             kl_bound::Float64)
    @assert length(θ_array) == ce_solver.num_samples
    cost_array = Vector{Float64}(undef, ce_solver.num_samples)
    for ii = 1 : ce_solver.num_samples
        ileqg_solver = ILEQGSolver(problem,
                                   μ_min=ce_solver.μ_min_ileqg,
                                   Δ_0=ce_solver.Δ_0_ileqg,
                                   λ=ce_solver.λ_ileqg,
                                   d=ce_solver.d_ileqg,
                                   iter_max=ce_solver.iter_max_ileqg,
                                   # β=ce_solver.β_ileqg,
                                   adaptive_ϵ_init=ce_solver.ϵ_init_auto_ileqg,
                                   ϵ_init=ce_solver.ϵ_init_ileqg,
                                   ϵ_min=ce_solver.ϵ_min_ileqg,
                                   f_returns_jacobian=ce_solver.f_returns_jacobian)
        #initialize!(ileqg_solver, problem, x, u_array, θ_array[ii])
        try
            cost_array[ii] = solve!(ileqg_solver, problem, x, u_array, θ=θ_array[ii], verbose=false)[4] +
                             kl_bound/θ_array[ii]
        catch
            cost_array[ii] = Inf
        end
    end
    return cost_array
end

"""
Get positive θ samples
"""
function get_positive_samples(μ::Float64, σ::Float64, num_samples::Int64, rng::AbstractRNG)
    θ_array = Float64[];
    d = Normal(μ, σ)
    while true
        θ_sampled = rand(rng, d);
        if θ_sampled > 0.0
            push!(θ_array, θ_sampled)
        end
        if length(θ_array) >= num_samples
            break
        end
    end
    return θ_array
end

"""
Single iteration of RATiLQR
"""
function step!(ce_solver::CrossEntropyBilevelOptimizationSolver,
               problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
               x::Vector{Float64},
               u_array::Vector{Vector{Float64}},
               kl_bound::Float64,
               rng::AbstractRNG,
               verbose=true, serial=false)
    ce_solver.iter_current += 1;
    if verbose
        println("**CE iteration $(ce_solver.iter_current)")
    end
    θ_array = Float64[];
    costs_array = Float64[];
    while true
        if ce_solver.iter_current == 1
            # draw from N(μ_init, σ_init)
            # if too few valid samples, then adjust μ_init, σ_init and redraw.
            # if all samples are valid, then increase μ_init, σ_init for the next iteration.
            if verbose
                println("****Drawing $(ce_solver.num_samples) positive samples of θ ~ N($(round(ce_solver.μ_init,digits=4)), $(round(ce_solver.σ_init,digits=4)))");
            end
            θ_array = get_positive_samples(ce_solver.μ_init, ce_solver.σ_init, ce_solver.num_samples, rng);
        else
            if verbose
                println("****Drawing $(ce_solver.num_samples) positive samples of θ ~ N($(round(ce_solver.μ,digits=4)), $(round(ce_solver.σ,digits=4)))");
            end
            θ_array = get_positive_samples(ce_solver.μ, ce_solver.σ, ce_solver.num_samples, rng);
        end
        if verbose
            println("****Evaluating costs of sampled points")
        end
        if !serial
            costs_array = compute_cost(ce_solver, problem, x, u_array, θ_array, kl_bound)
        else
            costs_array = compute_cost_serial(ce_solver, problem, x, u_array, θ_array, kl_bound)
        end
        if verbose
            println(costs_array)
        end
        num_inf = sum(isinf.(costs_array));
        num_valid = ce_solver.num_samples - num_inf;
        if ce_solver.iter_current == 1 && num_valid < max(ce_solver.num_elite, ce_solver.num_samples*ce_solver.λ)
            if verbose
                println("******$(num_inf)/$(ce_solver.num_samples) samples are Inf. Redrawing samples")
            end
            ce_solver.μ_init *= ce_solver.λ
            ce_solver.σ_init *= ce_solver.λ
        elseif ce_solver.iter_current == 1 && num_valid == ce_solver.num_samples
            ce_solver.μ_init /= ce_solver.λ
            ce_solver.σ_init /= ce_solver.λ
            if verbose
                println("******Increasing μ_init to $(ce_solver.μ_init) and σ_init to $(ce_solver.σ_init)")
            end
            break
        elseif num_valid >= max(ce_solver.num_elite, ce_solver.num_samples*ce_solver.λ)
            if verbose
                println("******$(num_valid)/$(ce_solver.num_samples) samples are valid")
            end
            break
        end
    end

    for ii = 1 : length(θ_array)
        if isinf(costs_array[ii])
            continue
        else
            if θ_array[ii] < ce_solver.θ_min
                ce_solver.θ_min = θ_array[ii]
            elseif θ_array[ii] > ce_solver.θ_max
                ce_solver.θ_max = θ_array[ii]
            end
        end
    end

    θ_cost_pair_array = collect(zip(θ_array, costs_array));
    θ_cost_pair_sorted_array = sort(θ_cost_pair_array, by=x->x[2]);
    θ_elite_array = [pair[1] for pair in θ_cost_pair_sorted_array[1 : ce_solver.num_elite]];
    μ_new = sum(θ_elite_array)/ce_solver.num_elite;
    σ_new = sqrt(sum((θ_elite_array .- μ_new).^2)/ce_solver.num_elite)
    if verbose
        println("****Updated with μ_new: $(round(μ_new,digits=4)) and σ_new: $(round(σ_new,digits=4))")
    end
    ce_solver.μ, ce_solver.σ = μ_new, σ_new;
end

"""
    solve!(ce_solver::CrossEntropyBilevelOptimizationSolver,
    problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
    x_0::Vector{Float64}, u_array::Vector{Vector{Float64}}, rng::AbstractRNG;
    kl_bound::Float64, verbose=true, serial=false)

Given `problem` and `ce_solver` (i.e. a RAT iLQR Solver), solve distributionally robust
control with current state `x_0` and nominal control schedule `u_array = [u_0, ..., u_{N-1}]`
under the KL divergence bound of `kl_bound` (>= 0).

# Return Values (Ordered)
- `θ_opt::Float64` -- optimal risk-sensitivity parameter.
- `x_array::Vector{Vector{Float64}}` -- nominal state trajectory `[x_0,...,x_N]`.
- `l_array::Vector{Vector{Float64}}` -- nominal control schedule `[l_0,...,l_{N-1}]`.
- `L_array::Vector{Matrix{Float64}}` -- feedback gain schedule `[L_0,...,L_{N-1}]`.
- `value::Float64` -- optimal cost-to-go (i.e. objective value) found by the solver.
- `θ_min::Float64` -- minimum feasible risk-sensitivity parameter found.
- `θ_max::Float64` -- maximum feasible risk-sensitivity parameter found.

# Notes
- Returns a time-varying affine state-feedback policy `π_k` of the form
  `π_k(x) = L_k(x - x_k) + l_k`.
- If `kl_bound` is 0.0, the solver reduces to iLQG.
- If `serial` is `true`, Monte Carlo sampling of the Cross Entropy method is serialized
  on a single process. If `false` it is distributed on all the available worker processes.
"""
function solve!(ce_solver::CrossEntropyBilevelOptimizationSolver,
                problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                x_0::Vector{Float64}, u_array::Vector{Vector{Float64}}, rng::AbstractRNG;
                kl_bound::Float64, verbose=true, serial=false)
    @assert kl_bound >= 0 "KL Divergence Bound must be non-negative"
    initialize!(ce_solver);
    if kl_bound > 0
        while ce_solver.iter_current < ce_solver.iter_max
            step!(ce_solver, problem, x_0, u_array, kl_bound, rng, verbose, serial)
        end
        θ_min, θ_max = ce_solver.θ_min, ce_solver.θ_max;
        if ce_solver.use_θ_max
            θ_opt = θ_max
            if verbose
                println("Using θ_max = $(round(θ_max, digits=4)) in place of μ = $(round(ce_solver.μ, digits=4))")
            end
        else
            θ_opt = ce_solver.μ
        end
        if verbose
            println("Maximum iteration number reached. θ_opt = $(round(θ_opt,digits=4)) from [$(round(θ_min,digits=4)), $(round(θ_max,digits=4))]")
        end
    else
        # kl_bound is 0. In this case the optimal controller is the risk-neutral (iLQG) controller.
        θ_opt = 0.0;
    end
    while true
        try
            ileqg_solver = ILEQGSolver(problem,
                                       μ_min=ce_solver.μ_min_ileqg,
                                       Δ_0=ce_solver.Δ_0_ileqg,
                                       λ=ce_solver.λ_ileqg,
                                       d=ce_solver.d_ileqg,
                                       iter_max=ce_solver.iter_max_ileqg,
                                       # β=ce_solver.β_ileqg,
                                       adaptive_ϵ_init=ce_solver.ϵ_init_auto_ileqg,
                                       ϵ_init=ce_solver.ϵ_init_ileqg,
                                       ϵ_min=ce_solver.ϵ_min_ileqg,
                                       f_returns_jacobian=ce_solver.f_returns_jacobian)
            #initialize!(ileqg_solver, problem, x_0, u_array, θ_opt)
            x_array, l_array, L_array, value, ~ = solve!(ileqg_solver, problem, x_0, u_array, θ=θ_opt, verbose=verbose)
            if kl_bound > 0
                return θ_opt, x_array, l_array, L_array, value + kl_bound/θ_opt, θ_min, θ_max
            else
                return θ_opt, x_array, l_array, L_array, value, 0.0, 0.0
            end
        catch
            @warn "θ_opt == $(θ_opt) resulted in neurotic breakdown. Re-trying with θ_opt == $(max(0.0, θ_opt - ce_solver.σ))"
            θ_opt = max(0.0, θ_opt - ce_solver.σ)
        end
    end
end


text:
using Printf

function looping_examples()
    println("--------------------------------------")
    println("Looping examples")

    println()

    println("We can use while loops to count to 20, displaying even numbers, and aborting when we get past 10")
    i = 1

    while i < 20
        if (i % 2) == 0
            println(i)
            i += 1
            continue # Skip the rest of the code in the loop
        end

        # Make the code use the global i
        i += 1

        if i > 10
            break # Jump out of the loop
        end
    end

    println("We can use for loops to go from one number to another, e.g. 1-5")
    for i = 1:5
        println(i)
    end

    println("We can use for loops to go from one number to another, e.g. 1-100, stepping at 10 using 1:10:100")
    for i = 1:10:100
        println(i)
    end

    println("We can use for loops to 'for-each' value in an array")
    for i in [2,4,6]
        println(i)
    end

    println("We can use for loops with multiple variables (a for-loop within a for-loop)")
    for i = 1:5, j = 1:10
        println((i, j))
    end
end

ext: jl

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
function by_single_type(t::Symbol)
    function single_type(model::ABM)
        ids = collect(allids(model))
        filter!(id -> model[id].type === t, ids)
        return ids
    end
    return single_type
end

function nearby_t(t::Symbol, id, model)
    ids = collect(nearby_ids(model[id].pos, model, model.properties.infection_distance))
    filter!(id -> model[id].type === t, ids)
    return ids
end

function random_without_bacterium(model::ABM{<:DiscreteSpace}, cutoff=0.998)
    function has_no_bacterium(pos, model)
        ids = ids_in_position(pos, model)
        return !any(id -> model[id].type === :bacterium, ids)
    end

    function no_bacterium_positions(model::ABM{<:DiscreteSpace})
        return Iterators.filter(i -> has_no_bacterium(i, model), positions(model))
    end

    if clamp(nagents(model) / prod(size(model.space.s)), 0.0, 1.0) < cutoff
        while true
            pos = random_position(model)
            has_no_bacterium(pos, model) && return pos
        end
    else
        no_bacterium = no_bacterium_positions(model)
        isempty(no_bacterium) && return ()
        return rand(model.rng, collect(no_bacterium))
    end
end

function random_without_bacterium(model::ABM{<:DiscreteSpace}, positions::Vector{Tuple{Int,Int}})
    function has_no_bacterium(pos, model)
        ids = ids_in_position(pos, model)
        return !any(id -> model[id].type === :bacterium, ids)
    end

    function no_bacterium_positions(model::ABM{<:DiscreteSpace}, positions)
        return Iterators.filter(i -> has_no_bacterium(i, model), positions)
    end

    no_bacterium = no_bacterium_positions(model, positions)
    isempty(no_bacterium) && return ()
    return rand(model.rng, collect(no_bacterium))
end

function add_bacterium_single!(
    agent::A,
    model::ABM{<:DiscreteSpace,A};
    cutoff=0.998
) where {A<:AbstractAgent}
    position = random_without_bacterium(model, cutoff)
    isempty(position) && return
    agent.pos = position
    add_agent_pos!(agent, model) >
    return agent
end

function move_bacterium_single!(
    agent::A,
    model::ABM{<:DiscreteSpace,A},
    positions
) where {A<:AbstractAgent}
    position = random_without_bacterium(model, positions)
    isempty(position) && return
    move_agent!(agent, position, model)
    return agent
end


text:
abstract type Scheme end

struct http <: Scheme end
struct https <: Scheme end
# struct ws <: Scheme end
# struct wss <: Scheme end

sockettype(::Type{http}) = TCPSocket
sockettype(::Type{https}) = TLS.SSLContext
schemetype(::Type{TCPSocket}) = http
schemetype(::Type{TLS.SSLContext}) = https

const Headers = Dict{String,String}

const Option{T} = Union{T, Void}
not(::Void) = true
not(x) = false
function get(value::T, name::Symbol, default::R)::R where {T, R}
    val = getfield(value, name)::Option{R}
    return not(val) ? default : val
end

"""
    RequestOptions(; chunksize=, connecttimeout=, readtimeout=, tlsconfig=, maxredirects=, allowredirects=)

A type to represent various http request options. Lives as a separate type so that options can be set
at the `HTTP.Client` level to be applied to every request sent. Options include:

  * `chunksize::Int`: if a request body is larger than `chunksize`, the "chunked-transfer" http mechanism will be used and chunks will be sent no larger than `chunksize`
  * `connecttimeout::Float64`: sets a timeout on how long to wait when trying to connect to a remote host; default = 10.0 seconds
  * `readtimeout::Float64`: sets a timeout on how long to wait when receiving a response from a remote host; default = 9.0 seconds
  * `tlsconfig::TLS.SSLConfig`: a valid `TLS.SSLConfig` which will be used to initialize every https connection
  * `maxredirects::Int`: the maximum number of redirects that will automatically be followed for an http request
  * `allowredirects::Bool`: whether redirects should be allowed to be followed at all; default = `true`
  * `forwardheaders::Bool`: whether user-provided headers should be forwarded on redirects; default = `false`
  * `retries::Int`: # of times a request will be tried before throwing an error; default = 3
  * `managecookies::Bool`: whether the request client should automatically store and add cookies from/to requests (following appropriate host-specific & expiration rules)
  * `statusraise::Bool`: whether an `HTTP.StatusError` should be raised on a non-2XX response status code
  * `insecure::Bool`: whether an "https" connection should allow insecure connections (no TLS verification); default = `false`
"""
mutable struct RequestOptions
    chunksize::Option{Int}
    gzip::Option{Bool}
    connecttimeout::Option{Float64}
    readtimeout::Option{Float64}
    tlsconfig::Option{TLS.SSLConfig}
    maxredirects::Option{Int}
    allowredirects::Option{Bool}
    forwardheaders::Option{Bool}
    retries::Option{Int}
    managecookies::Option{Bool}
    statusraise::Option{Bool}
    insecure::Option{Bool}
    RequestOptions(ch::Option{Int}, gzip::Option{Bool}, ct::Option{Float64}, rt::Option{Float64}, tls::Option{TLS.SSLConfig}, mr::Option{Int}, ar::Option{Bool}, fh::Option{Bool}, tr::Option{Int}, mc::Option{Bool}, sr::Option{Bool}, i::Option{Bool}) =
        new(ch, gzip, ct, rt, tls, mr, ar, fh, tr, mc, sr, i)
end

const RequestOptionsFieldTypes = Dict(:chunksize      => Int,
                                      :gzip           => Bool,
                                      :connecttimeout => Float64,
                                      :readtimeout    => Float64,
                                      :tlsconfig      => TLS.SSLConfig,
                                      :maxredirects   => Int,
                                      :allowredirects => Bool,
                                      :forwardheaders => Bool,
                                      :retries        => Int,
                                      :managecookies  => Bool,
                                      :statusraise    => Bool,
                                      :insecure       => Bool)

function RequestOptions(options::RequestOptions; kwargs...)
    for (k, v) in kwargs
        setfield!(options, k, convert(RequestOptionsFieldTypes[k], v))
    end
    return options
end

RequestOptions(chunk=nothing, gzip=nothing, ct=nothing, rt=nothing, tls=nothing, mr=nothing, ar=nothing, fh=nothing, tr=nothing, mc=nothing, sr=nothing, i=nothing; kwargs...) =
    RequestOptions(RequestOptions(chunk, gzip, ct, rt, tls, mr, ar, fh, tr, mc, sr, i); kwargs...)

function update!(opts1::RequestOptions, opts2::RequestOptions)
    for i = 1:nfields(RequestOptions)
        f = fieldname(RequestOptions, i)
        not(getfield(opts1, f)) && setfield!(opts1, f, getfield(opts2, f))
    end
    return opts1
end

# Request
"""
    Request()
    Request(method, uri, headers, body; options=RequestOptions())
    Request(; method=HTTP.GET, uri=HTTP.URI(""), major=1, minor=1, headers=HTTP.Headers(), body="")

A type representing an http request. `method` can be provided as a string or `HTTP.GET` type enum.
`uri` can be provided as an actual `HTTP.URI` or string. `headers` should be provided as a `Dict`.
`body` may be provided as string, byte vector, IO, or `HTTP.FIFOBuffer`.
`options` should be a `RequestOptions` type, see `?HTTP.RequestOptions` for details.

Accessor methods include:
  * `HTTP.method`: method for a request
  * `HTTP.major`: major http version for a request
  * `HTTP.minor`: minor http version for a request
  * `HTTP.uri`: uri for a request
  * `HTTP.headers`: headers for a request
  * `HTTP.body`: body for a request as a `HTTP.FIFOBuffer`

Two convenience methods are provided for accessing a request body:
  * `take!(r)`: consume the request body, returning it as a `Vector{UInt8}`
  * `String(r)`: consume the request body, returning it as a `String`
"""
mutable struct Request
    method::HTTP.Method
    major::Int16
    minor::Int16
    uri::URI
    headers::Headers # includes cookies
    body::Union{FIFOBuffer, Form}
end

# accessors
method(r::Request) = r.method
major(r::Request) = r.major
minor(r::Request) = r.minor
uri(r::Request) = r.uri
headers(r::Request) = r.headers
body(r::Request) = r.body

defaultheaders(::Type{Request}) = Headers(
    "User-Agent" => "HTTP.jl/0.0.0",
    "Accept" => "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8,application/json"
)

function Request(m::HTTP.Method, uri::URI, userheaders::Headers, b;
                    options::RequestOptions=RequestOptions(),
                    verbose::Bool=false,
                    logger::Option{IO}=STDOUT)
    if m != CONNECT
        headers = defaultheaders(Request)
        headers["Host"] = host(uri)
    else
        headers = Headers()
    end
    if !isempty(userinfo(uri)) && !haskey(headers, "Authorization")
        headers["Authorization"] = "Basic $(base64encode(userinfo(uri)))"
        @log "adding basic authentication header"
    end
    if isa(b, Dict) || isa(b, Form)
        # form data
        body = Form(b)
        headers["Content-Type"] = "multipart/form-data; boundary=$(body.boundary)"
    else
        body = FIFOBuffer(b)
    end
    if iscompressed(body) && length(body) > get(opts, :chunksize, 0)
        opts.chunksize = length(body) + 1
    end
    if !haskey(headers, "Content-Type") && length(body) > 0 && !isa(body, Form)
        sn = sniff(body)
        headers["Content-Type"] = sn
        @log "setting Content-Type header to: $sn"
    end
    return Request(m, Int16(1), Int16(1), uri, merge!(headers, userheaders), body)
end

Request(method, uri, h=Headers(), body=""; options::RequestOptions=RequestOptions(), io::IO=STDOUT, verbose::Bool=false) =
    Request(convert(HTTP.Method, method),
            isa(uri, String) ? URI(uri; isconnect=(method == "CONNECT" || method == CONNECT)) : uri,
            h, body; options=options, io=io, verbose=verbose)

Request(; method::Method=GET, major::Integer=Int16(1), minor::Integer=Int16(1), uri=URI(""), headers=Headers(), body=FIFOBuffer("")) =
    Request(method, major, minor, uri, headers, body)

==(a::Request,b::Request) = (a.method    == b.method)    &&
                            (a.major     == b.major)     &&
                            (a.minor     == b.minor)     &&
                            (a.uri       == b.uri)       &&
                            (a.headers   == b.headers)   &&
                            (a.body      == b.body)

Base.showcompact(io::IO, r::Request) = print(io, "Request(\"", resource(r.uri), "\", ",
                                        length(r.headers), " headers, ",
                                        length(r.body), " bytes in body)")

"""
    Response(status::Integer)
    Response(status::Integer, body::String)
    Response(status::Integer, headers, body)
    Response(; status=200, cookies=HTTP.Cookie[], headers=HTTP.Headers(), body="")

A type representing an http response. `status` represents the http status code for the response.
`headers` should be provided as a `Dict`. `body` can be provided as a string, byte vector, IO, or `HTTP.FIFOBuffer`.

Accessor methods include:
  * `HTTP.status`: status for a response
  * `HTTP.statustext`: statustext for a response
  * `HTTP.major`: major http version for a response
  * `HTTP.minor`: minor http version for a response
  * `HTTP.cookies`: cookies for a response, returned as a `Vector{HTTP.Cookie}`
  * `HTTP.headers`: headers for a response
  * `HTTP.request`: the `HTTP.Request` that resulted in this response
  * `HTTP.history`: history for a response if redirects were followed from an original request
  * `HTTP.body`: body for a response as a `HTTP.FIFOBuffer`

Two convenience methods are provided for accessing a response body:
  * `take!(r)`: consume the response body, returning it as a `Vector{UInt8}`
  * `String(r)`: consume the response body, returning it as a `String`
"""
mutable struct Response
    status::Int32
    major::Int16
    minor::Int16
    cookies::Vector{Cookie}
    headers::Headers
    body::FIFOBuffer
    request::Nullable{Request}
    history::Vector{Response}
end

# accessors
status(r::Response) = r.status
major(r::Response) = r.major
minor(r::Response) = r.minor
cookies(r::Response) = r.cookies
headers(r::Response) = r.headers
request(r::Response) = r.request
history(r::Response) = r.history
statustext(r::Response) = Base.get(STATUS_CODES, r.status, "Unknown Code")
body(r::Union{Request, Response}) = r.body
Base.take!(r::Union{Request, Response}) = readavailable(body(r))
Base.String(r::Union{Request, Response}) = String(body(r))

Response(; status::Int=200,
         cookies::Vector{Cookie}=Cookie[],
         headers::Headers=Headers(),
         body::FIFOBuffer=FIFOBuffer(""),
         request::Nullable{Request}=Nullable{Request}(),
         history::Vector{Response}=Response[]) =
    Response(status, Int16(1), Int16(1), cookies, headers, body, request, history)

Response(n::Integer, r::Request) = Response(; body=FIFOBuffer(n), request=Nullable(r))
Response(s::Integer) = Response(; status=s)
Response(s::Integer, msg) = Response(; status=s, body=FIFOBuffer(msg))
Response(b::Union{Vector{UInt8}, String}) = Response(; headers=defaultheaders(Response), body=FIFOBuffer(b))
Response(s::Integer, h::Headers, body) = Response(; status=s, headers=h, body=FIFOBuffer(body))

defaultheaders(::Type{Response}) = Headers(
    "Server"            => "Julia/$VERSION",
    "Content-Type"      => "text/html; charset=utf-8",
    "Content-Language"  => "en",
    "Date"              => Dates.format(now(Dates.UTC), Dates.RFC1123Format)
)

==(a::Response,b::Response) = (a.status  == b.status)  &&
                              (a.major   == b.major)   &&
                              (a.minor   == b.minor)   &&
                              (a.headers == b.headers) &&
                              (a.cookies == b.cookies) &&
                              (a.body    == b.body)

function Base.showcompact(io::IO, r::Response)
    print(io, "Response(", r.status, " ", Base.get(STATUS_CODES, r.status, "Unknown Code"), ", ",
          length(r.headers)," headers, ",
          length(r.body)," bytes in body)")
end

## Request & Response writing
# start lines
function startline(io::IO, r::Request)
    res = resource(uri(r); isconnect=r.method == CONNECT)
    res = ifelse(res == "", "/", res)
    write(io, "$(r.method) $res HTTP/$(r.major).$(r.minor)$CRLF")
end

function startline(io::IO, r::Response)
    write(io, "HTTP/$(r.major).$(r.minor) $(r.status) $(statustext(r))$CRLF")
end

# headers
function headers(io::IO, r::Union{Request, Response})
    for (k, v) in headers(r)
        write(io, "$k: $v$CRLF")
    end
    # write(io, CRLF); we let the body write this in case of chunked transfer
end

# body
# https://tools.ietf.org/html/rfc7230#section-3.3
function hasmessagebody(r::Response)
    if 100 <= status(r) < 200 || status(r) == 204 || status(r) == 304
        return false
    elseif !Base.isnull(request(r))
        req = Base.get(request(r))
        method(req) in (HEAD, CONNECT) && return false
    end
    return true
end
hasmessagebody(r::Request) = length(r.body) > 0 && !(r.method in (GET, HEAD, CONNECT))

function body(io::IO, r::Union{Request, Response}, opts)
    if !hasmessagebody(r)
        write(io, "$CRLF")
        return
    end
    chksz = get(opts, :chunksize, typemax(Int))
    pos = position(r.body)
    @sync begin
        @async begin
            chunked = false
            bytes = UInt8[]
            while !eof(r.body)
                bytes = read(r.body, chksz)
                eof(r.body) && !chunked && break
                if !chunked
                    write(io, "Transfer-Encoding: chunked$CRLF$CRLF")
                end
                chunked = true
                chunk = length(bytes)
                chunk == 0 && break
                write(io, "$(hex(chunk))$CRLF")
                write(io, bytes, CRLF)
            end
            if chunked
                write(io, "$(hex(0))$CRLF$CRLF")
            else
                write(io, "Content-Length: $(dec(length(bytes)))$CRLF$CRLF")
                write(io, bytes)
            end
        end
    end
    seek(r.body, pos)
    return
end

Base.write(io::IO, r::Union{Request, Response}, opts) = write(io, string(r))
function Base.string(r::Union{Request, Response}, opts=RequestOptions())
    i = IOBuffer()
    startline(i, r)
    headers(i, r)
    body(i, r, opts)
    return String(take!(i))
end

function Base.show(io::IO, r::Union{Request,Response}; opts=RequestOptions())
    println(io, typeof(r), ":")
    println(io, "\"\"\"")
    startline(io, r)
    headers(io, r)
    buf = IOBuffer()
    if isopen(r.body)
        println(io, "\n[open HTTP.FIFOBuffer with $(length(r.body)) bytes to read]")
    else
        body(buf, r, opts)
        b = take!(buf)
        if length(b) > 2
            contenttype = sniff(b)
            if contenttype in DISPLAYABLE_TYPES
                if length(b) > 750
                    println(io, "\n[$(typeof(r)) body of $(length(b)) bytes]")
                    println(io, String(b)[1:750])
                    println(io, "⋮")
                else
                    print(io, String(b))
                end
            else
                contenttype = Base.get(r.headers, "Content-Type", contenttype)
                encoding = Base.get(r.headers, "Content-Encoding", "")
                encodingtxt = encoding == "" ? "" : " with '$encoding' encoding"
                println(io, "\n[$(length(b)) bytes of '$contenttype' data$encodingtxt]")
            end
        else
            print(io, String(b))
        end
    end
    print(io, "\"\"\"")
end


text:
@testset "1317.convert-integer-to-the-sum-of-two-no-zero-integers.jl" begin
    @test get_no_zero_integers(1010) == [11, 999]
    @test get_no_zero_integers(10000) == [1, 9999]
end


text:
module ValueInferenceTools

export TaskStats, InferenceAccumulator, ValueAccumulator, Simulator
export accumulatefromdata, inferencesimulator, valuesimulator

using StatsFuns: logistic

include("probabilities.jl")
include("accumulate.jl")
include("agent.jl")

end # module


text:
using Clustering
using DataStructures

function filter_components(input, min_size)
    output = zeros(Int64, size(input))
    filter_components!(output, input, min_size)
end

function filter_components!(input, min_size)
    filter_components!(input, input, min_size)
end

function connected_components(input)
    output = zeros(UInt8, size(input))
    components!(output, input)
end

function components(input; radius=5.0, min_neighbors=8, min_cluster_size=12)
    blobs = Dict{Int, Array{Int}}()
    fg_points = hcat(map(x -> collect(convert.(Float64, ind2sub(input, x))), find(input))...)
    if size(fg_points, 2) < min_cluster_size
        return blobs
    end
    C = dbscan(fg_points, radius; min_neighbors=min_neighbors, min_cluster_size=min_cluster_size)
    for (label, x) in enumerate(C)
        blobs[label] = (convert.(Int, fg_points[:, x.core_indices]))'
    end
    return blobs
end

function connected_components!(output, input)
    (sets, sizes) = first_pass_label_components!(output, input)
    labels = Dict{Int, Int}(-1 => 0)
    blobs = Dict{Int, Array{Int}}()
    for idx in CartesianRange(size(input))
        if !is_fg(input[idx])
            continue
        end
        component = output[idx]
        root = find_root(sets, component)
        label = get_label!(labels, root)
        index_vec = [idx[n] for n in 1:length(idx)]
        if haskey(blobs, label)
            blobs[label] = append!(blobs[label], index_vec)
        else
            blobs[label] = index_vec
        end
    end
    for (key, blob) in blobs
        blobs[key] = reshape(blob, 2, :)'
    end
    return blobs
end

function filter_components!(output, input, min_size)
    (sets, sizes) = first_pass_label_components!(output, input)

    # avoid modifying counter while iterating over it
    size_list = [(key, count) for (key, count) in sizes]
    for (key, count) in size_list
        root = find_root(sets, key)
        push!(sizes, root, reset!(sizes, key))
    end
    labels = Dict{Int, Int}(-1 => 0)
    for idx in CartesianRange(size(input))
        if !is_fg(input[idx])
            continue
        end
        component = output[idx]
        root = find_root(sets, component)
        component_size = sizes[root]
        if component_size < min_size
            output[idx] = 0
        else
            label = get_label!(labels, root)
            output[idx] = label
        end
    end
    return output
end

function first_pass_label_components!(output, input)
    dims = length(size(input))
    sets = IntDisjointSets(0)
    sizes = counter(Int64)
    min_bound = CartesianIndex(ones(Int, dims)...)
    myEye = eye(Int64, dims)
    offsets = [CartesianIndex(myEye[x, :]...) for x in 1:dims]
    for idx in CartesianRange(size(input))
        if !is_fg(input[idx])
            continue
        end
        components = neighbor_components(idx, output, offsets, min_bound)
        if length(components) == 0
            component = push!(sets)
            output[idx] = component
            push!(sizes, component)
        else
            output[idx] = components[1]
            push!(sizes, components[1])
            if length(components) == 2
                union!(sets, components...)
            end
        end
    end
    return (sets, sizes)
end

is_fg(x) = x > zero(x)

function get_label!(labels, key)
    if haskey(labels, key)
        return labels[key]
    else
        labels[-1] += 1
        label = labels[-1]
        labels[key] = label
        return label
    end
end

function neighbor_components(idx, labels, offsets, min_bound)
    components = []
    for offset in offsets
        neighbor = idx - offset
        if max(neighbor, min_bound) == neighbor && labels[neighbor] > 0
            components = vcat(components, [labels[neighbor]])
        end
    end
    components
end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_weak_blast_wave

surface_flux = flux_ranocha
volume_flux = flux_ranocha
polydeg = 4
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=1.0,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral)

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

coordinates_min = (-1.0, -1.0)
coordinates_max = ( 1.0,  1.0)

trees_per_dimension = (4, 4)
mesh = P4estMesh(trees_per_dimension,
                 polydeg=4, initial_refinement_level=2,
                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                 periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_solution,
                        stepsize_callback)
###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# %%
using Plots
# %%
include("../src/Laser.jl")
# %%
function dfdt(f, dt)
    n = length(f)
    g = copy(f)
    g[1] = (f[2] - f[1]) / dt
    g[n] = (f[n] - f[n - 1]) / dt
    for i in 2:n - 1
        g[i] = (f[i + 1] - f[i - 1]) / 2dt
    end
    return g
end
 # %%
pulseEt = Laser.Pulse(1, 800, 0, :sin4, 10, :Et)
pulseAt = Laser.Pulse(1, 800, 0, :sin4, 10, :At)
pulse = pulseAt
time = Laser.Time(pulse, 1000, 10)
# %%
Et = Laser.getEt(time.t, pulse)
# %%
At = Laser.getAt(time.t, pulse)
# %%
dAdt = dfdt(At, time.dt)
# %%
p1 = plot(time.t, Et, label="Et", color=:black, linewidth=3.0)
p1 = plot!(time.t, -dAdt, label="-dAdt")
p2 = plot(time.t, At, label="At")
plot(p1, p2, layout=(2, 1))
# %%


text:
export hofbauerextension

"""
    HofbauerDomain

A `Domain` object annotated with a depth index for `HofbauerExtensions`.
It is not a `Domain` subtype itself as it should only be used in the context of Hofbauer extensions.
"""
struct HofbauerDomain{D<:Domain}
    domain::D
    depth::Int
end
show(io::IO, d::HofbauerDomain) = print(io, "$(d.domain) at depth $(d.depth)")
convert(::Type{Domain},d::HofbauerDomain) = d.domain

maxrangedomain(m::AbstractIntervalMap) =
    rangedomain(branches(m)[findmax(arclength.(rangedomain.(branches(m))))[2]])

"""
    HofbauerPoint(x, graphind)

Constructs a `HofbauerPoint`, representing point `x` on the `graphind`th domain
of a Hofbauer extension.
"""
struct HofbauerPoint{T}
    x::T
    graphind::Int
end
show(io::IO,hp::HofbauerPoint) = print(io,"$(hp.x) on domain #$(hp.graphind)")

# mapping Hofbauer extension

"""
    maphb(b::AbstractBranch, x)

Maps point `x` according to `b`. For `NeutralBranches` calculates the return map out of `domain(b)`.
"""
maphb(b::AbstractBranch,x) = b(x)
maphb(b::NeutralBranch,x) = _NI("maphb for NeutralBranch") # use domain(b) to fix.
for (M,MHB) = ((:mapD,:maphbD),(:mapP,:maphbP),(:mapinv,:maphbinv),
        (:mapinvD,:maphbinvD),(:mapinvP,:maphbinvP))
    @eval $MHB(b::AbstractBranch,x) = $M(b,x)
    @eval $MHB(b::NeutralBranch,x) = _NI("$MHB for NeutralBranch")
end

function _hofbauerdestination(he::HofbauerExtension, bind, x, graphind)
    graphinds = searchsorted(he.fedgelist[graphind], HofbauerEdge(0,0,bind,false), by=branchind)
    isempty(graphinds) && error("reached end of generated Hofbauer extension at depth $(hdomains(im.he)[graphind].depth)")
    for g in graphinds
        dst_ind = dst(im.he.fedgelist[graphind][g])
        if x ∈ convert(Domain,hdomains(he)[dst_ind])
            graphind = dst_ind
            break
        end
    end
    graphind
end
@inline _hofbauerdestination(he::HofbauerExtension, bind, p) =
    _hofbauerdestination(he, bind, p.x, p.graphind)

function (he::HofbauerExtension)(p::HofbauerPoint)
    bind = getbranchind(he.m, p.x)
    fx = maphb(branches(he.m)[bind], p.x)
    HofbauerPoint(x, _hofbauerdestination(he, bind, fx))
end
function mapP(he::HofbauerExtension, p::HofbauerPoint)
    bind = getbranchind(he.m, p.x)
    fx, dfx = maphbP(branches(he.m)[bind], p.x)
    HofbauerPoint(fx,graphind), dfx
end
mapD(he::HofbauerExtension, p::HofbauerPoint) = maphbD(he.m, p.x)

_hofbauersource(he::HofbauerExtension, bind, graphind) =
    searchsortedfirst(he.bedgelist[graphind], HofbauerEdge(0,0,bind,false), by=branchind)

# TODO: assert in rangedomain
mapinv(he::HofbauerExtension, i, p::HofbauerPoint) =
    HofbauerPoint(maphbinv(branches(he.m)[i],p.x), _hofbauersource(he,i,p.graphind))
mapinvD(he::HofbauerExtension, i, p::HofbauerPoint) =
        maphbinvD(branches(he.m)[i], p.x)
function mapinvP(he::HofbauerExtension, i, p::HofbauerPoint)
    (vx, dvx) = maphbinvP(branches(he.m)[i], p.x)
    HofbauerPoint(vx, _hofbauersource(he,i,p.graphind)), dvx
end

function mapinv(he::HofbauerExtension,e::HofbauerEdge,p::HofbauerPoint)
    @assert dst(e) == p.graphind
    HofbauerPoint(maphbinv(branches(he.m)[branchind(e)],p.x),src(e))
end
function mapinvD(he::HofbauerExtension,e::HofbauerEdge,p::HofbauerPoint)
    @assert dst(e) == p.graphind
    maphbinvD(branches(he.m)[branchind(e)],p.x)
end
function mapinvP(he::HofbauerExtension,e::HofbauerEdge,p::HofbauerPoint)
    @assert dst(e) == p.graphind
    vx,dvx = maphbinvP(branches(he.m)[branchind(e)],p.x)
    HofbauerPoint(vx,src(e)), dvx
end

"""
    hofbauerextension(m,basedomains;maxdepth=100,forcereturn=trues(length(basedomains)))

Generates a Hofbauer extension of `m` starting from `basedomains`, which may be
a single domain.

The keyword `maxdepth` says how deep the Hofbauer extension should go, and
`forcereturn` sets whether a return to given members of `basedomains` should
 forced if possible.

 For example, in the case of `f(x) = mod(2x,1)`, if `basedomains` is set to `Interval(0.,0.3)` then
`Interval(0.,0.5)` might map to `Interval(0.,1.)`` for a given map if `forcereturn=false` but would map to
Interval(0,0.3)` and `Interval(0.3,1)` if `forcereturn=true`.
"""
function hofbauerextension(m::AbstractIntervalMap,basedoms=(maxrangedomain(m),);maxdepth=100,forcereturn=trues(length(basedoms)))
    @assert domain(m) == rangedomain(m)

    G = HofbauerExtension{Int,HofbauerDomain}(m)
    for (n,bd) in enumerate(basedoms)
        basehdom = HofbauerDomain(convert(Domain,bd),1)
        add_vertex!(G,basehdom;forcereturn=forcereturn[n])
    end

    i = 1
    while nv(G) >= i && hdomains(G)[i].depth < maxdepth
        towerwalk!(G,i,hdomains(G)[i])
        i += 1
    end

    D = promote_type(collect(typeof(h.domain) for h in hdomains(G))...)
    # HofbauerExtension(ne(G),fedge(G),bedge(G),getmap(G),
    #     convert(Vector{HofbauerDomain{D}},hdomains(G)))
    G
end
hofbauerextension(m::AbstractIntervalMap,basedoms::Domain;maxdepth=100,forcereturn=trues(length(basedoms))) =
    hofbauerextension(m,(basedoms,);maxdepth=maxdepth,forcereturn=forcereturn)

# walks through domains in HofbauerExtension
function towerwalk!(G::HofbauerExtension, graphind, hdom)
    @assert hdomains(G)[graphind] == hdom # putting in hdom for type inference
    m = getmap(G)
    for branchind in eachbranchindex(m)
        b = branches(m)[branchind]
        bdom = domain(b)
        capdom = bdom ∩ convert(Domain,hdom)
        if ~isempty(capdom) && arclength(capdom) > 0 # hdom is partially contained in the domain of b
            if isa(b,NeutralBranch) #&& (neutralfixedpoint(b) ∈ capdom) # inducing at a neutral fixed point
                isnfp = true
                newdom = rangedomain(b) \ (atol_isapprox(bdom,capdom) ? domain(b) : b(capdom))
            else # normal bounce inducing
                isnfp = false
                newdom = atol_isapprox(bdom,capdom) ? rangedomain(b) : b(capdom)
            end
            for i in G.returnto
                forcedom = convert(Domain,hdomains(G)[i])
                if approx_issubset(forcedom,newdom)
                    towerextend!(G, graphind, hdom, forcedom, branchind, isnfp)
                    newdom = newdom \ forcedom
                end
            end
            arclength(newdom) > 0 && towerextend!(G, graphind, hdom, newdom, branchind, isnfp)
        end
    end
end

towerextend!(G::HofbauerExtension, graphind, hdom, newdom, branchind, isnfp) =
    isempty(newdom) ? true : actual_towerextend!(G::HofbauerExtension, graphind, hdom, newdom, branchind, isnfp)

function actual_towerextend!(G::HofbauerExtension, graphind, hdom, newdom, branchind, isnfp)
    arclength(newdom) < 300eps(arclength(convert(Domain,hdom))) && return true
    newgraphind = findfirst(collect(atol_isapprox(h.domain,newdom) for h in hdomains(G)))

    # add new domain if it doesn't already exist
    if newgraphind isa Nothing
        newhdom = HofbauerDomain(newdom,hdom.depth+1)
        add_vertex!(G,newhdom)
        newgraphind = nv(G)
    end

    # adding new edge
    newedge = HofbauerEdge(graphind, newgraphind, branchind, isnfp)
    add_edge!(G,newedge)
end
function actual_towerextend!(G::HofbauerExtension, graphind, hdom, newdom::UnionDomain, branchind, isnfp)
    for d in newdom.domains
        towerextend!(G, graphind, hdom, d, branchind, isnfp)
    end
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Compositions

using Nemo, IterTools, Triangles

export ModuleCompositions
export I097805, L097805, V097805, M097805

"""
* I097805, L097805, V097805, M097805
"""
const ModuleCompositions = ""

"""
Recurrence for `A097805`, the compositions of ``n`` with ``k`` parts.
"""
function R097805(n, k, prevrow::Function)
    k == 0 && return ZZ(k^n)
    prevrow(k - 1) + prevrow(k)
end

"""
Iterates over the first ``n`` rows of A097805.
"""
I097805(n) = RecTriangle(n, R097805)

"""
Lists the first ``n`` rows of A097805 by concatinating. This is the format for submissions to the OEIS.
"""
L097805(n) = vcat(I097805(n)...)

"""
Return the triangular array as a square matrix.
"""
M097805(dim) = fromΔ(L097805(dim))

"""
Return row ``n`` of A097805 based on the iteration I097805(n).
"""
V097805(n) = nth(I097805(n+1), n+1)

#START-TEST-########################################################

# Return row n of A097805 based on a closed formula.
VN097805(n) = n == 0 ? fmpz[1]   : [Nemo.binom(n-1, k-1)       for k in 0:n]
VJ097805(n) = n == 0 ? BigInt[1] : [binomial(BigInt(n-1), k-1) for k in 0:n]

function test()
end

function demo()
    println("\nIterates over the first n rows of A097805.")
    for r in I097805(9) println(r) end

    println("\nLists the first 9 rows of L097805 by concatinating.")
    println(L097805(9))

    #ShowAsMatrix(L097805(9))
    #println(M097805(9))

    println("\nReturns row 8 of A097805 based on the iteration I097805(8).")
    println(V097805(8))

    println("\nReturns row 8 of A097805 based on a closed formula.")
    println(VN097805(8))

    println("\nReturns row n of A097805 based on the iteration I097805(n).")
    for n in 0:8 println(n, ": ", V097805(n)) end

    println("\nBenchmark the construction of the first 500 rows of A097805 based on Iteration.")
    #  0.111066 seconds (253.51 k allocations: 7.795 MiB)
    GC.gc()
    @time L097805(500)

    println("\nBenchmark the construction of the first 500 rows of A097805 based on closed formula.")
    # Result: 0.479466 seconds (1.72 M allocations: 21.592 MiB)
    FList(n) = vcat([VN097805(k) for k in 0:n-1]...)
    GC.gc()
    @time FList(500)

    println("\nBenchmark the construction of the first 500 rows of A097805 based on closed formula.")
    # 0.476226 seconds (1.72 M allocations: 21.592 MiB)
    FList(n) = vcat([VJ097805(k) for k in 0:n-1]...)
    GC.gc()
    @time FList(500)

    println("\nResume: Construction based on Iteration is about 4 times faster")
    println("than based on the closed formula and uses much less memory.")
end

"""
L097805(500) :: 0.111066 seconds (253.51 k allocations: 7.795 MiB)
"""
function perf()
    GC.gc()
    @time L097805(500)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
export report

"""
    runCmd(cmd::Cmd, cmdOutput::IO)::IO
Run the external command `cmd`, which will write output to `cmdOutput`.
The stream that's returnd can be written to to provide inoput to the
command.  The second return value is the stderr stream.
"""
function runCmd(cmd::Cmd, cmdOutput::IO)
	cmdInput = Pipe()
	err = Pipe()
	proc = Base.run(pipeline(cmd,
			stdin=cmdInput,
			stdout=cmdOutput,
			stderr=err),
		wait=false)
	process_running(proc) || throw(Exception("Problem starting $cmd"))
	close(cmdInput.out)
	close(err.in)
	return cmdInput, err
end

"""
    report(::Searcher)
Generate an HTML fragment that provides a detailed report of
our cut search and results.
"""
function report(searcher::Searcher;
                includeCutDiagram=true,
                includeCutGraph=false)
    io = IOBuffer()
    function elt(f, io, tagname; attrs...)
        NativeSVG.element(f, tagname, io; attrs...)
    end
    elt(io, :div) do
        elt(io, :h2) do
            write(io, "Panel Cut Report")
        end
        elt(io, :p) do
            write(io, "Report of what stock to purchase and what" *
                " cuts to make to get panels of these aizes")
        end
        function th(io, heading)
            elt(io, :th) do
                write(io, heading)
            end
        end
        function td(io, val; attrs...)
            elt(io, :td; attrs...) do
		if val isa String
		    write(io, val)
		else
                    show(io, val)
		end
            end
        end
	# Table of wanted panels:
        elt(io, :table) do
            elt(io, :thread) do
                elt(io, :tr) do
                    for heading in ("Label", "Length", "Width", "Ok to Flip?")
                        th(io, heading)
                    end
                end
            end
            elt(io, :tbody) do
                for panel in searcher.wanted
                    elt(io, :tr) do
                        td(io, panel.label; align="center")
                        td(io, panel.length; align="right")
                        td(io, panel.width; align="right")
                        td(io, if panel isa FlippedPanel
                               "yes"
                           else
                               "no"
			   end;
			   align="center")
                    end
                end
            end
        end
        if searcher.cheapest == nothing
            elt(io, :p; style="font-weight: bold") do
                write(io, "No solution has been found!")
            end
            return
        else
            elt(io, :p) do
                write(io, "The best solution has a cost of " *
                    "$(searcher.cheapest.accumulated_cost).")
            end
        end
        # Table of panel areas
        elt(io, :div) do
            elt(io, :table) do
                elt(io, :thread) do
                    elt(io, :tr) do
                        th(io, "Group")
                        th(io, "Panel")
                        th(io, "Length")
                        th(io, "Width")
                        th(io, "Area")
                        th(io, "%")
                    end
                end
                elt(io, :tbody) do
                    bought_area = sum(area.(searcher.cheapest.bought))
                    function panel_group(group, panels)
                        for p in panels
                            elt(io, :tr; style=style) do
                                if p == first(panels)
                                    td(io, group;
                                       align="center",
                                       valign="top",
                                       rowspan=length(panels))
                                end
                                td(io, p.label; align="center")
                                td(io, p.length; align="right")
                                td(io, p.width; align="right")
                                td(io, area(p); align="right")
                                td(io, @sprintf("%.2f%%",
						100 * convert(Float64, area(p)/bought_area));
				   align="right")
                                if p == panels[1]
				    frac = convert(Float64, sum(area.(panels)) / bought_area)
				    td(io, @sprintf("%.2f%%", 100 * frac);
				       rowspan=length(panels),
				       valign="top",
				       align="right")
                                end
                            end
                        end
                    end
                    panel_group("bought", searcher.cheapest.bought)
                    panel_group("scrapped", searcher.cheapest.scrapped)
                    panel_group("left over", searcher.cheapest.working)
                    panel_group("finished", searcher.cheapest.finished)
                end
            end
        end
        if includeCutDiagram
            elt(io, :div) do
                toSVG(io, searcher.cheapest)
            end
        end
        if includeCutGraph
            elt(io, :div) do
		# Run the GraphViz dot command, inlining the SVG output
		# into the report:
                dot, err = runCmd(`dot -Tsvg`, io)
                dotgraph(dot, PanelCutGraph(searcher.cheapest))
                close(dot)
                err = read(err)
                if length(err) > 0
                    throw(Exception("Error running dot: $err"))
                end
            end
        end
    end
    return HTML(String(take!(io)))
end


text:
@testset "CHELSEA BioClim" begin
    using SimpleSDMDataSources: rasterurl

    @test rastername(CHELSA{BioClim}, 5) == "CHELSA_bio10_05.tif"

    bioclim_path = joinpath(ENV["ECODATASOURCES_PATH"], "CHELSA/BioClim")
    @test rasterpath(CHELSA{BioClim}) == bioclim_path
    @test rasterpath(CHELSA{BioClim}, 5) == joinpath(bioclim_path, "CHELSA_bio10_05.tif")

    @test rasterurl(CHELSA) == "ftp://envidatrepo.wsl.ch/uploads/chelsa/"
    @test rasterurl(CHELSA{BioClim}, 5) == "ftp://envidatrepo.wsl.ch/uploads/chelsa/chelsa_V1/climatologies/bio/CHELSA_bio10_05.tif"

    download_raster(CHELSA{BioClim}; layer=5)
    @test isfile(joinpath(bioclim_path, "CHELSA_bio10_05.tif"))
end


text:
"""
    struct ElectricPotential{T, N, S, AT} <: AbstractArray{T, N}

Electric potential of the simulation in units of volt (V).

## Parametric types
* `T`: Element type of `data`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{T, N}`: Array containing the values of the electric potential at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points for which the electric potential is determined.
"""
struct ElectricPotential{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S, AT}
end

@inline size(epot::ElectricPotential{T, N, S}) where {T, N, S} = size(epot.data)
@inline length(epot::ElectricPotential{T, N, S}) where {T, N, S} = length(epot.data)
@inline getindex(epot::ElectricPotential{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(epot.data, I...)
@inline getindex(epot::ElectricPotential{T, N, S}, i::Int) where {T, N, S} = getindex(epot.data, i)
@inline getindex(epot::ElectricPotential{T, N, S}, s::Symbol) where {T, N, S} = getindex(epot.grid, s)

# """
#     ElectricPotential(pss::PotentialSimulationSetup{T, 3, Cylindrical} ; kwargs...)::ElectricPotential{T, 3, Cylindrical}
#
# Extracts the electric potential from the `PotentialSimulationSetup`, `pss`, and extrapolate it to an 2π grid.
#
# For 2D grids (r and z) the user has to set the keyword `n_points_in_φ::Int`, e.g.: `n_points_in_φ = 36`.
# """
function ElectricPotential(pss::PotentialSimulationSetup{T, 3, Cylindrical} ; kwargs...)::ElectricPotential{T, 3, Cylindrical} where {T}
    return get_2π_potential(ElectricPotential{T, 3, Cylindrical}(pss.potential, pss.grid); kwargs...)
end

# """
#     ElectricPotential(pss::PotentialSimulationSetup{T, 3, Cartesian} ; kwargs...)::ElectricPotential{T, 3, Cartesian}
#
# Extracts the electric potential from the `PotentialSimulationSetup`, `pss`.
# """
function ElectricPotential(pss::PotentialSimulationSetup{T, 3, Cartesian} )::ElectricPotential{T, 3, Cartesian} where {T}
    return ElectricPotential{T, 3, Cartesian}(pss.potential, pss.grid)
end

function NamedTuple(epot::ElectricPotential{T, 3}) where {T}
    return (
        grid = NamedTuple(epot.grid),
        values = epot.data * u"V",
    )
end
Base.convert(T::Type{NamedTuple}, x::ElectricPotential) = T(x)

function ElectricPotential(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    ElectricPotential{T, N, S, typeof(grid.axes)}( ustrip.(uconvert.(u"V", nt.values)), grid)
end
Base.convert(T::Type{ElectricPotential}, x::NamedTuple) = T(x)


text:
module FunRefs

using Comm, Foldable, Functor, Funs, GIDs, Monad

import Base: show
import Base: eltype, getindex, setindex!, ==, isready, wait
import Base: serialize, deserialize
import Base: length, start, next, done, isempty

import GIDs: islocal

export FunRef
export proc, islocal
export set!, set_from_ref!   # reset!

immutable FunRefRemoteReady end
immutable FunRefRemoteUnready end

# TODO: Replace ready,cond by Maybe{Condition}

type FunRef{T}
    ready::Bool
    gid::GID                    # if ready
    obj::T                      # cache, if ready and local
    cond::Condition             # while not ready

    # Create unready
    function FunRef()
        ref = new(false)
        ref.cond = Condition()
        ref
    end
    # Create from local object
    function FunRef(obj::T)
        # gid = makegid(obj)
        gid = nullgid()
        ref = new(true, gid, obj)
        ref
    end
    # Create from ready remote reference
    function FunRef(::FunRefRemoteReady, gid::GID, prevgid::GID)
        @assert !isnull(gid) && !isnull(prevgid)
        ref = new(true, gid)
        finalizer(ref, finalize)
        if islocal(gid) ref.obj = getobj(gid) end
        take_ownership(ref, prevgid)
        ref
    end
    # Create from unready remote reference
    function FunRef(::FunRefRemoteUnready, prevgid::GID)
        @assert !isnull(prevgid)
        ref = new(false)
        ref.cond = Condition()
        wait_for_ref(ref, prevgid)
        ref
    end
end

FunRef{T}(obj::T) = FunRef{T}(obj)

# Do not free GIDs during garbage collection. Instead, put them on a to-free
# list, and have a worker thread periodically free the GIDs on this list.
global GIDLIST = GID[]
function mark_gid(gid::GID)
    global GIDLIST::Array{GID,1}
    push!(GIDLIST, gid)
end
function free_gids()
    global GIDLIST::Array{GID,1}
    gidlist = GID[]
    while true
        yield()
        gidlist, GIDLIST = GIDLIST, gidlist
        map(free, gidlist)
        empty!(gidlist)
    end
end
@schedule free_gids()

function finalize(ref::FunRef)
    # if !isready(ref) return end
    # if isnull(ref.gid) return end
    @assert isready(ref)
    @assert !isnull(ref.gid)
    # free(ref.gid)
    mark_gid(ref.gid)
end
function ensure_have_gid(ref::FunRef)
    if !isnull(ref.gid) return end
    ref.gid = makegid(ref.obj)
    finalizer(ref, finalize)
end
function take_ownership(ref::FunRef, prevgid::GID)
    @assert !isnull(ref.gid) && !isnull(prevgid)
    gid = ref.gid
    refgid = makegid(ref)
    rexec(gid.proc) do
        take_ownership2(gid, prevgid, refgid)
    end
end
function take_ownership2(gid::GID, prevgid::GID, refgid::GID)
    alloc(gid)
    free(prevgid)
    free(refgid)
end
function wait_for_ref(ref::FunRef, prevgid::GID)
    @assert !isnull(prevgid)
    refgid = makegid(ref)
    rexec(prevgid.proc) do
        wait_for_ref2(prevgid, refgid)
    end
end
function wait_for_ref2(prevgid::GID, refgid::GID)
    prev = getobj(prevgid)::FunRef
    wait(prev)
    ensure_have_gid(prev)
    rexec(refgid.proc) do
        wait_for_ref3(refgid, prev.gid, prevgid)
    end
end
function wait_for_ref3(refgid::GID, gid::GID, prevgid::GID)
    ref = getobj(refgid)::FunRef
    set!(ref, gid, prevgid)
end

function setindex!{T}(ref::FunRef{T}, obj)
    @assert !isready(ref)
    # gid = makegid(obj)
    gid = nullgid()
    ref.gid = gid
    ref.obj = obj
    ref.ready = true
    notify(ref.cond)
end
function set!{T}(ref::FunRef{T}, gid::GID, prevgid::GID)
    @assert !isready(ref)
    @assert !isnull(gid) && !isnull(prevgid)
    ref.gid = gid
    finalizer(ref, finalize)
    if islocal(gid) ref.obj = getobj(gid) end
    ref.ready = true
    notify(ref.cond)
    take_ownership(ref, prevgid)
end
# This could have a more generic signature
function set_from_ref!{T}(ref::FunRef{T}, other::FunRef{T})
    if islocal(other)
        ref[] = other[]
    else
        wait(other)
        othergid = makegid(other)
        set!(ref, other.gid, othergid)
    end
end
# Don't know how to unset ref.obj
# function reset!(ref::FunRef)
#     @assert isready(ref)
#     if islocal(ref.gid) ref.obj = nothing end
#     free(ref.gid)
#     ref.ready = false
#     ref.cond = Condition()
# end

function show(io::IO, ref::FunRef)
    if !isready(ref)
        print(io, "$(typeof(ref))($(ref.ready))")
    elseif !islocal(ref)
        print(io, "$(typeof(ref))($(ref.ready),$(ref.gid))")
    else
        print(io, "$(typeof(ref))($(ref.ready),$(ref.gid),$(ref.obj))")
    end
end

eltype{T}(ref::FunRef{T}) = T

function isready(ref::FunRef)
    ref.ready
end
function wait(ref::FunRef)
    if isready(ref) return end
    # Note: There must be no interruption between the "isready" and
    # the "wait"; in particular, there can be no I/O here.
    wait(ref.cond)
    # while !isready(ref)
    #     sleep(0.01)
    # end
end

function proc(ref::FunRef)
    wait(ref)
    if isnull(ref.gid) return Comm.myproc() end
    ref.gid.proc
end
function islocal(ref::FunRef)
    # proc(ref) == Comm.myproc()
    wait(ref)
    islocal(ref.gid)
end
function getindex{T}(ref::FunRef{T})
    wait(ref)
    @assert islocal(ref)
    ref.obj::T
end

# function =={T}(ref1::FunRef{T}, ref2::FunRef{T})
#     wait(ref1, ref2)
#     ref1.gid == ref2.gid
# end

function serialize{T}(s::AbstractSerializer, ref::FunRef{T})
    Base.serialize_type(s, FunRef{T})
    write(s.io, ref.ready)
    if ref.ready
        ensure_have_gid(ref)
        write(s.io, ref.gid)
        refgid = makegid(ref)
        write(s.io, refgid)
     else
        refgid = makegid(ref)
        write(s.io, refgid)
    end
end
function deserialize{T}(s::AbstractSerializer, ::Type{FunRef{T}})
    ready = read(s.io, Bool)
    if ready
        gid = read(s.io, GID)
        prevgid = read(s.io, GID)
        FunRef{T}(FunRefRemoteReady(), gid, prevgid)
    else
        prevgid = read(s.io, GID)
        FunRef{T}(FunRefRemoteUnready(), prevgid)
    end
end

length(ref::FunRef) = 1
start(ref::FunRef) = false
next(ref::FunRef, i) = ref[], true
done(ref::FunRef, i) = i
isempty(ref::FunRef) = false
getindex(ref::FunRef, i::Integer) = (@assert i==1; ref[])
setindex!(ref::FunRef, i::Integer, val) = (@assert i==1; ref[]=val)

end


text:
const EMPTY = UInt64(0)

struct ChessSet
    P::UInt64
    N::UInt64
    B::UInt64
    R::UInt64
    Q::UInt64
    K::UInt64
    friends::UInt64
end

function Base.show(io::IO, b::ChessSet)
end

struct Board
    white::ChessSet
    black::ChessSet

    taken::UInt64

    active::Bool
    castling::UInt8
    enpassant::UInt64
    halfmove::Int64
    fullmove::Int64
end

struct Piece
    type::Symbol
    square::UInt64
end
const NONE = Piece(:none, EMPTY)

struct Move
    type::Symbol
    from::UInt64
    to::UInt64
    take::Piece
    enpassant::UInt64
    promotion::Symbol
    castling::UInt8
end
const CK = UInt8(8)
const CQ = UInt8(4)
const Ck = UInt8(2)
const Cq = UInt8(1)
const NOCASTLING = UInt8(0)
const NOTAMOVE = Move(:none, EMPTY, EMPTY, NONE,
                      EMPTY, :none, NOCASTLING)

mutable struct Moves
    moves::Array{Move,1}
end

function Moves()
    return Moves(Array{Move,1}())
end

import Base.push!
function push!(list::Moves, move::Move)
    push!(list.moves, move)
end


text:
module TextureManager

export Texture

using CoreExtended
using MatrixMath
using FileManager
using JLGEngine

using ..GraphicsManager

""" TODO """
type Texture <: JLGEngine.IComponent
	id::Symbol
	enabled::Bool

	source::FileSource
	api::Any

	Texture(id=:_) = new(id,true,FileSource(),nothing)
end

API = nothing

""" TODO """
function init()
	global API = GRAPHICSDRIVER().TextureManager
end

""" TODO """
function init(this::Texture)
	this.api = API.create()
end

# preset for manager: includes code here
presetManager(Texture)

""" TODO """
load(this::Texture, path::String) = (this.source.path=path; API.load(this.api, this.source.path))

""" TODO """
bind(this::Texture) =	API.bind(api)

end #TextureManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using SafeTestsets

@time begin
    @time @safetestset "Check Position Tests" begin include("position_checking_test.jl") end
    @time @safetestset "IO Tests" begin include("io_test.jl") end
    @time @safetestset "Full Run Tests" begin include("full_run_test.jl") end
    @time @safetestset "Wall Repulsion Tests" begin include("wall_repulsion_test.jl") end
    @time @safetestset "Examples Test" begin include("examples_test.jl") end
    @time @safetestset "Exceptions Test" begin include("exceptions_test.jl") end
    @time @safetestset "Julia Logo Test" begin include("julia_logo_test.jl") end
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
# BGZF Compressor
# ===============

struct BGZFCompressor <: TranscodingStreams.Codec
    zstream::CodecZlib.ZStream
    level::Int
    windowbits::Int
end

function BGZFCompressor(;level::Integer=CodecZlib.Z_DEFAULT_COMPRESSION,
                         windowbits::Integer=CodecZlib.Z_DEFAULT_WINDOWBITS)
    return BGZFCompressor(CodecZlib.ZStream(), level, windowbits+16)
end

function Base.show(io::IO, codec::BGZFCompressor)
    print(io, summary(codec), "(level=$(codec.level),windowbits=$(codec.windowbits))")
end

const BGZFCompressorStream{S} = TranscodingStream{BGZFCompressor,S} where S<:IO

function BGZFCompressorStream(stream::IO)
    return TranscodingStream(BGZFCompressor(), stream; bufsize=BGZF_MAX_BLOCK_SIZE)
end

const EOF_BLOCK = [
    0x1f, 0x8b, 0x08, 0x04, 0x00, 0x00, 0x00, 0x00,
    0x00, 0xff, 0x06, 0x00, 0x42, 0x43, 0x02, 0x00,
    0x1b, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00,
    0x00, 0x00, 0x00, 0x00
]

function Base.close(stream::BGZFCompressorStream)
    mode = stream.state.mode
    if mode != :panic
        TranscodingStreams.changemode!(stream, :close)
    end
    if mode == :write
        @assert write(stream.stream, EOF_BLOCK) == 28
    end
    close(stream.stream)
    return nothing
end

# Methods
# -------

function TranscodingStreams.initialize(codec::BGZFCompressor)
    code = CodecZlib.deflate_init!(codec.zstream, codec.level, codec.windowbits)
    if code != CodecZlib.Z_OK
        CodecZlib.zerror(codec.zstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::BGZFCompressor)
    zstream = codec.zstream
    if zstream.state != C_NULL
        code = CodecZlib.deflate_end!(zstream)
        if code != CodecZlib.Z_OK
            CodecZlib.zerror(zstream, code)
        end
    end
    return
end

function TranscodingStreams.minoutsize(::BGZFCompressor, ::Memory)
    return Int(BGZF_MAX_BLOCK_SIZE)
end

function TranscodingStreams.process(codec::BGZFCompressor, input::Memory, output::Memory, error::Error)
    zstream = codec.zstream
    code = CodecZlib.deflate_reset!(zstream)
    if code != CodecZlib.Z_OK
        error[] = ErrorException(CodecZlib.zlib_error_message(zstream, code))
        return 0, 0, :error
    end
    avail_in = min(input.size, BGZF_SAFE_BLOCK_SIZE)
    zstream.next_in = input.ptr
    zstream.avail_in = avail_in
    avail_out = output.size - 8
    zstream.next_out = output.ptr + 8
    zstream.avail_out = avail_out
    code = CodecZlib.deflate!(zstream, CodecZlib.Z_FINISH)
    Δin = Int(avail_in - zstream.avail_in)
    Δout = Int(avail_out - zstream.avail_out) + 8
    if code == CodecZlib.Z_STREAM_END
        # BGZF header (compatible with gzip)
        output[ 1] = 0x1f  # ID1
        output[ 2] = 0x8b  # ID2
        output[ 3] = 0x08  # CM
        output[ 4] = 0x04  # FLG
        output[ 5] = 0x00  # MTIME
        output[ 6] = 0x00  # MTIME
        output[ 7] = 0x00  # MTIME
        output[ 8] = 0x00  # MTIME
        output[ 9] = 0x00  # XFL
        output[10] = 0xff  # OS
        output[11] = 0x06  # XLEN
        output[12] = 0x00  # XLEN
        output[13] = 0x42  # SI1
        output[14] = 0x43  # SI2
        unsafe_store!(Ptr{UInt16}(output.ptr+14), htol(0x0002))          # SLEN
        unsafe_store!(Ptr{UInt16}(output.ptr+16), htol(UInt16(Δout-1)))  # BSIZE
        return Δin, Δout, :end
    elseif code == Z_OK
        error[] = ErrorException("failed to deflate")
    else
        error[] = ErrorException(CodecZlib.zlib_error_message(zstream, code))
    end
    return Δin, Δout, :error
end


text:
#=
## Control System: Spaceship
As an example we want to build an IOSystem controlling the altitude of a spacecraft.
The spacecraft has mass ``m`` and can be controlled with thrusters which apply the force ``F(t)`` to the spacecraft. The altitude ``x(t)``
```math
\dot v(t) = \frac{ F(t) }{m}\\
\dot x(t) =  v(t)
```
in our model this system has the input ``F(t)``, the internal state ``v(t)`` (vertical velocity) and the output ``x(t)``.

```
       +------------+
F(t) --| spacecraft |-- x(t)
       | m, v(t)    |
       +------------+

```
=#

using BlockSystems
using ModelingToolkit
@parameters t M F(t)
@variables x(t) v(t)
D = Differential(t)

spacecraft = IOBlock([D(v) ~ F/M, D(x) ~ v], # define the equation
                     [F], # inputs of the system
                     [x], # outputs of the system
                     name = :spacecraft)

#=
We want to model a controller which takes a desired altitude as an input parameter and outputs the force for thrusters.

## Simple proportional controller
A proportional controller takes an input `i` and calculates the output proportional to the input.
```math
 o(t) = K\cdot i(t)
```
```
       +--------+
i(t) --| prop K |-- o(t)
       +--------+
```
=#
@parameters K i(t)
@variables o(t)

prop = IOBlock([o ~ K*i], [i], [o], name = :prop)
nothing # hide
#=
In order to make this useful as an controller, the input has to be the difference
between the reference and the system variable (negative feedback). We can model this
as an IOSystem where
```math
Δ = p - m\,.
```
```
       +--------+
p(t) --|  diff  |-- Δ(t)
m(t) --|        |
       +--------+
```
=#
@parameters p(t) m(t)
@variables Δ(t)
diff = IOBlock([Δ ~ p - m], [p, m], [Δ], name=:diff)
nothing # hide
#=
Now we can connect both of the defined models to create an proportional controller

```
             +----------------------------------------+
             | propc                                  |
             |         +--------+   +--------+        |
  target(t)--|--p(t) --|  diff  |---| prop K |--o(t)--|--o(t)
feedback(t)--|--m(t) --|        |   +--------+        |
             |         +--------+                     |
             +----------------------------------------+
```

If we don't provide additional information the system will try to promote all of the enclosed
variables to the new systemwide namespace.
=#
prop_c = IOSystem([diff.Δ => prop.i], # connect output of diff to input of prop
                  [diff, prop], # subsystems
                  name=:propc)

#=
For finer control, it is often preferred to give new names manually, this is
done with the `namespace_map` argument. Per default, all of the outputs of the
subsystems will become outputs of the connected system (in this case also the
output `diff.Δ`). We can prevent this by supplying the `outputs` argument
manually. Sub outputs which are not referenced here will become internal states
of the connected system.

The rhs of the namespace map can be given as a Variable/Parameter type from MTK.
For simple renaming one can also give the rhs as a `Symbol` type.
=#
prop_c = IOSystem([diff.Δ => prop.i], [diff, prop],
                  namespace_map = [prop.o => o,
                                   diff.p => :target,
                                   diff.m => :feedback],
                  outputs = [o],
                  name=:propc)

#=
Right now, the created object is a container for the two included systems. However,
it is possible to transform the object into a new `IOBlock` by calling the `connect_system`
function. The resulting is equivalent to
```
             +----------------------------------+
  target(t)--| prop_c_block                     |--o(t)
feedback(t)--| o(t)=K*(target(t) - feedback(t)) |
             +----------------------------------+
```
=#
prop_c_block = connect_system(prop_c)
nothing #hide

#=
Now we can hook our spaceship to this controller. It does not matter whether we use the
connected `IOBlock` version `prop_c` or the `IOSystem` version `prop_c_block`. We want to build
the connected system

```
           +--------------------------------------------+
           |  control system                            |
           |       +--------+   +------------+          |
target(t)--|-------| prop_c |---| spacecraft |-x(t)--+--|--altitude(t)
           |  +-fb-|        |   | m, v(t)    |       |  |
           |  |    +--------+   +------------+       |  |
           |  +--------------------------------------+  |
           +--------------------------------------------+
```
=#
@variables altitude(t)
space_controller = IOSystem([prop_c.o => spacecraft.F, spacecraft.x => prop_c.feedback],
                            [prop_c, spacecraft],
                            namespace_map = [spacecraft.x => altitude],
                            outputs = [altitude])
## we want to reduce the space_controller to a block
space_controller = connect_system(space_controller)
@info "Variables of space_controller" space_controller equations(space_controller.system)

#=
## Simulate System
In order to simulate the system we can have to build the Julia functions.
=#
gen = generate_io_function(space_controller)
nothing # hide
#=
By doing so we get access to a named tuple with the fileds
- `gen.f_ip` in-place function
- `gen.f_oop` out-of-place function
- `gen.massm` mass matrix of the system
- `gen.states` symbols of states (in order)
- `gen.inputs` symbols of inputs (in order)
- `gen.params` symbols of parameters (in order)
- (see docstring for full list)

The functions have the form
`f_ip(du, u, inputs, params, t)`
where `u` are all the states (outputs stacked on top of internal states) and `t` is the independent variable of the system.
The order of the inputs and states can be controlled.
=#
gen = generate_io_function(space_controller, f_states=[altitude, v], f_params=[K, M])
@info "Generated function" gen.massm gen.states gen.inputs gen.params
nothing # hide

#=
Well, let's see how our model is doing.
=#
using Plots
using OrdinaryDiffEq
targetfun(t) = t>1.0 ? 1.0 : 0
odefun(du, u, p, t) = gen.f_ip(du, u, [targetfun(t)], p, t)
p = [0.5, 1.0] # K, m
u0 = [0.0, 0.0] # altitude, v
tspan = (0.0, 30.0)
prob = ODEProblem(odefun, u0, tspan, p)
sol = solve(prob, Tsit5())
plot(t->sol(t)[1],tspan..., label="altitude", title="proportional control")
plot!(t->targetfun(t),tspan..., label="target")

#=
Well who could have thought, proportional control looks like an harmonic oscillator 🤷♂

## Defining a better controller
We might just add a damping term (a force proportional to the velocity of the spaceship).
If it works for a harmonic oscillator, it should work for our spaceship.
```
           +--------------------------------------------------------+
           |  control system                                        |
           |  +--------------------------------------------------+  |
           |  |    +--------+     +---+                          |  |
           |  +-v--| prop_v |-(-)-| d |     +------------+       |  |
           |       +--------+     | i |--F--| spacecraft |-v(t)--+  |
           |       +--------+     | f |     | m          |-x(t)--+--|--altitude(t)
target(t)--|-------| prop_c |-(+)-| f |     +------------+       |  |
           |  +-fb-|        |     +---+                          |  |
           |  |    +--------+                                    |  |
           |  +--------------------------------------------------+  |
           +--------------------------------------------------------+
```
In order to do so we have to slightly redefine the spaceship system: now the velocity `v(t)` is also an output and not and internal state.
=#

spacecraft = IOBlock([D(v) ~ F/M, D(x) ~ v],
                     [F],
                     [x,v],
                     name = :spacecraft)

# One can define new blocks based on previously defined blocks.
prop_v = IOBlock(prop, name=:prop_v)
fdiff = IOBlock(diff, name=:fdiff)

space_controller = IOSystem([spacecraft.v => prop_v.i,
                             spacecraft.x => prop_c.feedback,
                             prop_c.o => fdiff.p,
                             prop_v.o => fdiff.m,
                             fdiff.Δ => spacecraft.F],
                            [prop_v, prop_c, fdiff, spacecraft],
                            namespace_map = [spacecraft.x => altitude],
                            outputs = [altitude])

space_controller = connect_system(space_controller)
gen = generate_io_function(space_controller, f_states=[altitude, v], f_params=[prop_c.K, M, prop_v.K])

odefun(du, u, p, t) = gen.f_ip(du, u, [targetfun(t)], p, t)
p = [1.0, 1.0, 0.5] # propc, M, propv
u0 = [0.0, 0.0] # altitude, v
tspan = (0.0, 30.0)
prob = ODEProblem(odefun, u0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol, vars=(0,1), label="altitude", title="better control")
plot!(t->targetfun(t),tspan..., label="target")

#=
## Defining an PT1 controller
```
             +-----------------------------------------------+
             | pi_c                                   +---+  |
             |                           +------------|   |  |
             |        +----+  +--------+ | +-------+  |sum|--|--o(t)
  target(t)--|--p(t)--|diff|--| prop K |-+-| int T |--|   |  |
feedback(t)--|--m(t)--|    |  +--------+   +-------+  +---+  |
             |        +----+                                 |
             +-----------------------------------------------+
```
=#

@parameters T a1(t) a2(t)
@variables Σ(t) altitude(t)
int = IOBlock([D(o) ~ 1/T * i - o], [i], [o], name=:int)
adder = IOBlock([Σ ~ a1 + a2], [a1, a2], [Σ], name=:add)

pi_c = IOSystem([diff.Δ => prop.i,
                 prop.o => int.i,
                 prop.o => adder.a1,
                 int.o => adder.a2],
                [diff, prop, int, adder],
                namespace_map = [diff.p => :target,
                                 diff.m => :feedback,
                                 adder.Σ => o],
                outputs = [o],
                name=:pi_c)
nothing # hide

# as before we can close the loop and build the control circuit

space_controller = IOSystem([pi_c.o => spacecraft.F , spacecraft.x => pi_c.feedback],
                            [pi_c, spacecraft],
                            namespace_map = [spacecraft.x => altitude],
                            outputs = [altitude])
space_controller = connect_system(space_controller, verbose=false)
@info "Variables of space_controller" space_controller equations(space_controller.system)

# and we can simulate and plot the system
gen = generate_io_function(space_controller, f_states=[altitude], f_params=[K, T, M])

odefun(du, u, p, t) = gen.f_ip(du, u, [targetfun(t)], p, t)
p = [0.5, -1.5, 1.0] # K, T, m
u0 = [0.0, 0.0, 0.0] # altitude, int.o, v
tspan = (0.0, 50.0)
prob = ODEProblem(odefun, u0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol, vars=(0,[ 1,2 ]), label=["altitude" "integrator"], title="PT1 controller")
plot!(t->targetfun(t),tspan..., label="target")

# thank you for flying with us :)


text:
# This script is originally distributed by kdheepak under MIT License
# [kdheepak/JuliaFormatter\.vim](https://github.com/kdheepak/JuliaFormatter.vim)

using Pkg
Pkg.instantiate(verbose=false)

using JuliaFormatter
using Dates
using JSON

format_text("")

const logfile = open(joinpath(@__DIR__, "juliaformatter.log"), "w")

function log(msg; spacer = " ")
    write(logfile, "[$(Dates.now())]$spacer$msg\n")
    flush(logfile)
end

function main()

    server_state = "start"
    while server_state != "quit"
        text = readavailable(stdin)
        data = JSON.parse(String(text))
        if data["method"] == "exit"
            server_state = "quit"
        elseif data["method"] == "format"
            text = data["params"]["text"]
            output = text
            indent = typemax(Int64)
            for line in text
                if length(line) > 0
                    indent = min(length(line) - length(lstrip(line)), indent)
                end
            end
            log("Formatting: ")
            log(join(text, "\n"), spacer = "\n")
            try
                output = format_text(
                    join(text, "\n"),
                    indent = 4,
                    margin = 92,
                    always_for_in = false,
                    whitespace_typedefs = false,
                    whitespace_ops_in_indices = true,
                )
                data["status"] = "success"
            catch
                log("failed")
                output = join(text, "\n")
                data["status"] = "error"
            end
            log("\n---------------------------------------------------------------------\n")
            log(output, spacer = "\n")
            data["params"]["text"] =
                [rstrip(lpad(l, length(l) + indent)) for l in split(output, "\n")]
            print(stdout, JSON.json(data))
            log("Done.")
        end
    end

    log("exiting ...")
end

log("calling main ...")

main()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Exodus_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Exodus")
JLLWrappers.@generate_main_file("Exodus", UUID("fc7263ad-bce8-5229-97b1-c6ddabda41b1"))
end  # module Exodus_jll


text:
# Script to construct a biplot of compositional data:

using JLD2, DataFrames, CSV, Dates, Query
using LinearAlgebra
using Gadfly, Colors, Compose, ColorSchemes

include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

#load in the data:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"
#should return dataframe oligo_df_plus

## center the data:
data=convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)]) #convert columns of dataframe to just an array
y=closure(data)
cen=center(y) #calcuate center of dataset
centered=perturb(y,inverse(cen)) #perturb by inverse of this center

Z=clr(centered) #take clr transform
#sum(centered, dims=2) # should all be 1!

#Perform SVD:
bp = svd(clr(centered))
U, S, V = bp
#double checks:
#ztest = U * Diagonal(S) * V'
#Z ≈ U * Diagonal(S) * V' #this should be true!

#Construct the biplots:
#Z=FG^T
alpha=0 #α = 1, form biplot, α = 0; covariance biplot
F=U[:,1:2] * diagm(S[1:2].^alpha)
G=diagm(S[1:2].^(1-alpha)) * V'[1:2,:]
# simply scaled or unscaled singular values

#From CoDa book - doesn't say why scaling of sqrt(n) is there?
# nu=0
# a = sqrt(size(Z,1)) * U[:,1:2] * diagm(S[1:2].^(1-nu))
# b = (1/sqrt(size(Z,1))) * diagm(S[1:2].^nu) * V'[1:2,:]

#percentage of variance explained:
#values of S are the eigenvalues:
S[1].^2 / sum(S.^2) #PC 1
S[2].^2 / sum(S.^2) #PC 2
(S[2].^2 + S[1].^2) / sum(S.^2) #total var explained by first two PCs

rows=F
cols=G # / sqrt(128) #for covariance biplot, scale back

## Plot with scaled row and columns:
#  Rows are scaled by sqrt(n-1), column vectors are scaled by 1/sqrt(n-1):

Gadfly.set_default_plot_size(20cm,17cm)
s=7
coord = Coord.cartesian(xmin=-s, xmax=s, ymin=-s, ymax=s) # Geom.vector also needs scales for both axes
xsc  = Scale.x_continuous(minvalue=-s, maxvalue=s)
ysc  = Scale.y_continuous(minvalue=-s, maxvalue=s)

layer1 = layer(x=zeros(size(cols,2),1), y=zeros(size(cols,2),1), xend=(1/sqrt(128)).*cols[1,:], yend=(1/sqrt(128)).*cols[2,:], Geom.segment,style(default_color=RGB(0.4,0.4,0.4)))
layer2=layer(x=(1/sqrt(128)).*cols[1,:], y=(1/sqrt(128)).*cols[2,:],label=["O1","O2","O3","O4","O5","O6"],Geom.point,Geom.label,style(default_color=RGB(0.4,0.4,0.4)))
layer3=layer(x=sqrt(128).*rows[:,1], y=sqrt(128).*rows[:,2], color=dayofyear.(oligo_df_plus[!,:date]), Geom.point)

plt=plot(layer1,layer2,layer3,xsc, ysc,coord,
        Scale.color_continuous(minvalue=1,maxvalue=365,colormap=p -> get(ColorSchemes.jet1, p)),Guide.colorkey(title="Year Day"),
        Guide.xlabel("PCA1, 62.3%"),Guide.ylabel("PCA2, 29.5%"))

## gadfly magic:
#export:
plt |> PNG("/home/kristen/Documents/V6V8_analysis/paper_V6V8/V6V8_Syn_CoDa_paper/figures/Fig3_clrPCA.png"); #@async run(`eog /tmp/test.svg`)

##  export if needed somewhere else:
#CSV.write("/home/kristen/Documents/V6V8_analysis/scripts/rel_abnd_analysis_scripts/matlab_scripts/olgio_df_plus.csv", oligo_df_plus, delim = ',')


text:
using RecipesBase

struct T end

# This is all we define.  It uses a familiar signature, but strips it apart
# in order to add a custom definition to the internal method `RecipesBase.apply_recipe`
@recipe function plot(ap::AsymptoticPosterior, i::Int; customcolor = :green)
    # markershape --> :auto        # if markershape is unset, make it :auto
    markercolor :=  customcolor  # force markercolor to be customcolor
    # xrotation   --> 45           # if xrotation is unset, make it 45
    # zrotation   --> 90           # if zrotation is unset, make it 90
    ( x -> pdf(ap, x, i), ap.pl.map.θ[i] - 3ap.pl.map.std_estimates[i], ap.pl.map.θ[i] + 3ap.pl.map.std_estimates[i] )                   # return the arguments (input data) for the next recipe
end


text:
using Requests
using LightXML

immutable Question
    title::String
    url :: URI
    published :: DateTime
end

function Question(entry::XMLElement)
    Question(
         content(find_element(entry, "title")),
         URI(content(find_element(entry, "id"))),
         DateTime(content(find_element(entry, "published")), "yyyy-mm-ddTHH:MM:SSZ")
    )
end

"""
eg
siteurl = "http://stackoverflow.com"
tag = "julia-lang"
"""
function get_questions(siteurl, tag)
    siteuri = URI(siteurl * "/feeds/tag")
    feed_req = Requests.get(siteuri; query = Dict("tagnames" =>tag, "sort"=>"newest"))

    @assert(statuscode(feed_req)==200)
    feed_xml = parse_string(String(read(feed_req)))
    xroot = root(feed_xml)
    question_elements = get_elements_by_tagname(xroot, "entry")

    map(Question, question_elements)
end

function select_new_questions(questions, last_checked)
    questions = filter(pp -> pp.published > last_checked, questions)
    if length(questions) > 0
        last_checked =  maximum([pp.published for pp in questions]) #Use the most resent question as the one to look for questions after
    	println("updated last_checked: $(last_checked)")
    end

    questions, last_checked
end

function send_gitter_activity(webhook_url, message::String, level = "info")
    println("Sending: ", message)
    post(webhook_url; data = Dict("message" => message, "level" => level))
end

function format_question(question::Question, prefix="")
    "$(prefix): [$(question.title)]($(question.url))"
end

function send_gitter_question(webhook_url, question::Question, prefix)
    send_gitter_activity(webhook_url, format_question(question, prefix))
end

LASTCHECKED_FILENAME_SUFFIX = "_lastchecked.jsz"
function load_last_checked(prefix)
	try
		open(prefix*LASTCHECKED_FILENAME_SUFFIX, "r") do fp
			deserialize(fp)
		end
	catch ee
		warn("Failed to load Last checked. Defaulting: $ee")
		now() - Dates.Hour(24)
	end
end

function save_last_checked(prefix, last_checked)
	open(prefix*LASTCHECKED_FILENAME_SUFFIX, "w") do fp
		serialize(fp, last_checked)
	end
	println("Updated $prefix last check to $last_checked")
end

function main(prefix::String, siteurl::String, tag::String, gitter_webhook_url::URI)
	last_checked = load_last_checked(prefix)
	all_questions = get_questions(siteurl, tag)
	new_questions, new_last_checked = select_new_questions(all_questions, last_checked)

	for question in new_questions
		send_gitter_question(gitter_webhook_url, question, prefix)
	end

	save_last_checked(prefix, new_last_checked)
end

main(ARGS[1], ARGS[2], ARGS[3], URI(ARGS[4]))