text:
using StaticArrays, SparseArrays, LinearAlgebra

"""
    UniformBlockDiagonal{T}

Homogeneous block diagonal matrices.  `k` diagonal blocks each of size `m×m`
"""
struct UniformBlockDiagonal{T} <: AbstractMatrix{T}
    data::Array{T, 3}
    facevec::Vector{SubArray{T,2,Array{T,3}}}
end

function UniformBlockDiagonal(dat::Array{T,3}) where {T}
    UniformBlockDiagonal(dat,
        SubArray{T,2,Array{T,3}}[view(dat,:,:,i) for i in 1:size(dat, 3)])
end

Base.copy(A::UniformBlockDiagonal) = UniformBlockDiagonal(copy(A.data))

function Base.size(A::UniformBlockDiagonal)
    m, n, l = size(A.data)
    (l * m, l * n)
end

function Base.getindex(A::UniformBlockDiagonal{T}, i::Int, j::Int) where {T}
    Ad = A.data
    m, n, l = size(Ad)
    (0 < i ≤ l * m && 0 < j ≤ l * n) ||
        throw(IndexError("attempt to access $(l*m) × $(l*n) array at index [$i, $j]"))
    iblk, ioffset = divrem(i - 1, m)
    jblk, joffset = divrem(j - 1, n)
    iblk == jblk ? Ad[ioffset+1, joffset+1, iblk+1] : zero(T)
end

function LinearAlgebra.Matrix(A::UniformBlockDiagonal{T}) where {T}
    Ad = A.data
    m, n, l = size(Ad)
    mat = zeros(T, (m*l, n*l))
    @inbounds for k = 0:(l-1)
        kp1 = k + 1
        km = k * m
        kn = k * n
        for j = 1:n
            knpj = kn + j
            for i = 1:m
                mat[km + i, knpj] = Ad[i, j, kp1]
            end
        end
    end
    mat
end

"""
    RepeatedBlockDiagonal{T}

A block diagonal matrix consisting of `k` blocks each of which is the same `m×m` `Matrix{T}`.

This is the form of the `Λ` matrix from a `VectorFactorReTerm`.
"""
struct RepeatedBlockDiagonal{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
    data::S
    nblocks::Int

    function RepeatedBlockDiagonal{T,S}(data,nblocks) where {T,S<:AbstractMatrix{T}}
        new{T,S}(data, nblocks)
    end
end

function RepeatedBlockDiagonal(A::AbstractMatrix, nblocks::Integer)
    RepeatedBlockDiagonal{eltype(A), typeof(A)}(A, Int(nblocks))
end

function Base.size(A::RepeatedBlockDiagonal)
    m, n = size(A.data)
    nb = A.nblocks
    (m * nb, n * nb)
end

function Base.getindex(A::RepeatedBlockDiagonal{T}, i::Int, j::Int) where {T}
    m, n = size(A.data)
    nb = A.nblocks
    (0 < i ≤ nb * m && 0 < j ≤ nb * n) ||
        throw(IndexError("attempt to access $(nb*m) × $(nb*n) array at index [$i, $j]"))
    iblk, ioffset = divrem(i - 1, m)
    jblk, joffset = divrem(j - 1, n)
    iblk == jblk ? A.data[ioffset+1, joffset+1] : zero(T)
end

function LinearAlgebra.Matrix(A::RepeatedBlockDiagonal{T}) where T
    mat = zeros(T, size(A))
    Ad = A.data
    m, n = size(Ad)
    nb = A.nblocks
    for k = 0:(nb-1)
        km = k * m
        kn = k * n
        for j = 1:n
            knpj = kn + j
            for i = 1:m
                mat[km + i, knpj] = Ad[i, j]
            end
        end
    end
    mat
end

"""
    BlockedSparse{Tv, Ti}

A `SparseMatrixCSC` whose nonzeros form blocks of rows or columns or both.

# Members
* `cscmat`: `SparseMatrixCSC{Tv, Ti}` representation for general calculations
* `nzsasmat`: Matrix{Tv} `cscmat.nzval` as a matrix
* `rowblocks`: `Vector{Vector{SubArray{Tv,1,Vector{Tv}}}}` of row blocks of nonzeros
* `colblocks`: `Vector{StridedMatrix{Tv}}` of column blocks of nonzeros
"""
mutable struct BlockedSparse{Tv,Ti} <: AbstractMatrix{Tv}
    cscmat::SparseMatrixCSC{Tv,Ti}
    nzsasmat::Matrix{Tv}
    rowblocks::Vector{Vector{SubArray{Tv,1,Vector{Tv}}}}
    colblocks::Vector{StridedMatrix{Tv}}
end

Base.size(A::BlockedSparse) = size(A.cscmat)

Base.size(A::BlockedSparse, d) = size(A.cscmat, d)

Base.getindex(A::BlockedSparse, i::Integer, j::Integer) = getindex(A.cscmat, i, j)

LinearAlgebra.Matrix(A::BlockedSparse) = Matrix(A.cscmat)

SparseArrays.sparse(A::BlockedSparse) = A.cscmat

SparseArrays.nnz(A::BlockedSparse) = nnz(A.cscmat)

function Base.copyto!(L::BlockedSparse{T,I}, A::SparseMatrixCSC{T,I}) where {T,I}
    @argcheck(size(L) == size(A) && nnz(L) == nnz(A), DimensionMismatch)
    copyto!(nonzeros(L.cscmat), nonzeros(A))
    L
end

"""
    OptSummary

Summary of an `NLopt` optimization

# Fields

* `initial`: a copy of the initial parameter values in the optimization
* `lowerbd`: lower bounds on the parameter values
* `ftol_rel`: as in NLopt
* `ftol_abs`: as in NLopt
* `xtol_rel`: as in NLopt
* `xtol_abs`: as in NLopt
* `initial_step`: as in NLopt
* `maxfeval`: as in NLopt
* `final`: a copy of the final parameter values from the optimization
* `fmin`: the final value of the objective
* `feval`: the number of function evaluations
* `optimizer`: the name of the optimizer used, as a `Symbol`
* `returnvalue`: the return value, as a `Symbol`
* `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs
* `REML`: use the REML criterion for LMM fits

The latter field doesn't really belong here but it has to be in a mutable struct in case it is changed.
"""
mutable struct OptSummary{T <: AbstractFloat}
    initial::Vector{T}
    lowerbd::Vector{T}
    finitial::T
    ftol_rel::T
    ftol_abs::T
    xtol_rel::T
    xtol_abs::Vector{T}
    initial_step::Vector{T}
    maxfeval::Int
    final::Vector{T}
    fmin::T
    feval::Int
    optimizer::Symbol
    returnvalue::Symbol
    nAGQ::Integer           # doesn't really belong here but I needed some place to store it
    REML::Bool              # similarly, just needed a place to store this information
end
function OptSummary(initial::Vector{T}, lowerbd::Vector{T},
    optimizer::Symbol; ftol_rel::T=zero(T), ftol_abs::T=zero(T), xtol_rel::T=zero(T),
    initial_step::Vector{T}=T[]) where T <: AbstractFloat
    OptSummary(initial, lowerbd, T(Inf), ftol_rel, ftol_abs, xtol_rel, zero(initial),
        initial_step, -1, copy(initial), T(Inf), -1, optimizer, :FAILURE, 1, false)
end

function Base.show(io::IO, s::OptSummary)
    println(io, "Initial parameter vector: ", s.initial)
    println(io, "Initial objective value:  ", s.finitial)
    println(io)
    println(io, "Optimizer (from NLopt):   ", s.optimizer)
    println(io, "Lower bounds:             ", s.lowerbd)
    println(io, "ftol_rel:                 ", s.ftol_rel)
    println(io, "ftol_abs:                 ", s.ftol_abs)
    println(io, "xtol_rel:                 ", s.xtol_rel)
    println(io, "xtol_abs:                 ", s.xtol_abs)
    println(io, "initial_step:             ", s.initial_step)
    println(io, "maxfeval:                 ", s.maxfeval)
    println(io)
    println(io, "Function evaluations:     ", s.feval)
    println(io, "Final parameter vector:   ", s.final)
    println(io, "Final objective value:    ", s.fmin)
    println(io, "Return code:              ", s.returnvalue)
end

function NLopt.Opt(optsum::OptSummary)
    lb = optsum.lowerbd

    opt = NLopt.Opt(optsum.optimizer, length(lb))
    NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective
    NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective
    NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values
    if length(optsum.xtol_abs) == length(lb)  # not true for the second optimization in GLMM
        NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values
    end
    NLopt.lower_bounds!(opt, lb)
    NLopt.maxeval!(opt, optsum.maxfeval)
    if isempty(optsum.initial_step)
        optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb))
    else
        NLopt.initial_step!(opt, optsum.initial_step)
    end
    opt
end

abstract type MixedModel{T} <: RegressionModel end # model with fixed and random effects

"""
    LinearMixedModel

Linear mixed-effects model representation

## Fields

* `formula`: the formula for the model
* `trms`: a `Vector` of `AbstractTerm` types representing the model.  The last element is the response.
* `sqrtwts`: vector of square roots of the case weights.  Can be empty.
* `A`: an `nt × nt` symmetric `BlockMatrix` of matrices representing `hcat(Z,X,y)'hcat(Z,X,y)`
* `L`: a `nt × nt` `BlockMatrix` - the lower Cholesky factor of `Λ'AΛ+I`
* `optsum`: an [`OptSummary`](@ref) object

## Properties

* `θ` or `theta`: the covariance parameter vector used to form λ
* `β` or `beta`: the fixed-effects coefficient vector
* `λ` or `lambda`: a vector of lower triangular matrices repeated on the diagonal blocks of `Λ`
* `σ` or `sigma`: current value of the standard deviation of the per-observation noise
* `b`: random effects on the original scale, as a vector of matrices
* `u`: random effects on the orthogonal scale, as a vector of matrices
* `lowerbd`: lower bounds on the elements of θ
* `X`: the fixed-effects model matrix
* `y`: the response vector
"""
struct LinearMixedModel{T <: AbstractFloat} <: MixedModel{T}
    formula::Formula
    trms::Vector
    sqrtwts::Vector{T}
    A::BlockMatrix{T}            # cross-product blocks
    L::LowerTriangular{T,BlockArray{T,2,AbstractMatrix{T}}}
    optsum::OptSummary{T}
end

function normalized_variance_cumsum(Λ)
    vars = cumsum(abs2.(svdvals(Λ)))
    vars ./ vars[end]
end

function Base.getproperty(m::LinearMixedModel, s::Symbol)
    if s ∈ (:θ, :theta)
        getθ(m)
    elseif s ∈ (:β, :beta)
        fixef(m)
    elseif s ∈ (:λ, :lambda)
        getΛ(m)
    elseif s ∈ (:σ, :sigma)
        sdest(m)
    elseif s == :b
        ranef(m)
    elseif s == :u
        ranef(m, uscale = true)
    elseif s == :lowerbd
        m.optsum.lowerbd
    elseif s == :X
        m.trms[end - 1].x
    elseif s == :y
        vec(m.trms[end].x)
    elseif s == :rePCA
        normalized_variance_cumsum.(getΛ(m))
    else
        getfield(m, s)
    end
end

Base.setproperty!(m::LinearMixedModel, s::Symbol, y) = s == :θ ? setθ!(m, y) : setfield!(m, s, y)

Base.propertynames(m::LinearMixedModel, private=false) =
    (:formula, :trms, :A, :L, :optsum, :θ, :theta, :β, :beta, :λ, :lambda, :σ, :sigma, :b, :u, :lowerbd, :X, :y, :rePCA)

struct RaggedArray{T,I}
    vals::Vector{T}
    inds::Vector{I}
end
function Base.sum!(s::AbstractVector{T}, a::RaggedArray{T}) where T
    for (v, i) in zip(a.vals, a.inds)
        s[i] += v
    end
    s
end

