text:
##=================####==============##
## TEST SIMULATION ABM NETWORK-EPIDEMICS
##=================####==============##

##=================####==============##

using StatsBase     # For sample function
using SparseArrays  # For Adyacency Matrix
using Distributions # For probability distributions

using Plots         # plotting
using LaTeXStrings  # latex labels in plots

##=================####==============##

include("NetABM.jl") # Include file with functions and definitions

##=================####==============##
# PATH TO SAVE IMGS (MODIFY WITH YOUR OWN!)

img_out_path = joinpath(homedir(), "Code", "Collaborations", "networkEpidemiology_MX", "julia_ABM", "figures")

##=================####==============##
## SYSTEM'S PARAMETERS

params = Params()

params.p_infected_t0 = 0.01 # Initial fraction of infected agents
params.attack_rate   = 0.25 # Virus attack rate
params.repetition    = 1 # Repetition ID for ensemble statistics

params.μ_cop_agents  = 0.15 # Mean of cooperation probability
params.σ_cop_agents  = 0.1 # Std deviation of cooperation probability
params.μ_recovT_agents = 7 # Mean recovery time

# Link probability (Erdos-Renyi)
params.p_link = ((1+5exp10(-2))*log(params.num_agents))/params.num_agents

##=================####==============##
# INITIALIZE ADJACENCY MATRIX (ERDOS-RENYI)

adj_mat, rows = initialize_erdos_renyi(params);

##=================####==============##
# INITIALIZE MEETS, COOPERATION AND RECOVERY TIMES DISTRIBUTION

rt_μ = params.μ_recovT_agents # MEAN RECOVERY TIME
pc_μ = params.μ_cop_agents # MEAN COOPERATION PROBABILITY
pc_σ = params.σ_cop_agents # STD_DEV COOPERATION PROBABILITY

# meets_dist = Exponential()
meets_dist  = Geometric()
coop_dist   = truncated(Normal(pc_μ, pc_σ), 0.0, 1.0)
recovt_dist = truncated(Poisson(rt_μ), 7, 35)

##=================####==============##

# DEFINE AGENTS POPULATION (DEFAULT VALUES)
agents = [Agent(id) for id in 1:params.num_agents]

# MODIFY AGENT'S CHARACTERISTICS ACCORDING TO ITS DEMOGRAPHICS
initialize_demographics!(agents, params, coop_dist, meets_dist, recovt_dist)

# ASSGINS CONTACTS TO agent FROM adj_mat
for ag in agents
    assign_contacts!(ag, agents, adj_mat, rows)
end

##=================####==============##

time_lockdown = 10

tau = []
num_S = []
num_I = []
num_R = []

# INITIALIZE TIME AND POPULATIONS
params.now_t = 0
populations = get_populations(agents, params)

while populations[2] > 0.0 # WHILE THERE ARE INFECTED AGENTS

    # RECORD DATA
    push!(tau, params.now_t)
    push!(num_S, populations[1])
    push!(num_I, populations[2])
    push!(num_R, populations[3])

    if params.now_t == time_lockdown
        set_coop_agents!(agents, params)
        println("start lockdown")
    end

   # FINDS NEXT STATE AND UPDATES ALL AGENTS
    update_all_agents!(agents, params)

    # GET CURRENT POPULATIONS AND UPDATE TIME
    global populations = get_populations(agents, params);
    params.now_t += 1
end

println("END")

##=================####==============##
# POPULATIONS GRAPH

plt = plot(
        xlabel = L"t",
        ylabel = L"N(s)/N",
        guidefontsize = 16,
        tickfontsize = 10
    )

plot!(plt, tau, num_S, color = :blue,  lab = "S", lw = 1.25)
plot!(plt, tau, num_I, color = :red,   lab = "I", lw = 1.25)
plot!(plt, tau, num_R, color = :green, lab = "R", lw = 1.25)
# plot!(plt, tau, 1.0 .- num_R, color = :green, lab = "R", lw = 1.25)

vline!(plt, [time_lockdown], ls = :dash, color = :red, lab = "")

display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "SIR_N_$(params.num_agents)_AR_$(params.attack_rate).png"))
##=================####==============##
# Distribución de Números de Encuentros

plt = plot(
        xlabel = L"N_m",
        ylabel = L"N (N_m)",
        title = "Encounters per agent per time step distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

vals = 1 .+ rand(meets_dist, params.num_agents)
histogram!(plt, vals, lab = false, bar_width = 0.85)

# vals = rand(meets_dist, params.num_agents)
# histogram!(plt, vals, normalized = true, lab = false, bar_width = 0.85)

xticks!(plt, 1:maximum(vals))
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "meetsDist_N_$(params.num_agents).png"))
##=================####==============##
# Distribución de Tiempo de Recuperación

plt = plot(
        xlabel = L"R_t",
        ylabel = L"N(R_t)",
        title = "Agent's recovery times distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

vals = rand(recovt_dist, params.num_agents)
histogram!(plt, vals, lab = false)

# histogram!(plt, vals, normalized = true, lab = false)

# xticks!(plt, unique(vals))
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "recovtDist_N_$(params.num_agents).png"))
##=================####==============##
# Tiempo de Infección

plt = plot(
        xlabel = L"t_I",
        ylabel = L"p(t_I)",
        title = "Time of infection distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

histogram!(plt, [ag.infection_t for ag in agents], normalized = true, lab = false)
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "tiempoI_N_$(params.num_agents)_AR_$(params.attack_rate).png"))
##=================####==============##
# Distribución de Grado

plt = plot(
        xlabel = L"\kappa_i",
        ylabel = L"p\;(\kappa_i)",
        title = "Network's degree distribution ",
        guidefontsize = 16,
        tickfontsize = 10
    )

histogram!(plt, [ag.degree_t for ag in agents], normalized = true, lab = false)
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "degDist_N_$(params.num_agents)_ER.png"))
##=================####==============##


text:
using NestedEnvironments
using DifferentialEquations
using Test

## user-defined nested environments
struct SubSubEnv <: AbstractEnv
end
struct SubEnv <: AbstractEnv
    env21::SubSubEnv
    env22::SubSubEnv
end
struct Env <: AbstractEnv
    env1::SubSubEnv
    env2::SubEnv
end
# dynamics
function f(x, p, t)
    ẋ1 = -x.env1
    ẋ21 = -x.env2.env21
    ẋ22 = -x.env2.env22
    (; env1 = ẋ1, env2 = (; env21 = ẋ21, env22 = ẋ22))
end
# initial condition
function initial_condition(env::Env)
    (; env1 = 1, env2 = (; env21 = rand(4), env22 = rand(2, 7)))
end

# for register
__env = Env(SubSubEnv(), SubEnv(SubSubEnv(), SubSubEnv()))
__x0 = initial_condition(__env)
@reg_env __env __x0

## test
function test1()
    env = Env(SubSubEnv(), SubEnv(SubSubEnv(), SubSubEnv()))
    x0 = initial_condition(env)
    @test x0 == @readable @raw x0
end

function test2()
    env = Env(SubSubEnv(), SubEnv(SubSubEnv(), SubSubEnv()))
    x0 = initial_condition(env)
    _x0 = @raw env x0
    x0_new = @readable env _x0
    @test x0 == x0_new
end

function test3()
    env = Env(SubSubEnv(), SubEnv(SubSubEnv(), SubSubEnv()))
    x0 = initial_condition(env)
    t0 = 0.0
    tf = 100.0
    tspan = (t0, tf)
    ts = t0:0.01:tf
    prob = ODEProblem(env, f, x0, tspan)
    sol = solve(prob, Tsit5(), saveat=ts)
    @test sol.u[1] == @raw env x0
end

function test4()
    env1 = BaseEnv()
    @test NestedEnvironments.initial_condition(env1) == 0.0
    size_tuple = (1, 2)
    env2 = BaseEnv(size_tuple)
    @test NestedEnvironments.initial_condition(env2) == zeros(size_tuple)
    env3 = BaseEnv(size_tuple...)
    @test NestedEnvironments.initial_condition(env3) == zeros(size_tuple)
    x0 = reshape(collect(1:15), 3, 5)
    env4 = BaseEnv(x0)
    @test NestedEnvironments.initial_condition(env4) == x0
end

function test()
    test1()
    test2()
    test3()
    test4()
end

ext: jl


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

import ACE
import Base: ==, convert, Dict, vec
import JuLIP.FIO: read_dict, write_dict

# ---------

"""
`Bond1ParticleFcn` : represents the environment part of a single-particle
basis function.
"""
struct Bond1ParticleFcn{TI}
   kr::TI
   kθ::TI
   kz::TI
end

_inttype(b::Bond1ParticleFcn{TI}) where {TI} = TI

_tuple(b::Bond1ParticleFcn) = (b.kr, b.kθ, b.kz)
Base.vec(b::Bond1ParticleFcn) = [b.kr, b.kθ, b.kz]

function Bond1ParticleFcn(krθz::Union{AbstractVector, Tuple}, INT=Int16)
   @assert length(krθz) == 3
   return Bond1ParticleFcn(INT(krθz[1]), INT(krθz[2]), INT(krθz[3]))
end

Bond1ParticleFcn(b::Bond1ParticleFcn, INT=Int16) =
   Bond1ParticleFcn(INT(b.kr), INT(b.kθ), INT(b.kz))

"""
`BondAList` : datastructure for storing the density projection onto the
1-particle basis functions.
"""
struct BondAList{TI}
   i2krθz::Vector{Bond1ParticleFcn{TI}}
   krθz2i::Dict{Bond1ParticleFcn{TI}, TI}
   # firstz::Vector{TI}
end

_inttype(b::BondAList{TI}) where {TI} = TI
_inttype(b::NTuple{N, TI}) where {N, TI <: Integer} = TI
_inttype(b::StaticArray{DIMS, TI}) where {DIMS, TI <: Integer} = TI

# --------------(de-)serialisation----------------------------------------
write_dict(alist::BondAList{TI}) where {TI} =
      Dict( "__id__" => "ACE_BondAList",
            "TI"     => string(TI),
            "i2krθz" => vec.(alist.i2krθz))

read_dict(::Val{:ACE_BondAList}, D::Dict) =
      BondAList( Bond1ParticleFcn.(D["i2krθz"],
                 Meta.eval(Meta.parse(D["TI"])))  )

==(al1::BondAList, al2::BondAList) = (al1.i2krθz == al2.i2krθz)
# ------------------------------------------------------------------------

Base.length(alist::BondAList) = length(alist.i2krθz)
Base.getindex(alist::BondAList, i::Integer) = alist.i2krθz[i]
Base.getindex(alist::BondAList, krθz::Bond1ParticleFcn) = alist.krθz2i[krθz]

# TODO: rename this to alloc_B ???
alloc_A(alist::BondAList, T=Float64) = zeros(Complex{T}, length(alist))

function BondAList(krθzlist::AbstractVector)
   INT = _inttype(krθzlist[1])
   i2krθz = Bond1ParticleFcn.(krθzlist, INT)
   # create the inverse mapping
   krθz2i = Dict{Bond1ParticleFcn{INT}, INT}()
   for i = 1:length(i2krθz)
      krθz2i[i2krθz[i]] = INT(i)
   end
   # TODO: z dependence
   # # find the first index for each z
   # zmax = maximum( a.z for a in i2krθz )
   # firstz = [ findfirst([a.z == iz for a in i2krθz])
   #            for iz = 1:zmax ]
   return BondAList( i2krθz, krθz2i ) # , [firstz; length(i2krθz)+1] )
end

function BondAList(krθzlist::AbstractVector{<: Bond1ParticleFcn})
   INT = _inttype(krθzlist[1])
   i2krθz = collect(krθzlist)
   # create the inverse mapping
   krθz2i = Dict{Bond1ParticleFcn{INT}, INT}()
   for i = 1:length(i2krθz)
      krθz2i[i2krθz[i]] = INT(i)
   end
   # TODO: z dependence
   # # find the first index for each z
   # zmax = maximum( a.z for a in i2krθz )
   # firstz = [ findfirst([a.z == iz for a in i2krθz])
   #            for iz = 1:zmax ]
   return BondAList( i2krθz, krθz2i ) # , [firstz; length(i2krθz)+1] )
end

# ---------

"""
`BondBasisFcnIdx{N}` : represents a single multivariate basis function
in terms of polynomial degrees in each coordinate direction.
"""
struct BondBasisFcnIdx{N, TI}
   k0::TI
   kkrθz::NTuple{N, Bond1ParticleFcn{TI}}

   function BondBasisFcnIdx(k0::TI1, kkrθz::NTuple{N, Bond1ParticleFcn{TI2}}
            ) where {TI1 <: Integer, N, TI2 <: Integer}
      return new{N, TI2}(TI2(k0), kkrθz)
   end
   function BondBasisFcnIdx(k0::TI1, kkrθz::Tuple{}
            ) where {TI1 <: Integer}
      return new{0, TI1}(TI1(k0), kkrθz)
   end
end

Base.length(b::BondBasisFcnIdx) = length(b.kkrθz)

ACE.bodyorder(b::BondBasisFcnIdx{N}) where {N} = N

BondBasisFcnIdx(k0::Integer, ainds::AbstractVector) =
      BondBasisFcnIdx(k0, tuple(ainds...))

Base.getindex(b::BondBasisFcnIdx, n::Integer) = b.kkrθz[n]

vec(b::BondBasisFcnIdx{N}) where {N} =
   vcat([b.k0], [ vec(b.kkrθz[n]) for n = 1:N ]...)

