text:
module VanillaAtmosDiscretizations
using MPI

using ..AtmosDycore
AD = AtmosDycore
using ...Grids
using ...MPIStateArrays

export VanillaAtmosDiscretization

using ...ParametersType
using ...PlanetParameters: cp_d, cv_d, grav, MSLP
using CLIMA.MoistThermodynamics

# ASR Correction to Prandtl number from 7.1 to 0.71
@parameter prandtl 71//100 "Prandtl number: ratio of momentum diffusivity to thermal diffusivity"
@parameter λ_stokes  -2//3  "scaling for viscous effect associated with volume change"
@parameter k_μ cp_d/prandtl "thermal conductivity / dynamic viscosity"

const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E)
const _nstategrad = 15
const (_ρx, _ρy, _ρz, _ux, _uy, _uz, _vx, _vy, _vz, _wx, _wy, _wz,
       _Tx, _Ty, _Tz) = 1:_nstategrad

"""
    VanillaAtmosDiscretization{nmoist, ntrace}(grid; gravity = true,
    viscosity = 0)

Given a 'grid <: AbstractGrid' this construct all the data necessary to run a
vanilla discontinuous Galerkin discretization of the the compressible Euler
equations with `nmoist` moisture variables and `ntrace` tracer variables. If the
boolean keyword argument `gravity` is `true` then gravity is used otherwise it
is not. Isotropic viscosity can be used if `viscosity` is set to a positive
constant.
"""
struct VanillaAtmosDiscretization{T, dim, polynomialorder, numberofDOFs,
                                  DeviceArray, nmoist, ntrace, DASAT3,
                                  GT } <: AD.AbstractAtmosDiscretization
  "grid"
  grid::GT

  "gravitational acceleration (m/s^2)"
  gravity::T

  "viscosity constant"
  viscosity::T

  "storage for the grad"
  grad::DASAT3

  VanillaAtmosDiscretization(grid;
                             # How many tracer variables
                             ntrace=0,
                             # How many moisture variables
                             nmoist=0,
                             kw...) =
  VanillaAtmosDiscretization{nmoist, ntrace}(grid; kw...)

  function VanillaAtmosDiscretization{nmoist, ntrace
                                     }(grid::AbstractGrid{T, dim, N, Np, DA};
                                       # Use gravity?
                                       gravity = true,
                                       # viscosity constant
                                       viscosity = 0
                                      ) where {T, dim, N, Np, DA,
                                               nmoist, ntrace}
    topology = grid.topology

    ngrad = _nstategrad + 3*nmoist
    # FIXME: Remove after updating CUDA
    h_vgeo = Array(grid.vgeo)
    grad = MPIStateArray{Tuple{Np, ngrad}, T, DA}(topology.mpicomm,
                                                  length(topology.elems),
                                                  realelems=topology.realelems,
                                                  ghostelems=topology.ghostelems,
                                                  sendelems=topology.sendelems,
                                                  nabrtorank=topology.nabrtorank,
                                                  nabrtorecv=topology.nabrtorecv,
                                                  nabrtosend=topology.nabrtosend,
                                                  weights=view(h_vgeo, :, grid.Mid, :))

    GT = typeof(grid)
    DASAT3 = typeof(grad)

    new{T, dim, N, Np, DA, nmoist, ntrace, DASAT3, GT}(grid,
                                                       gravity ? grav : 0,
                                                       viscosity, grad)
  end
end

function Base.getproperty(X::VanillaAtmosDiscretization{T, dim, polynomialorder,
                                                        numberofDOFs,
                                                        DeviceArray, nmoist,
                                                        ntrace},
                          s::Symbol) where {T, dim, polynomialorder,
                                            numberofDOFs, DeviceArray, nmoist,
                                            ntrace}
  if s ∈ keys(stateid)
    stateid[s]
  elseif s == :nstate
    _nstate
  elseif s == :moistrange
    _nstate .+ (1:nmoist)
  elseif s == :tracerange
    _nstate+nmoist .+ (1:ntrace)
  else
    getfield(X, s)
  end
end

function Base.propertynames(X::VanillaAtmosDiscretization)
  (fieldnames(VanillaAtmosDiscretization)..., keys(stateid)...)
end

"""
    MPIStateArray(disc::VanillaAtmosDiscretization)

Given a discretization `disc` constructs an `MPIStateArrays` for holding a
solution state
"""
function MPIStateArrays.MPIStateArray(disc::VanillaAtmosDiscretization{
                                            T, dim, N, Np, DA, nmoist, ntrace}
                                      ) where {T, dim, N, Np, DA, nmoist,
                                               ntrace}
  topology = disc.grid.topology
  nvar = _nstate + nmoist + ntrace
  # FIXME: Remove after updating CUDA
  h_vgeo = Array(disc.grid.vgeo)
  MPIStateArray{Tuple{Np, nvar}, T, DA}(topology.mpicomm,
                                        length(topology.elems),
                                        realelems=topology.realelems,
                                        ghostelems=topology.ghostelems,
                                        sendelems=topology.sendelems,
                                        nabrtorank=topology.nabrtorank,
                                        nabrtorecv=topology.nabrtorecv,
                                        nabrtosend=topology.nabrtosend,
                                        weights=view(h_vgeo, :, disc.grid.Mid, :))
end

function MPIStateArrays.MPIStateArray(disc::VanillaAtmosDiscretization{
                                               T, dim, N, Np, DA, nmoist,
                                               ntrace}, ic::Function
                                       ) where {T, dim, N, Np, DA, nmoist,
                                                ntrace}
  Q = MPIStateArray(disc)

  nvar = _nstate + nmoist + ntrace
  G = disc.grid
  vgeo = G.vgeo

  # FIXME: GPUify me
  host_array = Array ∈ typeof(Q).parameters
  (h_vgeo, h_Q) = host_array ? (vgeo, Q) : (Array(vgeo), Array(Q))
  @inbounds for e = 1:size(Q, 3), i = 1:Np
    x, y, z = h_vgeo[i, G.xid, e], h_vgeo[i, G.yid, e], h_vgeo[i, G.zid, e]
    q0 = ic(x, y, z)

    # Assume that this will be compile time?

    @assert ((nmoist >  0 && nmoist == length(q0.Qmoist)) ||
             (nmoist == 0 && :Qmoist ∉ fieldnames(typeof(q0))))
    @assert ((ntrace >  0 && ntrace == length(q0.Qtrace)) ||
             (ntrace == 0 && :Qtrace ∉ fieldnames(typeof(q0))))

    h_Q[i, [_ρ, _U, _V, _W, _E], e] .= (q0.ρ, q0.U, q0.V, q0.W, q0.E)
    (nmoist > 0) && (h_Q[i, _nstate .+           (1:nmoist), e] .= q0.Qmoist)
    (ntrace > 0) && (h_Q[i, _nstate .+ nmoist .+ (1:ntrace), e] .= q0.Qtrace)
  end
  if !host_array
    Q .= h_Q
  end

  Q
end

MPIStateArrays.MPIStateArray(f::Function,
                             d::VanillaAtmosDiscretization
                            ) = MPIStateArray(d, f)

"""
    estimatedt(disc::VanillaAtmosDiscretization, Q::MPIStateArray)

Given a discretization `disc` and a state `Q` compute an estimate for the time
step

!!! todo

    This estimate is currently very conservative, needs to be revisited
"""
function estimatedt(disc::VanillaAtmosDiscretization{T, dim, N, Np, DA, nmoist},
                    Q::MPIStateArray) where {T, dim, N, Np, DA, nmoist}

  @assert T == eltype(Q)
  G = disc.grid
  vgeo = G.vgeo
  # FIXME: GPUify me
  host_array = Array ∈ typeof(Q).parameters
  (h_vgeo, h_Q) = host_array ? (vgeo, Q) : (Array(vgeo), Array(Q))
  estimatedt(Val(dim), Val(N), Val(nmoist), G, disc.gravity, h_Q, h_vgeo, G.topology.mpicomm)
end

# FIXME: This needs cleaning up
function estimatedt(::Val{dim}, ::Val{N}, ::Val{nmoist}, G, gravity, Q, vgeo,
                    mpicomm) where {dim, N, nmoist}

  DFloat = eltype(Q)

  Np = (N+1)^dim
  (_, _, nelem) = size(Q)

  dt = [floatmax(DFloat)]

  # Allocate 3 spaces for moist tracers qm, with a zero default value
  q_m = zeros(DFloat, max(3, nmoist))

  if dim == 2
    @inbounds for e = 1:nelem, n = 1:Np
      ρ, U, V = Q[n, _ρ, e], Q[n, _U, e], Q[n, _V, e]
      E = Q[n, _E, e]
      y = vgeo[n, G.yid, e]

      #compute temperature and internal energy
      #get moist variables from state vector
      for m = 1:nmoist
          s = _nstate + m
          q_m[m] = Q[n, s, e]/ρ
      end

      E_int = E - (U^2 + V^2)/(2*ρ) - ρ * gravity * y
      # get adjusted temperature and liquid and ice specific humidities
      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])

      ξx, ξy, ηx, ηy = vgeo[n, G.ξxid, e], vgeo[n, G.ξyid, e],
                       vgeo[n, G.ηxid, e], vgeo[n, G.ηyid, e]

       # Calculate local dt
      loc_dt = 2ρ / max(abs(U * ξx + V * ξy) + ρ * soundspeed_air(T),
                        abs(U * ηx + V * ηy) + ρ * soundspeed_air(T))
      dt[1] = min(dt[1], loc_dt)

    end
  end

  if dim == 3
    @inbounds for e = 1:nelem, n = 1:Np
      ρ, U, V, W = Q[n, _ρ, e], Q[n, _U, e], Q[n, _V, e], Q[n, _W, e]
      E = Q[n, _E, e]
      z = vgeo[n, G.zid, e]

      #Compute (Temperature) and (E_int per unit mass)
      E_int = E - (U^2 + V^2+ W^2)/(2*ρ) - ρ * gravity * z
      #Loop over moist variables and extract q_m where q_m[1] = q_t, q_m[2] = q_liq, q_m[3] = q_ice
      for m = 1:nmoist
          s = _nstate + m
          q_m[m] = Q[n, s, e] / ρ
      end

      # get adjusted temperature and liquid and ice specific humidities
      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])

      ξx, ξy, ξz = vgeo[n, G.ξxid, e], vgeo[n, G.ξyid, e], vgeo[n, G.ξzid, e]
      ηx, ηy, ηz = vgeo[n, G.ηxid, e], vgeo[n, G.ηyid, e], vgeo[n, G.ηzid, e]
      ζx, ζy, ζz = vgeo[n, G.ζxid, e], vgeo[n, G.ζyid, e], vgeo[n, G.ζzid, e]

      loc_dt = 2ρ / max(abs(U * ξx + V * ξy + W * ξz) + ρ * soundspeed_air(T),
                        abs(U * ηx + V * ηy + W * ηz) + ρ * soundspeed_air(T),
                        abs(U * ζx + V * ζy + W * ζz) + ρ * soundspeed_air(T))
      dt[1] = min(dt[1], loc_dt)

  end
  end

  MPI.Allreduce(dt[1], MPI.MIN, mpicomm) / N^√2
end

#{{{ rhs!

function AD.getrhsfunction(disc::VanillaAtmosDiscretization)
  (x...) -> rhs!(x..., disc)
end

rhs!(dQ, Q, p::Nothing, t, sd::VanillaAtmosDiscretization) = rhs!(dQ, Q, t, sd)

function rhs!(dQ::MPIStateArray{S, T}, Q::MPIStateArray{S, T}, t::T,
              disc::VanillaAtmosDiscretization{T, dim, N, Np, DA, nmoist,
                                               ntrace}
             ) where {S, T, dim, N, Np, DA, nmoist, ntrace}
  grid = disc.grid
  topology = grid.topology

  gravity = disc.gravity

  grad = disc.grad

  vgeo = grid.vgeo
  sgeo = grid.sgeo
  Dmat = grid.D
  vmapM = grid.vmapM
  vmapP = grid.vmapP
  elemtobndy = grid.elemtobndy

  DFloat = eltype(Q)
  @assert DFloat == eltype(Q)
  @assert DFloat == eltype(vgeo)
  @assert DFloat == eltype(sgeo)
  @assert DFloat == eltype(Dmat)
  @assert DFloat == eltype(grad)
  @assert DFloat == eltype(dQ)

  ########################
  # Gradient Computation #
  ########################
  MPIStateArrays.start_ghost_exchange!(Q)

  volumegrad!(Val(dim), Val(N), Val(nmoist), Val(ntrace), grad.Q, Q.Q, vgeo,
              gravity, Dmat, topology.realelems)

  MPIStateArrays.finish_ghost_exchange!(Q)

  facegrad!(Val(dim), Val(N), Val(nmoist), Val(ntrace), grad.Q, Q.Q, vgeo,
            sgeo, gravity, topology.realelems, vmapM, vmapP, elemtobndy)

  ###################
  # RHS Computation #
  ###################

  viscosity::DFloat = disc.viscosity

  MPIStateArrays.start_ghost_exchange!(grad)

  volumerhs!(Val(dim), Val(N), Val(nmoist), Val(ntrace), dQ.Q, Q.Q, grad.Q,
             vgeo, gravity, viscosity, Dmat, topology.realelems)

  MPIStateArrays.finish_ghost_exchange!(grad)

  facerhs!(Val(dim), Val(N), Val(nmoist), Val(ntrace), dQ.Q, Q.Q, grad.Q,
           vgeo, sgeo, gravity, viscosity, topology.realelems, vmapM, vmapP,
           elemtobndy)
end
# }}}

# {{{ FIXME: remove this after we've figure out how to pass through to kernel
const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _MJ, _MJI,
       _x, _y, _z = 1:_nvgeo

const _nsgeo = 5
const _nx, _ny, _nz, _sMJ, _vMJI = 1:_nsgeo
# }}}

