text:
#= src/sketch.jl

References:

  N. Halko, P.G. Martinsson, J.A. Tropp. Finding structure with randomness:
    Probabilistic algorithms for constructing approximate matrix
    decompositions. SIAM Rev. 53 (2): 217-288, 2011.

  F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert. A fast randomized algorithm for
    the approximation of matrices. Appl. Comput. Harmon. Anal. 25: 335-366,
    2008.
=#

abstract type SketchMatrix end

size(A::SketchMatrix, dim::Integer...) = A.k

  function *(A::SketchMatrix, B::AbstractMatOrLinOp{T}) where T
    C = Array{T}(undef, A.k, size(B,2))
    mul!(C, A, B)
  end
  function *(A::SketchMatrix, Bc::Adjoint{T,<:AbstractMatOrLinOp{T}}) where T
    C = Array{T}(undef, A.k, size(Bc,2))
    mul!(C, A, Bc)
  end
  function *(A::AbstractMatOrLinOp{T}, B::SketchMatrix) where T
    C = Array{T}(undef, size(A,1), B.k)
    mul!(C, A, B)
  end
  function *(Ac::Adjoint{T,<:AbstractMatOrLinOp{T}}, B::SketchMatrix) where T
    C = Array{T}(undef, size(Ac,1), B.k)
    mul!(C, Ac, B)
  end

function sketch(
    side::Symbol, trans::Symbol, A::AbstractMatOrLinOp{T}, order::Integer,
    opts::LRAOptions=LRAOptions(T); args...) where T
  sketch_chkargs(side, trans, order)
  opts = copy(opts; args...)
  sketch_chkopts!(opts, A)
  if     opts.sketch == :randn  return sketch_randn(side, trans, A, order, opts)
  elseif opts.sketch == :sprn   return  sketch_sprn(side, trans, A, order, opts)
  elseif opts.sketch == :srft   return  sketch_srft(side, trans, A, order, opts)
  elseif opts.sketch == :sub    return   sketch_sub(side, trans, A, order, opts)
  end
end
sketch(side::Symbol, trans::Symbol, A, order::Integer, args...; kwargs...) =
  sketch(side, trans, LinOp(A), order, args...; kwargs...)
sketch(A, order::Integer, args...; kwargs...) =
  sketch(:left, :n, A, order, args...; kwargs...)

function sketchfact(
    side::Symbol, trans::Symbol, A::AbstractMatOrLinOp{T},
    opts::LRAOptions=LRAOptions(T); args...) where T
  sketchfact_chkargs(side, trans)
  opts = copy(opts; args...)
  sketch_chkopts!(opts, A)
  if     opts.sketch == :randn  return sketchfact_randn(side, trans, A, opts)
  elseif opts.sketch == :sprn   return  sketchfact_sprn(side, trans, A, opts)
  elseif opts.sketch == :srft   return  sketchfact_srft(side, trans, A, opts)
  elseif opts.sketch == :sub    return   sketchfact_sub(side, trans, A, opts)
  end
end
sketchfact(side::Symbol, trans::Symbol, A, args...; kwargs...) =
  sketchfact(side, trans, LinOp(A), args...; kwargs...)
sketchfact(A, args...; kwargs...) = sketchfact(:left, :n, A, args...; kwargs...)

function sketch_chkopts!(opts::LRAOptions, A)
  chkopts!(opts)
  if opts.sketch == :none
    @warn "invalid sketch method; using \"randn\""
    opts.sketch = :randn
  end
end

function sketch_chkargs(side::Symbol, trans::Symbol, order::Integer)
  sketchfact_chkargs(side, trans)
  order >= 0 || throw(ArgumentError("order"))
end
function sketchfact_chkargs(side::Symbol, trans::Symbol)
  side in (:left, :right) || throw(ArgumentError("side"))
  trans in (:n, :c) || throw(ArgumentError("trans"))
end

# RandomGaussian

mutable struct RandomGaussian <: SketchMatrix
  k::Int
end

function mul!(C, A::RandomGaussian, B::AbstractMatOrLinOp{T}) where T
  S = crandn(T, A.k, size(B,1))
  mul!(C, S, B)
end
function mul!(C, A::AbstractMatOrLinOp{T}, B::RandomGaussian) where T
  S = crandn(T, size(A,2), B.k)
  mul!(C, A, S)
end

  function mul!(C, A::RandomGaussian, Bc::Adjoint{T,<:AbstractMatOrLinOp{T}}) where T
    S = crandn(T, A.k, size(Bc,1))
    mul!(C, S, Bc)
  end

  function mul!(C, Ac::Adjoint{T,<:AbstractMatOrLinOp{T}}, B::RandomGaussian) where T
    S = crandn(T, size(Ac,2), B.k)
    mul!(C, Ac, S)
  end

## sketch interface

function sketch_randn(
    side::Symbol, trans::Symbol, A::AbstractMatOrLinOp, order::Integer,
    opts::LRAOptions)
  if side == :left
    if trans == :n  return sketch_randn_ln(A, order, opts)
    else            return sketch_randn_lc(A, order, opts)
    end
  else
    if trans == :n  return sketch_randn_rn(A, order, opts)
    else            return sketch_randn_rc(A, order, opts)
    end
  end