_tuple(b::BondBasisFcnIdx) = tuple( vec(b)... )
BondBasisFcnIdx(t::Tuple) =
      BondBasisFcnIdx( t[1],
                       [ Bond1ParticleFcn((t[i], t[i+1], t[i+3]))
                         for i = 2:3:length(t) ] )

"""
`BondAAList` : represents a basis functions for an EnvPairPot
"""
mutable struct BondAAList{TI}
   alist::BondAList{TI}
   i2Aidx::Matrix{TI}    # where in A can we find the ith basis function
   i2k0::Vector{TI}      # where in P0 can we find the ith basis function
   len::Vector{TI}       # body-order of ith basis function
   kkrθz2i::Dict{BondBasisFcnIdx, Int}   # inverse mapping
end

kkrθz2i(aalist::BondAAList, kkrθz::BondBasisFcnIdx) = aalist.kkrθz2i[kkrθz]

i2k0(aalist::BondAAList, i) = aalist.i2k0[i]

# dispatch this to a version with the individual arguments so we can also
# use it for de-serialization
i2kkrθz(aalist::BondAAList, i::Integer) =
      i2kkrθz(alist, aalist.i2Aidx, aalist.len, i)

i2kkrθz(alist::BondAList, i2Aidx::AbstractMatrix, len::AbstractVector,
        i::Integer) =
   [ alist.i2krθz[ i2Aidx[i, j] ] for j = 1:len[i] ]

Base.length(aalist::BondAAList) = length(aalist.len)

function BondAAList(atuples, aatuples::AbstractVector{<: NTuple{N}}) where {N}
   INT = _inttype(atuples[1])
   len = length(aatuples)
   alist = BondAList(atuples)
   # assemble the aalist elements
   i2Aidx = zeros(INT, length(aatuples), maximum(length, aatuples))
   len = zeros(INT, length(aatuples))
   kkrθz2i = Dict{BondBasisFcnIdx, Int}()
   for (i, aa) in enumerate(aatuples)
      i2Aidx[i, 1:length(aa)] .= aa
      len[i] = length(aa)
      b1 = BondBasisFcnIdx(0, ntuple(i -> alist[aa[i]], N))
      kkrθz2i[b1] = i
   end
   i2k0 =  zeros(INT, length(aatuples))
   # wrap up...
   return BondAAList(alist, i2Aidx, i2k0, len, kkrθz2i)
end

function BondAAList( Abasis::AbstractVector{<: Bond1ParticleFcn},
                    AAbasis::AbstractVector{<: BondBasisFcnIdx} )
   INT = _inttype(Abasis[1])
   alist = BondAList(Abasis)
   # assemble the aalist elements
   i2Aidx = zeros(INT, length(AAbasis), maximum(length, AAbasis))
   len = zeros(INT, length(AAbasis))
   kkrθz2i = Dict{BondBasisFcnIdx, Int}()
   i2k0 =  zeros(INT, length(AAbasis))
   for (i, AA) in enumerate(AAbasis)
      kkrθz2i[AA] = i
      len[i] = length(AA)
      i2k0[i] = AA.k0
      for j = 1:length(AA)
         i2Aidx[i, j] = alist.krθz2i[AA.kkrθz[j]]
      end
   end
   # wrap up...
   return BondAAList(alist, i2Aidx, i2k0, len, kkrθz2i)

end

# --------------(de-)serialisation----------------------------------------

write_dict(aalist::BondAAList) =
      Dict( "__id__" => "ACE_BondAAList",
            "alist"  => write_dict(aalist.alist),
            "i2Aidx" => write_dict(aalist.i2Aidx),
            "i2k0"   => aalist.i2k0,
            "len"    => aalist.len )

function read_dict(::Val{:ACE_BondAAList}, D::Dict)
   alist = read_dict(D["alist"])
   TI = _inttype(alist)
   i2Aidx = read_dict(D["i2Aidx"])::Matrix{TI}
   i2k0 = convert(Vector{TI}, D["i2k0"])
   len = convert(Vector{TI}, D["len"])
   kkrθz2i = Dict{BondBasisFcnIdx, Int}()
   for i = 1:length(len)
      kkrθz = i2kkrθz(alist, i2Aidx, len, i)
      bfcn = BondBasisFcnIdx(i2k0[i], kkrθz)
      kkrθz2i[bfcn] = Int(i)
   end
   return BondAAList(alist, i2Aidx, i2k0, len, kkrθz2i)
end

==(aal1::BondAAList, aal2::BondAAList) = (
         (aal1.i2Aidx == aal2.i2Aidx) &&
         (aal1.i2k0 == aal2.i2k0)     &&
         (aal1.alist == aal2.alist)
      )


text:
a = 0;
b = 1;
c = 0;
d = 1;
n = 9;

# set up exact solution and f
u_exact = (x,y) ->  x*(1-x)*y*(1-y)
f = (x,y)->2*y*(1-y)+2*x*(1-x)

# set up problem
problem = FDPoisson2DProblem(a,b,c,d,n,f);
# assemble system
assemble_system!(problem);
# solve
U = direct_solve_poisson2d(problem);
u = reshape(U,n,n);
# compute error
err = @. u - [u_exact(x_,y_) for x_ in problem.x, y_ in problem.y];
norm(err,Inf)<1e-14;


text:
function maximize_elbo(rng, logp, dists, ndraws, executor)
    EE = Core.Compiler.return_type(
        elbo_and_samples, Tuple{typeof(rng),typeof(logp),eltype(dists),Int}
    )
    estimates = similar(dists, EE)
    isempty(estimates) && return 0, estimates
    Folds.map!(estimates, dists, executor) do dist
        return elbo_and_samples(rng, logp, dist, ndraws)
    end
    _, iteration_opt = _findmax(estimates |> Transducers.Map(est -> est.value))
    return iteration_opt, estimates
end

function elbo_and_samples(rng, logp, dist, ndraws)
    ϕ, logqϕ = rand_and_logpdf(rng, dist, ndraws)
    logpϕ = similar(logqϕ)
    logpϕ .= logp.(eachcol(ϕ))
    logr = logpϕ - logqϕ
    elbo = Statistics.mean(logr)
    elbo_se = sqrt(Statistics.var(logr) / length(logr))
    return ELBOEstimate(elbo, elbo_se, ϕ, logpϕ, logqϕ, logr)
end

struct ELBOEstimate{T,P,L<:AbstractVector{T}}
    value::T
    std_err::T
    draws::P
    log_densities_target::L
    log_densities_fit::L
    log_density_ratios::L
end

function Base.show(io::IO, ::MIME"text/plain", elbo::ELBOEstimate)
    print(io, "ELBO estimate: ", _to_string(elbo))
    return nothing
end

function _to_string(est::ELBOEstimate; digits=2)
    return "$(round(est.value; digits)) ± $(round(est.std_err; digits))"
end


text:
export _At_mul_B

# default tolerance for matrix condition number (see 'isinvertible')
const DEFAULT_COND_TOL = 1e6

# matrix-matrix multiplication and division
@inline _At_mul_B(A, B) = transpose(A) * B
@inline _At_ldiv_B(A, B) = transpose(A) \ B

"""
    dot_zero(x::AbstractVector{N}, y::AbstractVector{N}) where{N<:Real}

Dot product with preference for zero value in the presence of infinity values.

### Input

- `x` -- first vector
- `y` -- second vector

### Output

The dot product of `x` and `y`, but with the rule that `0 * Inf == 0`.
"""
function dot_zero(x::AbstractVector{N}, y::AbstractVector{N}) where{N<:Real}
    res = zero(N)
    for i in 1:length(x)
        if !iszero(x[i]) && !iszero(y[i])
            res += x[i] * y[i]
        end
    end
    return res
end

"""
    sign_cadlag(x::N)::N where {N<:Real}

This function works like the sign function but is ``1`` for input ``0``.

### Input

- `x` -- real scalar

### Output

``1`` if ``x ≥ 0``, ``-1`` otherwise.

### Notes

This is the sign function right-continuous at zero (see
[càdlàg function](https://en.wikipedia.org/wiki/C%C3%A0dl%C3%A0g)).
It can be used with vector-valued arguments via the dot operator.

### Examples

```jldoctest
julia> LazySets.sign_cadlag.([-0.6, 1.3, 0.0])
3-element Array{Float64,1}:
 -1.0
  1.0
  1.0
```
"""
function sign_cadlag(x::N)::N where {N<:Real}
    return x < zero(x) ? -one(x) : one(x)
end

"""
    ispermutation(u::AbstractVector{T}, v::AbstractVector)::Bool where {T}

Check that two vectors contain the same elements up to reordering.

### Input

- `u` -- first vector
- `v` -- second vector

### Output

`true` iff the vectors are identical up to reordering.

### Examples

```jldoctest
julia> LazySets.ispermutation([1, 2, 2], [2, 2, 1])
true

julia> LazySets.ispermutation([1, 2, 2], [1, 1, 2])
false

```
"""
function ispermutation(u::AbstractVector{T}, v::AbstractVector)::Bool where {T}
    if length(u) != length(v)
        return false
    end
    occurrence_map = Dict{T, Int}()
    has_duplicates = false
    for e in u
        if e ∉ v
            return false
        end
        if haskey(occurrence_map, e)
            occurrence_map[e] += 1
            has_duplicates = true
        else
            occurrence_map[e] = 1
        end
    end
    if has_duplicates
        for e in v
            if !haskey(occurrence_map, e) || occurrence_map[e] == 0
                return false
            end
            occurrence_map[e] -= 1
        end
    end
    return true
end

"""
    issquare(M::AbstractMatrix)::Bool

Check whether a matrix is square.

### Input

- `M` -- matrix

### Output

`true` iff the matrix is square.
"""
function issquare(M::AbstractMatrix)::Bool
    m, n = size(M)
    return m == n
end

"""
    isinvertible(M::Matrix; [cond_tol]::Number=DEFAULT_COND_TOL)

A sufficient check of a matrix being invertible (or nonsingular).

### Input

- `M`        -- matrix
- `cond_tol` -- (optional, default: `DEFAULT_COND_TOL`) tolerance of matrix
                condition

### Output

If the result is `true`, `M` is invertible.
If the result is `false`, the matrix is non-square or this function could not
conclude.

### Algorithm

We check whether the matrix is square and whether the
[matrix condition number](https://en.wikipedia.org/wiki/Condition_number#Matrices)
`cond(M)` is below some prescribed tolerance.
"""
function isinvertible(M::Matrix; cond_tol::Number=DEFAULT_COND_TOL)
    return issquare(M) && cond(M) < cond_tol
end

# cond is not available for sparse matrices; see JuliaLang#6485 and related issues
function isinvertible(M::AbstractSparseMatrix;
                      cond_tol::Number=DEFAULT_COND_TOL)
    return issquare(M) && isinvertible(Matrix(M), cond_tol=cond_tol)
end

function isinvertible(M::Diagonal; cond_tol=nothing)
    return !any(iszero, diag(M))
end

"""
    remove_duplicates_sorted!(v::AbstractVector)

Remove duplicate entries in a sorted vector.

### Input

- `v` -- sorted vector

### Output

The input vector without duplicates.
"""
function remove_duplicates_sorted!(v::AbstractVector)
    for i in length(v)-1:-1:1
        if v[i] == v[i+1]
            splice!(v, i+1)
        end
    end
    return v
end

"""
    samedir(u::AbstractVector{N},
            v::AbstractVector{N})::Tuple{Bool, Real} where {N<:Real}

Check whether two vectors point in the same direction.

### Input

- `u` -- first vector
- `v` -- second vector

### Output

`(true, k)` iff the vectors are identical up to a positive scaling factor `k`,
and `(false, 0)` otherwise.

### Examples

```jldoctest
julia> using LazySets: samedir

julia> samedir([1, 2, 3], [2, 4, 6])
(true, 0.5)

julia> samedir([1, 2, 3], [3, 2, 1])
(false, 0)

julia> samedir([1, 2, 3], [-1, -2, -3])
(false, 0)

```
"""
function samedir(u::AbstractVector{N},
                 v::AbstractVector{N}
                )::Tuple{Bool, Real} where {N<:Real}
    @assert length(u) == length(v) "wrong dimension"
    no_factor = true
    factor = 0
    @inbounds for i in 1:length(u)
        if u[i] == 0
            if v[i] != 0
                return (false, 0)
            end
            continue
        elseif v[i] == 0
            return (false, 0)
        end
        if no_factor
            no_factor = false
            factor = u[i] / v[i]
            if factor < 0
                return (false, 0)
            end
        elseif factor != u[i] / v[i]
            return (false, 0)
        end
    end
    if no_factor
        # both vectors are zero
        return (true, 0)
    end
    return (true, factor)
end

"""
    nonzero_indices(v::AbstractVector{N})::Vector{Int} where {N<:Real}

Return the indices in which a vector is non-zero.

### Input

- `v` -- vector

### Output

A vector of ascending indices `i` such that the vector is non-zero in dimension
`i`.
"""
function nonzero_indices(v::AbstractVector{N})::Vector{Int} where {N<:Real}
    n = length(v)
    res = Vector{Int}()
    sizehint!(res, n)
    for i in 1:n
        if v[i] != zero(N)
            push!(res, i)
        end
    end
    return res
end

function nonzero_indices(v::SparseVector{N})::Vector{Int} where {N<:Real}
    return v.nzind
end

"""
    substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector

### Output

A fresh vector corresponding to `x` after `substitution` was applied.
"""
function substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(substitution, copy(x))
end

