text:
## initialize!

@muladd function initialize!(nlsolver::NLSolver{<:NLNewton,false}, integrator)
  @unpack dt = integrator
  @unpack cache = nlsolver

  cache.invγdt = inv(dt * nlsolver.γ)
  cache.tstep = integrator.t + nlsolver.c * dt

  nothing
end

@muladd function initialize!(nlsolver::NLSolver{<:NLNewton,true}, integrator)
  @unpack u,uprev,t,dt,opts = integrator
  @unpack cache = nlsolver
  @unpack weight = cache

  cache.invγdt = inv(dt * nlsolver.γ)
  cache.tstep = integrator.t + nlsolver.c * dt
  calculate_residuals!(weight, fill!(weight, one(eltype(u))), uprev, u,
                       opts.abstol, opts.reltol, opts.internalnorm, t)

  nothing
end

## compute_step!

"""
    compute_step!(nlsolver::NLSolver{<:NLNewton}, integrator)

Compute next iterate of numerically stable modified Newton iteration
that is specialized for implicit methods.

Please check
https://github.com/SciML/DiffEqDevMaterials/blob/master/newton/output/main.pdf
for more details.

# References

M.E.Hoseaa and L.F.Shampine, "Analysis and implementation of TR-BDF2",
Applied Numerical Mathematics, Volume 20, Issues 1–2, February 1996, Pages
21-37.
[doi:10.1016/0168-9274(95)00115-8](https://doi.org/10.1016/0168-9274(95)00115-8).

Ernst Hairer and Gerhard Wanner, "Solving Ordinary Differential
Equations II, Springer Series in Computational Mathematics. ISBN
978-3-642-05221-7. Section IV.8.
[doi:10.1007/978-3-642-05221-7](https://doi.org/10.1007/978-3-642-05221-7).
"""
@muladd function compute_step!(nlsolver::NLSolver{<:NLNewton,false}, integrator)
  @unpack uprev,t,p,dt,opts = integrator
  @unpack z,tmp,γ,α,cache = nlsolver
  @unpack tstep,W,invγdt = cache

  f = nlsolve_f(integrator)
  isdae = f isa DAEFunction

  if isdae
    # not all predictors are uprev, for other forms of predictors, defined in u₀
    if isdefined(integrator.cache, :u₀)
      ustep = @.. integrator.cache.u₀ + z
    else
      ustep = @.. uprev + z
    end
    dustep = @. (tmp + α * z) * invγdt
    ztmp = f(dustep, ustep, p, t)
  else
    mass_matrix = integrator.f.mass_matrix
    if nlsolver.method === COEFFICIENT_MULTISTEP
      ustep = z
      # tmp = outertmp ./ hγ
      if mass_matrix === I
        ztmp = tmp .+ f(z, p, tstep) .- (α * invγdt) .* z
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        ztmp = tmp .+ f(z, p, tstep) .- (mass_matrix * z) .* (α * invγdt)
      end
    else
      ustep = @. tmp + γ * z
      if mass_matrix === I
        ztmp = (dt .* f(ustep, p, tstep) .- z) .* invγdt
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        ztmp = (dt .* f(ustep, p, tstep) .- mass_matrix * z) .* invγdt
      end
    end
  end

  if DiffEqBase.has_destats(integrator)
    integrator.destats.nf += 1
  end

  # update W
  if W isa DiffEqBase.AbstractDiffEqLinearOperator
    W = update_coefficients!(W, ustep, p, tstep)
  end

  dz = _reshape(W \ _vec(ztmp), axes(ztmp))
  if DiffEqBase.has_destats(integrator)
    integrator.destats.nsolve += 1
  end

  atmp = calculate_residuals(dz, uprev, ustep, opts.abstol, opts.reltol, opts.internalnorm, t)
  ndz = opts.internalnorm(atmp, t)
  # NDF and BDF are special because the truncation error is directly
  # propertional to the total displacement.
  if integrator.alg isa QNDF
    ndz *= error_constant(integrator, alg_order(integrator.alg))
  end

  # compute next iterate
  nlsolver.ztmp = z .- dz

  ndz
end

@muladd function compute_step!(nlsolver::NLSolver{<:NLNewton,true}, integrator)
  @unpack uprev,t,p,dt,opts = integrator
  @unpack z,tmp,ztmp,γ,α,iter,cache = nlsolver
  @unpack W_γdt,ustep,tstep,k,atmp,dz,W,new_W,invγdt,linsolve,weight = cache

  f = nlsolve_f(integrator)
  isdae = f isa DAEFunction

  if DiffEqBase.has_destats(integrator)
    integrator.destats.nf += 1
  end

  if isdae
    @.. ztmp = (tmp + α * z) * invγdt
    # not all predictors are uprev, for other forms of predictors, defined in u₀
    if isdefined(integrator.cache, :u₀)
      @.. ustep = integrator.cache.u₀ + z
    else
      @.. ustep = uprev + z
    end
    f(k, ztmp, ustep, p, tstep)
    b = _vec(k)
  else
    mass_matrix = integrator.f.mass_matrix
    if nlsolver.method === COEFFICIENT_MULTISTEP
      ustep = z
      f(k, z, p, tstep)
      if mass_matrix === I
        @.. ztmp = tmp + k - (α * invγdt) * z
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        mul!(_vec(ztmp), mass_matrix, _vec(z))
        @.. ztmp = tmp + k - (α * invγdt) * ztmp
      end
    else
      @.. ustep = tmp + γ * z
      f(k, ustep, p, tstep)
      if mass_matrix === I
        @.. ztmp = (dt * k - z) * invγdt
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        mul!(_vec(ztmp), mass_matrix, _vec(z))
        @.. ztmp = (dt * k - ztmp) * invγdt
      end
    end
    b = _vec(ztmp)
  end

  # update W
  if W isa DiffEqBase.AbstractDiffEqLinearOperator
    update_coefficients!(W, ustep, p, tstep)
  end

  if integrator.opts.adaptive
    reltol = integrator.opts.reltol
  else
    reltol = eps(eltype(dz))
  end

  if iter == 1 && new_W
    linsolve = LinearSolve.set_A(linsolve,W)
  end
  linres = dolinsolve(integrator, linsolve; b = _vec(b), u = _vec(dz), reltol = reltol)
  cache.linsolve = linres.cache

  if DiffEqBase.has_destats(integrator)
    integrator.destats.nsolve += 1
  end

  # relaxed Newton
  # Diagonally Implicit Runge-Kutta Methods for Ordinary Differential
  # Equations. A Review, by Christopher A. Kennedy and Mark H. Carpenter
  # page 54.
  if isdae
    γdt = α * invγdt
  else
    γdt = γ * dt
  end

  !(W_γdt ≈ γdt) && (rmul!(dz, 2/(1 + γdt / W_γdt)))

  calculate_residuals!(atmp, dz, uprev, ustep, opts.abstol, opts.reltol, opts.internalnorm, t)
  ndz = opts.internalnorm(atmp, t)
  # NDF and BDF are special because the truncation error is directly
  # propertional to the total displacement.
  if integrator.alg isa QNDF
    ndz *= error_constant(integrator, alg_order(integrator.alg))
  end

  # compute next iterate
  @.. ztmp = z - dz

  ndz
end

@muladd function compute_step!(nlsolver::NLSolver{<:NLNewton,true,<:Array}, integrator)
  @unpack uprev,t,p,dt,opts = integrator
  @unpack z,tmp,ztmp,γ,α,iter,cache = nlsolver
  @unpack W_γdt,ustep,tstep,k,atmp,dz,W,new_W,invγdt,linsolve,weight = cache
  f = nlsolve_f(integrator)
  isdae = f isa DAEFunction

  if DiffEqBase.has_destats(integrator)
    integrator.destats.nf += 1
  end

  if isdae
    @inbounds @simd ivdep for i in eachindex(z)
      ztmp[i] = (tmp[i] + α * z[i]) * invγdt
    end
    if isdefined(integrator.cache, :u₀)
      @inbounds @simd ivdep for i in eachindex(z)
        ustep[i] = integrator.cache.u₀[i] + z[i]
      end
      #@.. ustep = integrator.cache.u₀ + z
    else
      @inbounds @simd ivdep for i in eachindex(z)
        ustep[i] = uprev[i] + z[i]
      end
    end
    f(k, ztmp, ustep, p, tstep)
    b = _vec(k)
  else
    mass_matrix = integrator.f.mass_matrix
    if nlsolver.method === COEFFICIENT_MULTISTEP
      ustep = z
      f(k, z, p, tstep)
      if mass_matrix === I
        @inbounds @simd ivdep for i in eachindex(z)
          ztmp[i] = tmp[i] + k[i] - (α * invγdt) * z[i]
        end
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        mul!(_vec(ztmp), mass_matrix, _vec(z))

        @inbounds @simd ivdep for i in eachindex(z)
          ztmp[i] = tmp[i] + k[i] - (α * invγdt) * ztmp[i]
        end
      end
    else

      @inbounds @simd ivdep for i in eachindex(z)
        ustep[i] = tmp[i] + γ * z[i]
      end
      f(k, ustep, p, tstep)
      if mass_matrix === I
        @inbounds @simd ivdep for i in eachindex(z)
          ztmp[i] = (dt * k[i] - z[i]) * invγdt
        end
      else
        update_coefficients!(mass_matrix, ustep, p, tstep)
        mul!(_vec(ztmp), mass_matrix, _vec(z))
        @inbounds @simd ivdep for i in eachindex(z)
          ztmp[i] = (dt * k[i] - ztmp[i]) * invγdt
        end
      end
    end
    b = _vec(ztmp)
  end

  # update W
  if W isa DiffEqBase.AbstractDiffEqLinearOperator
    update_coefficients!(W, ustep, p, tstep)
  end

  if integrator.opts.adaptive
    reltol = integrator.opts.reltol
  else
    reltol = eps(eltype(dz))
  end

  if iter == 1 && new_W
    linsolve = LinearSolve.set_A(linsolve,W)
  end
  linres = dolinsolve(integrator, linsolve; b = _vec(b), u = _vec(dz), reltol = reltol)
  cache.linsolve = linres.cache

  if DiffEqBase.has_destats(integrator)
    integrator.destats.nsolve += 1
  end

  # relaxed Newton
  # Diagonally Implicit Runge-Kutta Methods for Ordinary Differential
  # Equations. A Review, by Christopher A. Kennedy and Mark H. Carpenter
  # page 54.
  if isdae
    γdt = α * invγdt
  else
    γdt = γ * dt
  end

  !(W_γdt ≈ γdt) && (rmul!(dz, 2/(1 + γdt / W_γdt)))

  calculate_residuals!(atmp, dz, uprev, ustep, opts.abstol, opts.reltol, opts.internalnorm, t)
  ndz = opts.internalnorm(atmp, t)
  # NDF and BDF are special because the truncation error is directly
  # propertional to the total displacement.
  if integrator.alg isa QNDF
    ndz *= error_constant(integrator, alg_order(integrator.alg))
  end

  # compute next iterate
  @inbounds @simd ivdep for i in eachindex(z)
    ztmp[i] = z[i] - dz[i]
  end

  ndz
