text:
"""
    stepforward!(prob)

Step forward `prob` for one timestep.
"""
function stepforward!(prob::Problem)
  stepforward!(prob.state, prob.ts, prob.eqn, prob.vars, prob.params, prob.grid)
  prob.t = prob.state.t
  prob.step = prob.state.step
  nothing
end

"""
    stepforward!(prob, nsteps)

Step forward `prob` for `nsteps`.
"""
function stepforward!(prob::Problem, nsteps)
  for step = 1:nsteps
    stepforward!(prob)
  end
  nothing
end

"""
    stepforward!(prob, diags, nsteps)

Step forward `prob` for `nsteps`, incrementing diagnostics in the array `diags` along the way.
"""
function stepforward!(prob::Problem, diags::AbstractArray, nsteps)
  for step = 1:nsteps
    stepforward!(prob)
    for diag in diags
      if (prob.step+1) % diag.freq == 0
        increment!(diag)
      end
    end
  end
  nothing
end

"""
    getetdcoeffs(dt, LC; ncirc=32, rcirc=1)

Calculate ETDRK4 coefficients associated with the (diagonal) linear coefficient
LC by integrating over a small circle in complex space.
Note: arbitrary-precision arithmetic might provide a more robust method for
calculating these coefficients.
"""
function getetdcoeffs(dt, LC; ncirc=32, rcirc=1)

  shape = Tuple(cat(ncirc, ones(Int, ndims(LC)), dims=1))

  circ = zeros(cxeltype(LC), shape)
  circ .= rcirc * exp.(2π*im/ncirc*(0.5:1:(ncirc-0.5)))
  circ = permutedims(circ, ndims(circ):-1:1)

  zc = dt*LC .+ circ
  M = ndims(LC)+1

  # Four coefficients: ζ, α, β, Γ
  ζc = @.            ( exp(zc/2)-1 ) / zc
  αc = @. ( -4 - zc + exp(zc)*(4 - 3zc + zc^2) ) / zc^3
  βc = @.    ( 2  + zc + exp(zc)*(-2 + zc) ) / zc^3
  Γc = @. ( -4 - 3zc - zc^2 + exp(zc)*(4 - zc) ) / zc^3

  if eltype(LC) <: Real
    ζ = dt*real.(dropdims(mean(ζc, dims=M), dims=M))
    α = dt*real.(dropdims(mean(αc, dims=M), dims=M))
    β = dt*real.(dropdims(mean(βc, dims=M), dims=M))
    Γ = dt*real.(dropdims(mean(Γc, dims=M), dims=M))
  else
    ζ = dt*dropdims(mean(ζc, dims=M), dims=M)
    α = dt*dropdims(mean(αc, dims=M), dims=M)
    β = dt*dropdims(mean(βc, dims=M), dims=M)
    Γ = dt*dropdims(mean(Γc, dims=M), dims=M)
  end

  ζ, α, β, Γ
end

# ------------
# Timesteppers
# ------------

abstract type AbstractForwardEulerTimeStepper <: AbstractTimeStepper end
abstract type AbstractFilteredForwardEulerTimeStepper <: AbstractTimeStepper end

abstract type AbstractETDRK4TimeStepper <: AbstractTimeStepper end
abstract type AbstractFilteredETDRK4TimeStepper <: AbstractTimeStepper end

abstract type AbstractRK4TimeStepper <: AbstractTimeStepper end
abstract type AbstractFilteredRK4TimeStepper <: AbstractTimeStepper end

abstract type AbstractDualRK4TimeStepper <: AbstractTimeStepper end
abstract type AbstractDualFilteredRK4TimeStepper <: AbstractTimeStepper end

abstract type AbstractDualETDRK4TimeStepper <: AbstractTimeStepper end
abstract type AbstractDualFilteredETDRK4TimeStepper <: AbstractTimeStepper end

abstract type AbstractAB3TimeStepper <: AbstractTimeStepper end
abstract type AbstractFilteredAB3TimeStepper <: AbstractTimeStepper end

# -------------
# Forward Euler
# -------------

"""
    ForwardEulerTimeStepper(dt, LC)

Initialize a forward Euler timestepper. The forward Euler method is
the simplest time-stepping method in the books and is explicit and 1st-order accurate.
"""
struct ForwardEulerTimeStepper{T,dim} <: AbstractForwardEulerTimeStepper
  dt::Float64
  N::Array{T,dim}    # Explicit linear and nonlinear terms
  ForwardEulerTimeStepper(dt, N::Array{T,dim}) where {T,dim} = (
    new{cxeltype(N),ndims(N)}(dt, zeros(cxeltype(N), size(N))))
end

function stepforward!(s, ts::AbstractForwardEulerTimeStepper, eq, v, p, g)
  eq.calcN!(ts.N, s.sol, s.t, s, v, p, g)
  @. s.sol += ts.dt*(ts.N + eq.LC*s.sol)
  s.t += ts.dt
  s.step += 1
  nothing
end

# ----------------------
# Filtered Forward Euler
# ----------------------

"""
    FilteredForwardEulerTimeStepper(dt, LC, g; filterkwargs...)

Initialize a forward Euler timestepper with spectral filtering. The forward Euler method is
the simplest time-stepping method in the books and is explicit and 1st-order accurate.
"""
struct FilteredForwardEulerTimeStepper{dim} <: AbstractFilteredForwardEulerTimeStepper
  dt::Float64
  N::Array{Complex{Float64},dim}    # Explicit linear and nonlinear terms
  filter::Array{Float64,dim}        # Filter for solution
end

function FilteredForwardEulerTimeStepper(dt, LC, g; filterkwargs...)
  @createarrays eltype(LC) size(LC) N
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  FilteredForwardEulerTimeStepper{ndims(N)}(dt, N, filter)
end

function stepforward!(s, ts::AbstractFilteredForwardEulerTimeStepper, eq, v, p, g)
  eq.calcN!(ts.N, s.sol, s.t, s, v, p, g)
  @. s.sol = ts.filter*(s.sol + ts.dt*(ts.N + eq.LC*s.sol) )
  s.t += ts.dt
  s.step += 1
  nothing
end

# ------
# ETDRK4
# ------

"""
    ETDRK4TimeStepper(dt, LC)

The Rolls-Royce of time-stepping. Exact treatment of the implicit linear part
of the equation, explicit and 4th-order accurate integration of nonlinear
parts of equation.
"""
struct ETDRK4TimeStepper{T,dim} <: AbstractETDRK4TimeStepper
  dt::Float64
  LC::Array{T,dim}          # Linear coefficient
  # ETDRK4 coefficents
  ζ::Array{T,dim}
  α::Array{T,dim}
  β::Array{T,dim}
  Γ::Array{T,dim}
  expLCdt::Array{T,dim}     # Precomputed exp(LC*dt)
  expLCdt2::Array{T,dim}    # Precomputed exp(LC*dt/2)
  # Intermediate times, solutions, and nonlinear evaluations
  sol₁::Array{T,dim}
  sol₂::Array{T,dim}
  N₁::Array{T,dim}
  N₂::Array{T,dim}
  N₃::Array{T,dim}
  N₄::Array{T,dim}
end

function ETDRK4TimeStepper(dt, LC; cxsol=true)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  expLCdt  = exp.(dt*LC)
  expLCdt2 = exp.(0.5*dt*LC)
  ζ, α, β, Γ = getetdcoeffs(dt, LC)
  @createarrays T size(LC) sol₁ sol₂ N₁ N₂ N₃ N₄
  ETDRK4TimeStepper{T,ndims(LC)}(dt, LC, ζ, α, β, Γ, expLCdt, expLCdt2, sol₁, sol₂, N₁, N₂, N₃, N₄)
end

function stepforward!(s, ts::AbstractETDRK4TimeStepper, eq, v, p, g)
  # Substep 1
  eq.calcN!(ts.N₁, s.sol, s.t, s, v, p, g)
  @. ts.sol₁ = ts.expLCdt2*s.sol + ts.ζ*ts.N₁
  # Substep 2
  t2 = s.t + 0.5*ts.dt
  eq.calcN!(ts.N₂, ts.sol₁, t2, s, v, p, g)
  @. ts.sol₂ = ts.expLCdt2*s.sol + ts.ζ*ts.N₂
  # Substep 3
  eq.calcN!(ts.N₃, ts.sol₂, t2, s, v, p, g)
  @. ts.sol₂ = ts.expLCdt2*ts.sol₁ + ts.ζ*(2*ts.N₃ - ts.N₁)
  # Substep 4
  t3 = s.t + ts.dt
  eq.calcN!(ts.N₄, ts.sol₂, t3, s, v, p, g)

  # Update
  @. s.sol = (ts.expLCdt.*s.sol +   ts.α * ts.N₁
                                + 2*ts.β * (ts.N₂ + ts.N₃)
                                +   ts.Γ * ts.N₄ )
  s.t += ts.dt
  s.step += 1

  nothing
end

struct DualETDRK4TimeStepper{Tc,Tr,dimc,dimr} <: AbstractDualETDRK4TimeStepper
  dt::Float64
  c::ETDRK4TimeStepper{Tc,dimc}
  r::ETDRK4TimeStepper{Tr,dimr}
end

function DualETDRK4TimeStepper(dt, LCc, LCr)
  c = ETDRK4TimeStepper(dt, LCc)
  r = ETDRK4TimeStepper(dt, LCr)
  DualETDRK4TimeStepper{cxeltype(LCc),cxeltype(LCr),ndims(LCc),ndims(LCr)}(dt, c, r)
end
ETDRK4TimeStepper(dt, LCc, LCr) = DualETDRK4TimeStepper(dt, LCc, LCr)