"""
    substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

### Input

- `substitution` -- substitution (a mapping from an index to a new value)
- `x`            -- vector (modified in this function)

### Output

The same (but see the Notes below) vector `x` but after `substitution` was
applied.

### Notes

The vector `x` is modified in-place if it has type `Vector` or `SparseVector`.
Otherwise, we first create a new `Vector` from it.
"""
function substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}
    return substitute!(Vector(x), substitution)
end

function substitute!(substitution::Dict{Int, T},
                     x::Union{Vector{T}, SparseVector{T}}) where {T}
    for (index, value) in substitution
        x[index] = value
    end
    return x
end

"""
    reseed(rng::AbstractRNG, seed::Union{Int, Nothing})::AbstractRNG

Reset the RNG seed if the seed argument is a number.

### Input

- `rng`  -- random number generator
- `seed` -- seed for reseeding

### Output

The input RNG if the seed is `nothing`, and a reseeded RNG otherwise.
"""
function reseed(rng::AbstractRNG, seed::Union{Int, Nothing})::AbstractRNG
    if seed != nothing
        return Random.seed!(rng, seed)
    end
    return rng
end

"""
    cross_product(M::AbstractMatrix{N}) where {N<:Real}

Compute the high-dimensional cross product of ``n-1`` ``n``-dimensional vectors.

### Input

- `M` -- ``n × n - 1``-dimensional matrix

### Output

A vector.

### Algorithm

The cross product is defined as follows:

```math
\\left[ \\dots, (-1)^{n+1} \\det(M^{[i]}), \\dots \\right]^T
```
where ``M^{[i]}`` is defined as ``M`` with the ``i``-th row removed.
See *Althoff, Stursberg, Buss: Computing Reachable Sets of Hybrid Systems Using
a Combination of Zonotopes and Polytopes. 2009.*
"""
function cross_product(M::AbstractMatrix{N})::Vector{N} where {N<:Real}
    n = size(M, 1)
    @assert 1 < n == size(M, 2) + 1 "the matrix must be n x (n-1) dimensional"

    v = Vector{N}(undef, n)
    minus = false
    for i in 1:n
        Mi = view(M, 1:n .!= i, :)  # remove i-th row
        d = det(Mi)
        if minus
            v[i] = -d
            minus = false
        else
            v[i] = d
            minus = true
        end
    end
    return v
end

"""
    StrictlyIncreasingIndices

Iterator over the vectors of `m` strictly increasing indices from 1 to `n`.

### Fields

- `n` -- size of the index domain
- `m` -- number of indices to choose (resp. length of the vectors)

### Notes

The vectors are modified in-place.

The iterator ranges over ``\\binom{n}{m}`` (`n` choose `m`) possible vectors.

This implementation results in a lexicographic order with the last index growing
first.

### Examples

```jldoctest
julia> for v in LazySets.StrictlyIncreasingIndices(4, 2)
           println(v)
       end
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]
```
"""
struct StrictlyIncreasingIndices
    n::Int
    m::Int

    function StrictlyIncreasingIndices(n::Int, m::Int)
        @assert n >= m > 0 "require n >= m > 0"
        new(n, m)
    end
end

Base.eltype(::Type{StrictlyIncreasingIndices}) = Vector{Int}
Base.length(sii::StrictlyIncreasingIndices) = binomial(sii.n, sii.m)

# initialization
function Base.iterate(sii::StrictlyIncreasingIndices)
    v = [1:sii.m;]
    return (v, v)
end

# normal iteration
function Base.iterate(sii::StrictlyIncreasingIndices, state::AbstractVector{Int})
    v = state
    i = sii.m
    diff = sii.n
    if i == diff
        return nothing
    end
    while v[i] == diff
        i -= 1
        diff -= 1
    end
    # update vector
    v[i] += 1
    for j in i+1:sii.m
        v[j] = v[j-1] + 1
    end
    # detect termination: first index has maximum value
    if i == 1 && v[1] == (sii.n - sii.m + 1)
        return (v, nothing)
    end
    return (v, v)
end

# termination
function Base.iterate(sii::StrictlyIncreasingIndices, state::Nothing)
    return nothing
end

"""
    subtypes(interface, concrete::Bool)

Return the concrete subtypes of a given interface.

### Input

- `interface` -- an abstract type, usually a set interface
- `concrete`  -- if `true`, seek further the inner abstract subtypes of the given
                 interface, otherwise return only the direct subtypes of `interface`

### Output

A list with the subtypes of the abstract type `interface`, sorted alphabetically.

### Examples

Consider the `AbstractPolytope` interface. If we include the abstract subtypes
of this interface,

```jldoctest subtypes
julia> using LazySets: subtypes

julia> subtypes(AbstractPolytope, false)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope
```

We can use this function to obtain the concrete subtypes of
`AbstractCentrallySymmetricPolytope` and `AbstractPolygon` (further until all
concrete types are obtained), using the `concrete` flag:

```jldoctest subtypes
julia> subtypes(AbstractPolytope, true)
14-element Array{Type,1}:
 Ball1
 BallInf
 HPolygon
 HPolygonOpt
 HPolytope
 Hyperrectangle
 Interval
 LineSegment
 Singleton
 SymmetricIntervalHull
 VPolygon
 VPolytope
 ZeroSet
 Zonotope
```
"""
function subtypes(interface, concrete::Bool)

    subtypes_to_test = subtypes(interface)

    # do not seek the concrete subtypes further
    if !concrete
        return sort(subtypes_to_test, by=string)
    end

    result = Vector{Type}()
    i = 0
    while i < length(subtypes_to_test)
        i += 1
        subtype = subtypes_to_test[i]
        new_subtypes = subtypes(subtype)
        if isempty(new_subtypes)
            # base type found
            push!(result, subtype)
        else
            # yet another interface layer
            append!(subtypes_to_test, new_subtypes)
        end
    end
    return sort(result, by=string)
end

"""
    delete_zero_columns(A::AbstractMatrix)

Remove all columns that only contain zeros from a given matrix.

### Input

- `A`    -- matrix
- `copy` -- (optional, default: `false`) flag to copy the matrix

### Output

A matrix.

If the input matrix `A` does not contain any zero column, we return `A` unless
the option `copy` is set.
If the input matrix contains zero columns, we always return a copy if the option
`copy` is set and otherwise a `SubArray` via `@view`.
"""
function delete_zero_columns(A::AbstractMatrix, copy::Bool=false)
    n = size(A, 2)
    nonzero_columns = Vector{Int}()
    sizehint!(nonzero_columns, n)
    for i in 1:n
        if !iszero(A[:, i])
            push!(nonzero_columns, i)
        end
    end
    if copy
        if length(nonzero_columns) == n
            return copy(A)
        else
            return A[:, nonzero_columns]
        end
    else
        if length(nonzero_columns) == n
            return A
        else
            return @view A[:, nonzero_columns]
        end
    end
end

"""
    inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
         ) where {N}

Compute the inner product ``xᵀ A y``.

### Input

- `x` -- vector on the left
- `A` -- matrix
- `y` -- vector on the right

### Output

The (scalar) result of the multiplication.
"""
function inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
              ) where {N}
    dot(x, A * y)
end

"""
    rectify(x::AbstractVector{N}) where {N<:Real}

Rectify a vector, i.e., take the element-wise maximum with zero.

### Input

- `x` -- vector

### Output

A copy of the vector where each negative entry is replaced by zero.
"""
function rectify(x::AbstractVector{N}) where {N<:Real}
    return map(xi -> max(xi, zero(N)), x)
end


text:
    export array_put
    """

        julia> array_put(aws,bucket,path,array;[level])

    Stores an array to AWS S3 bucket.

    # Signature

        function array_put{AT<:Number}(aws::AWSCore.AWSConfig,
            bucket::String,path::String,array::DenseArray{AT};
            level::Int=1,max_size::Int=2000)

    # Arguments

    - `aws`: aws config created by AWSCore.aws_config
    - `bucket`: name of AWS S3 bucket
    - `path`: file key/path name
    - `array`: dense numeric array
    - `level`: Blosc compression level (0-9); 1 is typically OK, anything above 5 is typically an over-kill
    - `max_size`: maximum array size (MB<=2000) before going into multi-part mode

    # Examples

    - `array_put(aws,"slimbucket","tmp/test/small",a)`: put array `a` into bucket `slimbucket` under path `tmp/test/small`

    # Notes:

    - use `array_delete` to delete object created with `array_put`

    """
    function array_put(aws::AWSCore.AWSConfig,bucket::String,path::String,a::DenseArray{AT};
            level::Int=1,max_size::Int=2000) where {AT<:Number}
        max_size > 2000 && warn("AWSS3/array_put: given max_size > 2000; using default 2000")
        cmp_max=min(2000*1024^2,max_size*1024^2) #blosc compression max 2147483631 bytes < (2*1024^3)
        szs=size(a)
        dims=length(szs)
        if sizeof(a)<cmp_max # single file
            tags=Dict("creator"=>"S3-SLIM","type"=>"Array")
                tags["eltype"]="$(AT)"
                tags["dims"]="$(dims)"
                tags["ns"]=join(map(i->(@sprintf "%d" szs[i]),1:dims),":")
            ac=Blosc.compress(vec(a);level=level);
            s3_put(aws, bucket, path, ac, tags=tags);
            return nothing
        else # multi-part files
            #warn("AWSS3/array_put: large array - going into multi-part mode";key="AWS S3 array_put",once=true)
            (nfiles,nelmts,parts,idxs,idxe)=file_parts(a,cmp_max)
                tags=Dict("creator"=>"S3-SLIM","type"=>"metaArray")
                tags["nfiles"]="$(nfiles)"
                tags["nelmts"]="$(nelmts)"
                tags["eltype"]="$(AT)"
                tags["dims"]="$(dims)"
                tags["ns"]=join(map(i->(@sprintf "%d" szs[i]),1:dims),":")
            #println((nfiles,nelmts,parts,idxs,idxe))
            #for i=1:nfiles println((i,parts[i],idxs[i],idxe[i])) end
            av=vec(a)
            for i=1:nfiles
                ppath=@sprintf("%s-parts/%6.6d",path,i)
                #println(ppath)
                pc=Blosc.compress(av[idxs[i]:idxe[i]];level=level);
                s3_put(aws, bucket, ppath, pc);
            end
            d=[parts idxs idxe]
            dc=Blosc.compress(vec(d);level=level);
            s3_put(aws, bucket, path, dc, tags=tags);
            return nothing
        end
    end

    export array_get
    """

        julia> array_get(aws,bucket,path;delete=false)

    Reads from AWS S3 bucket an array stored by `array_put`.

    # Signature

        function array_get(aws::AWSCore.AWSConfig,
            bucket::String,path::String;delete::Bool=false)

    # Arguments

    - `aws`: aws config created by AWSCore.aws_config
    - `bucket`: name of AWS S3 bucket
    - `path`: file key/path name
    - `delete`: delete file key/path after reading

    # Examples

    - `A=array_get(aws,"slimbucket","tmp/test/small")`: gets array `a` from bucket `slimbucket` and path `tmp/test/small`

    # Notes:

    - use `array_delete` to delete object created with `array_put`

    """
    function array_get(aws::AWSCore.AWSConfig,bucket::String,path::String;delete::Bool=false)
        s3_exists(aws, bucket, path) || error("AWSS3/array_get: file $path does not exist in $bucket.")
        tags=s3_get_tags(aws, bucket, path);
        (haskey(tags,"creator")&&tags["creator"]=="S3-SLIM") || error("AWSS3/array_get: file $path in $bucket is unknown.")
        haskey(tags,"type") || error("AWSS3/array_get: file $path in $bucket does not have known type.")
        if tags["type"]=="Array" # single file
            eval(Meta.parse("edt=$(tags["eltype"])"))
            eval(Meta.parse("dims=$(tags["dims"])"))
            szs=parse.(Int,split(tags["ns"],":"))
            ac=s3_get(aws, bucket, path);
            a=reshape(Blosc.decompress(edt,ac),(szs...,));
            delete && array_delete(aws, bucket, path);
            return a
        elseif tags["type"]=="metaArray" # multi-part files
            eval(Meta.parse("nfiles=$(tags["nfiles"])"))
            eval(Meta.parse("nelmts=$(tags["nelmts"])"))
            eval(Meta.parse("edt=$(tags["eltype"])"))
            eval(Meta.parse("dims=$(tags["dims"])"))
            szs=parse.(Int,split(tags["ns"],":"))
            dc=s3_get(aws, bucket, path);
            d=reshape(Blosc.decompress(Int,dc),(nfiles,3));
            parts=d[:,1]; idxs=d[:,2]; idxe=d[:,3];
            #println((nfiles,nelmts,parts,idxs,idxe))
            #for i=1:nfiles println((i,parts[i],idxs[i],idxe[i])) end
            av=Vector{edt}(undef,nelmts)
            for i=1:nfiles
                ppath=@sprintf("%s-parts/%6.6d",path,i)
                #println(ppath)
                pc=s3_get(aws, bucket, ppath);
                av[idxs[i]:idxe[i]]=Blosc.decompress(edt,pc)
            end
            a=reshape(av,(szs...,))
            delete && array_delete(aws, bucket, path);
            return a
        else
            error("AWSS3/array_get: file $path in $bucket is not an array stored with array_put.")
        end
    end

    export array_delete
    """

        julia> array_delete(aws,bucket,path)

    Deletes from AWS S3 bucket an array stored by `array_put`.

    # Signature

        function array_delete(aws::AWSCore.AWSConfig,
            bucket::String,path::String)

    # Arguments

    - `aws`: aws config created by AWSCore.aws_config
    - `bucket`: name of AWS S3 bucket
    - `path`: file key/path name

    # Examples

    - `array_delete(aws,"slimbucket","tmp/test/small")`: deletes array stored by `array_put` as `tmp/test/small` and bucket `slimbucket`

    """
    function array_delete(aws::AWSCore.AWSConfig,bucket::String,path::String)
        s3_exists(aws, bucket, path) || error("AWSS3/array_delete: file $path does not exist in $bucket.")
        tags=s3_get_tags(aws, bucket, path);
        (haskey(tags,"creator")&&tags["creator"]=="S3-SLIM") || error("AWSS3/array_delete: file $path in $bucket is unknown.")
        haskey(tags,"type") || error("AWSS3/array_delete: file $path in $bucket does not have known type.")
        if tags["type"]=="Array" # single file
            s3_delete(aws,bucket,path)
        elseif tags["type"]=="metaArray" # multi-part files
            eval(Meta.parse("nfiles=$(tags["nfiles"])"))
            for i=1:nfiles
                ppath=@sprintf("%s-parts/%6.6d",path,i)
                s3_delete(aws,bucket,ppath)
            end
            s3_delete(aws,bucket,path)
        else
            error("AWSS3/array_delete: file $path in $bucket is not an array.")
        end
        return nothing
    end