end

## resize!

function Base.resize!(nlcache::NLNewtonCache, ::AbstractNLSolver, integrator, i::Int)
  resize!(nlcache.ustep, i)
  resize!(nlcache.k, i)
  resize!(nlcache.atmp, i)
  resize!(nlcache.dz, i)
  resize!(nlcache.du1, i)
  if nlcache.jac_config !== nothing
    resize_jac_config!(nlcache.jac_config, i)
  end
  resize!(nlcache.weight, i)

  # resize J and W (or rather create new ones of appropriate size and type)
  resize_J_W!(nlcache, integrator, i)

  nothing
end


text:
abstract type Component end
abstract type CurrentComponent <:Component end

struct Mapper{V,F<:Function}
  W::V
  b::V
  initial_params::F
  paramlength::Int
end
function Mapper(in::Integer)
  W = ones(Float32,in)
  b = zeros(Float32,in)
  p = vcat(W,b)
  initial_params() = p
  Mapper(W, b, initial_params, length(p))
end

function (m::Mapper{<:Vector{T}})(x::AbstractVecOrMat{T}, p) where T
  Wl = length(m.W)
  W = @view p[1 : Wl]
  b = @view p[Wl + 1 : end]
  W .* x .+ b
end
Base.show(io::IO, m::Mapper) = print(io, "Mapper(", length(m.W), ")")

initial_params(m::Mapper) = m.initial_params()
paramlength(m::Mapper) = m.paramlength

import Base: +, -
struct ZeroSynapse <: CurrentComponent end
ZeroSynapse(args...) = ZeroSynapse()
(m::ZeroSynapse)(h,x,p=nothing) = zeros(eltype(x),size(x,1))#fill(Zeros(), size(x,1))##Flux.Zeros()#zeros(I,size(h,2))
+(::Zeros, b::Real) = b
+(a::Real, ::Zeros) = a
-(::Zeros, b::Real) = -b
-(a::Real, ::Zeros) = -a
initial_params(m::ZeroSynapse) = Float32[]
paramlength(m::ZeroSynapse) = 0

struct SigmoidSynapse{V,F<:Function} <: CurrentComponent
  μ::V
  σ::V
  G::V
  E::V
  initial_params::F
end
function SigmoidSynapse(;
             μr=Float32.((0.3, 0.8)), σr=Float32.((3, 8)), Gr=Float32.((0.001, 1)), Er=Float32.((-0.3, 0.3)))
  μ = rand_uniform(eltype(μr),μr..., 1, 1)[1]#[1,:]
  σ = rand_uniform(eltype(σr),σr..., 1, 1)[1]#[1,:]
  G = rand_uniform(eltype(Gr),Gr..., 1, 1)[1]#[1,:]
  E = rand_uniform(eltype(Er),Er..., 1, 1)[1]#[1,:]
  p = vcat(μ,σ,G,E)
  initial_params() = p
  SigmoidSynapse(μ,σ,G,E,initial_params)
end
function (m::SigmoidSynapse{V})(h::AbstractVecOrMat, x::AbstractVecOrMat{T}, p) where {V,T}
  μ, σ, G, E = p
  @. G * sigmoid((x - μ) * σ) * (h - E)
end

Base.show(io::IO, m::SigmoidSynapse) = print(io, "SigmoidSynapse")
initial_params(m::SigmoidSynapse) = m.initial_params()
paramlength(m::SigmoidSynapse) = 4

struct LeakChannel{V,F<:Function} <: CurrentComponent
  G::V
  E::V
  initial_params::F
  paramlength::Int
end
function LeakChannel(n::Int; Gr=Float32.((0.001,1)), Er=Float32.((-0.3,0.3)))
  G = rand_uniform(eltype(Gr),Gr..., n)
  E = rand_uniform(eltype(Er),Er..., n)
  p = vcat(G,E)
  initial_params() = p
  LeakChannel(G, E, initial_params, length(p))
end
function (m::LeakChannel{V})(h::AbstractVecOrMat{T}, p) where {V,T}
  Gl = Int(size(p,1) / 2)
  G = @view p[1:Gl]
  E = @view p[Gl+1:end]
  @. G * (h - E)
end
Base.show(io::IO, m::LeakChannel) = print(io, "LeakChannel(", size(m.G,1), ")")
initial_params(m::LeakChannel) = m.initial_params()
paramlength(m::LeakChannel) = m.paramlength

struct ComponentContainer{C,F<:Function}
  components::C
  initial_params::F
  paramlength::Int
end
function ComponentContainer(components)
  p = reduce(vcat, DiffEqFlux.initial_params.(components))
  initial_params() = p
  ComponentContainer(components, initial_params, length(p))
end
# (m::ComponentContainer{C})(x::AbstractVecOrMat, I::AbstractVecOrMat) where C = doComponentContainer(m.components, x, I)
# doComponentContainer(components, x, I) = reshape(mapreduce(dst -> mapreduce(src -> components[src,dst](x[dst,:],I[src,:]), +, 1:size(I,1)), vcat, 1:size(x,1)), size(x))
# Flux.@functor ComponentContainer
function (m::ComponentContainer{C})(x::AbstractVecOrMat, I::AbstractVecOrMat, p) where C
  out = Zygote.Buffer(x, size(x,1), size(x,2))
  for i in eachindex(out)
    out[i] = 0
  end
  comps = m.components
  i = 1
  @inbounds for dst in 1:size(x,1)

    buf = Zygote.Buffer(x, size(comps,1), size(x,2))
    for j in eachindex(buf)
      buf[j] = 0
    end

    for src in 1:size(comps,1)
      c = comps[src,dst]
      cpl = paramlength(c)
      p_comp = @view p[i:i-1+cpl]
      c_out = c(x[dst,:],I[src,:],p_comp)
      buf[src,:] = c_out
      i = i + cpl
    end

    bb = copy(buf)
    out[dst,:] = sum(bb,dims=1)
  end
  return copy(out)
end

initial_params(m::ComponentContainer) = m.initial_params()
paramlength(m::ComponentContainer) = m.paramlength

struct LTCCell{W<:Wiring,SE<:ComponentContainer,SY<:ComponentContainer,LE<:LeakChannel,V<:AbstractArray,S<:AbstractMatrix,SOLVER,SENSEALG,F<:Function}
  wiring::W
  sens::SE
  syns::SY
  leaks::LE
  # cc::CC
  # cc_p::CCP
  # cc_re::CCRE
  cm::V
  state0::S
  solver::SOLVER
  sensealg::SENSEALG
  initial_params::F
  paramlength::Int
end

function LTCCell(wiring, solver, sensealg; cmr=Float32.((1.6,2.5)), state0r=Float32.((0.01)))
  n_in = wiring.n_in
  out = wiring.out
  n_total = wiring.n_total
  cm    = rand_uniform(eltype(cmr),cmr..., n_total)
  state0 = fill(state0r[1], n_total, 1)

  sens = Union{ZeroSynapse,SigmoidSynapse}[]
  syns = Union{ZeroSynapse,SigmoidSynapse}[]
  # sens = SigmoidSynapse[]
  # syns = SigmoidSynapse[]
  # sens = []
  # syns = []

  for dst in 1:size(wiring.sens_mask,2)
    for src in 1:size(wiring.sens_mask,1)
      toadd = wiring.sens_mask[src,dst]
      s = toadd == 0 ? ZeroSynapse() : SigmoidSynapse()
      push!(sens, s)
    end
  end
  for dst in 1:size(wiring.syn_mask,2)
    for src in 1:size(wiring.syn_mask,1)
      toadd = wiring.syn_mask[src,dst]
      s = toadd == 0 ? ZeroSynapse() : SigmoidSynapse()
      push!(syns, s)
    end
  end

  sens = reshape(sens, size(wiring.sens_mask))
  syns = reshape(syns, size(wiring.syn_mask))

  #sens = StructArray(sens)
  #syns = StructArray(syns)

  sensc = ComponentContainer(sens)
  synsc = ComponentContainer(syns)
  leaks = LeakChannel(n_total)

  # cc = Flux.Parallel(+, sensc, synsc, leaks)
  # cc_p, cc_re = Flux.destructure(cc)

  #p = vcat(DiffEqFlux.initial_params(sensc),DiffEqFlux.initial_params(synsc),DiffEqFlux.initial_params(leaks),cm,vec(state0))
  p = vcat(DiffEqFlux.initial_params(sensc),DiffEqFlux.initial_params(synsc),DiffEqFlux.initial_params(leaks),cm)
  initial_params() = p

  LTCCell(wiring, sensc, synsc, leaks, cm, state0, solver, sensealg, initial_params, length(p))
end
# Flux.@functor LTCCell
# Flux.trainable(m::LTCCell) = (m.cc_p, m.cm, m.state0,)
Base.show(io::IO, m::LTCCell) = print(io, "LTCCell(", m.wiring.n_sensory, ",", m.wiring.n_inter, ",", m.wiring.n_command, ",", m.wiring.n_motor, ")")
initial_params(m::LTCCell) = m.initial_params()
paramlength(m::LTCCell) = m.paramlength

#function (m::LTCCell{W,SE,SY,LE,CC,CCP,CCRE,V,<:AbstractMatrix{T},SOLVER,SENSEALG})(h::AbstractVecOrMat{T}, x::AbstractVecOrMat{T}, p) where {W,SE,SY,LE,CC,CCP,CCRE,V,T,SOLVER,SENSEALG}
function (m::LTCCell)(h::AbstractVecOrMat, x::AbstractVecOrMat, p)
  h = repeat(h, 1, size(x,2)-size(h,2)+1)
  # p_ode = @view p[1:end-length(m.state0)]
  p_ode = p
  h = solve_ode(m,h,x, p_ode)
  h, h