"""
    GeneralizedLinearMixedModel

Generalized linear mixed-effects model representation

# Fields
- `LMM`: a [`LinearMixedModel`](@ref) - the local approximation to the GLMM.
- `β`: the fixed-effects vector
- `β₀`: similar to `β`. Used in the PIRLS algorithm if step-halving is needed.
- `θ`: covariance parameter vector
- `b`: similar to `u`, equivalent to `broadcast!(*, b, LMM.Λ, u)`
- `u`: a vector of matrices of random effects
- `u₀`: similar to `u`.  Used in the PIRLS algorithm if step-halving is needed.
- `resp`: a `GlmResp` object
- `η`: the linear predictor
- `wt`: vector of prior case weights, a value of `T[]` indicates equal weights.
The following fields are used in adaptive Gauss-Hermite quadrature, which applies
only to models with a single random-effects term, in which case their lengths are
the number of levels in the grouping factor for that term.  Otherwise they are
zero-length vectors.
- `devc`: vector of deviance components
- `devc0`: vector of deviance components at offset of zero
- `sd`: approximate standard deviation of the conditional density
- `mult`: multiplier

# Properties

In addition to the fieldnames, the following names are also accessible through the `.` extractor

- `theta`: synonym for `θ`
- `beta`: synonym for `β`
- `σ` or `sigma`: common scale parameter (value is `NaN` for distributions without a scale parameter)
- `lowerbd`: vector of lower bounds on the combined elements of `β` and `θ`
- `formula`, `trms`, `A`, `L`, and `optsum`: fields of the `LMM` field
- `X`: fixed-effects model matrix
- `y`: response vector

"""
struct GeneralizedLinearMixedModel{T <: AbstractFloat} <: MixedModel{T}
    LMM::LinearMixedModel{T}
    β::Vector{T}
    β₀::Vector{T}
    θ::Vector{T}
    b::Vector{Matrix{T}}
    u::Vector{Matrix{T}}
    u₀::Vector{Matrix{T}}
    resp::GlmResp
    η::Vector{T}
    wt::Vector{T}
    devc::Vector{T}
    devc0::Vector{T}
    sd::Vector{T}
    mult::Vector{T}
end

function Base.getproperty(m::GeneralizedLinearMixedModel, s::Symbol)
    if s == :theta
        m.θ
    elseif s == :beta
        m.β
    elseif s ∈ (:λ, :lambda)
        getΛ(m)
    elseif s ∈ (:σ, :sigma)
        sdest(m)
    elseif s == :lowerbd
        m.LMM.optsum.lowerbd
    elseif s ∈ (:formula, :trms, :A, :L, :optsum)
        getfield(m.LMM, s)
    elseif s == :X
        m.LMM.trms[end - 1].x
    elseif s == :y
        vec(m.LMM.trms[end].x)
    else
        getfield(m, s)
    end
end

function Base.setproperty!(m::GeneralizedLinearMixedModel, s::Symbol, y)
    if s ∈ (:θ, :theta)
        setθ!(m, y)
    elseif s ∈ (:β, :beta)
        setβ!(m, y)
    elseif s ∈ (:βθ, :betatheta)
        setβθ!(m, y)
    else
        setfield!(m, s, y)
    end
end

Base.propertynames(m::GeneralizedLinearMixedModel, private=false) =
    (:theta, :beta, :λ, :lambda, :σ, :sigma, :X, :y, :lowerbd, fieldnames(typeof(m))..., fieldnames(typeof(m.LMM))...)

"""
    VarCorr

An encapsulation of information on the fitted random-effects
variance-covariance matrices.

# Members
* `σ`: a `Vector{Vector{T}}` of unscaled standard deviations
* `ρ`: a `Vector{Matrix{T}}` of correlation matrices
* `fnms`: a `Vector{Symbol}` of grouping factor names
* `cnms`: a `Vector{Vector{String}}` of column names
* `s`: the estimate of σ, the standard deviation of the per-observation noise.  When there is no scaling factor this value is `NaN`

The main purpose of defining this type is to isolate the logic in the show method.
"""
struct VarCorr{T}
    σ::Vector{Vector{T}}
    ρ::Vector{Matrix{T}}
    fnms::Vector{Symbol}
    cnms::Vector{Vector{String}}
    s::T
end
function VarCorr(m::MixedModel{T}) where T
    fnms = Symbol[]
    cnms = Vector{String}[]
    σ = Vector{T}[]
    ρ = Matrix{T}[]
    for trm in reterms(m)
        σi, ρi = stddevcor(trm)
        push!(σ, σi)
        push!(ρ, ρi)
        push!(fnms, trm.fnm)
        push!(cnms, trm.cnms)
    end
    VarCorr(σ, ρ, fnms, cnms, sdest(m))
end

cpad(s::String, n::Integer) = rpad(lpad(s, (n + textwidth(s)) >> 1), n)

function Base.show(io::IO, vc::VarCorr)
    # FIXME: Do this one term at a time
    fnms = copy(vc.fnms)
    stdm = copy(vc.σ)
    cor = vc.ρ
    cnms = reduce(append!, vc.cnms, init=String[])
    if isfinite(vc.s)
        push!(fnms, :Residual)
        push!(stdm, [1.])
        rmul!(stdm, vc.s)
        push!(cnms, "")
    end
    nmwd = maximum(map(textwidth, string.(fnms))) + 1
    write(io, "Variance components:\n")
    cnmwd = max(6, maximum(map(textwidth, cnms))) + 1
    tt = vcat(stdm...)
    vars = showoff(abs2.(tt), :plain)
    stds = showoff(tt, :plain)
    varwd = 1 + max(length("Variance"), maximum(map(textwidth, vars)))
    stdwd = 1 + max(length("Std.Dev."), maximum(map(textwidth, stds)))
    write(io, " "^(2+nmwd))
    write(io, cpad("Column", cnmwd))
    write(io, cpad("Variance", varwd))
    write(io, cpad("Std.Dev.", stdwd))
    any(s -> length(s) > 1, stdm) && write(io,"  Corr.")
    println(io)
    ind = 1
    for i in 1:length(fnms)
        stdmi = stdm[i]
        write(io, ' ')
        write(io, rpad(fnms[i], nmwd))
        write(io, rpad(cnms[ind], cnmwd))
        write(io, lpad(vars[ind], varwd))
        write(io, lpad(stds[ind], stdwd))
        ind += 1
        println(io)
        for j in 2:length(stdmi)
            write(io, " "^(1 + nmwd))
            write(io, rpad(cnms[ind], cnmwd))
            write(io, lpad(vars[ind], varwd))
            write(io, lpad(stds[ind], stdwd))
            ind += 1
            for k in 1:(j-1)
                @printf(io, "%6.2f", cor[i][j, k])
            end
            println(io)
        end
    end
end


text:
using TeaEM
using Documenter

DocMeta.setdocmeta!(TeaEM, :DocTestSetup, :(using TeaEM); recursive=true)

makedocs(;
    modules=[TeaEM],
    authors="Chen Huang <chen1huang2@gmail.com> and contributors",
    repo="https://github.com/chenspc/TeaEM.jl/blob/{commit}{path}#{line}",
    sitename="TeaEM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chenspc.github.io/TeaEM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/chenspc/TeaEM.jl",
)

ext: jl


text:
baremodule AgentClass
    const MOTORCYCLE = 1
    const CAR        = 2
    const TRUCK      = 3
    const PEDESTRIAN = 4
end

"""
    Vehicle definition which contains a class and a bounding box.
"""
struct VehicleDef
    class::Int64 # ∈ AgentClass
    length::Float64
    width::Float64
end
function VehicleDef(;
    class::Int=AgentClass.CAR,
    length::Float64=4.0,
    width::Float64=1.8,
    )

    VehicleDef(class, length, width)
end

const NULL_VEHICLEDEF = VehicleDef(AgentClass.CAR, NaN, NaN)

function Base.show(io::IO, d::VehicleDef)
    class = d.class == AgentClass.CAR ? "CAR" :
          d.class == AgentClass.MOTORCYCLE ? "MOTORCYCLE" :
          d.class == AgentClass.TRUCK ? "TRUCK" :
          d.class == AgentClass.PEDESTRIAN ? "PEDESTRIAN" :
          "UNKNOWN"
    @printf(io, "VehicleDef(%s, %.3f, %.3f)", class, d.length, d.width)
end
Base.write(io::IO, ::MIME"text/plain", def::VehicleDef) = @printf(io, "%d %.16e %.16e", def.class, def.length, def.width)
function Base.read(io::IO, ::MIME"text/plain", ::Type{VehicleDef})
    tokens = split(strip(readline(io)), ' ')
    class = parse(Int, tokens[1])
    length = parse(Float64, tokens[2])
    width = parse(Float64, tokens[3])
    return VehicleDef(class, length, width)
end

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
# Display
DISPLAY_MEAN = true
MEAN_STYLE = [nothing, :stroke, :fill, :gradient]
# 1. Horizontal (overlayed horizontal bars with opacity)
# 2. Vertical
#   - column to sort: Sample.start; all
#   - sort column by: sample mean, cumulative sum of previous values
#   - align at: beginning, ending, middle
# 3. Mean
#   - align at: middle
#   - source (location): cumulative sum of previous values (beginning) vs. sample (ending)
#   - visualization for each source:
#       - box,
#       - vertical histogram,
#       - horizontal histogram,
#       - violin plot (probably with opacity)


text:
export ExtremeValueDecision, extreme_value_decision, extreme_value_decision_one, demeaned, pref_scale, draw_gumbel_shocks

const EulerConst = 0.5772;

"""
   $(SIGNATURES)

Solve a decision problem with type I extreme value shocks

# Arguments
- value_ixM
      values of alternatives
      rows are types; columns are alternatives
      works with a single type
- prefScale
      scale of type I extreme value shocks
"""
struct ExtremeValueDecision{F <: AbstractFloat}
    prefScale :: F
    # Preference shocks demeaned?
    demeaned :: Bool
    dbg :: Bool
end

function ExtremeValueDecision(prefScale :: F, dbg :: Bool = false) where
    F <: AbstractFloat

    return ExtremeValueDecision(prefScale, true, dbg)
end

"""
	$(SIGNATURES)

Use demeaned Gumbel preference shocks?
"""
demeaned(d :: ExtremeValueDecision{F}) where F = d.demeaned;

"""
	$(SIGNATURES)

Returns the scale parameter of the Gumbel preference shocks.
"""
pref_scale(d :: ExtremeValueDecision{F}) where F = d.prefScale;

## Decision: One alternative
function extreme_value_decision(d :: ExtremeValueDecision{F},
   value_iV :: AbstractVector{F}) where F

    return extreme_value_decision(value_iV, pref_scale(d); demeaned = demeaned(d));
end

function extreme_value_decision(value_iV :: AbstractVector{F}, prefScale :: F;
   demeaned :: Bool = true) where F

    nTypes = length(value_iV);
    probV = ones(F, nTypes);
    valueV = copy(value_iV);
    if !demeaned
        valueV .+= (prefScale * F(EulerConst));
    end
    return probV, valueV
end

"""
   $(SIGNATURES)

Decision: Multiple alternatives.
Always returns a vector for `eVal_iV` (even with one agent).

OUT
   `prob_ixM`
      probability of choosing each option
   `eVal_iV`
      expected value for each type

TEST
   by simulation in `extreme_value_decision_test`
"""
function extreme_value_decision(d :: ExtremeValueDecision{F},
   value_ixM :: AbstractMatrix{F}) where  F <: Real

   return extreme_value_decision(value_ixM, pref_scale(d);
      demeaned = demeaned(d))
end

function extreme_value_decision(value_ixM :: AbstractMatrix{F}, prefScale :: F;
   demeaned :: Bool = true) where F <: Real

   nTypes, nx = size(value_ixM);
   prob_ixM = similar(value_ixM);
   eVal_iV = Vector{F}(undef, nTypes);
   for j = 1 : nTypes
      prob_ixM[j,:], eVal_iV[j] = extreme_value_decision_one(vec(value_ixM[j,:]),
         prefScale; demeaned = demeaned);
   end
    return prob_ixM, eVal_iV
end

# Input is a Tuple with Vectors for each type (across all alternatives).
# Output is a Vector containing a Vector of decision probabilities for each type.
# test this +++++
function extreme_value_decision(value_ixV :: NTuple{nTypes, Vector{F}},
   prefScale :: F;  demeaned :: Bool = true) where {nTypes, F}

   # nx = length(first(value_ixV));
   prob_ixV = Vector{Vector{F}}(undef, nTypes);
   eVal_iV = Vector{F}(undef, nTypes);
   for j = 1 : nTypes
      prob_ixV[j], eVal_iV[j] = extreme_value_decision_one(value_ixV[j],
         prefScale; demeaned = demeaned);
   end
    return prob_ixV, eVal_iV

end