text:
function process_assets(request::HTTP.Request, state::HandlerState; file_path::AbstractString)
    app = state.app
    filename = joinpath(get_assets_path(app), file_path)
    try
        headers = Pair{String,String}[]
        mimetype = mime_by_path(filename)
        !isnothing(mimetype) && push!(headers, "Content-Type" => mimetype)
        file_contents = read(filename)
        return HTTP.Response(200, headers;body = file_contents)
    catch
        return HTTP.Response(404)
    end

end


text:
export DefaultRecipe, qsfit, qsfit_multi

abstract type DefaultRecipe <: AbstractRecipe end

function default_options(::Type{T}) where T <: DefaultRecipe
    out = OrderedDict{Symbol, Any}()
    out[:wavelength_range] = [1215, 7.3e3]
    out[:min_spectral_coverage] = Dict(:default => 0.6,
                                       :ironuv  => 0.3,
                                       :ironopt => 0.3)
    out[:skip_lines] = Symbol[]

    out[:host_template] = Dict(:library=>"swire", :template=>"Ell5")
    out[:use_host_template] = true
    out[:host_template_range] = [4000., 7000.]

    out[:use_balmer] = true
    out[:use_ironuv] = true;      out[:ironuv_fwhm]    = 3000.
    out[:use_ironopt] = true;     out[:ironoptbr_fwhm] = 3000.;  out[:ironoptna_fwhm] =  500.

    out[:line_profiles] = :gauss
    out[:line_broadening] = true
    out[:iron_broadening] = true

    out[:n_unk] = 10
    out[:unk_avoid] = [4863 .+ [-1,1] .* 50, 6565 .+ [-1,1] .* 150]

    return out
end

function known_spectral_lines(source::QSO{T}) where T <: DefaultRecipe
    list = [
        MultiCompLine(                      :Lyb                     , [BroadLine, NarrowLine]),
        # NarrowLine( custom_transition(tid=:OV_1213     , 1213.8  )),  # Ferland+92, Shields+95
        MultiCompLine(                      :Lya                     , [BroadLine, NarrowLine]),
        # NarrowLine( custom_transition(tid=:OV_1218     , 1218.3  )),  # Ferland+92, Shields+95
        NarrowLine(                         :NV_1241                ),
        BroadLine(                          :OI_1306                ),
        BroadLine(                          :CII_1335               ),
        BroadLine(                          :SiIV_1400              ),
        MultiCompLine(                      :CIV_1549                , [BroadLine, NarrowLine]),
        BroadLine(                          :HeII_1640              ),
        BroadLine(                          :OIII_1664              ),
        BroadLine(                          :AlIII_1858             ),
        BroadLine(                          :CIII_1909              ),
        BroadLine(                          :CII_2326               ),
        BroadLine(    custom_transition(tid=:F2420       , 2420.0  )),
        MultiCompLine(                      :MgII_2798               , [BroadLine, NarrowLine]),
        NarrowLine(                         :NeV_3345               ),
        NarrowLine(                         :NeV_3426               ),
        NarrowLine(                         :OII_3727               ),
        NarrowLine(                         :NeIII_3869             ),
        BroadLine(                          :Hd                     ),
        BroadLine(                          :Hg                     ),
        NarrowLine(                         :OIII_4363              ),
        BroadLine(                          :HeII_4686              ),
        MultiCompLine(                      :Hb                      , [BroadLine, NarrowLine]),
        NarrowLine(                         :OIII_4959              ),
        NarrowLine(                         :OIII_5007              ),
        NarrowLine(   cname=:OIII_5007_bw,  :OIII_5007              ),
        BroadLine(                          :HeI_5876               ),
        NarrowLine(                         :OI_6300                ),
        NarrowLine(                         :OI_6364                ),
        NarrowLine(                         :NII_6549               ),
        MultiCompLine(                      :Ha                      , [BroadLine, NarrowLine, BroadBaseLine]),
        NarrowLine(                         :NII_6583               ),
        NarrowLine(                         :SII_6716               ),
        NarrowLine(                         :SII_6731               )]
    return list
end

function LineComponent(source::QSO{T}, line::GenericLine, multicomp::Bool) where T <: DefaultRecipe
    tt = transition(line.tid)
    λ = tt.LAMBDA_VAC_ANG
    if length(λ) > 1
        println(logio(source), "Considering average wavelength for the $(line.tid) multiplet: " * string(mean(λ)) * "Å")
        λ = mean(λ)  # average lambda o multiplets
    else
        λ = λ[1]
    end
    @assert source.options[:line_profiles] in [:gauss, :lorentz, :voigt]
    if source.options[:line_profiles] == :gauss
        comp = SpecLineGauss(λ)
    elseif source.options[:line_profiles] == :lorentz
        comp = SpecLineLorentz(λ)
    elseif source.options[:line_profiles] == :voigt
        comp = SpecLineVoigt(λ)
    else
        error("options[:line_profiles] must be one of: :gauss, :lorentz, :voigt.")
    end
    return LineComponent(line, comp, multicomp)
end

function LineComponent(source::QSO{T}, line::BroadLine, multicomp::Bool) where T <: DefaultRecipe
    comp = LineComponent(source, GenericLine(line.tid), multicomp).comp
    comp.fwhm.val  = 5e3
    comp.fwhm.low  = 900
    comp.fwhm.high = 1.5e4
    comp.voff.low  = -3e3
    comp.voff.high =  3e3

    if line.tid == :MgII_2798
        comp.voff.low  = -1e3
        comp.voff.high =  1e3
    end
    return LineComponent(line, comp, multicomp)
end

function LineComponent(source::QSO{T}, line::NarrowLine, multicomp::Bool) where T <: DefaultRecipe
    comp = LineComponent(source, GenericLine(line.tid), multicomp).comp
    comp.fwhm.val  = 5e2
    comp.fwhm.low  = 100
    comp.fwhm.high = (multicomp  ?  1e3  :  2e3)
    comp.voff.low  = -1e3
    comp.voff.high =  1e3

    if line.cname == :OIII_5007_bw
        comp.fwhm.val  = 500
        comp.fwhm.high = 1e3
        comp.voff.low  = 0
        comp.voff.high = 2e3
    end
    return LineComponent(line, comp, multicomp)
end

function LineComponent(source::QSO{T}, line::BroadBaseLine, multicomp::Bool) where T <: DefaultRecipe
    comp = LineComponent(source, GenericLine(line.tid), multicomp).comp
    comp.fwhm.val  = 2e4
    comp.fwhm.low  = 1e4
    comp.fwhm.high = 3e4
    comp.voff.fixed = true
    return LineComponent(line, comp, multicomp)
end

function PreparedSpectrum(source::QSO{T}; id=1) where T <: DefaultRecipe
    data = deepcopy(source.specs[id])
    println(logio(source), "Spectrum: " * data.label)
    println(logio(source), "  good fraction:: ", goodfraction(data))
    if goodfraction(data) < 0.5
        error("Good fraction < 0.5")
    end
    println(logio(source), "  resolution: ", @sprintf("%.4g", data.resolution), " km / s (FWHM)")

    λ = data.λ ./ (1 + source.z)
    data.good[findall(λ .< source.options[:wavelength_range][1])] .= false
    data.good[findall(λ .> source.options[:wavelength_range][2])] .= false

    #= Emission line are localized features whose parameter can be
    reliably estimated only if there are sufficient samples to
    constrain the corresponding parameters.  If data coverage is
    not sufficient the component should not be added to the model,
    and corresponding spectral samples should be ignored to avoid
    worsening the fit due to missing model components. =#
    println(logio(source), "Good samples before line coverage filter: ", length(findall(data.good)))

    # Collect LineComponent objects
    lcs = collect_LineComponent(source)

    # Sort lnames according to the expected line FWHM.  This is
    # necessary to avoid cases where a single line has sufficient
    # coverage, but the check for a subsequent broad lines has not
    # (dropping the "good" flag also for the narrow one).
    lnames = collect(keys(lcs))[reverse(sortperm(getfield.(getfield.(getfield.(values(lcs), :comp), :fwhm), :val)))]
    for lname in lnames
        lc = lcs[lname]
        (λmin, λmax, coverage) = spectral_coverage(λ .* data.good, data.resolution, lc.comp)
        coverage = round(coverage * 1e3) / 1e3  # keep just 3 significant digits...
        threshold = get(source.options[:min_spectral_coverage], lname, source.options[:min_spectral_coverage][:default])
        print(logio(source), @sprintf("Line %-15s coverage: %4.2f (threshold: %4.2f)", lname, coverage, threshold))
        if coverage < threshold
            print(logio(source), @sprintf("  neglecting range: %10.5g < λ < %10.5g", λmin, λmax))
            ii = findall(λmin .<= λ .< λmax)
            data.good[ii] .= false
            delete!(lcs, lname)
        end
        println(logio(source))
    end
    println(logio(source), "Good samples after line coverage filtering: ", length(findall(data.good)))

    # De-reddening
    dered = ccm_unred([1450, 3000, 5100.], source.mw_ebv)
    println(logio(source), "Dereddening factors @ 1450, 3000, 5100 AA: ", dered)
    dered = ccm_unred(data.λ, source.mw_ebv)

    ii = findall(data.good)
    dom = Domain(data.λ[ii] ./ (1 + source.z))
    lum = Measures(dom,
                   data.flux[ii] .* dered[ii] .* source.flux2lum .* (1 + source.z),
                   data.err[ ii] .* dered[ii] .* source.flux2lum .* (1 + source.z))

    return PreparedSpectrum(id, data, dom, lum, lcs)
end

function minimizer(source::QSO{T}) where T <: DefaultRecipe
    mzer = GFit.cmpfit()
    mzer.Δfitstat_threshold = 1.e-5
    return mzer
end

function fit!(source::QSO{T}, model::Model, pspec::PreparedSpectrum) where T <: DefaultRecipe
    mzer = minimizer(source)
    fitres = fit!(model, pspec.data, minimizer=mzer)
    show(logio(source), model)
    show(logio(source), fitres)
    # @gp :QSFit pspec.data model
    # printstyled(color=:blink, "Press ENTER to continue..."); readline()
    return fitres
end

function fit!(source::QSO{T}, multi::MultiModel, pspecs::Vector{PreparedSpectrum}) where T <: DefaultRecipe
    mzer = minimizer(source)
    fitres = fit!(multi, getfield.(pspecs, :data), minimizer=mzer)
    show(logio(source), multi)
    show(logio(source), fitres)
    # for id in 1:length(multi)
    #     @gp Symbol("QSFit$(id)") pspecs[id].data multi[id]
    # end
    # printstyled(color=:blink, "Press ENTER to continue..."); readline()
    return fitres
end

function add_qso_continuum!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    λ = domain(model)[:]

    comp = QSFit.powerlaw(3000)
    comp.x0.val = median(λ)
    comp.norm.val = Spline1D(λ, pspec.data.val, k=1, bc="error")(comp.x0.val)
    comp.norm.low = comp.norm.val / 1000.  # ensure contiuum remains positive (needed to estimate EWs)
    comp.alpha.val  = -1.5
    comp.alpha.low  = -3
    comp.alpha.high =  1

    model[:qso_cont] = comp
    push!(model[:Continuum].list, :qso_cont)
    evaluate!(model)
end

function add_host_galaxy!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_host_template]  &&
        (source.options[:host_template_range][1] .< maximum(λ))  &&
        (source.options[:host_template_range][2] .> minimum(λ))
        model[:galaxy] = QSFit.hostgalaxy(source.options[:host_template])
        push!(model[:Continuum].list, :galaxy)

        # Split total flux between continuum and host galaxy
        vv = Spline1D(λ, pspec.data.val, k=1, bc="extrapolate")(5500.)
        model[:galaxy].norm.val    = 1/2 * vv
        model[:qso_cont].norm.val *= 1/2 * vv / Spline1D(λ, model(:qso_cont), k=1, bc="extrapolate")(5500.)
        evaluate!(model)
    end
end