function stepforward!(s::DualState, ts::AbstractDualETDRK4TimeStepper, eq, v, p, g)
  # Substep 1
  eq.calcN!(ts.c.N₁, ts.r.N₁, s.solc, s.solr, s.t, s, v, p, g)
  @. ts.c.sol₁ = ts.c.expLCdt2*s.solc + ts.c.ζ*ts.c.N₁
  @. ts.r.sol₁ = ts.r.expLCdt2*s.solr + ts.r.ζ*ts.r.N₁
  # Substep 2
  t2 = s.t + 0.5*ts.c.dt
  eq.calcN!(ts.c.N₂, ts.r.N₂, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.sol₂ = ts.c.expLCdt2*s.solc + ts.c.ζ*ts.c.N₂
  @. ts.r.sol₂ = ts.r.expLCdt2*s.solr + ts.r.ζ*ts.r.N₂
  # Substep 3
  eq.calcN!(ts.c.N₃, ts.r.N₃, ts.c.sol₂, ts.r.sol₂, t2, s, v, p, g)
  @. ts.c.sol₂ = (ts.c.expLCdt2*ts.c.sol₁ + ts.c.ζ*(2*ts.c.N₃ - ts.c.N₁))
  @. ts.r.sol₂ = (ts.r.expLCdt2*ts.r.sol₁ + ts.r.ζ*(2*ts.r.N₃ - ts.r.N₁))
  # Substep 4
  t3 = s.t + ts.c.dt
  eq.calcN!(ts.c.N₄, ts.r.N₄, ts.c.sol₂, ts.r.sol₂, t3, s, v, p, g)

  # Update
  @. s.solc = (ts.c.expLCdt.*s.solc +   ts.c.α * ts.c.N₁
                                    + 2*ts.c.β * (ts.c.N₂ + ts.c.N₃)
                                    +   ts.c.Γ * ts.c.N₄ )
  @. s.solr = (ts.r.expLCdt.*s.solr +   ts.r.α * ts.r.N₁
                                    + 2*ts.r.β * (ts.r.N₂ + ts.r.N₃)
                                    +   ts.r.Γ * ts.r.N₄ )
  s.t += ts.dt
  s.step += 1

  nothing
end

# ---------------
# Filtered ETDRK4
# ---------------

struct FilteredETDRK4TimeStepper{T,dim} <: AbstractFilteredETDRK4TimeStepper
  dt::Float64
  LC::Array{T,dim}          # Linear coefficient
  # ETDRK4 coefficents
  ζ::Array{T,dim}
  α::Array{T,dim}
  β::Array{T,dim}
  Γ::Array{T,dim}
  expLCdt::Array{T,dim}     # Precomputed exp(LC*dt)
  expLCdt2::Array{T,dim}    # Precomputed exp(LC*dt/2)
  # Intermediate times, solutions, and nonlinear evaluations
  sol₁::Array{T,dim}
  sol₂::Array{T,dim}
  N₁::Array{T,dim}
  N₂::Array{T,dim}
  N₃::Array{T,dim}
  N₄::Array{T,dim}
  filter::Array{T,dim}    # Filter for solution
end

function FilteredETDRK4TimeStepper(dt, LC, g; cxsol=true, filterkwargs...)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  expLCdt  = exp.(dt*LC)
  expLCdt2 = exp.(0.5*dt*LC)
  ζ, α, β, Γ = getetdcoeffs(dt, LC)
  @createarrays T size(LC) sol₁ sol₂ N₁ N₂ N₃ N₄
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  FilteredETDRK4TimeStepper{T,ndims(LC)}(dt, LC, ζ, α, β, Γ, expLCdt, expLCdt2, sol₁, sol₂, N₁, N₂, N₃, N₄, filter)
end

function stepforward!(s, ts::AbstractFilteredETDRK4TimeStepper, eq, v, p, g)
  # Substep 1
  eq.calcN!(ts.N₁, s.sol, s.t, s, v, p, g)
  @. ts.sol₁ = ts.expLCdt2*s.sol + ts.ζ*ts.N₁
  # Substep 2
  t2 = s.t + 0.5*ts.dt
  eq.calcN!(ts.N₂, ts.sol₁, t2, s, v, p, g)
  @. ts.sol₂ = ts.expLCdt2*s.sol + ts.ζ*ts.N₂
  # Substep 3
  eq.calcN!(ts.N₃, ts.sol₂, t2, s, v, p, g)
  @. ts.sol₂ = ts.expLCdt2*ts.sol₁ + ts.ζ*(2*ts.N₃ - ts.N₁)
  # Substep 4
  t3 = s.t + ts.dt
  eq.calcN!(ts.N₄, ts.sol₂, t3, s, v, p, g)

  # Update
  @. s.sol = ts.filter*(ts.expLCdt*s.sol +   ts.α * ts.N₁
                                         + 2*ts.β * (ts.N₂ + ts.N₃)
                                         +   ts.Γ * ts.N₄ )
  s.t += ts.dt
  s.step += 1

  nothing
end

struct DualFilteredETDRK4TimeStepper{Tc,Tr,dimc,dimr} <: AbstractDualFilteredETDRK4TimeStepper
  dt::Float64
  c::FilteredETDRK4TimeStepper{Tc,dimc}
  r::FilteredETDRK4TimeStepper{Tr,dimr}
end

function DualFilteredETDRK4TimeStepper(dt, LCc, LCr, g)
  c = FilteredETDRK4TimeStepper(dt, LCc, g)
  r = FilteredETDRK4TimeStepper(dt, LCr, g)
  DualFilteredETDRK4TimeStepper{cxeltype(LCc),cxeltype(LCr),ndims(LCc),ndims(LCr)}(dt, c, r)
end
FilteredETDRK4TimeStepper(dt, LCc, LCr, g) = DualFilteredETDRK4TimeStepper(dt, LCc, LCr, g)

function stepforward!(s::DualState, ts::AbstractDualFilteredETDRK4TimeStepper, eq, v, p, g)
  # Substep 1
  eq.calcN!(ts.c.N₁, ts.r.N₁, s.solc, s.solr, s.t, s, v, p, g)
  @. ts.c.sol₁ = ts.c.expLCdt2*s.solc + ts.c.ζ*ts.c.N₁
  @. ts.r.sol₁ = ts.r.expLCdt2*s.solr + ts.r.ζ*ts.r.N₁
  # Substep 2
  t2 = s.t + 0.5*ts.c.dt
  eq.calcN!(ts.c.N₂, ts.r.N₂, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.sol₂ = ts.c.expLCdt2*s.solc + ts.c.ζ*ts.c.N₂
  @. ts.r.sol₂ = ts.r.expLCdt2*s.solr + ts.r.ζ*ts.r.N₂
  # Substep 3
  eq.calcN!(ts.c.N₃, ts.r.N₃, ts.c.sol₂, ts.r.sol₂, t2, s, v, p, g)
  @. ts.c.sol₂ = (ts.c.expLCdt2*ts.c.sol₁ + ts.c.ζ*(2*ts.c.N₃ - ts.c.N₁))
  @. ts.r.sol₂ = (ts.r.expLCdt2*ts.r.sol₁ + ts.r.ζ*(2*ts.r.N₃ - ts.r.N₁))
  # Substep 4
  t3 = s.t + ts.c.dt
  eq.calcN!(ts.c.N₄, ts.r.N₄, ts.c.sol₂, ts.r.sol₂, t3, s, v, p, g)

  # Update
  @. s.solc = ts.c.filter*(ts.c.expLCdt.*s.solc +   ts.c.α * ts.c.N₁
                                                + 2*ts.c.β * (ts.c.N₂ + ts.c.N₃)
                                                +   ts.c.Γ * ts.c.N₄ )
  @. s.solr = ts.r.filter*(ts.r.expLCdt.*s.solr +   ts.r.α * ts.r.N₁
                                                + 2*ts.r.β * (ts.r.N₂ + ts.r.N₃)
                                                +   ts.r.Γ * ts.r.N₄ )
  s.t += ts.dt
  s.step += 1

  nothing
end

# ---
# RK4
# ---

"""
    RK4TimeStepper(dt, LC)

RK4 is the classical explicit 4th-order Runge-Kutta time-stepping
method. It uses a series of substeps/estimators to achieve 4th-order
accuracy over each individual time-step, at the cost of requiring
relatively more evaluations of the nonlinear right hand side.
It is described, among other places, in Bewley's Numerical
Renaissance.
"""
struct RK4TimeStepper{T,dim} <: AbstractRK4TimeStepper
  dt::Float64
  sol₁::Array{T,dim}
  RHS₁::Array{T,dim}
  RHS₂::Array{T,dim}
  RHS₃::Array{T,dim}
  RHS₄::Array{T,dim}
end

function RK4TimeStepper(dt, LC; cxsol=true)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  @createarrays T size(LC) sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  RK4TimeStepper{T,ndims(LC)}(dt, sol₁, RHS₁, RHS₂, RHS₃, RHS₄)
end

function stepforward!(s, ts::AbstractRK4TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS₁, s.sol, s.t, s, v, p, g)
  @. ts.RHS₁ += eq.LC*s.sol
  # Substep 1
  t2 = s.t + 0.5*ts.dt
  @. ts.sol₁ = s.sol + 0.5*ts.dt*ts.RHS₁
  eq.calcN!(ts.RHS₂, ts.sol₁, t2, s, v, p, g)
  @. ts.RHS₂ += eq.LC*ts.sol₁
  # Substep 2
  @. ts.sol₁ = s.sol + 0.5*ts.dt*ts.RHS₂
  eq.calcN!(ts.RHS₃, ts.sol₁, t2, s, v, p, g)
  @. ts.RHS₃ += eq.LC*ts.sol₁
  # Substep 3
  t3 = s.t + ts.dt
  @. ts.sol₁ = s.sol + ts.dt*ts.RHS₃
  eq.calcN!(ts.RHS₄, ts.sol₁, t3, s, v, p, g)
  @. ts.RHS₄ += eq.LC*ts.sol₁

  # Substep 4 and final step
  @. s.sol += ts.dt*(1/6*ts.RHS₁ + 1/3*ts.RHS₂ + 1/3*ts.RHS₃ + 1/6*ts.RHS₄)
  s.t += ts.dt
  s.step += 1
  nothing
end

struct DualRK4TimeStepper{Tc,Tr,dimc,dimr} <: AbstractDualRK4TimeStepper
  dt::Float64
  c::RK4TimeStepper{Tc,dimc}
  r::RK4TimeStepper{Tr,dimr}
end

function DualRK4TimeStepper(dt, LCc, LCr)
  c = RK4TimeStepper(dt, LCc)
  r = RK4TimeStepper(dt, LCr)
  DualRK4TimeStepper{cxeltype(LCc),cxeltype(LCr),ndims(LCc),ndims(LCr)}(dt, c, r)
end

RK4TimeStepper(dt, LCc, LCr) = DualRK4TimeStepper(dt, LCc, LCr)

function stepforward!(s, ts::AbstractDualRK4TimeStepper, eq, v, p, g)
  eq.calcN!(ts.c.RHS₁, ts.r.RHS₁, s.solc, s.solr, s.t, s, v, p, g)
  @. ts.c.RHS₁ += eq.LCc*s.solc
  @. ts.r.RHS₁ += eq.LCr*s.solr
  # Substep 1
  t2 = s.t + 0.5*ts.c.dt
  @. ts.c.sol₁ = s.solc + 0.5*ts.c.dt*ts.c.RHS₁
  @. ts.r.sol₁ = s.solr + 0.5*ts.r.dt*ts.r.RHS₁
  eq.calcN!(ts.c.RHS₂, ts.r.RHS₂, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.RHS₂ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₂ += eq.LCr*ts.r.sol₁
  # Substep 2
  @. ts.c.sol₁ = s.solc + 0.5*ts.c.dt*ts.c.RHS₂
  @. ts.r.sol₁ = s.solr + 0.5*ts.r.dt*ts.r.RHS₂
  eq.calcN!(ts.c.RHS₃, ts.r.RHS₃, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.RHS₃ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₃ += eq.LCr*ts.r.sol₁
  # Substep 3
  t3 = s.t + ts.c.dt
  @. ts.c.sol₁ = s.solc + ts.c.dt*ts.c.RHS₃
  @. ts.r.sol₁ = s.solr + ts.r.dt*ts.r.RHS₃
  eq.calcN!(ts.c.RHS₄, ts.r.RHS₄, ts.c.sol₁, ts.r.sol₁, t3, s, v, p, g)
  @. ts.c.RHS₄ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₄ += eq.LCr*ts.r.sol₁

  # Substep 4 and final step
  @. s.solc += ts.c.dt*(1/6*ts.c.RHS₁ + 1/3*ts.c.RHS₂ + 1/3*ts.c.RHS₃ + 1/6*ts.c.RHS₄)
  @. s.solr += ts.r.dt*(1/6*ts.r.RHS₁ + 1/3*ts.r.RHS₂ + 1/3*ts.r.RHS₃ + 1/6*ts.r.RHS₄)
  s.t += ts.c.dt
  s.step += 1
  nothing
end

# ------------
# Filtered RK4
# ------------

struct FilteredRK4TimeStepper{T,dim} <: AbstractFilteredRK4TimeStepper
  dt::Float64
  sol₁::Array{T,dim}
  RHS₁::Array{T,dim}
  RHS₂::Array{T,dim}
  RHS₃::Array{T,dim}
  RHS₄::Array{T,dim}
  filter::Array{T,dim}    # Filter for solution
end

function FilteredRK4TimeStepper(dt, LC, g; cxsol=true, filterkwargs...)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  @createarrays T size(LC) sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  FilteredRK4TimeStepper{T,ndims(LC)}(dt, sol₁, RHS₁, RHS₂, RHS₃, RHS₄, filter)
end

function stepforward!(s, ts::AbstractFilteredRK4TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS₁, s.sol, s.t, s, v, p, g)
  @. ts.RHS₁ += eq.LC*s.sol
  # Substep 1
  t2 = s.t + 0.5*ts.dt
  @. ts.sol₁ = s.sol + 0.5*ts.dt*ts.RHS₁
  eq.calcN!(ts.RHS₂, ts.sol₁, t2, s, v, p, g)
  @. ts.RHS₂ += eq.LC*ts.sol₁
  # Substep 2
  @. ts.sol₁ = s.sol + 0.5*ts.dt*ts.RHS₂
  eq.calcN!(ts.RHS₃, ts.sol₁, t2, s, v, p, g)
  @. ts.RHS₃ += eq.LC*ts.sol₁
  # Substep 3
  t3 = s.t + ts.dt
  @. ts.sol₁ = s.sol + ts.dt*ts.RHS₃
  eq.calcN!(ts.RHS₄, ts.sol₁, t3, s, v, p, g)
  @. ts.RHS₄ += eq.LC*ts.sol₁

  # Substep 4 and final step
  @. s.sol = ts.filter*(s.sol + ts.dt*(1/6*ts.RHS₁ + 1/3*ts.RHS₂ + 1/3*ts.RHS₃ + 1/6*ts.RHS₄))
  s.t += ts.dt
  s.step += 1

  nothing
end

struct DualFilteredRK4TimeStepper{Tc,Tr,dimc,dimr} <: AbstractDualFilteredRK4TimeStepper
  dt::Float64
  c::FilteredRK4TimeStepper{Tc,dimc}
  r::FilteredRK4TimeStepper{Tr,dimr}
end

function DualFilteredRK4TimeStepper(dt, LCc, LCr, g)
  c = FilteredRK4TimeStepper(dt, LCc, g)
  r = FilteredRK4TimeStepper(dt, LCr, g)
  DualFilteredRK4TimeStepper{cxeltype(LCc),cxeltype(LCr),ndims(LCc),ndims(LCr)}(dt, c, r)
end
FilteredRK4TimeStepper(dt, LCc, LCr, g) = DualFilteredRK4TimeStepper(dt, LCc, LCr, g)

function stepforward!(s, ts::AbstractDualFilteredRK4TimeStepper, eq, v, p, g)
  eq.calcN!(ts.c.RHS₁, ts.r.RHS₁, s.solc, s.solr, s.t, s, v, p, g)
  @. ts.c.RHS₁ += eq.LCc*s.solc
  @. ts.r.RHS₁ += eq.LCr*s.solr
  # Substep 1
  t2 = s.t + 0.5*ts.c.dt
  @. ts.c.sol₁ = s.solc + 0.5*ts.c.dt*ts.c.RHS₁
  @. ts.r.sol₁ = s.solr + 0.5*ts.r.dt*ts.r.RHS₁
  eq.calcN!(ts.c.RHS₂, ts.r.RHS₂, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.RHS₂ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₂ += eq.LCr*ts.r.sol₁
  # Substep 2
  @. ts.c.sol₁ = s.solc + 0.5*ts.c.dt*ts.c.RHS₂
  @. ts.r.sol₁ = s.solr + 0.5*ts.r.dt*ts.r.RHS₂
  eq.calcN!(ts.c.RHS₃, ts.r.RHS₃, ts.c.sol₁, ts.r.sol₁, t2, s, v, p, g)
  @. ts.c.RHS₃ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₃ += eq.LCr*ts.r.sol₁
  # Substep 3
  t3 = s.t + ts.c.dt
  @. ts.c.sol₁ = s.solc + ts.c.dt*ts.c.RHS₃
  @. ts.r.sol₁ = s.solr + ts.r.dt*ts.r.RHS₃
  eq.calcN!(ts.c.RHS₄, ts.r.RHS₄, ts.c.sol₁, ts.r.sol₁, t3, s, v, p, g)
  @. ts.c.RHS₄ += eq.LCc*ts.c.sol₁
  @. ts.r.RHS₄ += eq.LCr*ts.r.sol₁

  # Substep 4 and final step
  @. s.solc = ts.c.filter*(s.solc + ts.c.dt*(1/6*ts.c.RHS₁ + 1/3*ts.c.RHS₂ + 1/3*ts.c.RHS₃ + 1/6*ts.c.RHS₄))
  @. s.solr = ts.r.filter*(s.solr + ts.r.dt*(1/6*ts.r.RHS₁ + 1/3*ts.r.RHS₂ + 1/3*ts.r.RHS₃ + 1/6*ts.r.RHS₄))
  s.t += ts.c.dt
  s.step += 1
  nothing
end

# ---
# AB3
# ---

"""
    AB3TimeStepper(dt, LC)

3rd order Adams-Bashforth time stepping is an explicit scheme that uses
solutions from two previous time-steps to achieve 3rd order accuracy.
"""
struct AB3TimeStepper{T,dim} <: AbstractAB3TimeStepper
  dt::Float64
  RHS::Array{T,dim}
  RHS₋₁::Array{T,dim}
  RHS₋₂::Array{T,dim}
end

function AB3TimeStepper(dt, LC; cxsol=true)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  @createarrays T size(LC) RHS RHS₋₁ RHS₋₂
  AB3TimeStepper{T,ndims(LC)}(dt, RHS, RHS₋₁, RHS₋₂)
end

function stepforward!(s, ts::AbstractAB3TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS, s.sol, s.t, s, v, p, g)
  @. ts.RHS += eq.LC.*s.sol   # Add linear term to RHS

  if s.step < 3  # forward Euler steps to initialize AB3
    @. s.sol += ts.dt*ts.RHS    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. s.sol += ts.dt*(23/12*ts.RHS - 16/12*ts.RHS₋₁ + 5/12*ts.RHS₋₂)
  end

  s.t += ts.dt
  s.step += 1
  ts.RHS₋₂ .= ts.RHS₋₁          # Store
  ts.RHS₋₁ .= ts.RHS            # ... previous values of RHS
  nothing
end

# ------------
# Filtered AB3
# ------------

"""
    FilteredAB3TimeStepper(dt, LC, g; filterkwargs...)

3rd order Adams-Bashforth time stepping is an explicit scheme that uses
solutions from two previous time-steps to achieve 3rd order accuracy.
"""
struct FilteredAB3TimeStepper{T,dim} <: AbstractFilteredAB3TimeStepper
  dt::Float64
  RHS::Array{T,dim}
  RHS₋₁::Array{T,dim}
  RHS₋₂::Array{T,dim}
  filter::Array{T,dim}    # Filter for solution
end

function FilteredAB3TimeStepper(dt, LC, g; cxsol=true, filterkwargs...)
  T = cxsol ? cxeltype(LC) : eltype(LC)
  @createarrays T size(LC) RHS RHS₋₁ RHS₋₂
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  FilteredAB3TimeStepper{T,ndims(LC)}(dt, RHS, RHS₋₁, RHS₋₂, filter)
end

function stepforward!(s, ts::AbstractFilteredAB3TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS, s.sol, s.t, s, v, p, g)
  @. ts.RHS += eq.LC.*s.sol   # Add linear term to RHS

  if s.step < 3  # forward Euler steps to initialize AB3
    @. s.sol = ts.filter*(s.sol + ts.dt*ts.RHS)    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. s.sol = ts.filter*(s.sol + ts.dt*(23/12*ts.RHS - 16/12*ts.RHS₋₁ + 5/12*ts.RHS₋₂))
  end

  s.t += ts.dt
  s.step += 1
  ts.RHS₋₂ .= ts.RHS₋₁          # Store
  ts.RHS₋₁ .= ts.RHS            # ... previous values of RHS
  nothing
end


text:
using ImageMorphology
using ImageCore
using Test
using OffsetArrays
using ImageMetadata

@test isempty(detect_ambiguities(ImageMorphology))

using Documenter
Base.VERSION >= v"1.6" && doctest(ImageMorphology, manual = false)

@testset "ImageMorphology" begin
    include("convexhull.jl")
    include("connected.jl")
    include("dilation_and_erosion.jl")
    include("thinning.jl")
    include("imfill.jl")
    include("maxtree.jl")
    include("feature_transform.jl")
    include("clearborder.jl")
end

# multithreaded.jl runs on CI

ext: jl


text:
using GadgetIO, Test, DelimitedFiles, Downloads

function filter_dummy(filename::String)
    mtop = read_subfind(filename, "MTOP")
    return findall(mtop .> 7.0)
end

function pass_all(snap_file)
    h = read_header(snap_file)
    return collect(1:h.npart[1])
end

@info "downloading test data..."
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_sedov", "./snap_sedov")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/pos_sedov.dat", "./pos_sedov.dat")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/sub_002.0", "./sub_002.0")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/sub_002.1", "./sub_002.1")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/sub_002.2", "./sub_002.2")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/sub_002.3", "./sub_002.3")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.0", "./snap_002.0")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.1", "./snap_002.1")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.2", "./snap_002.2")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.3", "./snap_002.3")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.0.key", "./snap_002.0.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.1.key", "./snap_002.1.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.2.key", "./snap_002.2.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.3.key", "./snap_002.3.key")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_002.key.index", "./snap_002.key.index")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.0", "./snap_144.0")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.1", "./snap_144.1")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.2", "./snap_144.2")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.3", "./snap_144.3")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.0.key", "./snap_144.0.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.1.key", "./snap_144.1.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.2.key", "./snap_144.2.key")
Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.3.key", "./snap_144.3.key")

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.key.index", "./snap_144.key.index")