end

function solve_ode(m,h,x,p)#::Matrix{Float32} where T

  function dltcdt!(dh,h,p,t)
    #cc_p = @view p[1 : cc_p_l]
    sens_p  = @view p[1 : sens_pl]
    syns_p  = @view p[sens_pl + 1 : sens_pl + syns_pl]
    leaks_p = @view p[sens_pl + syns_pl + 1 : sens_pl + syns_pl + leaks_pl]
    cm      = @view p[sens_pl + syns_pl + leaks_pl + 1 : end]
    #cc = cc_re(cc_p)
    #argsv = ((h,x), (h,h), (h,h))
    #I_components = mapreduce(i -> cc.layers[i](argsv[i][1],argsv[i][2]), cc.connection, 1:length(cc.layers))
    I_components = m.sens(h,x,sens_p) .+ m.syns(h,h,syns_p) .+ m.leaks(h,leaks_p)
    @. dh = - cm * I_components
    nothing
  end
  function oop(h,p,t)
    dh = similar(h)
    dltcdt!(dh,h,p,t)
    dh
  end

  sens_pl = paramlength(m.sens)
  syns_pl = paramlength(m.syns)
  leaks_pl = paramlength(m.leaks)

  prob = ODEProblem{true}(dltcdt!,h,Float32.((0,1)),p)
  # prob = ODEProblem{false}(oop,h,Float32.((0,1)),p)

  # de = ModelingToolkit.modelingtoolkitize(prob)
  # jac = eval(ModelingToolkit.generate_jacobian(de)[2])
  # f = ODEFunction((dh,h,p,t)->dltcdt!(dh,h,p,t,x,cc_p_l), jac=jac)
  # prob_jac = ODEProblem(f,h,tspan,p)

  solve(prob,m.solver; sensealg=m.sensealg, save_everystep=false, save_start=false, abstol=1e-3, reltol=1e-3)[:,:,end]
end

mutable struct LTCNet{MI<:Mapper,MO<:Mapper,T<:LTCCell,S,F<:Function}
  mapin::MI
  mapout::MO
  cell::T
  state::S
  initial_params::F
  paramlength::Int
  #LTCNet(mapin,mapout,cell,state) = new{typeof(mapin),typeof(mapout),typeof(cell),typeof(state)}(mapin,mapout,cell,state)
end
function LTCNet(wiring,solver,sensealg)
  mapin = Mapper(wiring.n_in)
  mapout = Mapper(wiring.n_total)
  cell = LTCCell(wiring,solver,sensealg)

  p = vcat(DiffEqFlux.initial_params(mapin), DiffEqFlux.initial_params(cell), DiffEqFlux.initial_params(mapout))
  initial_params() = p

  LTCNet(mapin,mapout,cell,cell.state0,initial_params,length(p))
end

function (m::LTCNet{MI,MO,T,<:AbstractMatrix{T2}})(x::AbstractVecOrMat{T2}, p) where {MI,MO,T,T2}
  mapin_pl = paramlength(m.mapin)
  cell_pl = paramlength(m.cell)
  mapout_pl = paramlength(m.mapout)
  p_mapin  = @view p[1 : mapin_pl]
  p_cell   = @view p[mapin_pl + 1 : mapin_pl + cell_pl]
  p_mapout = @view p[mapin_pl + cell_pl + 1 : mapin_pl + cell_pl + mapout_pl]

  x = m.mapin(x, p_mapin)
  m.state, y = m.cell(m.state, x, p_cell)
  m.mapout(y, p_mapout)
end

reset!(m::LTCNet) = (m.state = m.cell.state0)
reset_state!(m::LTCNet, p) = (m.state = reshape(p[end-length(m.cell.state0)+1:end],:,1))

initial_params(m::LTCNet) = m.initial_params()
paramlength(m::LTCNet) = m.paramlength

Base.show(io::IO, m::LTCNet) = print(io, "LTCNet(", m.mapin, ",", m.mapout, ",", m.cell, ")")

function get_bounds(m::LeakChannel)
end
function get_bounds(m::SigmoidSynapse)
  [0.1, 1, 0.001, -1] |> f32, [0.9, 10, 1, 1] |> f32
end
function get_bounds(m::Mapper)
  lb = [[-20.1 for i in 1:length(m.W)]...,
        [-20.1 for i in 1:length(m.b)]...] |> f32

  ub = [[20.1 for i in 1:length(m.W)]...,
        [20.1 for i in 1:length(m.b)]...] |> f32
  return lb, ub
end

function get_bounds(m::LTCNet)
  lower = [
      get_bounds(m.mapin)[1]...,
      get_bounds(m.mapout)[1]...,
      reduce(vcat,[get_bounds(s)[1] for s in [m.cell.sens.components;m.cell.syns.components]])...,
      [0.01 for _ in m.cell.state0]...,            #
      [-2 for _ in m.cell.state0]...,
      [0.001 for _ in m.cell.cm]...,
      [-50 for _ in m.cell.state0]...,
  ] |> f32
  upper = [
    get_bounds(m.mapin)[2]...,
    get_bounds(m.mapout)[2]...,
    reduce(vcat,[get_bounds(s)[1] for s in [m.cell.sens.components;m.cell.syns.components]])...,
    [10.0 for _ in m.cell.state0]...,
    [-0.1 for _ in m.cell.state0]...,
    [5 for _ in m.cell.cm]...,
    [20 for _ in m.cell.state0]...,
  ] |> f32

  lower, upper
end

ext: jl


text:
mutable struct Value
	value :: Int
	owner# :: Union{Missing, Node} # circular dependency on each other so we'll have to leave it as Any
end

mutable struct Node
	left :: Union{Value, Node}
	right :: Union{Value, Node}
	parent :: Union{Missing, Nothing, Node}
end

function makenode(numbers :: Union{Vector, Int})
	if numbers isa Vector
		left, right = numbers
		return Node(makenode(left), makenode(right), missing)
	elseif numbers isa Int
		return Value(numbers, missing)
	end
end

function assignparents!(node :: Node)
	if node.left isa Node
		node.left.parent = node
		assignparents!(node.left)
	else
		node.left.owner = node
	end

	if node.right isa Node
		node.right.parent = node
		assignparents!(node.right)
	else
		node.right.owner = node
	end
end

function getorder!(current :: Node, result :: Vector{Value})
	for next in (current.left, current.right)
		if next isa Node
			getorder!(next, result)
		else
			push!(result, next)
		end
	end
end

function findexplodeable(values :: Vector{Value})
	eligibles = Value[]
	for value in values
		current = value.owner
		depth = 1
		while !isnothing(current.parent)
			current = current.parent
			depth += 1
		end

		if depth >= 5
			push!(eligibles, value)
		end
	end

	for index in 1:length(eligibles) - 1
		if eligibles[index].owner == eligibles[index + 1].owner
			return eligibles[index].owner
		end
	end
end

magnitude(n :: Node) = 3 * magnitude(n.left) + 2 * magnitude(n.right)
magnitude(v :: Value) = v.value

function reducetree(node1 :: Node, node2 :: Node)
	newtotal = Node(node1, node2, nothing)
	newtotal.left.parent = newtotal
	newtotal.right.parent = newtotal

	while true
		valuesinorder = Value[]
		getorder!(newtotal, valuesinorder)

		explodeable = findexplodeable(valuesinorder)

		splitindex = findfirst(v -> v.value >= 10, valuesinorder)

		if !isnothing(explodeable)
			left = explodeable.left
			right = explodeable.right

			leftindex = findfirst(==(left), valuesinorder)
			rightindex = findfirst(==(right), valuesinorder)
			if leftindex - 1 >= 1
				valuesinorder[leftindex - 1].value += left.value
			end

			if rightindex + 1 <= length(valuesinorder)
				valuesinorder[rightindex + 1].value += right.value
			end

			parent = explodeable.parent
			if parent.left == explodeable
				parent.left = Value(0, parent)
			else
				parent.right = Value(0, parent)
			end

		elseif !isnothing(splitindex)
			value = valuesinorder[splitindex]
			parent = value.owner

			lower = Value(Int(floor(value.value / 2)), missing)
			upper = Value(Int(ceil(value.value / 2)), missing)
			newnode = Node(lower, upper, parent)
			lower.owner = newnode
			upper.owner = newnode
			if parent.left == value
				parent.left = newnode
			else
				parent.right = newnode
			end

		else
			break # cannot be further reduced
		end
	end

	return newtotal
end

function part1()
	lines = readlines("inputD18.txt")
	snailnumbers = makenode.(eval.(Meta.parse.(lines)))

	assignparents!.(snailnumbers)

	total = first(snailnumbers)

	for index in 2:length(snailnumbers)
		total = reducetree(total, snailnumbers[index])
	end

	return magnitude(total)
end

function part2()
	lines = readlines("inputD18.txt")
	peakmag = 0
	for i in 1:100
		for j in 1:100
			if i != j
				println((i, j, peakmag))
				snailnumbers = makenode.(eval.(Meta.parse.(lines)))

				assignparents!.(snailnumbers)

				total = reducetree(snailnumbers[i], snailnumbers[j])

				mag = magnitude(total)
				if mag > peakmag
					peakmag = mag
				end
			end
		end
	end

	return peakmag
end

println(part1())
println(part2()) # warning: very slow


text:
pareto_dominates_clear(u::Vector{T}, v::Vector{T}) where {T <: Real} =
    all(u .<= v) && any(u .< v)

function pareto_dominates_fast(u::Vector{T}, v::Vector{T}) where {T <: Real}
    anylt = false
    @inbounds for i in eachindex(u)
        if u[i] > v[i]
            return false
        elseif u[i] < v[i]
            anylt = true
        end
    end
    return anylt
end

pareto_dominates(v1::Vector{T}, v2::Vector{T}) where {T <: Real} =
    pareto_dominates_fast(v1, v2)

function pareto_dominates_hat(u::Vector{T}, v::Vector{T}) where {T <: Real}
    anygt = anylt = false
    @inbounds for i in eachindex(u)
        if u[i] > v[i]
            anygt = true
        elseif u[i] < v[i]
            anylt = true
        end
    end
    if anygt
        return anylt ? 0 : 1
    else
        return anylt ? -1 : 0
    end