function add_balmer_cont!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    if source.options[:use_balmer]
        model[:balmer] = QSFit.balmercont(0.1, 0.5)
        push!(model[:Continuum].list, :balmer)
        c = model[:balmer]
        c.norm.val  = 0.1
        c.norm.fixed = false
        c.norm.high = 0.5
        c.ratio.val = 0.5
        c.ratio.fixed = false
        c.ratio.low  = 0.1
        c.ratio.high = 1
        @try_patch! model[:balmer].norm *= model[:qso_cont].norm
        evaluate!(model)
    end
end

function renorm_cont!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    freeze(model, :qso_cont)
    c = model[:qso_cont]
    initialnorm = c.norm.val
    if c.norm.val > 0
        println(logio(source), "Cont. norm. (before): ", c.norm.val)
        while true
            residuals = (model() - pspec.data.val) ./ pspec.data.unc
            ratio = count(residuals .< 0) / length(residuals)
            (ratio > 0.9)  &&  break
            (c.norm.val < (initialnorm / 5))  &&  break # give up
            c.norm.val *= 0.99
            evaluate!(model)
        end
        println(logio(source), "Cont. norm. (after) : ", c.norm.val)
    else
        println(logio(source), "Skipping cont. renormalization")
    end
    evaluate!(model)
    # @gp (domain(model), pspec.data) model
    # printstyled(color=:blink, "Press ENTER to continue..."); readline()
end

function guess_norm_factor!(pspec::PreparedSpectrum, model::Model, name::Symbol; quantile=0.95)
    @assert model[name].norm.val != 0
    m = model(name)
    c = cumsum(m)
    @assert maximum(c) != 0. "Model for $name evaluates to zero over the whole domain"
    c ./= maximum(c)
    i1 = findfirst(c .> ((1 - quantile)/2))
    i2 = findlast( c .< ((1 + quantile)/2))
    if i1 >= i2
        return #Can't calculate normalization for component
    end
    resid = pspec.data.val - model()
    ratio = model[name].norm.val / sum(m[i1:i2])
    off = sum(resid[i1:i2]) * ratio
    model[name].norm.val += off
    @assert !isnan(off) "Norm. offset is NaN for $name"
    if model[name].norm.val < 0  # ensure line has positive normalization
        model[name].norm.val = abs(off)
    end
end

function add_iron_uv!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_ironuv]
        fwhm = source.options[:ironuv_fwhm]
        if source.options[:iron_broadening]
            fwhm = sqrt(fwhm^2 + pspec.orig.resolution^2)
        end
        comp = QSFit.ironuv(fwhm)
        (_1, _2, coverage) = spectral_coverage(λ, pspec.orig.resolution, comp)
        threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])
        if coverage >= threshold
            model[:ironuv] = comp
            model[:ironuv].norm.val = 1.
            push!(model[:Iron].list, :ironuv)
            evaluate!(model)
            QSFit.guess_norm_factor!(pspec, model, :ironuv)
            evaluate!(model)
        else
            println(logio(source), "Ignoring ironuv component (threshold: $threshold)")
        end
    end
end

function add_iron_opt!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_ironopt]
        fwhm = source.options[:ironoptbr_fwhm]
        if source.options[:iron_broadening]
            fwhm = sqrt(fwhm^2 + pspec.orig.resolution^2)
        end
        comp = QSFit.ironopt_broad(fwhm)
        (_1, _2, coverage) = spectral_coverage(λ, pspec.orig.resolution, comp)
        threshold = get(source.options[:min_spectral_coverage], :ironopt, source.options[:min_spectral_coverage][:default])
        if coverage >= threshold
            model[:ironoptbr] = comp
            model[:ironoptbr].norm.val = 1 # TODO: guess a sensible value
            fwhm = source.options[:ironoptna_fwhm]
            if source.options[:iron_broadening]
                fwhm = sqrt(fwhm^2 + pspec.orig.resolution^2)
            end
            model[:ironoptna] = QSFit.ironopt_narrow(fwhm)
            model[:ironoptna].norm.val = 1 # TODO: guess a sensible value
            model[:ironoptna].norm.fixed = false
            push!(model[:Iron].list, :ironoptbr)
            push!(model[:Iron].list, :ironoptna)
            evaluate!(model)
            QSFit.guess_norm_factor!(pspec, model, :ironoptbr)
            evaluate!(model)
        else
            println(logio(source), "Ignoring ironopt component (threshold: $threshold)")
        end
    end
end

function add_emission_lines!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    groups = OrderedDict{Symbol, Vector{Symbol}}()
    for (cname, lc) in pspec.lcs

        # All QSFit line progiles take spectral resolution into
        # account.  This is significantly faster than convolving the
        # whole model with an instrument response but has some
        # limitations:
        # - works only with built-in profiles (Gaussian, Lorentzian, Voigt);
        # - all components must be additive (i.e. no absorptions);
        # - further narrow components (besides known emission lines)
        #   will not be corrected for instrumental resolution.
        if source.options[:line_broadening]
            lc.comp.resolution = pspec.orig.resolution
        end

        model[cname] = lc.comp
        grp = group(lc.line)
        haskey(groups, grp)  ||  (groups[grp] = Vector{Symbol}())
        push!(groups[grp], cname)
    end
    for (group, lnames) in groups
        model[group] = SumReducer(lnames)
    end
    evaluate!(model)
end

function guess_emission_lines!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    groups_to_go = unique([QSFit.group(v.line) for (k,v) in pspec.lcs])
    for grp in [:BroadLines, :NarrowLines, :BroadBaseLines, :AsymmTailLines]  # Note: order is important
        found = false
        for (cname, lc) in pspec.lcs
            (grp == group(lc.line))  ||  continue
            QSFit.guess_norm_factor!(pspec, model, cname)
            found = true
        end
        if found
            push!(model[:main].list, grp)
            deleteat!(groups_to_go, findfirst(groups_to_go .== grp))
        end
        evaluate!(model)
    end

    # Ensure all groups have been considered
    @assert length(groups_to_go) == 0 "The following line groups have not been considered: $groups_to_go"
end

function add_patch_functs!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    # Patch parameters
    @try_patch! begin
        # model[:OIII_4959].norm = model[:OIII_5007].norm / 3
        model[:OIII_4959].voff = model[:OIII_5007].voff
    end
    @try_patch! begin
        model[:OIII_5007_bw].voff += model[:OIII_5007].voff
        model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
    end
    @try_patch! begin
        # model[:OI_6300].norm = model[:OI_6364].norm / 3
        model[:OI_6300].voff = model[:OI_6364].voff
    end
    @try_patch! begin
        # model[:NII_6549].norm = model[:NII_6583].norm / 3
        model[:NII_6549].voff = model[:NII_6583].voff
    end
    @try_patch! begin
        # model[:SII_6716].norm = model[:SII_6731].norm / 1.5
        model[:SII_6716].voff = model[:SII_6731].voff
    end

    @try_patch! model[:Hb_na].voff = model[:Ha_na].voff

    # The following are required to avoid degeneracy with iron
    # template
    @try_patch! begin
        model[:Hg].voff = model[:Hb_br].voff
        model[:Hg].fwhm = model[:Hb_br].fwhm
    end
    @try_patch! begin
        model[:Hg_br].voff = model[:Hb_br].voff
        model[:Hg_br].fwhm = model[:Hb_br].fwhm
    end
    @try_patch! begin
        model[:Hg_na].voff = model[:Hb_na].voff
        model[:Hg_na].fwhm = model[:Hb_na].fwhm
    end

    # Ensure luminosity at peak of the broad base component is
    # smaller than the associated broad component:
    if  haskey(model, :Hb_br)  &&
        haskey(model, :Hb_bb)
        model[:Hb_bb].norm.high = 1
        model[:Hb_bb].norm.val  = 0.5
        @try_patch! model[:Hb_bb].norm *= model[:Hb_br].norm / model[:Hb_br].fwhm * model[:Hb_bb].fwhm
    end
    if  haskey(model, :Ha_br)  &&
        haskey(model, :Ha_bb)
        model[:Ha_bb].norm.high = 1
        model[:Ha_bb].norm.val  = 0.5
        @try_patch! model[:Ha_bb].norm *= model[:Ha_br].norm / model[:Ha_br].fwhm * model[:Ha_bb].fwhm
    end
end

function default_unk_line(source::QSO{T}) where T <: DefaultRecipe
    comp = SpecLineGauss(5e3)
    comp.norm.val = 0.
    comp.center.fixed = false
    comp.center.low = 0
    comp.center.high = Inf
    comp.fwhm.val  = 5e3
    comp.fwhm.low  = 600
    comp.fwhm.high = 1e4
    comp.voff.fixed = true
    return comp
end

function add_unknown_lines!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    (source.options[:n_unk] > 0)  ||  (return nothing)
    λ = domain(model)[:]
    for i in 1:source.options[:n_unk]
        model[Symbol(:unk, i)] = default_unk_line(source)
    end
    model[:UnkLines] = SumReducer([Symbol(:unk, i) for i in 1:source.options[:n_unk]])
    push!(model[:main].list, :UnkLines)
    evaluate!(model)
    for j in 1:source.options[:n_unk]
        freeze(model, Symbol(:unk, j))
    end
    evaluate!(model)

    # Set "unknown" line center wavelength where there is a maximum in
    # the fit residuals, and re-run a fit.
    λunk = Vector{Float64}()
    while true
        (length(λunk) >= source.options[:n_unk])  &&  break
        evaluate!(model)
        Δ = (pspec.data.val - model()) ./ pspec.data.unc

        # Avoid considering again the same region (within 1A) TODO: within resolution
        for l in λunk
            Δ[findall(abs.(l .- λ) .< 1)] .= 0.
        end

        # Avoidance regions
        for rr in source.options[:unk_avoid]
            Δ[findall(rr[1] .< λ .< rr[2])] .= 0.
        end

        # Do not add lines close to from the edges since these may
        # affect qso_cont fitting
        Δ[findall((λ .< minimum(λ)*1.02)  .|
                  (λ .> maximum(λ)*0.98))] .= 0.
        iadd = argmax(Δ)
        (Δ[iadd] <= 0)  &&  break  # No residual is greater than 0, skip further residuals....
        push!(λunk, λ[iadd])

        cname = Symbol(:unk, length(λunk))
        model[cname].norm.val = 1.
        model[cname].center.val  = λ[iadd]
        model[cname].center.low  = λ[iadd] - λ[iadd]/10. # allow to shift 10%
        model[cname].center.high = λ[iadd] + λ[iadd]/10.

        thaw(model, cname)
        fitres = fit!(source, model, pspec)
        freeze(model, cname)
    end
    evaluate!(model)
end

function neglect_weak_features!(source::QSO{T}, pspec::PreparedSpectrum, model::Model, fitres::GFit.FitResult) where T <: DefaultRecipe
    # Disable "unknown" lines whose normalization uncertainty is larger
    # than X times the normalization
    needs_fitting = false
    for ii in 1:source.options[:n_unk]
        cname = Symbol(:unk, ii)
        isfixed(model, cname)  &&  continue
        if model[cname].norm.val == 0.
            freeze(model, cname)
            needs_fitting = true
            println(logio(source), "Disabling $cname (norm. = 0)")
        elseif model[cname].norm.unc / model[cname].norm.val > 3
            model[cname].norm.val = 0.
            freeze(model, cname)
            needs_fitting = true
            println(logio(source), "Disabling $cname (unc. / norm. > 3)")
        end
    end
    return needs_fitting
end

include("DefaultRecipe_single.jl")
include("DefaultRecipe_multi.jl")


text:
@kwdef struct CarX{R} <: FieldVector{5,R}
    "x position"
    x::R
    "y position"
    y::R
    "orientation"
    θ::R
    "linear speed"
    v::R = 0.0
    "steering angle"
    ψ::R = 0.0
end

get_pos(s::CarX) = @SVector[s.x, s.y]

@kwdef struct CarU{R} <: FieldVector{2,R}
    "desired linear speed"
    v̂::R = 0.0
    "desired steering angle"
    ψ̂::R = 0.0
end

# TODO: fine-tune these parameters
@kwdef struct CarDynamics{NF} <: SysDynamics
    "control time interval in seconds"
    delta_t::ℝ
    "maximal linear speed"
    max_v::ℝ = 5.0
    "maximal steering angle"
    max_ψ::ℝ = 30°
    "wheelbase between the front and rear wheels"
    l::ℝ = 0.1
    "rate of convergence for v to converge to v̂"
    k_v::ℝ = 5.0
    "rate of convergence for ψ to converge to ψ̂"
    k_ψ::ℝ = 10.0
    "add_noise(x, t) -> x′"
    add_noise::NF = (x, t) -> x
    integrator_samples::ℕ = 1
end

ψ_from_v_ω(v, ω, l) = (abs(v) < 1e-4) ? 0.0 : atan(ω * l / v)
ω_from_v_ψ(v, ψ, l) = tan(ψ) * v / l

function u_from_v_ω(v, ω, dy::CarDynamics)
    v = clamp(v, -dy.max_v, dy.max_v)
    ψ = ψ_from_v_ω(v, ω, dy.l)
    ψ = clamp(ψ, -dy.max_ψ, dy.max_ψ)
    CarU(v, ψ)
end

# Simple clamping limit
function OneVision.limit_control(dy::CarDynamics, u::U, x, t)::U where {U}
    v̂, ψ̂ = u
    v̂1 = clamp(v̂, -dy.max_v, dy.max_v)
    ψ̂1 = clamp(ψ̂, -dy.max_ψ, dy.max_ψ)
    if v̂1 == v̂ && ψ̂1 == ψ̂
        u
    else
        U(v̂, ψ̂)
    end