@info "done!"

@testset "GadgetIO" begin

    @testset "Objects" begin
        @test_nowarn SnapshotHeader()
        @test_nowarn InfoLine()
    end

    @testset "Read Snapshot" begin

        snap_file = "snap_sedov"

        @testset "Read blocks" begin
            @test_nowarn read_snap(snap_file)
            @test_nowarn read_header(snap_file)
            @test_nowarn read_info(snap_file)

            d = read_snap(snap_file, "POS", 0)

            ideal_file = joinpath(dirname(@__FILE__), "pos_sedov.dat")
            d_ideal = copy(transpose(Float32.(readdlm(ideal_file))))

            @test d == d_ideal

            # check if read to dict works
            @test_nowarn d = read_snap(snap_file)
            @test d["PartType0"]["POS"] == d_ideal

            snap_base = "snap_002"

            blocks = ["POS", "RHO"]
            data = read_blocks_over_all_files(snap_base, blocks, filter_function=pass_all, parttype=0)

            rho = read_block(snap_base, "RHO", parttype=0)

            @test data["RHO"] == rho

            pos = read_block(snap_base, "POS", parttype=0)

            @test data["POS"] == pos
        end

        @testset "Read particles in box" begin

            @testset "Brute-Force" begin
                center = Float32[3978.9688, -95.40625, -8845.25]
                rvir   = 118.76352
                pos = read_particles_in_volume("snap_002", "POS", center, rvir, use_keys=false, parttype=1)

                @test pos[:,1] ≈ Float32[3882.5537, -20.574343, -8768.669]

                pos = read_particles_in_box("snap_002", "POS", center .- rvir, center .+ rvir,
                                            use_keys=false, parttype=1, verbose=false)

                @test pos[:,1] ≈ Float32[3882.5537, -20.574343, -8768.669]
            end

            @testset "Peano-Hilbert reading" begin

                center = Float32[32572.607, 33825.95, 26314.158]
                rvir   = 1138.885

                pos = read_particles_in_volume("snap_144", "POS", center, rvir, use_keys=true, parttype=1)

                # check if correct number of particles were read
                @test size(pos,2) == 11981

                # check if the positions are the same as in Klaus' IDL code
                @test pos[:,1] ≈ Float32[ 31486.8, 33696.7, 25475.5 ]
                @test pos[:,2] ≈ Float32[ 32027.1, 33818.6, 25249.8 ]
                @test pos[:,3] ≈ Float32[ 31978.6, 33792.4, 25175.6 ]
            end

        end

        @testset "Read particles in geometry" begin
            center = [3978.9688, -95.40625, -8845.25]
            rvir   = 118.76352

            @testset "Cube" begin

                cube = GadgetCube(center .- rvir, center .+ rvir)
                pos  = read_particles_in_geometry("snap_002", "POS", cube, use_keys=false, parttype=1)

                @test pos["POS"][:,1] ≈ Float32[3882.5537, -20.574343, -8768.669]
            end

            @testset "Sphere" begin
                sphere = GadgetSphere(center, rvir)

                @test_nowarn read_particles_in_geometry("snap_002", "POS", sphere, use_keys=false, parttype=1)
            end

            @testset "Cylinder" begin
                cylinder = GadgetCylinder(center .- 0.5rvir, center .+ 0.5rvir,
                                        0.5rvir)

                @test_nowarn read_particles_in_geometry("snap_002", "POS", cylinder, use_keys=false, parttype=1)
            end

            # to do: use key files!

        end

        @testset "Read particles in halo" begin
            pos = read_particles_in_halo("snap_002", "POS", "sub_002", HaloID(0,4), use_keys=false)

            @test pos[:,1] ≈ Float32[3909.1545, -189.9392, -8845.135]

            ids = UInt32[0x000028fc, 0x00002594, 0x00002963, 0x00002681, 0x00001af4, 0x00001ff1, 0x000022d7, 0x00002267, 0x000029c0, 0x0000277b]
            pos = read_particles_by_id("snap_002", ids, "POS")

            @test pos ≈ copy(transpose(Float32[-692.6776 -5005.1025 1474.2584; -734.53326 -4894.864 1665.7646; -756.7661 -4985.657 1942.4185; -801.0376 -4920.4683 1884.446; -907.67645 -4945.71 1895.1641; -939.883 -4893.6753 1874.1469; -932.33136 -4891.3984 1109.0826; -819.5988 -5004.6147 1254.0176; -644.03674 -4939.248 1164.3943; -667.2112 -5048.75 995.14856]))
        end

        @testset "Read positions" begin
            center = Float32[3978.9688, -95.40625, -8845.25]
            rvir   = 118.76352

            # with filter function
            ff(filename) = filter_cube(filename, center .- rvir, center .+ rvir, parttype=1)
            read_positions = find_read_positions("snap_002", ff)

            @test read_positions["N_part"] == 87

            @test read_positions[0]["index"][1] == 2441
            @test read_positions[0]["index"][5] == 3966

            @test read_positions[0]["n_to_read"][1] == 2
            @test read_positions[0]["n_to_read"][4] == 74

            # with gadget geometry
            cube = GadgetCube(center .- rvir, center .+ rvir)
            read_positions_geo = find_read_positions("snap_002", cube, parttype=1)

            @test read_positions_geo["N_part"] == 87

            @test read_positions_geo[0]["index"][1] == 2441
            @test read_positions_geo[0]["index"][5] == 3966

            @test read_positions_geo[0]["n_to_read"][1] == 2
            @test read_positions_geo[0]["n_to_read"][4] == 74

            # test IO
            save_read_positions("dummy.bin", read_positions)
            loaded_read_positions = load_read_positions("dummy.bin")
            delete!(loaded_read_positions, "N_part")

            @test read_positions == loaded_read_positions

            # delete dummy file
            rm("dummy.bin")
        end

        @testset "Error Handling" begin
            @test_throws ErrorException("Please specify particle type!") read_block("snap_002.0", "POS")
            @test_throws ErrorException("Particle Type 5 not present in simulation!") read_block("snap_002.0", "POS", parttype=5, h=SnapshotHeader())
            @test_throws ErrorException("Block ABCD not present!") read_block("snap_002.0", "ABCD", parttype=0)
            @test_throws ErrorException("Requested block ABCD not present!") GadgetIO.check_block_position("snap_002.0", "ABCD")
            @test_throws ErrorException("Please provide either a dictionary with read positions or a filter function!") read_blocks_over_all_files("snap_002", ["POS"])
        end
    end

    @testset "Read subfind" begin

        subbase = "sub_002"
        subfile = "./sub_002.0"

        @testset "standard read" begin
            # check if standard reading works
            mtop = read_subfind("sub_002.0", "MTOP")
            @test mtop ≈ Float32[6.793532, 6.0309854, 7.230924, 8.415218]
        end

        @testset "Filter Subfind" begin
           # check if filter works
            find_mass_gt_7(M) = ( (M > 7.0) ? true : false )
            dummy = filter_subfind("sub_002", "MTOP", find_mass_gt_7)
            @test dummy[1] == HaloID(0, 3)
            @test dummy[2] == HaloID(0, 4)

            dummy2 = filter_subfind("sub_002", filter_dummy)
            @test dummy == dummy2

            # find the most massive halo in the sample subfind output
            center, rvir, haloid = find_most_massive_halo("sub_002", 4)
            @test center ≈ Float32[3978.9688, -95.40625, -8845.25]
            @test rvir   ≈ 118.76352
            @test haloid == HaloID(0, 4)
        end

        @testset "Read halo props" begin
            prop, haloid = read_halo_prop_and_id("sub_002", 4, "MTOP", 4)

            @test prop ≈ 5.431016
            @test haloid == HaloID(1, 1)

            prop2 = read_halo_prop("sub_002", haloid, "MTOP")

            @test prop2 == prop
        end

        @testset "Error handling" begin
            # check error handling
            @test_throws ErrorException("Block MVIR not present!") read_subfind(subfile, "MVIR")
            @test_throws ErrorException("Halo at index 1000 does not exist!") read_halo_prop_and_id(subfile, 1000, "MTOP", verbose=false)
        end

    end

    @testset "Snapshot utility" begin

        ref_file = joinpath(dirname(@__FILE__), "snap_sedov")

        @test_nowarn print_blocks(ref_file)

        # Check if a block is present
        present = block_present(ref_file, "POS")
        @test present == true
        present = block_present(ref_file, "BHMA")
        @test present == false

        blocks = ["POS", "VEL", "MASS"]
        blocks_checked, no_mass_block = GadgetIO.check_blocks(ref_file, blocks)

        @test blocks_checked == ["POS", "VEL"]
        @test no_mass_block == true

        @test_nowarn GadgetIO.get_block_positions(ref_file)

        filename = GadgetIO.select_file(ref_file, 0)

        @test filename == ref_file

        info = read_info(ref_file)
        d = GadgetIO.allocate_data_dict(["POS", "ID"], 10, info, true)
        @test haskey(d, "MASS")

        @test_throws ErrorException("File dummy_snap.0 not present!") GadgetIO.select_file("dummy_snap", 0)

        # shift across box border
        boxsize, boxsize_half = 10, 5
        @test GadgetIO.shift_across_box_border(1, 2, boxsize, boxsize_half) == 1
        @test GadgetIO.shift_across_box_border(8, 2, boxsize, boxsize_half) == -2
        @test GadgetIO.shift_across_box_border(2, 8, boxsize, boxsize_half) == 12
    end

    @testset "Peano-Hilbert" begin
        h_key = GadgetIO.read_keyheader("snap_002.0.key")

        @test h_key.nkeys_file[1] == Int32(597)
        @test h_key.bits == 10

        low_list, high_list, file_list = GadgetIO.read_key_index("snap_002.key.index")
        @test file_list == UInt32[0x00000000, 0x00000001, 0x00000002, 0x00000003, 0x00000000, 0x00000001, 0x00000002, 0x00000003, 0x00000000, 0x00000001, 0x00000002, 0x00000003, 0x00000000]

        @test GadgetIO.peano_hilbert_key(h_key.bits, 0, 0, 0) == 0
        @test GadgetIO.peano_hilbert_key(h_key.bits, 0, 0, 1) == 1
        @test GadgetIO.peano_hilbert_key(h_key.bits, 0, 1, 1) == 6

        @test GadgetIO.get_int_pos( 1000.5, h_key.domain_corners[1], h_key.domain_fac ) == 1

        @testset "Get Index Bounds" begin
            low_bounds, high_bounds = [0, 4, 7, 10], [1, 5, 8, 16]

            @test GadgetIO.get_index_bounds([1, 2, 3], low_bounds, high_bounds) == [1]
            @test GadgetIO.get_index_bounds([1, 2, 3, 10], low_bounds, high_bounds) == [1, 4]
            @test GadgetIO.get_index_bounds([1, 2, 10, 17], low_bounds, high_bounds) == [1, 4]
            @test GadgetIO.get_index_bounds([4, 5, 8, 12, 15], low_bounds, high_bounds) == [2, 3, 4]
        end

        keylist_ideal = UInt64[0x0000000000000000, 0x0000000000000001, 0x0000000000000002, 0x0000000000000003, 0x0000000000000004, 0x0000000000000005, 0x0000000000000006, 0x0000000000000007]

        h_key = GadgetIO.read_keyheader("snap_002.0.key")
        x0 = [-100.0, -100.0, -100.0]
        x1 = [1_000.0, 1_000.0, 1_000.0]

        keylist = GadgetIO.get_keylist(h_key, x0, x1)

        @test keylist == keylist_ideal
    end

    @testset "Write Snapshot" begin

        # read in reference file
        ref_file = joinpath(dirname(@__FILE__), "snap_sedov")
        head = head_to_obj(ref_file)
        x = read_snap(ref_file, "POS", 0)

        # specify output file for testing
        output_file = joinpath(dirname(@__FILE__), "write_test.dat")
        f = open(output_file, "w")

        @test_nowarn write_header(f, head)
        @test_nowarn write_block(f, x, "POS")
        close(f)

        pos_info = InfoLine("POS", Float32, 3, [1, 0, 0, 0, 0, 0])
        x_check = read_block(output_file, "POS", info=pos_info, parttype=0)

        # check if we read the same thing we wrote
        @test x_check == x

        # test info block write
        @testset "Write INFO block" begin
            ref_file = joinpath(dirname(@__FILE__), "snap_sedov")
            info_snap = read_info(ref_file)

            output_file = joinpath(dirname(@__FILE__), "write_test.dat")
            f = open(output_file, "w")
            write_info_block(f, info_snap)
            close(f)

            info_file = read_info(output_file)

            for i = 1:length(info_snap)

                @test info_file[i].block_name == info_snap[i].block_name
                @test info_file[i].data_type  == info_snap[i].data_type
                @test info_file[i].n_dim      == info_snap[i].n_dim
                @test info_file[i].is_present == info_snap[i].is_present

            end
        end
        # snap format 1
        f = open(output_file, "w")
        @test_nowarn write_block(f, x, "", snap_format=1)

        @test_throws ErrorException("Please specify blockname!") write_block(f, x, "")
        close(f)

        rm(output_file)
    end

    @testset "Find Index Locations" begin

        # create find and check arrays
        list_to_find = [2, 56, 354, 254, 653, 452, 7523, 45, 42, 742, 5423, 942, 105, 425, 815, 7821]
        list_to_check = [523, 9, 254, 653, 452, 2, 923, 815, 7821, 742, 354, 543, 942, 15, 25]

        # create sorted versions for forward search
        list_to_find_sorted = sort(list_to_find)
        list_to_check_sorted = sort(list_to_check)

        # correct indices of unsorted and sorted arrays
        indices = [3, 4, 5, 6, 8, 9, 10, 11, 13]
        indices_sorted = [1, 5, 6, 7, 10, 11, 12, 14, 15]

        # check methods
        @test issetequal(GadgetIO.get_index_list_arr(list_to_find_sorted, list_to_check_sorted), indices_sorted)
        @test issetequal(GadgetIO.get_index_list_dict(list_to_find, list_to_check), indices)
        @test issetequal(GadgetIO.get_index_list_dict(list_to_find_sorted, list_to_check_sorted), indices_sorted)
        @test issetequal(GadgetIO.get_index_list_set(list_to_find, list_to_check), indices)
        @test issetequal(GadgetIO.get_index_list_set(list_to_find_sorted, list_to_check_sorted), indices_sorted)

        # check that right methods are called
        @test get_index_list(list_to_check, list_to_find) == GadgetIO.get_index_list_set(list_to_check, list_to_find)
        @test get_index_list(list_to_check_sorted, list_to_find_sorted) == GadgetIO.get_index_list_arr(list_to_check_sorted, list_to_find_sorted)
    end