end

pareto_dominates(v1::Vector{T}, v2::Vector{T}) where {T <: Real} =
    pareto_dominates_fast(v1, v2)

# FIXME
pareto_dominates(f1::NewFitness, f2::NewFitness) =
    pareto_dominates(fitnessvalues(f1), fitnessvalues(f2))


text:
function add_input!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer,
    config::Config
)
    graph, agent_list = state
    # neighbors of neighbors
    input_candidates = Integer[]
    for neighbor in inneighbors(graph, agent_idx)
        append!(input_candidates, setdiff(inneighbors(graph, neighbor), inneighbors(graph, agent_idx)))
    end
    shuffle!(input_candidates)
    # Order neighbors by frequency of occurence in input_candidates descending
    input_queue = first.(sort(collect(countmap(input_candidates)), by=last, rev=true))
    # add edges
    if (length(input_queue) - config.network.new_follows) < 0
        new_input_count = length(input_queue)
    else
        new_input_count = config.network.new_follows
    end
    for _ in 1:new_input_count
        new_Neighbor = popfirst!(input_queue)
        add_edge!(graph, new_Neighbor, agent_idx)
        if (abs(agent_list[agent_idx].opinion - agent_list[new_Neighbor].opinion) < config.opinion_treshs.follow
            && indegree(graph, agent_idx) > indegree(graph, new_Neighbor))
            add_edge!(graph, agent_idx, new_Neighbor)
        end
    end
    return state
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using Test, DependenciesParser

@testset "Base Test" begin
    for pkg ∈ DependenciesParser.data
        @inferred installable(pkg)
        deps = installable(pkg)
        @test last(installable(DependenciesParser.data[1], direct = true)) ⊆ last(deps)
        first(deps) && break
    end
end


text:
using LinearAlgebra: LinearAlgebra, Symmetric, SymTridiagonal, I

using NQCBase: NQCBase, PeriodicCell, InfiniteCell
using RingPolymerArrays: RingPolymerArrays, NormalModeTransformation

struct RingPolymerParameters{T<:AbstractFloat}
    n_beads::Int
    ω_n::T
    ω_n²::T
    springs::Symmetric{T,Matrix{T}}
    normal_mode_springs::Vector{T}
    transformation::NormalModeTransformation{T}
    quantum_atoms::Vector{Int}
    function RingPolymerParameters{T}(
        n_beads::Integer, temperature::Real,
        quantum_atoms::Vector{<:Integer}) where {T<:AbstractFloat}

        ω_n = n_beads * temperature
        new(n_beads, ω_n, ω_n^2,
            get_spring_matrix(n_beads, ω_n),
            get_normal_mode_springs(n_beads, ω_n),
            NormalModeTransformation{T}(n_beads),
            quantum_atoms,
        )
    end
end

Base.broadcastable(beads::RingPolymerParameters) = Ref(beads)

"""Constructor for choosing specific elements to be quantum."""
function RingPolymerParameters{T}(n_beads::Integer, temperature::Real, atom_types::AbstractVector{Symbol}, quantum_nuclei::Vector{Symbol}) where {T}
    quantum_atoms = findall(in(quantum_nuclei), atom_types)
    RingPolymerParameters{T}(n_beads, temperature, quantum_atoms)
end

"""Constructor for the case where all nuclei are quantum."""
function RingPolymerParameters{T}(n_beads::Integer, temperature::Real, n_atoms::Integer) where {T}
    RingPolymerParameters{T}(n_beads, temperature, collect(1:n_atoms))
end

"""
    get_L(n_beads, mass, ω_n)

Return the Circulant symmetric matrix for the ring polymer springs.
"""
function get_spring_matrix(n_beads::Integer, ω_n::Real)::Symmetric
    if n_beads == 1
        spring = zeros(1, 1)
    elseif n_beads == 2
        spring = [2 -2; -2 2]
    else
        spring = SymTridiagonal(fill(2, n_beads), fill(-1, n_beads-1))
        spring = convert(Matrix, spring)
        spring[end,1] = spring[1, end] = -1
    end
    Symmetric(spring .*  ω_n^2 / 2)
end

get_normal_mode_springs(n_beads::Integer, ω_n::Real) = get_matsubara_frequencies(n_beads, ω_n) .^2 / 2
get_matsubara_frequencies(n::Integer, ω_n::Real) = 2ω_n*sin.((0:n-1)*π/n)

nbeads(beads::RingPolymerParameters) = beads.n_beads
Base.length(beads::RingPolymerParameters) = nbeads(beads)
Base.range(beads::RingPolymerParameters) = range(1; length=length(beads))

"""
    cayley_propagator(beads::RingPolymerParameters{T}, dt::Real; half::Bool=true) where {T}

J. Chem. Phys. 151, 124103 (2019); doi: 10.1063/1.5120282
"""
function cayley_propagator(beads::RingPolymerParameters{T}, dt::Real; half::Bool=true) where {T}

    cay(dtA::Matrix)::Matrix = LinearAlgebra.inv(I - dtA/2) * (I + dtA/2)

    ω_k = get_matsubara_frequencies(length(beads), beads.ω_n)
    prop = [Array{T}(undef, 2, 2) for i=1:length(beads)]
    for (i, ω) in enumerate(ω_k)
        A = [0 1; -ω^2 0]
        prop[i] .= half ? real.(sqrt(cay(dt.*A))) : cay(dt.*A)
    end
    prop
end

"""
    get_spring_energy(beads::RingPolymerParameters, masses, R)

Calculate the ring polymer spring potential.
"""
function get_spring_energy(beads::RingPolymerParameters, masses, R)
    E = zero(eltype(R))
    for bead=1:nbeads(beads)-1
        for i in beads.quantum_atoms # Only for quantum nuclei
            for j in axes(R, 1)
                E += masses[i] * (R[j, i, bead] - R[j, i, bead+1])^2
            end
        end
    end
    for i in beads.quantum_atoms
        for j in axes(R, 1)
            E += masses[i] * (R[j, i, end] - R[j, i, 1])^2
        end
    end
    E * beads.ω_n^2/2
end

function NQCBase.apply_cell_boundaries!(cell::PeriodicCell, R::AbstractArray{T,3}, beads::RingPolymerParameters) where {T}
    RingPolymerArrays.transform_to_normal_modes!(R, beads.transformation)
    R[:,beads.quantum_atoms,1] ./= sqrt(length(beads))
    @views NQCBase.apply_cell_boundaries!(cell, R[:,:,1])
    R[:,beads.quantum_atoms,1] .*= sqrt(length(beads))
    RingPolymerArrays.transform_from_normal_modes!(R, beads.transformation)
end
NQCBase.apply_cell_boundaries!(::InfiniteCell, ::AbstractArray{T,3}, ::RingPolymerParameters) where {T} = nothing


text:
export  asdf_rqml,
        asdf_waux,
        asdf_wqml,
        read_asdf_evt,
        read_hdf5,
        read_hdf5!,
        scan_hdf5,
        write_hdf5


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.7.1
#   kernelspec:
#     display_name: Julia 1.5.3
#     language: julia
#     name: julia-1.5
# ---

using Plots, Statistics
using Random
using Pkg
Pkg.activate(@__DIR__)
using Revise

# +
using ParticleInCell

alpha = 0.1
kx = 0.5
nx = 50
np = 10000 * nx
xmin, xmax = 0, 4π
mesh = OneDGrid( xmin, xmax, nx)
particles = ParticleGroup{1,1}(
        np;
        charge = 1.0,
        mass = 1.0,
        n_weights = 1)

sampler = LandauDamping( alpha, kx )
sample!( particles, mesh, sampler)
# -

p = plot(layout=2)
histogram!(p[1], particles.array[1,:], normalize=true, bins=50)
histogram!(p[2], particles.array[2,:], normalize=true, bins=50)

# +
function compute_rho!( rho, pm, particles)

    fill!(rho, 0.0)

    for i_part = 1:particles.n_particles
       xi = particles.array[1, i_part]
       wi = particles.array[3, i_part]
       add_charge!(rho, pm, xi, wi)
    end

end
# -

rho = zeros(Float64, nx)
compute_rho!(rho, pm, particles)

# +
maxwell = Maxwell1DFEM(mesh, spline_degree);

ex = zeros(Float64, nx)
x = LinRange(xmin, xmax, nx+1)[1:end-1]
compute_e_from_rho!( ex, maxwell, rho )
sval = eval_uniform_periodic_spline_curve(spline_degree-1, ex);
# -

p = plot(layout=2)
plot!(p[1], x, rho)
plot!(p[1], x, 1 .+ alpha * cos.(kx * x))
plot!(p[2], x, sval)
plot!(p[2], x, alpha/kx * sin.(kx * x))


text:
#
# images.jl -
#
# Manipulation of Tk images.
#

import ColorTypes, FixedPointNumbers

TkImage(obj::TkObject, kind::Name, args...; kwds...) =
    TkImage(getinterp(obj), kind, args; kwds...)

TkImage(interp::TclInterp, kind::Name; kwds...) =
    __wrapimage(interp, kind,
                Tcl.eval(interp, "image", "create", kind; kwds...))

function TkImage(interp::TclInterp, kind::AbstractString,
                 name::AbstractString; kwds...)
    if Tcl.eval(TclStatus, interp, "image", "type", name) == TCL_OK
        # Image already exists.
        if kind != getresult(interp)
            Tcl.error("image exists with a different type")
        end
        if length(kwds) > 0
            Tcl.eval(interp, name, "configure"; kwds...)
        end
    else
        name = Tcl.eval(interp, "image", "create", kind, name; kwds...)
    end
    return __wrapimage(interp, kind, name)
end

function __wrapimage(interp::TclInterp, kind::AbstractString,
                     name::AbstractString)
    T = Symbol(string(lowercase(kind[1]), lowercase(kind[2:end])))
    return TkImage{T}(interp, name)
end

(img::TkImage)(args...; kdws...) =
    Tcl.eval(img.interp, img, args...; kwds...)

getinterp(img::TkImage) = img.interp

getpath(img::TkImage) = img.path

# FIXME: extend __objptr and add a TclObj field to the TkImage structure and
# replace getpath(img) below by just img.
TclObj(img::T) where {T<:TkImage} = TclObj{T}(__newobj(getpath(img)))