"""
	$(SIGNATURES)

Extreme value decision for one individual.

This is probably not efficient. See efficient implementation for `logsumexp` in `StatsFuns.jl`. But here we also need the scaled exp(valueV) directly for the probabilities.

# Arguments
- valueV: can be a Vector or a Tuple with a value for each alternative.
"""
function extreme_value_decision_one(valueV, prefScale :: F;
   demeaned :: Bool = true) where F <: Real

   # Decision probability is log(sum(exp(V / prefScale)))
   # This needs to be nicely scaled to avoid overflow
   vMax = maximum(valueV) ./ prefScale .- F(4.0);

   # This uses a loop to ensure that probV is a Vector, even if valueV is a Tuple.
   # The following line is expensive
   probV = zeros(length(valueV));
   for (j, value) in enumerate(valueV)
      probV[j] = exp(value / prefScale - vMax);
   end

   # For each type: sum over alternatives
   expSum = sum(probV);

   # Prob of each choice
   probV ./= expSum;

   # Expected value
   eVal = prefScale * (vMax + log(expSum));

   if !demeaned
      eVal += prefScale * F(EulerConst);
   end

   return probV, eVal
end

# """
# Same with vector input (b/c a vector is not the same as an array with one row)
# """
# function extreme_value_decision(valueV :: Array{T,1},  prefScale :: T,
#     dbg :: Bool)   where T <: Real
#
#     return extreme_value_decision(collect(valueV'), prefScale, dbg)
# end

"""
	$(SIGNATURES)

Draw Gumbel shocks, optionally demeaned. Returns vector, even if `sizeV == 1`.
"""
function draw_gumbel_shocks(rng, prefScale :: F1, sizeV; demeaned :: Bool = true) where F1 <: Real

   drawM = prefScale .* rand(rng, Gumbel(zero(F1)),  sizeV...);
   if demeaned
       drawM .-= prefScale * F1(EconLH.EulerConst);
   end
   return drawM
end

# ----------------


text:
mutable struct StatObj

    vars   :: Dict
    ws     :: Dict

    function StatObj(
        vars :: Union{Nothing, Dict},
    )
        _vars = Dict()

        ws = Dict()

        for (k, v) in vars
            if typeof(v) <: Float64
                _vars[k] = 0.0
            elseif typeof(v) <: AbstractArray{Float64}
                _vars[k] = zeros(Float64, size(v)...)
            else
                throw(ErrorException("Not a subtype of Number or AbstractArray{Number}"))
            end

            ws[k] = 0.0
        end

        return new(_vars, ws)
    end

end

function addStatObj!(sobj::StatObj, vars::Dict; w::Float64=1.0)
    for (k, v) in vars
        if typeof(v) <: Number
            sobj.vars[k] += w * v
        else
            sobj.vars[k] .+= w .* v
        end
        sobj.ws[k] += w
    end
end

function normStatObj!(sobj::StatObj)
    for (k, v) in sobj.vars
        if typeof(v) <: Number
            sobj.vars[k] /= sobj.ws[k]
        else
            sobj.vars[k] ./= sobj.ws[k]
        end
        sobj.ws[k] = 1.0
    end

end

function zeroStatObj!(sobj::StatObj)
    for (k, v) in sobj.vars
        if typeof(v) <: Number
            sobj.vars[k] = 0.0
        else
            sobj.vars[k] .= 0.0
        end
        sobj.ws[k] = 0.0
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module HyperbolicPlane

using LinearFractionalTransformations, Plots, SimpleDrawing, AbstractLattices

export HObject, HPlane, equality_threshold

import Base: getindex, setindex!, isequal, length, ==, show, adjoint, -, in, +
import Base: angle, in, issubset
import AbstractLattices: ∨, ∧, dist