end

@inline function sys_derivates(dy::CarDynamics, x::X, u)::X where X
    x, y, θ, v, ψ = x
    v̂, ψ̂ = u

    ẋ = cos(θ) * v
    ẏ = sin(θ) * v
    θ̇ = ω_from_v_ψ(v, ψ, dy.l)
    v̇ = dy.k_v * (v̂ - v)
    ψ̇ = dy.k_ψ * (ψ̂ - ψ)

    X(ẋ, ẏ, θ̇, v̇, ψ̇)
end

function OneVision.sys_forward(dy::CarDynamics, x::X, u, t::𝕋)::X where X
    N = dy.integrator_samples
    dt = dy.delta_t / N
    @inline f(x) = sys_derivates(dy, x, u)

    x′ = integrate_forward_invariant(f, x, dt, RK38, N)
    dy.add_noise(x′,t)
end

"""
Should implement `track_ref`.
"""
abstract type TrackingControl end

@kwdef struct RefPointTrackControl <: TrackingControl
    dy::CarDynamics
    "The distance between the reference point and rear axis, positive means forward"
    ref_pos::ℝ
    ctrl_interval::ℝ
    "Propotional gain"
    kp::ℝ
    "(Discrete) Integral gain"
    ki::ℝ = 0.0
    "(Discrete) Derivative gain"
    kd::ℝ = 0.0
end

function ref_point(ref_pos, s::CarX)
    x = s.x + ref_pos * cos(s.θ)
    y = s.y + ref_pos * sin(s.θ)
    @SVector[x, y]
end

ref_point(K::RefPointTrackControl, s::CarX) = ref_point(K.ref_pos, s)

function ref_point_v(ref_pos, car_l, s::CarX)
    d = ref_pos
    @unpack v, θ, ψ = s
    ω = ω_from_v_ψ(v, ψ, car_l)
    v̂_x = v * cos(θ) - sin(θ) * ω * d
    v̂_y = v * sin(θ) + cos(θ) * ω * d
    @SVector[v̂_x, v̂_y]
end

ref_point_v(K::RefPointTrackControl, s::CarX) = ref_point_v(K.ref_pos, K.dy.l, s)

function track_ref(
    K::RefPointTrackControl, ξ::SymbolMap, ŝ::CarX{R}, s::CarX{R}, t
)::CarU{R} where R
    # compute the desired ref point velocity
    p̂ = ref_point(K, ŝ)
    v_p̂ = ref_point_v(K, ŝ)

    track_refpoint(K, ξ, (p̂, v_p̂), s, t)
end

function track_refpoint(
    K::RefPointTrackControl, ξ::SymbolMap, (p̂, v_p̂), s::CarX{R}, t
)::CarU{R} where R
    ξ = submap(ξ, :track_refpoint)
    p = ref_point(K, s)
    v_p = let
        Δt = K.ctrl_interval
        ∫edt = K.ki == 0 ? zero(p̂) : K.ki * Δt * integral!(ξ, :integral, t, p̂ - p)
        dedt = K.kd == 0 ? zero(p̂) : K.kd / Δt * derivative!(ξ, :derivative, t, p̂ - p)
        K.kp * (p̂ - p) + ∫edt + v_p̂
    end
    # convert `v_p` back into the control `CarU`
    d = K.ref_pos
    θ = s.θ
    v, v_y = rotation2D(-θ) * v_p
    ω = v_y / d
    (v < 0) && (ω *= -1)
    u_from_v_ω(v, ω, K.dy)
end

@kwdef struct ConfigTrackControl <: TrackingControl
    dy::CarDynamics
    k1::ℝ
    k2::ℝ
    k3::ℝ
end

ref_point(::ConfigTrackControl, s::CarX) = ref_point(0.0, s)
ref_point_v(K::ConfigTrackControl, s::CarX) = ref_point_v(0.0, K.dy.l, s)

"""
Full configuration tracking control for a 2D Car.
"""
function track_ref(
    K::ConfigTrackControl, ξ, x̂::CarX{R}, x::CarX{R}, t
)::CarU{R} where R
    θ_d, x_d, y_d, v_d = x̂.θ, x̂.x, x̂.y, x̂.v
    w_d = ω_from_v_ψ(x̂.v, x̂.ψ, K.dy.l)
    θ_e = x.θ - θ_d
    Δx = x.x - x_d
    Δy = x.y - y_d
    x_e =  cos(θ_d) * Δx + sin(θ_d) * Δy
    y_e = -sin(θ_d) * Δx + cos(θ_d) * Δy

    co = cos(θ_e)
    ta = tan(θ_e)
    v̂ = (v_d - K.k1 * abs(v_d) * (x_e + y_e * ta)) / co
    ŵ = w_d - (K.k2 * v_d * y_e + K.k3 * abs(v_d) * ta) * co * co
    u_from_v_ω(v̂, ŵ, K.dy)
end

struct RefTrackCentralControl{TC <: TrackingControl,Tr} <: CentralControl{CarU{ℝ},SymbolMap}
    K::TC
    trajectories::FuncT{Tuple{ℕ,𝕋},CarX{ℝ},Tr}
end

function OneVision.control_one(
    ctrl::RefTrackCentralControl, ξ::SymbolMap, xs, zs, t::𝕋, id::ℕ
)
    x = xs[id]
    x̂ = ctrl.trajectories((id, t))
    ξ = submap(ξ, Symbol(id))
    track_ref(ctrl.K, ξ, x̂, x, t)::CarU
end

# === formation driving example ====
const Formation{R} = Vector{CarX{R}}

function rotate_formation(form::Formation{R}, α)::Formation{R} where R
    mat = rotation2D(α)
    function f(s::CarX{R})
        @unpack x, y, θ, v, ψ = s
        x, y = mat * @SVector([x, y])
        θ += α
        CarX(; x, y, θ, v, ψ)
    end
    f.(form)
end

@kwdef struct CollisionAvoidance
    "Force scale"
    scale::ℝ
    "Force goes to infinity below this distance"
    min_r::ℝ
    "Force drops to zero beyound this distance"
    max_r::ℝ
end

"""
Leader-follower formation control
"""
struct FormationControl{TC <: TrackingControl, F} <: CentralControl{CarU{ℝ}, SymbolMap}
    "formation(xs, t)::Formation{ℝ}"
    formation::F
    "TrackingControl"
    K::TC
    dy::CarDynamics
    avoidance::CollisionAvoidance
end

const repl_rotation = rotation2D(10°)
function repel_force(avoidance::CollisionAvoidance, center, pos)
    @unpack scale, min_r, max_r = avoidance
    r = pos - center
    d = norm(r)
    (d > max_r) && return zero(center)

    α = max(eps(), (d - min_r) / (max_r - min_r))
    mag = scale * (1/α)
    repl_rotation * (r / d) * mag
end

function to_formation_frame(ctrl::FormationControl, s_leader)
    K = ctrl.K
    rot = rotation2D(s_leader.θ)
    pos_offset = ref_point(K, s_leader)
    v_offset = ref_point_v(K, s_leader)
    ω = ω_from_v_ψ(s_leader.v, s_leader.ψ, ctrl.dy.l)

    function formpoint_to_refpoint(fp)
        pos = rot * fp + pos_offset
        x, y = fp
        v = rot * @SVector[-ω*y, ω*x] + v_offset
        (pos = pos, v = v)
    end

    formpoint_to_refpoint
end

function formation_controller(ctrl::FormationControl{RefPointTrackControl}, ξ, xs, zs, t)
    formpoint_to_refpoint = to_formation_frame(ctrl, xs[1])
    form = ctrl.formation(xs, zs, t)

    N = length(xs)
    function avoid_collision(i)
        ca = ctrl.avoidance
        pos = get_pos(xs[i]) # ref_point(ctrl.K, xs[i])
        sum(repel_force(ca, get_pos(xs[j]), pos) for j in 1:N if j != i)
    end

    function action(id)::CarU
        (id == 1) && return zs[1].c

        s = form[id]
        (p, v_p) = formpoint_to_refpoint(@SVector[s.x, s.y])
        ξi = submap(ξ, Symbol(id))
        v_o = avoid_collision(id)
        track_refpoint(ctrl.K, ξi, (p, v_p + v_o), xs[id], t)
    end

    action
end

"""
Compute the `CarX` that rotates around the given `center` at the specified angular
velocity `ω`
"""
function rotate_around(center, ω, l, pos)
    x, y = pos
    dir = pos - center
    θ = atan(dir[2], dir[1]) + sign(ω) * 90°
    r = norm(dir)
    v = r * abs(ω)
    ψ = ψ_from_v_ω(v, ω, l)
    CarX(x, y, θ, v, ψ)
end

function center_of_rotation(s::CarX, l)
    ω = ω_from_v_ψ(s.v, s.ψ, l)
    r = ω ≈ 0 ? sign(ω) * 1e10 : s.v / ω
    θ = s.θ
    center = get_pos(s) + @SVector[-sin(θ), cos(θ)] * r
    center, ω
end

function formation_controller(ctrl::FormationControl{ConfigTrackControl}, ξ, xs, zs, t)
    s_leader = xs[1]
    @unpack l = ctrl.dy
    θ = s_leader.θ
    center, ω = center_of_rotation(s_leader, l)

    form = ctrl.formation(xs, zs, t)
    s′ = form[1]
    rot = rotation2D(s_leader.θ - s′.θ)
    offset = get_pos(s_leader) - get_pos(s′)

    function action(id)::CarU
        (id == 1) && return zs[1].c

        s = form[id]
        pos = rot * get_pos(s) + offset
        ŝ = rotate_around(center, ω, l, pos)
        ξi = submap(ξ, Symbol(id))
        track_ref(ctrl.K, ξi, ŝ, xs[id], t)
    end

    action
end

function OneVision.control_one(ctrl::FormationControl, ξ::SymbolMap, xs, zs, t::𝕋, id::ℕ)
    formation_controller(ctrl, ξ, xs, zs, t)(id)::CarU
end

function OneVision.control_all(
    ctrl::FormationControl, ξ::SymbolMap, xs, zs, t::𝕋, ids:: AbstractVector{ℕ}
)
    f = formation_controller(ctrl, ξ, xs, zs, t)
    f.(ids)
end

function car_triangle(x, y, θ; len = 0.2, width = 0.06)
    base = Point2f0(x, y)
    dir = Point2f0(cos(θ), sin(θ))
    left = Point2f0(sin(θ), -cos(θ))
    p1 = base + dir * len
    p2 = base + left * width
    p3 = base - left * width
    [p1,p2,p3]
end

function with_alpha(color, alpha)
    @unpack r, g, b = color
    RGBAf0(r, g, b, alpha)
end


text:
## Laplacian

function _imfilter_inbounds!(r::AbstractResource, out, A::AbstractArray, L::Laplacian, border::NoPad, inds)
    TT = eltype(out) # accumtype(eltype(out), eltype(A))
    n = 2*length(L.offsets)
    R = CartesianRange(inds)
    @unsafe for I in R
        tmp = convert(TT, - n * A[I])
        for J in L.offsets
            tmp += A[I+J]
            tmp += A[I-J]
        end
        out[I] = tmp
    end
    out
end

## imgradients
"""
```
imgradients(img, [points], [method], [border])
```

Performs edge detection filtering in the N-dimensional array `img`.
Gradients are computed at specified `points` (or indexes) in the
array or everywhere.

Available methods for 2D images: `"sobel"`, `"prewitt"`, `"ando3"`, `"ando4"`,
                                 `"ando5"`, `"ando4_sep"`, `"ando5_sep"`.

Available methods for ND images: `"sobel"`, `"prewitt"`, `"ando3"`, `"ando4"`.

Border options:`"replicate"`, `"circular"`, `"reflect"`, `"symmetric"`.

If `points` is specified, returns a 2D array `G` with the
gradients as rows. The number of rows is the number of
points at which the gradient was computed and the number
of columns is the dimensionality of the array.

If `points` is ommitted, returns a tuple of arrays, each
of the same size of the input image: (gradx, grady, ...)
"""
function imgradients{T,N}(img::AbstractArray{T,N}, points::AbstractVector;
                          method::AbstractString="ando3", border::AbstractString="replicate")
    extent = size(img)
    ndirs = length(extent)
    npoints = length(points)

    # pad input image only on appropriate directions
    imgpad = _gradientpad(img, border)

    # gradient matrix
    Tret = typeof(zero(T)*zero(Float64))
    G = zeros(Tret, npoints, ndirs)

    for dir in 1:ndirs
        # kernel = centered difference + perpendicular smoothing
        if extent[dir] > 1
            kern = _directional_kernel(dir, extent, method)

            # compute gradient at specified points
            A = zeros(kern)
            shape = size(kern)
            for (k, p) in enumerate(points)
                icenter = CartesianIndex(ind2sub(extent, p))
                i1 = CartesianIndex(tuple(ones(Int, ndirs)...))
                for ii in CartesianRange(shape)
                    A[ii] = imgpad[ii + icenter - i1]
                end

                G[k,dir] = sum(kern .* A)
            end
        end
    end

    G
end

function imgradients{T,N}(img::AbstractArray{T,N}, kernelfun::Function=KernelFactors.ando3, border=Pad{:replicate}())
    extended = map(isextended, indices(img))
    # kern1 = kernelfun(extended, 1)
    # S = filter_type(img, kern1)
    _imgradients((), img, kernelfun, extended, border)