using Requires

include("VanillaAtmosDiscretizations_kernels.jl")

end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
include("solvers/annealing.jl")

function main_annealing(args)
    ln1("="^70)
    ln1("Début de l'action annealing")
    inst = Instance(args[:infile])

    # Construction de la solution initiale :
    sol = Solution(inst)
    ln1("Solution correspondant à l'ordre de l'instance")
    ln1(to_s(sol))

    initial_sort!(sol)
    ln1("Solution initiale envoyée au solver")
    ln1(to_s(sol))

    sv = AnnealingSolver(
        inst;
        temp_init_rate=0.8,
        step_size=1,
        startsol=sol,
        temp_coef=0.985,
    )
    ln1(get_stats(sv))

    ms_start = ms()
    solve(sv, swap_close_planes!, durationmax=15*60)
    ms_stop = ms()

    bestsol = sv.bestsol
    print_sol(bestsol)

    nb_calls = bestsol.solver.nb_calls
    nb_infeasable = bestsol.solver.nb_infeasable
    nb_sec = round(ms_stop - ms_start, digits=3)
    nb_call_per_sec = round(nb_calls/nb_sec, digits=3)
    println("Performance: ")
    println("  nb_calls=$nb_calls")
    println("  nb_infeasable=$nb_infeasable")
    println("  nb_sec=$nb_sec")
    println("  => nb_call_per_sec = $nb_call_per_sec call/sec")

    ln1("Fin de l'action annealing")
end

main_annealing(Args.args)

ext: jl


text:
module Julia2Tikz

include("./TikzTypes.jl")
include("./TikzAxis.jl")
include("./TikzFigure.jl")
include("./PlotStyles/PlotStyles.jl")

# Mapping between the plot_type argument and underlying function when calling plot()
const plot_type_map = Dict{Symbol, Function}(
										     :scatter => scatter
									        )

# Plotting routines
function plot(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, attributes::Dict{String,String}; plotType::Symbol=:scatter)
	plt = TikzFigure()
	defaultAxis = TikzAxis()
	push!(plt.axes, defaultAxis)
	plot!(plt, x, y, attributes; plotType)
	return plt
end

function plot!(figure::AbstractTikzFigure, x::AbstractVector{<: Real}, y::AbstractVector{<: Real}, attributes::Dict{String,String}; plotType::Symbol=:scatter)
    plot!(figure.axes[end], x, y, attributes; plotType)
end

function plot!(axis::AbstractTikzAxis, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, attributes::Dict{String,String}; plotType::Symbol=:scatter)
	plotType in keys(plot_type_map) || ArgumentError("Invalid plot type $plotType")
	thisPlot = plot_type_map[plotType](x, y, attributes)
	push!(axis.plots, thisPlot)
end

function saveTikz(figure::TikzFigure, fileName::String)
	open(fileName, "w") do f
		write(f, "\\begin{tikzpicture}\n")
		for ax in figure.axes
			write(f, "\\begin{axis}[\n]\n")
			for line in ax.plots
				write(f, ToTikzString(line))
			end
			write(f, "\\end{axis}\n")
		end
		write(f, "\\end{tikzpicture}\n")
	end
end

end


text:
using Optim
using Distributions

function loglik(x, μ, logσ)
    nd = Normal(μ, exp(logσ))
    -sum(logpdf(nd, v) for v in x)
end

function testoptim(x)
    res = optimize(par -> loglik(x, par[1], par[2]), zeros(2))
    display(res)
    res.minimizer[1], exp(res.minimizer[2])
end

max_stars_repo_path: Chapter09/03. Estimation using maximum likelihood/opt.jl


text:
"""
Data and priors for estimating R using Kalman filter.
"""
struct RtModel{C,P,V,T, T0, IT,TT}
  "Δlog(Δcases)"
  dlogk::C
  "Covariates shifting mean of R[t]"
  X::T
  "Covariates shifting mean of R[0]"
  X0::T0
  priors::NamedTuple{P,V}
  MA1::Bool
  "Identifier"
  id::IT
  t::TT
end

"""
    RtModel(d::Array{Vector{T}}, p::NamedTuple{P,V})

Create RtModel with `d = Δlog(Δcases)`, prior `p`, and no covariates.
"""
function RtModel(d::Array{Vector{T}}, p::NamedTuple{P,V}) where {T, P , V}
  RtModel(d, [ones(length(ds),1) for ds in d],
          fill(ones(1), length(d)), p, true, 1:length(d),
          [1:length(di) for di in d])
end

RtModel(d, x, x0, p) = RtModel(d,x,x0,p,false, 1:length(d), [1:length(di) for di in d])

"""
    RtModel(df::AbstractDataFrame, casevar::Symbol,
            xvars::Vector{Symbol}, x0vars::Vector{Symbol},
            priors::NamedTuple; L1=1, L2=1,
            time0=r->(r[casevar].>=5), t=:date, id=:fips)

Construct `RtModel` from DataFrame `df.

# Arguments
- `df`
- `casevar` variable in `df` that contains cumulative cases
- `xvars` variables in `df` shifting mean of `R[t]`
- `x0vars` variables in `df` shifting mean of `R[0]`
- `priors`
- `L1=1` lags to use to calculate `Δcases`
- `L2=1` lags to use calculage `Δlog(Δcases)`
- `time0=r->(r[casevar].>=5)` criteria for model time 0. The defaults sets t=0 when cases are first 5 or more.
- `t` time variable in `df`
- `id` id varaible in `df`
- `minvalforlog` instead of `log(Δcases)`, we take `log(min(minvalforlog, Δcases))` to avoid log(0)

"""
function RtModel(df::DataFrames.AbstractDataFrame, casevar::Symbol,
                 xvars::Vector{Symbol}, x0vars::Vector{Symbol},
                 priors::NamedTuple; L1=1, L2=1,
                 time0=r->(r[casevar].>=5), t=:date, id=:state,
                 minvalforlog=0.1)

  sdf = copy(df[!,unique(vcat(id, t, casevar, xvars..., x0vars...))])
  sort!(sdf, (id, t))

  sdf[!,:Δcases] = by(sdf, :state, Δcases = casevar => x->lagdiff(x, L1))[!,:Δcases]./L1
  sdf[!,:ΔlogΔcases] = by(sdf, :state, ΔlogΔcases =
                          :Δcases => x->lagdiff(log.(max.(x,minvalforlog)), L2))[!,:ΔlogΔcases]./L2

  sort!(sdf, (id, t))
  ids = unique(sdf[!,id])
  dlogk = Vector{Vector{Float64}}(undef, length(ids))
  dates = Vector{Vector{eltype(sdf[!,t])}}(undef, length(ids))
  X = Vector{Matrix{Float64}}(undef, length(ids))
  X0 = Vector{Vector{Float64}}(undef, length(ids))
  for (i, d) in enumerate(ids)
    gdf = filter(x->x[id].==d, sdf)
    start = findfirst(time0(gdf) .& .!(ismissing.(gdf.ΔlogΔcases)))
    dlogk[i] = Vector{Float64}(gdf[start:end, :ΔlogΔcases])
    X[i] = Matrix{Float64}(gdf[start:end, xvars])
    X0[i] = Vector{Float64}(gdf[start, x0vars])
    dates[i] = Vector{eltype(sdf[!,t])}(gdf[start:end, t])
  end

  return(RtModel(dlogk, X, X0, priors, false, ids, dates))
end

"""
    function (m::RtModel)(param)

Loglikelihood for RtModel
"""
function (m::RtModel)(param)
  Parameters.@unpack dlogk, X, X0, priors = m
  Parameters.@unpack σR, σk, σR0, γ, ρ, α, α0 = param
  logp = try
    logp = logpdf(priors.σR, σR) +
      logpdf(priors.γ, γ) +
      logpdf(priors.σR0, σR0) +
      #logpdf(priors.μR0, μR0) +
      logpdf(priors.σk, σk) +
      logpdf(priors.α, α) +
      logpdf(priors.α0, α0) +
      logpdf(priors.ρ,ρ)

    # Note that this not the usual Kalman filter because the noise in
    # dlogk is an MA(1) instead of independent. We use
    # https://arxiv.org/pdf/1909.10582.pdf
    # and adopt their notation
    A = m.MA1 ? SMatrix{2,2}([2*σk^2  -σk^2; -σk^2  2*σk^2]) : SMatrix{2,2}([2*σk^2 0; 0 2*σk^2])
    S = length(dlogk)
    T = maximum(length.(dlogk))
    meanR = Vector{typeof(σR)}(undef, T)
    zhat = Vector{typeof(σR)}(undef, T)
    varR, varRprior, K, varz, zcoef = ma1kalman_variance(T, σR0^2, ρ, σR^2, γ, A)
    for s in 1:S
      z = dlogk[s] .+ γ .- γ* (X[s]*α)
      μR0 = dot(X0[s],α0)
      ma1kalman_mean!(meanR, zhat, z, ρ, γ, μR0, K, zcoef) #, X[s]*α)
      for t in 1:length(dlogk[s])
        μ = zhat[t]
        σ = sqrt(varz[t])
        logp += logpdf(Normal(μ, σ), z[t])
      end
    end
    logp
  catch err
    @warn "Error while evaluating log likelihood. Returning -Inf.\n"*
    "The error message was:\n"*
    err.msg

    logp = -Inf*σR
  end
  return(logp)
end

"""
     mcmc(mdl::RtModel;
              θ0=(γ = 1/7, σR0=1.0, α0 = zeros(length(mdl.X0[1])),
                  σR=0.25, σk=2.0, ρ=1.0, α=zeros(size(mdl.X[1],2))),
              iterations = 2000,
              warmup=default_warmup_stages(local_optimization=nothing,
                                           stepsize_search=nothing,
                                           init_steps=100, middle_steps=100,
                                           terminating_steps=2*100,
                                           doubling_stages=3, M=Symmetric))

Estimate `RtModel` using `DynamicHMC`
"""
function mcmc(mdl::RtModel;
              θ0=(γ = 1/7, σR0=1.0, α0 = zeros(length(mdl.X0[1])),
                  σR=0.25, σk=2.0, ρ=1.0, α=zeros(size(mdl.X[1],2))),
              iterations = 2000,
              warmup=default_warmup_stages(local_optimization=nothing,
                                           stepsize_search=nothing,
                                           init_steps=100, middle_steps=100,
                                           terminating_steps=2*100,
                                           doubling_stages=3, M=Symmetric))

  K = size(mdl.X[1], 2)
  rlo = minimum(mdl.priors.ρ)
  rhi = maximum(mdl.priors.ρ)
  trans = as( (γ = asℝ₊, σR0 = asℝ₊, α0 = as(Array, length(mdl.X0[1])),
               σR = asℝ₊, σk = asℝ₊, ρ=as(Real, rlo, rhi),
               α = as(Array, K)) )
  P = TransformedLogDensity(trans, mdl)
  ∇P = ADgradient(:ForwardDiff, P)
  x0 = inverse(trans,θ0)

  rng = Random.MersenneTwister()
  res = DynamicHMC.mcmc_keep_warmup(rng, ∇P, 2000;initialization = (q = x0, ϵ=0.1),
                                    reporter = LogProgressReport(nothing, 25, 15),
                                    warmup_stages =warmup);
  post = transform.(trans,res.inference.chain)
  return(post)

end

"""
      MCMCChain(post::Array{NamedTuple,1}, xvars::Vector{Symbol}, x0vars::Vector{Symbol})

   Convert posterior from RtModel.mcmc to an MCMCChain.Chains
"""
function MCMCChain(post::Array{NamedTuple{N,V},1}, xvars::Vector{Symbol}, x0vars::Vector{Symbol}) where {N, V}
  p = post[1]
  vals = hcat([vcat([length(v)==1 ? v : vec(v) for v in values(p)]...) for p in post]...)'
  vals = reshape(vals, size(vals)..., 1)
  names = vcat([length(p[s])==1 ? String(s) : String.(s).*"[".*string.(1:length(p[s])).*"]"
                for s in keys(p)]...)
  for (i,m) in enumerate(match.(r"α0\[(\d+)\]", names))
    if (m !=nothing)
      j = parse(Int, m.captures[1])
      names[i] = "α0[$(x0vars[j])]"
    end
  end
  for (i,m) in enumerate(match.(r"α\[(\d+)\]", names))
    if (m !=nothing)
      j = parse(Int, m.captures[1])
      names[i] = "α[$(xvars[j])]"
    end
  end
  names .= replace.(names, r"\[" => s"(")
  names .= replace.(names, r"\]" => s")")
  names .= replace.(names, r"_" => s" ")
  names .= replace.(names, r"\." => s" ")
  return(MCMCChains.Chains(vals, names))
end

"""
    function plotpostr(dates, dlogk, post, X, X0; Δ=1)

Plots estimated posterior of `R[t]` on `dates` given posterior draws
of the parameters `post`, and covariates `X` and `X0`.

`post`: should be an array of named tuples of parameters, as returned
by `DynamicHMC.mcmc_with_warmup`
"""
function plotpostr(dates, dlogk, post, X, X0; Δ=1)
  k = [kalman(dlogk, p.σR, p.σk, p.σR0, dot(X0,p.α0), p.γ, p.ρ, X*p.α) for p in post];
  γ = [p.γ for p in post];
  Xa = hcat([X*p.α for p in post]...)
  meanR = hcat([x[1] for x in k]...)./Δ;
  varR = hcat([x[2] for x in k]...)./Δ^2;
  zhat = hcat([x[3] for x in k]...);
  c = "black"
  figr = plot(dates, mean(meanR, dims=2), ribbon=1.64*mean(sqrt.(varR),dims=2),
              color=c, fillalpha=0.2, linewidth=1.5, legend=:none, ylab="Rₜ")
  r=([quantile(meanR[t,:] - 1.64*sqrt.(varR[t,:]), 0.05) for t in 1:size(meanR,1)],
     [quantile(meanR[t,:] + 1.64*sqrt.(varR[t,:]), 0.95) for t in 1:size(meanR,1)])
  figr = plot!(figr, dates, zeros(length(r[1])), ribbon=(-r[1], r[2]), color=c,
               linewidth=0, ylim=nothing, fillalpha=0.2)
  T = length(dlogk)
  rt = [mean(dlogk[t]./(Δ*γ) .+ 1) for t in 1:T]
  lo = [quantile(dlogk[t]./(Δ*γ) .+ 1, 0.05) for t in 1:T]
  hi = [quantile(dlogk[t]./(Δ*γ) .+ 1, 0.95) for t in 1:T]
  figr = scatter!(figr, dates, rt, yerror=(rt-lo, hi-rt))
  figr = hline!(figr, [1.], color="red", linewidth=1.5, linestyle=:dash)
  return(figr)