end

  function sketch_randn_ln(A::AbstractMatOrLinOp{T}, order::Integer, opts::LRAOptions) where T
    S = RandomGaussian(order)
    m, n = size(A)
    isherm = ishermitian(A)
    Bp = Array{T}(undef, order, n)
    if opts.sketch_randn_niter > 0
      Bq   = Array{T}(undef, order, m)
      tau  = Array{T}(undef, 1)
      work = Array{T}(undef, 1)
    end
    mul!(Bp, S, A)
    for i = 1:opts.sketch_randn_niter
      Bp, tau, work = orthrows!(Bp, tau, work, thin=false)
      mul!(Bq, Bp, A')
      if isherm
        Bp, Bq = Bq, Bp
      else
        Bq, tau, work = orthrows!(Bq, tau, work, thin=false)
        mul!(Bp, Bq, A)
      end
    end
    Bp
  end
  function sketch_randn_lc(A::AbstractMatOrLinOp{T}, order::Integer, opts::LRAOptions) where T
    S = RandomGaussian(order)
    m, n = size(A)
    isherm = ishermitian(A)
    Bp = Array{T}(undef, order, m)
    if opts.sketch_randn_niter > 0
      Bq   = Array{T}(undef, order, n)
      tau  = Array{T}(undef, 1)
      work = Array{T}(undef, 1)
    end
    mul!(Bp, S, A')
    for i = 1:opts.sketch_randn_niter
      Bp, tau, work = orthrows!(Bp, tau, work, thin=false)
      mul!(Bq, Bp, A)
      if isherm
        Bp, Bq = Bq, Bp
      else
        Bq, tau, work = orthrows!(Bq, tau, work, thin=false)
        mul!(Bp, Bq, A')
      end
    end
    Bp
  end
  function sketch_randn_rn(A::AbstractMatOrLinOp{T}, order::Integer, opts::LRAOptions) where T
    S = RandomGaussian(order)
    m, n = size(A)
    isherm = ishermitian(A)
    Bp = Array{T}(undef, m, order)
    if opts.sketch_randn_niter > 0
      Bq   = Array{T}(undef, n, order)
      tau  = Array{T}(undef, 1)
      work = Array{T}(undef, 1)
    end
    mul!(Bp, A, S)
    for i = 1:opts.sketch_randn_niter
      Bp, tau, work = orthcols!(Bp, tau, work, thin=false)
      mul!(Bq, A', Bp)
      if isherm
        Bp, Bq = Bq, Bp
      else
        Bq, tau, work = orthcols!(Bq, tau, work, thin=false)
        mul!(Bp, A, Bq)
      end
    end
    Bp
  end
  function sketch_randn_rc(A::AbstractMatOrLinOp{T}, order::Integer, opts::LRAOptions) where T
    S = RandomGaussian(order)
    m, n = size(A)
    isherm = ishermitian(A)
    Bp = Array{T}(undef, n, order)
    if opts.sketch_randn_niter > 0
      Bq   = Array{T}(undef, m, order)
      tau  = Array{T}(undef, 1)
      work = Array{T}(undef, 1)
    end
    mul!(Bp, A', S)
    for i = 1:opts.sketch_randn_niter
      Bp, tau, work = orthcols!(Bp, tau, work, thin=false)
      mul!(Bq, A, Bp)
      if isherm
        Bp, Bq = Bq, Bp
      else
        Bq, tau, work = orthcols!(Bq, tau, work, thin=false)
        mul!(Bp, A', Bq)
      end
    end
    Bp
  end

function sketchfact_randn(
    side::Symbol, trans::Symbol, A::AbstractMatOrLinOp, opts::LRAOptions)
  if opts.sketchfact_adap || opts.rank < 0
    n = opts.nb
    opts_ = copy(opts, maxdet_tol=-1.)
    while true
      order = opts.sketchfact_randn_samp(n)
      B = sketch_randn(side, trans, A, order, opts)
      p, tau, k = geqp3_adap!(B, opts_)
      k < n && return pqrback_postproc(B, p, tau, k, opts)
      n *= 2
    end
  else
    order = opts.sketchfact_randn_samp(opts.rank)
    B = sketch_randn(side, trans, A, order, opts)
    return pqrfact_backend!(B, opts)
  end
end

# RandomSubset

mutable struct RandomSubset <: SketchMatrix
  k::Int
end

function mul!(C, A::RandomSubset, B::AbstractMatrix)
  k = A.k
  m, n = size(B)
  size(C) == (k, n) || throw(DimensionMismatch)
  r = rand(1:m, k)
  for i = 1:k
    copyto!(view(C,[i],:), view(B,r[i:i],:))
  end
  C
end

function mul!(C, A::AbstractMatrix, B::RandomSubset)
  k = B.k
  m, n = size(A)
  size(C) == (m, k) || throw(DimensionMismatch)
  r = rand(1:n, k)
  for i = 1:k
    copyto!(view(C,:,[i]), view(A,:,r[i:i]))
  end
  C
end

  function mul!(C, A::RandomSubset, Bc::Adjoint{T,<:AbstractMatrix{T}}) where T
    B = parent(Bc)
    k = A.k
    m, n = size(B)
    size(C) == (k, m) || throw(DimensionMismatch)
    r = rand(1:n, k)
    for i = 1:k
      adjoint!(view(C,[i],:), view(B,:,r[i:i]))
    end
    C
  end

  function mul!(C, Ac::Adjoint{T,<:AbstractMatrix{T}}, B::RandomSubset) where T
    A = parent(Ac)
    k = B.k
    m, n = size(A)
    size(C) == (n, k) || throw(DimensionMismatch)
    r = rand(1:m, k)
    for i = 1:k
      adjoint!(view(C,:,[i]), view(A,r[i:i],:))
    end
    C
  end

## sketch interface

function sketch_sub(
    side::Symbol, trans::Symbol, A::AbstractMatrix, order::Integer,
    opts::LRAOptions)
  S = RandomSubset(order)
  if side == :left
    if trans == :n  return S*A
    else            return S*A'
    end
  else
    if trans == :n  return A *S
    else            return A'*S
    end
  end
end

function sketchfact_sub(
    side::Symbol, trans::Symbol, A::AbstractMatrix, opts::LRAOptions)
  if opts.sketchfact_adap || opts.rank < 0
    n = opts.nb
    opts_ = copy(opts, maxdet_tol=-1.)
    while true
      order = opts.sketchfact_sub_samp(n)
      B = sketch_sub(side, trans, A, order, opts)
      p, tau, k = geqp3_adap!(B, opts_)
      k < n && return pqrback_postproc(B, p, tau, k, opts)
      n *= 2
    end
  else
    order = opts.sketchfact_sub_samp(opts.rank)
    B = sketch_sub(side, trans, A, order, opts)
    return pqrfact_backend!(B, opts)
  end
end

# SRFT

mutable struct SRFT <: SketchMatrix
  k::Int
end

function srft_rand(::Type{T}, n::Integer) where T<:Real
  x = rand(n)
  @simd for i = 1:n
    @inbounds x[i] = 2*(x[i] > 0.5) - 1
  end
  x
end
function srft_rand(::Type{T}, n::Integer) where T<:Complex
  x = crandn(T, n)
  @inbounds for i = 1:n
    x[i] /= abs(x[i])
  end
  x
end

function srft_init(::Type{T}, n::Integer, k::Integer) where T<:Real
  l = k
  while n % l > 0
    l -= 1
  end
  m = div(n, l)
  X = Array{T}(undef, l, m)
  d = srft_rand(T, n)
  idx = rand(1:n, k)
  r2rplan! = plan_r2r!(X, R2HC, 1)
  X, d, idx, r2rplan!
end
function srft_init(::Type{T}, n::Integer, k::Integer) where T<:Complex
  l = k
  while n % l > 0
    l -= 1
  end
  m = div(n, l)
  X = Array{T}(undef, l, m)
  d = srft_rand(T, n)
  idx = rand(1:n, k)
  fftplan! = plan_fft!(X, 1)
  X, d, idx, fftplan!
end

function srft_reshape!(X::AbstractMatrix, d::AbstractVector, x::AbstractVecOrMat)
  l, m = size(X)
  i = 0
  @inbounds for j = 1:l, k = 1:m
    i += 1
    X[j,k] = d[i]*x[i]
  end
end
function srft_reshape_conj!(
    X::AbstractMatrix, d::AbstractVector, x::AbstractVecOrMat)
  l, m = size(X)
  i = 0
  @inbounds for j = 1:l, k = 1:m
    i += 1
    X[j,k] = d[i]*conj(x[i])
  end
end

function srft_apply!(
    y::AbstractVecOrMat{T}, X::AbstractMatrix{T}, idx::AbstractVector,
    r2rplan!::r2rFFTWPlan) where T<:Real
  l, m = size(X)
  n = l*m
  k = length(idx)
  mul!(X, r2rplan!, X)
  wn = exp(-2im*pi/n)
  wm = exp(-2im*pi/m)
  nnyq = div(n, 2)
  cnyq = div(l, 2)
  i = 1
  @inbounds while i <= k
    idx_ = idx[i] - 1
    row = fld(idx_, l) + 1
    col = rem(idx_, l) + 1
    w = wm^(row - 1)*wn^(col - 1)

    # find indices of real/imag parts
    col_ = col - 1
    cswap = col_ > cnyq
    ia = cswap ? l - col_ : col_
    ib = col_ > 0 ? l - ia : 0

    # initialze next entry to fill
    y[i] = 0
    s = one(T)

    # compute only one entry if purely real or no more space
    if in == 0 || in == nnyq || i == k
      for j = 1:m
        a = X[ia+1,j]
        b = ib == 0 || ib == ia ? zero(T) : X[ib+1,j]
        b = cswap ? -b : b
        y[i] += real(s*(a + b*im))
        s *= w
      end

    # else compute one entry each for real/imag parts
    else
      y[i+1] = 0
      for j = 1:m
        a = X[ia+1,j]
        b = ib == 0 || ib == ia ? zero(T) : X[ib+1,j]
        b = cswap ? -b : b
        z = s*(a + b*im)
        y[i  ] += real(z)
        y[i+1] += imag(z)
        s *= w
      end
      i += 1
    end
    i += 1
  end
end

function srft_apply!(
    y::AbstractVecOrMat{T}, X::AbstractMatrix, idx::AbstractVector,
    fftplan!::FFTWPlan) where T<:Complex
  l, m = size(X)
  n = l*m
  k = length(idx)
  mul!(X, fftplan!, X)
  wn = exp(-2im*pi/n)
  wm = exp(-2im*pi/m)
  @inbounds for i = 1:k
    row = fld(idx[i] - 1, l) + 1
    col = rem(idx[i] - 1, l) + 1
    w = wm^(row - 1)*wn^(col - 1)
    y[i] = 0
    s = one(T)
    for j = 1:m
      y[i] += s*X[col,j]
      s *= w
    end
  end
end

function mul!(C, A::SRFT, B::AbstractMatrix{T}) where T
  m, n = size(B)
  k = A.k
  size(C) == (k, n) || throw(DimensionMismatch)
  X, d, idx, fftplan! = srft_init(T, m, k)
  for i = 1:n
    srft_reshape!(X, d, view(B,:,i))
    srft_apply!(view(C,:,i), X, idx, fftplan!)
  end
  C
end

function mul!(C, A::AbstractMatrix{T}, B::SRFT) where T
  m, n = size(A)
  k = B.k
  size(C) == (m, k) || throw(DimensionMismatch)
  X, d, idx, fftplan! = srft_init(T, n, k)
  for i = 1:m
    srft_reshape!(X, d, view(A,i,:))
    srft_apply!(view(C,i,:), X, idx, fftplan!)
  end
  C
end

  function mul!(C, A::SRFT, Bc::Adjoint{T,<:AbstractMatrix{T}}) where T
    B = parent(Bc)
    m, n = size(B)
    k = A.k
    size(C) == (k, m) || throw(DimensionMismatch)
    X, d, idx, fftplan! = srft_init(T, n, k)
    for i = 1:m
      srft_reshape_conj!(X, d, view(B,i,:))
      srft_apply!(view(C,:,i), X, idx, fftplan!)
    end
    C
  end
  function mul!(C, Ac::Adjoint{T,<:AbstractMatrix{T}}, B::SRFT) where T
    A = parent(Ac)
    m, n = size(A)
    k = B.k
    size(C) == (n, k) || throw(DimensionMismatch)
    X, d, idx, fftplan! = srft_init(T, m, k)
    for i = 1:n
      srft_reshape_conj!(X, d, view(A,:,i))
      srft_apply!(view(C,i,:), X, idx, fftplan!)
    end
    C
  end

## sketch interface

function sketch_srft(
    side::Symbol, trans::Symbol, A::AbstractMatrix, order::Integer,
    opts::LRAOptions)
  S = SRFT(order)
  if side == :left
    if trans == :n  return S*A
    else            return S*A'
    end
  else
    if trans == :n  return A *S
    else            return A'*S
    end
  end
end

function sketchfact_srft(
    side::Symbol, trans::Symbol, A::AbstractMatrix, opts::LRAOptions)
  if opts.sketchfact_adap || opts.rank < 0
    n = opts.nb
    opts_ = copy(opts, maxdet_tol=-1.)
    while true
      order = opts.sketchfact_srft_samp(n)
      B = sketch_srft(side, trans, A, order, opts)
      p, tau, k = geqp3_adap!(B, opts_)
      k < n && return pqrback_postproc(B, p, tau, k, opts)
      n *= 2
    end
  else
    order = opts.sketchfact_srft_samp(opts.rank)
    B = sketch_srft(side, trans, A, order, opts)
    return pqrfact_backend!(B, opts)
  end
end

# SparseRandomGaussian

mutable struct SparseRandomGaussian <: SketchMatrix
  k::Int
end
const SparseRandGauss = SparseRandomGaussian

function mul!(C, A::SparseRandGauss, B::AbstractMatrix{T}) where T
  k = A.k
  m, n = size(B)
  size(C) == (k, n) || throw(DimensionMismatch)
  r = randperm(m)
  idx = 0
  @inbounds for i = 1:k
    p = fld(m - i, k) + 1
    s = crandn(T, p)
    for j = 1:n
      C[i,j] = 0
      for l = 1:p
        C[i,j] += s[l]*B[r[idx+l],j]
      end
    end
    idx += p
  end
  C
end

function mul!(C, A::AbstractMatrix{T}, B::SparseRandGauss) where T
  k = B.k
  m, n = size(A)
  size(C) == (m, k) || throw(DimensionMismatch)
  r = randperm(n)
  idx = 0
  @inbounds for j = 1:k
    p = fld(n - j, k) + 1
    s = crandn(T, p)
    for i = 1:m
      C[i,j] = 0
      for l = 1:p
        C[i,j] += A[i,r[idx+l]]*s[l]
      end
    end
    idx += p
  end
  C
end

  function mul!(C, A::SparseRandGauss, Bc::Adjoint{T,<:AbstractMatrix{T}}) where T
    B = parent(Bc)
    k = A.k
    m, n = size(B)
    size(C) == (k, m) || throw(DimensionMismatch)
    r = randperm(n)
    idx = 0
    @inbounds for i = 1:k
      p = fld(n - i, k) + 1
      s = crandn(T, p)
      for j = 1:m
        C[i,j] = 0
        for l = 1:p
          C[i,j] += s[l]*conj(B[j,r[idx+l]])
        end
      end
      idx += p
    end
    C
  end
  function mul!(C, Ac::Adjoint{T,<:AbstractMatrix{T}}, B::SparseRandGauss) where T
    A = parent(Ac)
    k = B.k
    m, n = size(A)
    size(C) == (n, k) || throw(DimensionMismatch)
    r = randperm(m)
    idx = 0
    @inbounds for j = 1:k
      p = fld(m - j, k) + 1
      s = crandn(T, p)
      for i = 1:n
        C[i,j] = 0
        for l = 1:p
          C[i,j] += conj(A[r[idx+l],i])*s[l]
        end
      end
      idx += p
    end
    C
  end

## sketch interface

function sketch_sprn(
    side::Symbol, trans::Symbol, A::AbstractMatrix, order::Integer,
    opts::LRAOptions)
  S = SparseRandGauss(order)
  if side == :left
    if trans == :n  return S*A
    else            return S*A'
    end
  else
    if trans == :n  return A *S
    else            return A'*S
    end
  end
end

function sketchfact_sprn(
    side::Symbol, trans::Symbol, A::AbstractMatrix, opts::LRAOptions)
  if opts.sketchfact_adap || opts.rank < 0
    n = opts.nb
    opts_ = copy(opts, maxdet_tol=-1.)
    while true
      B = sketch_sprn(side, trans, A, n, opts)
      p, tau, k = geqp3_adap!(B, opts_)
      k < n && return pqrback_postproc(B, p, tau, k, opts)
      n *= 2
    end
  else
    order = opts.rank
    B = sketch_sprn(side, trans, A, order, opts)
    return pqrfact_backend!(B, opts)
  end
end


text:
using TimetableSolver

C9 = Division(9, 3)
D10 = Division(10, 4)

subjectcounts_C9 = SubjectCounts(
    "MATH"=>2, "ENG"=>2, "SPORTS"=>1
)
subjectcounts_D10 = SubjectCounts(
    "MATH"=>2, "ENG"=>1, "SCIENCE"=>2
)

mark = Teacher("Mark Andrews",  "MARK",     ["MATH","SCIENCE"], [9, 10])
john = Teacher("John Curry",    "JOHN",     ["MATH"],           [9])
geeta = Teacher("Geeta Gupta",  "GGA",      ["ENG"],            [9, 10])
paul = Teacher("Paul Reid",     "PAUL",     ["SPORTS"],         [9])
jane = Teacher("Jane Murphy",   "JANE",     ["SCIENCE"],        [10])

tt_C9 = Timetable(
    [3, 2],
    subjectcounts_C9,
    [mark, john, geeta, paul],
    C9
)
tt_D10 = Timetable(
    [3, 2],
    subjectcounts_D10,
    [mark, geeta, jane],
    D10
)

schedule = Schedule(tt_C9, tt_D10)

@time solution, status = solve!(schedule)
schedule

ext: jl


text:
using DataDrivenDiffEq
using DataDrivenDiffEq.Optimize
using ModelingToolkit

using LinearAlgebra
using SafeTestsets
using Random

@info "Loading OrdinaryDiffEq"
using OrdinaryDiffEq

using Test
@info "Finished loading packages"

const GROUP = get(ENV, "GROUP", "All")

@time begin
    if GROUP == "All" || GROUP == "DataDrivenDiffEq" || GROUP == "Standard"
        @testset "Basis" begin include("./basis/basis.jl") end
        @testset "Basis Generators" begin include("./basis/generators.jl") end
        @testset "DataDrivenProblem" begin include("./problem.jl") end

        @testset "Sparse Identification" begin
            @testset "Pendulum" begin include("./sindy/pendulum.jl") end
            @testset "Michaelis Menten" begin include("./sindy/michaelis_menten.jl") end
            @testset "Cartpole" begin include("./sindy/cartpole.jl") end
        end

        @testset "Koopman" begin
            @testset "Linear Autonomous" begin include("./dmd/linear_autonomous.jl") end
            @testset "Linear Forced" begin include("./dmd/linear_forced.jl") end
            @testset "Nonlinear Autonomous" begin include("./dmd/nonlinear_autonomous.jl") end
            @testset "Nonlinear Forced" begin include("./dmd/nonlinear_forced.jl") end
        end

        include("./utils.jl")
    end

    if GROUP == "All" || GROUP == "Optional"

        @info "Loading Flux"
        using Flux
        @info "Loading Symbolic Regression"
        using SymbolicRegression

        @testset "Symbolic Regression" begin
            @testset "OccamNet" begin include("./symbolic_regression/occamnet.jl") end
            @testset "SymbolicRegression" begin include("./symbolic_regression/symbolic_regression.jl") end
        end
    end

    # These are excluded right now, until the deps are figured out
    #if GROUP == "Integration" || GROUP == "All"
    #    @safetestset "Partial Lotka Volterra Discovery " begin include("./applications/partial_lotka_volterra.jl") end
    #end
end


text:
# # Exploring model variants with `screen`

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# This notebooks demonstrates a simple use of [`screen`](@ref) to explore large
# number of model variants. On the toy *E. Coli* model, we try to map the
# impact of knocking out single reactions and 2-reaction combinations.

# First, let's download the data and load the packages and the model
!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json");

using COBREXA, Tulip

model = load_model(StandardModel, "e_coli_core.json")

# ## Preparing the functions
#
# While we could use the [`with_changed_bound`](@ref) to limit the reaction rates,
# but we will make a slightly more precise, usage-tailored modification. This
# is a straightforward modification of the [`with_changed_bound`](@ref) that does
# not set bounds "outside" of the original bounds:

with_limited_rate(reaction_id::String, limit) =
    model::StandardModel -> begin
        m = copy(model)
        m.reactions = copy(model.reactions)
        r = m.reactions[reaction_id] = copy(model.reactions[reaction_id])
        if -limit > r.lb
            r.lb = -limit
        end
        if limit < r.ub
            r.ub = limit
        end
        m
    end

# ## Knocking out single reactions
#
# This can be applied to see how much biomass can the model produce with
# certain reactions limited to almost zero:

get_biomass(x) = isnothing(x) ? 0 : x["BIOMASS_Ecoli_core_w_GAM"]

productions = screen_variants(
    model,
    [[with_limited_rate(rxn, 0.1)] for rxn in reactions(model)],
    model -> get_biomass(flux_balance_analysis_dict(model, Tulip.Optimizer)),
)

# This can be nicely plotted to give a more comprehensive overview of which
# reactions are critical or not:

using CairoMakie

disp_rxns = rand(1:95, 20) # only display 20 random fluxes to save space
barplot(
    productions[disp_rxns],
    direction = :x,
    axis = (
        yticks = (1:20, reactions(model)[disp_rxns]),
        xlabel = "Flux [mmol/gDW/h]",
        ylabel = "Reaction ID",
    ),
)

# ## Knocking out reaction combinations
#
# It is very easy to prepare matrices of biomass productions from all possible
# two-reaction knockouts. To make it more interesting, we will restrict one of
# the reactions of the pair a bit less, to see more possible outcomes.

# We do not process all reactions here to make the notebook rendering
# efficient, but you can easily remove the restriction, and speed the process
# up using parallel execution, by specifying `workers` argument (see
# documentation of [`screen`](@ref) for details)

rxns = reactions(model)

productions = screen_variants(
    model,
    [
        [with_limited_rate(rxn1, 3), with_limited_rate(rxn2, 0.1)] for rxn1 in rxns,
        rxn2 in rxns
    ],
    model -> get_biomass(flux_balance_analysis_dict(model, Tulip.Optimizer)),
)

f = Figure()
ax = Axis(
    f[1, 1],
    ylabel = "Reaction KO",
    xlabel = "Reaction KO",
    yticks = (1:20, reactions(model)[disp_rxns]),
    xticks = (1:20, reactions(model)[disp_rxns]),
    xticklabelrotation = -pi / 4,
)
hm = heatmap!(ax, productions[disp_rxns, disp_rxns])
Colorbar(f[1, 2], hm, label = "Flux [mmol/gDW/h]")
f


text:
immutable DatasTrends{norder} <: TimeIntergral_onestep
  n :: Int64
  function DatasTrends()
    @assert eltype(norder) <: Int64
    @assert norder >0
    new(norder)
  end
end
export DatasTrends

function eval!(tf :: DatasTrends{2}, p :: ModelParams,
  data :: Tuple, trends :: Tuple,
  t :: Float64, dt :: Float64)
  ## Trend on y(n)
  ytrend0 = fTrend(data[1],p)
  ## Trend on y(n-1)
  ytrendp = trends
  (y0,yp) = data
  yn = map((1:length(y0)...)) do ii
    correction1 = y0[ii] - yp[ii] - ytrend0[ii]*dt
    correction2 = dt*(ytrend0[ii] - ytrendp[ii])
    y0[ii] + ytrend0[ii]*dt - 5*correction1 + 2*correction2
  end
  return (p,(yn,y0),ytrend0)
end


text:
function out_neighbors(cg::CloudGraphs, nodeid::Int64)
  q = "match (n)-[:RELATION]-(neig) where id(n) = $(nodeid) return neig"

end

function out_neighbors(cg::CloudGraphs, nodesymb::ASCIIString)
  q = "match (n)-[:RELATION]-(neig) where n.symbol = $(nodesymb) return neig"

end


text:
module Chemostat_Human1

    import Chemostat
    import CSV
    import DataFrames: DataFrame

    include("Utils/Utils.jl")

    function __init__()
        _make_dirs()
        _load_NCI60_exchanges_data()
    end
end


text:
using AutonomousMerging
using Random
using AutomotiveDrivingModels
using POMDPs
using POMDPPolicies
using POMDPSimulators
using Test

function test_state(mdp::GenerativeMergingMDP, s,v=0.0, acc=0.0)
    ego = Vehicle(vehicle_state(35.0, merge_lane(mdp.env), 5.0, mdp.env.roadway), VehicleDef(), EGO_ID)
    veh1 = Vehicle(vehicle_state(s, main_lane(mdp.env), 4.9, mdp.env.roadway), VehicleDef(), EGO_ID + 1)
    scene = Scene()
    push!(scene, ego)
    push!(scene, veh1)
    return AugScene(scene, (acc=acc,))
end

@testset "Environment" begin
    env = MergingEnvironment(main_lane_angle = 0.0, merge_lane_angle = pi/12)
    main = main_lane(env)
    merge = merge_lane(env)
    @test main == env.roadway[LaneTag(MAIN_LANE_ID, 1)]
    @test main.curve[end].s == env.main_lane_length + env.after_merge_length
    @test merge.curve[end].s == env.merge_lane_length
end

@testset "GenerativeMDP" begin

    rng = MersenneTwister(1)
    mdp = GenerativeMergingMDP(max_cars=12,
                               min_cars=10,
                               driver_type = :random,
                               observe_cooperation = true,
                               initial_ego_velocity=0.0)
    s0 = initialstate(mdp, rng)
    policy = RandomPolicy(mdp, rng=rng)

    hr = HistoryRecorder(rng = rng, max_steps=100)
    hist = POMDPSimulators.simulate(hr, mdp, policy, s0)

    s = s0
    svec = convert_s(Vector{Float64}, s, mdp)
    srec = convert_s(AugScene, svec, mdp)
    @test get_by_id(srec.scene, EGO_ID).state.posG ≈ get_by_id(s.scene, EGO_ID).state.posG

    s = state_hist(hist)[end]
    svec = convert_s(Vector{Float64}, s, mdp)
    srec = convert_s(AugScene, svec, mdp)
    @test get_by_id(srec.scene, EGO_ID).state.posG ≈ get_by_id(s.scene, EGO_ID).state.posG
end

@testset "CooperativeIDM" begin
    rng = MersenneTwister(1)

    mdp = GenerativeMergingMDP(random_n_cars=true, dt=0.5)

    mdp.driver_models[2] = CooperativeIDM(c=1.0)
    set_desired_speed!(mdp.driver_models[2], 5.0)

    s0 = test_state(mdp, 85.0)
    policy = FunctionPolicy(s->7)
    hr = HistoryRecorder(rng = rng, max_steps=40)
    hist = POMDPSimulators.simulate(hr, mdp, policy, s0)

end


text:
module Writer

export write_pddl, write_domain, write_problem
export save_domain, save_problem

using Julog
using ..PDDL:
    IMPLIED_REQUIREMENTS, Domain, Problem, Action,
    get_name, get_requirements, get_typetree, get_constants, get_constypes,
    get_predicates, get_functions, get_actions, get_axioms,
    get_domain_name, get_objects, get_objtypes,
    get_init_terms, get_goal, get_metric, get_constraints,
    get_argvars, get_argtypes, get_precond, get_effect

"Write list of typed formulae in PDDL syntax."
function write_typed_list(formulae::Vector{<:Term}, types::Vector{Symbol})
    if all(x -> x == :object, types)
        return join(write_subformula.(formulae), " ")
    end
    str = ""
    for (i, (f, t)) in enumerate(zip(formulae, types))
        str *= "$(write_subformula(f)) "
        if i == length(formulae) || t != types[i+1]
            str *= "- $t "
        end
    end
    return str[1:end-1]
end

function write_typed_list(formulae::Vector{<:Term})
    formulae, types = zip([(t.args[1], t.name) for t in formulae]...)
    return write_typed_list(collect(formulae), collect(types))
end

function write_typed_list(formulae::Vector{<:Term}, types::Dict{<:Term,Symbol})
    types = Symbol[types[f] for f in formulae]
    return write_typed_list(formulae, types)
end

"Indent list of PDDL constants."
function indent_const_list(str::String, indent::Int, maxchars::Int=80)
    if length(str) < maxchars - indent return str end
    lines = String[]
    while length(str) > 0
        l_end = min(maxchars - indent, length(str))
        if l_end != length(str) l_end = findlast(' ', str[1:l_end]) end
        if isnothing(l_end) l_end = findfirst(' ', str) end
        if isnothing(l_end) l_end = length(str) end
        push!(lines, str[1:l_end])
        str = str[l_end+1:end]
    end
    return join(lines, "\n" * ' '^indent)
end

"Indent list of typed PDDL constants."
function indent_typed_list(str::String, indent::Int, maxchars::Int=80)
    if length(str) < maxchars - indent return str end
    if !occursin(" - ", str) return indent_const_list(strs[1]) end
    substrs = String[]
    while length(str) > 0
        idxs = findnext(r" - (\S+)", str, 1)
        if isnothing(idxs) idxs = [length(str)] end
        push!(substrs, str[1:idxs[end]])
        str = strip(str[idxs[end]+1:end])
    end
    substrs = indent_const_list.(substrs, indent, maxchars)
    return join(substrs, "\n" * ' '^indent)
end

"Write formula in PDDL syntax."
function write_formula(f::Compound)
    if f.name in [:exists, :forall]
        typecond, body = f.args
        var_str = write_typed_list(flatten_conjs(typecond))
        body_str = write_formula(body)
        return "($(f.name) ($var_str) $body_str)"
    elseif f.name in [:and, :or, :not, :when, :imply,
                      :(==), :>, :<, :!=, :>=, :<=, :+, :-, :*, :/,
                      :assign, :increase, :decrease,
                      Symbol("scale-up"), Symbol("scale-down")]
        name = f.name == :(==) ? "=" : string(f.name)
        args = join([write_formula(a) for a in f.args], " ")
        return "($name $args)"
    else
        args = join([write_subformula(a) for a in f.args], " ")
        return "($(f.name) $args)"
    end
end
write_formula(f::Var) = "?" * lowercasefirst(repr(f))
write_formula(f::Const) = "(" * repr(f) * ")"
write_formula(::Nothing) = ""

write_subformula(f::Compound) = write_formula(f)
write_subformula(f::Var) = "?" * lowercasefirst(repr(f))
write_subformula(f::Const) = repr(f)

"Write to string in PDDL syntax."
write_pddl(f::Term) = write_formula(f)

"Write domain in PDDL syntax."
function write_domain(domain::Domain, indent::Int=2)
    strs = Dict{Symbol,String}()
    fields = [:requirements, :types, :constants, :predicates, :functions]
    strs[:requirements] = write_requirements(get_requirements(domain))
    strs[:types] = write_typetree(get_typetree(domain))
    strs[:constants] = indent_typed_list(
        write_typed_list(get_constants(domain), get_constypes(domain)), 14)
    strs[:predicates] = write_signatures(get_predicates(domain))
    strs[:functions] = write_signatures(get_functions(domain))
    strs = ["($(repr(k)) $(strs[k]))" for k in fields
            if haskey(strs, k) && length(strs[k]) > 0]
    append!(strs, write_axiom.(values(get_axioms(domain))))
    append!(strs, write_action.(values(get_actions(domain)), indent=3))
    pushfirst!(strs, "(define (domain $(domain.name))")
    return join(strs, "\n" * ' '^indent) * "\n)"
end
write_pddl(domain::Domain) = write_domain(domain)

"Write domain requirements."
function write_requirements(requirements::Dict{Symbol,Bool})
    reqs = Set(k for (k, v) in requirements if v)
    for (k, implied) in IMPLIED_REQUIREMENTS
        if haskey(requirements, k)
            setdiff!(reqs, implied)
        end
    end
    return join([":$r" for r in reqs], " ")
end

"Write domain typetree."
function write_typetree(typetree)
    strs = Dict{Symbol,String}()
    for (type, subtypes) in pairs(typetree)
        if isempty(subtypes) || type == :object continue end
        subtype_str = join(subtypes, " ")
        strs[type] = "$subtype_str - $type"
    end
    maxtypes = collect(setdiff(get(typetree, :object, Symbol[]), keys(strs)))
    strs[:object] = join(maxtypes, " ")
    return strip(join(values(strs), " "))
end

"Write domain predicates or functions."
function write_signatures(signatures)
    strs = String[]
    for (name, sig) in pairs(signatures)
        args = collect(Term, sig.args)
        types = collect(Symbol, sig.argtypes)
        args_str = write_typed_list(args, types)
        sig_str = isempty(args_str) ? "($name)" : "($name $args_str)"
        if sig.type != :boolean && sig.type != :numeric
            sig_str = "$sig_str - $(sig.type)"
        end
        push!(strs, sig_str)
    end
    return join(strs, " ")
end

"Write PDDL axiom / derived predicate."
function write_axiom(c::Clause; key=":derived")
    head_str = write_formula(c.head)
    body_str = length(c.body) == 1 ? write_formula(c.body[1]) :
        "(and " * join(write_formula.(c.body), " ") * ")"
    return "($key $head_str $body_str)"
end

"Write action in PDDL syntax."
function write_action(action::Action; indent::Int=1)
    strs = Dict{Symbol,String}()
    fields = [:action, :parameters, :precondition, :effect]
    strs[:action] = string(get_name(action))
    strs[:parameters] =
        "(" * write_typed_list(get_argvars(action), get_argtypes(action)) * ")"
    strs[:precondition] = write_formula(get_precond(action))
    strs[:effect] = write_formula(get_effect(action))
    strs = ["$(repr(k)) $(strs[k])" for k in fields]
    return "(" * join(strs, "\n" * ' '^indent) * ")"
end
write_pddl(action::Action) = write_action(action)

"Write problem in PDDL syntax."
function write_problem(problem::Problem, indent::Int=2)
    strs = Dict{Symbol,String}()
    fields = [:domain, :objects, :init, :goal, :metric, :constraints]
    strs[:domain] = string(problem.domain)
    strs[:objects] = indent_typed_list(
        write_typed_list(get_objects(problem), get_objtypes(problem)), 12)
    strs[:init] = write_init(get_init_terms(problem), 9)
    strs[:goal] = write_formula(get_goal(problem))
    strs[:metric] = write_formula(get_metric(problem))[2:end-1]
    strs[:constraints] = write_formula(get_constraints(problem))
    strs = ["($(repr(k)) $(strs[k]))" for k in fields
            if haskey(strs, k) && length(strs[k]) > 0]
    pushfirst!(strs, "(define (problem $(get_name(problem)))")
    return join(strs, "\n" * ' '^indent) * "\n)"
end
write_pddl(problem::Problem) = write_problem(problem)

"Write initial problem formulae in PDDL syntax."
function write_init(init::Vector{<:Term}, indent::Int=2, maxchars::Int=80)
    strs = write_formula.(init)
    if sum(length.(strs)) + length("(:init )") < maxchars
        return join(strs, " ") end
    return join(strs, "\n" * ' '^indent)
end

"Save PDDL domain to specified path."
function save_domain(path::String, domain::Domain)
    open(f->write(f, write_domain(domain)), path, "w")
    return path
end

"Save PDDL problem to specified path."
function save_problem(path::String, problem::Problem)
    open(f->write(f, write_problem(problem)), path, "w")
    return path
end

Base.show(io::IO, ::MIME"text/pddl", domain::Domain) =
    print(io, write_pddl(domain))
Base.show(io::IO, ::MIME"text/pddl", problem::Problem) =
    print(io, write_pddl(problem))
Base.show(io::IO, ::MIME"text/pddl", action::Action) =
    print(io, write_pddl(action))
Base.show(io::IO, ::MIME"text/pddl", term::Term) =
    print(io, write_pddl(term))

end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type RandomVariable end
abstract type RandomVector end

struct RawStoch{N}
    name::String
    random_variables::NTuple{N, Vector{RandomVariable}}
    random_vectors::NTuple{N, Vector{RandomVector}}

    function RawStoch(name::String, ran_vars::NTuple{N, Vector{RandomVariable}}, ran_vectors::NTuple{N, Vector{RandomVector}}) where N
        return new{N}(name, ran_vars, ran_vectors)
    end
end

struct IndepDiscrete{T <: AbstractFloat} <: RandomVariable
    rowcol::RowCol
    inclusion::InclusionType
    support::Vector{T}
    probabilities::Vector{T}

    function IndepDiscrete(::Type{T}, rowcol::RowCol, inclusion::InclusionType) where T <: AbstractFloat
        inclusion in INCLUSIONS || error("Unknown inclusion $inclusion")
        return new{T}(rowcol, inclusion, Vector{T}(), Vector{T}())
    end
end

struct IndepDistribution{T <: AbstractFloat} <: RandomVariable
    rowcol::RowCol
    inclusion::InclusionType
    distribution::DistributionType
    parameters::Pair{T}

    function IndepDistribution(rowcol::RowCol, inclusion::InclusionType, distribution::DistributionType, parameters::Pair{T}) where T <: AbstractFloat
        inclusion in INCLUSIONS || error("Unknown inclusion $inclusion")
        distribution in DISTRIBUTIONS || error("Unknown distribution $distribution")
        return new{T}(rowcol, inclusion, distribution, parameters)
    end
end

struct BlockDiscrete{T <: AbstractFloat} <: RandomVector
    inclusion::InclusionType
    support::Vector{Dict{RowCol,T}}
    probabilities::Vector{T}

    function BlockDiscrete(::Type{T}, inclusion::InclusionType) where T
        inclusion in INCLUSIONS || error("Unknown inclusion $inclusion")
        return new{T}(inclusion, Vector{Dict{RowCol,T}}(), Vector{T}())
    end
end

function parse_sto(::Type{T}, tim::RawTime, cor::RawCor, filename::AbstractString; N::Integer = 2) where T <: AbstractFloat
    # Initialize auxiliary variables
    name      = "SLP"
    mode      = INDEP
    dist      = DISCRETE
    inclusion = REPLACE
    N         = num_stages(tim)
    stage     = 2
    v₁        = zero(T)
    v₂        = zero(T)
    block     = 1
    # Prepare sto data
    random_vars = ntuple(Val(N-1)) do i
        Vector{RandomVariable}()
    end
    random_vecs = ntuple(Val(N-1)) do i
        Vector{RandomVector}()
    end
    # Bookkeeping
    ranvar_indices = ntuple(Val(N-1)) do i
        Dict{RowCol, Int}()
    end
    ranvec_indices = ntuple(Val(N-1)) do i
        Dict{Int, Int}()
    end
    # Parse the file
    data = open(filename) do io
        firstline = split(readline(io))
        if Symbol(firstline[1]) == :STOCH
            name = join(firstline[2:end], " ")
        else
            throw(ArgumentError("`STOCH` field is expected on the first line."))
        end
        for line in eachline(io)
            if mode == END
                # Parse finished
                break
            end
            words = split(line)
            first_word = Symbol(words[1])
            if first_word in STO_MODES
                mode = first_word
                if length(words) == 2
                    dist = DistributionType(words[2])
                    dist in DISTRIBUTIONS || error("Unknown distribution $dist.")
                elseif length(words) == 3
                    dist = DistributionType(words[2])
                    dist in DISTRIBUTIONS || error("Unknown distribution $dist.")
                    inclusion = InclusionType(words[3])
                    inclusion in INCLUSIONS || error("Unknown inclusion $inclusion.")
                end
                 continue
            end
            if mode == INDEP
                length(words) == 4 || length(words) == 5 || error("Malformed sto file at line $words.")
                col = Col(words[1])
                col = col == cor.rhsname ? RHS : col
                row = Row(words[2])
                if length(words) == 4
                    v₁ = parse(Float64, words[3])
                    v₂ = parse(Float64, words[4])
                elseif length(words) == 5
                    v₁ = parse(Float64, words[3])
                    v₂ = parse(Float64, words[5])
                    stage_name = Symbol(words[4])
                    haskey(tim.stages, stage_name) || error("Stage name $stage_name not specified in .tim file.")
                    stage = tim.stages[stage_name]
                else
                    error("Malformed sto file at line $words")
                end
                if dist == DISCRETE
                    # Check if random variable exists already
                    idx = get(ranvar_indices[stage - 1], (row, col), 0)
                    ran_var = if idx > 0
                        ran_var = random_vars[stage - 1][idx]
                        # Sanity check
                        ran_var isa IndepDiscrete || error("Random variable at $((row, col)) has more than one specified distribution.")
                        ran_var
                    else
                        # Create new INDEP random variable
                        ran_var = IndepDiscrete(T, (row,col), inclusion)
                        push!(random_vars[stage - 1], ran_var)
                        # Bookkeep
                        ranvar_indices[stage - 1][(row, col)] = length(random_vars[stage - 1])
                        ran_var
                    end
                    # Update values
                    push!(ran_var.support, v₁)
                    push!(ran_var.probabilities, v₂)
                else
                    haskey(ranvar_indices[stage - 1], (row, col)) && error("Random variable at $((row, col)) has more than one specified distribution.")
                    push!(random_vars[stage - 1], IndepDistribution((row, col), inclusion, dist, Pair(v₁, v₂)))
                    # Bookkeep
                    ranvar_indices[(row, col)] = length(random_vars[stage - 1])
                end
            elseif mode == BLOCKS
                if dist == DISCRETE
                    if words[1] == "BL"
                        length(words) == 4 || error("Malformed sto file at line $words.")
                        stage_name = Symbol(words[3])
                        haskey(tim.stages, stage_name) || error("Stage name $stage_name not specified in .tim file.")
                        stage = tim.stages[stage_name]
                        prob = parse(Float64, words[4])
                        # Check if random vector exists already
                        idx = get(ranvec_indices[stage - 1], block, 0)
                        ran_vec = if idx > 0
                            ran_vec = random_vecs[stage - 1][idx]
                            # Sanity check
                            ran_vec isa BlockDiscrete || error("Random vector at $block has more than one specified distribution.")
                            ran_vec
                        else
                            # Create new DISCRETE random vector
                            ran_vec = BlockDiscrete(T, inclusion)
                            push!(random_vecs[stage - 1], ran_vec)
                            # Bookkeep
                            block += 1
                            ranvec_indices[stage - 1][block] = length(random_vecs[stage - 1])
                            ran_vec
                        end
                        # Update probability
                        push!(ran_vec.support, Dict{RowCol,T}())
                        push!(ran_vec.probabilities, prob)
                        continue
                    else
                        length(words) == 3 || error("Malformed sto file at line $words.")
                        col = Col(words[1])
                        col = col == cor.rhsname ? RHS : col
                        row = Row(words[2])
                        val = parse(Float64, words[3])
                        idx = ranvec_indices[stage - 1][block]
                        ran_vec = random_vecs[stage - 1][idx]
                        ran_vec.support[end][(row, col)] = val
                    end
                else
                    error("Distribution $dist is not supported.")
                end
            elseif mode == SCENARIOS
                if words[1] == "SC"
                    length(words) == 5 || error("Malformed+ sto file at line $words.")
                    words[3] == "ROOT" || error("SCEN format currently only supported for two-stage models.")
                    stage_name = Symbol(words[5])
                    haskey(tim.stages, stage_name) || error("Stage name $stage_name not specified in .tim file.")
                    stage = tim.stages[stage_name]
                    prob = parse(Float64, words[4])
                    # Check if random vector exists already
                    idx = get(ranvec_indices[stage - 1], block, 0)
                    ran_vec = if idx > 0
                        ran_vec = random_vecs[stage - 1][idx]
                        # Sanity check
                        ran_vec isa BlockDiscrete || error("Random vector at $block has more than one specified distribution.")
                        ran_vec
                    else
                        # Create new DISCRETE random vector
                        ran_vec = BlockDiscrete(T, inclusion)
                        push!(random_vecs[stage - 1], ran_vec)
                        # Bookkeep
                        block += 1
                        ranvec_indices[stage - 1][block] = length(random_vecs[stage - 1])
                        ran_vec
                    end
                    # Update probability
                    push!(ran_vec.support, Dict{RowCol,T}())
                    push!(ran_vec.probabilities, prob)
                    continue
                else
                    length(words) == 3 || error("Malformed sto file at line $words.")
                    col = Col(words[1])
                    col = col == cor.rhsname ? RHS : col
                    row = Row(words[2])
                    val = parse(Float64, words[3])
                    idx = ranvec_indices[stage - 1][block]
                    ran_vec = random_vecs[stage - 1][idx]
                    ran_vec.support[end][(row, col)] = val
                end
            else
                throw(ArgumentError("$(mode) is not a valid sto file mode."))
            end
        end
    end
    # Return raw data
    return RawStoch(name, random_vars, random_vecs)
end
parse_sto(tim::RawTime, cor::RawCor, filename::AbstractString) = parse_sto(Float64, tim, cor, filename)

function uncertainty_template(sto::RawStoch{N}, cor::RawCor, model::LPData{T, Matrix{T}}, stage::Integer) where {N, T}
    1 <= stage <= N || error("$(stage.id + 1) not in range 2 to $(N + 1).")
    map = model.indexmap
    # Prepare data
    Δc  = zero(model.c₁)
    ΔA  = zero(model.A)
    Δb  = zero(model.b)
    Δd₁ = zero(model.d₁)
    ΔC  = zero(model.C)
    Δd₂ = zero(model.d₂)
    # Collect rows, columns, and inclusions
    uncertainty_structure = Dict{RowCol, InclusionType}()
    for ran_var in sto.random_variables[stage]
        uncertainty_structure[ran_var.rowcol] = ran_var.inclusion
    end
    for ran_vec in sto.random_vectors[stage]
        isempty(ran_vec.support[1]) && error("Block $block has empty support.")
        for (rowcol, val) in ran_vec.support[1]
            rowcol = (row, col)
            uncertainty_structure[rowcol] = ran_vec.inclusion
        end
    end
    # Loop over all scenarios
    for (rowcol, inclusion) in uncertainty_structure
        (row, col) = rowcol
        rowsense = cor.rows[row][3]
        if row == cor.objname
            # Objective
            j = map[rowcol][2]
            if inclusion == REPLACE
                Δc[j] = -model.c[j]
            else
                Δc[j] = model.c[j]
            end
        elseif col == RHS
            # RHS
            i = map[rowcol][1]
            if rowsense == EQ
                if inclusion == REPLACE
                    Δb[i] = -model.b[i]
                else
                    Δb[i] = model.b[i]
                end
            elseif rowsense == LEQ
                if inclusion == REPLACE
                    Δd₂[i] = -model.d₂[i]
                    haskey(cor.ranges, row) && (Δd₁[i] = -model.d₁[i])
                else
                    Δd₂[i] = model.d₂[i]
                    haskey(cor.ranges, row) && (Δd₁[i] = model.d₁[i])
                end
            elseif rowsense == GEQ
                if inclusion == REPLACE
                    Δd₁[i] = -model.d₁[i]
                    haskey(cor.ranges, row) && (Δd₂[i] = -model.d₂[i])
                else
                    Δd₁[i] = model.d₁[i]
                    haskey(cor.ranges, row) && (Δd₂[i] = model.d₂[i])
                end
            end
        else
            # Matrix coeffs
            (i, j, _) = map[rowcol]
            if rowsense == EQ
                if inclusion == REPLACE
                    ΔA[i,j] = -model.A[i,j]
                else
                    ΔA[i] = model.A[i,j]
                end
            else
                if inclusion == REPLACE
                    ΔC[i,j] = -model.C[i,j]
                else
                    ΔC[i,j] = model.C[i,j]
                end
            end
        end
    end
    # Return template
    return LPData(T, Δc, zero(T), ΔA, Δb, Δd₁, ΔC, Δd₂, model.lb, model.ub, model.is_binary, model.is_integer, model.indexmap)
end

function uncertainty_template(sto::RawStoch{N}, cor::RawCor, model::LPData{T, SparseMatrixCSC{T,Int}}, stage::Integer) where {N, T}
    1 <= stage <= N || error("$(stage.id + 1) not in range 2 to $(N + 1).")
    map = model.indexmap
    # Get sizes
    n   = size(model.A, 2)
    m₁  = size(model.A, 1)
    m₂  = size(model.C, 1)
    # Prepare vectors
    Δc  = zero(model.c₁)
    Δb  = zero(model.b)
    Δd₁ = zero(model.d₁)
    Δd₂ = zero(model.d₂)
    # Prepare matrices
    ΔAᵢ = Vector{Int}()
    ΔAⱼ = Vector{Int}()
    ΔAᵥ = Vector{T}()
    ΔCᵢ = Vector{Int}()
    ΔCⱼ = Vector{Int}()
    ΔCᵥ = Vector{T}()
    # Collect rows, columns, and inclusions
    uncertainty_structure = Dict{RowCol, InclusionType}()
    for ran_var in sto.random_variables[stage]
        uncertainty_structure[ran_var.rowcol] = ran_var.inclusion
    end
    for ran_vec in sto.random_vectors[stage]
        isempty(ran_vec.support[1]) && error("Block $block has empty support.")
        for (rowcol, val) in ran_vec.support[1]
            uncertainty_structure[rowcol] = ran_vec.inclusion
        end
    end
    # Loop over all scenarios
    for (rowcol, inclusion) in uncertainty_structure
        (row, col) = rowcol
        rowsense = cor.rows[row][3]
        if row == cor.objname
            # Objective
            j = map[rowcol][2]
            if inclusion == REPLACE
                Δc[j] = -model.c₁[j]
            else
                Δc[j] = model.c₁[j]
            end
        elseif col == RHS
            # RHS
            i = map[rowcol][1]
            if rowsense == EQ
                if inclusion == REPLACE
                    Δb[i] = -model.b[i]
                else
                    Δb[i] = model.b[i]
                end
            elseif rowsense == LEQ
                if inclusion == REPLACE
                    Δd₂[i] = -model.d₂[i]
                    haskey(cor.ranges, row) && (Δd₁[i] = -model.d₂[i])
                else
                    Δd₂[i] = model.d₂[i]
                    haskey(cor.ranges, row) && (Δd₁[i] = model.d₂[i])
                end
            elseif rowsense == GEQ
                if inclusion == REPLACE
                    Δd₁[i] = -model.d₁[i]
                    haskey(cor.ranges, row) && (Δd₂[i] = -model.d₁[i])
                else
                    Δd₁[i] = model.d₁[i]
                    haskey(cor.ranges, row) && (Δd₂[i] = model.d₁[i])
                end
            end
        else
            # Matrix coeffs
            (i, j, _) = map[rowcol]
            if rowsense == EQ
                push!(ΔAᵢ, i)
                push!(ΔAⱼ, j)
                if inclusion == REPLACE
                    push!(ΔAᵥ, -model.A[i,j])
                else
                    push!(ΔAᵥ, model.A[i,j])
                end
            else
                push!(ΔCᵢ, i)
                push!(ΔCⱼ, j)
                if inclusion == REPLACE
                    push!(ΔCᵥ, -model.C[i,j])
                else
                    push!(ΔCᵥ, model.C[i,j])
                end
            end
        end
    end
    # Create sparse structures
    ΔA  = sparse(ΔAᵢ, ΔAⱼ, ΔAᵥ, m₁, n)
    ΔC  = sparse(ΔCᵢ, ΔCⱼ, ΔCᵥ, m₂, n)
    # Return template
    return LPData(T, Δc, zero(T), ΔA, Δb, Δd₁, ΔC, Δd₂, model.lb, model.ub, model.is_binary, model.is_integer, model.indexmap)
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum__DataType <: ProtoEnum
    DT_INVALID::Int32
    DT_FLOAT::Int32
    DT_DOUBLE::Int32
    DT_INT32::Int32
    DT_UINT8::Int32
    DT_INT16::Int32
    DT_INT8::Int32
    DT_STRING::Int32
    DT_COMPLEX64::Int32
    DT_INT64::Int32
    DT_BOOL::Int32
    DT_QINT8::Int32
    DT_QUINT8::Int32
    DT_QINT32::Int32
    DT_BFLOAT16::Int32
    DT_QINT16::Int32
    DT_QUINT16::Int32
    DT_UINT16::Int32
    DT_COMPLEX128::Int32
    DT_HALF::Int32
    DT_RESOURCE::Int32
    DT_VARIANT::Int32
    DT_UINT32::Int32
    DT_UINT64::Int32
    DT_FLOAT_REF::Int32
    DT_DOUBLE_REF::Int32
    DT_INT32_REF::Int32
    DT_UINT8_REF::Int32
    DT_INT16_REF::Int32
    DT_INT8_REF::Int32
    DT_STRING_REF::Int32
    DT_COMPLEX64_REF::Int32
    DT_INT64_REF::Int32
    DT_BOOL_REF::Int32
    DT_QINT8_REF::Int32
    DT_QUINT8_REF::Int32
    DT_QINT32_REF::Int32
    DT_BFLOAT16_REF::Int32
    DT_QINT16_REF::Int32
    DT_QUINT16_REF::Int32
    DT_UINT16_REF::Int32
    DT_COMPLEX128_REF::Int32
    DT_HALF_REF::Int32
    DT_RESOURCE_REF::Int32
    DT_VARIANT_REF::Int32
    DT_UINT32_REF::Int32
    DT_UINT64_REF::Int32
    __enum__DataType() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123)
end #struct __enum__DataType
const _DataType = __enum__DataType()

export _DataType


text:
using Primes

export SlotEncoding, ScalarEncoding

struct ScalarEncoding{T, R} <: AbstractArray{T, 0}
    plain::R
end

ScalarEncoding{T}(plain::R) where {T,R} =
    ScalarEncoding{T,R}(plain)

Base.axes(a::ScalarEncoding) = ()
Base.size(a::ScalarEncoding) = 1
Base.getindex(a::ScalarEncoding) = NTT.coeffs_primal(a.plain)[0]
Base.setindex!(a::ScalarEncoding, v) =  NTT.coeffs_primal(a.plain)[0] = v

function Base.convert(T::Type{<:NTT.RingElement}, s::ScalarEncoding)
    convert(T, s.plain)
end

function Base.showarg(io::IO, v::ScalarEncoding{T, R}, toplevel) where {T,R}
    # In the default case where we're just encoding into a ring
    # element, no need to make special note of it.
    if R <: NTT.RingElement
        print(io, "ScalarEncoding{", T, "}")
    else
        print(io, typeof(v))
    end
end

struct SlotEncoding{T, S, R} <: AbstractVector{T}
    plain::R
end

function SlotEncoding(r::R) where R<:NTT.RingElement{ℛ} where ℛ
    N = degree(ℛ)
    p = modulus(eltype(R))
    if Primes.isprime(p) && gcd(p-1, 2N) == 2N
        # TODO: Should allow CRTEncoded here also
        return SlotEncoding{eltype(R), eltype(R), R}(R(NTT.coeffs_dual(r), nothing))
    end
    error("Not supported yet")
end

Base.axes(s::SlotEncoding) = Base.axes(s.plain)
Base.size(s::SlotEncoding) = Base.size(s.plain)
Base.getindex(s::SlotEncoding, idxs...) = getindex(s.plain, idxs...)
Base.setindex!(s::SlotEncoding, args...) = setindex!(s.plain, args...)

function Base.convert(::Type{NTT.RingElement}, s::SlotEncoding)
    typeof(s.plain)(nothing, NTT.coeffs_primal(s.plain))
end

function Base.convert(T::Type{<:NTT.RingElement}, s::SlotEncoding)
    T(nothing, NTT.coeffs_primal(s.plain))
end

function SlotEncoding(r::R) where R<:NmodPolyRing
    f = modulus(r)
    p = modulus(base_ring(r))
end


text:
push!(LOAD_PATH,"../src/") # load module from local directory

using TheDancer, Documenter

makedocs()


text:
"""
    signals = init_zoom_rubberband(canvas::GtkReactive.Canvas,
                                   zr::Signal{ZoomRegion},
                                   initiate = btn->(btn.button == 1 && btn.clicktype == BUTTON_PRESS && btn.modifiers == CONTROL),
                                   reset = btn->(btn.button == 1 && btn.clicktype == DOUBLE_BUTTON_PRESS && btn.modifiers == CONTROL),
                                   minpixels = 2)

Initialize rubber-band selection that updates `zr`. `signals` is a
dictionary holding the Reactive.jl signals needed for rubber-banding;
you can push `true/false` to `signals["enabled"]` to turn rubber
banding on and off, respectively. Your application is responsible for
making sure that `signals` does not get garbage-collected (which would
turn off rubberbanding).

`initiate(btn)` returns `true` when the condition for starting a
rubber-band selection has been met (by default, clicking mouse button
1). The argument `btn` is a [`MouseButton`](@ref) event. `reset(btn)`
returns true when restoring the full view (by default, double-clicking
mouse button 1). `minpixels` can be used for aborting rubber-band
selections smaller than some threshold.
"""
function init_zoom_rubberband(canvas::Canvas{U},
                              zr::Signal{ZoomRegion{T}},
                              initiate::Function = zrb_init_default,
                              reset::Function = zrb_reset_default,
                              minpixels::Integer = 2) where {U,T}
    enabled = Signal(true)
    active = Signal(false)
    function update_zr(widget, bb)
        push!(active, false)
        fv = value(zr).fullview
        push!(zr, ZoomRegion(fv, XY(interior(bb.xmin..bb.xmax, fv.x),
                                    interior(bb.ymin..bb.ymax, fv.y))
                             ))
        nothing
    end
    rb = RubberBand(XY{U}(-1,-1), XY{U}(-1,-1), false, minpixels)
    dummybtn = MouseButton{U}()
    local ctxcopy
    init = map(filterwhen(enabled, dummybtn, canvas.mouse.buttonpress)) do btn
        if initiate(btn)
            push!(active, true)
            ctxcopy = copy(getgc(canvas))
            rb.pos1 = rb.pos2 = btn.position
        elseif reset(btn)
            push!(active, false)  # double-clicks need to cancel the previous single-click
            push!(zr, GtkReactive.reset(value(zr)))
        end
        nothing
    end
    drag = map(filterwhen(active, dummybtn, canvas.mouse.motion)) do btn
        btn.button == 0 && return nothing
        rubberband_move(canvas, rb, btn, ctxcopy)
    end
    finish = map(filterwhen(active, dummybtn, canvas.mouse.buttonrelease)) do btn
        btn.button == 0 && return nothing
        push!(active, false)
        rubberband_stop(canvas, rb, btn, ctxcopy, update_zr)
    end
    Dict("enabled"=>enabled, "active"=>active, "init"=>init, "drag"=>drag, "finish"=>finish)
end

zrb_init_default(btn) = btn.button == 1 && btn.clicktype == BUTTON_PRESS && (btn.modifiers & 0x0f) == CONTROL
zrb_reset_default(btn) = btn.button == 1 && btn.clicktype == DOUBLE_BUTTON_PRESS && (btn.modifiers & 0x0f) == CONTROL

# For rubberband, we draw the selection region on the front canvas, and repair
# by copying from the back.
mutable struct RubberBand{U}
    pos1::XY{U}
    pos2::XY{U}
    moved::Bool
    minpixels::Int
end

const dash   = Float64[3.0,3.0]
const nodash = Float64[]

function rb_erase(r::GraphicsContext, rb::RubberBand, ctxcopy)
    # Erase the previous rubberband by copying from back surface to front
    rb_set(r, rb)
    # Because line widths are expressed in pixels, let's go to device units
    save(r)
    reset_transform(r)
    save(ctxcopy)
    reset_transform(ctxcopy)
    set_source(r, ctxcopy)
    set_line_width(r, 3)
    set_dash(r, nodash)
    stroke(r)
    restore(r)
    restore(ctxcopy)
end

function rb_draw(r::GraphicsContext, rb::RubberBand)
    rb_set(r, rb)
    save(r)
    reset_transform(r)
    set_line_width(r, 1)
    set_dash(r, dash, 3.0)
    set_source_rgb(r, 1, 1, 1)
    stroke_preserve(r)
    set_dash(r, dash, 0.0)
    set_source_rgb(r, 0, 0, 0)
    stroke(r)
    restore(r)
end

function rb_set(r::GraphicsContext, rb::RubberBand)
    x1, y1 = rb.pos1.x, rb.pos1.y
    x2, y2 = rb.pos2.x, rb.pos2.y
    rectangle(r, x1, y1, x2 - x1, y2 - y1)
end

function rubberband_move(c::Canvas, rb::RubberBand, btn, ctxcopy)
    if btn.button == 0
        return nothing
    end
    r = getgc(c)
    if rb.moved
        rb_erase(r, rb, ctxcopy)
    end
    rb.moved = true
    # Draw the new rubberband
    rb.pos2 = btn.position
    rb_draw(r, rb)
    reveal(c, false)
    nothing
end

function rubberband_stop(c::GtkReactive.Canvas, rb::RubberBand, btn, ctxcopy, callback_done)
    if !rb.moved
        return nothing
    end
    r = getgc(c)
    rb_set(r, rb)
    rb_erase(r, rb, ctxcopy)
    reveal(c, false)
    pos = btn.position
    x, y = pos.x, pos.y
    x1, y1 = rb.pos1.x, rb.pos1.y
    xd, yd = convertunits(DeviceUnit, r, x, y)
    x1d, y1d = convertunits(DeviceUnit, r, x1, y1)
    if abs(x1d-xd) > rb.minpixels || abs(y1d-yd) > rb.minpixels
        # It moved sufficiently, let's execute the callback
        bb = BoundingBox(min(x1,x), max(x1,x), min(y1,y), max(y1,y))
        callback_done(c, bb)
    end
    nothing
end


text:
function create(config::Config; dir::AbstractString="", name::AbstractString="",
    julia::AbstractString="", reqs_file::AbstractString="", reqs_type::Symbol=:REQUIRE)

    init(config)

    pg = PlaygroundConfig(config, dir, name)
    create_paths(pg)

    info("Playground folders created")

    if julia != ""
        mklink(joinpath(config.dir.bin, julia), pg.julia_path)
    else
        sys_julia_path = abspath(readchomp(`which julia`))
        mklink(sys_julia_path, pg.julia_path)
    end

    if dir != "" && name != ""
        mklink(pg.root_path, abspath(joinpath(config.dir.store, name)))
    end

    ENV["JULIA_PKGDIR"] = pg.pkg_path

    if reqs_file != "" && ispath(reqs_file)
        if basename(reqs_file) == "DECLARE" || reqs_type == :DECLARE
            if ispath(DECLARATIVE_PACKAGES_DIR)
                ENV["DECLARE"] = reqs_file
                ENV["JULIA_PKGDIR"] = joinpath(pg.pkg_path, ".declare_packages")
                mkpath(ENV["JULIA_PKGDIR"])
                copy(reqs_file, joinpath(ENV["JULIA_PKGDIR"], "DECLARE"))

                old_folders = readdir(pg.pkg_path)
                info("Installing packages from DECLARE file $reqs_file...")
                passed = false
                try
                    run(`$(pg.julia_path) $(DECLARATIVE_PACKAGES_DIR)/src/installpackages.jl`)
                    passed = true
                catch
                    passed = false
                    warn("Failed to install from DECLARE file. Perhaps there is something wrong with your DECLARE file or you need to update DeclarativePackages.jl")
                end

                if passed
                    # Despite a declarative packages JULIA_PKGDIR being provided
                    # DeclarativePackages will just create a temporary folder
                    # in the parent directory. So we compare the pg.pkg_path before
                    # and after running installpackages.jl. Once, we find the generated folder
                    # that isn't just symlinks we copy the version folder over.
                    new_folders = readdir(pg.pkg_path)
                    to_delete = []
                    for f in new_folders
                        if !(f in old_folders)
                            declared_folder = joinpath(pg.pkg_path, f)
                            push!(to_delete, declared_folder)
                            if !islink(declared_folder) && isdir(declared_folder)
                                for v in readdir(declared_folder)
                                    vfolder = joinpath(declared_folder, v)
                                    copy(vfolder, joinpath(pg.pkg_path, v))
                                end
                            end
                        end
                    end
                    for d in to_delete
                        try
                            rm(d, recursive=true)
                        catch
                            @unix_only begin
                                run(`chmod -R 755 $d`)
                                run(`rm -rf $d`)
                            end
                        end
                    end
                end
            else
                warn("DeclarativePackages isn't installed")
            end
        elseif basename(reqs_file) == "REQUIRE" || reqs_type == :REQUIRE
            info("Installing packages from REQUIRE file $reqs_file...")
            run(`$(pg.julia_path) -e Pkg.init()`)
            for v in readdir(pg.pkg_path)
                copy(reqs_file, joinpath(pg.pkg_path, v, "REQUIRE"))
                try
                    run(`$(pg.julia_path) -e Pkg.resolve()`)
                catch
                    warn("Failed to resolve requirements. Perhaps there is something wrong with your REQUIRE file.")
                end
            end
        else
            error("Unknown requirements file type $reqs_type for file $reqs_file")
        end
    else
        run(`$(pg.julia_path) -e Pkg.init()`)
    end
end


text:
#demo_cktanalysis.jl: Run sample code for CircuitAnalysis tools
#-------------------------------------------------------------------------------

using CMDimCircuits
using CMDimCircuits.CircuitAnalysis

#==Input data
===============================================================================#
f = collect(1:5)*1e9

#==Sample code
===============================================================================#
@info("Define circuit elements:")
#@show r = resistance(20.5)
#@show g = conductance(1/50)
@show l = inductance(20e-9)
@show c = capacitance(5e-15)

@info("Compute single-frequency values:")
@show ycap = admittance(c, f=1e9)
@show zcap = impedance(c, f=1e9)

@info("Compute frequency-dependent values:")
@show f
@show ycap = admittance(c, f=f)
@show zcap = impedance(c, f=f)

#println("\nOperations:")
#TODO: fix operations

@info("Perform unit conversions:")
@show dB10(2), dB20(2)
@show dB(2,:Wratio), dB(2,:Vratio), dB(2,:Iratio)
@show dBm(2,:W), dBm(2,:VRMS), dBm(2,:Vpk)
@show dBW(2,:W), dBW(2,:VRMS), dBW(2,:Vpk)
@show Vpk(2,:W; R=50), Vpk(2,:VRMS)
@show Ipk(2,:W; R=50), Ipk(2,:IRMS)
@show VRMS(2,:W; R=50), VRMS(2,:Vpk)
@show IRMS(2,:W; R=50), IRMS(2,:Ipk)

:SampleCode_Executed


text:
function central_diff(Qbase,cellxmax,cellymax,D,vecAx,vecAy,volume)
    E_hat = zeros(cellxmax+1,cellymax)

    Threads.@threads for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            T_av = 0.5*(Qbase[i,j]+Qbase[i-1,j])
            dTdxi = (Qbase[i,j]-Qbase[i-1,j])
            dTdeta = 0.25*((Qbase[i,j+1]-Qbase[i,j-1]) + (Qbase[i-1,j+1]-Qbase[i-1,j-1]))

            vecAy_xav    = 0.25*( vecAy[i,j,1] + vecAy[i-1,j,1] + vecAy[i,j+1,1] + vecAy[i-1,j+1,1] )
            vecAy_yav    = 0.25*( vecAy[i,j,2] + vecAy[i-1,j,2] + vecAy[i,j+1,2] + vecAy[i-1,j+1,2] )

            dTdx = vecAx[i,j,1]*dTdxi + vecAy_xav*dTdeta
            dTdy = vecAx[i,j,2]*dTdxi + vecAy_yav*dTdeta

            volume_temp = 0.5*(volume[i,j]+volume[i-1,j])
            E_hat[i,j] = D*(vecAx[i,j,1]*dTdx + vecAx[i,j,2]*dTdy)/volume_temp
        end
    end

    F_hat = zeros(cellxmax,cellymax+1)

    Threads.@threads for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1
			T_av = 0.5*(Qbase[i,j]+Qbase[i,j-1])
			dTdxi = 0.25*((Qbase[i+1,j]-Qbase[i-1,j]) + (Qbase[i+1,j-1]-Qbase[i-1,j-1]))
			dTdeta = (Qbase[i,j]-Qbase[i,j-1])

            vecAx_xav    = 0.25*( vecAx[i,j,1] + vecAx[i,j-1,1] + vecAx[i+1,j,1] + vecAx[i+1,j-1,1] )
            vecAx_yav    = 0.25*( vecAx[i,j,2] + vecAx[i,j-1,2] + vecAx[i+1,j,2] + vecAx[i+1,j-1,2] )

            dTdx = vecAx_xav*dTdxi + vecAy[i,j,1]*dTdeta
            dTdy = vecAx_yav*dTdxi + vecAy[i,j,2]*dTdeta

            volume_temp = 0.5*(volume[i,j]+volume[i,j-1])
            F_hat[i,j] = D*(vecAy[i,j,1]*dTdx + vecAy[i,j,2]*dTdy)/volume_temp
        end
    end

    return E_hat,F_hat
end

function setup_RHS(cellxmax,cellymax,E_hat, F_hat)
	RHS = zeros(cellxmax,cellymax)

    Threads.@threads for i in 2:cellxmax-1
		for j in 2:cellymax-1
			RHS[i,j] = -((E_hat[i,j]-E_hat[i+1,j]) + (F_hat[i,j]-F_hat[i,j+1]))
		end
	end

	return RHS
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
gamma = 5/3
equations = IdealGlmMhdEquations1D(gamma)

"""
    initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Torrilhon's shock tube test case for one dimensional ideal MHD equations.
- Torrilhon (2003)
  Uniqueness conditions for Riemann problems of ideal magnetohydrodynamics
  [DOI: 10.1017/S0022377803002186](https://doi.org/10.1017/S0022377803002186)
"""
function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-1, 1.5], γ = 5/3, final time = 0.4
  rho = x[1] <= 0 ? 3.0 : 1.0
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = x[1] <= 0 ? 3.0 : 1.0
  B1 = 1.5
  B2 = x[1] <= 0 ? 1.0 : cos(1.5)
  B3 = x[1] <= 0 ? 0.0 : sin(1.5)
  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end
initial_condition = initial_condition_torrilhon_shock_tube

boundary_conditions = BoundaryConditionDirichlet(initial_condition)

surface_flux = flux_lax_friedrichs
volume_flux  = flux_central
basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = -1.0
coordinates_max =  1.5
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=7,
                n_cells_max=10_000,
                periodicity=false)

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

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

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

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_errors=(:l2_error_primitive,
                                                            :linf_error_primitive))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

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

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

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


text:
include("./util/Corrv3.jl")

parsed_args = parseCommandLine()

dict = init_dict(parsed_args)
main(dict; test=parsed_args["test"])

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text: #export


text:
# Abbreviation to help keep code short!
const AMSA = AbstractMultiScaleArray

Base.map!(f::F, m::AMSA, A0::AbstractArray, As::AbstractArray...) where {F} =
        broadcast!(f, m, A0, As...)
Base.map!(f::F, m::AMSA, A0, As...) where {F} =
        broadcast!(f, m, A0, As...)

struct AMSAStyle <: Broadcast.AbstractArrayStyle{Any} end
Broadcast.BroadcastStyle(::AMSAStyle, ::Broadcast.DefaultArrayStyle{0}) = AMSAStyle()
Broadcast.BroadcastStyle(::Broadcast.DefaultArrayStyle{0}, ::AMSAStyle) = AMSAStyle()
Broadcast.BroadcastStyle(::AMSAStyle, ::Broadcast.DefaultArrayStyle{N}) where N = Broadcast.DefaultArrayStyle{N}()
Broadcast.BroadcastStyle(::Type{<:AMSA}) = AMSAStyle()

@inline function Base.copy(bc::Broadcast.Broadcasted{<:AMSAStyle})
    first_amsa = find_amsa(bc)

    out = similar(first_amsa,Base.Broadcast._broadcast_getindex_eltype(bc))

    #=
    ElType = Base.Broadcast.combine_eltypes(bc.f, bc.args)
    if Base.isconcretetype(ElType)
        # We can trust it and defer to the simpler `copyto!`
        return copyto!(similar(bc, ElType), bc)
    end
    =#

    copyto!(out,bc)
    out
end

@inline function Base.copyto!(dest::AMSA, bc::Broadcast.Broadcasted{<:AMSAStyle})
    N = length(dest.nodes)
    for i in 1:N
        copyto!(dest.nodes[i], unpack(bc, i))
    end
    copyto!(dest.values,unpack(bc, nothing))
    dest
end

@inline function Base.copyto!(dest::AbstractMultiScaleArrayLeaf, bc::Broadcast.Broadcasted{<:AMSAStyle})
    copyto!(dest.values,unpack(bc,nothing))
    dest
end

# drop axes because it is easier to recompute
@inline unpack(bc::Broadcast.Broadcasted, i) = Broadcast.Broadcasted(bc.f, unpack_args(i, bc.args))
unpack(x,::Any) = x
unpack(x::AMSA, i) = x.nodes[i]
unpack(x::AMSA, ::Nothing) = x.values

@inline unpack_args(i, args::Tuple) = (unpack(args[1], i), unpack_args(i, Base.tail(args))...)
unpack_args(i, args::Tuple{Any}) = (unpack(args[1], i),)
unpack_args(::Any, args::Tuple{}) = ()

nnodes(A) = 0
nnodes(A::AMSA) = length(A.nodes)
nnodes(bc::Broadcast.Broadcasted) = _nnodes(bc.args)
nnodes(A, Bs...) = common_number(nnodes(A), _nnodes(Bs))

@inline _nnodes(args::Tuple) = common_number(nnodes(args[1]), _nnodes(Base.tail(args)))
_nnodes(args::Tuple{Any}) = nnodes(args[1])
_nnodes(args::Tuple{}) = 0

"`A = find_amsa(As)` returns the first AMSA among the arguments."
find_amsa(bc::Base.Broadcast.Broadcasted) = find_amsa(bc.args)
find_amsa(args::Tuple) = !isempty(args) && find_amsa(find_amsa(args[1]), Base.tail(args))
find_amsa(x) = x
find_amsa(a::AMSA, rest) = a
find_amsa(::Any, rest) = find_amsa(rest)

any_non_amsa(bc::Base.Broadcast.Broadcasted) = any_non_amsa(bc.args)
any_non_amsa(args::Tuple) = any_non_amsa(any_non_amsa(args[1]), Base.tail(args))
any_non_amsa(x::AMSA) = false
any_non_amsa(x::Number) = false
any_non_amsa(x::Any) = true
any_non_amsa(x::AbstractArray) = true
any_non_amsa(x::Bool, rest) = isempty(rest) ? x : x || any_non_amsa(rest)

## utils
common_number(a, b) =
    a == 0 ? b :
    (b == 0 ? a :
     (a == b ? a :
      throw(DimensionMismatch("number of nodes must be equal"))))

## Linear Algebra

function LinearAlgebra.ldiv!(A::LinearAlgebra.LU,b::AMSA)
    x = Array(b)
    ldiv!(A,x)
    b .= x
end
function LinearAlgebra.ldiv!(A::LinearAlgebra.QR,b::AMSA)
    x = Array(b)
    ldiv!(A,x)
    b .= x
end
function LinearAlgebra.ldiv!(A::LinearAlgebra.SVD,b::AMSA)
    x = Array(b)
    ldiv!(A,x)
    b .= x
end


text:
using Plots

# Horizon
T = 101
Tm = 51

# Time step
tf = 1.0
h = tf / (T - 1)

function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z
	z̄ = 0.0
	# x = range(x_start, stop = x_goal, length = T)
	x = circular_projection_range(x_start, stop = x_goal, length = T)
	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))
	return x, z