end

@info "delete test data..."
rm("snap_sedov")
rm("pos_sedov.dat")
rm("sub_002.0")
rm("sub_002.1")
rm("sub_002.2")
rm("sub_002.3")

rm("snap_002.0")
rm("snap_002.1")
rm("snap_002.2")
rm("snap_002.3")
rm("snap_002.0.key")
rm("snap_002.1.key")
rm("snap_002.2.key")
rm("snap_002.3.key")
rm("snap_002.key.index")

rm("snap_144.0")
rm("snap_144.1")
rm("snap_144.2")
rm("snap_144.3")
rm("snap_144.0.key")
rm("snap_144.1.key")
rm("snap_144.2.key")
rm("snap_144.3.key")
rm("snap_144.key.index")

@info "done!"


text:
macro except_str(expr, err_type)
    return quote
        let err = nothing
            try
                $(esc(expr))
            catch err
            end
            err === nothing && error("expected failure, but no exception thrown")
            @test typeof(err) === $(esc(err_type))
            buf = IOBuffer()
            showerror(buf, err)
            String(take!(buf))
        end
    end
end


text:
"""
Authors: Kara Ignatenko, 2021

Contains useful UI elements for Pluto.jl visualizations.

ui.jl
2021-09-22 v.0.1: First modular version
2021-09-24 v.0.2: Add support for arbitrary steps in the day picker
2021-10-07 v.0.3: Added the booleanButton widget
"""
module VizUI