end


text:
abstract type AbstractParameter end

"""
    value(x)

Return the "value" of an object.
For `AbstractParameter`s this typically applies some transformation to some data
contained in the parameter, and returns a plain data type.
It might, for example, return a transformation of some internal data, the result of which
is guaranteed to satisfy some constraint.
"""
value(x)

# Various basic `value` definitions.
value(x::Number) = x
value(x::AbstractArray{<:Number}) = x
value(x::AbstractArray) = map(value, x)
value(x::Tuple) = map(value, x)
value(x::NamedTuple) = map(value, x)
value(x::Dict) = Dict(k => value(v) for (k, v) in x)


text:
srand(1)

context("neighbors.jl") do
  landscapes = [
    NK.NKLandscape(4, 1),
    NK.NKqLandscape(4, 1, 2),
    NK.NKpLandscape(4, 1, 0.90),
  ]

  for ls = landscapes
    # The ith locus is linked to the locus i places to its right, assuming a
    # periodic boundary.
    ls.links = [
      0b1000,
      0b0100,
      0b0010,
      0b0001
    ]

    g = NK.Genotype(0b0000, ls)
    actual = [0b0001, 0b0010, 0b0100, 0b1000]

    context("$(ls)") do
      context("number_neighbors(...)") do
        @fact NK.number_neighbors(g) --> 4
      end

      context("all_neighbors(...)") do
        nbrs = NK.all_neighbors(g)
        @fact length(nbrs) --> 4
        @fact nbrs[1].alleles --> anyof(actual...)
        @fact nbrs[1].landscape --> ls
        @fact nbrs[2].alleles --> anyof(actual...)
        @fact nbrs[2].landscape --> ls
        @fact nbrs[3].alleles --> anyof(actual...)
        @fact nbrs[3].landscape --> ls
        @fact nbrs[4].alleles --> anyof(actual...)
        @fact nbrs[4].landscape --> ls
      end

      context("random_neighbor(...)") do
        rnbr = NK.random_neighbor(g)
        @fact rnbr.alleles --> anyof(actual...)
        @fact rnbr.landscape --> ls
      end

      context("neutral_neighbors(...)") do
        # TODO: Fake out `contribs` to guarantee neighbors.
        nbrs = NK.neutral_neighbors(g)
        for nbr = nbrs
          @fact NK.fitness(nbr) --> NK.fitness(g)
          @fact nbr.alleles --> anyof(actual...)
          @fact nbr.landscape --> ls
        end
      end

      context("fitter_neighbors(...)") do
        # TODO: Fake out `contribs` to guarantee neighbors.
        nbrs = NK.fitter_neighbors(g)
        for nbr = nbrs
          @fact NK.fitness(nbr) --> greater_than(fitness(g))
          @fact nbr.alleles --> anyof(actual...)
          @fact nbr.landscape --> ls
        end
      end

      function test_fitness_range_neighbors(lb::Float64, ub::Float64)
        nbrs = NK.fitness_range_neighbors(g, lb, ub)
        if lb == 0.0 && ub > 1.0
          @fact length(nbrs) --> 4
        end
        for nbr = nbrs
          @fact NK.fitness(nbr) --> greater_than_or_equal(lb)
          @fact NK.fitness(nbr) --> less_than(ub)
          @fact nbr.alleles --> anyof(actual...)
          @fact nbr.landscape --> ls
        end
      end

      context("fitness_range_neighbors(...)") do
        test_fitness_range_neighbors(0.0, 1.0 + eps(Float64))
        test_fitness_range_neighbors(0.25, 0.75)
      end

      context("fittest_neighbors(...)") do
        nbrs = NK.fittest_neighbors(g, 2)
        allfits = [NK.fitness(x) for x = NK.all_neighbors(g)] |> sort
        @fact length(nbrs) --> 2
        for nbr = nbrs
          @fact NK.fitness(nbr) --> anyof(allfits[(end - 1):end]...)
          @fact nbr.alleles --> anyof(actual...)
          @fact nbr.landscape --> ls
        end
      end

      context("fittest_neighbor(...)") do
        fnbr = NK.fittest_neighbor(g)
        allfits = [NK.fitness(x) for x = NK.all_neighbors(g)]
        @fact NK.fitness(fnbr) --> maximum(allfits)
        @fact fnbr.alleles --> anyof(actual...)
        @fact fnbr.landscape --> ls
      end
    end
  end
end


text:
# -------- Auxiliary functions ----------------------------------------------- #

dict(t :: Symbol; kwargs...) = Dict{Symbol, Any}(:type => t, kwargs...)
compose(d :: Dict{Symbol, Any}) = compose(Val(d[:type]), d)

# -------- Basic Conversions ------------------------------------------------- #

decompose(p :: Bool)    = p
decompose(p :: Int)     = p
decompose(p :: Float32) = p
decompose(p :: Float64) = p
decompose(p :: String)  = p
decompose(p :: Symbol)          = p

compose(p :: Bool)    = p
compose(p :: Int)     = p
compose(p :: Float32) = p
compose(p :: Float64) = p
compose(p :: String)  = p
compose(p :: Symbol)  = p

decompose(p :: Array)           = dict(:array, value = p)
decompose(p :: Tuple{Int, Int}) = dict(:pair, a = p[1], b = p[2])
decompose(p :: Nothing)         = dict(:nothing)

compose(:: Val{:array}, d)   = d[:value]
compose(:: Val{:pair}, d)    = (d[:a], d[:b])
compose(:: Val{:nothing}, d) = nothing

# -------- Function Conversion ----------------------------------------------- #

decompose(:: typeof(identity)) = dict(:function, name = :identity)
decompose(:: typeof(softmax))  = dict(:function, name = :softmax)
decompose(:: typeof(relu))     = dict(:function, name = :relu)
decompose(:: typeof(tanh))     = dict(:function, name = :tanh)
decompose(:: typeof(sigm))     = dict(:function, name = :sigm)
decompose(:: typeof(elu))      = dict(:function, name = :elu)
decompose(:: typeof(zeros))    = dict(:function, name = :zeros)
decompose(:: typeof(ones))     = dict(:function, name = :ones)

function compose(:: Val{:function}, d)
  if d[:name] in [ :identity, :softmax, :relu
                 , :tanh, :sigm, :elu, :zeros, :ones ]
    eval(d[:name])
  else
    error("Cannot compose function $(d[:name])")
  end
end

# -------- Conversion Macros ------------------------------------------------- #

macro decompose(t, instance)
  names = :(Base.fieldnames(typeof($instance)))
  pairs = :([n => decompose(Base.getfield($instance, n)) for n in $names])
  esc(:(dict($t; $pairs...)))
end

macro compose(t, dict)
  names = :(Base.fieldnames($t))
  args = :([compose($dict[entry]) for entry in $names])
  esc(:($t($args...)))
end

# -------- Knet Conversions -------------------------------------------------- #

decompose(p :: Knet.Param)      = dict(:parameter, value = Knet.value(p))
decompose(bm :: Knet.BNMoments) = @decompose :bnmoments bm

compose(:: Val{:parameter}, d) = Knet.Param(d[:value])
compose(:: Val{:bnmoments}, d) = @compose Knet.BNMoments d

# -------- Layer Conversions ------------------------------------------------- #

decompose(v :: Vector{Layer})   = dict(:layers, value = decompose.(v))
compose(:: Val{:layers}, d)   = Layer[compose(l) for l in d[:value]]

# -------- Feature Conversions ----------------------------------------------- #

decompose(f :: ConstantFeature) = @decompose :constantfeature f
compose(:: Val{:constantfeature}, d) = @compose ConstantFeature d

decompose(f :: Vector{Feature}) = dict(:features, features = decompose.(f))
compose(:: Val{:features}, d) = compose.(d[:features])

# -------- Game Type Conversions --------------------------------------------- #

function decompose(:: Type{G}) where {G <: Game}
  name = nameof(G)
  params = G.parameters
  dict(:gametype, name = nameof(G), params = collect(G.parameters))
end

function compose(:: Val{:gametype}, d)
  name   = compose(d[:name])
  params = compose.(d[:params]) |> collect

  @assert isa(name, Symbol) "Cannot compose gametype"
  @assert all(isbits, params) "Cannot compose gametype"

  if isempty(d[:params])
    eval(name)
  else
    eval(Expr(:curly, name, params...))
  end
end

# TODO: The :gametype and :features decompositions are not transparent!
# But maybe this is not a problem, as they are not needed for remote
# reconstruction of playing ability?

# -------- Layers ------------------------------------------------------------ #

decompose(l :: Pointwise{false}) = @decompose :pointwise l
decompose(l :: Dense{false})     = @decompose :dense l
decompose(l :: Conv{false})      = @decompose :conv l
decompose(l :: Deconv{false})    = @decompose :deconv l
decompose(l :: Dropout{false})   = @decompose :dropout l
decompose(l :: Pool{false})      = @decompose :pool l
decompose(l :: Batchnorm{false}) = @decompose :batchnorm l
decompose(l :: Chain{false})     = @decompose :chain l
decompose(l :: Stack{false})     = @decompose :stack l
decompose(l :: Residual{false})  = @decompose :residual l

compose(:: Val{:pointwise}, d) = @compose Pointwise{false} d
compose(:: Val{:dense}, d)     = @compose Dense{false} d
compose(:: Val{:conv}, d)      = @compose Conv{false} d
compose(:: Val{:deconv}, d)    = @compose Deconv{false} d
compose(:: Val{:dropout}, d)   = @compose Dropout{false} d
compose(:: Val{:pool}, d)      = @compose Pool{false} d
compose(:: Val{:batchnorm}, d) = @compose Batchnorm{false} d
compose(:: Val{:chain}, d)     = @compose Chain{false} d
compose(:: Val{:stack}, d)     = @compose Stack{false} d
compose(:: Val{:residual}, d)  = @compose Residual{false} d

# -------- Base Models ------------------------------------------------------- #

decompose(m :: DummyModel) = dict(:dummy)
decompose(m :: RandomModel) = dict(:random)
decompose(m :: RolloutModel) = dict(:rollout)

compose(:: Val{:dummy}, d) = DummyModel()
compose(:: Val{:random}, d) = RandomModel()
compose(:: Val{:rollout}, d) = RolloutModel()

# -------- Neural Models ----------------------------------------------------- #

function decompose(m :: NeuralModel{G, false}) where {G <: Game}
  d = @decompose :model m
  push!(d, :gametype => decompose(G))
  d
end

function compose(:: Val{:model}, d)
  G = compose(d[:gametype])
  @compose NeuralModel{G, false} d
end

# -------- Async Models ------------------------------------------------------ #

function decompose(m :: Async{G}) where {G <: Game}
  @assert !on_gpu(training_model(m)) "Cannot decompose GPU models"
  dict( :async
      , model = decompose(m.model)
      , max_batchsize = decompose(m.max_batchsize)
      , buffersize = decompose(m.buffersize) )
end

function compose(:: Val{:async}, d)
  model = compose(d[:model])
  max_batchsize = compose(d[:max_batchsize])
  buffersize = compose(d[:buffersize])

  Async(model, max_batchsize = max_batchsize, buffersize = buffersize)
end

# -------- Saving Models ----------------------------------------------------- #

"""
    save_model(name, model)

Save `model` under the filename `name` with automatically appended extension
".jtm". Note that the model is first converted to a saveable format, i.e., it is
moved to the CPU and `training_model(model)` is extracted.
"""
function save_model(fname :: String, model :: Model{G}) where {G}

  # Decompose the training_model after bringing it to the cpu
  dict = model |> training_model |> to_cpu |> decompose

  # Save the decomposed model
  BSON.bson(fname * ".jtm", model = dict)

end

"""
    load_model(name)

Load a model from file `name`, where the extension ".jtm" is automatically
appended.
"""
load_model(fname :: String) = BSON.load(fname * ".jtm")[:model] |> compose

# -------- Saving and Loading Datasets --------------------------------------- #

"""
    save_dataset(name, dataset)

Save `dataset` under filename `name` with automatically appended extension
".jtd". Dataset caches are not saved.
"""
function save_dataset(fname :: String, d :: DataSet)

  # Save the file
  BSON.bson(fname * ".jtd", dataset = d)

end

"""
    load_dataset(name)

Load a dataset from file "name", where the extension ".jtd" is automatically
appended.
"""
load_dataset(fname :: String) = BSON.load(fname * ".jtd")[:dataset]

max_stars_repo_name: roSievers/jtac


text:
function Line(arg0::Line)
    return Line((Line,), arg0)
end

function Line(arg0::Vector2D, arg1::Vector2D, arg2::jdouble)
    return Line((Vector2D, Vector2D, jdouble), arg0, arg1, arg2)
end

function Line(arg0::Vector2D, arg1::jdouble, arg2::jdouble)
    return Line((Vector2D, jdouble, jdouble), arg0, arg1, arg2)
end

function contains(obj::Line, arg0::Vector2D)
    return jcall(obj, "contains", jboolean, (Vector2D,), arg0)