delete(img::TkImage) =
    Tcl.eval(getinterp(img), "delete", getpath(img))

getwidth(img::TkImage) =
    Tcl.eval(Int, getinterp(img), "width", getpath(img))

getheight(img::TkImage) =
    Tcl.eval(Int, getinterp(img), "height", getpath(img))

exists(img::TkImage) =
    Tcl.eval(TclStatus, getinterp(img), "image", "inuse", getpath(img)) == TCL_OK

Base.resize!(img::TkImage{:Photo}, args...) =
    setphotosize!(getinterp(img), getpath(img), args...)

Base.size(img::TkImage{:Photo}) =
    getphotosize(getinterp(img), getpath(img))

function Base.size(img::TkImage{:Photo}, i::Integer)
    i ≥ 1 || throw(BoundsError("out of bounds dimension index"))
    return (i ≤ 2 ? size(img)[i] : 1)
end

Base.size(img::TkImage) =
    getwidth(img), getheight(img)

function Base.size(img::TkImage, i::Integer)
    i ≥ 1 || throw(BoundsError("out of bounds dimension index"))
    return (i == 1 ? getwidth(img) :
            i == 2 ? getheight(img) :
            1)
end

#------------------------------------------------------------------------------
# Apply a "color" map to an array of gray levels.

colorize(arr::Array{T,N}, lut::AbstractVector{C}; kwds...) where {T<:Real,N,C} =
    colorize!(Array{C,N}(undef, size(arr)), arr, lut; kwds...)

function colorize!(dst::Array{C,N}, arr::Array{T,N},
                   lut::AbstractVector{C};
                   cmin::Union{Real,Nothing} = nothing,
                   cmax::Union{Real,Nothing} = nothing) where {T<:Real,N,C}
    # Get the clipping values (not needed if number of color levels smaller
    # than 2).
    local _cmin::T, _cmax::T
    if length(lut) < 2
        @assert length(arr) ≥ 1
        _cmin = _cmax = arr[1]
    elseif cmin === nothing
        if cmax === nothing
            _cmin, _cmax = extrema(arr)
        else
            _cmin, _cmax = minimum(arr), cmax
        end
    else
        if cmax === nothing
            _cmin, _cmax = cmin, maximum(arr)
        else
            _cmin, _cmax = cmin, cmax
        end
    end
    colorize!(dst, arr, _cmin, _cmax, lut)
end

function threshold!(dst::AbstractArray{C,N},
                    src::AbstractArray{T,N},
                    lvl::T, lo::C, mid::C, hi::C) where {T<:Real,N,C}
    @assert size(dst) == size(src)
    @inbounds for i in eachindex(dst, src)
        val = src[i]
        dst[i] = (val < lvl ? lo :
                  val > lvl ? hi :
                  mid)
    end
    return dst
end

function colorize!(dst::Array{C,N},
                   arr::Array{T,N},
                   cmin::T, cmax::T,
                   lut::AbstractVector{C})  :: Array{C,N} where {T<:Real,N,C}
    @assert size(dst) == size(arr)
    n = length(lut)
    @assert length(arr) ≥ 1
    @assert 1 ≤ n

    if n < 2
        # Just fill the result with the same look-up table entry.
        fill!(dst, lut[1])
    elseif cmin == cmax
        # Perform a thresholding.
        threshold!(dst, src, cmin, lut[1], lut[div(n,2)+1], lut[n])
    else
        # Use at least Float32 for the computations.
        R = promote_type(Float32, T)
        kmin = 1
        kmax = n
        scl = R(kmax - kmin)/(R(cmax) - R(cmin))
        off = R(cmin*kmax - cmax*kmin)/R(kmax - kmin)
        # FIXME: adjust cmin and cmax to speedup
        # adj = (cmax - cmin)/(2*(kmax - kmin))
        # cmin += zadj
        # cmax -= zadj
        if scl > zero(R)
            @inbounds for i in eachindex(dst, src)
                val = src[i]
                k = (val ≤ cmin ? kmin :
                     val ≥ cmax ? kmax :
                     round(Int, (R(val) - off)*scl)) :: Int
                dst[i] = lut[k]
            end
        else
            @inbounds for i in eachindex(dst, src)
                val = src[i]
                k = (val ≥ cmin ? kmin :
                     val ≤ cmax ? kmax :
                     round(Int, (R(val) - off)*scl)) :: Int
                dst[i] = lut[k]
            end
        end
    end
    return dst
end

#------------------------------------------------------------------------------
# Implement reading/writing of Tk "photo" images.

mutable struct TkPhotoImageBlock
    # Pointer to the first pixel.
    ptr::Ptr{UInt8}

    # Width of block, in pixels.
    width::Cint

    # Height of block, in pixels.
    height::Cint

    # Address difference between corresponding pixels in successive lines.
    pitch::Cint

    # Address difference between successive pixels in the same line.
    pixelsize::Cint

    # Address differences between the red, green, blue and alpha components of
    # the pixel and the pixel as a whole.
    red::Cint
    green::Cint
    blue::Cint
    alpha::Cint

    TkPhotoImageBlock() = new(C_NULL,0,0,0,0,0,0,0,0)
end

findphoto(img::TkImage) = findphoto(getinterp(img), getpath(img))

function findphoto(interp::TclInterp, name::AbstractString)
    imgptr = ccall((:Tk_FindPhoto, libtk), Ptr{Cvoid},
                   (TclInterpPtr, Ptr{UInt8}), interp.ptr, name)
    if imgptr == C_NULL
        Tcl.error("invalid image name")
    end
    return imgptr
end

getpixels(img::TkImage, args...) =
    getpixels(getinterp(img), getpath(img), args...)

getpixels(name::Name, args...) =
    getpixels(getinterp(), name, args...)

getpixels(interp::TclInterp, name::Symbol, args...) =
    getpixels(interp, string(name), args...)

function getpixels(interp::TclInterp, name::AbstractString,
                   colormode::Type{T} = Val{:gray}) where {T <: Val}
    # Get photo image data.
    imgptr = findphoto(interp, name)
    block = TkPhotoImageBlock()
    code = ccall((:Tk_PhotoGetImage, libtk), Cint,
                 (Ptr{Cvoid}, Ref{TkPhotoImageBlock}),
                 imgptr, block)
    if code != 1
        error("unexpected returned code")
    end
    width     = Int(block.width)
    height    = Int(block.height)
    pixelsize = Int(block.pixelsize)
    pitch     = Int(block.pitch)
    @assert pitch ≥ width*pixelsize
    @assert rem(pitch, pixelsize) == 0
    src = unsafe_wrap(Array, block.ptr,
                      (pixelsize, div(pitch, pixelsize), height), false)
    r = 1 + Int(block.red)
    g = 1 + Int(block.green)
    b = 1 + Int(block.blue)
    a = 1 + Int(block.alpha)

    return __getpixels(T, src, width, height, r, g, b, a)
end