end

function circular_projection_range(start; stop=1.0, length=10)
	dist = stop - start
	θr = range(π, stop=0, length=length)
	r = start .+ dist * ((1 .+ cos.(θr))./2)
	return r
end

function initial_configuration(model;
    l_torso = 0.19,
    w_torso = 0.1025,
    offset = 0.025)
	q_body = zeros(6)

	# position
	q_body[3] = body_height

	# orientation
	mrp = MRP(RotZ(0.0))
	q_body[4:6] = [mrp.x; mrp.y; mrp.z]

	# feet positions (in body frame)
	q_f1 = [l_torso + offset; w_torso; 0.0]
	q_f2 = [l_torso - offset; -w_torso; 0.0]
	q_f3 = [-l_torso - offset; w_torso; 0.0]
	q_f4 = [-l_torso + offset; -w_torso; 0.0]

	return q_body, q_f1, q_f2, q_f3, q_f4
end

q_body, q_f1, q_f2, q_f3, q_f4 = initial_configuration(model)
visualize!(vis, model, [q_body], [q_f1], [q_f2], [q_f3], [q_f4])

# feet positions
pf1 = copy(q_f1)
pf2 = copy(q_f2)

pf3 = copy(q_f3)
pf4 = copy(q_f4)

strd = 2 * (pf1 - pf2)[1]