# this is faster than `abs(z)`
_mag(z::Number)::Real = real(z*z')

abstract type HObject end

THRESHOLD = 100   # Global default threshold = THRESHOLD*eps(1.0)

"""
`equality_threshold(val)` sets the "sloppiness" for equality checking.
The default value is 100. In general, two objects are equal if their
stored values are within `val * eps(1.0)`.

Calling `equality_threshold()` with no arguments returns the current
value.
"""
function equality_threshold()
    return THRESHOLD
end
function equality_threshold(newthresh::Real)
    global THRESHOLD = abs(newthresh)
end

setindex!(X::T, val, key::Symbol) where T<:HObject = setindex!(X.attr, val, key)
getindex(X::T, key) where T<:HObject = getindex(X.attr, key)

struct HPlane <: HObject
    attr::Dict{Symbol,Any}
    function HPlane()
        HP = new(Dict{Symbol,Any}())
        set_color(HP)
        set_thickness(HP)
        set_line_style(HP,:dot)
        return HP
    end
end

HPlane(X::HObject) = HPlane()   # copy constructor

(==)(A::HPlane, B::HPlane) = true
show(io::IO, A::HPlane) = print(io,"HPlane()")

include("point.jl")
include("segment.jl")
include("line.jl")
include("ray.jl")
include("triangle.jl")
include("polygon.jl")
include("circle.jl")
include("regular.jl")
include("horocycle.jl")

# straight things
HLinear = Union{HSegment,HLine,HRay}
# Round, fillable things
HRound  = Union{HPlane,HCircle,Horocycle}

export HLinear, HRound

include("meet.jl")
include("triangulate.jl")
include("container.jl")
include("tesselation.jl")

include("isometry.jl")
include("attributes.jl")
include("drawing.jl")

end #end of module


text:
@testset "Generation tools" begin

    t = (1,2)
    @test ϕ_tuple(t::NTuple) ≈ atan(2, 1)
    t = (1.23,129.23)
    @test ϕ_tuple(t::NTuple) ≈ atan(t[2], t[1])

    @test idx_to_dim([(1, 2) (3, 4); (5, 6) (7, 8)]) == [1 3; 5 7;;; 2 4; 6 8]

    @test idx_to_dim([(1, 2), (3, 4)]) == [1 2; 3 4]

end


text:
function optype(atype, btype)
    #If atype is signed, optype must be signed and at least big enough.
    if atype <: Integer || btype <: Integer
        if atype <: Signed || btype <: Signed
            p = promote_type(atype, btype)
            if p <: Integer
                return signed(p)
            else
                return p
            end
        else
            return promote_type(atype, btype)
        end
    else
        return promote_type(atype, btype)
    end
end

optype(::GBArray{T}, ::GBArray{U}) where {T, U} = optype(T, U)

function inferoutputtype(::GBArray{T}, ::GBArray{U}, op) where {T, U}
    t = optype(T, U)
    if op isa Tuple
        op = Semiring(op...)
    else
        op = BinaryOp(op)
    end
    return ztype(op, t)
end

function inferoutputtype(::GBArray{T}, ::GBArray{U}, op::AbstractOp) where {T, U}
    t = optype(T, U)
    return ztype(op, t)
end
function inferoutputtype(::GBArray{T}, op::AbstractOp) where {T}
    return ztype(op, T)
end
function inferoutputtype(::GBArray{T}, op) where {T}
    return ztype(UnaryOp(op), T)
end
function inferoutputtype(::GBArray{T}, ::AbstractTypedOp{Z}) where {T, Z}
    return Z
end
function inferoutputtype(::GBArray{T}, ::GBArray{U}, ::AbstractTypedOp{Z}) where {T, U, Z}
    return Z
end

function _handlenothings(kwargs...)
    return (x === nothing ? C_NULL : x for x in kwargs)
end
"""
    xtype(op::GrBOp)::DataType

Determine type of the first argument to a typed operator.
"""
function xtype end

"""
    ytype(op::GrBOp)::DataType

Determine type of the second argument to a typed operator.
"""
function ytype end

"""
    ytype(op::GrBOp)::DataType

Determine type of the output of a typed operator.
"""
function ztype end


text:
pushopt(options::Dict,expr::Expr) = Base.Meta.isexpr(expr,:(=)) && (options[expr.args[1]] = expr.args[2])

type MarkupInput
    codestart::Regex
    codeend::Regex
end

type ScriptInput
  doc_line::Regex
  doc_start::Regex
  opt_line::Regex
  opt_start::Regex
end

const input_formats = Dict{AbstractString, Any}(
        "noweb" => MarkupInput(r"^<<(.*?)>>=\s*$",
                    r"^@\s*$"),
        "markdown" => MarkupInput(
                      r"^[`~]{3,}(?:\{|\{\.|)julia(?:;|)\s*(.*?)(\}|\s*)$",
                      r"^[`~]{3,}\s*$"),
        "script" => ScriptInput(
          r"(^#'.*)|(^#%%.*)|(^# %%.*)",
          r"(^#')|(^#%%)|(^# %%)",
          r"(^#\+.*$)|(^#%%\+.*$)|(^# %%\+.*$)",
          r"(^#\+)|(^#%%\+)|(^# %%\+)")
        )

"""Detect the input format based on file extension"""
function detect_informat(source::AbstractString)
  ext = lowercase(splitext(source)[2])

  ext == ".jl" && return "script"
  ext == ".jmd" && return "markdown"
  return "noweb"
end

"""Read and parse input document"""
function read_doc(source::AbstractString, format=:auto)
    format == :auto && (format = detect_informat(source))
    document = @compat readstring(source)
    parsed = parse_doc(document, format)
    doc = WeaveDoc(source, parsed)
end

function parse_doc(document::AbstractString, format="noweb"::AbstractString)
  return parse_doc(document, input_formats[format])
end

"""Parse documents with Weave.jl markup"""
function parse_doc(document::AbstractString, format::MarkupInput)
  lines = split(document, "\n")

  codestart = format.codestart
  codeend = format.codeend
  state = "doc"

  docno = 1
  codeno = 1
  content = ""
  start_line = 0

  options = Dict()
  optionString = ""
  parsed = Any[]
  for lineno in 1:length(lines)
    line = lines[lineno]
    if (m = match(codestart, line)) != nothing && state=="doc"
      state = "code"
      if m.captures[1] == nothing
          optionString = ""
      else
          optionString=strip(m.captures[1])
      end
      #@show optionAbstractString
      options = Dict{Symbol,Any}()
      if length(optionString) > 0
          expr = parse(optionString)
          Base.Meta.isexpr(expr,:(=)) && (options[expr.args[1]] = expr.args[2])
          Base.Meta.isexpr(expr,:toplevel) && map(pushopt,fill(options,length(expr.args)),expr.args)
      end
      haskey(options, :label) && (options[:name] = options[:label])
      haskey(options, :name) || (options[:name] = nothing)
      #options = merge(rcParams[:chunk_defaults], options)
      #@show options
      chunk = DocChunk(content, docno, start_line)
      #chunk = @compat Dict{Symbol,Any}(:type => "doc", :content => content,
      #                                 :number => docno,:start_line => start_line)
      docno += 1
      start_line = lineno
      push!(parsed, chunk)
      content = ""
      continue
    end
    if ismatch(codeend, line) && state=="code"

      chunk = CodeChunk(content, codeno, start_line, optionString, options)
      #chunk = @compat Dict{Symbol,Any}(:type => "code", :content => content,
#                                   :number => codeno, :options => options,
#                                       :optionAbstractString => optionAbstractString,
#                                       :start_line => start_line)

      codeno+=1
      start_line = lineno
      content = ""
      state = "doc"
      push!(parsed, chunk)
      continue
    end

    if lineno == 1
      content *= line
    else
      content *= "\n" * line
    end
  end

  #Remember the last chunk
  if strip(content) != ""
    chunk = DocChunk(content, docno, start_line)
    #chunk = @compat Dict{Symbol,Any}(:type => "doc", :content => content,
    #                                 :number =>  docno, :start_line => start_line)
    push!(parsed, chunk)
  end
  return parsed
end

"""Parse .jl scripts with Weave.jl markup"""
function parse_doc(document::AbstractString, format::ScriptInput)
  lines = split(document, "\n")

  doc_line = format.doc_line
  doc_start = format.doc_start
  opt_line = format.opt_line
  opt_start = format.opt_start

  read = ""
  chunks = []
  docno = 1
  codeno = 1
  content = ""
  start_line = 1
  options = Dict{Symbol,Any}()
  optionString = ""
  parsed = Any[]
  state = "code"
  lineno = 1
  n_emptylines = 0

  for lineno in 1:length(lines)
    line = lines[lineno]
    if (m = match(doc_line, line)) != nothing && (m = match(opt_line, line)) == nothing
          line = replace(line, doc_start, "", 1)
      if startswith(line, " ")
          line = replace(line, " ", "", 1)
      end
      if state == "code"  && strip(read) != ""
          chunk = CodeChunk("\n" * rstrip(read), codeno, start_line, optionString, options)
          push!(parsed, chunk)
          codeno +=1
          read = ""
          start_line = lineno
      end
      state = "doc"
    elseif (m = match(opt_line, line)) != nothing
      start_line = lineno
      if state == "code" && strip(read) !=""
          chunk = CodeChunk("\n" * rstrip(read), codeno, start_line, optionString, options)
          push!(parsed, chunk)
          read = ""
          codeno +=1
      end
      if state == "doc" && strip(read) != ""
          (docno > 1) && (read = "\n" * read) # Add whitespace to doc chunk. Needed for markdown output
          chunk = DocChunk(read, docno, start_line)
          push!(parsed, chunk)
          read = ""
          docno += 1
      end

      optionString = replace(line, opt_start, "", 1)
      #Get options
      options = Dict{Symbol,Any}()
      if length(optionString) > 0
          expr = parse(optionString)
          Base.Meta.isexpr(expr,:(=)) && (options[expr.args[1]] = expr.args[2])
          Base.Meta.isexpr(expr,:toplevel) && map(pushopt,fill(options,length(expr.args)),expr.args)
      end
      haskey(options, :label) && (options[:name] = options[:label])
      haskey(options, :name) || (options[:name] = nothing)

      state = "code"
      continue
    elseif state == "doc" && strip(line) != "" && strip(read) != ""
      state = "code"
      (docno > 1) && (read = "\n" * read) # Add whitespace to doc chunk. Needed for markdown output
      chunk = DocChunk(read, docno, start_line)
      push!(parsed, chunk)
      options = Dict{Symbol,Any}()
      start_line = lineno
      read = ""
      docno += 1
    end
    read *= line * "\n"

    if strip(line) == ""
      n_emptylines += 1
    else
      n_emptylines = 0
    end
  end

  # Handle the last chunk
  if state == "code"
    chunk = CodeChunk("\n" * rstrip(read), codeno, start_line, optionString, options)
    push!(parsed, chunk)
  else
    chunk = DocChunk(read, docno, start_line)
    push!(parsed, chunk)
  end

  return parsed
end


text:
module Snappy

export
    compress,
    uncompress

@static if VERSION < v"0.6"
    ⊻ = $
end

include("varint.jl")
include("internal.jl")

"""
    compress(input::Vector{UInt8})

Compress a byte array using the snappy compression library.
"""
function compress(input::Vector{UInt8})
    length(input) > typemax(UInt32) && error("Input too large.")
    local sourcelen::UInt32 = length(input)

    # preallocate the output buffer, and resize down afterwards.
    output = zeros(UInt8, maxlength_compressed(sourcelen))
    outputindex = encode32!(output, 1, sourcelen)
    table = alloc_hashtable(sourcelen)

    for i in 0:K_BLOCK_SIZE:sourcelen
        fill!(table, 0xffff) # see compress_fragment! for explanation of using 0xffff
        outputindex = compress_fragment!(output, outputindex,
            input, i+1, min(i+K_BLOCK_SIZE, endof(input)), table)
    end

    return resize!(output, outputindex-1)
end
# for convenience
compress(input::String) = compress(Vector{UInt8}(input))

"""
    uncompress(input::Vector{UInt8})

Uncompress a snappy encoded byte array.
"""
function uncompress(input::Vector{UInt8})
    output_size, offset = length_uncompressed(input)
    output = zeros(UInt8, output_size)
    uncompressed_len = decompress_all_tags!(output, input, offset)
    (output_size != uncompressed_len-1) && error("Invalid input.")
    return output
end

# Compressed data can be defined as:
#    compressed := item* literal*
#    item       := literal* copy
#
# The trailing literal sequence has a space blowup of at most 62/60
# since a literal of length 60 needs one tag byte + one extra byte
# for length information.
#
# Item blowup is trickier to measure.  Suppose the "copy" op copies
# 4 bytes of data.  Because of a special check in the encoding code,
# we produce a 4-byte copy only if the offset is < 65536.  Therefore
# the copy op takes 3 bytes to encode, and this type of item leads
# to at most the 62/60 blowup for representing literals.
#
# Suppose the "copy" op copies 5 bytes of data.  If the offset is big
# enough, it will take 5 bytes to encode the copy op.  Therefore the
# worst case here is a one-byte literal followed by a five-byte copy.
# I.e., 6 bytes of input turn into 7 bytes of "compressed" data.
#
# This last factor dominates the blowup, so the final estimate is:
"""
    maxlength_compressed(sourcelen::Integer)

Return the maximum possible output size (in bytes) for a given input.
"""
function maxlength_compressed(sourcelen::Integer)
    return 32 + sourcelen + (sourcelen ÷ 6) # this is division, that op sure looks like plus though
end

"""
    length_uncompressed(input::Vector{UInt8})

Read uncompressed length encoded as a varint32 from the start of the array.
"""
function length_uncompressed(input::Vector{UInt8})
    return parse32(input, start(input))
end

end


text:
###################################
## Simulation models ##
###################################

function count_active_agents(agents::Vector{Agent})
    r = 0
    @simd for i = eachindex(agents)
        @inbounds r += agents[i].active
    end
    r
end

"""
`simulation_run` run traffic simulation in two modes.

**Input parameters**
* `mode` : simulation model switch
    * `:base` : simulation with regular agents only
    * `:smart` : simulation with smart agents enabled
* `OSMmap` : OpenStreetMapX MapData object with road network data
* `inAgents` : set of agents created with generate_agents function
* `density_factor` : road length reserved for one vehicle
* `k_routes_dict` : dictionary with multiple shortest paths (values) between vertices (keys)
* `U` : period of weights updates
* `T` : distribution parameter in k-shortest path rerouting
* `k` : number of fastest routes generated in rerouting function
* `debug` : debug messages switch
* `track_avg_speeds` : return average speeds on edges if true
"""
function simulation_run(mode::Symbol,
                        OSMmap::MapData,
                        inAgents::Vector{Agent},
                        density_factor::Float64 = 5.0,
                        k_routes_dict::Dict{Tuple{Int,Int},Array{Vector{Int}}}=
                        Dict{Tuple{Int,Int},Array{Vector{Int}}}(),
                        U::Int = 100,
                        T::Float64 = 1.0,
                        k::Int = 3,
                        debug::Bool = false;
                        track_avg_speeds::Bool = false)
    if !in(mode, [:base,:smart]) error("Wrong mode specified.") end
    Agents = deepcopy(inAgents) #Creating working copy of agents
    Agents[1].start_time = 0.0
    Agents[1].active = true
    max_densities, max_speeds = traffic_constants(OSMmap, density_factor) #Traffic characteristic constants
    densities, speeds = init_traffic_variables(OSMmap, Agents) #Traffic characteristic variables
    update_weights!(speeds, densities, max_densities, max_speeds) #Initial speeds update
    curr_edges = [a.edge for a in Agents] #Initialize edges list
    if track_avg_speeds
        avg_speeds = deepcopy(speeds)
        tick = 1
    end
    #Initialize simulation variables
    simtime = 0.0
    steps = 0
    active = length(Agents)
    #Calculate initial time to nearest junction for all agents
    times_to_edge = next_edge(Agents, speeds, OSMmap.w)
    times_to_start = [a.start_time for a in Agents] #Get agents starting times
    times_to_start[1] = Inf
    event_update = Inf
    #Loop until all agents are deactivated
    while active != 0 || !all(times_to_start.==Inf)
        steps += 1
        #Calculate next event time and update sim clock
        event_edge, IDe = findmin(times_to_edge)
        event_start, IDs = findmin(times_to_start)
        if mode == :smart
            event_update = next_update = (simtime ÷ U + 1) * U - simtime
        end
        min_event, event_type = findmin([event_edge, #1
                                        event_start, #2
                                        event_update]) #3
        simtime += min_event #Increase total simulation time
        times_to_edge .-= min_event #Move agents forward to event time
        times_to_start .-= min_event
        if event_type == 1 #Switching edges
            vAgent = Agents[IDe] #Pick agent for which event is occuring
            #Process agent for which event is occuring
            changed_edges = update_event_agent!(vAgent, simtime, densities, OSMmap.v)
            length(changed_edges) == 1 ? curr_edges[IDe] = (0,0) : curr_edges[IDe] = changed_edges[2]
            #Update speeds and correct events time
            update_weights_and_events!(Agents, curr_edges, times_to_edge, speeds,
                changed_edges, densities, max_densities, max_speeds)
            vAgent.active ? times_to_edge[IDe] = next_edge(vAgent, speeds, OSMmap.w) : times_to_edge[IDe] = Inf
        elseif event_type == 2 #New agent joins traffic
            Agents[IDs].active = true
            times_to_start[IDs] = Inf
            times_to_edge[IDs] = next_edge(Agents[IDs], speeds, OSMmap.w)
            c_edge = Agents[IDs].edge
            haskey(densities, c_edge) ? densities[c_edge] += 1 : densities[c_edge] = 1
        elseif event_type == 3 #Velocities update
            if debug
                update = Int(simtime ÷ U)
                smartAgents = sum([a.smart*a.active for a in Agents])
                print("Update: $(update) | Smart agents: $(smartAgents) ")
            end
            #Process rerouting for smart agents
            for a in Agents
                a.smart && a.active &&
                k_shortest_path_rerouting!(OSMmap, k_routes_dict, a, speeds, max_speeds, k, T, U)
            end
            debug && println("Finished")
        end
        #Update average speeds every 30 secs
        if track_avg_speeds && simtime > tick*30
            #Iterative mean
            tick +=1
            avg_speeds = ((tick-1)/tick)*avg_speeds+(1/tick)*speeds
        end
        active = count_active_agents(Agents)
    end
    times = getfield.(Agents,:travel_time)
    if track_avg_speeds
        output_tuple = (Steps = steps,
                    Simtime = simtime,
                    TravelTimes = times,
                    AvgSpeeds = avg_speeds)
    else
        output_tuple = (Steps = steps,
                    Simtime = simtime,
                    TravelTimes = times)
    end
    return output_tuple
end

"""
`gather_statistics` return various measures describing TMS model quality

**Input parameters**
* `smart_ind` : vector with smart agent indicators
* `times_base` : agents travelling time in base scenario
* `times_smart` : agents travelling time in smart scenario
"""
function gather_statistics(smart_ind::BitArray{1},
                            times_base::Vector{Float64},
                            times_smart::Vector{Float64})
    overall_time = (sum(times_base) - sum(times_smart))/sum(times_base)
    smart_time = (sum(times_base[smart_ind]) - sum(times_smart[smart_ind]))/sum(times_base[smart_ind])
    other_time = (sum(times_base[.!smart_ind]) - sum(times_smart[.!smart_ind]))/sum(times_base[.!smart_ind])
    statistics_tuple = (
        overall_time = round(overall_time, digits=3),
        smart_time = round(smart_time, digits=3),
        other_time = round(other_time, digits=3),
        avg_base = round(StatsBase.mean(times_base), digits=2),
        avg_overall_V2I = round(StatsBase.mean(times_smart), digits=2),
        avg_smart_V2I = round(StatsBase.mean(times_smart[smart_ind]), digits=2),
        avg_regular_V2I = round(StatsBase.mean(times_smart[.!smart_ind]), digits=2)
                        )
    return statistics_tuple
end

"""
`run_parameter_analysis` run base and V2I simulation scenario to compare agents performance, returns comparison statistics

**Input parameters**
* `GridElement` : element from parameters grid to be used in simulation
* `ParamGrid` : array with parameters values
* `Start` : vector of areas from which agents randomly pick starting point
* `End` : vector of areas from which agents randomly pick ending point
* `density_factor` : road length reserved for one vehicle
* `K_Paths_Dict` : dictionary with multiple shortest paths (values) between vertices (keys)
* `mapdata` : OpenStreetMapX MapData object with road network data
"""
function run_parameter_analysis(GridElement::Int,
                                ParamGrid::Vector{Tuple{Int,Float64,Int,Int,Float64,Int,Float64}},
                                Start::Vector{Rect}, End::Vector{Rect}, density_factor::Float64,
                                K_Paths_Dict::Dict{Tuple{Int,Int},Array{Vector{Int}}},
                                mapdata::OpenStreetMapX.MapData)
  startTime = time()
  p = ParamGrid[GridElement]
  #Generating agents
  Agents = generate_agents(mapdata, p[3], Start, End, p[2], p[6], p[5], p[7], K_Paths_Dict)
  #Running base simulation - no V2I system
  BaseOutput = simulation_run(:base, mapdata, Agents, density_factor)
  #Running simulation with smart agents - V2I system enabled
  SmartOutput = simulation_run(:smart, mapdata, Agents, density_factor,
                              K_Paths_Dict, p[4], p[5], p[6], false)
  step_stat = gather_statistics(getfield.(Agents,:smart),
                                      BaseOutput.TravelTimes,
                                      SmartOutput.TravelTimes)
  runtime = time()-startTime
  print("$(GridElement),$(p[1]),$(p[2]),$(p[3]),$(p[4]),$(p[5]),$(p[6]),$(p[7]),")
  print("$(round(runtime,digits = 2)),$(step_stat[1]),$(step_stat[2]),$(step_stat[3]),")
  println("$(step_stat[4]),$(step_stat[5]),$(step_stat[6]),$(step_stat[7])")
end

"""
`run_parameter_analysis` run base and V2I simulation scenario to compare agents performance, returns comparison statistics

**Input parameters**
* `GridElement` : element from parameters grid to be used in simulation
* `ParamGrid` : array with parameters values
* `AgentsVec` : vector with agents used in simulation
* `density_factor` : road length reserved for one vehicle
* `K_Paths_Dict` : dictionary with multiple shortest paths (values) between vertices (keys)
* `mapdata` : OpenStreetMapX MapData object with road network data
"""
function run_parameter_analysis(GridElement::Int,
                                ParamGrid::Vector{Tuple{Int,Float64,Int,Int,Float64,Int,Float64}},
                                AgentsVec::Vector{Agent}, density_factor::Float64,
                                K_Paths_Dict::Dict{Tuple{Int,Int},Array{Vector{Int}}},
                                mapdata::OpenStreetMapX.MapData)
  startTime = time()
  p = ParamGrid[GridElement]
  #Sample N agents from supplied vector and adjust smart agents number
  Agents = deepcopy(sample(AgentsVec, p[3], replace = false))
  N_int= Int(ceil(p[3]*p[2]))
  smart_ind = [trues(N_int); falses(p[3]-N_int)]
  for i in 1:p[3] Agents[i].smart = smart_ind[i] end
  #Running base simulation - no V2I system
  BaseOutput = simulation_run(:base, mapdata, Agents, density_factor)
  #Running simulation with smart agents - V2I system enabled
  SmartOutput = simulation_run(:smart, mapdata, Agents, density_factor,
                              K_Paths_Dict, p[4], p[5], p[6], false)
  step_stat = gather_statistics(getfield.(Agents,:smart),
                                      BaseOutput.TravelTimes,
                                      SmartOutput.TravelTimes)
  runtime = time()-startTime
  print("$(GridElement),$(p[1]),$(p[2]),$(p[3]),$(p[4]),$(p[5]),$(p[6]),$(p[7]),")
  print("$(round(runtime,digits = 2)),$(step_stat[1]),$(step_stat[2]),$(step_stat[3]),")
  println("$(step_stat[4]),$(step_stat[5]),$(step_stat[6]),$(step_stat[7])")
end


text:
module JuNix

# Write your package code here.

using Pkg: Pkg, Types, Operations, Registry
using Pkg.MiniProgressBars: MiniProgressBar, print_progress_bottom, show_progress, end_progress
using Pkg.Types: Context, EnvCache
using Artifacts: Artifacts
using Base: SHA1, UUID
using Base.BinaryPlatforms: AbstractPlatform, HostPlatform
using TOML
# using Git
using LibGit2
using Dates
# using ArgParse
using Printf
using SHA
using JSON3
using StructTypes: StructTypes, StructType
using Downloads: download

include("types.jl")
include("api.jl")

const GENERATED_FILENAME = "Depot.nix"

function collect_packages(ctx::Context, platforms::Vector{<:AbstractPlatform}, include_lazy::Bool)
    pkgs = PackageInfo[]
    deps = Pkg.dependencies(ctx.env)
    pkg_server = Pkg.pkg_server()
    for (uuid, dep) in deps
        (; name, version, tree_hash, is_tracking_path, is_tracking_repo, is_tracking_registry, source) = dep

        is_tracking_path && error("$name ($uuid) because it is tracking a local path: $source")
        tree_hash !== nothing || continue # stdlib

        pkg = PackageInfo(;
            name, uuid, version, tree_hash, source, is_tracking_path, is_tracking_repo, is_tracking_registry
        )
        load_repos!(pkg, ctx)
        load_archives!(pkg, ctx, pkg_server)
        load_artifacts!(pkg, ctx, pkg_server)
        push!(pkgs, pkg)
    end
    return pkgs
end

function load_repos!(pkg::PackageInfo, ctx::Context)
    if (pkg.is_tracking_registry || pkg.is_tracking_repo) && pkg.source !== nothing
        for url in Operations.find_urls(ctx.registries, pkg.uuid)
            @debug "Adding repo $repo_url to dependency $(pkg.name) ($(pkg.uuid))"
            push!(pkg.repos, Repo(; url))
        end
    end
    return pkg
end

function load_archives!(pkg::PackageInfo, ctx::Context, pkg_server)
    if pkg_server !== nothing
        url = "$pkg_server/package/$(pkg.uuid)/$(pkg.tree_hash)"
        @debug "Adding Pkg server archive $url to dependency $(pkg.name) ($(pkg.uuid))"
        push!(pkg.archives, Archive(; url))
    end
    for repo in pkg.repos
        url = Operations.get_archive_url_for_version(repo.url, pkg.tree_hash)
        url !== nothing && push!(pkg.archives, Archive(; url))
    end
    return pkg
end

function load_artifacts!(pkg::PackageInfo, ctx::Context, pkg_server)
    for (artifacts_toml, artifacts) in Operations.collect_artifacts(pkg.source)
        for (name, artifact) in artifacts
            archives = Archive[]
            if pkg_server !== nothing
                archive = Archive(; url="$pkg_server/artifact/$(pkg.tree_hash)")
                push!(archives, archive)
            end
            if haskey(artifact, "download")
                for d in artifact["download"]
                    archive = Archive(; url=d["url"], sha256=d["sha256"])
                    push!(archives, archive)
                end
            else
                @warn "Artifact missing 'download' section: $(artifact.name) $(artifact.tree_hash)"
            end
            augmented = Artifact(;
                name,
                tree_hash=artifact["git-tree-sha1"],
                lazy=get(artifact, "lazy", false),
                arch=get(artifact, "arch", nothing),
                os=get(artifact, "os", nothing),
                libc=get(artifact, "libc", nothing),
                libstdcxx_version=get(artifact, "libstdcxx_version", nothing),
                cxxstring_abi=get(artifact, "cxxstring_abi", nothing),
                julia_version=get(artifact, "julia_version", nothing),
                extra=Dict{String,Any}(
                    k => v for (k, v) in artifact if !(
                        k in (
                            "lazy",
                            "arch",
                            "os",
                            "libc",
                            "libstdcxx_version",
                            "cxxstring_abi",
                            "jucxxstring_abilia_version",
                            "download",
                            "git-tree-sha1",
                        )
                    )
                ),
                archives,
            )
            push!(pkg.artifacts, augmented)
        end
    end
end

function generate_depot(pkgs::Vector{PackageInfo}, ctx::Context)
    depot = Depot()
    (; registries, packages, artifacts) = depot
    for pkg in pkgs
        packages[pkg.tree_hash] = pkg
        for artifact in pkg.artifacts
            if haskey(artifacts, artifact.tree_hash)
                a = artifacts[artifact.tree_hash]
                b = artifact
                @warn "Artifact conflict detected. Attempting merge." a=artifacts[artifact.tree_hash].name b=artifact.name tree_hash=artifact.tree_hash
                for n in (:arch, :os, :libc, :libstdcxx_version, :cxxstring_abi, :julia_version)
                    if getfield(a, n) != getfield(b, n)
                        @warn "Merge failed" key=n a=getfield(a, n) b=getfield(b, n)
                        continue
                    end
                end
                a.lazy &= b.lazy
                for d in get(b, "download", [])
                    archive = Archive(; url=d["url"], sha256=d["sha256"])
                    push!(a.archives, archive)
                end
            else
                artifacts[artifact.tree_hash] = artifact
            end
        end
    end
    return depot
end

function resolve_depot(depot::Depot, ctx::Context)
    (; registries, packages, artifacts) = depot

    jobs = Channel{Union{Artifact,PackageInfo}}(ctx.num_concurrent_downloads)
    results = Channel(ctx.num_concurrent_downloads)

    to_download = eltype(jobs)[]
    foreach(pkg -> push!(to_download, pkg), values(packages))
    foreach(artifact -> push!(to_download, artifact), values(artifacts))

    @sync begin
        @async begin
            foreach(job -> put!(jobs, job), to_download)
        end
        for i in 1:ctx.num_concurrent_downloads
            @async begin
                for job in jobs
                    try
                        found = false
                        allerrors = CompositeException()
                        if job isa Artifact
                            res = resolve_archive!(job.archives)
                            res isa Exception && (append!(allerrors.errors, errors) && break)
                            put!(results, job => res)
                        elseif job isa PackageInfo
                            res = resolve_archive!(job.archives)
                            res isa Exception && (append!(allerrors.errors, errors) && break)
                            put!(results, job => res)
                            # for repo in job.repos
                            #     result = get_repo(repo)
                            #     result isa Exception && continue
                            #     found = true
                            #     put!(results, artifact => repo)
                            # end
                        end
                        isempty(allerrors) || put!(results, job => CompositeException(errors))
                    catch e
                        put!(results, job => e)
                    end
                end
            end
        end

        resolved = Depot()
        bar = MiniProgressBar(; indent=2, header = "Progress", color = Base.info_color(),
                                  percentage=false, always_reprint=true)
        bar.max = length(to_download)
        fancyprint = Pkg.can_fancyprint(ctx.io)

        try
            for i=1:length(to_download)
                job, res = take!(results)
                bar.current = i
                # fancyprint && print_progress_bottom(ctx.io)
                res isa Exception && @error "Error downloading $(job.name) $(job.tree_hash)" exception=res
                sha256 = res.sha256 === nothing ? nothing : to_nixhash(res.sha256)
                entry = res isa Archive ? (;type="archive", res.url, sha256, origsha256=res.sha256) : (;type="git", res.url, res.tree_hash)
                if job isa PackageInfo
                    entry = merge(entry, (;job.name, job.uuid , job.version))
                    resolved.packages[job.depot_path] = entry
                elseif job isa Artifact
                    entry = merge(entry, (;job.name, job.lazy, job.arch, job.os, job.libc, job.libstdcxx_version, job.cxxstring_abi, job.julia_version))
                    resolved.artifacts[job.depot_path] = entry
                end
                fancyprint && show_progress(ctx.io, bar)
            end
        finally
            fancyprint && end_progress(ctx.io, bar)
            close(jobs)
            return resolved
        end
    end
end

function resolve_archive!(archives::Vector{Archive})
    errors = Exception[]
    for archive in archives
        io = IOBuffer()
        try
            # download(archive.url, io)
            # archive.sha256 === nothing && (archive.sha256 = bytes2hex(sha256(io)))
            file = download(archive.url)
            archive.sha256 === nothing && (archive.sha256 = nix_hashfile(file))
            return archive
        catch e
            e isa InterruptException ? rethrow() : push!(errors, e)
            # push!(errors, CapturedException(e, catch_backtrace()))
        end
    end
    return CompositeException(errors)
end

function to_nixhash(hash::String, type="sha256")
    read(`nix hash to-sri --type sha256 $hash`, String)
end

function nix_hashfile(path::String)
    read(`nix hash file --sri $path`, String)
end

end # module


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

isym=txt2sym("openfst/sym.txt")

S = String
fst = txt2fst("openfst/A.fst", isym)
#println(fst)
@test size(fst) == (5,4)
fst = txt2fst("openfst/B.fst", "openfst/sym.txt")
#println(fst)
@test size(fst) == (4,3)
fst = txt2fst("openfst/C_trivial.fst", "openfst/sym.txt")
#println(fst)
@test size(fst) == (8,11)
fst = txt2fst("openfst/L.fst", "openfst/chars.txt", "openfst/words.txt")
@test size(fst) == (14,18)
#println(fst)
## TODO add many more tests

dot = fst2dot(fst)
# fst2pdf needs dot to be installed
if success(`which dot`)
  fst2pdf(fst,"fst.pdf")
end


text:
using StatsPlots
using Distributions
using Plots
plotlyjs()

mu = 3.14;
sigma = 1.0;
x = [mu - 5 * sigma:0.1:mu + 5 * sigma;];
N = Normal(mu, sigma^2);
p = plot(x, pdf.(N, x));
display(p)

max_stars_repo_name: jtschneider/plutonotebooks

max_issues_repo_name: jtschneider/plutonotebooks

max_forks_repo_name: jtschneider/plutonotebooks


text:
#=
Plots consumption, income and debt for the simple infinite horizon LQ
permanent income model with Gaussian iid income.

@author : Spencer Lyon

@date: 07/09/2014
=#

using PyPlot

const r = 0.05
const beta = 1.0 / (1.0 + r)
const T = 60
const sigma = 0.15
const mu = 1.0

function time_path()
    w = randn(T+1)
    w[1] =  0.0
    b = Array(Float64, T+1)
    for t=2:T+1
        b[t] = sum(w[1:t])
    end
    b .*= -sigma
    c = mu + (1.0 - beta) .* (sigma .* w .- b)
    return w, b, c
end

# == Figure showing a typical realization == #
function single_realization()
    fig, ax = subplots()

    ax[:grid]()
    ax[:set_xlabel]("Time")
    bbox = [0.0, 1.02, 1.0, 0.102]
    p_args = {:lw=> 2, :alpha => 0.7}

    w, b, c = time_path()
    ax[:plot](0:T, mu + sigma .* w, "g-", label="non-financial income"; p_args...)
    ax[:plot](0:T, c, "k-", label="consumption"; p_args...)
    ax[:plot](0:T, b, "b-", label="debt"; p_args...)
    ax[:legend](ncol=3, bbox_to_anchor=bbox, loc="upper left", mode="expand")

    return Void
end

# == Figure showing multiple consumption paths == #
function consumption_paths(n=250)  # n is number of paths
    fix, ax = subplots()
    ax[:grid]()
    ax[:set_xlabel]("Time")
    ax[:set_ylabel]("Consumption")
    colors = ["c-", "g-", "k-", "b-"]
    for i=1:n
        ax[:plot](0:T, time_path()[3], colors[rand(1:4)], lw=0.8, alpha=0.7)
    end
    return Void
end


text:
#
# @description Implements the RC4 stream cipher algorithm.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export Rc4Cipher

export RC4_LEGAL_KEY_SIZES
export key_size, reset!, encrypt!, decrypt!

# @description Defines the RC4 cipher data struture.
mutable struct Rc4Cipher <: StreamCipher
  key::Vector{UInt8}
  s::Vector{UInt8}
  i::UInt16
  j::UInt16

  # @description Initialize the RC4 cipher.
  #
  # @param {Vector{UInt8}} key The symmetric key.
  function Rc4Cipher(key::Vector{UInt8})
    return new(key, compute_rc4_ksa(key), 0, 0)
  end
end

# @description The DES cipher legal key sizes in bytes.
const RC4_LEGAL_KEY_SIZES = UInt64[ x for x in 1 : 32 ]

# @description Gets the cipher key size.
#
# @param {Rc4Cipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::Rc4Cipher)
  return length(self.key)
end

# @description Resets the state of the RC4 data structure.
#
# @param {Rc4Cipher} self The hash algorithm.
#
# @return The reference to the resetted hash algorithm.
function reset!(self::Rc4Cipher)
  self.s = compute_rc4_ksa(self.key)
  self.i = 0
  self.j = 0
end

# @description Encrypts the given plaintext block with the key set at initialization.
#
# @param {Rc4Cipher}     self      The cipher data struture.
# @param {Vector{UInt8}} plaintext The message to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext block.
function encrypt!(self::Rc4Cipher, plaintext::Vector{UInt8})
  return transform_rc4_bytes!(self, plaintext)
end

# @description Decrypts the given ciphertext block with the key set at initialization.
#
# @param {Rc4Cipher}     self       The cipher data struture.
# @param {Vector{UInt8}} ciphertext The message to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function decrypt!(self::Rc4Cipher, ciphertext::Vector{UInt8})
  return transform_rc4_bytes!(self, ciphertext)
end

# PRIVATE IMPLEMENTATION #######################################################

# @description Computes the RC4 key-scheduling algorithm (KSA).
#
# @param {Vector{UInt8}} key The symmetric key.
#
# @return {Vector{UInt8}} The keystream 'S'.
function compute_rc4_ksa(key::Vector{UInt8})
  if length(key) ∉ RC4_LEGAL_KEY_SIZES
    throw(ArgumentError("Invalid key length. $(length(key) * 8) bits provided."))
  end

  s = UInt8[ x for x in 0 : 255 ]
  j = 0

  for i in 0 : 255
    k = mod(i, length(key))
    j = mod(j + s[i + 1] + key[k + 1], 256)
    s[i + 1], s[j + 1] = s[j + 1], s[i + 1]
  end

  return s
end

# @description Computes the RC4 Pseudo-random generation algorithm (PRGA).
#
# @param {Rc4Cipher} self The cipher data struture.
#
# @return {UInt8} The random byte.
function compute_next_rc4_number!(self::Rc4Cipher)
  self.i = mod(self.i + 1, 256)
  self.j = mod(self.j + self.s[self.i + 1], 256)
  index  = mod(self.s[self.i + 1] + self.s[self.j + 1], 256)

  self.s[self.i + 1], self.s[self.j + 1] = self.s[self.j + 1], self.s[self.i + 1]

  return self.s[index + 1]
end

# @description Transforms a RC4 byte's stream.
#
# @param {Rc4Cipher}     self The cipher data struture.
# @param {Vector{UInt8}} data The message to transform.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function transform_rc4_bytes!(self::Rc4Cipher, data::Vector{UInt8})
  result = UInt8[]

  for x in data
    push!(result, x ⊻ compute_next_rc4_number!(self))
  end

  return result
end


text:
function test_diagnosticsteps(dev::Device=CPU(); nsteps=100, freq=1, ndata=ceil(Int, (nsteps+1)/freq))
  prob = Problem(nx=6, Lx=2π, dev=dev)
  getone(prob) = 1
  d = Diagnostic(getone, prob, nsteps=nsteps, freq=freq, ndata=ndata)
  stepforward!(prob, d, nsteps)
  expectedsteps = cat([0], freq:freq:nsteps, dims=1)
  d.steps[1:d.i] == expectedsteps
end

function test_scalardiagnostics(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2, nsteps=100, freq=1, ndata=ceil(Int, (nsteps+1)/freq))
  k1 = 2π/Lx
   τ = 1/(kappa*k1^2) # time-scale for diffusive decay
  dt = 1e-4 * τ # dynamics are resolved.

  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="ETDRK4", dev=dev)
  g = prob.grid

  ct(t) = exp(-kappa*k1^2*t)

  t1 = dt*nsteps
  c0 = @. sin(k1*g.x)
  c1 = @. ct(t1) * c0 # analytical solution

  set_c!(prob, c0)

  variance(prob) = parsevalsum2(prob.sol, prob.grid) / Lx
  d = Diagnostic(variance, prob, nsteps=nsteps, freq=freq, ndata=ndata)

  stepforward!(prob, d, nsteps)

  da = @. 0.5*ct(d[:t])^2
  isapprox(d[:data], da)