end

function copy_self(obj::Line)
    return jcall(obj, "copySelf", Line, ())
end

function distance(obj::Line, arg0::Vector2D)
    return jcall(obj, "distance", jdouble, (Vector2D,), arg0)
end

function empty_hyperplane(obj::Line)
    return jcall(obj, "emptyHyperplane", SubLine, ())
end

function get_angle(obj::Line)
    return jcall(obj, "getAngle", jdouble, ())
end

function get_offset(obj::Line, arg0::JVector)
    return jcall(obj, "getOffset", jdouble, (JVector,), arg0)
end

function get_offset(obj::Line, arg0::Line)
    return jcall(obj, "getOffset", jdouble, (Line,), arg0)
end

function get_offset(obj::Line, arg0::Point)
    return jcall(obj, "getOffset", jdouble, (Point,), arg0)
end

function get_origin_offset(obj::Line)
    return jcall(obj, "getOriginOffset", jdouble, ())
end

function get_point_at(obj::Line, arg0::Vector1D, arg1::jdouble)
    return jcall(obj, "getPointAt", Vector2D, (Vector1D, jdouble), arg0, arg1)
end

function get_reverse(obj::Line)
    return jcall(obj, "getReverse", Line, ())
end

function get_tolerance(obj::Line)
    return jcall(obj, "getTolerance", jdouble, ())
end