# q_shift1 = zeros(model.nq)
# q_shift1[1] = 0.5 * strd
# q_shift1[10] = strd
# q_shift1[13] = strd

body_shift = [0.5 * strd; 0.0; 0.0; 0.0; 0.0; 0.0]
foot_shift = [1.0 * strd; 0.0; 0.0]

qM_body = q_body + body_shift
qM_f1 = copy(q_f1)
qM_f2 = q_f2 + foot_shift
qM_f3 = q_f3 + foot_shift
qM_f4 = copy(q_f4)

qT_body = qM_body + body_shift
qT_f1 = copy(qM_f1) + foot_shift
qT_f2 = copy(qM_f2)
qT_f3 = copy(qM_f3)
qT_f4 = copy(qM_f4) + foot_shift

q_body_ref = linear_interpolation(q_body, qT_body, T)#[q_body, linear_interpolation(q_body, qM_body, Tm - 1)..., linear_interpolation(qM_body, qT_body, Tm - 1)...]
q_f1_ref = [q_f1, linear_interpolation(q_f1, qM_f1, Tm - 1)..., linear_interpolation(qM_f1, qT_f1, Tm - 1)...]
q_f2_ref = [q_f2, linear_interpolation(q_f2, qM_f2, Tm - 1)..., linear_interpolation(qM_f2, qT_f2, Tm - 1)...]
q_f3_ref = [q_f3, linear_interpolation(q_f3, qM_f3, Tm - 1)..., linear_interpolation(qM_f3, qT_f3, Tm - 1)...]
q_f4_ref = [q_f4, linear_interpolation(q_f4, qM_f4, Tm - 1)..., linear_interpolation(qM_f4, qT_f4, Tm - 1)...]