end

function test_basicdiagnostics(dev::Device=CPU(); nx=6, Lx=2π, kappa=1e-2)
  k1 = 2π/Lx
   τ = 1/(kappa*k1^2) # time-scale for diffusive decay
  dt = 1e-2 * τ # dynamics are resolved.

  prob = Problem(nx=nx, Lx=Lx, kappa=kappa, dt=dt, stepper="ETDRK4", dev=dev)
  g = prob.grid

  c0(x) = sin(k1*x)
  set_c!(prob, c0)

  getsol(prob) = prob.sol
  soldiag = Diagnostic(getsol, prob, nsteps=1)

  stepforward!(prob, soldiag, 1)

  isapprox(soldiag.data[soldiag.i], prob.sol)
end


text:
"""
    AbstractAllocator

Types implementing the `AbstractAllocator` interface need to define the following two
methods:

    allocate(allocator, bytes::Integer) -> (Ptr{Nothing}, Token)

Allocate `bytes` from `allocator`. Return type is a a tuple whose first element is a
pointer with type `Ptr{Nothing}` to the beginning of the allocated data and whose second
element is some token that can be passed to `free`.

    free(ptr::Ptr{Nothing}, token)

Complement to `allocate`: Takes the returned pointer and token and takes any steps
necessary to free the data.

In otherwords, the following should and clean up all allocated resources:

    free(allocate(allocator, bytes)...)

Example
=======
```julia-repl
julia> allocator = CachedArrays.AlignedAllocator()
CachedArrays.AlignedAllocator()

julia> ptr, token = CachedArrays.allocate(allocator, 100)
(Ptr{Nothing} @0x00005587e09cdf40, CachedArrays.AlignedAllocator())

julia> CachedArrays.free(ptr, token)
```

Note: If `token` will automatically handle resource reclamation using a finalizer, then
the type `NoopWrapper` can be wrapped around the token to turn `free` into a no-op.
"""
abstract type AbstractAllocator end