using HypertextLiteral
using VegaLite

export dayPicker, vegaEmbed

"""
Create an interactive day picker widget for Pluto.jl.
The widget accepts integer values between `MIN` and `MAX`.
In addition to the textbox, buttons are provided to jump to the
first and last values, as well as an arbitrary number of steps
forwards or backwards as given in `steps`.
"""
function dayPicker(min::Int, max::Int, steps::Vector{Int}=[1])
    # create arrays of classes for backwards and forwards step buttons
    btnClasses = vcat([ ".prev-" * string(step) for step in reverse(steps) ],
                      [ ".next-" * string(step) for step in steps ])

    @htl("""

        <button class="first-btn"> first </button>
        $((
            @htl("""<button class="prev-$step"> -$step </button>""") for step in reverse(steps)
        ))
        <input/>
        $((
            @htl("""<button class="next-$step"> +$step </button>""") for step in steps
        ))
        <button class="last-btn"> last </button>

    """)
end

"""
Create the HTML for a Vega-Embed instance that will render
the VegaLite spec `spec`.
"""
function vegaEmbed(spec::VegaLite.VLSpec)
    # transform spec into JSON
    json = sprint(VegaLite.our_json_print, spec)

    # generate random number to allow for more than one visualization
    uid = abs(rand(Int))

    @htl("""

    """)

end

"""
Create a Pluto.jl button that will set a boolean value to `true`,
and then reset it back to false after `timeout` milliseconds.
"""
function booleanButton(text::String, timeout::Int=1000)
    @htl("""

        <button> $(text) </button>

    """)
end

end


text:
facts("Create MachineEnv structures") do
  vm = activate(RudeOilTests.machine)
  @fact length(vm.config) => greater_than(0)
  @fact vm.config => x -> ismatch(r"--tls", x)
end


text:
# same as `sum.jl`, but reusing the Julia compiler to compile and execute the IR

using Test

using LLVM
using LLVM.Interop

if length(ARGS) == 2
    x, y = parse.([Int32], ARGS[1:2])
else
    x = Int32(1)
    y = Int32(2)
end

JuliaContext() do ctx
    param_types = [LLVM.Int32Type(ctx), LLVM.Int32Type(ctx)]
    ret_type = LLVM.Int32Type(ctx)
    sum, _ = create_function(ret_type, param_types)

    # generate IR
    Builder(ctx) do builder
        entry = BasicBlock(sum, "entry", ctx)
        position!(builder, entry)

        tmp = add!(builder, parameters(sum)[1], parameters(sum)[2], "tmp")
        ret!(builder, tmp)
    end

    # make Julia compile and execute the function
    push!(function_attributes(sum), EnumAttribute("alwaysinline"))
    @eval call_sum(x, y) = $(call_function(sum, Int32, Tuple{Int32, Int32}, :(x, y)))
end

@test call_sum(x, y) == x + y


text:
export random, random_range, random_seed, set_random_seed

const random_seed = Parameter(12345)

set_random_seed(v::Int) = random_seed(v)

next_random(previous_random::Int) =
  let test = 16807*rem(previous_random,127773) - 2836*div(previous_random,127773)
    if test > 0
      if test > 2147483647
        test - 2147483647
      else
        test
      end
    else
      test + 2147483647
    end
  end

next_random!() = random_seed(next_random(random_seed()))

random(x::Int) = rem(next_random!(), x)

random(x::Real) = x*next_random!()/2147483647.0

random_range(x0, x1) =
  if x0 == x1
    x0
  else
    x0 + random(x1 - x0)
  end

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

export RGB, rgb, rgba, rgb_radiance
const rgb = RGB
const rgba = RGBA

rgb_radiance(c::RGB) = 0.265*red(c)+0.67*green(c)+0.065*blue(c)

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

required() = error("Required parameter")

#########################################
export division, map_division

division(t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    collect(include_last ? iter : take(iter, n))
  end

# Generic objects are processed using map_division
division(obj::Any, n::Real) = map_division(identity, obj, n)

map_division(f, t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    map(f, include_last ? iter : take(iter, n))
  end

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv),
               u0, u1, nu, include_last_u)

map_division(f, u0, u1, nu::Real, v0, v1, nv::Real, include_last_v::Bool=true) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu)

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real, include_last_v::Bool) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu, include_last_u)

########################################
# Grasshopper compatibility

export series, crossref, remap, cull

series(start::Real, step::Real, count::Int) =
  range(start, step=step, length=count)

export crossref_holistic
crossref_holistic(arr1, arr2) =
  vcat([arr1[i] for i in range(1, stop=length(arr1)) for j in arr2]...),
  vcat([arr2 for i in arr1]...)

crossref(as, bs) = [(a, b) for a in as, b in bs]

remap(in, (min_in, max_in), (min_out, max_out)) =
  min_out + (max_out-min_out)/(max_in-min_in)*(in-min_in)

cull(template, as) =
  [a for (a, t) in zip(as, cycle(template)) if t]

# To create paths from paths
export path_replace_suffix
path_replace_suffix(path::String, suffix::String) =
  let (base, old_suffix) = splitext(path)
    base * suffix
  end

#

replace_in(expr::Expr, replacements) =
    if expr.head == :.
        Expr(expr.head,
             replace_in(expr.args[1], replacements), expr.args[2])
    elseif expr.head == :quote
        expr
    else
        Expr(expr.head,
             map(arg -> replace_in(arg, replacements), expr.args) ...)
    end
replace_in(expr::Symbol, replacements) =
    get(replacements, expr, esc(expr))
replace_in(expr::Any, replacements) =
    expr

function process_named_params(def)
  call, body = def.args[1], def.args[2]
  name, params = call.args[1], call.args[2:end]
  idx = findfirst(p->p.head==:(::), params)
  mand, opts = isnothing(idx) ? ([], params) : (esc.(params[1:idx]), params[idx+1:end])
  opt_names = map(opt -> opt.args[1].args[1], opts)
  opt_types = map(opt -> esc(opt.args[1].args[2]), opts)
  opt_inits = map(opt -> opt.args[2], opts)
  opt_renames = map(Symbol ∘ string, opt_names)
  opt_replacements = Dict(zip(opt_names, opt_renames))
  mk_param(name,typ,init) = Expr(:kw, Expr(:(::), name, typ), init)
  opt_params = map(mk_param, opt_renames, opt_types, map(init -> replace_in(init, opt_replacements), opt_inits))
  key_params = map(mk_param, opt_names, opt_types, opt_renames)
  func_name = esc(name)
  :($(func_name)($(mand...), $(opt_params...); $(key_params...)) = $(esc(body)))
end

macro named_params(def)
  process_named_params(def)
end

export reverse_dict
reverse_dict(dict) =
  let rev_dict = Dict()
    for k in keys(dict)
      if dict[k] in keys(rev_dict)
         push!(rev_dict[dict[k]], k)
      else
         rev_dict[dict[k]] = [k]
      end
    end
  rev_dict
end

# To compute the sun's altitude and azimuth
export sun_pos
# Based on "How to compute planetary positions" by Paul Schlyter, Stockholm, Sweden
function sun_pos(year, month, day, hour, minute, Lstm, latitude, longitude)
  if abs(longitude-Lstm)>30
     @info("Longitude $(longitude) differs by more than 30 degrees from timezone meridian $(Lstm).")
  end
  # Calculate universal time (utime)
  utime = hour+(minute/60)-Lstm/15;
  if utime < 0
     day = day-1;
     utime = utime+24;
  end
  if utime > 24
     day = day+1;
     utime = utime-24;
  end
  int(x) = floor(Int, x)
  degrees(x) = 180*x/pi
  # Amount of days to, or from, the year 2000
  d = 367*year-int((7*int((year+int((month+9))/12)))/4)+int((275*month)/9)+day-730530+utime/24;
  # Longitude of perihelion (w), eccentricity (e)
  w = 282.9404+4.70935E-5*d;
  e = 0.016709-1.151E-9*d;
  mean_anomaly = 356.0470+0.9856002585*d;
  sun_longitude = w + mean_anomaly;
  # Obliquity of the ecliptic, eccentric anomaly (E)
  oblecl = 23.4393-3.563E-7*d;
  E = mean_anomaly+(180/pi)*e*sind(mean_anomaly)*(1+e*cosd(mean_anomaly));
  # Sun's rectangular coordinates in the plane of ecliptic (A,B)
  A = cosd(E)-e;
  B = sind(E)*sqrt(1-e*e);
  # Distance (r), longitude of the sun (lon)
  r = sqrt(A*A+B*B);
  true_anomaly = degrees(atan(B,A));
  lon = true_anomaly + w;
  # Calculate declination and right ascension
  decl = asin(sind(oblecl)*sind(lon));
  RA = degrees(atan(sind(lon)*cosd(oblecl),cosd(lon)))/15;
  # Greenwich meridian siderial time at 00:00 (GMST0),siderial time (SIDTIME), hour angle (HA)
  GMST0 = sun_longitude/15+12;
  SIDTIME = GMST0+utime+longitude/15;
  HA = (SIDTIME-RA)*15;
  # This is what we're looking for: Altitude & Azimuth
  Al = degrees(asin(sind(latitude)*sin(decl)+cosd(latitude)*cos(decl)*cosd(HA)));
  Az = degrees(atan(sind(HA),cosd(HA)*sind(latitude)-tan(decl)*cosd(latitude)))+180;
  Al, Az