function __getpixels(::Type{Val{:gray}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    # Below is an approximation to:
    #	GRAY = 0.30*RED + 0.59*GREEN + 0.11*BLUE
    # rounded to nearest integer.
    cr = UInt16( 77)
    cg = UInt16(151)
    cb = UInt16( 28)
    c0 = UInt16(128)
    dst = Array{UInt8,2}(undef, width, height)
    @inbounds for y in 1:height
        @simd for x in 1:width
            dst[x,y] = ((cr*src[r,x,y] + cg*src[g,x,y] +
                         cb*src[b,x,y] + c0) >> 8)
        end
    end
    return dst
end

function __getpixels(::Type{Val{:red}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    return __getpixels(src, width, height, r)
end

function __getpixels(::Type{Val{:green}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    return __getpixels(src, width, height, g)
end

function __getpixels(::Type{Val{:blue}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    return __getpixels(src, width, height, b)
end

function __getpixels(::Type{Val{:alpha}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    return __getpixels(src, width, height, a)
end

function __getpixels(src::Array{UInt8,3}, width::Int, height::Int, c::Int)
    dst = Array{UInt8,3}(undef, width, height)
    @inbounds for y in 1:height
        @simd for x in 1:width
            dst[x,y] = src[c,x,y]
        end
    end
    return dst
end

function __getpixels(::Type{Val{:rgb}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    dst = Array{UInt8,3}(undef, 3, width, height)
    @inbounds for y in 1:height
        @simd for x in 1:width
            dst[1,x,y] = src[r,x,y]
            dst[2,x,y] = src[g,x,y]
            dst[3,x,y] = src[b,x,y]
        end
    end
    return dst
end

function __getpixels(::Type{Val{:rgba}}, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    dst = Array{UInt8,3}(undef, 4, width, height)
    @inbounds for y in 1:height
        @simd for x in 1:width
            dst[1,x,y] = src[r,x,y]
            dst[2,x,y] = src[g,x,y]
            dst[3,x,y] = src[b,x,y]
            dst[4,x,y] = src[a,x,y]
        end
    end
    return dst
end

function __getpixels(::Type, src::Array{UInt8,3},
                     width::Int, height::Int, r::Int, g::Int, b::Int, a::Int)
    error("invalid color mode")
end

getphotosize(img::TkImage) =
    getphotosize(getinterp(img), getpath(img))

function getphotosize(interp::TclInterp, name::AbstractString)
    w, h = __getphotosize(findphoto(interp, name))
    return (Int(w), Int(h))
end

function setphotosize!(interp::TclInterp, name::AbstractString,
                       width::Integer, height::Integer)
    __setphotosize(interp, findphoto(interp, name), Cint(width), Cint(height))
end

function __getphotosize(imgptr::Ptr{Cvoid})
    width, height = Ref{Cint}(0), Ref{Cint}(0)
    if imgptr != C_NULL
        ccall((:Tk_PhotoGetSize, libtk), Cvoid,
              (Ptr{Cvoid}, Ref{Cint}, Ref{Cint}),
              imgptr, width, height)
    end
    return (width[], height[])
end

function __setphotosize(interp::TclInterp, imgptr::Ptr{Cvoid},
                        width::Cint, height::Cint)
    code = ccall((:Tk_PhotoSetSize, libtk), Cint,
                 (TclInterpPtr, Ptr{Cvoid}, Cint, Cint),
                 interp.ptr, imgptr, width, height)
    code == TCL_OK || Tcl.error(tclresult(interp))
    return nothing
end

function __expandphotosize(interp::TclInterp, imgptr::Ptr{Cvoid},
                           width::Cint, height::Cint)
    code = ccall((:Tk_PhotoExpand, libtk), Cint,
                 (TclInterpPtr, Ptr{Cvoid}, Cint, Cint),
                 interp.ptr, imgptr, width, height)
    code == TCL_OK || Tcl.error(tclresult(interp))
    return nothing
end

function __setpixels(interp::TclInterp, name::AbstractString,
                     block::TkPhotoImageBlock,
                     x::Cint = Cint(0), y::Cint = Cint(0),
                     composite::Cint = TK_PHOTO_COMPOSITE_SET)
    # Get photo image.
    imgptr = findphoto(interp, name)
    width, height = __getphotosize(imgptr)
    println("image size: $width x $height")

    # Resize the image if it is too small.
    if width < x + block.width || height < y + block.height
        # Not clear (from Tcl/Tk doc.) why the following should be done and I
        # had to dive into the source code TkImgPhoto.c to figure out how to
        # actually resize the image (just calling Tk_PhotoSetSize with the
        # correct size yields segmentation fault).
        width = max(width, x + block.width)
        height = max(height, y + block.height)
        __setphotosize(interp, imgptr, zero(Cint), zero(Cint))
        __expandphotosize(interp, imgptr, width, height)
    end

    width, height = __getphotosize(imgptr)
    println("image size: $width x $height")

    # Assume (TCL_MAJOR_VERSION == 8) && (TCL_MINOR_VERSION >= 5), for older
    # versions, the interpreter argument is missing in Tk_PhotoPutBlock.
    code = ccall((:Tk_PhotoPutBlock, libtk), Cint,
                 (TclInterpPtr, Ptr{Cvoid}, Ref{TkPhotoImageBlock},
                  Cint, Cint, Cint, Cint, Cint),
                 interp.ptr, imgptr, block, x, y,
                 block.width, block.height, composite)
    code == TCL_OK || Tcl.error(tclresult(interp))

    if false
        # Notify that image has changed (FIXME: not needed as it is done by
        # Tk_PhotoPutBlock).
        ccall((:Tk_ImageChanged, libtk), Cvoid,
              (Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cint, Cint),
              imgptr, x, y, block.width, block.height,
              width, height)
    end

    return nothing
end

setpixels!(img::TkImage{:Photo}, args...) =
    setpixels!(getinterp(img), getpath(img), args...)

setpixels!(name::Name, args...) = setpixels!(getinterp(), name, args...)

setpixels!(interp::TclInterp, name::Symbol, args...) =
    setpixels!(interp, string(name), args...)

function setpixels!(interp::TclInterp, name::AbstractString,
                    src::DenseArray{UInt8,3})
    block = TkPhotoImageBlock()
    block.ptr       = pointer(src)
    block.pixelsize = size(src, 1)
    block.width     = size(src, 2)
    block.height    = size(src, 3)
    block.pitch     = block.pixelsize*block.width
    if block.pixelsize == 3
        block.red   = 0;
        block.green = 1;
        block.blue  = 2;
        block.alpha = 0;
    elseif block.pixelsize == 4
        block.red   = 0;
        block.green = 1;
        block.blue  = 2;
        block.alpha = 3;
    else
        error("invalid first dimension")
    end
    __setpixels(interp, name, block)
end

const Normed8 = FixedPointNumbers.Normed{UInt8,8}
const Gray8   = Union{UInt8,TkGray{UInt8},ColorTypes.Gray{Normed8}}
const RGB24   = Union{TkRGB{UInt8},ColorTypes.RGB{Normed8}}
const BGR24   = Union{TkBGR{UInt8},ColorTypes.BGR{Normed8}}
const RGBA32  = Union{TkRGBA{UInt8},ColorTypes.RGBA{Normed8}}
const BGRA32  = Union{TkBGRA{UInt8},ColorTypes.BGRA{Normed8}}
const ARGB32  = Union{TkARGB{UInt8},ColorTypes.ARGB{Normed8}}
const ABGR32  = Union{TkABGR{UInt8},ColorTypes.ABGR{Normed8}}

for (T, r, g, b, a) in ((:Gray8,  0, 0, 0, 0),
                        (:RGB24,  0, 1, 2, 0),
                        (:BGR24,  2, 1, 0, 0),
                        (:RGBA32, 0, 1, 2, 3),
                        (:BGRA32, 2, 1, 0, 3),
                        (:ARGB32, 1, 2, 3, 0),
                        (:ABGR32, 3, 2, 1, 0))
    @eval function setpixels!(interp::TclInterp, name::AbstractString,
                              A::DenseArray{T,2}) where {T<:$T}
        block = TkPhotoImageBlock()
        block.ptr       = pointer(A)
        block.pixelsize = sizeof(T)
        block.width     = size(A, 1)
        block.height    = size(A, 2)
        block.pitch     = block.pixelsize*block.width
        block.red       = $r
        block.green     = $g
        block.blue      = $b
        block.alpha     = $a
        __setpixels(interp, name, block)
    end
end


text:
using Base.Test
using ThermodynamicsTable
include("maintest.jl")
include("libtest.jl")


text:
using Test
using MappedArrays

# check for ambiguities
refambs = detect_ambiguities(Base, Core)
using MLLabelUtils
ambs = detect_ambiguities(MLLabelUtils, Base, Core)
#@test length(setdiff(ambs, refambs)) <= 24 # these are fine for now

tests = [
    "tst_labelencoding.jl"
    "tst_classify.jl"
    "tst_convertlabel.jl"
    "tst_labelmap.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
export TabularRandomPolicy

using StatsBase:Weights, sample

"""
    TabularRandomPolicy(prob::Array{Float64, 2})

`prob` describes the distribution of actions for each state.
"""
struct TabularRandomPolicy <: AbstractPolicy
    prob::Array{Float64,2}
end

(π::TabularRandomPolicy)(s::Int) = sample(Weights(get_prob(π, s)))
(π::TabularRandomPolicy)(obs) = π(get_state(obs))

RLBase.get_prob(π::TabularRandomPolicy, s) = @view π.prob[:, s]
RLBase.get_prob(π::TabularRandomPolicy, s, a) = π.prob[a, s]


text:
using SimplexSplitting
using Base.Test

# Only run these if the refinement tests fail
#include("test_tensordecomposition.jl")
#include("test_embed.jl")
#include("test_triangulate.jl")

include("triangulate.jl")
#@show example_triangulation(n_simplices = 100)

# If these succeed, then the above doesn't need testing
#include("test_refine_triangulation.jl")
#include("test_refine_triangulation_with_images.jl")
#include("test_refine_recursive.jl")
#include("test_refine_recursive_with_images.jl")

#include("test_refine_variable_k.jl")

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
module Bot

using Revise

using Discord

function bot()
    c = Client("BOT_KEY"; prefix = ">")

    function reply_bot_ping(c, m::Message)
        bot_ping = heartbeat_ping(c)
        reply(c, m, "My ping is $bot_ping.")
    end
    function reply_bot_info(c, m::Message)
        bot_guilds = get_guild(c)
        bot_members =
        reply(c, m, "\`\`\` Bot Version: 1.0.0
        \n Using: Discord.jl
        \n Current Guilds = $bot_guilds
        \n Current Members = $bot_members
        \`\`\`")
    end

    add_handler!(c, compile=true, MessageCreate, (_, e) -> println("Received message: $(e.message.content)"))
    add_command!(c, compile=true, :echo; pattern=r"^echo (.+)", help="Echos what you typed in") do c, msg, noprefix
        reply(c, msg, noprefix)
    end
    add_command!(c, compile=true, :roll; pattern=r"^roll (.+)", help="Rolls a dice with the number of sides given") do c, msg, value
        println("Number recieved: $value")
        isnum = isa(value, Number)
        println("Is num: $isnum")
        roll_num = int(rand(int32(0):int32(value)))
        reply(c, msg, "You got $roll_num.")
        # if isa(value, Number)
        #     println("Is number")
        #     roll_num = int(rand(int32(0):int32(value)))
        #     reply(c, msg, "You got $roll_num.")
        # else
        #     println("Not number")
        #     reply(c, msg, "Please enter a number!")
        # end
    end
    add_command!(c, compile=true, :botinfo, reply_bot_info; help = "Shows the stats of the bot")
    add_command!(c, compile=true, :ping, reply_bot_ping; help="Shows the ping from the bot to Discord servers")
    add_help!(c; pattern=r"^help(?: (.+))?", help="Show this help message", nohelp="No help provided", nocmd="Command not found")

    open(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Bot.bot()
    wait(c)
end


text:
module ManeuversWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractConstantThrustPropulsionModel = @jimport org.orekit.forces.maneuvers.propulsion.AbstractConstantThrustPropulsionModel
const AbstractDetector = @jimport org.orekit.propagation.events.AbstractDetector
const AbstractForceModel = @jimport org.orekit.forces.AbstractForceModel
const Action = @jimport org.hipparchus.ode.events.Action
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const Class = @jimport java.lang.Class
const ConfigurableLowThrustManeuver = @jimport org.orekit.forces.maneuvers.ConfigurableLowThrustManeuver
const ConstantThrustManeuver = @jimport org.orekit.forces.maneuvers.ConstantThrustManeuver
const DateBasedManeuverTriggers = @jimport org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const EventHandler = @jimport org.orekit.propagation.events.handlers.EventHandler
const Field = @jimport org.hipparchus.Field
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeDerivativesEquations = @jimport org.orekit.propagation.numerical.FieldTimeDerivativesEquations
const FieldVector3D = @jimport org.hipparchus.geometry.euclidean.threed.FieldVector3D
const ForceModel = @jimport org.orekit.forces.ForceModel
const Frame = @jimport org.orekit.frames.Frame
const ImpulseManeuver = @jimport org.orekit.forces.maneuvers.ImpulseManeuver
const JString = @jimport java.lang.String
const Maneuver = @jimport org.orekit.forces.maneuvers.Maneuver
const ManeuverTriggers = @jimport org.orekit.forces.maneuvers.trigger.ManeuverTriggers
const Object = @jimport java.lang.Object
const Orbit = @jimport org.orekit.orbits.Orbit
const ParameterDriver = @jimport org.orekit.utils.ParameterDriver
const PositionAngle = @jimport org.orekit.orbits.PositionAngle
const PropulsionModel = @jimport org.orekit.forces.maneuvers.propulsion.PropulsionModel
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SmallManeuverAnalyticalModel = @jimport org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const ThrustDirectionAndAttitudeProvider = @jimport org.orekit.forces.maneuvers.propulsion.ThrustDirectionAndAttitudeProvider
const TimeDerivativesEquations = @jimport org.orekit.propagation.numerical.TimeDerivativesEquations
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export ConfigurableLowThrustManeuver
export ConstantThrustManeuver
export ImpulseManeuver
export Maneuver
export SmallManeuverAnalyticalModel
export acceleration
export add_contribution
export apply
export depends_on_position_only
export equals
export event_occurred
export g
export get_attitude_override
export get_class
export get_date
export get_delta_v_sat
export get_direction
export get_duration
export get_end_date
export get_events_detectors
export get_field_events_detectors
export get_flow_rate
export get_handler
export get_inertial_dv
export get_inertial_frame
export get_isp
export get_jacobian
export get_maneuver_triggers
export get_max_check_interval
export get_max_iteration_count
export get_name
export get_parameter_driver
export get_parameters
export get_parameters_drivers
export get_propulsion_model
export get_start_date
export get_threshold
export get_thrust
export get_thrust_direction_provider
export get_thrust_vector
export get_trigger
export hash_code
export init
export is_firing
export is_forward
export is_supported
export notify
export notify_all
export reset_state
export to_string
export update_mass
export wait
export with_handler
export with_max_check
export with_max_iter
export with_threshold

include("PropulsionWrapper/PropulsionWrapper.jl")
include("TriggerWrapper/TriggerWrapper.jl")
include("configurable_low_thrust_maneuver.jl")
include("constant_thrust_maneuver.jl")
include("impulse_maneuver.jl")
include("maneuver.jl")
include("small_maneuver_analytical_model.jl")

end


text:
function ikegamida(m_prefs, f_prefs)

    m_num, f_num = size(m_prefs)[2], size(f_prefs)[2]
    m_pool, f_pool = collect(1:m_num), collect(1:f_num)
    match_or_unmatch = Array(Bool, m_num)
    for i in 1:m_num
        match_or_unmatch[i] = true
    end
    m_matched, f_matched = ones(Int64, m_num) + f_num, ones(Int64, f_num) + m_num

    m_rank = Array(Int64, (f_num + 1, m_num))
    f_rank = Array(Int64, (m_num + 1, f_num))

    sorting = Array(Int64, (f_num + 1, 2))
    sorting[:, 1] = collect(1:(f_num + 1))

    for i in 1:m_num
        sorting[:, 2] = m_prefs[:, i]
        m_rank[:,i] = sortrows(sorting, by = x->(x[2]))[:,1]
        sorting = sortrows(sorting, by = x->(x[1]))
    end

    sorting2 = Array(Int64, (m_num + 1, 2))
    sorting2[:, 1] = collect(1:(m_num + 1))

    for i in 1:f_num
        sorting2[:, 2] = f_prefs[:, i]
        f_rank[:,i] = sortrows(sorting2, by = x->(x[2]))[:,1]
        sorting2 = sortrows(sorting2, by = x->(x[1]))
    end

    while sum(match_or_unmatch) > 0
        for i in 1:m_num
            if match_or_unmatch[i] == true

                for j in m_prefs[:, i]

                    if j == 0
                        m_matched[i] = 0
                        match_or_unmatch[i] = false
                        break

                    else
                        no_marriage = f_rank[1, j]

                        if f_rank[i+1, j] < no_marriage

                            if f_matched[j] == m_num + 1
                                f_matched[j] = i
                                m_matched[i] = j
                                match_or_unmatch[i] = false
                                break

                            else

                                if f_rank[f_matched[j]+1, j] > f_rank[i+1, j]
                                    match_or_unmatch[f_matched[j]] = true
                                    f_matched[j] = i
                                    m_matched[i] = j
                                    match_or_unmatch[i] = false
                                    break
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    for (i,t) in enumerate(f_matched)
        if t == m_num + 1
            f_matched[i] = 0
        end
    end

     for (i,t) in enumerate(m_matched)
        if t == f_num + 1
            m_matched[i] = 0
        end
    end

    return m_matched, f_matched
end


text:
using DataArrays,DataFrames
using PyCall

@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.figure as pyfigure

function main()

  path = "../ml-100k/u.data"
  df = readrel(path,'\t')

  usergroups= sort(aggregate(df,:user_id,length), cols = :rating_length, rev=true)

  user_ids = usergroups[:user_id]

  indexes = 1:length(user_ids)

  num_ratings = usergroups[:rating_length]
  plt.figure(figsize=(40,10))
  plt.xlabel("USER_IDS")
  plt.ylabel("NUMBER_OF_RATINGS")
  plt.bar(indexes,num_ratings)
  plt.xticks(user_ids,rotation=45)

  plt.show()

end

# data =

#df[:user_id]

# describe(d)

# row = data[:]

# println(row)

# usergroups = groupby(x -> x[:x1],data)

# ndims(usergroups)

# for i in usergroups
#   println(length(i))
# end

# for i in usergroups
#   println(i[1])
# end

# digits = imap(i -> parseint(string(i[1][1])), groupedbyuser)

# println(digits)

# for i in usergroups

#   print("user_id: ")
#   print(i[1][1])
# #   print("  length: ")
# #   println(length(i))
# end

function readrel(relpath,delimiter=',')

  base = dirname(@__FILE__)

  if relpath[1] == '/'
    relpathwslash = relpath
  else
    relpathwslash = string("/",relpath)
  end

  path = string(base,relpathwslash)

  df = readtable(path,separator = delimiter)

  names!(df,[:user_id,:item_id,:rating,:timestamp])

  df

#   readdlm(path,delimiter)

end

# length([1,2,3])

# main()

# readrel("../ml-100k/u.data",'\t')

# str = "/foobar"

# str[1] == '/'

main()


text:
@testset "MDP Simulation" begin
    pomdp = SingleCovidPOMDP()

    b0 = initialstate(pomdp)
    s0 = rand(b0)
    a = rand(actions(pomdp))

    hist = simulate(pomdp, s0)
    validate_history(hist)
end

@testset "POMDP Simulation" begin
    testing_props = [0,0.5,1.0]
    pomdp = SingleCovidPOMDP(actions = CovidPOMDPs.CovidActionSpace(testing_props))
    A = actions(pomdp)
    @assert length(A) == length(testing_props)
    @assert all(isa.(A, CovidAction))

    sol = CovidPOMDPs.ProportionalControlSolver()
    planner = solve(sol, pomdp)

    b0 = initialstate(pomdp, Distributions.Uniform(1,10_000))
    s0 = rand(b0)
    sir = CovidPOMDPs.SIR(s0)
    @test all(≥(0), sir)
    @test sum(sir) == pomdp.N

    a = rand(actions(pomdp))

    hist = simulate(
        pomdp,
        b0,
        planner;
        T = 50,
        upd = BootstrapFilter,
        n_p = 1_000,
        progress = true
    )
    validate_history(hist)
end


text:
# Example showing the use of waveform triggering and segmented capture
using AcqSynth, PlotlyJS, Statistics
const acq = AcqSynth

# As in the first example, we setup the board
boardnum = 0
acq.init_board(boardnum)

clock_ratio = 20
num_seg=50000
seg_len=32000
numblocks = cld(2*num_seg*seg_len, 1024^2)
num_seg = numblocks*1024^2÷(2*seg_len)

# Use a helper function to configure for ttl triggering, with EXT clock and CHAN 1
acq.configure_for_ttl_triggering(boardnum,0,1;num_seg=num_seg,seg_len=seg_len)

# Check the frequency, in single channel mode, this will be twice the EXT clock
acq.get_frequency(boardnum) |> Int

# Trigger measurement and download the data from the card
@time data = acq.read_seg_waveforms_ddc(boardnum,numblocks,clock_ratio,seg_len)

# Plot the stuff
t = (1:size(data,1))./2E8 # time vector, we are sampling at 2GSPS and decimate by 10x
layout = Layout(yaxis_range=[0, 0.02])
plot(scatter(x=t,y=abs.(vec(mean(data,dims=2)))),layout)

# Taking many segments and downconverting is slow, especially if we are just going to average them
# instead, average in hardware (up to 64x)!
acq.set_segmented_capture(boardnum,0,seg_len)
acq.set_averager(boardnum,64,seg_len)
@time data = acq.read_seg_waveforms_ddc(boardnum,cld(numblocks,64),clock_ratio,seg_len)

# Plot the stuff
t = (1:size(data,1))./2E8 # time vector, we are sampling at 2GSPS and decimate by 10x
layout = Layout(yaxis_range=[0, 0.02])
plot(scatter(x=t,y=abs.(vec(mean(data,dims=2)))),layout)


text:
################################################################################
# Control Input
################################################################################

function input_impulse!(joint::Rotational{T}, pbody::Node, cbody::Node,
    timestep::T, clear::Bool) where T

    orientation_offset = joint.orientation_offset
    τ = joint.input * timestep
    xa, qa = current_configuration(pbody.state)
    xb, qb = current_configuration(cbody.state)

    pbody.state.Jτ2 += vector_rotate(-τ, orientation_offset)
    cbody.state.Jτ2 += vector_rotate(τ, inv(qb) * qa * orientation_offset)
    clear && (joint.input = szeros(T,3))
    return
end

################################################################################
# Control Jacobian
################################################################################

function input_jacobian_control(relative::Symbol,
    joint::Rotational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    timestep) where T

    orientation_offset = joint.orientation_offset
    if relative == :parent
        BFa = szeros(T, 3, 3)
        Bτa = - rotation_matrix(orientation_offset)
        return [BFa; Bτa] * timestep
    elseif relative == :child
        BFb = szeros(T, 3, 3)
        Bτb = rotation_matrix(inv(qb) * qa * orientation_offset)
        return [BFb; Bτb] * timestep
    end
end

function input_jacobian_configuration(relative::Symbol,
    joint::Rotational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion) where T

    orientation_offset = joint.orientation_offset
    τ = joint.input

    if relative == :parent
        FaXa = szeros(T,3,3)
        FaQa = szeros(T,3,4)
        τaXa = szeros(T,3,3)
        τaQa = szeros(T,3,4)
        FbXa = szeros(T,3,3)
        FbQa = szeros(T,3,4)
        τbXa = szeros(T,3,3)
        τbQa = rotation_matrix(inv(qb)) * ∂rotation_matrix∂q(qa * orientation_offset, τ)#*LVᵀmat(qa)
        return FaXa, FaQa, τaXa, τaQa, FbXa, FbQa, τbXa, τbQa
    elseif relative == :child
        FaXb = szeros(T,3,3)
        FaQb = szeros(T,3,4)
        τaXb = szeros(T,3,3)
        τaQb = szeros(T,3,4)
        FbXb = szeros(T,3,3)
        FbQb = szeros(T,3,4)
        τbXb = szeros(T,3,3)
        τbQb = ∂rotation_matrix_inv∂q(qb, vector_rotate(τ, qa * orientation_offset))#*LVᵀmat(qb)
        return FaXb, FaQb, τaXb, τaQb, FbXb, FbQb, τbXb, τbQb
    end
end


text:
module TreatmentPanels

include("TreatmentPanel.jl")
include("show_and_plot.jl")

# Export main types
export TreatmentPanel
export BalancedPanel, UnbalancedPanel

# Export treatment description types
export UnitTreatmentType
export SingleUnitTreatment, MultiUnitTreatment
export TreatmentTimingType
export Staggered, Simultaneous
export TreatmentDurationType
export Continuous, Discontinuous

# Export utility functions
export treated_ids, treated_labels, first_treated_period_ids, first_treated_period_labels, length_T₀, length_T₁
export get_y₀₀, get_y₀₁, get_y₁₀, get_y₁₁, decompose_y

end


text:
using Clang

# Only wrap VEDA for now since VERA requires C++
INCLUDE_PATH = "/usr/local/ve/veda/include"

HEADERS = [
    "veda.h",
    "veda_types.h",
    "veda_enums.h",
    "veda_device.h",
    "veda_device_omp.h",
    "veda_ptr.h",
]
HEADERS = map(h->joinpath(INCLUDE_PATH, h), HEADERS)
NAME = "veda"

wc = init(; headers = HEADERS,
            output_file = joinpath(@__DIR__, "lib$(NAME).jl"),
            common_file = joinpath(@__DIR__, "lib$(NAME)_common.jl"),
            clang_includes = vcat(INCLUDE_PATH, CLANG_INCLUDE),
            clang_args = ["-I", INCLUDE_PATH],
            header_wrapped = (root, current)->  current == root,
            header_library = x->"lib$(NAME)",
            clang_diagnostics = true,
            )
run(wc)


text:
mutable struct SDEIntegrator{algType,IIP,uType,uEltype,tType,tdirType,P,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,rateType,solType,cacheType,F4,F5,OType,noiseType,EventErrorType,CallbackCacheType} <: AbstractSDEIntegrator{algType,IIP,uType,tType}
  f::F4
  g::F5
  noise::noiseType
  uprev::uType
  tprev::tType
  t::tType
  u::uType
  p::P
  dt::tType
  dtnew::tType
  dtpropose::tType
  dtcache::tType
  T::tType
  tdir::tdirType
  just_hit_tstop::Bool
  isout::Bool
  event_last_time::Int
  vector_event_last_time::Int
  last_event_error::EventErrorType
  accept_step::Bool
  last_stepfail::Bool
  force_stepfail::Bool
  dtchangeable::Bool
  u_modified::Bool
  saveiter::Int
  alg::algType
  sol::solType
  cache::cacheType
  callback_cache::CallbackCacheType
  sqdt::tType
  W::randType
  opts::OType
  iter::Int
  success_iter::Int
  eigen_est::eigenType
  EEst::tTypeNoUnits
  q::tTypeNoUnits
  qold::tTypeNoUnits
  q11::tTypeNoUnits
end


text:
struct Executor
  net_sums::Array{<:NetworkSum}
  feed_dict::Dict
end

NetworkSum() = NetworkSum([])

function Executor(networks::Vector{<:Vector{<:AbstractTensor}})
  nodes, node_dict = generate_einsum_expr(networks; optimize=true)
  net_sums = [NetworkSum([n]) for n in nodes]
  return Executor(net_sums, node_dict)
end

@non_differentiable Executor(networks::Vector{<:Vector{<:AbstractTensor}})

function Executor(trees::Vector{SubNetwork}; optimize=true)
  nodes, node_dict = generate_einsum_expr(trees; optimize=optimize)
  net_sums = [NetworkSum([n]) for n in nodes]
  return Executor(net_sums, node_dict)
end

@non_differentiable Executor(trees::Vector{SubNetwork})

function Executor(networks::Vector{<:Vector{<:AbstractTensor}}, cache::NetworkCache)
  node_dict = Dict()
  for (node, index) in cache.index_dict
    i, j = index
    node_dict[node] = networks[i][j]
  end
  return Executor(cache.net_sums, node_dict)
end

@non_differentiable Executor(
  networks::Vector{<:Vector{<:AbstractTensor}}, cache::NetworkCache
)

function Executor(trees::Vector{SubNetwork}, cache::NetworkCache)
  return Executor(get_leaves(trees), cache)
end

@non_differentiable Executor(trees::Vector{SubNetwork}, cache::NetworkCache)

function run(executor::Executor; kwargs...)
  return [run(net_sum, executor.feed_dict; kwargs...) for net_sum in executor.net_sums]
end

Base.length(executor::Executor) = Base.length(executor.net_sums)


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))

using Oceananigans
using Oceananigans.Units
using Oceananigans.TimeSteppers: time_step!
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity

using Profile
using StatProfilerHTML
using BenchmarkTools

grid = RegularRectilinearGrid(size=64, z=(-64, 0), topology=(Flat, Flat, Bounded))

closure = CATKEVerticalDiffusivity()

Qᵇ = 1e-8
Qᵘ = - 1e-4
Qᵛ = 0.0

u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵛ))
b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ))

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    tracers = (:b, :e),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = FPlane(f=1e-4),
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    closure = closure)

N² = 1e-5
bᵢ(x, y, z) = N² * z
set!(model, b = bᵢ)

function nsteps!(model, n)
    for _ = 1:n
        time_step!(model, 1e-9)
    end
    return nothing
end

@info "Running one time-step..."
nsteps!(model, 1) # warmup

@info "Benchmarking..."
@btime time_step!(model, 1e-9)

n = 100
@info "Profiling $n time-steps..."
@profile nsteps!(model, n)

statprofilehtml()


text:
# ---
# title: 738. Monotone Increasing Digits
# id: problem738
# author: Indigo
# date: 2020-12-16
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/monotone-increasing-digits/description/>
# hidden: true
# ---
#
# Given a non-negative integer `N`, find the largest number that is less than or
# equal to `N` with monotone increasing digits.
#
# (Recall that an integer has _monotone increasing digits_ if and only if each
# pair of adjacent digits `x` and `y` satisfy `x <= y`.)
#
# **Example 1:**
#
#
#
#     Input: N = 10
#     Output: 9
#
#
# **Example 2:**
#
#
#
#     Input: N = 1234
#     Output: 1234
#
#
# **Example 3:**
#
#
#
#     Input: N = 332
#     Output: 299
#
#
# **Note:** `N` is an integer in the range `[0, 10^9]`.
#
#
## @lc code=start
using LeetCode

function monotone_increasing_digits(N::Int)::Int
    some_ones = 111111111
    res = 0
    for _ in 1:9
        while res + some_ones > N
            some_ones ÷= 10
        end
        res += some_ones
    end
    return res
end

## @lc code=end


text:
# Julia wrapper for header: /usr/local/include/libavcodec/qsv.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_qsv_alloc_context

function av_qsv_alloc_context()
    ccall((:av_qsv_alloc_context,libavcodec),Ptr{AVQSVContext},())
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
# read_samples

"""

Read sample output files created by StanSample.jl and return the requested `output_format`.
The default output_format is :namedtuple. Optionally the list of parameter symbols can be returned.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `model`                     : SampleModel
```

### Optional arguments
```julia
* `output_format=:namedtuple` : Requested format for samples
* `include_internals=false`   : Include internal Stan paramenters
* `return_parameters=false`   : Return a tuple of (output_format, parameter_symbols)
* `chains=1:m.n_chains[1]`    : Chains to be included in output
* `start=1`                   : First sample to be included
* `kwargs...`                 : Capture all other keyword arguments
```

Currently supported formats are:

1. :array (3d array format - [samples, parameters, chains])
2. :namedtuple (DEFAULT: NamedTuple object, all chains appended)
3. :namedtuples (Vector{NamedTuple} object, individual chains)
4. :table (Tables object, all chains appended)
5. :tables (Vector{Tables} object, individual chains)
6. :dataframe (DataFrames.DataFrame object, all chains appended)
7. :dataframes (Vector{DataFrames.DataFrame} object, individual chains)
8. :particles (Dict{MonteCarloMeasurements.Particles})
9. :mcmcchains (MCMCChains.Chains object)

Basically chains can be returned as a NamedTuple, a StanTable, a DataFrame,
a Particles or an MCMCChains.Chains object.

For NamedTuple, StanTable and DataFrame all chains are appended or can be returned
as a Vector{...} for each chain.

By default all chains will be read in. With the optional keyword argument `chains`
a subset of chains can be included, e.g. `chains = [2, 4]`.

The optional keyword argument `start` specifies which initial (warm-up) samples
should be removed.

Notes:
1. Use of the Stan `thinning` option will interfere with the value of start.
2. Start is the first sample included, e.g. with 1000 warm-up samples, start
should be set to 1001.

The NamedTuple output-format will extract and combine parameter vectors, e.g.
if Stan's cmdstan returns `a.1, a.2, a.3` the NamedTuple will just contain `a`.

For Tables object you can use the `select_block()` function to create an object
that conforms to the Tables interface:
```
stantable = read+samples(m10.4s; output_format=:table)
atable = select_block(stantable, "a")
```

The glue code for option 9 is enabled by Requires.jl if MCMCChains is loaded,.
"""
function read_samples(model::SampleModel;
  output_format=:namedtuple,
  include_internals=false,
  return_parameters=false,
  chains=1:model.n_chains[1],
  start=1,
  kwargs...)

  (res, names) = read_csv_files(model::SampleModel, output_format;
    include_internals=include_internals, start=start, chains=chains,
    kwargs...
  )

  if return_parameters
    return( (res, names) )
  else
    return(res)
  end

end