struct NoopWrapper{T}
    data::T
end
free(::Ptr{Nothing}, ::NoopWrapper) = nothing

#####
##### Memory Map Allocator
#####

struct MmapAllocator <: AbstractAllocator
    path::String
    persistent::Bool
end

MmapAllocator(dir::AbstractString) = MmapAllocator(dir, false)
PersistentMmapAllocator(path::AbstractString) = MmapAllocator(path, true)

function allocate(allocator::MmapAllocator, sz)
    # Decide if we're going to create a temporary mapping or use an existing mapping.
    path = allocator.path
    path = allocator.persistent ? path : tempname(path; cleanup = false)
    A = open(path; create = true, read = true, write = true) do io
        Mmap.mmap(io, Vector{UInt8}, sz)
    end

    # Remove the file - OS will keep it around as long as our process is running and
    # the cleanup once the process terminates.
    allocator.persistent || rm(path)
    ptr = Ptr{Nothing}(pointer(A))
    return ptr, NoopWrapper(A)
end

#####
##### DRAM Allocator
#####

struct AlignedAllocator <: AbstractAllocator end

allocate(::AlignedAllocator, sz) = aligned_alloc(sz)
free(ptr::Ptr{Nothing}, ::AlignedAllocator) = Libc.free(ptr)

function aligned_alloc(sz::Integer; alignment = UInt(64))
    ptr_ref = Ref(Ptr{Nothing}())
    ret = ccall(
        :posix_memalign,
        Cint,
        (Ptr{Ptr{Cvoid}}, Csize_t, Csize_t),
        ptr_ref,
        alignment,
        sz
    )
    @assert ret == 0
    return ptr_ref[], AlignedAllocator()
end


text:
function load_rdd_data(dataset::Symbol)
    if dataset == :lee08
        lee08_path = joinpath(dirname(@__FILE__), "..", "data", "lee08.feather")

        lee08 = Feather.read(lee08_path)
        ZsR = RunningVariable(lee08.margin ./ 100; cutoff = 0.0, treated = :≥)
        Ys = lee08.voteshare ./ 100
        dataset = RDData(Ys, ZsR)
    end

    return dataset
end


text:
using ModernGL
using Makie

function run_experiment(escreen::PsychoScreen)
    #disable renderloop
    Makie.opengl_renderloop[] = (_screen)->nothing

    stop = false
    #sceen =  get the screen for the scene
    while !stop
        GLFW.PollEvents()
        Makie.render_frame(screen)
        GLFW.SwapBuffers(Makie.to_native(screen))
        glFinish()
    end
    Makie.opengl_renderloop[] = Makie.renderloop # restore previous loop
    Makie.destroy!(screen)
end


text:
#==============================================================================#
# KinesisVideoMedia.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/kinesis-video-media-2017-09-30.normal.json
#==============================================================================#

__precompile__()

module KinesisVideoMedia

using AWSCore

"""
    using AWSSDK.KinesisVideoMedia.get_media
    get_media([::AWSConfig], arguments::Dict)
    get_media([::AWSConfig]; StartSelector=, <keyword arguments>)

    using AWSCore.Services.kinesis_video_media
    kinesis_video_media([::AWSConfig], "POST", "/getMedia", arguments::Dict)
    kinesis_video_media([::AWSConfig], "POST", "/getMedia", StartSelector=, <keyword arguments>)

# GetMedia Operation

Use this API to retrieve media content from a Kinesis video stream. In the request, you identify stream name or stream Amazon Resource Name (ARN), and the starting chunk. Kinesis Video Streams then returns a stream of chunks in order by fragment number.

**Note**
> You must first call the `GetDataEndpoint` API to get an endpoint to which you can then send the `GetMedia` requests.

When you put media data (fragments) on a stream, Kinesis Video Streams stores each incoming fragment and related metadata in what is called a "chunk." For more information, see . The `GetMedia` API returns a stream of these chunks starting from the chunk that you specify in the request.

The following limits apply when using the `GetMedia` API:

*   A client can call `GetMedia` up to five times per second per stream.

*   Kinesis Video Streams sends media data at a rate of up to 25 megabytes per second (or 200 megabits per second) during a `GetMedia` session.

# Arguments

## `StreamName = ::String`
The Kinesis video stream name from where you want to get the media content. If you don't specify the `streamName`, you must specify the `streamARN`.

## `StreamARN = ::String`
The ARN of the stream from where you want to get the media content. If you don't specify the `streamARN`, you must specify the `streamName`.

## `StartSelector = [ ... ]` -- *Required*
Identifies the starting chunk to get from the specified stream.
```
 StartSelector = [
        "StartSelectorType" => <required> "FRAGMENT_NUMBER", "SERVER_TIMESTAMP", "PRODUCER_TIMESTAMP", "NOW", "EARLIEST" or "CONTINUATION_TOKEN",
        "AfterFragmentNumber" =>  ::String,
        "StartTimestamp" =>  timestamp,
        "ContinuationToken" =>  ::String
    ]
```

# Returns

`GetMediaOutput`

# Exceptions

`ResourceNotFoundException`, `NotAuthorizedException`, `InvalidEndpointException`, `ClientLimitExceededException`, `ConnectionLimitExceededException` or `InvalidArgumentException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/kinesis-video-media-2017-09-30/GetMedia)
"""
@inline get_media(aws::AWSConfig=default_aws_config(); args...) = get_media(aws, args)