end

sun_pos(date, timezone, latitude, longitude) =
  sun_pos(year(date), month(date), day(date), hour(date), minute(date), timezone, latitude, longitude)


text:
using SQLite
using NCDatasets
using DataStructures
using ProgressMeter
using Statistics

NC_FILE = "time_series_info.nc"
DB_FILE = "socat.sqlite"

"""
Create a netCDF containing a count of the unique days with measurements in
each grid cell from the SOCAT database

Run `make_socat_database.jl` first to extract data from the raw SOCAT data file
into a database for easier interrogation.
"""
function main()
  println("Initialising netCDF...")
  nc = initnetcdf(NC_FILE)

  println("Connecting to database...")
  db = SQLite.DB(DB_FILE)

  rowcount = 0
  countquery = DBInterface.execute(db, "SELECT COUNT(*) AS c FROM socat")
  for cr in countquery
    rowcount = cr[:c]
  end

  timeseriescount(nc, db, rowcount)
  measurementcount(nc, db, rowcount)
  meanR̅(nc, db, rowcount)

  close(nc)
  close(db)
end

"""
Initialise a netCDF file with lon/lat dimensions
"""
function initnetcdf(file::String)::NCDataset
  if ispath(file)
    rm(file)
  end

  nc = Dataset(file, "c")
  nc.attrib["title"] = "Summary of diel biases in the SOCAT dataset"

  defDim(nc, "longitude", 360)
  defDim(nc, "latitude", 180)

  lonvar = defVar(nc, "longitude", Float32, ["longitude"],
    attrib = OrderedDict("units" => "degrees_east"))

  lonvar[:] = collect(0.5:1:359.5)

  latvar = defVar(nc, "latitude", Float32, ["latitude"],
    attrib = OrderedDict("units" => "degrees_north"))

  latvar[:] = collect(-89.5:1:89.5)

  return nc
end

"""
Add a variable to a NetCDF file that has lon/lat dimensions

"""
function addvar(nc::NCDataset, name::String, type::Type)::NCDatasets.CFVariable
  return defVar(nc, name, type, ["longitude", "latitude"],
    attrib = OrderedDict("_FillValue" => convert(type, -999)))
end

"""
Create the time series count stats
"""
function timeseriescount(nc::NCDataset, db::SQLite.DB, rowcount::Int64)
  println("Days Per Cell")
  println("-------------")

  countvar = addvar(nc, "days_with_measurements", Int32)

  println("Searching...")
  rows = DBInterface.execute(db,
    "SELECT lonindex, latindex, year||dayofyear AS ydoy FROM socat " *
    "ORDER BY lonindex, latindex, ydoy"
  )

  prog = Progress(rowcount, "Reading:")

  currentlon = -999
  currentlat = -999
  seriescount = 0
  currentydoy = nothing

  for row in rows

    rowlon = row[:lonindex]
    rowlat = row[:latindex]
    rowydoy = row[:ydoy]

    if rowlon != currentlon || rowlat != currentlat
      if currentlon != -999
        countvar[currentlon, currentlat] = seriescount
      end

      currentlon = rowlon
      currentlat = rowlat
      seriescount = 0
      currentydoy = nothing
    end

    if rowydoy != currentydoy
      seriescount += 1
      currentydoy = rowydoy
    end

    next!(prog)
  end

  # Write last count
  countvar[currentlon, currentlat] = seriescount

  finish!(prog)

  return nothing
end

"""
Create the measurement count stats
"""
function measurementcount(nc::NCDataset, db::SQLite.DB, rowcount::Int64)
  println("Measurements Per Cell")
  println("---------------------")

  measurementcountvar = addvar(nc, "measurement_count", Int32)

  println("Searching...")
  rows = DBInterface.execute(db,
    "SELECT lonindex, latindex FROM socat ORDER BY lonindex, latindex"
  )

  prog = Progress(rowcount, "Reading:")

  currentlon = -999
  currentlat = -999
  measurements = 0

  for row in rows
    rowlon = row[:lonindex]
    rowlat = row[:latindex]

    if rowlon != currentlon || rowlat != currentlat
      if currentlon != -999
        measurementcountvar[currentlon, currentlat] = measurements
      end

      currentlon = rowlon
      currentlat = rowlat
      measurements = 0
    end

    measurements += 1

    next!(prog)
  end

  measurementcountvar[currentlon, currentlat] = measurements

  finish!(prog)
end

"""
Make the mean and R̅ stats
"""
function meanR̅(nc::NCDataset, db::SQLite.DB, rowcount::Int64)

  println("Minute Mean and R̅")
  println("-----------------")

  meanvar = addvar(nc, "minute_mean", Float64)
  R̅var = addvar(nc, "minute_R̅", Float64)

  println("Searching...")
  rows = DBInterface.execute(db,
    "SELECT lonindex, latindex, minuteofday FROM socat " *
    "ORDER BY lonindex, latindex"
  )

  prog = Progress(rowcount, "Reading:")

  currentlon = -999
  currentlat = -999
  cellminutes = Int64[]

  for row in rows

    rowlon = row[:lonindex]
    rowlat = row[:latindex]

    if rowlon != currentlon || rowlat != currentlat
      if currentlon != -999
        celldegrees = cellminutes .* 0.25
        meanvar[currentlon, currentlat] = meandegree(celldegrees) * 4
        R̅var[currentlon, currentlat] = R̅(celldegrees)
      end

      currentlon = rowlon
      currentlat = rowlat
      cellminutes = Int64[]
    end

    push!(cellminutes, row[:minuteofday])

    next!(prog)
  end

  # Write last count
  celldegrees = cellminutes .* 0.25
  meanvar[currentlon, currentlat] = meandegree(celldegrees) * 4
  R̅var[currentlon, currentlat] = R̅(celldegrees)

  finish!(prog)
end

"""
Calculate the mean minute of day from a series of minutes
"""
function meandegree(degrees::Vector{Float64})::Float64

  meandegree = rad2deg(atan(mean(sind.(degrees)),
    mean(cosd.(degrees))))

  if meandegree < 0
    meandegree = 360 - abs(meandegree)
  end

  return meandegree
end

"""
Calculate the standard deviation of a series of minutes
"""
function R̅(degrees::Vector{Float64})::Float64

  if length(degrees) == 1
    return 1
  else
    C = sum(cosd.(degrees)) / length(degrees)
    S = sum(sind.(degrees)) / length(degrees)
    return sqrt(C^2 + S^2)
  end
end

main()


text:
@testset "Interior Points" begin

    @testset "Optimal" begin

        A = [2 1 1 0; 1 2 0 1]
        b = [4, 4]
        c = [4, 3, 0, 0]
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        @test isapprox(output.x, [1.3333333333333335, 1.3333333333333333, 0.0, 0.0])
        @test isapprox(output.z, 9.333333333333334)
        @test output.termination_status == 1
    end

    @testset "Unbound" begin
        A = [-2 1 1 0; 1 -1 0 1]
        b = [2, 2]
        c = [1, 1, 0, 0]
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        @test isapprox(output.z, Inf)
        @test output.termination_status == 2
    end

    @testset "10 variables" begin
        Random.seed!(123)
        A_prime = rand(1:20,10,10)
        A = hcat(A_prime,Matrix(I,10,10))
        c = zeros(20)
        c[1:10] = rand(1:10,10)
        b = rand(1:30,10)
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        @test isapprox(output.x, [0.021665538253215974, 0.0, 0.07752200406228843, 0.0, 0.0, 0.0, 0.0, 0.04570074475287744, 0.0, 0.0, 8.20514556533514, 6.321259309410968, 21.437711577522006, 19.794515910629656, 20.410291130670277, 0.0, 13.41198375084631, 0.0, 0.3706838185511172, 0.0])
        @test isapprox(output.z, 1.052132701421801)
        @test output.termination_status == 1
    end

    @testset "100 variables" begin
        Random.seed!(123)
        A_prime = rand(1:200,100,100)
        A = hcat(A_prime,Matrix(I,100,100))
        c = zeros(200)
        c[1:100] = rand(1:10,100)
        b = rand(1:30,100)
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        # @test isapprox(output.x, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.005497498395887944, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.004566856299882416, 0.0, 0.0, 0.0034876174209097503, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0009317676782381164, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01570455795057045, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 17.40184009409551, 9.05933662252928, 2.4542994260563264, 17.887884835159845, 13.397392259170859, 9.497814553143392, 6.627345442541626, 24.510258929360113, 2.3023649010590277, 17.66006648141052, 6.994701120611391, 21.028703812751868, 4.793191254800283, 15.975326321340914, 15.130598356793328, 6.877495170470716, 7.635531583782399, 14.082784291797406, 28.537004213026215, 1.4083399576531552, 0.621640090069343, 5.852332912027959, 4.407094245045586, 0.0, 15.56328896075071, 1.2900329328036855, 4.539929122243899, 16.554462942593233, 21.631972955912353, 10.721262510221951, 22.19213851528068, 24.03888250539417, 22.857064994868264, 15.790811698611146, 27.1573999252441, 15.127882386214093, 20.345108081507963, 11.74172961111829, 1.4040609247162312, 19.348874070921234, 0.3098945926466158, 0.5185442843180255, 19.125269632747475, 2.112958913675148, 13.723058512535461, 2.890858614109204, 1.1990316170939093, 3.708410919794461, 6.438440493149732, 9.350162823884915, 19.16521657760707, 23.4226279552358, 21.831087891692743, 5.49088378229085, 21.080597467637197, 20.054570534101273, 12.289641088070667, 22.092125462963345, 0.0, 2.404603452588636, 2.2766266207447607, 7.4193128055816215, 3.952940853610769, 26.867919347410837, 11.395965262366085, 17.735575814931213, 19.070638090576377, 21.20402949545797, 7.293782612475125, 20.49986144945302, 17.341505314609986, 5.035809923511386, 22.35329451250835, 13.096585037952156, 0.0, 16.04799246352496, 6.240022717228514, 5.468877370887461, 15.15083430400292, 0.0, 25.053640364028315, 28.77180679303822, 27.109831323288216, 27.094220203936434, 25.565645772847475, 18.540856592257754, 4.006456577225859, 21.67052590950605, 24.85558428213155, 11.978660346389486, 20.4855200520266, 16.62394413711698, 0.0, 10.101361626820168, 14.425391420934153, 20.494717093391312, 20.146004124590927, 1.061291729338838, 9.893427272040896, 23.90678197299795];atol=1E-2)
        @test isapprox(output.z, 0.2832476238901244;atol=1E-2)
        @test output.termination_status == 1
    end

    @testset "Phase 1" begin
        A = [2 1 1 0 0; 1 2 0 1 0; -1 -1 0 0 1]
        b = [4, 4, -1]
        c = [4, 3, 0, 0, 0]
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        @test isapprox(output.x, [1.3333333333333333, 1.3333333333333335, 0.0, 0.0, 1.6666666666666667])
        @test isapprox(output.z, 9.333333333333332)
        @test output.termination_status == 1
    end

    @testset "Infeasible" begin
        A = [2 1 1 0 0; 1 2 0 1 0; -1 -1 0 0 1]
        b = [4, 4, -5]
        c = [4, 3, 0, 0, 0]
        input = create(A, b, c, tol=1E-10,solver=1,verbose=false)
        output = solve(input)

        @test output.termination_status == 3

    end

    @testset "Cicle" begin
        using LinearOptimizationSolver
        A = [1 2 1 0 0; -3 1 0 1 0; 1 -1 0 0 1]
        b = [14, 0, 2]
        c = [3, 4, 0, 0, 0]
        input = create(A, b, c, tol=1E-10, solver = 1,verbose=false)
        output = solve(input)

        @test isapprox(output.x, [6.0, 4.0, 0.0, 14.0, 0.0])
        @test isapprox(output.z, 34.0)
        @test output.termination_status == 1

    end