end
isextended(ind) = length(ind) > 1

# When all N gradients have been calculated, return the result
_imgradients{T,N}(G::NTuple{N}, img::AbstractArray{T,N}, kernelfun::Function, extent, border) = G

# Add the next dimension to G
function _imgradients{T,M,N}(G::NTuple{M}, img::AbstractArray{T,N}, kernelfun::Function, extended, border)
    d = M+1  # the dimension we're working on now
    kern = kernelfun(extended, d)
    _imgradients((G..., imfilter(img, kern, border)), img, kernelfun, extended, border)
end


text:
using Random

function earthquake()
    # explosion detector (using spiral coorindate system) # run this file from julia by typing:]
    # julia> using PyPlot
    # julia> include("earthquake.jl")
    # julia> earthquake()
    # define the coordinate system:
    S=5000 # number of points on the spiral
    x=zeros(S); y=zeros(S)
    for s=1:S
        theta=50*2*pi*s/S; r=s/S
        x[s]=r*cos(theta); y[s]=r*sin(theta)
    end
    plot(x,y,".")
    # define the locations of the detection stations on the surface
    # Also define what value on each sensor would be generated by an explostion at internal location s
    N=10 # number of stations
    x_sensor=zeros(N); y_sensor=zeros(N)
    v=zeros(S,N)
    for sensor=1:N
        theta_sensor=2*pi*sensor/N
        x_sensor[sensor]=cos(theta_sensor); y_sensor[sensor]=sin(theta_sensor)
        for s=1:S
            v[s,sensor]=value(x[s],y[s],x_sensor[sensor],y_sensor[sensor]) # explosion value
        end
    end
    sd=1 # standard deviation of the Gaussian noise
    # Make the explosion data:
    true_s=rand(1:S) # true location of the explosion
    theta=50*2*pi*true_s/S
    r=true_s/S
    x_true=r*cos(theta); y_true=r*sin(theta)
    # Get the noisy sensor values that will be observed for this explosion:
    val=zeros(N)
    val_clean=zeros(N) # unknown clean values (just for interest)
    for sensor=1:N
        val_clean[sensor]=value(x_true,y_true,x_sensor[sensor],y_sensor[sensor])
        val[sensor]=val_clean[sensor]+sd*randn()
    end
    figure()
    plot(1:N,val,label="noisy observed sensor measurements")
    plot(1:N,val_clean,label="clean (unknown) sensor measurements")
    legend()
    # Perform inference p(location|observed sensor values) given these sensor values
    logp=zeros(S)
    for s=1:S
        for sensor=1:N
            logp[s] += -0.5*(val[sensor]-v[s,sensor])^2/(sd^2) # Gaussian distribution
        end
    end
    p=exp.(logp.-maximum(logp)) # do exponentiation (and avoid over/underflow)
    p=p/sum(p) # normalise
    # plot the posterior and most likely location of the explosion:
    maxp,maxind =findmax(p)
    figure()
    for s=1:S
        plot(x[s],y[s],".",color=(1-(p[s]/maxp))*[1,1,1])
    end
    for theta=0:0.01:2*pi
        plot(cos(theta),sin(theta),".",color=[0,0,0])
    end
    for sensor=1:N
        plot(x_sensor[sensor],y_sensor[sensor],"o",color=[1,0,0])
    end
    plot(x_true,y_true,"rx",markersize=20,label="true")
    plot(x[maxind],y[maxind],"m+",markersize=20,label="estimated (most likely)")
    plot(sum(p.*x),sum(p.*y),"go",markersize=10,label="estimated (average)")
    legend()
end

function value(x_true,y_true,x_sensor,y_sensor)
    return 1/(0.1+ (x_true-x_sensor)^2 + (y_true-y_sensor)^2)
end

max_stars_repo_path: Earthqauke - Assignment 1/.ipynb_checkpoints/earthquake-checkpoint.jl


text:
"""
    ANLI

For details, see the [GitHub page](https://github.com/facebookresearch/anli) or read the [2019 paper](https://arxiv.org/pdf/1910.14599.pdf).

Available data:
```julia
ANLI.R1_train_jsonl()
ANLI.R1_dev_jsonl()
ANLI.R1_test_jsonl()
ANLI.R2_train_jsonl()
ANLI.R2_dev_jsonl()
ANLI.R2_test_jsonl()
ANLI.R3_train_jsonl()
ANLI.R3_dev_jsonl()
ANLI.R3_test_jsonl()
```
"""
module ANLI

import ...register_data
using DataDeps: @datadep_str, unpack

anli_file(file, files...) = joinpath(datadep"ANLI", "anli_v0.1", file, files...)

R1_train_jsonl() = anli_file("R1", "train.jsonl")
R1_dev_jsonl()   = anli_file("R1", "dev.jsonl")
R1_test_jsonl()  = anli_file("R1", "test.jsonl")
R2_train_jsonl() = anli_file("R2", "train.jsonl")
R2_dev_jsonl()   = anli_file("R2", "dev.jsonl")
R2_test_jsonl()  = anli_file("R2", "test.jsonl")
R3_train_jsonl() = anli_file("R3", "train.jsonl")
R3_dev_jsonl()   = anli_file("R3", "dev.jsonl")
R3_test_jsonl()  = anli_file("R3", "test.jsonl")

function __init__()
    register_data(
        "ANLI",
        """
        ANLI is the Adversarial Natural Language Inference Benchmark.

        ANLI is licensed under Creative Commons-Non Commercial 4.0.
        """,
        "https://dl.fbaipublicfiles.com/anli/anli_v0.1.zip",
        "16ac929a7e90ecf9093deaec89cc81fe86a379265a5320a150028efe50c5cde8",
        postfetch=unpack
    )
end

end # module


text:
include(pwd()*raw"\GRASPEL.jl")

vars = matread(pwd()*raw"\USPS.mat") #read in matfile with USPS data
fea = vars["fea"] #read in matfile matrix as a Vector
gnd = vec(vars["gnd"].-1) #Read in matfile gnd, labels in gnd are offset by 1?? Also convert to vector
#Graspel = 1 #do graspel iters? unused currently
kNN = 2 #for setting up initial kNNT
samplePerc = 1/100 #percentage of sorted nodes as candidate edges
edge_iter = .01/100 #edge budget per iteration (percentage of nodes)
sig = 1e3 #feature variance
embedTol = 10 #embedding distortion tolerance (should be > 1)
maxIterations = 100 #max number of GRASPEL maxIterations

Subset = [0:1:9;]
num_cluster = length(Subset)
num_eigs = 5*num_cluster
num_test = 10 #number of kmean clustering runs
num_res_test = 1000 #number of resistance TestResistance

#order dataset into subsets
num_node, dim = size(fea)
p = sortperm(gnd)
fea = fea[p,:]
#feaVec = [vec(fea[i,:]) for i = 1:num_node]#fea must be a vector of vectors
gnd = gnd[p]
#preprocessing
##center images by average value
for i = 1:num_node
    fea[i,:] = fea[i,:].-mean(fea[i,:])
end
##normalize
##Julia norm(A,2) returns Frobenius Norm... not the same as Matlab's default 2-norm
##max(svd(A)), in this case svdvals(A)[1] , should eventually test this with larger data sets
fea = (fea/svdvals(fea)[1])*sqrt(dim)

#Do GRASPEL
G, varRatio, costVec, numIter = FiedlerConstruct(fea, samplePerc, edge_iter, sig, maxIterations, kNN, embedTol, num_node, dim)

num_edge=(length(edges(G)));
density=num_edge/num_node
L_J = laplacian_matrix(G)
L_J = L_J + (1/(sig^2))*sparse(LinearAlgebra.I,L_J.m,L_J.n)
A_J = adjacency_matrix(G)
mmwrite(pwd()*raw"\Julia_L.mtx",L_J)
mmwrite(pwd()*raw"\Julia_A.mtx",A_J)
@info "Density: $density"

TestResistance = 1 #resistance for Reff, check against data Distances

max_stars_repo_name: janticev/GRASPEL

max_issues_repo_name: janticev/GRASPEL

max_forks_repo_name: janticev/GRASPEL


text:
module ModuleConstruct
using MultipleViewGeometry
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using MultipleViewGeometry.ModuleProjection
using StaticArrays, LinearAlgebra
export construct

include("construct_projectionmatrix.jl")
include("construct_fundamentalmatrix.jl")
include("construct_homographymatrix.jl")
include("construct_essentialmatrix.jl")
end


text:
## File with graph search algorithms (BFS, path, etc...)

# BFS algorithm for Network
function bfsPath(g::Network, flows = zeros(length(g.links)))
  n = length(g.tails) - 1
  colors = zeros(Int,n) # unvisited = 0; open = 1; closed = 2
  paths = zeros(Int,n,n)

  # Init
  colors[1] = 1
  queue = [(g.source,typemax(Float64),1)]

  while !isempty(queue)
    v, flow, endofpath = pop!(queue)
    paths[v,endofpath] = v
    if v == g.target
          return (paths[v,1:endofpath], flow)
    end
    for (id,e) in enumerate(g.links[g.tails[v]:(g.tails[v+1]-1)])
      x = e.cap - flows[g.tails[v] + id - 1]
      if (x > 0) && (colors[e.head] == 0)
        pathflow = min(flow,x)
        q = (e.head, pathflow, endofpath + 1)
        push!(queue,q)
        colors[e.head] = 1
        paths[e.head,:] = paths[v,:]
      end
    end
    colors[v] = 2
  end
  return ([0],0.)
end

function bfsMinCut(g::Network, flows = zeros(length(g.links)))
  n = length(g.tails) - 1
  colors = zeros(Int,n) # unvisited = 0; open = 1; closed = 2
  cutedges = Vector{Int}()

  # Init
  colors[g.source] = 1
  queue = [g.source]

  while !isempty(queue)
    v = pop!(queue)
#    println(g.links[g.tails[v]:(g.tails[v+1]-1)])
    for (id,e) in enumerate(g.links[g.tails[v]:(g.tails[v+1]-1)])
#      println("colors[e.head]",colors[e.head])
#      println("e.cap > flows[id]",e.cap > flows[g.tails[v] + id - 1])
      if (e.cap > flows[g.tails[v] + id - 1]) && (colors[e.head] == 0)
        colors[v]=1
        q = e.head
        push!(queue,q)
      end
    end
    colors[v] = 2
  end

#  println(colors)
#  printflow(g,flows)

  for (tail,v) in enumerate(g.tails[1:n])
    for (id,a) in enumerate(g.links[g.tails[tail]:(g.tails[tail+1]-1)])
      if (colors[tail] - colors[a.head] == 2)
#        push!(cutedges,g.links[id+v-1].sym)
       push!(cutedges,id+v-1)
      end
    end
  end
  return cutedges
end

function bfs(g::Network, request::Symbol, flows = zeros(length(g.links)))
  if request == :Path
    bfsPath(g,flows)
  elseif request == :MinCut
    bfsMinCut(g, flows)
  end
end


text:
"""
Variational inference for Bayesian logistic regression.
"""
module LogisticRegression
using Distributions

export sigmoid, sample_data, VI

function sigmoid(x)
    return 1.0 / (1.0 + exp.(-x[1]))
end

function bern_sample(mu)
    i = rand(Bernoulli(mu))
    val = zeros(2)
    val[i+1] = 1
    return val
end