v_body_ref = (qT_body - q_body) / tf

visualize!(vis, model, q_body_ref, q_f1_ref, q_f2_ref, q_f3_ref, q_f4_ref,
    Δt = h)

T_fix = 5
zh = 0.05
xf1_el, zf1_el = ellipse_traj(pf1[1], pf1[1] + strd, zh, Tm - T_fix)
xf1 = [[pf1[1] for t = 1:Tm + T_fix]..., xf1_el[2:end]...]
zf1 = [[pf1[3] for t = 1:Tm + T_fix]..., zf1_el[2:end]...]
pf1_ref = [[xf1[t]; pf1[2];  zf1[t]] for t = 1:T]
f1_contact = [t < (Tm + T_fix) ? 1 : 0 for t = 1:T]

xf4_el, zf4_el = ellipse_traj(pf4[1], pf4[1] + strd, zh, Tm - T_fix)
xf4 = [[pf4[1] for t = 1:Tm + T_fix]..., xf4_el[2:end]...]
zf4 = [[pf4[3] for t = 1:Tm + T_fix]..., zf4_el[2:end]...]
pf4_ref = [[xf4[t]; pf4[2]; zf4[t]] for t = 1:T]
f4_contact = [t < (Tm + T_fix) ? 1 : 0 for t = 1:T]