end


text:
# Model-space preconditioners for LSRTM
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: December 2016
#

export model_topmute, opTopmute, findWaterBottom, depth_scaling, opDepthScaling, laplace

function model_topmute(n::Tuple{Int64,Int64},mute_end::Array{Integer,1},length::Int64,x)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
	x = reshape(x,n)
	for j=1:n[1]
		mute_start = mute_end[j] - length
		filter = zeros(Float32,n[2])
		filter[1:mute_start-1] = 0f0
		filter[mute_end[j]+1:end] = 1f0
		taper_length = mute_end[j] - mute_start + 1
		taper = (1f0 + sin.((pi*(0:taper_length-1))/(taper_length - 1)-pi/2f0))/2f0
		filter[mute_start:mute_end[j]] = taper
		M = spdiagm(filter)
		x[j,:] = x[j,:].*filter
	end
	return vec(x)
end

function model_topmute(n::Tuple{Int64,Int64},mute_end::Int64,length::Int64,x)
# Model domain topmute for a velocity model of dimensions n = [nx, nz].
	x = reshape(x,n)
	mute_start = mute_end - length
	filter = zeros(Float32,n[2])
	filter[1:mute_start-1] = 0f0
	filter[mute_end+1:end] = 1f0
	taper_length = mute_end - mute_start + 1
	taper = (1f0 + sin.((pi*(0:taper_length-1))/(taper_length - 1)-pi/2f0))/2f0
	filter[mute_start:mute_end] = taper
	M = spdiagm(filter)
	for j=1:n[1]
		x[j,:] = x[j,:].*filter
	end
	return vec(x)
end

function opTopmute(n,mute_start,length)
	# JOLI wrapper for model domain topmute
	N = prod(n)
	T = joLinearFunctionFwdT(N,N,
 							 v -> model_topmute(n,mute_start,length,v),
 							 w -> model_topmute(n,mute_start,length,w),
							 Float32,Float32,name="Model topmute")
	return T
end

function findWaterBottom(m)
	#return the indices of the water bottom of a seismic image
	n = size(m)
	idx = zeros(Integer,n[1])
	eps = 1e-4
	for j=1:n[1]
		k=1
		while true
			if abs(m[j,k]) > eps
				idx[j] = k
				break
			end
			k += 1
		end
	end
	return idx
end

function depth_scaling(m,model)
# Linear depth scaling function for seismic images
	m = reshape(m,model.n)
	filter = sqrt(0f0:model.d[2]:model.d[2]*(model.n[2]-1))
	F = spdiagm(filter)
	for j=1:model.n[1]
		m[j,:] = F*m[j,:]
	end
	return vec(m)
end

function opDepthScaling(model)
# JOLI wrapper for the linear depth scaling function
	N = prod(model.n)
	D = joLinearFunctionFwdT(N, N,
							 v -> depth_scaling(v,model),
 							 w -> depth_scaling(w,model),
							 Float32,Float32,name="Depth scaling")
end

function laplace(model::TimeModeling.Model)
# 2D Laplace operator

	# 2nd derivative in x direction
	d1 = ones(Float32,model.n[1]-1)*1f0/(model.d[1]^2)
	d2 = ones(Float32,model.n[1])*-2f0/(model.d[1]^2)
	d3 = d1
	d = (d1, d2, d3)
	position = (-1, 0, 1)

	Dx = spdiagm(d,position,model.n[1],model.n[1])
	Ix = speye(Float32,model.n[1])

	# 2nd derivative in z direction
	d1 = ones(Float32,model.n[2]-1)*1f0/(model.d[2]^2)
	d2 = ones(Float32,model.n[2])*-2f0/(model.d[2]^2)
	d3 = d1
	d = (d1, d2, d3)
	position = (-1, 0, 1)

	Dz = spdiagm(d,position,model.n[2],model.n[2])
	Iz = speye(Float32,model.n[2])

	# 2D Laplace operator
	D = kron(Dz,Ix) + kron(Iz,Dx)
end

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
using Random
using RecipesBase
using LaTeXStrings
# restricted to 2-player game

"""
    `train_both!(player1, player2, N::Int)`

Find Nash equilibrium by training two players against each other

## Example
```
RPS = MatrixGame([
    (0,0) (-1,1) (1,-1);
    (1,-1) (0,0) (-1,1);
    (-1,1) (1,-1) (0,0)
])
init_strategy = [0.1,0.3,0.6]
player1 = MatrixPlayer(RPS, 1, copy(init_strategy))
player2 = MatrixPlayer(RPS, 2, copy(init_strategy))
train_both!(player1, player2, 1000)
```
"""
function train_both! end

"""
    `train_one!(player1, player2, N::Int)`

Train only player 1 while keeping strategy of player 2 constant.
Yields strategy for player 1 that is maximally exploitative of player 2.

```
RPS = MatrixGame([
    (0,0) (-1,1) (1,-1);
    (1,-1) (0,0) (-1,1);
    (-1,1) (1,-1) (0,0)
])
init_strategy = [0.1,0.3,0.6]
player1 = MatrixPlayer(RPS, 1, copy(init_strategy))
player2 = MatrixPlayer(RPS, 2, copy(init_strategy))
train_one!(player1, player2, 1000)
```
"""
function train_one! end

"""
Matrix Game
- Solved with regret matching
- Takes reward matrix as input

    `MatrixGame(R::Matrix{NTuple{2,T}})`

## Example
```
RPS = MatrixGame([
    (0,0) (-1,1) (1,-1);
    (1,-1) (0,0) (-1,1);
    (-1,1) (1,-1) (0,0)
])
```
"""
struct MatrixGame{T}
    R::Matrix{NTuple{2,T}}
end

"""
Matrix Game Player

Instantiate with `MatrixPlayer(game, id[, initial_strategy])`
Default to uniform strategy

## Example
```
init_strategy = [0.1,0.3,0.6]
player1 = MatrixPlayer(game, 1, init_strategy)
```
"""
struct MatrixPlayer{T}
    id::Int
    game::MatrixGame{T}
    strategy::Vector{Float64}
    hist::Vector{Vector{Float64}}
    regret_sum::Vector{Float64}
    strat_sum::Vector{Float64}
end

player(g::MatrixGame, id::Int) = MatrixPlayer(g, id)
player(g::MatrixGame, id::Int, s::Vector{Float64}) = MatrixPlayer(g, id, s)

function MatrixPlayer(game::MatrixGame, id::Int)
    n_actions = size(game.R, id)
    strategy = fill(1/n_actions, n_actions)
    MatrixPlayer(
        id,
        game,
        strategy,
        [copy(strategy)],
        zeros(n_actions),
        copy(strategy)
    )
end

function clear!(p::MatrixPlayer)
    resize!(p.hist, 1)
    p.regret_sum .= 0.0
    p.strat_sum .= p.strategy
    p.hist[1] .= p.strategy
end

function MatrixPlayer(game::MatrixGame, id::Int, strategy::Vector{Float64})
    n_actions = size(game.R, id)
    @assert length(strategy) == n_actions
    MatrixPlayer(
        id,
        game,
        copy(strategy),
        [copy(strategy)],
        zeros(n_actions),
        copy(strategy)
    )
end

function regret(game::MatrixGame, i::Int, a1::Int, a2::Int)
    u = game.R[a1,a2][i]
    if i === 1
        return [
            game.R[a,a2][i] - u
            for a in 1:size(game.R, i)
            ]
    elseif i === 2
        return [
            game.R[a1,a][i] - u
            for a in 1:size(game.R, i)
            ]
    else
        error("i must be 1 or 2")
    end
end

function gen_action(rng::AbstractRNG, p::MatrixPlayer)
    σ = p.strategy
    t = rand(rng)
    i = 1
    cw = σ[1]
    while cw < t && i < length(σ)
        i += 1
        @inbounds cw += σ[i]
    end
    return i
end

gen_action(p::MatrixPlayer) = gen_action(Random.GLOBAL_RNG, p)

function update_regret!(p::MatrixPlayer, a1::Int, a2::Int)
    p.regret_sum .+= regret(p.game, p.id, a1, a2)
end

function fill_normed_regret!(v::Vector{Float64}, r::Vector)
    s = 0.0
    for (i,k) in enumerate(r)
        if k > 0
            s += k
            v[i] = k
        else
            v[i] = 0.0
        end
    end
    if s == 0
        fill!(v, 1/length(v))
    else
        v ./= s
    end
end

function update_strategy!(p::MatrixPlayer)
    fill_normed_regret!(p.strategy, p.regret_sum)
    σ = p.strategy
    p.strat_sum .+= σ
    s = copy(p.strat_sum)
    s ./= sum(s)
    push!(p.hist, s)
end

function finalize_strategy!(p::MatrixPlayer)
    σ = p.strategy .= 0.0
    for σ_i in p.hist
        σ .+= σ_i
    end
    σ ./= sum(σ)
end

function evaluate(p1::MatrixPlayer, p2::MatrixPlayer)
    # strategies assumed already finalized
    game = p1.game
    σ1 = p1.strategy
    σ2 = p2.strategy
    R = game.R
    s1,s2 = size(R)
    p1_eval = 0.0
    p2_eval = 0.0
    for i in 1:s1, j in 1:s2
        prob = σ1[i]*σ2[j]
        p1_eval += prob*R[i,j][1]
        p2_eval += prob*R[i,j][2]
    end
    return p1_eval, p2_eval
end

function train_both!(p1::MatrixPlayer, p2::MatrixPlayer, N::Int; show_progress::Bool=false)
    prog = Progress(N; enabled=show_progress)
    for i in 1:N
        a1, a2 = gen_action(p1), gen_action(p2)

        update_regret!(p1, a1, a2)
        update_regret!(p2, a1, a2)

        update_strategy!(p1)
        update_strategy!(p2)

        next!(prog)
    end
    finalize_strategy!(p1)
    finalize_strategy!(p2)
end

function train_one!(p1::MatrixPlayer, p2::MatrixPlayer, N::Int; show_progress::Bool=false)
    prog = Progress(N; enabled=show_progress)
    for i in 1:N
        a1, a2 = gen_action(p1), gen_action(p2)

        update_regret!(p1, a1, a2)

        update_strategy!(p1)
        next!(prog)
    end
    finalize_strategy!(p1)
end

function cumulative_strategies(p::MatrixPlayer)
    mat = Matrix{Float64}(undef, length(p.hist), length(p.strategy))
    σ = zeros(Float64, length(p.strategy))

    for (i,σ_i) in enumerate(p.hist)
        σ = σ + (σ_i - σ)/i
        mat[i,:] .= σ
    end
    return mat
end