function get_transform(::Type{Line}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return jcall(Line, "getTransform", Transform, (jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function intersection(obj::Line, arg0::Line)
    return jcall(obj, "intersection", Vector2D, (Line,), arg0)
end

function is_parallel_to(obj::Line, arg0::Line)
    return jcall(obj, "isParallelTo", jboolean, (Line,), arg0)
end

function project(obj::Line, arg0::Point)
    return jcall(obj, "project", Point, (Point,), arg0)
end

function reset(obj::Line, arg0::Vector2D, arg1::Vector2D)
    return jcall(obj, "reset", void, (Vector2D, Vector2D), arg0, arg1)
end

function reset(obj::Line, arg0::Vector2D, arg1::jdouble)
    return jcall(obj, "reset", void, (Vector2D, jdouble), arg0, arg1)
end

function revert_self(obj::Line)
    return jcall(obj, "revertSelf", void, ())
end

function same_orientation_as(obj::Line, arg0::Hyperplane)
    return jcall(obj, "sameOrientationAs", jboolean, (Hyperplane,), arg0)
end

function set_angle(obj::Line, arg0::jdouble)
    return jcall(obj, "setAngle", void, (jdouble,), arg0)
end

function set_origin_offset(obj::Line, arg0::jdouble)
    return jcall(obj, "setOriginOffset", void, (jdouble,), arg0)
end

function to_space(obj::Line, arg0::JVector)
    return jcall(obj, "toSpace", Vector2D, (JVector,), arg0)
end

function to_space(obj::Line, arg0::Point)
    return jcall(obj, "toSpace", Vector2D, (Point,), arg0)
end

function to_sub_space(obj::Line, arg0::JVector)
    return jcall(obj, "toSubSpace", Vector1D, (JVector,), arg0)
end

function to_sub_space(obj::Line, arg0::Point)
    return jcall(obj, "toSubSpace", Vector1D, (Point,), arg0)
end

function translate_to_point(obj::Line, arg0::Vector2D)
    return jcall(obj, "translateToPoint", void, (Vector2D,), arg0)
end

function whole_hyperplane(obj::Line)
    return jcall(obj, "wholeHyperplane", SubLine, ())
end

function whole_space(obj::Line)
    return jcall(obj, "wholeSpace", PolygonsSet, ())
end


text:
abstract type AbstractProblemMinimallyAugmented end
abstract type AbstractCodim2EigenSolver <: AbstractEigenSolver end

@inline getLens(pb::AbstractProblemMinimallyAugmented) = pb.lens
getsolver(eig::AbstractCodim2EigenSolver) = eig.eigsolver

for op in (:FoldProblemMinimallyAugmented, :HopfProblemMinimallyAugmented)
	@eval begin
		"""
		$(TYPEDEF)

		Structure to encode Fold / Hopf functional based on a Minimally Augmented formulation.

		# Fields

		$(FIELDS)
		"""
		struct $op{TF, TJ, TJa, Td2f, Tl <: Lens, vectype, S <: AbstractLinearSolver, Sa <: AbstractLinearSolver, Sbd <: AbstractBorderedLinearSolver, Sbda <: AbstractBorderedLinearSolver} <: AbstractProblemMinimallyAugmented
			"Function F(x, p) = 0"
			F::TF
			"Jacobian of F w.r.t. x"
			J::TJ
			"Adjoint of the Jacobian of F"
			Jᵗ::TJa
			"Hessian of F"
			d2F::Td2f
			"parameter axis for the codim 2 point"
			lens::Tl
			"close to null vector of Jᵗ"
			a::vectype
			"close to null vector of J"
			b::vectype
			"vector zero, to avoid allocating it many times"
			zero::vectype
			"linear solver. Used to invert the jacobian of MA functional"
			linsolver::S
			"linear solver for the jacobian adjoint"
			linsolverAdjoint::Sa
			"bordered linear solver"
			linbdsolver::Sbd
			"linear bordered solver for the jacobian adjoint"
			linbdsolverAdjoint::Sbda
			"whether the Jacobian is Symmetric, avoid computing Jᵗ"
			issymmetric::Bool
		end

		@inline hasHessian(pb::$op{TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd, Sbda}) where {TF, TJ, TJa, Td2f, Tp, Tl, vectype, S, Sa, Sbd, Sbda} = Td2f != Nothing

		@inline hasAdjoint(pb::$op{TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd, Sbda}) where {TF, TJ, TJa, Td2f, Tl, vectype, S, Sa, Sbd, Sbda} = TJa != Nothing
	end
end

@inline issymmetric(pb::AbstractProblemMinimallyAugmented) = pb.issymmetric

function applyJacobian(pb::AbstractProblemMinimallyAugmented, x, par, dx, transposeJac = false)
	if issymmetric(pb)
		return apply(pb.J(x, par), dx)
	else
		if transposeJac == false
			return apply(pb.J(x, par), dx)
		else
			if hasAdjoint(pb)
				return apply(pb.Jᵗ(x, par), dx)
			else
				return apply(transpose(pb.J(x, par)), dx)
			end
		end
	end
end

function detectCodim2Parameters(detectCodim2Bifurcation, options_cont; kwargs...)
	if detectCodim2Bifurcation > 0
		if get(kwargs, :updateMinAugEveryStep, 0) == 0
			@error "You ask for detection of codim 2 bifurcations but passed the option `updateMinAugEveryStep = 0`. The bifurcation detection algorithm may not work faithfully. Please use `updateMinAugEveryStep > 0`."
		end
		return setproperties(options_cont; detectBifurcation = 0, detectEvent = detectCodim2Bifurcation, detectFold = false)
	else
		return options_cont
	end
end

"""
$(SIGNATURES)

This function turns an initial guess for a Fold/Hopf point into a solution to the Fold/Hopf problem based on a Minimally Augmented formulation. The arguments are as follows
- `F  = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `J  = (x, p) -> d_xF(x, p)` associated jacobian
- `br` results returned after a call to [continuation](@ref Library-Continuation)
- `ind_bif` bifurcation index in `br`
- `lens` parameter axis used to locate the Fold/Hopf point.
- `options::NewtonPar`

# Optional arguments:
- `issymmetric` whether the Jacobian is Symmetric (for Fold)
- `Jᵗ = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of Matrix / Sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoids recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1,v2]`.
- `normN = norm`
- `options` You can pass newton parameters different from the ones stored in `br` by using this argument `options`.
- `bdlinsolver` bordered linear solver for the constraint equation
- `startWithEigen = false` whether to start the Minimally Augmented problem with information from eigen elements
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver

!!! tip "ODE problems"
    For ODE problems, it is more efficient to use the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`
"""
function newton(F, J, br::AbstractBranchResult, ind_bif::Int64; Jᵗ = nothing, d2F = nothing, normN = norm, options = br.contparams.newtonOptions, startWithEigen = false, issymmetric::Bool = false, kwargs...)
	@assert length(br.specialpoint) > 0 "The branch does not contain bifurcation points"
	if br.specialpoint[ind_bif].type == :hopf
		d2Fc = isnothing(d2F) ? nothing : (x,p,dx1,dx2) -> BilinearMap((_dx1, _dx2) -> d2F(x,p,_dx1,_dx2))(dx1,dx2)
		return newtonHopf(F, J, br, ind_bif; Jᵗ = Jᵗ, d2F = d2Fc, normN = normN, options = options, startWithEigen = startWithEigen, kwargs...)
	else
		return newtonFold(F, J, br, ind_bif; issymmetric = issymmetric, Jᵗ = Jᵗ, d2F = d2F, normN = normN, options = options, startWithEigen = startWithEigen, kwargs...)
	end
end

"""
$(SIGNATURES)

Codimension 2 continuation of Fold / Hopf points. This function turns an initial guess for a Fold/Hopf point into a curve of Fold/Hopf points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `br` results returned after a call to [continuation](@ref Library-Continuation)
- `ind_bif` bifurcation index in `br`
- `lens2` parameters used for the vector field
- `options_cont = br.contparams` arguments to be passed to the regular [continuation](@ref Library-Continuation)

# Optional arguments:
- `issymmetric` whether the Jacobian is Symmetric (for Fold)
- `Jᵗ = (x, p) -> transpose(d_xF(x, p))` associated jacobian transpose
- `d2F = (x, p, v1, v2) -> d2F(x, p, v1, v2)` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`. This helps solving the linear problem associated to the minimally augmented formulation.
- `d3F = (x, p, v1, v2, v3) -> d3F(x, p, v1, v2, v3)` this is the third derivative of `F` computed at `(x, p)` and evaluated at `(v1, v2, v3)`. This is used to detect **Bautin** bifurcation.
- `bdlinsolver` bordered linear solver for the constraint equation
- `updateMinAugEveryStep` update vectors `a,b` in Minimally Formulation every `updateMinAugEveryStep` steps
- `startWithEigen = false` whether to start the Minimally Augmented problem with information from eigen elements
- `detectCodim2Bifurcation ∈ {0,1,2}` whether to detect Bogdanov-Takens, Bautin and Cusp. If equals `1` non precise detection is used. If equals `2`, a bisection method is used to locate the bifurcations.
- `kwargs` keywords arguments to be passed to the regular [continuation](@ref Library-Continuation)

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Hopf point in `br` you want to refine.

!!! tip "Jacobian transpose"
    The adjoint of the jacobian `J` is computed internally when `Jᵗ = nothing` by using `transpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jᵗ = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! tip "ODE problems"
    For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option `bdlinsolver = MatrixBLS()`
"""
function continuation(F, J,
				br::AbstractBranchResult, ind_bif::Int64,
				lens2::Lens, options_cont::ContinuationPar = br.contparams ;
				startWithEigen = false,
				issymmetric::Bool = false,
				Jᵗ = nothing,
				d2F = nothing,
				d3F = nothing,
				detectCodim2Bifurcation::Int = 0,
				kwargs...)
	@assert length(br.specialpoint) > 0 "The branch does not contain bifurcation points"
	# options to detect codim2 bifurcations
	computeEigenElements = options_cont.detectBifurcation > 0
	_options_cont = detectCodim2Parameters(detectCodim2Bifurcation, options_cont; kwargs...)

	if br.specialpoint[ind_bif].type == :hopf
		# redefine the multilinear form to accept complex arguments
		d2Fc = isnothing(d2F) ? nothing : (x,p,dx1,dx2) -> BilinearMap((_dx1, _dx2) -> d2F(x,p,_dx1,_dx2))(dx1,dx2)
		d3Fc = isnothing(d3F) ? nothing : (x,p,dx1,dx2,dx3) -> TrilinearMap((_dx1, _dx2, _dx3) -> d3F(x,p,_dx1,_dx2,_dx3))(dx1,dx2,dx3)
		return continuationHopf(F, J, br, ind_bif, lens2, _options_cont; Jᵗ = Jᵗ, d2F = d2Fc, d3F = d3Fc, startWithEigen = startWithEigen, computeEigenElements = computeEigenElements, kwargs...)
	else
		return continuationFold(F, J, br, ind_bif, lens2, _options_cont;
			issymmetric = issymmetric,
			Jᵗ = Jᵗ,
			d2F = d2F,
			startWithEigen = startWithEigen,
			computeEigenElements = computeEigenElements,
			kwargs...)
	end
end
####################################################################################################
# branch switching at bt
function continuation(F, dF, d2F, d3F,
			br::ContResult{Ta, Teigvals, Teigvecbr, Biftype, Ts, Tparc, Tfunc, Tpar, Tl},
			ind_bif::Int, options_cont::ContinuationPar = br.contparams;
			Jᵗ = nothing,
			δ::Real = 1e-8, δp = nothing, ampfactor::Real = 1,
			nev = options_cont.nev,
			issymmetric = false,
			detectCodim2Bifurcation::Int = 0,
			Teigvec = getvectortype(br),
			scaleζ = norm,
			startWithEigen = false,
			autodiff = false,
			kwargs...) where {Ta, Teigvals, Teigvecbr, Biftype, Ts, Tparc, Tfunc <: AbstractProblemMinimallyAugmented, Tpar, Tl <: Lens}

		verbose = get(kwargs, :verbosity, 0) > 0 ? true : false
		verbose && println("--> Considering bifurcation point:"); _show(stdout, br.specialpoint[ind_bif], ind_bif)

		@assert br.specialpoint[ind_bif].type in (:bt,:zh) "Only branching from Bogdanov-Takens and Zero-Hopf (for now)"

		# functional
		prob = br.functional

		# continuation parameters
		computeEigenElements = options_cont.detectBifurcation > 0
		optionsCont = detectCodim2Parameters(detectCodim2Bifurcation, options_cont; kwargs...)

		# higher order differentials
		d2Fc = isnothing(d2F) ? nothing : (x,p,dx1,dx2) -> BilinearMap((_dx1, _dx2) -> d2F(x,p,_dx1,_dx2))(dx1,dx2)
		d3Fc = isnothing(d3F) ? nothing : (x,p,dx1,dx2,dx3) -> TrilinearMap((_dx1, _dx2, _dx3) -> d3F(x,p,_dx1,_dx2,_dx3))(dx1,dx2,dx3)

		# scalar type
		Ty = eltype(Teigvec)

		# compute the normal form of the bifurcation point
		nf = computeNormalForm(F, dF, d2F, d3F, br, ind_bif; Jᵗ = Jᵗ, δ = δ, nev = nev, verbose = verbose, issymmetric = issymmetric, Teigvec = Teigvec, scaleζ = scaleζ, autodiff = autodiff)

		# compute predictor for point on new branch
		ds = isnothing(δp) ? optionsCont.ds : δp

		if prob isa FoldProblemMinimallyAugmented
			# define guess for the first Hopf point on the branch
			pred = predictor(nf, Val(:HopfCurve), ds)

			# new continuation parameters
			parcont = pred.hopf(ds)

			# new full parameters
			params = set(set(nf.params, nf.lens[2], parcont[2]), nf.lens[1], parcont[1])

			# guess for the Hopf point
			hopfpt = BorderedArray(nf.x0 .+ pred.x0(ds), [parcont[1], pred.ω(ds)])

			# estimates for eigenvectors for ±iω
			ζ = pred.EigenVec(ds)
			ζstar = pred.EigenVecAd(ds)

			# put back original options
			@set! optionsCont.newtonOptions.eigsolver =
								getsolver(optionsCont.newtonOptions.eigsolver)
			@set! optionsCont.newtonOptions.linsolver = prob.linsolver

			branch, u, τ = continuationHopf(F, dF, hopfpt, params,
					nf.lens...,
					ζ, ζstar,
					optionsCont;
					d2F = d2Fc, d3F = d3Fc,
					Jᵗ = Jᵗ,
					bdlinsolver = prob.linbdsolver,
					startWithEigen = startWithEigen,
					computeEigenElements = computeEigenElements,
					kwargs...
					)
			return Branch(branch, nf), u, τ

		else
			@assert prob isa HopfProblemMinimallyAugmented
			pred = predictor(nf, Val(:FoldCurve), 0.)
			# new continuation parameters
			parcont = pred.fold(ds)

			# new full parameters
			params = set(set(nf.params, nf.lens[2], parcont[2]), nf.lens[1], parcont[1])

			# guess for the fold point
			foldpt = BorderedArray(nf.x0 .+ 0 .* pred.x0(ds), parcont[1])

			# estimates for null eigenvectors
			ζ = pred.EigenVec(ds)
			ζstar = pred.EigenVecAd(ds)

			# put back original options
			@set! optionsCont.newtonOptions.eigsolver =
								getsolver(optionsCont.newtonOptions.eigsolver)
			@set! optionsCont.newtonOptions.linsolver = prob.linsolver
			# @set! optionsCont.detectBifurcation = 0
			# @set! optionsCont.detectEvent = 0

			branch, u, τ = continuationFold(F, dF, foldpt, params,
					nf.lens...,
					ζstar, ζ,
					optionsCont;
					d2F = d2Fc, #d3F = d3Fc,
					Jᵗ = Jᵗ,
					bdlinsolver = prob.linbdsolver,
					startWithEigen = startWithEigen,
					computeEigenElements = computeEigenElements,
					kwargs...
					)
			return Branch(branch, nf), u, τ
		end
end

"""
$(SIGNATURES)

This function uses information in the branch to detect codim 2 bifurcations like BT, ZH and Cusp.
"""
function correctBifurcation(contres::ContResult)
	if contres.functional isa AbstractProblemMinimallyAugmented == false
		return contres
	end
	if contres.functional isa FoldProblemMinimallyAugmented
		conversion = Dict(:bp => :bt, :hopf => :zh, :fold => :cusp, :nd => :nd)
	elseif contres.functional isa HopfProblemMinimallyAugmented
		conversion = Dict(:bp => :zh, :hopf => :hh, :fold => :nd, :nd => :nd, :ghbt => :bt)
	else
		throw("Error! this should not occur. Please open an issue on the website of BifurcationKit.jl")
	end
	for (ind, bp) in pairs(contres.specialpoint)
		if bp.type in keys(conversion)
			@set! contres.specialpoint[ind].type = conversion[bp.type]
		end
	end
	return contres
end


text:
const WSInRoom = Set{WebSocket}()
const RoomPass = Dict{String, String}()
const RoomDict = Dict{String, Set{WebSocket}}()
const Usernames = Dict{String, Set{String}}()
const getRoomfromWS = Dict{WebSocket, String}()
const getHandlefromWS = Dict{WebSocket, String}()

function makeroom(currws, logdetails)
    !validJSON(logdetails, ["handle", "roomname", "roompass"]) && return
    handle = logdetails["handle"]
    roomname = logdetails["roomname"]
    roompass = logdetails["roompass"]

    if currws in WSInRoom
        # WebSocket already in Room
        return
    elseif !validName(handle)
        # Whitespace only Handle
        return errormsg(currws, "login-invalid-handle")
    elseif !validName(roomname)
        # Whitespace only Name
        return errormsg(currws, "login-invalid-roomname")
    elseif haskey(RoomPass, roomname)
        # Room Already Exists
        errormsg(currws, "login-room-exists")
        return
    end

    push!(WSInRoom, currws)
    getRoomfromWS[currws] = roomname
    getHandlefromWS[currws] = handle

    RoomPass[roomname] = roompass
    Usernames[roomname] = Set{String}()
    push!(Usernames[roomname], handle)
    RoomDict[roomname] = Set{WebSocket}()
    push!(RoomDict[roomname], currws)

    roomsuccess(currws)
end

function joinroom(currws, logdetails)
    !validJSON(logdetails, ["handle", "roomname", "roompass"]) && return
        handle = logdetails["handle"]
        roomname = logdetails["roomname"]
        roompass = logdetails["roompass"]

    if currws in WSInRoom
        # WebSocket already in Room
        return
    elseif !validName(handle)
        # Whitespace only Handle
        return errormsg(currws, "login-invalid-handle")
    elseif !validName(roomname)
        # Whitespace only Name
        return errormsg(currws, "login-invalid-roomname")
    elseif !haskey(RoomPass, roomname)
        # Room Doesn't Exist
        return errormsg(currws, "login-room-missing")
    elseif RoomPass[roomname] != roompass
        # Wrong Password
        return errormsg(currws, "login-wrong-password")
    elseif handle in Usernames[roomname]
        # Duplicated Username
        return errormsg(currws, "login-duplicate-username")
    end

    push!(WSInRoom, currws)
    getRoomfromWS[currws] = roomname
    getHandlefromWS[currws] = handle

    push!(Usernames[roomname], handle)
    push!(RoomDict[roomname], currws)
    roomsuccess(currws)
end

roomsize(roomname) = length(Usernames[roomname])

# WebSocket disconnected.
function removereferences(ws)
    if ws in WSInRoom
        roomname = getRoomfromWS[ws]
        handle = getHandlefromWS[ws]

        pop!(WSInRoom, ws)
        pop!(getRoomfromWS, ws)
        pop!(getHandlefromWS, ws)

        if handle in Usernames[roomname]
            pop!(Usernames[roomname], handle)
        end
        if ws in RoomDict[roomname]
            # Player Disconnected from Room
            pop!(RoomDict[roomname], ws)
            playerdisconnected(RoomDict[roomname], handle)
        end

        if roomsize(roomname) == 0
            # Clean-up Empty Room (Free-up Password)
            pop!(RoomPass, roomname)
            delete!(RoomGames, roomname)
            delete!(SpecDict, roomname)
        end

        dcplayerleft(handle, roomname)
    end

    nothing
end

function leaveroom(currws)
    removereferences(currws)
    msg = Dict{String, Any}()
    msg["responsetype"] = "room-left-room"
    msg["server-time"] = Dates.now()

    writeguarded(currws, JSON.json(msg))

    nothing
end

function roomsuccess(currws)
    # Send join-room
    msg = Dict{String, Any}()
    msg["responsetype"] = "room-joined-room"

    roomname = getRoomfromWS[currws]
    handle = getHandlefromWS[currws]

    # Includes yourself in player list
    msg["users"] = getUsers(roomname)
    msg["game-ongoing"] = haskey(RoomGames, roomname)
    msg["room-name"] = roomname
    msg["server-time"] = Dates.now()
    # Only include online players in spectation-list
    msg["spectation-list"] = filter(x -> x in getUsers(roomname), getSpectators(roomname))
    if haskey(RoomGames, roomname)
        # Send Game State to Player
        msg["game-status"] = gameState(RoomGames[roomname])
    end

    dcplayerjoined(handle, roomname)
    writeguarded(currws, JSON.json(msg))

    # Send the rest of players a player-join message
    playerjoined(RoomDict[roomname], handle, roomname)
    nothing
end

function playerjoined(listofws, handle, roomname)
    msg = Dict{String, Any}()
    msg["responsetype"] = "room-player-joined"
    msg["user"] = handle
    msg["is-spectator"] = handle in getSpectators(roomname)
    msg["server-time"] = Dates.now()
    return broadcastmsg(listofws, JSON.json(msg),
                        cond = x -> (getHandlefromWS[x] != handle))
end

function playerdisconnected(listofws, handle)
    msg = Dict{String, Any}()
    msg["responsetype"] = "room-player-disconnected"
    msg["user"] = handle
    msg["server-time"] = Dates.now()

    return broadcastmsg(listofws, JSON.json(msg),
                        cond = x -> (getHandlefromWS[x] != handle))
end


text:
using Katz
using Base.Test
using LightGraphs
using GraphMatrices
import GraphMatrices: KatzAdjacency
a = 0.15
g = LightGraphs.Datasets.BullGraph()
function test1(g, a::Number)
    n = nv(g)
    A = sparse(CombinatorialAdjacency(g))
    v = stdbasis(n, 1)
    k = katz(A, a, v)
    k2 = katz(A, a, ones(n))
    @test_approx_eq_eps k2 k 1e-5
end

function test_iterative(g, a)
    K = KatzAdjacency(CombinatorialAdjacency(g), a)
    n = nv(g)
    y = K*ones(n)
    z, convhist = ikatz(K,ones(n))
    x = K.A*z
    @show convhist.isconverged
    x_lu = katz(sparse(K.A), a, ones(n))
    @test norm(x_lu - x) < 1e-6
    @show ρ(convhist)
end
datasets = [:BullGraph,
            :TutteGraph,
            :DesarguesGraph]
for a in [0.1,0.01,0.001]
for s in datasets
    g = eval(:(LightGraphs.Datasets.$s()))
    println("test1 on $s, $a")
    test1(g, a)
    test_iterative(g,a)
end
end

for s in ["Newman/football"]
    g = MDGraph(s,:read)
    test1(g, α)
    test_iterative(g,α)
end


text:
function AngularCoordinates()
    return AngularCoordinates(())
end

function AngularCoordinates(arg0::FieldRotation)
    return AngularCoordinates((FieldRotation,), arg0)
end

function AngularCoordinates(arg0::PVCoordinates, arg1::PVCoordinates)
    return AngularCoordinates((PVCoordinates, PVCoordinates), arg0, arg1)
end

function AngularCoordinates(arg0::PVCoordinates, arg1::PVCoordinates, arg2::PVCoordinates, arg3::PVCoordinates, arg4::jdouble)
    return AngularCoordinates((PVCoordinates, PVCoordinates, PVCoordinates, PVCoordinates, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function AngularCoordinates(arg0::Rotation, arg1::Vector3D)
    return AngularCoordinates((Rotation, Vector3D), arg0, arg1)
end

function AngularCoordinates(arg0::Rotation, arg1::Vector3D, arg2::Vector3D)
    return AngularCoordinates((Rotation, Vector3D, Vector3D), arg0, arg1, arg2)
end

function add_offset(obj::AngularCoordinates, arg0::AngularCoordinates)
    return jcall(obj, "addOffset", AngularCoordinates, (AngularCoordinates,), arg0)
end

function apply_to(obj::AngularCoordinates, arg0::FieldPVCoordinates)
    return jcall(obj, "applyTo", FieldPVCoordinates, (FieldPVCoordinates,), arg0)
end

function apply_to(obj::AngularCoordinates, arg0::PVCoordinates)
    return jcall(obj, "applyTo", PVCoordinates, (PVCoordinates,), arg0)
end

function apply_to(obj::AngularCoordinates, arg0::TimeStampedFieldPVCoordinates)
    return jcall(obj, "applyTo", TimeStampedFieldPVCoordinates, (TimeStampedFieldPVCoordinates,), arg0)
end

function apply_to(obj::AngularCoordinates, arg0::TimeStampedPVCoordinates)
    return jcall(obj, "applyTo", TimeStampedPVCoordinates, (TimeStampedPVCoordinates,), arg0)
end

function create_from_modified_rodrigues(::Type{AngularCoordinates}, arg0::Vector{Vector{jdouble}})
    return jcall(AngularCoordinates, "createFromModifiedRodrigues", AngularCoordinates, (Vector{Vector{jdouble}},), arg0)
end

function estimate_rate(::Type{AngularCoordinates}, arg0::Rotation, arg1::Rotation, arg2::jdouble)
    return jcall(AngularCoordinates, "estimateRate", Vector3D, (Rotation, Rotation, jdouble), arg0, arg1, arg2)
end

function get_modified_rodrigues(obj::AngularCoordinates, arg0::jdouble)
    return jcall(obj, "getModifiedRodrigues", Vector{Vector{jdouble}}, (jdouble,), arg0)
end

function get_rotation(obj::AngularCoordinates)
    return jcall(obj, "getRotation", Rotation, ())
end

function get_rotation_acceleration(obj::AngularCoordinates)
    return jcall(obj, "getRotationAcceleration", Vector3D, ())
end

function get_rotation_rate(obj::AngularCoordinates)
    return jcall(obj, "getRotationRate", Vector3D, ())
end

function revert(obj::AngularCoordinates)
    return jcall(obj, "revert", AngularCoordinates, ())
end

function shifted_by(obj::AngularCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", AngularCoordinates, (jdouble,), arg0)
end

function subtract_offset(obj::AngularCoordinates, arg0::AngularCoordinates)
    return jcall(obj, "subtractOffset", AngularCoordinates, (AngularCoordinates,), arg0)
end

function to_derivative_structure_rotation(obj::AngularCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructureRotation", FieldRotation, (jint,), arg0)
end

function to_univariate_derivative1_rotation(obj::AngularCoordinates)
    return jcall(obj, "toUnivariateDerivative1Rotation", FieldRotation, ())
end


text:
const EnvPair = Pair{String,Union{String,Nothing}}

struct CLISpec
    options::Cmd
    env::Vector{EnvPair}
end

function Base.show(io::IO, ::MIME"text/plain", spec::CLISpec)
    get(io, :typeinfo, Any) == CLISpec || print(io, "CLISpec: ")
    for (k, v) in spec.env
        print(io, k, '=')
        if v === nothing
            printstyled(io, '❌'; color = :red)
        else
            show(io, v)
        end
        print(io, ' ')
    end
    printstyled(io, raw"$julia"; color = :blue)
    isempty(spec.options) || print(io, ' ')
    join(io, spec.options, ' ')
end

CLISpec(spec::AbstractDict) = CLISpec(``, collect(EnvPair, spec))

function CLISpec(spec::AbstractVector)
    env = EnvPair[]
    options = ``
    for s in spec
        if s isa Cmd
            options = `$options $s`
        else
            push!(env, s)
        end
    end
    return CLISpec(options, env)
end

function env_from_spec(spec::CLISpec)
    env = copy(ENV)
    for (k, v) in spec.env
        if v === nothing
            pop!(env, k)
        else
            env[k] = v
        end
    end
    return env
end

if VERSION < v"1.1"
    function _readboth(cmd)
        pipe = Pipe()
        try
            reader = @async read(pipe)
            proc = run(pipeline(cmd, stdout = pipe, stderr = pipe, stdin = devnull))
            close(pipe)
            return proc, fetch(reader)
        finally
            close(pipe)
        end
    end
else
    function _readboth(cmd)
        io = IOBuffer()
        proc = run(pipeline(cmd, stdout = io, stderr = io, stdin = devnull))
        return proc, take!(io)
    end
end

struct IncludeResult
    script::String
    spec::CLISpec
    output::String
    success::Bool
    proc::Base.Process
end

function Base.show(io::IO, ::MIME"text/plain", result::IncludeResult)
    color = result.success ? :green : :red

    printstyled(io, "┌"; color = color)
    printstyled(io, " Test result: ", bold = true)
    if result.success
        printstyled(io, "All success"; color = color, bold = true)
    else
        printstyled(io, "At least one failure"; color = color, bold = true)
    end
    println(io)

    printstyled(io, '│'; color = color)
    printstyled(io, " Script : "; color = :blue)
    printstyled(io, result.script; bold = true)
    println(io)

    printstyled(io, '│'; color = color)
    printstyled(io, " Command: "; color = :blue)
    show(IOContext(io, :typeinfo => typeof(result.spec)), MIME"text/plain"(), result.spec)
    println(io)

    for line in eachline(IOBuffer(result.output))
        printstyled(io, '│'; color = color)
        println(io, ' ', line)
    end
    printstyled(io, "└"; color = color)
    println(io)
end

"""
    PerformanceTestTools.include_foreach(script, speclist)

Like [`@include_foreach`](@ref) but relative path `script` is resolved
with respect to the current working directory, instead of the file in
which this function is called.
"""
function include_foreach(script, speclist0; parallel::Bool = false, __include = _include)
    speclist = CLISpec[]
    hasnothing = false
    for x in speclist0
        if x === nothing
            hasnothing = true
        else
            push!(speclist, CLISpec(x))
        end
    end

    script = abspath(script)
    if parallel
        _map = (f, xs) -> @async asyncmap(f, xs; ntasks = Sys.CPU_THREADS)
    else
        _map = map
    end
    results = _map(speclist) do spec
        env = env_from_spec(spec)
        code = """
        $(Base.load_path_setup_code())
        include($(repr(script)))
        """
        cmd = setenv(`$(_julia_cmd()) -e $code $(spec.options)`, env)
        @info "Running `$script` in a subprocess..." spec
        if parallel
            proc, output = _readboth(ignorestatus(cmd))
            output = String(output)
        else
            proc = run(pipeline(
                ignorestatus(cmd);
                stdout = stdout,
                stderr = stderr,
                stdin = devnull,
            ))
            output = ""
        end
        @info "Running `$script` in a subprocess...DONE"
        return IncludeResult(script, spec, output, success(proc), proc)
    end
    try
        if hasnothing
            @info "Including `$script` in this process..."
            __include(script)
            @info "Including `$script` in this process...DONE"
        end
    finally
        if parallel
            istaskdone(results) || @info "Waiting for parallel tasks..."
            results = fetch(results)
            for r in results
                show(stdout, "text/plain", r)
            end
        end
        @test all(x -> x.success, results)
    end
end

@doc read(joinpath(@__DIR__, "foreach.md"), String) :(@include_foreach)

macro include_foreach(script, args...)
    dir = dirname(string(__source__.file))
    args = map(args) do ex
        isexpr(ex, :(=)) ? Expr(:kw, ex.args...) : ex
    end
    esc(:($include_foreach($joinpath($dir, $script), __include = include, $(args...))))
end


text:
function get_estimated_measurement(obj::EstimationsProvider, arg0::jint)
    return jcall(obj, "getEstimatedMeasurement", EstimatedMeasurement, (jint,), arg0)
end

function get_number(obj::EstimationsProvider)
    return jcall(obj, "getNumber", jint, ())
end


text:
"""

```
mddf_linkedcells(trajectory::Trajectory, options::Options,samples::Samples, RNG, mddf_compute!)
```

Computes the MDDF using linked cells, serial version.

"""
function mddf_linkedcells(
    trajectory::Trajectory,
    options::Options,
    samples::Samples,
    RNG,
    mddf_compute!,
)

    # Initializing the structure that carries all results
    R = Result(trajectory, options)

    # Data structure to be passed to mddf_frame
    framedata = FrameData(trajectory, R)

    # Print some information about this run
    options.silent || title(R, trajectory.solute, trajectory.solvent)

    # Computing all minimum-distances
    if !options.silent
        progress = Progress(R.nframes_read, 1)
    end
    for iframe = 1:R.lastframe_read

        # reading coordinates of next frame
        nextframe!(trajectory)
        if iframe < options.firstframe
            continue
        end
        if iframe % options.stride != 0
            continue
        end
        mddf_compute!(iframe, framedata, options, RNG, R)

        options.silent || next!(progress)
    end # frames
    closetraj(trajectory)

    # Setup the final data structure with final values averaged over the number of frames,
    # sampling, etc, and computes final distributions and integrals
    finalresults!(R, options, trajectory, samples)
    options.silent || println(bars)

    return R

end


text:
function zero!(arr::A; kwargs...) where A <: AbstractArray{T, N} where T where N
    arr[:] .= zero(T)
    return arr
end

Base.zero(::Type{Any}) = 0

Base.zero(::Type{T}) where T = 0

Base.zero(::Type{T}) where T <: AbstractString = ""

Base.zero(::Type{T}) where T <: Nothing = nothing

Base.zero(::Type{T}) where T <: Missing = missing


text:
"""
    HMMApproximator

Sampleable approximator object.
"""
mutable struct HMMApproximator{T<:HMM}
  hmm::T
  obs::Vector
  zi::Vector
  llfw::Vector
  z::Vector
  n::Int
  nk::Int
  k::Int
end

"""
    sample(approx::HMMApproximator,n=1)

Sample `n` sequences from `approx` using generalized Backward procedure.
"""
function sample(approx::HMMApproximator)
  @unpack hmm, zi, llfw, k, n, z, nk = approx
  @unpack logP = hmm

  p = Vector{Float64}(undef, nk)
  s = Vector{Int}(undef, n)
  r = k - 1

  ll_n = 0

  # last (k-1) elements
  idx = rand(Categorical(exp.(z)))
  s[n-r+1:n] = decode(idx, nk, r)
  ll_n += zi[end][idx]

  # intermediate elements
  for i = n-r:-1:2
    for s1 = 1:nk
      s[i] = s1
      p[s1] = logP[s[i+r-1], s[i+r]] + llfw[i][encode(s[i:i+r], nk)] .- zi[i][encode(s[i+1:i+r], nk)]
    end
    idx = rand(Categorical(exp.(p)))
    s[i] = idx
    ll_n += p[idx]
  end

  # first element
  for s1 = 1:nk
    s[1] = s1
    p[s1] = llfw[1][encode(s[1:k], nk)] - zi[1][encode(s[2:k], nk)]
  end
  idx = rand(Categorical(exp.(p)))
  s[1] = idx
  ll_n += p[idx]

  s, ll_n
end

function sample(approx::HMMApproximator, nsamples::Int)
  samples = Vector{Vector{Int}}()
  lls     = Vector{Float64}()

  for i in 1:nsamples
    s, ll = sample(approx)
    push!(samples, s)
    push!(lls, ll)
  end

  samples, lls
end

function viterbi(approx::HMMApproximator)
  @unpack hmm, zi, llfw, k = approx
  @unpack logP = hmm

  nk = size(logP, 1)
  r = k - 1
  n = length(zi) + r
  z = zi[end] .- logsumexp(zi[end])

  p = Vector{Float64}(undef, nk)
  s = Vector{Int}(undef, n)
  ll = 0

  # last (k-1) elements
  idx = argmax(z)
  s[n-r+1:n] = decode(idx, nk, r)
  ll += zi[end][idx]

  @inbounds for i = n-r:-1:2
    for s1 = 1:nk
      s[i] = s1
      p[s1] = logP[s[i+r-1], s[i+r]] + llfw[i][encode(s[i:i+r], nk)] - zi[i][encode(s[i+1:i+r], nk)]
    end
    idx = argmax(p)
    s[i] = idx
    ll += p[idx]
  end

  # first element
  @inbounds for s1 = 1:nk
    s[1] = s1
    p[s1] = llfw[1][encode(s[1:k], nk)] - zi[1][encode(s[2:k], nk)]
  end
  idx = argmax(p)
  s[1] = idx
  ll += p[idx]

  s, ll
end

function mcmc!(samples::Vector, approx::HMMApproximator, nIter=1000)
  @unpack hmm, obs = approx
  n = length(obs)

  @assert length(samples) == nIter "dimension mismatch"
  @assert length(first(samples)) == n "dimension mismatch"

  s_prev, p_a_prev = sample(approx)
  p_prev = logpdf(hmm, s_prev, obs)
  α = 0

  for j in 1:nIter
    s, p_a = sample(approx)
    p = logpdf(hmm, s, obs)
    αⱼ = min(1, exp(p - p_a + p_a_prev - p_prev))
    s_prev = rand() < αⱼ ? s : s_prev
    samples[j] = s_prev
    α = ((j-1)*α + αⱼ)/j
  end

  samples, α
end

function mcmc(approx::HMMApproximator, nIter=1000)
  n = length(approx.obs)
  samples = [Vector{Int}(undef, n) for _ in 1:nIter]
  mcmc!(samples, approx, nIter)
end

function Base.show(io::IO, h::HMMApproximator)
    print(io, "HMMApproximator(k=$(h.k))")
end


text:
module TestLocks

using Base.Experimental: @sync
using Reagents
using Test

include("../../../examples/locks.jl")

end  # module


text:
# This file is part of the Taylor1Series.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

"""
    AbstractSeries{T<:Number} <: Number

Parameterized abstract type for [`Taylor1`](@ref),
[`HomogeneousPolynomial`](@ref) and [`TaylorN`](@ref).
"""
abstract type AbstractSeries{T<:Number} <: Number end

## Constructors ##

######################### Taylor1
"""
    Taylor1{T<:Number} <: AbstractSeries{T}

DataType for polynomial expansions in one independent variable.

**Fields:**

- `coeffs :: Array{T,1}` Expansion coefficients; the ``i``-th
    component is the coefficient of degree ``i-1`` of the expansion.
- `order  :: Int` Maximum order (degree) of the polynomial.

Note that `Taylor1` variables are callable. For more information, see
[`evaluate`](@ref).
"""
struct Taylor1{T<:Number} <: AbstractSeries{T}
    coeffs :: Array{T,1}
    order :: Int

    ## Inner constructor ##
    function Taylor1{T}(coeffs::Array{T,1}, order::Int) where T<:Number
        resize_coeffs1!(coeffs, order)
        return new{T}(coeffs, order)
    end
end

## Outer constructors ##
Taylor1(x::Taylor1{T}) where {T<:Number} = x
Taylor1(coeffs::Array{T,1}, order::Int) where {T<:Number} = Taylor1{T}(coeffs, order)
Taylor1(coeffs::Array{T,1}) where {T<:Number} = Taylor1(coeffs, length(coeffs)-1)
function Taylor1(x::T, order::Int) where {T<:Number}
    v = fill(zero(x), order+1)
    v[1] = x
    return Taylor1(v, order)
end

# Methods using 1-d views to create Taylor1's
Taylor1(a::SubArray{T,1}, order::Int) where {T<:Number} = Taylor1(a.parent[a.indices...], order)
Taylor1(a::SubArray{T,1}) where {T<:Number} = Taylor1(a.parent[a.indices...])

# Shortcut to define Taylor1 independent variables
"""
    Taylor1([T::Type=Float64], order::Int)

Shortcut to define the independent variable of a `Taylor1{T}` polynomial of
given `order`. The default type for `T` is `Float64`.

```julia
julia> Taylor1(16)
 1.0 t + 𝒪(t¹⁷)

julia> Taylor1(Rational{Int}, 4)
 1//1 t + 𝒪(t⁵)
```
"""
Taylor1(::Type{T}, order::Int) where {T<:Number} = Taylor1( [zero(T), one(T)], order)
Taylor1(order::Int) = Taylor1(Float64, order)

######################### HomogeneousPolynomial
"""
    HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}

DataType for homogenous polynomials in many (>1) independent variables.

**Fields:**

- `coeffs  :: Array{T,1}` Expansion coefficients of the homogeneous
polynomial; the ``i``-th component is related to a monomial, where the degrees
of the independent variables are specified by `coeff_table[order+1][i]`.
- `order   :: Int` order (degree) of the homogenous polynomial.

Note that `HomogeneousPolynomial` variables are callable. For more information,
see [`evaluate`](@ref).
"""
struct HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}
    coeffs  :: Array{T,1}
    order   :: Int

    function HomogeneousPolynomial{T}(coeffs::Array{T,1}, order::Int) where T<:Number
        resize_coeffsHP!(coeffs, order)
        return new{T}(coeffs, order)
    end
end

HomogeneousPolynomial(x::HomogeneousPolynomial{T}) where {T<:Number} = x
HomogeneousPolynomial(coeffs::Array{T,1}, order::Int) where {T<:Number} =
    HomogeneousPolynomial{T}(coeffs, order)
HomogeneousPolynomial(coeffs::Array{T,1}) where {T<:Number} =
    HomogeneousPolynomial(coeffs, orderH(coeffs))
HomogeneousPolynomial(x::T, order::Int) where {T<:Number} =
    HomogeneousPolynomial([x], order)

# Shortcut to define HomogeneousPolynomial independent variable
"""
    HomogeneousPolynomial([T::Type=Float64], nv::Int])

Shortcut to define the `nv`-th independent `HomogeneousPolynomial{T}`.
The default type for `T` is `Float64`.

```julia
julia> HomogeneousPolynomial(1)
 1.0 x₁

julia> HomogeneousPolynomial(Rational{Int}, 2)
 1//1 x₂
```
"""
function HomogeneousPolynomial(::Type{T}, nv::Int) where {T<:Number}
    @assert 0 < nv ≤ get_numvars()
    v = zeros(T, get_numvars())
    @inbounds v[nv] = one(T)
    return HomogeneousPolynomial(v, 1)
end
HomogeneousPolynomial(nv::Int) = HomogeneousPolynomial(Float64, nv)

######################### TaylorN
"""
    TaylorN{T<:Number} <: AbstractSeries{T}

DataType for polynomial expansions in many (>1) independent variables.

**Fields:**

- `coeffs  :: Array{HomogeneousPolynomial{T},1}` Vector containing the
`HomogeneousPolynomial` entries. The ``i``-th component corresponds to the
homogeneous polynomial of degree ``i-1``.
- `order   :: Int`  maximum order of the polynomial expansion.

Note that `TaylorN` variables are callable. For more information, see
[`evaluate`](@ref).
"""
struct TaylorN{T<:Number} <: AbstractSeries{T}
    coeffs  :: Array{HomogeneousPolynomial{T},1}
    order   :: Int

    function TaylorN{T}(v::Array{HomogeneousPolynomial{T},1}, order::Int) where T<:Number
        coeffs = zeros(HomogeneousPolynomial{T}, order)
        @inbounds for i in eachindex(v)
            ord = v[i].order
            if ord ≤ order
                coeffs[ord+1] += v[i]
            end
        end
        new{T}(coeffs, order)
    end
end

TaylorN(x::TaylorN{T}) where {T<:Number} = x
function TaylorN(x::Array{HomogeneousPolynomial{T},1}, order::Int) where {T<:Number}
    if order == 0
        order = maxorderH(x)
    end
    return TaylorN{T}(x, order)
end
TaylorN(x::Array{HomogeneousPolynomial{T},1}) where {T<:Number} =
    TaylorN(x, maxorderH(x))
TaylorN(x::HomogeneousPolynomial{T}, order::Int) where {T<:Number} =
    TaylorN( [x], order )
TaylorN(x::HomogeneousPolynomial{T}) where {T<:Number} = TaylorN(x, x.order)
TaylorN(x::T, order::Int) where {T<:Number} =
    TaylorN(HomogeneousPolynomial([x], 0), order)

# Shortcut to define TaylorN independent variables
"""
    TaylorN([T::Type=Float64], nv::Int; [order::Int=get_order()])

Shortcut to define the `nv`-th independent `TaylorN{T}` variable as a
polynomial. The order is defined through the keyword parameter `order`,
whose default corresponds to `get_order()`. The default of type for
`T` is `Float64`.

```julia
julia> TaylorN(1)
 1.0 x₁ + 𝒪(‖x‖⁷)

julia> TaylorN(Rational{Int},2)
 1//1 x₂ + 𝒪(‖x‖⁷)
```
"""
TaylorN(::Type{T}, nv::Int; order::Int=get_order()) where {T<:Number} =
    TaylorN( HomogeneousPolynomial(T, nv), order )
TaylorN(nv::Int; order::Int=get_order()) = TaylorN(Float64, nv, order=order)

# A `Number` which is not an `AbstractSeries`
const NumberNotSeries = Union{Real,Complex}

# A `Number` which is not `TaylorN` nor a `HomogeneousPolynomial`
const NumberNotSeriesN = Union{Real,Complex,Taylor1}

## Additional Taylor1 outer constructor ##
Taylor1{T}(x::S) where {T<:Number,S<:NumberNotSeries} = Taylor1([convert(T,x)], 0)


text:
#=
This code provides the framework to stich together two separate tables
(either concatenating them horizontally or vertically).
=#

mutable struct IndexedTable{N, M} <: TexTable
    columns::Vector
    row_index::Index{N}
    col_index::Index{M}
end

IndexedTable(t::TableCol) = begin
    columns     = [t]
    row_index   = keys(t.data) |> collect |> sort
    col_index   = [t.header]
    return IndexedTable(columns, row_index, col_index)
end

convert(::Type{IndexedTable}, t::TexTable) = IndexedTable(t)
convert(::Type{IndexedTable}, t::IndexedTable) = t

################################################################################
#################### Merging and Concatenating #################################
################################################################################

function vcat(t1::IndexedTable, t2::IndexedTable)

    # Promote to the same dimensions
    t1, t2 = deepcopy.(promote(t1, t2))

    # Row Indices stay the same except within the highest group, where they need
    # to be shifted up in order to keep the index unique
    shift       =   maximum(get_idx(t1.row_index, 1)) -
                    minimum(get_idx(t2.row_index, 1)) + 1
    new_index   = map(t2.row_index) do idx
        i1 = idx.idx[1] + shift
        new_idx = tuple(i1, idx.idx[2:end]...)
        return update_index(idx, new_idx)
    end

    new_columns = deepcopy(t2.columns)
    for col in new_columns
        for (idx, new_idx) in zip(t2.row_index, new_index)
            if haskey(col.data, idx)
                col[new_idx] = pop!(col.data, idx)
            end
        end
    end

    row_index = vcat(t1.row_index, new_index)

    # Columns
    col_index = deepcopy(t1.col_index)
    columns   = deepcopy(t1.columns)
    for (i, idx) in enumerate(t2.col_index)

        # Figure out where to insert the column
        new_idx, s = insert_index!(col_index, idx)

        # It might be a new column
        if s > length(columns)
            push!(columns, new_columns[i])
        # If not, we need to move all the data over
        else
            for (key, value) in new_columns[i].data
                columns[s].data[key] = value
            end
        end
    end

    return IndexedTable(columns, row_index, col_index)

end

function hcat(t1::IndexedTable, t2::IndexedTable)

    # Promote to the same dimensions
    t1, t2 = deepcopy.(promote(t1, t2))

    # Column Indices stay the same except within the highest group,
    # where they need to be shifted up in order to keep the index unique
    shift       =   maximum(get_idx(t1.col_index, 1)) -
                    minimum(get_idx(t2.col_index, 1)) + 1
    new_index   = map(t2.col_index) do idx
        i1 = idx.idx[1] + shift
        return update_index(idx, tuple(i1, idx.idx[2:end]...))
    end
    col_index   = vcat(t1.col_index, new_index)

    # Row indices are merged in (inserted) one at a time, maintaining
    # strict insertion order in all index levels but the lowest one
    new_columns = deepcopy(t2.columns)
    row_index   = t1.row_index
    for idx in t2.row_index

        # Insert the index and recover the new_index and the required
        # insertion point
        new_idx, s = insert_index!(row_index, idx)

        # Rename the old indexes to the new ones
        for col in new_columns
            if haskey(col.data, idx)
                val          = pop!(col.data, idx)
                col[new_idx] = val
            end
        end
    end

    # Remap the internal column headers to keep them consistent
    old_new     = Dict(Pair.(t2.col_index, new_index))
    for col in new_columns
        col.header = old_new[col.header]
    end

    # Now, we're ready to append the columns together.
    columns     = vcat(t1.columns, new_columns)
    return IndexedTable(columns, row_index, col_index)
end

hcat(tables::Vararg{TexTable}) = reduce(hcat, tables)
vcat(tables::Vararg{TexTable}) = reduce(vcat, tables)

function hvcat(rows::Tuple{Vararg{Int}}, as::Vararg{TexTable})
    nbr = length(rows)  # number of block rows
    rs = Array{Any,1}(undef, nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(as[a:a-1+rows[i]]...)
        a += rows[i]
    end
    vcat(rs...)
end

# Make vcat and hcat work for all TexTables
vcat(t1::TexTable, t2::TexTable) = vcat(convert.(IndexedTable,
                                                 (t1, t2))...)
hcat(t1::TexTable, t2::TexTable) = hcat(convert.(IndexedTable,
                                                 (t1, t2))...)

join_table(t1::IndexedTable) = t1

function join_table(t1::TexTable, t2::TexTable)

    # Promote to the same dimensions
    t1, t2 = promote(convert.(IndexedTable, (t1, t2))...)

    t1_new = add_col_level(t1, 1)
    t2_new = add_col_level(t2, 2)

    return hcat(t1_new, t2_new)
end

function join_table(t1::TexTable, t2::TexTable,
                    t3::TexTable, args...)
    return join_table(join_table(t1,t2), t3, args...)
end

# Joining on Pairs
function join_table(p1::Pair{P1,T1}) where {P1 <: Printable,
                                      T1 <: TexTable}
    t1      = convert(IndexedTable, p1.second)
    t1_new  = add_col_level(p1.second, 1, p1.first)
end

function join_table(p1::Pair{P1,T1}, p2::Pair{P2,T2}) where
    {P1 <: Printable, P2 <: Printable, T1 <: TexTable, T2<:TexTable}

    t1, t2 = promote(convert.(IndexedTable, (p1.second, p2.second))...)

    t1_new = add_col_level(t1, 1, p1.first)
    t2_new = add_col_level(t2, 2, p2.first)

    return hcat(t1_new, t2_new)
end

function join_table(p1::Pair{P1,T1},
              p2::Pair{P2,T2},
              p3::Pair{P3,T3}, args...) where
                {P1 <: Printable, P2 <: Printable, P3<:Printable,
                T1 <: TexTable, T2<:TexTable, T3<:TexTable}

    return join_table(join_table(p1, p2), p3, args...)
end

join_table(t1::IndexedTable, p2::Pair{P2,T2}) where {P2, T2} = begin
    join_table(t1, join_table(p2))
end

join_table(p2::Pair{P2,T2},t1::IndexedTable) where {P2, T2} = begin
    join_table(join_table(p2), t1)
end

# Appending
append_table(t1::TexTable) = t1
function append_table(t1::TexTable, t2::TexTable)

    # Promote to the same dimensions
    t1, t2 = promote(t1, t2)

    t1_new = add_row_level(t1, 1)
    t2_new = add_row_level(t2, 2)

    return vcat(t1_new, t2_new)
end

function append_table(t1::TexTable, t2::TexTable, t3::TexTable, args...)
    return append_table(append_table(t1,t2), t3, args...)
end

# Appending on Pairs
function append_table(p1::Pair{P1,T1}) where {P1 <: Printable,
                                      T1 <: TexTable}
    t1      = convert(IndexedTable, p1.second)
    t1_new  = add_row_level(p1.second, 1, p1.first)
end

function append_table(p1::Pair{P1,T1}, p2::Pair{P2,T2}) where
    {P1 <: Printable, P2 <: Printable, T1 <: TexTable, T2<:TexTable}

    t1, t2 = promote(convert.(IndexedTable, (p1.second, p2.second))...)
    t1_new = add_row_level(t1, 1, p1.first)
    t2_new = add_row_level(t2, 2, p2.first)

    return vcat(t1_new, t2_new)
end

function append_table(p1::Pair{P1,T1},
              p2::Pair{P2,T2},
              p3::Pair{P3,T3}, args...) where
                {P1 <: Printable, P2 <: Printable, P3<:Printable,
                T1 <: TexTable, T2<:TexTable, T3<:TexTable}

    return append_table(append_table(p1, p2), p3, args...)
end

append_table(t1::IndexedTable, p2::Pair{P2,T2}) where {P2, T2} = begin
    append_table(t1, append_table(p2))
end

append_table(p2::Pair{P2,T2},t1::IndexedTable) where {P2, T2} = begin
    append_table(append_table(p2), t1)
end

append_table(t1::IndexedTable, p2::Pair{P2,T2}, args...) where {P2, T2} = begin
    append_table(append_table(t1, p2), args...)
end

################################################################################
#################### Conversion Between Dimensions #############################
################################################################################

function promote_rule(::Type{IndexedTable{N1,M1}},
                      ::Type{IndexedTable{N2,M2}}) where
                      {N1, M1, N2, M2}
    N = max(N1, N2)
    M = max(M1, M2)
    return IndexedTable{N, M}
end

function convert(::Type{IndexedTable{N, M}}, t::IndexedTable{N0, M0}) where
    {N,M,N0,M0}
    if (N0 > N) | (M0 > M)
        msg = """
        Cannot convert IndexedTable{$N0,$M0} to IndexedTable{$N,$M}
        """
        throw(error(msg))
    else
        for i=1:N-N0
            t = add_row_level(t, 1)
        end

        for i=1:M-M0
            t = add_col_level(t, 1)
        end
    end
    return t
end

function promote_rule(::Type{T1}, ::Type{T2}) where
    {T1 <: IndexedTable, T2 <: TexTable}
    return IndexedTable
end

################################################################################
#################### General Indexing ##########################################
################################################################################

function insert_index!(index::Index{N}, idx::TableIndex{N}) where N

    range = searchsorted(index, idx, lt=isless_group)

    # If it's empty, insert it in the right position
    if isempty(range)
        insert!(index, range.start, idx)
        return idx, range.start

    # Otherwise, check to see whether or not the last level matches already
    else

        N_index = get_idx(index[range], N)
        N_names = get_name(index[range], N)

        # If it does, then we don't have to do anything except check that the
        # strings are right
        if idx.name[N] in N_names
            loc = findall(N_names .== idx.name[N])[1]

            # Here's the new index
            new_idx = update_index(idx, tuple(idx.idx[1:N-1]..., loc))
            return new_idx, range.start + loc - 1
        else
            # Otherwise, it's not there so we need to insert it into the index,
            # and its last integer level should be one higher than all the
            # others
            new_idx = update_index(idx, tuple(idx.idx[1:N-1]...,
                                              maximum(N_index)+1))

            insert!(index, range.stop+1, new_idx)
            return new_idx, range.stop + 1
        end
    end
end

get_idx(index)               = map(x->x.idx, index)
get_idx(index, level::Int)   = map(x->x.idx[level], index)
get_name(index)              = map(x->x.name, index)
get_name(index, level::Int)  = map(x->x.name[level], index)

function find_level(index::Index{N}, idx::Idx{N}, level::Int) where N
    range = searchsorted(get_level(index, level), idx[level])
    return range
end

function add_level(index::Vector{TableIndex{N}}, level,
                   name::Printable="") where N
    return map(index) do idx
        return TableIndex(tuple(level, idx.idx...),
                          tuple(Symbol(name), idx.name...))
    end
end

"""
```
add_row_level(t::IndexedTable, level::Int, name::$Printable="")
```
Add's a new level to the row index with the given `level` for the integer
component of the index, and `name` for the symbol component
"""
function add_row_level(t::IndexedTable{N,M}, level::Int,
                       name::Printable="") where {N,M}

    new_rows = add_level(t.row_index, level, name)

    old_new  = Dict(Pair.(t.row_index, new_rows)...)

    new_columns = []
    for col in t.columns
        data = TableDict{N+1, FormattedNumber}()
        for (key, value) in col.data
            data[old_new[key]] = value
        end
        push!(new_columns, TableCol(col.header, data))
    end

    return IndexedTable(new_columns, new_rows, t.col_index)
end

"""
```
add_col_level(t::IndexedTable, level::Int, name::$Printable="")
```
Add's a new level to the column index with the given `level` for the integer
component of the index, and `name` for the symbol component
"""
function add_col_level(t::IndexedTable{N,M},
                       level::Int, name::Printable="") where {N,M}

    new_cols = add_level(t.col_index, level, name)
    old_new  = Dict(Pair.(t.col_index, new_cols))

    new_columns = []
    for col in t.columns
        push!(new_columns, TableCol(old_new[col.header],
                                    col.data))
    end

    return IndexedTable(new_columns, t.row_index, new_cols)
end

add_row_level(t::TexTable, args...) = add_row_level(IndexedTable(t), args...)
add_col_level(t::TexTable, args...) = add_col_level(IndexedTable(t), args...)

################################################################################
#################### Access Methods ############################################
################################################################################

Indexable{N}  = Union{TableIndex{N}, Tuple}
Indexable1D   = Union{Printable, Integer}

function row_loc(t::IndexedTable{N,M}, idx::Indexable{N}) where {N,M}
    locate(t.row_index, idx)
end

function col_loc(t::IndexedTable{N,M}, idx::Indexable{N}) where {N,M}
    locate(t.col_index, idx)
end

function loc(t::IndexedTable{N,M}, ridx::Indexable{N},
             cidx::Indexable{M}) where {N,M}

    rloc = locate(t.row_index, ridx)
    cloc = locate(t.col_index, cidx)

    if isempty(rloc) | isempty(cloc)
        throw(KeyError("key ($ridx, $cidx) not found"))
    elseif length(rloc) > 1
        throw(KeyError("$ridx does not uniquely identify a row"))
    elseif length(cloc) > 1
        throw(KeyError("$cidx does not uniquely identify a column"))
    else
        return rloc[1], cloc[1]
    end
end

function locate(index::Vector{TableIndex{N}}, idx::TableIndex{N}) where N
    return findall(index .== Ref(idx))
end

function locate(index::Vector{TableIndex{N}}, idx) where N
    length(idx) == N || throw(ArgumentError("$idx does not have dimension $N"))
    return findall(index) do x
        for i=1:N
            match_index(x, idx[i], i) || return false
        end
        return true
    end
end

function match_index(index::TableIndex{N}, idx::Printable, level::Int) where N
    return index.name[level] == Symbol(idx)
end

function match_index(index::TableIndex{N}, idx::Int, level::Int) where N
    return index.idx[level] == idx
end

function getindex(t::IndexedTable{N,M}, row::Indexable{N},
                  col::Indexable{M}) where {N,M}
    rloc, cloc = loc(t, row, col)
    return t.columns[cloc][t.row_index[rloc]]
end

function setindex!(t::IndexedTable, args...)
    throw(error("setindex! not implemented yet"))
end

# Fallback Methods
function getindex(t::IndexedTable, row::Indexable1D, col::Indexable1D)
    return t[tuple(row), tuple(col)]
end

function getindex(t::IndexedTable, row::Indexable, col::Indexable1D)
    return t[row, tuple(col)]
end

function getindex(t::IndexedTable, row::Indexable1D, col::Indexable)
    return t[tuple(row), col]
end

# Getvals
function get_vals(t::IndexedTable, row, col)
    get_vals(t[row, col])
end


text:
using Test
using Mods

@testset "Constructors" begin
    @test one(Mod{17}) == Mod{17}(1)
    @test oneunit(Mod{17}) == Mod{17}(1)
    @test zero(Mod{17}) == 0
    @test iszero(zero(Mod{17}))
    @test Mod{17}(1) == Mod{17}(1, 0)
    @test Mod{17}(1, 2) == 1 + 2im
    a = Mod{17}(3)
    @test typeof(a) == Mod{17,Int}
    a = Mod{17}(3, 2)
    @test typeof(a) == GaussMod{17,Int}
    a = zero(Mod{17})
    @test typeof(a) == Mod{17,Int}
    a = Mod{17}(1 // 2 + (3 // 4)im)
    @test typeof(a) == GaussMod{17,Int}
end

@testset "Mod arithmetic" begin
    p = 23
    a = Mod{p}(2)
    b = Mod{p}(25)
    @test a == b
    @test a == 2
    @test a == -21

    b = Mod{p}(20)
    @test a + b == 22
    @test a - b == -18
    @test a + a == 2a
    @test 0 - a == -a

    @test a * b == Mod{p}(17)
    @test (a / b) * b == a
    @test (b // a) * (2 // 1) == b
    @test a * (2 // 3) == (2a) * inv(Mod{p}(3))

    @test is_invertible(a)
    @test !is_invertible(Mod{10}(4))

    @test a^(p - 1) == 1
    @test a^(-1) == inv(a)

    @test Mod{p}(3 // 7) == 3 // 7
    @test isequal(3 // 7, Mod{p}(3 // 7))

    @test -Mod{13,Int}(typemin(Int)) == -Mod{13}(mod(typemin(Int), 13))
end

@testset "Mod arithmetic with UInt" begin
    p = UInt(23)
    a = Mod{p,UInt}(2)
    b = Mod{p,UInt}(25)
    @test a == b
    @test a == 2
    @test a == 25

    b = Mod{p,UInt}(20)
    @test a + b == 22
    @test a - b == 28
    @test a + a == 2a
    @test 0 - a == -a

    @test a * b == Mod{p,UInt}(17)
    @test (a / b) * b == a
    @test (b // a) * (2 // 1) == b
    @test a * (2 // 3) == (2a) * inv(Mod{p,UInt}(3))

    @test is_invertible(a)
    @test !is_invertible(Mod{10}(4))

    @test a^(p - 1) == 1
    @test a^(-1) == inv(a)
end

@testset "GaussMod arithmetic" begin
    @test one(GaussMod{6}) == Mod{6}(1 + 0im)
    @test zero(GaussMod{6}) == Mod{6}(0 + 0im)
    @test GaussMod{6}(2 + 3im) == Mod{6}(2 + 3im)
    @test GaussMod{6,Int}(2 + 3im) == Mod{6}(2 + 3im)
    @test rand(GaussMod{6}) isa GaussMod
    @test eltype(rand(GaussMod{6}, 4, 4)) <: GaussMod
    p = 23
    a = GaussMod{p}(3 - im)
    b = Mod{p}(5 + 5im)

    @test a + b == 8 + 4im
    @test a + Mod{p}(11) == Mod{p}(14, 22)
    @test -a == 20 + im
    @test a - b == Mod{p}(3 - im - 5 - 5im)

    @test a * b == Mod{p}((3 - im) * (5 + 5im))
    @test a / b == Mod{p}((3 - im) // (5 + 5im))

    @test a^(p * p - 1) == 1
    @test is_invertible(a)
    @test a * inv(a) == 1

    @test a / (1 + im) == a / Mod{p}(1 + im)
    @test imag(a * a') == 0
end

@testset "Large Modulus" begin

    p = 9223372036854775783   # This is a large prime
    x = Mod{p}(-2)
    @test x * x == 4
    @test x + x == -4
    @test x / x == 1
    @test x / 3 == x / Mod{p}(3)
    @test (x / 3) * (3 // x) == 1
    @test x // x == value(x) / x
    @test x^4 == 16
    @test x^(p - 1) == 1   # Fermat Little Theorem test
    @test 2x == x + x
    @test x - x == 0
    y = inv(x)
    @test x * y == 1
    @test x + p == x
    @test x * p == 0
    @test p - x == x - 2x

    p = 9223372036854775783   # This is a large prime
    x = Mod{p}(-2 + 0im)
    @test x * x == 4
    @test x + x == -4
    @test x / x == 1
    @test x / 3 == x / Mod{p}(3)
    @test (x / 3) * (3 // x) == 1
    @test x // x == value(x) / x
    @test x^4 == 16
    @test x^(p - 1) == 1   # Fermat Little Theorem test
    @test 2x == x + x
    @test x - x == 0
    y = inv(x)
    @test x * y == 1
    @test x + p == x
    @test x * p == 0
    @test p - x == x - 2x

    @test 0 <= value(rand(Mod{p})) < p

end

@testset "CRT" begin
    p = 23
    q = 91
    a = Mod{p}(17)
    b = Mod{q}(32)

    @test CRT(Int32, [], []) === Int32(0)
    x = CRT(a, b)
    @test typeof(x) <: BigInt
    @test a == mod(x, p)
    @test b == mod(x, q)

    c = Mod{101}(86)
    x = CRT(Int, a, b, c)
    @test typeof(x) <: Int

    @test a == mod(x, p)
    @test b == mod(x, q)
    @test c == mod(x, 101)

    @test CRT(
        BigInt,
        Mod{9223372036854775783}(9223372036854775782),
        Mod{9223372036854775643}(9223372036854775642),
    ) == 85070591730234614113402964855534653468
end

@testset "Matrices" begin
    M = zeros(Mod{11}, 3, 3)
    @test sum(M) == 0

    M = ones(Mod{11}, 5, 5)
    @test sum(M) == 3

    M = rand(Mod{11}, 5, 6)
    @test size(M) == (5, 6)

    A = rand(Mod{17}, 5, 5)
    X = values.(A)
    @test sum(X) == sum(Mod{17}.(A))
end

@testset "Hashing/Iterating" begin
    x = Mod{17}(11)
    y = x + 0im
    @test x == y
    @test hash(x) == hash(y)
    @test typeof(x) !== typeof(y)

    A = Set([x, y])
    @test length(A) == 1

    v = [Mod{10}(t) for t = 1:15]
    w = [Mod{10}(t + 0im) for t = 1:15]
    S = Set(v)
    T = Set(w)
    @test length(S) == 10
    @test S == T
    @test union(S, T) == intersect(S, T)
end

@testset "Iteration" begin
    @test sum(k for k in Mod{7}) == zero(Mod{7})
    @test collect(Mod{7}) == Mod{7}.(0:6)
    @test prod(k for k in Mod{7} if k != 0) == -1  # Wilson's theorem

    @test sum(GaussMod{7}) == 0
    @test length(GaussMod{5}) == 25
end

@testset "Linear system solution" begin
    A = Mod{13}.([1 2; 3 -4])
    b = Mod{13}.([3, 4])
    x = A \ b
    @test A * x == b
    @test inv(A) * b == x
end