"""
Sample data & parameter given covariance Sigma_w and inputs X.
"""
function sample_data(X, Sigma_w)
    N = size(X, 2)
    M = size(Sigma_w, 1)

    # sample parameters
    W = rand(MvNormal(zeros(M), Sigma_w))

    # sample data
    Y = [rand(Bernoulli(sigmoid(W'*X[:, n]))) for n in 1 : N]
    return Y, W
end

"""
Compute variational parameters.
"""
function VI(Y, X, M, Sigma_w, alpha, max_iter)
    function rho2sig(rho)
        return log.(1 + exp.(rho))
    end

    function compute_df_dw(Y, X, Sigma_w, mu, rho, W)
        M, N = size(X)
        term1 = (mu - W) ./ rho2sig(rho).^2
        term2 = inv(Sigma_w)*W
        term3 = 0
        for n in 1 : N
            term3 += -(Y[n] - sigmoid(W'*X[:,n])) * X[:,n]
        end
        return term1 + term2 + term3
    end

    function compute_df_dmu(mu, rho, W)
        return (W - mu) ./ rho2sig(rho).^2
    end

    function compute_df_drho(Y, X, Sigma_w, mu, rho, W)
        return -0.5*((W - mu).^2 - rho2sig(rho).^2) .* compute_dprec_drho(rho)
    end

    function compute_dprec_drho(rho)
        return 2 * rho2sig(rho) .^ (-3) .* (1 ./ (1+exp.(rho))).^2 .* (1 ./ (1+exp.(-rho)))
    end

    # diag gaussian for approximate posterior
    mu = randn(M)
    rho = randn(M) # sigma = log.(1 + exp.(rho))

    for i in 1 : max_iter
        # sample epsilon
        ep = rand(M)
        W_tmp = mu + log.(1 + exp.(rho)) .* ep

        # calculate gradient
        df_dw = compute_df_dw(Y, X, Sigma_w, mu, rho, W_tmp)
        df_dmu = compute_df_dmu(mu, rho, W_tmp)
        df_drho = compute_df_drho(Y, X, Sigma_w, mu, rho, W_tmp)
        d_mu = df_dw + df_dmu
        d_rho = df_dw .* (ep ./ (1+exp.(-rho))) + df_drho

        # update variational parameters
        mu = mu - alpha * d_mu
        rho = rho - alpha * d_rho
    end
    return mu, rho
end

end


text:
"""
meansbands_to_matrix(m, input_type, cond_type, output_vars;
    forecast_string = "", verbose = :low)

Reformat `MeansBands` object into matrices, and save to individual files.

### Inputs

- `m::AbstractModel`: model object
- `input_type::Symbol`: same as input into `forecast_one`
- `cond_type::Symbol`: same as input into `forecast_one`
- `output_vars::Symbol`: same as input into `forecast_one`

### Keyword Arguments

- `forecast_string::String`: identifies the forecast (if
  desired). Required if `input_type == :subset`
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`
"""
function meansbands_to_matrix(m::AbstractModel, input_type::Symbol,
                              cond_type::Symbol, output_vars::Vector{Symbol};
                              forecast_string::String = "", verbose::Symbol = :low)

    # Determine full set of output_vars necessary for plotting desired results
    output_vars = add_requisite_output_vars(output_vars)
    output_dir  = workpath(m, "forecast")

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println()
        info("Converting means and bands to matrices for input_type = $input_type, cond_type = $cond_type...")
        println("Start time: $(now())")
        println("Means and bands matrices will be saved in $output_dir")
    end

    for output_var in output_vars
        meansbands_to_matrix(m, input_type, cond_type, output_var;
                             forecast_string = forecast_string, verbose = verbose)
    end

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        println("\nConversion of means and bands complete: $(now())")
    end
end

function meansbands_to_matrix(m::AbstractModel, input_type::Symbol,
                              cond_type::Symbol, output_var::Symbol;
                              forecast_string::String = "", verbose::Symbol = :low)

    mb = read_mb(m, input_type, cond_type, output_var,
                 forecast_string = forecast_string)

    prod     = get_product(mb)   # history? forecast? shockdec?
    cl       = get_class(mb)     # pseudo? obs? state?
    condtype = get_cond_type(mb)

    ## Get name of file to write
    outfile = get_forecast_filename(m, input_type, cond_type,
                                    Symbol("mb_matrix_", output_var);
                                    pathfcn = workpath,
                                    forecast_string = forecast_string,
                                    fileformat = :h5)

    # Convert MeansBands objects to matrices
    means, bands = meansbands_to_matrix(mb)

    # Save to file
    h5open(outfile, "w") do file
        write(file, "means", means)
        write(file, "bands", bands)
    end

    if VERBOSITY[verbose] >= VERBOSITY[:high]
        println(" * Wrote $(basename(outfile))")
    end
end

"""
```
meansbands_to_matrix(mb::MeansBands)
```

Convert a `MeansBands` object to matrix form.
"""
function meansbands_to_matrix(mb::MeansBands)

    # extract useful info
    vars       = get_vars_means(mb)             # get names of variables
    nvars      = length(vars)                   # get number of variables

    tmp        = setdiff(names(mb.means), [:date])[1]
    T          = eltype(mb.means[tmp])          # type of elements of mb structure

    nperiods   = n_periods_means(mb)            # get number of periods
    prod       = get_product(mb)                # history? forecast? shockdec?
    inds       = mb.metadata[:indices]          # mapping from names of vars to indices
    bands_list = which_density_bands(mb)        # which bands are stored
    nbands     = length(bands_list)             # how many bands are stored

    # extract  matrices from MeansBands structure
    if prod in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                :bddforecast, :bddforecastut, :bddforecast4q,
                :dettrend, :trend]

        # construct means and bands arrays
        means = Array{T,2}(nvars, nperiods)
        bands = Array{T}(nbands, nvars, nperiods)

        # extract each series and place in arrays
        for series in setdiff(names(mb.means), [:date])
            ind = inds[series]
            means[ind,:] = convert(Array{T},mb.means[series])

            for (i,band) in enumerate(bands_list)  # these are ordered properly already
                bands[i,ind,:] = convert(Array{T},mb.bands[series][Symbol(band)])
            end
        end

    elseif prod in [:shockdec, :irf]

        shock_inds = mb.metadata[:shock_indices]
        nshocks = length(shock_inds)

        # construct means and bands arrays
        means = Array{T}(nvars, nperiods, nshocks)
        bands = Array{T}(nbands, nvars, nperiods, nshocks)

        for series in setdiff(names(mb.means), [:date])

            var, shock = DSGE.parse_mb_colname(series)
            ind = inds[var]

            shock_ind = shock_inds[shock]
            means[ind,:,shock_ind] = convert(Array{T}, mb.means[series])

            for (band_ind, band) in enumerate(bands_list)
                bands[band_ind, ind, :, shock_ind] =
                    convert(Array{T}, mb.bands[series][Symbol(band)])
            end
        end
    end

    # return matrix
    return means, bands
end


text:
using BenchmarkTools

"""
    lookandsay(n::String)

Perform the look and say algorithm once on the given string.
"""
function lookandsay(n::String)
    sequence = IOBuffer()
    i, c = 1, 1
    while i ≤ length(n)
        if i != length(n) && n[i] == n[i+1]
            c += 1
        else
            print(sequence, c, n[i])
            c = 1
        end
        i += 1
    end
    return newstring(take!(sequence))
end

"""
    lookandsayseq(n::Integer, startsequence::String)

Iteratively construct the nth look and say sequence.
"""
function lookandsayseq(n::Integer, startsequence::String)
    sequence = Array{String}(n)
    sequence[1] = startsequence
    if n > 1
        for i in 2:n
            sequence[i] = lookandsay(sequence[i-1])
        end
    end
    return sequence
end

N = parse(Int, ARGS[1])

println("--------------------------")
println(@sprintf "Look and say sequence: %d" N)
println("--------------------------")

@btime begin
    r = lookandsayseq(N, "1223334444")
end
    #println(r[N])


text:
module HypertextTools

# Write your package code here.

	using Weave, Markdown
    include("weavehugo.jl")
	export lazyhugo, cp2content, getdoc, defolder

    using InteractiveUtils: clipboard # this is the old df2latex.jl
    using DataFrames
    include("latextable.jl")
    export latextable

	include("renamefile.jl");
	export renamefile
end

max_stars_repo_name: CGRG-lab/HypertextTools

max_issues_repo_name: CGRG-lab/HypertextTools

max_forks_repo_name: CGRG-lab/HypertextTools


text:
using HPAT

#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.CGen.setCreateMain(true)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(3)

@acc hpat function test1(file_name)
    t1 = DataSource(DataTable{:userid=Int64, :val2=Float64}, HDF5, file_name)
    t1 = t1[:userid>2]
    return t1[:userid], t1[:val2]
end

using HDF5
using MPI
if MPI.Comm_rank(MPI.COMM_WORLD)==0
    h5write("test1_1.hdf5", "/userid", [1, 2, 3, 1, 2])
    h5write("test1_1.hdf5", "/val2", [1.1, 2.1, 3.1, 3.2, 1.9])
end
a,b = test1("test1_1.hdf5")
println(a)
println(b)

using Base.Test
@test a==[3]
@test_approx_eq b [3.1]

if MPI.Comm_rank(MPI.COMM_WORLD)==0
    rm("test1_1.hdf5")
end


text:
module LyceumExamples

import IJulia

function notebooks()
    IJulia.notebook(dir = joinpath(@__DIR__, "../notebooks"), detached = true)
end

end # module


text:
# import GitCommand: git
import GitHub

function git(f)
  return f("git")
end

function install(workflow::Workflow,
                 org::AbstractString;
                 token::Union{AbstractString, Nothing},
                 cc::AbstractVector{<:AbstractString})
    if token === nothing
        auth = nothing
    else
        auth = GitHub.authenticate(token)
    end
    if auth === nothing
        orgrepos, page_data = GitHub.repos(org)
    else
        orgrepos, page_data = GitHub.repos(org; auth = auth)
    end
    pkgs = Vector{String}(undef, 0)
    for r in orgrepos
        name = r.name
        endswith(name, ".jl") && push!(pkgs, name[1:(end-3)])
    end
    unique!(pkgs)
    sort!(pkgs)
    install(workflow,
            org,
            pkgs;
            cc = cc,
            token = token)
    return nothing
end

function install(workflow::Workflow,
                 org::AbstractString,
                 pkgs::AbstractVector{<:AbstractString};
                 token,
                 cc::AbstractVector{<:AbstractString})
    if token === nothing
        auth = nothing
    else
        auth = GitHub.authenticate(token)
    end
    cc_string = string("cc: ", join(string.("@",
                                            strip.(strip.(strip.(cc),
                                                          '@'))), " "))
    my_pr_branch_name = "massinstallaction/set-up-$(workflow.name)"
    my_pr_title = "MassInstallAction: Install the $(workflow.name) workflow on this repository"
    my_pr_body = "This pull request sets up the $(workflow.name) workflow on this repository. $(cc_string)"
    for pkg in pkgs
        if auth === nothing
            pkgrepo = GitHub.repo("$(org)/$(pkg).jl")
        else
            pkgrepo = GitHub.repo("$(org)/$(pkg).jl";
                                  auth = auth)
        end
        pkg_url_with_auth = _get_repo_url_with_auth(org, pkg; token = token)
        with_temp_dir() do tmp_dir
            git() do git
                cd(tmp_dir)
                run(`$(git) clone $(pkg_url_with_auth) REPO`)
                cd("REPO")
                run(`$(git) checkout -B $(my_pr_branch_name)`)
                workflows_directory = joinpath(pwd(), ".github", "workflows")
                mkpath(workflows_directory)
                cd(workflows_directory)
                for filename in workflow.files_to_delete
                    rm(filename; force = true, recursive = true)
                end
                for filename in keys(workflow.files_to_create)
                    file_content = workflow.files_to_create[filename]
                    rm(filename; force = true, recursive = true)
                    open(filename, "w") do io
                        print(io, file_content)
                    end
                end
                run(`$(git) add -A`)
                try
                    run(`$(git) commit -m "Automated commit made by MassInstallAction.jl"`)
                catch
                end
                try
                    run(`$(git) push --force origin $(my_pr_branch_name)`)
                catch
                end
                try
                    run(`$(git) push --force origin $(my_pr_branch_name)`)
                catch
                end
                params = Dict{String, String}()
                params["title"] = my_pr_title
                params["head"] = my_pr_branch_name
                params["base"] = pkgrepo.default_branch
                params["body"] = my_pr_body
                if auth === nothing
                    try
                        GitHub.create_pull_request(pkgrepo;
                                                   params = params)
                    catch
                    end
                else
                    try
                        GitHub.create_pull_request(pkgrepo;
                                                   params = params,
                                                   auth = auth)
                    catch
                    end
                end
            end
        end
    end
    return nothing
end


text:
export dvec, bow
import SparseArrays: sparsevec, sparse

function dvec(x::AbstractSparseVector)
    #DVEC{Symbol,Float64}(model.id2token[x.nzind[i]] => x.nzval[i] for i in nonzeroinds(x))
    DVEC{eltype(x.nzind), eltype(x.nzval)}(t => v for (t, v) in zip(x.nzind, x.nzval))
end

sparse2dvec(x) = dvec

function sparsevec(vec::DVEC{Ti,Tv}, m=0) where {Ti<:Integer,Tv<:Number}
    I = Ti[]
    F = Tv[]

    for (t, weight) in vec
        push!(I, t)
        push!(F, weight)
    end

    if m == 0
        sparsevec(I, F)
    else
        sparsevec(I, F, m)
    end
end

function sparse(cols::AbstractVector{S}, m=0) where S<:DVEC{Ti,Tv} where {Ti<:Integer,Tv<:Number}
    I = Ti[]
    J = Ti[]
    F = Tv[]

    for j in eachindex(cols)
        for (t, weight) in cols[j]
            push!(I, t)
            push!(J, j)
            push!(F, weight)
        end
    end
    if m == 0
        sparse(I, J, F)
    else
        sparse(I, J, F, m, length(cols))
    end
end

function bow(model::Model, x::AbstractSparseVector)
    DVEC{Symbol,eltype{x.nzval}}(model.id2token[t] => v for (t, v) in zip(x.nzind, x.nzval))
end

function bow(model::Model, x::DVEC{Ti,Tv}) where {Ti<:Integer,Tv<:Number}
    DVEC{Symbol,Tv}(model.id2token[t] => v for (t, v) in x)
end

function dvec(model::Model, x::DVEC{Symbol,Tv}, Ti=Int) where Tv<:Number
    DVEC{Ti,Tv}(model.tokens[t].id => v for (t, v) in zip(x.nzind, x.nzval))
end

function sparsevec(model::VectorModel, bow::DVEC{Symbol,Tv}, Ti=Int) where {Tv<:Number}
    I = Ti[]
    F = Tv[]

    for (sym, weight) in bow
        idfreq = get(model.tokens, sym, nothing)
        if idfreq === nothing
            continue
        end

        push!(I, idfreq.id)
        push!(F, weight)
    end

    sparsevec(I, F, model.m)
end

function sparsevec(model::EntModel, bow::DVEC{Symbol,Tv}, Ti=Int) where {Tv<:Number}
    I = Ti[]
    F = Tv[]

    for (sym, weight) in bow
        idweight = get(model.tokens, sym, nothing)
        if idweight === nothing
            continue
        end

        push!(I, idweight.id)
        push!(F, weight)
    end

    sparsevec(I, F, model.m)
end