xf2_el, zf2_el = ellipse_traj(pf2[1], pf2[1] + strd, zh, Tm - T_fix)
xf2 = [[xf2_el[1] for t = 1:T_fix]..., xf2_el..., [xf2_el[end] for t = 1:Tm-1 + T_fix]...]
zf2 = [[zf2_el[1] for t = 1:T_fix]..., zf2_el..., [zf2_el[end] for t = 1:Tm-1 + T_fix]...]
pf2_ref = [[xf2[t]; pf2[2]; zf2[t]] for t = 1:T]
f2_contact = [t <= T_fix ? 1 : (t > Tm ? 1 : 0) for t = 1:T]

xf3_el, zf3_el = ellipse_traj(pf3[1], pf3[1] + strd, zh, Tm - T_fix)
xf3 = [[xf3_el[1] for t = 1:T_fix]..., xf3_el..., [xf3_el[end] for t = 1:Tm-1]...]
zf3 = [[zf3_el[1] for t = 1:T_fix]..., zf3_el..., [zf3_el[end] for t = 1:Tm-1]...]
pf3_ref = [[xf3[t]; pf3[2]; zf3[t]] for t = 1:T]
f3_contact = [t <= T_fix ? 1 : (t > Tm ? 1 : 0) for t = 1:T]