@recipe function f(p1::MatrixPlayer, p2::MatrixPlayer)
    layout --> 2
    link := :both
    framestyle := [:axes :axes]

    xlabel := "Training Steps"

    L1 = length(p1.strategy)
    labels1 = Matrix{String}(undef, 1, L1)
    for i in eachindex(labels1); labels1[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 1
        ylabel := "Strategy"
        title := "Player 1"
        labels := labels1
        reduce(hcat,p1.hist)'
    end

    L2 = length(p2.strategy)
    labels2 = Matrix{String}(undef, 1, L2)
    for i in eachindex(labels2); labels2[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 2
        title := "Player 2"
        labels := labels2
        reduce(hcat,p2.hist)'
    end
end

@recipe function f(p::MatrixPlayer)

    xlabel := "Training Steps"

    L = length(p.strategy)
    labels = Matrix{String}(undef, 1, L)
    for i in eachindex(labels); labels[i] = L"a_{%$(i)}"; end

    @series begin
        subplot := 1
        ylabel := "Strategy"
        title := "Player 1"
        labels := labels
        cumulative_strategies(p)
    end
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
include("plots.jl")

function get_max(coordinates, data)
    max, max_coor = findmax(data)
    return max, coordinates[:, max_coor]
end

function get_metrics(scan::ScanData, params::ScanParameters)
    voltage = scan.waveform
    pressure = voltage * params.calibration_factor

    isppa = calc_intensity_sppa(pressure, params.medium, params.excitation)
    ispta = calc_intensity_spta(pressure, params.medium, params.excitation)
    mi = calc_mechanical_index(pressure, params.excitation)

    isppa_max =  get_max(scan.coordinates, isppa.val)
    ispta_max = get_max(scan.coordinates, ispta.val)
    mi_max = get_max(scan.coordinates, mi.val)

    return ScanMetric(
        params,
        pressure,
        isppa,
        isppa_max,
        ispta,
        ispta_max,
        mi,
        mi_max,
     )
end

function compute_metrics(scan::ScanData, params::ScanParameters)
    return ScanData(
        scan.version,
        scan.axes,
        scan.scope_info,
        scan.time,
        scan.waveform,
        scan.coordinates,
        get_metrics(scan, params),
    )
end

function calc_intensity_sppa(pressure::PressureArray{N}, medium, excitation) where {N}
    return Metric{ISPPA, N-1}(IntensityMetrics.intensity_sppa(pressure, medium, excitation))
end

function calc_intensity_spta(pressure::PressureArray{N}, medium, excitation) where {N}
    return Metric{ISPTA, N-1}(IntensityMetrics.intensity_spta(pressure, medium, excitation))
end

function calc_mechanical_index(pressure::PressureArray{N}, excitation) where {N}
    return Metric{MI, N-1}(IntensityMetrics.mechanical_index(pressure, excitation))
end


text:
module SimpleED

export fermionicOperators, fermionicOperatorsFunc, fermionicHamiltonianBuilder, measure

function fermionicOperators(sites)

    f = [[1 0]; [0 -1.]]
    ii = [[1 0]; [0 1.]]
    c = [[0 0]; [1. 0]]
    a = [[0 1.]; [0 0]]
    n = c*a

    ci = zeros(Float64, (sites,2^sites,2^sites))
    ai = zeros(Float64, (sites,2^sites,2^sites))
    ni = zeros(Float64, (sites,2^sites,2^sites))
    fi = zeros(Float64, (sites,2^sites,2^sites))

    if sites == 1
        ci[1,:,:] = c
        ai[1,:,:] = a
        ni[1,:,:] = n
        fi[1,:,:] = f
    else
        for i in 1:sites
            ci[i,:,:] = kron((f for j in 1:i-1)...,c,(ii for j in i:sites-1)...)
            ai[i,:,:] = kron((f for j in 1:i-1)...,a,(ii for j in i:sites-1)...)
            ni[i,:,:] = kron((ii for j in 1:i-1)...,n,(ii for j in i:sites-1)...)
            fi[i,:,:] = kron((f for j in 1:i)...,(ii for j in i:sites-1)...)
        end
    end

    return ci,ai,ni,fi
    end;

function fermionicOperatorsFunc(sites)
    ci,ai,ni,fi = fermionicOperators(sites)
    return (x::Int -> ci[x,:,:], x::Int -> ai[x,:,:], x::Int -> ni[x,:,:], x::Int ->fi[x,:,:])
end

function fermionicHamiltonianBuilder(e,t,U)

    N = length(e)

    c,a,n,f = fermionicOperatorsFunc(N)

    H = zeros(typeof(e[1]*U[1][1]*t[1][1]),(2^N, 2^N))

    for i in 1:N
        H .+= e[i]*n(i)
        for j in 1:N
            H .+= t[i,j][1] * c(i) * a(j)
            H .+= U[i,j][1] * n(i) * n(j)
        end
    end

    return H

end;

function measure(
        observable,
        probabilities,
        ES
    )
    """"""

    obsTilde = transpose(ES.vectors)*observable*ES.vectors

    output ::typeof(probabilities[1]*observable[1]*ES.vectors[1]) = 0

    for (i,P) in enumerate(probabilities)
        output += P * obsTilde[i,i]
    end

    return output

end;

#function spfermionicOperators(sites)
#

#    f = sparse([[1 0]; [0 -1.]])
#    ii = sparse([[1 0]; [0 1.]])
#    c = sparse([[0 0]; [1. 0]])
#    a = sparse([[0 1.]; [0 0]])
#    n = c*a

#    ci = [kron((f for j in 1:i-1)...,c,(ii for j in i:sites-1)...) for i in 1:sites]
#    ai = [kron((f for j in 1:i-1)...,a,(ii for j in i:sites-1)...) for i in 1:sites]
#    ni = [kron((ii for j in 1:i-1)...,n,(ii for j in i:sites-1)...) for i in 1:sites]
#    fi = [kron((f for j in 1:i)...,(ii for j in i:sites-1)...) for i in 1:sites]

#    return ci,ai,ni,fi
#    end;

#function spHamiltonianBuilder(e,t,U)

#    N = length(e)

#    c,a,n,f = spfermionicOperators(N)

#    H = spzeros(Complex{Float64},2^N,2^N)

#    for i in 1:N
#        H .+= e[i]*n[i]
#        for j in 1:N
#            H .+= t[i,j][1] * c[i] * a[j]
#            H .+= U[i,j][1] * n[i] * n[j]
#        end
#    end

#    return H

#end;

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Variables (ModiaMath/Variables/_module.jl)
#

using Base.Meta:quot, isexpr

"""
    @signal signalName(arguments) begin ... end - Generate a function that instantiates a new signal

`signalName(arguments) begin ... end` is treated as the constructor function of a new
(mutable) struct that consists of all left-hand-side (scalar or vector) symbols
present between `begin ... end`.

# Examples
```julia
using ModiaMath

```
"""
macro signal(head, top_ex)
   # @signal Name(arguments) begin ... end
   #    is mapped to:
   # mutable struct Name <: ModiaMath.AbstractComponentWithVariables
   #    < declarations >
   #    function Name(<arguments>)
   #       < equations top_ex >
   #       _internal  = ModiaMath.ComponentInternal(:Name, nothing, computeVariables!)
   #       _signal = new(_internal, <declarations>)
   #       ModiaMath.initComponent(_internal, name1, :name1)
   #       ModiaMath.initComponent(_internal, name2, :name2)
   #         ...
   #       return _signal
   #    end
   # end

   if typeof(head) == Symbol
      signal_symbol = head
      fc = Expr(:call, head)
   else
      @assert isexpr(head, :call)
      signal_symbol = head.args[1]::Symbol
      fc = head
   end

   # Build:  internalComponentName = ModiaMath.Component(component_name)
   # internalComponentName = gensym()
   # push!(eq, :( a = ModiaMath.Component($component_name)) )

   # Inspect top_ex find
   #   symbol = < statements >
   # and store "symbol" in "names"
   this = gensym("this")
   eqInternal = []
   eq      = []
   push!(eqInternal,:( $this = new( ModiaMath.ComponentInternal($(quot(signal_symbol))) ) ))
   namesSet = DataStructures.OrderedSet{Symbol}() # Set of left-hand side symbols
   helpArray = []

   # Add keyword arguments to namesSet
   if length(fc.args) <= 1
      # No keyword arguments
   elseif typeof(fc.args[2]) == Symbol
      # Positional arguments present, but no keyword arguments
   elseif fc.args[2].head == :parameters
      # Keyword arguments present
      kwargs = fc.args[2].args
      # println("... keyword arguments = ", kwargs)
      for kw in kwargs
         # println("... kw = ", kw, ", kw.head = ", kw.head, ", kw.args = ", kw.args, ", name = ", kw.args[1])
         symbol = kw.args[1]
         if in(symbol, namesSet)
            error("\nfrom @signal ", signal_symbol, ": keyword argument",
                  "\n   ", kw,
                  "\nis defined twice.")
         end
         push!(namesSet, symbol)
         push!(helpArray, :($this.$symbol = $symbol))
      end
      # push!(fcp.args[2].args, Expr(:kw, :sceneOptions, :nothing))
   else
      error("Unknown error in signals.jl: typeof(fc.args[2]) = ", typeof(fc.args[2]))
   end

   for ex in top_ex.args
      # dump(ex)
      push!(eq, ex)
      if isexpr(ex, :(=), 2) && typeof(ex.args[1]) == Symbol
         symbol = ex.args[1]
         if in(symbol, namesSet)
            error("\nfrom @signal ", signal_symbol, ": Left hand side symbol in",
                  "\n   ", ex,
                  "\nis defined twice.")
         end
         push!(namesSet, symbol)
         push!(eq, :( ModiaMath.initComponent!($this,$symbol,$(quot(symbol)) ) ) )
     end
   end
   names = collect(namesSet)

   # Build:
   #   mutable struct signal_symbol <: ModiaMath.AbstractComponentWithVariables
   #      _internal::ModiaMath.ComponentInternal
   #      names[1]
   #      names[2]
   #       ...
   #      function signal_symbol(...)
   #         $(eq...)
   #         ModiaMath.initComponentInternal!( new(...) )
   #      end
   #   end
   code = :(mutable struct $signal_symbol <: Modia3D.AbstractSignal; end)
   # println("code.args = ", code.args)
   code.args[3].args = vcat( :( _internal::ModiaMath.ComponentInternal ), names )

   code2 = code.args[3].args
   push!(code2, :(function fc(); end) )
   i = length(code2)
   code2[i].args[1] = fc

   # Generate constructor function
   code2[i].args[2] = quote
      $(eqInternal...)
      $(helpArray...)
      $(eq...)
      return $this
   end

   # Print generated function (without "#....\n")
   # println( replace(sprint(showcompact,code) , r"# .*\n", "\n") )
   # println( sprint(showcompact,bcode) )

   return esc(code)
end


text:
using RandomFunctions
using Test

@testset "collatz.jl" begin
    @test collatz_steps(10) == [5, 16, 8, 4, 2, 1]
    @test collatz_steps_02(10) == [5, 16, 8, 4, 2, 1]
    @test max_stop_time(10) == (19,9)
    @test max_stop_time_02(10^7) == (8400511, 685)
    @test max_stop_time_03(10^7) == (8400511, 685)
end

@testset "edit_distance.jl" begin
    @test edit_distance("ATCTCGT", "ACTCCTC")[1] == 3.0
end

@testset "sieve.jl" begin
    @test sieve_01(7) == [2, 3, 5, 7]
    @test sieve_02(7) == [2, 3, 5, 7]
    @test sieve_03(7) == [2, 3, 5, 7]
end