@inline get_media(aws::AWSConfig, args) = AWSCore.Services.kinesis_video_media(aws, "POST", "/getMedia", args)

@inline get_media(args) = get_media(default_aws_config(), args)

end # module KinesisVideoMedia

#==============================================================================#
# End of file
#==============================================================================#


text:
using QuantumESPRESSO
using Documenter

DocMeta.setdocmeta!(QuantumESPRESSO, :DocTestSetup, :(using QuantumESPRESSO); recursive=true)

makedocs(;
    modules=[QuantumESPRESSO],
    authors="Qi Zhang <singularitti@outlook.com>",
    repo="https://github.com/MineralsCloud/QuantumESPRESSO.jl/blob/{commit}{path}#{line}",
    sitename="QuantumESPRESSO.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MineralsCloud.github.io/QuantumESPRESSO.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MineralsCloud/QuantumESPRESSO.jl",
)


text:
immutable RK4Stepper
    propagation_step
end

immutable RK4StepperState
    rhs_function
    previous_pdim_value
    current_pdim_value
    previous_field
    current_field
    previous_step_times_deriv
    current_step_times_deriv
end

pdim(state::RK4StepperState) = state.current_pdim_value

function (stepper::RK4Stepper)(
        problem::IntegrationProblem, initial_state::AbstractArray, initial_pdim::Number)

    rhs = replace_differentials(problem.rhs)

    rhs_afunc = compile_wrapped_array_function(as_wrapped_array_function(
        rhs,
        [problem.unknown_field, problem.propagation_dimension],
        problem.transverse_dimensions))

    RK4StepperState(
        rhs_afunc,
        initial_pdim,
        initial_pdim,
        initial_state,
        initial_state,
        nothing,
        nothing)
end

function make_step(stepper::RK4Stepper, state::RK4StepperState)
    rhs = state.rhs_function
    field = state.current_field
    pdim = state.current_pdim_value
    step = stepper.propagation_step

    step_times_deriv_prev = step * rhs(field, pdim)
    k1 = step_times_deriv_prev
    k2 = step * rhs(field + k1 / 2, pdim + step / 2)
    k3 = step * rhs(field + k2 / 2, pdim + step / 2)
    k4 = step * rhs(field + k3, pdim + step)
    step_times_deriv = k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
    new_field = field + step_times_deriv

    RK4StepperState(
        state.rhs_function,
        pdim,
        pdim + stepper.propagation_step,
        field,
        new_field,
        step_times_deriv_prev,
        step_times_deriv)
end

function _interpolate_at(pdim, pdim1, pdim2, field1, field2, f1, f2)
    h = pdim2 - pdim1
    t = (pdim - pdim1) / h

    # Third-order approximation
    (
        (1 - t) * field1 + t * field2
        + t * (t - 1) * ((1 - 2 * t) * (field2 - field1) + (t - 1) * f1 + t * f2))
end

function interpolate_at(stepper::RK4Stepper, state::RK4StepperState, pdim_val)
    @assert state.previous_pdim_value <= pdim_val <= state.current_pdim_value

    if state.previous_pdim_value == pdim_val
        return state.previous_field
    end

    if state.current_pdim_value == pdim_val
        return state.current_field
    end

    _interpolate_at(
        pdim_val, state.previous_pdim_value, state.current_pdim_value,
        state.previous_field, state.current_field,
        state.previous_step_times_deriv, state.current_step_times_deriv)
end


text:
module LKH

using Random, LinearAlgebra

# Write your package code here.
include("../deps/deps.jl") # const LKH_EXECUTABLE
include("utils.jl")
include("solver.jl")

export solve_tsp

end


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

#Worked with Ahmed (In Class) and Waleed
function PS5()
# read in function to create state transitions for dynamic model
include("create_grids.jl")

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: reshaping the data
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# load in the data

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)
df
df1 = CSV.read(HTTP.get(url).body,DataFrame)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: estimate a static version of the model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#ul(x1t,b)=00+01x1t
gamma = glm(@formula(Y ~ Branded + Odometer), df_long, Binomial(), LogitLink())
println(gamma)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3a: read in data for dynamic model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)

Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
#[print(string(":"), x, ", ")) for x in names(df)]
set1 = [:Odo1, :Odo2, :Odo3, :Odo4, :Odo5, :Odo6, :Odo7, :Odo8, :Odo9, :Odo10, :Odo11, :Odo12, :Odo13, :Odo14, :Odo15, :Odo16, :Odo17, :Odo18, :Odo19, :Odo20]
Odo = Matrix(df[:,set1])

set2 = [:Xst1, :Xst2, :Xst3, :Xst4, :Xst5, :Xst6, :Xst7, :Xst8, :Xst9, :Xst10, :Xst11, :Xst12, :Xst13, :Xst14, :Xst15, :Xst16, :Xst17, :Xst18, :Xst19, :Xst20]
Xst = Matrix(df[:,set2])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3b: generate state transition matrices
#:::::::::::::::::::::::::::::::::::::::::::::::::::
zval,zbin,xval,xbin,xtran = create_grids()

rowIndex = size(xtran)[1]
outIndx = length(unique(df.Branded))
TT = size(Y)[2]

# future value
FV = zeros(rowIndex, outIndx, TT+1)

            @views @inbounds function myfun(alpha)

                    # Backwards loop from T+1 to 1 over T
                    for t in TT:-1:1
                        #Loop over the two possible brand state {0, 1}
                        for b in 0:1:1
                            #Loop over the possible permanent route usage states (i.e. from 1 to zbin)
                            for z in 1:1:zbin
                                #Loop over the possible permanent route usage states(i.e. from 1 to xbin)
                                for x in 1:1:xbin
                                    row = x + (z-1) * xbin
                                    v1 = gamma[1] + gamma[2]*xval[x] + gamma[3]*b + xtran[row,:]'*FV[(z-1)*xbin+1:z*xbin,b+1,t+1]
                                    v0 = xtran[1+(z-1)*xbin,:]' * FV[(z-1)*xbin+1:z*xbin,b+1,t+1]

                                    FV[row,b+1,t]=0.9*log(exp(v1) + exp(v0))
                                end
                            end
                        end
                    end

#3d Construct the log likelihood using the future value terms from the previous step andonly using the observed states in the data. This will entail a for loop over buses and time periods.
                    Bus = size(Y)[1]
                    lglk = 0
                                for i in Bus
                                    for t in 1:TT
                                        row1 = 1 + (df.Zst[i]-1)*xbin
                            row0 = Xst[i,t] + (df.Zst[i]-1)*xbin

                            util = coef(gamma_hat_glm)[1] + coef(gamma_hat_glm)[2]*Xst[i,t] + coef(gamma_hat_glm)[3]*df.Branded[i] +
                                        (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,df.Branded[i]+1,t+1]
                            P = exp(util)/(1+exp(util))
                            lglk += Y[i,t]log(P)+(Y[i,t]-1)log(1-P)
                                    end
                                end
                                return -lglk
            end
#I had complications running lines 135 and 136.
        optim = optimize(gamma -> myfun(gamma),coef(gamma) , LBFGS(), Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true))
    println(optim.minimizer)

end

PS5()


text:
using FMMLIB2D

using LinearAlgebra
using Random
using Test

Random.seed!(0)

# TODO: Compute gradient and hessian

function direct_self(source, charge, dipvec, dipstr)
    pot = zero(charge)
    N = size(source, 2)
    for i=1:N
        for j=1:N
            if i == j
                continue
            end
            r1 = source[1,j]-source[1,i]
            r2 = source[2,j]-source[2,i]
            rsq = r1^2 + r2^2
            rdotq = r1*dipvec[1,j] + r2*dipvec[2,j]
            pot[i] += charge[j]*log(rsq)/2 + dipstr[j]*rdotq/rsq
        end
    end
    return pot
end
function direct_targ(source, charge, dipvec, dipstr, target)
    M = size(target, 2)
    N = size(source, 2)
    pot = zeros(typeof(charge[1]), M)
    for i=1:M
        for j=1:N
            r1 = source[1,j]-target[1,i]
            r2 = source[2,j]-target[2,i]
            rsq = r1^2 + r2^2
            rdotq = r1*dipvec[1,j] + r2*dipvec[2,j]
            pot[i] += charge[j]*log(rsq)/2 + dipstr[j]*rdotq/rsq
        end
    end
    return pot
end

@testset "Laplace FMM 2D (complex)" begin
    iprec = 5
    N = 40
    M = 30
    source = rand(2, N)
    charge = rand(N) + 1im*rand(N)
    dipstr = rand(N) + 1im*rand(N)
    dipvec = rand(2, N)
    target = rand(2, M)
    ifcharge = 1
    ifdipole = 1
    ifpot = 1
    ifgrad = 1
    ifhess = 1
    ifpottarg = 1
    ifgradtarg = 1
    ifhesstarg = 1
    refpottarg = direct_targ(source, charge, dipvec, dipstr, target)
    refpot = direct_self(source, charge, dipvec, dipstr)
    @testset "Direct interface: lfmm2dpartself" begin
        pot, grad, hess =
            lfmm2dpartself(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess)
        pot_relerr = norm(refpot-pot, Inf) / norm(refpot, Inf)
        @test pot_relerr < 1e-15
    end

    @testset "Direct interface: lfmm2dparttarg" begin
        pot, grad, hess, pottarg, gradtarg, hesstarg =
            lfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
                           target, ifpottarg, ifgradtarg, ifhesstarg)
        pot_relerr = norm(refpot-pot, Inf) / norm(refpot, Inf)
        pottarg_relerr = norm(refpottarg-pottarg, Inf) / norm(refpottarg, Inf)
        @test pot_relerr < 1e-15
        @test pottarg_relerr < 1e-15
    end
    @testset "Keyword interface" begin
        U = lfmm2d(source=source, charge=charge, dipstr=dipstr, dipvec=dipvec)
        @test norm(U.pot-refpot, Inf) / norm(refpot, Inf) < 1e-14
        U = lfmm2d(source=source, charge=charge, dipstr=dipstr, dipvec=dipvec, target=target)
        @test norm(U.pottarg-refpottarg, Inf) / norm(refpot, Inf) < 1e-14
    end
    @testset "Tolerance test" begin
        N = 2000
        source = rand(2, N)
        charge = rand(N) + 1im*rand(N)
        dipstr = rand(N) + 1im*rand(N)
        dipvec = rand(2, N)
        refpot = direct_self(source, charge, dipvec, dipstr)
        for tolexp=-14:-1
            TOL = 10.0^tolexp
            U = lfmm2d(source=source, charge=charge, tol=TOL, dipstr=dipstr, dipvec=dipvec)
            relerr = norm(U.pot-refpot) / norm(refpot)
            @test relerr < TOL
        end
    end
end