tr = range(0, stop = tf, length = T)
plot(tr, hcat(pf1_ref...)', labels = "")
plot!(tr, hcat(pf4_ref...)', labels = "")

plot!(tr, hcat(pf2_ref...)', labels = "")
plot!(tr, hcat(pf3_ref...)', labels = "")

visualize!(vis, model, q_body_ref,
    pf1_ref, pf2_ref, pf3_ref, pf4_ref,
    Δt = h)

p_fr = plot(tr, f1_contact, linetype = :steppost, color = :black, label = "FR", xaxis = false, x_ticks = false, yaxis = false, yticks = false)
p_rl = plot(tr, f4_contact, linetype = :steppost, color = :black, label = "RL", xaxis = false, x_ticks = false, yaxis = false, yticks = false)

p_fl = plot(tr, f2_contact, linetype = :steppost, color = :black, label = "FL", xaxis = false, x_ticks = false, yaxis = false, yticks = false)
p_rr = plot(tr, f3_contact, linetype = :steppost, color = :black, label = "RR", yaxis = false, yticks = false, xaxis = "time (s)")

plot(p_fr, p_rl, p_fl, p_rr, layout = (4, 1))

contact_modes = [[f1_contact[t];
                  f2_contact[t];
                  f3_contact[t];
                  f4_contact[t]] for t = 1:T]

contact_modes_2 = [contact_modes..., contact_modes[2:end]...]

plot(hcat(contact_modes_2...)', linetype = :steppost, label = "")

T2 = 2 * T - 1
tf2 = 2 * tf
tr2 = range(0, stop = tf2, length = T2)
pf1_ref_2 = [pf1_ref..., [p + foot_shift for p in pf1_ref[2:end]]...]
pf2_ref_2 = [pf2_ref..., [p + foot_shift for p in pf2_ref[2:end]]...]
pf3_ref_2 = [pf3_ref..., [p + foot_shift for p in pf3_ref[2:end]]...]
pf4_ref_2 = [pf4_ref..., [p + foot_shift for p in pf4_ref[2:end]]...]
qT_body_2 = q_body + 4 * body_shift
q_body_ref_2 = linear_interpolation(q_body, qT_body_2, T2)#[q_body, linear_interpolation(q_body, qM_body, Tm - 1)..., linear_interpolation(qM_body, qT_body, Tm - 1)...]

plot(tr2, hcat(pf1_ref_2...)', labels = "")
plot!(tr2, hcat(pf2_ref_2...)', labels = "")
plot!(tr2, hcat(pf3_ref_2...)', labels = "")
plot!(tr2, hcat(pf4_ref_2...)', labels = "")

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
`SqExponentialKernel([ρ=1.0])`

The squared exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ²‖x-y‖²)
```
See also [`ExponentialKernel`](@ref) for a
related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
"""
struct SqExponentialKernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::SqEuclidean
    function SqExponentialKernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
        return new{T,Tr}(transform,SqEuclidean())
    end
end

@inline kappa(κ::SqExponentialKernel, d²::Real) = exp(-d²)
@inline iskroncompatible(::SqExponentialKernel) = true

## Aliases ##
const RBFKernel = SqExponentialKernel
const GaussianKernel = SqExponentialKernel

"""
`ExponentialKernel([ρ=1.0])`
The exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ‖x-y‖)
```
"""
struct ExponentialKernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::Euclidean
    function ExponentialKernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
        return new{T,Tr}(transform,Euclidean())
    end
end

@inline kappa(κ::ExponentialKernel, d::Real) = exp(-d)
@inline iskroncompatible(::ExponentialKernel) = true

## Alias ##
const LaplacianKernel = ExponentialKernel

"""
`GammaExponentialKernel([ρ=1.0,[gamma=2.0]])`
The γ-exponential kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = exp(-ρ^(2γ)‖x-y‖^(2γ))
```
"""
struct GammaExponentialKernel{T,Tr,Tᵧ<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::SqEuclidean
    γ::Tᵧ
    function GammaExponentialKernel{T,Tr,Tᵧ}(transform::Tr,γ::Tᵧ) where {T,Tr<:Transform,Tᵧ<:Real}
        return new{T,Tr,Tᵧ}(transform,SqEuclidean(),γ)
    end
end

params(k::GammaExponentialKernel) = (params(transform),γ)
opt_params(k::GammaExponentialKernel) = (opt_params(transform),γ)

function GammaExponentialKernel(ρ::T₁=1.0,gamma::T₂=2.0) where {T₁<:Real,T₂<:Real}
    @check_args(GammaExponentialKernel, gamma, gamma >= zero(T₂), "gamma > 0")
    GammaExponentialKernel{T₁,ScaleTransform{T₁},T₂}(ScaleTransform(ρ),gamma)
end

function GammaExponentialKernel(ρ::AbstractVector{T₁},gamma::T₂=2.0) where {T₁<:Real,T₂<:Real}
    @check_args(GammaExponentialKernel, gamma, gamma >= zero(T₂), "gamma > 0")
    GammaExponentialKernel{T₁,ARDTransform{T₁,length(ρ)},T₂}(ARDTransform(ρ),gamma)
end

function GammaExponentialKernel(t::Tr,gamma::T₁=2.0) where {Tr<:Transform,T₁<:Real}
    @check_args(GammaExponentialKernel, gamma, gamma >= zero(T₁), "gamma > 0")
    GammaExponentialKernel{eltype(Tr),Tr,T₁}(t,gamma)
end

@inline kappa(κ::GammaExponentialKernel, d²::Real) = exp(-d²^κ.γ)
@inline iskroncompatible(::GammaExponentialKernel) = true


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors, FileIO
if !isdefined(:runtests)
    window = glscreen()
end

description = """
Demonstration of loding different image formats. The blank images and the ones
showing checkmarks should all be green. As you can see, you can also quite easily
load and display animated images.
"""

# a few helper functions to generate images

const NColor{N, T} = Colorant{T, N}
fillcolor(::Type{T}) where {T <: NColor{4}} = T(0,1,0,1)
fillcolor(::Type{T}) where {T <: NColor{3}} = T(0,1,0)

# create different images with different color types (not an exhaustive list of supported types)
arrays = map((RGBA{N0f8}, RGBA{Float32}, RGB{N0f8}, RGB{Float32}, BGRA{N0f8}, BGR{Float32})) do C
     C[fillcolor(C) for x=1:45,y=1:45]
 end
# load a few images from the asset folder with FileIO.load (that's what loadasset calls)
loaded_imgs = map(x->loadasset("test_images", x), readdir(assetpath("test_images")))

# combine them all into one array and add an animated gif and a few other images

# backward compatible imconvert
imconvert(im) = im
if isdefined(:Image)
    function imconvert(im::Image)
        if ndims(im) == 2
            convert(Array{eltype(im), ndims(im)}, im)
        else
            permutedims(im.data, (2, 1, 3))
        end
    end
end
x = Any[
    arrays..., loaded_imgs...,
    loadasset("kittens-look.gif"),
    loadasset("mario", "stand", "right.png"),
    loadasset("mario", "jump", "left.gif"),
]
x = map(imconvert, x)

# visualize all images and convert the array to be a vector of element type context
# This shouldn't be necessary, but it seems map is not able to infer the type alone
images = convert(Vector{Context}, map(visualize, x))
# make it a grid
images = reshape(images, (4, 4))
# GLVisualize offers a few helpers to visualize arrays of render objects
# spaced out as the underlying array. So this will create a grid whereas every
# item is 128x128x128 pixels big
img_vis = visualize(images, scale = Vec3f0(128))
_view(img_vis, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
module STABLE

## %% Load packages
using Dieter
# import Dieter: initialise_data_file_dict!, check_files_exist
import Dieter: parse_file, parse_nodes!, parse_base_technologies!, parse_storages!, parse_load!, parse_availibility!
import Dieter: initialise_set_relation_data!, parse_set_relations!,parse_arcs!, calc_base_parameters!, parse_extensions!
import Dieter: dvalmatch, dkeymatch, split_df_tuple

using JuMP
# import MathOptInterface
# import CPLEX

using DataFrames
using DataFramesMeta
import DBInterface
import SQLite
import CSV
import TimeSeries
import Dates
import Serialization
# using Tables
# import XLSX

import OrderedCollections: OrderedDict

include("structures.jl")
include("util.jl")
# include("datamodels.jl")

end


text:
###
# SphereTrav
###

"""
    SphereTrav(A::AbstractMatrix)

is an anlogue of `DiagTrav` but for coefficients stored according to
FastTransforms.jl spherical harmonics layout
"""
struct SphereTrav{T, AA<:AbstractMatrix{T}} <: AbstractBlockVector{T}
    matrix::AA
    function SphereTrav{T, AA}(matrix::AA) where {T,AA<:AbstractMatrix{T}}
        n,m = size(matrix)
        m == 2n-1 || throw(ArgumentError("size must match"))
        new{T,AA}(matrix)
    end
end

SphereTrav{T}(matrix::AbstractMatrix{T}) where T = SphereTrav{T,typeof(matrix)}(matrix)
SphereTrav(matrix::AbstractMatrix{T}) where T = SphereTrav{T}(matrix)

axes(A::SphereTrav) = (blockedrange(range(1; step=2, length=size(A.matrix,1))),)

function getindex(A::SphereTrav, K::Block{1})
    k = Int(K)
    m = size(A.matrix,1)
    st = stride(A.matrix,2)
    # nonnegative terms
    p = A.matrix[range(k; step=2*st-1, length=k)]
    k == 1 && return p
    # negative terms
    n = A.matrix[range(k+st-1; step=2*st-1, length=k-1)]
    [reverse!(n); p]
end

getindex(A::SphereTrav, k::Int) = A[findblockindex(axes(A,1), k)]

"""
    RealSphereTrav(A::AbstractMatrix)

    takes coefficients as provided by the spherical harmonics layout of FastTransforms.jl and
    makes them accessible sorted such that in each block the m=0 entries are always in first place,
    followed by alternating sin and cos terms of increasing |m|.
"""
struct RealSphereTrav{T, AA<:AbstractMatrix{T}} <: AbstractBlockVector{T}
    matrix::AA
    function RealSphereTrav{T, AA}(matrix::AA) where {T,AA<:AbstractMatrix{T}}
        n,m = size(matrix)
        m == 2n-1 || throw(ArgumentError("size must match"))
        new{T,AA}(matrix)
    end
end

RealSphereTrav{T}(matrix::AbstractMatrix{T}) where T = RealSphereTrav{T,typeof(matrix)}(matrix)
RealSphereTrav(matrix::AbstractMatrix{T}) where T = RealSphereTrav{T}(matrix)

axes(A::RealSphereTrav) = (blockedrange(range(1; step=2, length=size(A.matrix,1))),)

function getindex(A::RealSphereTrav, K::Block{1})
    k = Int(K)
    m = size(A.matrix,1)
    st = stride(A.matrix,2)
    # nonnegative terms
    p = A.matrix[range(k; step=2*st-1, length=k)]
    k == 1 && return p
    # negative terms
    n = A.matrix[range(k+st-1; step=2*st-1, length=k-1)]
    interlace(p,n)
end

getindex(A::RealSphereTrav, k::Int) = A[findblockindex(axes(A,1), k)]

for Typ in (:SphereTrav,:RealSphereTrav)
    @eval function resize!(A::$Typ, K::Block{1})
        k = Int(K)
        $Typ(A.matrix[1:k, 1:2k-1])
    end
end


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

ArgParse.parse_item(::Type{Vector{String}}, x::AbstractString) = split(x, ",")
ArgParse.parse_item(::Type{Union{Int, Nothing}}, x::AbstractString) = x == "nothing" ? nothing : parse(Int, x)

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--min_nv"
            help = "Minimum number of vertices a network had to have to be processed."
            arg_type = Union{Int, Nothing}
            default = nothing
        "--max_nv"
            help = "Maximum number of vertices a network had to have to be processed."
            arg_type = Union{Int, Nothing}
            default = nothing
        "--checks"
            help = "Name of the checks to apply on the database."
            arg_type = Vector{String}
            default = ["chordality", "connectivity", "self_loops", "index_clique", "source_graph", "serialize_graph", "serialize_network", "basic_feature"]
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    check_sanity(db, parsed_args["checks"];
                 min_nv=parsed_args["min_nv"], max_nv=parsed_args["max_nv"])
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
# Draw the Mandelbrot set
#
# Inspired by the gtk-server fractal demo, by Peter van Eerten
#
# Requires the ColorSchemes package

using GtkServer, ColorSchemes

fc(c, z) = z^2 + c

function mandelbrot(c, maxiters, radius)
    iters = 0
    z = 0
    r = 0.0
    while r < radius
        iters += 1
        z = fc(c, z)
        r = abs(z)^2
        iters == maxiters && break
    end
    return iters
end

# The x and y range to plot can be specified as tuples. Also, the maximum number of iterations,
# the radius at which the point is assumed to have escaped, and the palette can also be specified.
function demo(; x = (-2.3, 1), y = (-1.2, 1.2), maxiters = 50, radius = 5.0, palette = :dense)
    # Initialize GtkServer
    gtk_init()
    # Initialize GTK itself
    gtk_send("gtk_init NULL NULL")
    # Create the top-level window, give it a title, size it, center it, and make it non-resizeable
    win = gtk_send("gtk_window_new 0")
    gtk_send("gtk_window_set_title $win 'Julia GTK Mandelbrot Generator'")
    gtk_send("gtk_widget_set_size_request $win 800 650")  # width, height
    gtk_send("gtk_window_set_position $win 1")
    gtk_send("gtk_window_set_resizable $win 0")
    # Create widget to display image
    image = gtk_send("gtk_image_new")
    ebox = gtk_send("gtk_event_box_new")
    gtk_send("gtk_container_add $ebox $image")
    # Create a separator
    sep = gtk_send("gtk_hseparator_new")
    # Create buttons
    draw_btn = gtk_send("gtk_button_new_with_label Draw")
    gtk_send("gtk_widget_set_size_request $draw_btn 75 30")
    clear_btn = gtk_send("gtk_button_new_with_label Clear")
    gtk_send("gtk_widget_set_size_request $clear_btn 75 30")
    exit_btn = gtk_send("gtk_button_new_with_label Exit")
    gtk_send("gtk_widget_set_size_request $exit_btn 75 30")
    # place widgets using boxes
    hbox = gtk_send("gtk_hbox_new 0 0")
    gtk_send("gtk_box_pack_start $hbox $draw_btn 0 0 1")
    gtk_send("gtk_box_pack_start $hbox $clear_btn 0 0 1")
    gtk_send("gtk_box_pack_end $hbox $exit_btn 0 0 1")
    vbox = gtk_send("gtk_vbox_new 0 0")
    gtk_send("gtk_box_pack_start $vbox $ebox 0 0 1")
    gtk_send("gtk_box_pack_start $vbox $sep 0 0 1")
    gtk_send("gtk_box_pack_end $vbox $hbox 0 0 1")
    gtk_send("gtk_container_add $win $vbox")
    # Show all widgets
    gtk_send("gtk_widget_show_all $win")
    # Create an image to display the fractal
    gdkwin = gtk_send("gtk_widget_get_parent_window $image")
    pix = gtk_send("gdk_pixmap_new $gdkwin 799 600 -1")
    gcn = gtk_send("gdk_gc_new $pix")
    gtk_send("gtk_image_set_from_pixmap $image $pix NULL")
    color = gtk_send("gtk_frame_new NULL")
    gtk_send("gdk_color_parse #ffffff $color")
    gtk_send("gdk_gc_set_rgb_bg_color $gcn $color")
    gtk_send("gdk_gc_set_rgb_fg_color $gcn $color")
    gtk_send("gdk_draw_rectangle $pix $gcn 1 0 0 799 600")
    gtk_send("gdk_color_parse #000000 $color")
    gtk_send("gdk_gc_set_rgb_fg_color $gcn $color")
    gtk_send("gtk_widget_queue_draw $image")
    # show all widgets
    gtk_send("gtk_widget_show_all $win")
    gtk_send("gtk_server_callback update")

    # Mandelbrot variables
    X = range(x[1], x[2], length = 799)
    Y = range(y[1], y[2], length = 600)
    maxiters = 50
    radius = 5.0

    # choose palette
    cm = colorschemes[palette]

    # cache colors
    colorcache = Dict{Int, String}()

    # main loop
    loop = true
    px = 0
    py = 0
    while loop
        result = gtk_send("gtk_server_callback wait")
        if result == draw_btn
            for yy in Y
                px = 0
                py += 1
                loop == false && break
                for xx in X
                    px += 1
                    iters = mandelbrot(Complex(xx, yy), maxiters, radius)
                    if haskey(colorcache, iters)
                        pcolor = colorcache[iters]
                    else
                        c = get(cm, iters/maxiters)
                        cr = string(floor(Int,256*c.r), base=16, pad=2)
                        cg = string(floor(Int,256*c.g), base=16, pad=2)
                        cb = string(floor(Int,256*c.b), base=16, pad=2)
                        pcolor = "#$cr$cg$cb"
                        push!(colorcache, iters => pcolor)
                    end
                    gtk_send("gdk_color_parse $pcolor $color")
                    gtk_send("gdk_gc_set_rgb_fg_color $gcn $color")
                    gtk_send("gdk_draw_point $pix $gcn $px $py")
                    result = gtk_send("gtk_server_callback update")
                    if result == exit_btn || result == win
                        loop = false
                        break
                    end
                end
                gtk_send("gtk_widget_queue_draw $image")
                gtk_send("gtk_server_callback update")
            end
        elseif result == clear_btn
            gtk_send("gdk_color_parse #ffffff $color")
            gtk_send("gdk_gc_set_rgb_fg_color $gcn $color")
            gtk_send("gdk_draw_rectangle $pix $gcn 1 0 0 799 600")
            gtk_send("gdk_color_parse #000000 $color")
            gtk_send("gdk_gc_set_rgb_fg_color $gcn $color")
            gtk_send("gtk_widget_queue_draw $image")
        elseif result == exit_btn || result == win
            break
        end
    end
    gtk_exit()
end


text:
function regulargrid2d(mins, maxs, ns, dz)
	linearindex = (i1, i2)->i2 + ns[2] * (i1 - 1)
	coords = Array{Float64}(undef, 2, prod(ns))
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	dx = xs[2] - xs[1]
	dy = ys[2] - ys[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, 2 * prod(ns) - ns[1] - ns[2])
	areasoverlengths = Array{Float64}(undef, 2 * prod(ns) - ns[1] - ns[2])
	volumes = Array{Float64}(undef, 0)
	for i1 = 1:ns[1]
		for i2 = 1:ns[2]
			push!(volumes, dx * dy * dz)
			coords[1, linearindex(i1, i2)] = xs[i1]
			coords[2, linearindex(i1, i2)] = ys[i2]
			if i1 < ns[1]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1 + 1, i2)
				areasoverlengths[j] = dy * dz / dx
				j += 1
			end
			if i2 < ns[2]
				neighbors[j] = linearindex(i1, i2)=>linearindex(i1, i2 + 1)
				areasoverlengths[j] = dx * dz / dy
				j += 1
			end
		end
	end
	return coords, neighbors, areasoverlengths, volumes
end

function regulargrid3d(mins, maxs, ns)
	linearindex = (i1, i2, i3)->i3 + ns[3] * (i2 - 1) + ns[3] * ns[2] * (i1 - 1)
	coords = Array{Float64}(undef, 3, prod(ns))
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	zs = range(mins[3]; stop=maxs[3], length=ns[3])
	dx = xs[2] - xs[1]
	dy = ys[2] - ys[1]
	dz = zs[2] - zs[1]
	j = 1
	neighbors = Array{Pair{Int, Int}}(undef, 3 * prod(ns) - ns[1] * ns[2] - ns[1] * ns[3] - ns[2] * ns[3])
	areasoverlengths = Array{Float64}(undef, 3 * prod(ns) - ns[1] * ns[2] - ns[1] * ns[3] - ns[2] * ns[3])
	volumes = fill(dx * dy * dz, size(coords, 2))
	for i1 = 1:ns[1]
		for i2 = 1:ns[2]
			for i3 = 1:ns[3]
				coords[1, linearindex(i1, i2, i3)] = xs[i1]
				coords[2, linearindex(i1, i2, i3)] = ys[i2]
				coords[3, linearindex(i1, i2, i3)] = zs[i3]
				if i1 < ns[1]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1 + 1, i2, i3)
					areasoverlengths[j] = dy * dz / dx
					j += 1
				end
				if i2 < ns[2]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1, i2 + 1, i3)
					areasoverlengths[j] = dx * dz / dy
					j += 1
				end
				if i3 < ns[3]
					neighbors[j] = linearindex(i1, i2, i3)=>linearindex(i1, i2, i3 + 1)
					areasoverlengths[j] = dx * dy / dz
					j += 1
				end
			end
		end
	end
	return coords, neighbors, areasoverlengths, volumes
end

function darcy_velocity(h, Ks, coords, ns)
	@assert length(ns) == 2
	@show size(Ks)
	@show size(ns)
	hr = reshape(h, ns...)
	vr = zeros(2, ns...)
	@show size(vr)
	@show ns
	coordsr = reshape(coords, 2, reverse(ns)...)
	@show size(coordsr)
	@show size(hr)
	for i = 2:ns[1] - 1
		for j = 2:ns[2] - 1
			#=
			if coordsr[1, j, i + 1, j] - coordsr[1, i - 1, j] == 0
				@show i, j
				@show coordsr[1, i + 1, j], coordsr[1, i - 1, j]
			end
			=#
			vr[1, i, j] = Ks[i, j] * (hr[i + 1, j] - hr[i - 1, j]) / (coordsr[1, j, i + 1] - coordsr[1, j, i - 1])
			vr[2, i, j] = Ks[i, j] * (hr[i, j + 1] - hr[i, j - 1]) / (coordsr[2, j + 1, i] - coordsr[2, j - 1, i])
		end
	end
	return vr
end