@testset "Laplace FMM 2D (real)" begin
    iprec = 5
    N = 40
    M = 30
    source = rand(2, N)
    charge = rand(N)
    dipstr = rand(N)
    dipvec = rand(2, N)
    target = rand(2, M)
    ifcharge = 1
    ifdipole = 1
    ifpot = 1
    ifgrad = 1
    ifhess = 1
    ifpottarg = 1
    ifgradtarg = 1
    ifhesstarg = 1
    refpottarg = direct_targ(source, charge, dipvec, dipstr, target)
    refpot = direct_self(source, charge, dipvec, dipstr)
    @testset "Direct interface: rfmm2dpartself" begin
        pot, grad, hess =
            rfmm2dpartself(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess)
        pot_relerr = norm(refpot-pot, Inf) / norm(refpot, Inf)
        @test pot_relerr < 1e-15
    end
    @testset "Direct interface: rfmm2dparttarg" begin
        pot, grad, hess, pottarg, gradtarg, hesstarg =
            rfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
                           target, ifpottarg, ifgradtarg, ifhesstarg)
        pot_relerr = norm(refpot-pot, Inf) / norm(refpot, Inf)
        pottarg_relerr = norm(refpottarg-pottarg, Inf) / norm(refpottarg, Inf)
        @test pot_relerr < 1e-15
        @test pottarg_relerr < 1e-15
    end
    @testset "Keyword interface" begin
        U = rfmm2d(source=source, charge=charge, dipstr=dipstr, dipvec=dipvec)
        @test norm(U.pot-refpot, Inf) / norm(refpot, Inf) < 1e-14
        U = rfmm2d(source=source, charge=charge, dipstr=dipstr, dipvec=dipvec, target=target)
        @test norm(U.pottarg-refpottarg, Inf) / norm(refpot, Inf) < 1e-14
    end
    @testset "Tolerance test" begin
        N = 2000
        source = rand(2, N)
        charge = rand(N)
        dipstr = rand(N)
        dipvec = rand(2, N)
        refpot = direct_self(source, charge, dipvec, dipstr)
        for tolexp=-14:-1
            TOL = 10.0^tolexp
            U = rfmm2d(source=source, charge=charge, tol=TOL, dipstr=dipstr, dipvec=dipvec)
            relerr = norm(U.pot-refpot) / norm(refpot)
            @test relerr < TOL
        end
    end
end


text:
# Function to create matrix from array of arrays. This version is faster than vcat
function fast_mat(S)
    X = Array{Float32, 2}(undef, length(S[1]), length(S))
    for i=1:length(S)
        X[:,i] = S[i]
    end
    X
end

# This function samples episodes from the policy using importance sampling to weight them
function sample_episodes(mdp, policy::Policy, Neps::Int64; max_steps::Int64 = 1000,
                         a_and_p = action_and_probability)
    S, G, avg_ep_return = [], Float32[], 0.
    for i=1:Neps
        s, steps = initialstate(mdp), 0
        Si, Ri, ρi = [], [], []
        while !isterminal(mdp, s)
            push!(Si, convert_s(Array{Float64,1}, s, mdp))
            a, prob = a_and_p(policy, s)
            push!(ρi, policy.action_probability(mdp, s, a) / prob)
            s, r = gen(mdp, s, a)
            push!(Ri, r)
            steps += 1
            steps >= max_steps && break
        end
        steps >= max_steps && println("Episode timeout at ", max_steps, " steps")
        weighted_returns = reverse(cumsum(reverse(Ri))) .* reverse(cumprod(reverse(ρi)))
        push!(S, Si...)
        push!(G, weighted_returns...)
        avg_ep_return += sum(Ri)
    end
    S, G, avg_ep_return = fast_mat(S), G', avg_ep_return/Neps
    @assert size(S,2) == size(G,2)
    @assert size(G,1) == 1
    S, G, avg_ep_return
end


text:
using NonparEconometricsTool
using Test

@testset "NonparEconometricsTool.jl" begin
    # Write your own tests here.
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    polymat(xs, d)

Return the matrix of monomials for the iterator `xs`, i.e.
for each item `x = (x₁, x₂,…, xₙ)` in `xs`, evaluate
the monomial terms of the expansion `(x₁ + x₂ + ⋯ + xₙ)ᵈ`
for a given degree `d`.

The resulting matrix has a number of rows that is equal
to the number of items in the iterator `xs`. The number
of columns is a function of the degree. For `d=0`, a
single column of ones is returned that corresponds to
the constant term `x₁⁰⋅x₂⁰⋅⋯⋅xₙ⁰` for all items in `xs`.
"""
function polymat(xs, d)
  x, _ = iterate(xs)
  n    = length(x)
  exps = Iterators.flatten(multiexponents(n, d) for d in 0:d)
  m = map(exps) do e
    map(xs) do x
      prod(x .^ e)
    end
  end
  reduce(hcat, m)
end

"""
    trend(data, vars; degree=1)

Return the deterministic spatial trend for the variables `vars`
in the spatial `data`. Approximate the trend with a polynomial
of given `degree`.
"""
function trend(data, vars::AbstractVector{Symbol}; degree=1)
  𝒯 = values(data)
  𝒟 = domain(data)

  # build LHS of linear system
  xs = (coordinates(centroid(𝒟, i)) for i in 1:nelements(𝒟))
  X  = polymat(xs, degree)

  # solve for each variable
  ŷs = map(vars) do v
    y  = Tables.getcolumn(𝒯, v)
    θ  = X'*X \ X'*y
    ŷ  = X*θ
  end

  ctor  = Tables.materializer(𝒯)
  table = ctor((; zip(vars, ŷs)...))

  georef(table, 𝒟)
end

trend(data, var::Symbol; kwargs...) = trend(data, [var]; kwargs...)


text:
function __solve end
function __init end
function solve! end

NO_TSPAN_PROBS = Union{AbstractSteadyStateProblem,AbstractJumpProblem}

function init(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    __init(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    __init(_prob,args...;kwargs...)
  else
    __init(_prob,args...;kwargs...)
  end
end

function solve(prob::DEProblem,args...;kwargs...)
  _prob = get_concrete_problem(prob,kwargs)
  if haskey(kwargs,:alg) && (isempty(args) || args[1] === nothing)
    alg = kwargs[:alg]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    __solve(_prob,alg,args...;kwargs...)
  elseif !isempty(args) && typeof(args[1]) <: DEAlgorithm
    alg = args[1]
    isadaptive(alg) &&
    !(typeof(prob) <: NO_TSPAN_PROBS) &&
    adaptive_warn(_prob.u0,_prob.tspan)
    __solve(_prob,args...;kwargs...)
  else
    __solve(_prob,args...;kwargs...)
  end
end

function get_concrete_problem(prob::AbstractJumpProblem,kwargs)
  prob
end

function get_concrete_problem(prob::AbstractSteadyStateProblem,kwargs)
  if typeof(prob.u0) <: Function
    _u0 = prob.u0(prob.p,Inf)
  else
    _u0 = prob.u0
  end

  __u0 = handle_distribution_u0(_u0)

  remake(prob;u0=__u0)
end

function get_concrete_problem(prob,kwargs)
  if typeof(prob.tspan) <: Function
    _tspan = prob.tspan(prob.p)
  elseif prob.tspan == (nothing,nothing)
    if haskey(kwargs,:tspan)
      _tspan = kwargs.tspan
    else
      error("No tspan is set in the problem or chosen in the init/solve call")
    end
  else
    _tspan = prob.tspan
  end

  if typeof(prob.u0) <: Function
    _u0 = prob.u0(prob.p,_tspan[1])
  else
    _u0 = prob.u0
  end

  __u0 = handle_distribution_u0(_u0)

  remake(prob;u0=__u0,tspan=_tspan)
end

handle_distribution_u0(_u0) = _u0

function adaptive_warn(u0,tspan)
  adaptive_integer_warn(tspan)
end

function adaptive_integer_warn(tspan)
  if eltype(tspan) <: Integer
    @warn("Integer time values are incompatible with adaptive integrators. Utilize floating point numbers instead of integers in this case, i.e. (0.0,1.0) instead of (0,1).")
  end
end


text:
using Plots, RecipesBase, Manifolds, ManifoldsBase, Documenter, PyPlot
# required for loading the manifold tests functios
using Test, ForwardDiff, ReverseDiff, FiniteDifferences
ENV["GKSwstype"] = "100"

generated_path = joinpath(@__DIR__, "src", "misc")
base_url = "https://github.com/JuliaManifolds/Manifolds.jl/blob/master/"
isdir(generated_path) || mkdir(generated_path)

open(joinpath(generated_path, "contributing.md"), "w") do io
    # Point to source license file
    println(
        io,
        """
        ```@meta
        EditURL = "$(base_url)CONTRIBUTING.md"
        ```
        """,
    )
    # Write the contents out below the meta block
    for line in eachline(joinpath(dirname(@__DIR__), "CONTRIBUTING.md"))
        println(io, line)
    end
end

makedocs(
    # for development, we disable prettyurls
    format=Documenter.HTML(prettyurls=false, assets=["assets/favicon.ico"]),
    modules=[Manifolds, ManifoldsBase],
    authors="Seth Axen, Mateusz Baran, Ronny Bergmann, and contributors.",
    sitename="Manifolds.jl",
    pages=[
        "Home" => "index.md",
        "ManifoldsBase.jl" => "interface.md",
        "Examples" => ["How to implement a Manifold" => "examples/manifold.md"],
        "Manifolds" => [
            "Basic manifolds" => [
                "Centered matrices" => "manifolds/centeredmatrices.md",
                "Cholesky space" => "manifolds/choleskyspace.md",
                "Circle" => "manifolds/circle.md",
                "Elliptope" => "manifolds/elliptope.md",
                "Essential manifold" => "manifolds/essentialmanifold.md",
                "Euclidean" => "manifolds/euclidean.md",
                "Fixed-rank matrices" => "manifolds/fixedrankmatrices.md",
                "Generalized Stiefel" => "manifolds/generalizedstiefel.md",
                "Generalized Grassmann" => "manifolds/generalizedgrassmann.md",
                "Grassmann" => "manifolds/grassmann.md",
                "Hyperbolic space" => "manifolds/hyperbolic.md",
                "Lorentzian manifold" => "manifolds/lorentz.md",
                "Multinomial doubly stochastic matrices" => "manifolds/multinomialdoublystochastic.md",
                "Multinomial matrices" => "manifolds/multinomial.md",
                "Multinomial symmetric matrices" => "manifolds/multinomialsymmetric.md",
                "Oblique manifold" => "manifolds/oblique.md",
                "Probability simplex" => "manifolds/probabilitysimplex.md",
                "Positive numbers" => "manifolds/positivenumbers.md",
                "Projective space" => "manifolds/projectivespace.md",
                "Rotations" => "manifolds/rotations.md",
                "Skew-Hermitian matrices" => "manifolds/skewhermitian.md",
                "Spectrahedron" => "manifolds/spectrahedron.md",
                "Sphere" => "manifolds/sphere.md",
                "Stiefel" => "manifolds/stiefel.md",
                "Symmetric matrices" => "manifolds/symmetric.md",
                "Symmetric positive definite" => "manifolds/symmetricpositivedefinite.md",
                "Symmetric positive semidefinite fixed rank" => "manifolds/symmetricpsdfixedrank.md",
                "Torus" => "manifolds/torus.md",
                "Tucker" => "manifolds/tucker.md",
                "Unit-norm symmetric matrices" => "manifolds/spheresymmetricmatrices.md",
            ],
            "Combined manifolds" => [
                "Graph manifold" => "manifolds/graph.md",
                "Power manifold" => "manifolds/power.md",
                "Product manifold" => "manifolds/product.md",
                "Vector bundle" => "manifolds/vector_bundle.md",
            ],
            "Manifold decorators" => [
                "Connection manifold" => "manifolds/connection.md",
                "Group manifold" => "manifolds/group.md",
                "Metric manifold" => "manifolds/metric.md",
            ],
        ],
        "Features on Manifolds" => [
            "Atlases and charts" => "features/atlases.md",
            "Differentiation" => "features/differentiation.md",
            "Distributions" => "features/distributions.md",
            "Statistics" => "features/statistics.md",
            "Testing" => "features/testing.md",
            "Utilities" => "features/utilities.md",
        ],
        "Miscellanea" => [
            "About" => "misc/about.md",
            "Contributing" => "misc/contributing.md",
            "Internals" => "misc/internals.md",
            "Notation" => "misc/notation.md",
        ],
    ],
)
deploydocs(repo="github.com/JuliaManifolds/Manifolds.jl.git", push_preview=true)


text:
# setup project paths -
const _PATH_TO_SRC = dirname(pathof(@__MODULE__))
const _PATH_TO_BASE = joinpath(_PATH_TO_SRC, "base")

# import packages that we will use -
using JSON
using Logging
using DataFrames
using Plots
using BSON
using Distributions
using Random
using StatsBase
using Statistics
using Test
using DataStructures

# set the seed -
Random.seed!(Random.make_seed())

# include our codes -
include(joinpath(_PATH_TO_BASE, "VLTypes.jl"))
include(joinpath(_PATH_TO_BASE, "VLGameEngine.jl"))
include(joinpath(_PATH_TO_BASE, "VLGameFactory.jl"))
include(joinpath(_PATH_TO_BASE, "VLGameAnalysis.jl"))