text:
function initialize!(integrator,cache::SSPRK22ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK22ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator

  # u1 -> stored as u
  u = uprev + dt*integrator.fsalfirst
  k = f(u, p, t+dt)
  # u
  u = (uprev + u + dt*k) / 2

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK22Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK22Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,stage_limiter!,step_limiter! = cache

  # u1 -> stored as u
  @. u = uprev + dt*fsalfirst
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  # u
  @. u = (uprev + u + dt*k) / 2
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK33ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK33ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator

  # u1
  u = uprev + dt*integrator.fsalfirst
  k = f(u, p, t+dt)
  # u2
  u = (3*uprev + u + dt*k) / 4
  k = f(u,p,t+dt/2)
  # u
  u = (uprev + 2*u + 2*dt*k) / 3

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK33Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK33Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,stage_limiter!,step_limiter! = cache

  # u1
  @. u = uprev + dt*fsalfirst
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  # u2
  @. u = (3*uprev + u + dt*k) / 4
  stage_limiter!(u, f, t+dt/2)
  f(k,u,p,t+dt/2)
  # u
  @. u = (uprev + 2*u + 2*dt*k) / 3
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK53ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK53ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache

  # u1
  tmp = uprev + β10 * dt * integrator.fsalfirst
  k = f(tmp, p, t+c1*dt)
  # u2 -> stored as u
  u = tmp + β21 * dt * k
  k = f(u, p, t+c2*dt)
  # u3
  tmp = α30 * uprev + α32 * u + β32 * dt * k
  k = f(tmp, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α43 * tmp + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u
  u = α52 * u + α54 * tmp + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK53Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK53Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4 = cache.tab

  # u1
  @. tmp = uprev + β10 * dt * fsalfirst
  stage_limiter!(tmp, f, t+c1*dt)
  f( k,  tmp, p, t+c1*dt)
  # u2 -> stored as u
  @. u = tmp + β21 * dt * k
  stage_limiter!(u, f, t+c2*dt)
  f( k,  u, p, t+c2*dt)
  # u3
  @. tmp = α30 * uprev + α32 * u + β32 * dt * k
  stage_limiter!(tmp, f, t+c3*dt)
  f( k,  tmp, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α43 * tmp + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u
  @. u = α52 * u + α54 * tmp + β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK63ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK63ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5 = cache

  # u1 -> stored as u
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2
  u₂ = u + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u3
  tmp = u₂ + β32 * dt * k
  k = f(tmp, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α41 * u + α43 * tmp + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u
  u = α62 * u₂ + α65 * tmp + β65 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK63Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK63Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₂,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5 = cache.tab

  # u1 -> stored as u
  @. u = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2
  @. u₂ = u + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u3
  @. tmp = u₂ + β32 * dt * k
  stage_limiter!(tmp, f, t+c3*dt)
  f( k,  tmp, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α41 * u + α43 * tmp + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u
  @. u = α62 * u₂ + α65 * tmp + β65 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK73ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK73ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6 = cache

  # u1
  u₁ = uprev + β10 * dt * integrator.fsalfirst
  k = f(u₁,p,t+c1*dt)
  # u2
  tmp = u₁ + β21 * dt * k
  k = f(tmp, p, t+c2*dt)
  # u3 -> stored as u
  u = tmp + β32 * dt * k
  k = f(u, p, t+c3*dt)
  # u4
  tmp = α40 * uprev + α43 * u + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = α50 * uprev + α51 * u₁ + α54 * tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u6
  tmp = tmp + β65 * dt * k
  k = f(tmp, p, t+c6*dt)
  # u
  u = α73 * u + α76 * tmp + β76 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK73Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK73Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₁,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6 = cache.tab

  # u1
  @. u₁ = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u₁, f, t+c1*dt)
  f(k,u₁,p,t+c1*dt)
  # u2
  @. tmp = u₁ + β21 * dt * k
  stage_limiter!(tmp, f, t+c2*dt)
  f( k,  tmp, p, t+c2*dt)
  # u3 -> stored as u
  @. u = tmp + β32 * dt * k
  stage_limiter!(u, f, t+c3*dt)
  f( k,  u, p, t+c3*dt)
  # u4
  @. tmp = α40 * uprev + α43 * u + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = α50 * uprev + α51 * u₁ + α54 * tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u6
  @. tmp = tmp + β65 * dt * k
  stage_limiter!(tmp, f, t+c6*dt)
  f( k,  tmp, p, t+c6*dt)
  # u
  @. u = α73 * u + α76 * tmp + β76 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK83ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK83ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7 = cache

  # u1 -> save as u
  u = uprev + β10 * dt * integrator.fsalfirst
  k = f(u, p, t+c1*dt)
  # u2
  u₂ = u + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u3
  u₃ = u₂ + β32 * dt * k
  k = f(u₃,p,t+c3*dt)
  # u4
  tmp = u₃ + β43 * dt * k
  k = f(tmp, p, t+c4*dt)
  # u5
  tmp = α50 * uprev + α51 * u + α54 * tmp + β54 * dt * k
  k = f(tmp, p, t+c5*dt)
  # u6
  tmp = α61 * u + α65 * tmp + β65 * dt * k
  k = f(tmp, p, t+c6*dt)
  # u7
  tmp = α72 * u₂ + α73 * u₃ + α76 * tmp + β76 * dt * k
  k = f(tmp, p, t+c7*dt)
  # u
  u = tmp + β87 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK83Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK83Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,tmp,u₂,u₃,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7 = cache.tab

  # u1 -> save as u
  @. u = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u, f, t+c1*dt)
  f( k,  u, p, t+c1*dt)
  # u2
  @. u₂ = u + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u3
  @. u₃ = u₂ + β32 * dt * k
  stage_limiter!(u₃, f, t+c3*dt)
  f(k,u₃,p,t+c3*dt)
  # u4
  @. tmp = u₃ + β43 * dt * k
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u5
  @. tmp = α50 * uprev + α51 * u + α54 * tmp + β54 * dt * k
  stage_limiter!(tmp, f, t+c5*dt)
  f( k,  tmp, p, t+c5*dt)
  # u6
  @. tmp = α61 * u + α65 * tmp + β65 * dt * k
  stage_limiter!(tmp, f, t+c6*dt)
  f( k,  tmp, p, t+c6*dt)
  # u7
  @. tmp = α72 * u₂ + α73 * u₃ + α76 * tmp + β76 * dt * k
  stage_limiter!(tmp, f, t+c7*dt)
  f( k,  tmp, p, t+c7*dt)
  # u
  @. u = tmp + β87 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK432ConstantCache)
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK432ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_2 = dt / 2

  # u1
  u = uprev + dt_2*integrator.fsalfirst
  k = f(u, p, t+dt_2)
  # u2
  u = u + dt_2*k
  k = f(u, p, t+dt)
  u = u + dt_2*k
  if integrator.opts.adaptive
    utilde = (uprev + 2*u) / 3
  end
  # u3
  u = (2*uprev + u) / 3
  k = f(u, p, t+dt_2)
  # u
  u = u + dt_2*k

  integrator.fsallast = f(u, p, t+dt)
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm(@. (utilde-u)/(integrator.opts.abstol+max(integrator.opts.internalnorm(uprev),integrator.opts.internalnorm(u))*integrator.opts.reltol))
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK432Cache)
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator,cache::SSPRK432Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,utilde,atmp,stage_limiter!,step_limiter! = cache
  dt_2 = dt / 2

  # u1
  @. u = uprev + dt_2*fsalfirst
  stage_limiter!(u, f, t+dt_2)
  f( k,  u, p, t+dt_2)
  # u2
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
  #
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt+dt_2)
  if integrator.opts.adaptive
    @. utilde = (uprev + 2*u) / 3
  end
  # u3
  @. u = (2*uprev + u) / 3
  f( k,  u, p, t+dt_2)
  #
  @. u = u + dt_2*k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)

  if integrator.opts.adaptive
    @. atmp = (utilde-u)/(integrator.opts.abstol+max(integrator.opts.internalnorm(uprev),integrator.opts.internalnorm(u))*integrator.opts.reltol)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK932ConstantCache)
  integrator.kshortsize = 1
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
end

@muladd function perform_step!(integrator,cache::SSPRK932ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_6 = dt / 6
  dt_3 = dt / 3
  dt_2 = dt / 2

  # u1
  u = uprev + dt_6*integrator.fsalfirst
  k = f(u, p, t+dt_6)
  # u2
  u = u + dt_6*k
  k = f(u, p, t+dt_3)
  # u3
  u = u + dt_6*k
  k = f(u, p, t+dt_2)
  # u4
  u = u + dt_6*k
  k = f(u, p, t+2*dt_3)
  # u5
  u = u + dt_6*k
  k = f(u, p, t+5*dt_6)
  # u6
  u = u + dt_6*k
  if integrator.opts.adaptive
    k = f(u, p, t+dt)
    utilde = (uprev + 6*u + 6*dt*k) / 7
  end
  # u6*
  u = (3*uprev + dt_2*integrator.fsalfirst + 2*u) / 5
  k = f(u, p, t+dt_2)
  # u7*
  u = u + dt_6*k
  k = f(u, p, t+2*dt_3)
  # u8*
  u = u + dt_6*k
  k = f(u, p, t+5*dt_6)
  # u
  u = u + dt_6*k

  integrator.fsallast = f(u, p, t+dt)
  if integrator.opts.adaptive
    integrator.EEst = integrator.opts.internalnorm(@. (utilde-u)/(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK932Cache)
  integrator.kshortsize = 1
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator,cache::SSPRK932Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,fsalfirst,utilde,atmp,stage_limiter!,step_limiter! = cache
  dt_6 = dt / 6
  dt_3 = dt / 3
  dt_2 = dt / 2

  # u1
  @. u = uprev + dt_6*fsalfirst
  stage_limiter!(u, f, t+dt_6)
  f( k,  u, p, t+dt_6)
  # u2
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt_3)
  f( k,  u, p, t+dt_3)
  # u3
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt_2)
  f( k,  u, p, t+dt_2)
  # u4
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+2*dt_3)
  f( k,  u, p, t+2*dt_3)
  # u5
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+5*dt_6)
  f( k,  u, p, t+5*dt_6)
  # u6
  @. u = u + dt_6*k
  if integrator.opts.adaptive
    stage_limiter!(u, f, t+dt)
    f( k,  u, p, t+dt)
    @. utilde = (uprev + 6*u + 6*dt*k) / 7
    stage_limiter!(utilde, f, t+dt)
  end
  # u6*
  @. u = (3*uprev + dt_2*fsalfirst + 2*u) / 5
  stage_limiter!(u, f, t+dt_6)
  f( k,  u, p, t+dt_2)
  # u7*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+2*dt_3)
  f( k,  u, p, t+2*dt_3)
  # u8*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+5*dt_6)
  f( k,  u, p, t+5*dt_6)
  # u9*
  @. u = u + dt_6*k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)

  if integrator.opts.adaptive
    @. atmp = (utilde-u)/(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK54ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::SSPRK54ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4 = cache

  # u₁
  u₂ = uprev + β10 * dt * integrator.fsalfirst
  k = f(u₂,p,t+c1*dt)
  # u₂
  u₂ = α20 * uprev + α21 * u₂ + β21 * dt * k
  k = f(u₂,p,t+c2*dt)
  # u₃
  u₃ = α30 * uprev + α32 * u₂ + β32 * dt * k
  k₃ = f(u₃,p,t+c3*dt)
  # u₄ -> stored as tmp
  tmp = α40 * uprev + α43 * u₃ + β43 * dt * k₃
  k = f(tmp, p, t+c4*dt)
  # u
  u = α52 * u₂ + α53 * u₃ + β53 * dt * k₃ + α54 * tmp + β54 * dt * k

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK54Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK54Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,k₃,u₂,u₃,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  @unpack β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4 = cache.tab

  # u₁
  @. u₂ = uprev + β10 * dt * integrator.fsalfirst
  stage_limiter!(u₂, f, t+c1*dt)
  f(k,u₂,p,t+c1*dt)
  # u₂
  @. u₂ = α20 * uprev + α21 * u₂ + β21 * dt * k
  stage_limiter!(u₂, f, t+c2*dt)
  f(k,u₂,p,t+c2*dt)
  # u₃
  @. u₃ = α30 * uprev + α32 * u₂ + β32 * dt * k
  stage_limiter!(u₃, f, t+c3*dt)
  f(k₃,u₃,p,t+c3*dt)
  # u₄ -> stored as tmp
  @. tmp = α40 * uprev + α43 * u₃ + β43 * dt * k₃
  stage_limiter!(tmp, f, t+c4*dt)
  f( k,  tmp, p, t+c4*dt)
  # u
  @. u = α52 * u₂ + α53 * u₃ + β53 * dt * k₃ + α54 * tmp + β54 * dt * k
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f( k,  u, p, t+dt)
end

function initialize!(integrator,cache::SSPRK104ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev,integrator.p,integrator.t) # Pre-start fsal
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::SSPRK104ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  dt_6 = dt/6
  dt_3 = dt/3
  dt_2 = dt/2

  tmp = uprev + dt_6 * integrator.fsalfirst # u₁
  k = f(tmp, p, t+dt_6)
  tmp = tmp + dt_6 * k # u₂
  k = f(tmp, p, t+dt_3)
  tmp = tmp + dt_6 * k # u₃
  k = f(tmp, p, t+dt_2)
  u = tmp + dt_6 * k # u₄
  k₄ = f(u, p, t+2*dt_3)
  tmp = (3*uprev + 2*u + dt_3 * k₄) / 5 # u₅
  k = f(tmp, p, t+dt_3)
  tmp = tmp + dt_6 * k # u₆
  k = f(tmp, p, t+dt_2)
  tmp = tmp + dt_6 * k # u₇
  k = f(tmp, p, t+2*dt_3)
  tmp = tmp + dt_6 * k # u₈
  k = f(tmp, p, t+5*dt_6)
  tmp = tmp + dt_6 * k # u₉
  k = f(tmp, p, t+dt)
  u = (uprev + 9*(u + dt_6*k₄) + 15*(tmp + dt_6*k)) / 25

  integrator.fsallast = f(u, p, t+dt) # For interpolation, then FSAL'd
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator,cache::SSPRK104Cache)
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst,integrator.uprev,integrator.p,integrator.t) # FSAL for interpolation
end

@muladd function perform_step!(integrator,cache::SSPRK104Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack k,k₄,tmp,fsalfirst,stage_limiter!,step_limiter! = cache
  dt_6 = dt/6
  dt_3 = dt/3
  dt_2 = dt/2

  @. tmp = uprev + dt_6 * integrator.fsalfirst
  stage_limiter!(tmp, f, t+dt_6)
  f( k,  tmp, p, t+dt_6)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_3)
  f( k,  tmp, p, t+dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_2)
  f( k,  tmp, p, t+dt_2)
  @. u = tmp + dt_6 * k
  stage_limiter!(u, f, t+2*dt_3)
  f(k₄,u,p,t+2*dt_3)
  @. tmp = (3*uprev + 2*u + 2*dt_6 * k₄) / 5
  stage_limiter!(tmp, f, t+dt_3)
  f( k,  tmp, p, t+dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt_2)
  f( k,  tmp, p, t+dt_2)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+2*dt_3)
  f( k,  tmp, p, t+2*dt_3)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+5*dt_6)
  f( k,  tmp, p, t+5*dt_6)
  @. tmp = tmp + dt_6 * k
  stage_limiter!(tmp, f, t+dt)
  f( k,  tmp, p, t+dt)
  @. u = (uprev + 9*(u + dt_6*k₄) + 15*(tmp + dt_6*k)) / 25
  stage_limiter!(u, f, t+dt)
  step_limiter!(u, f, t+dt)
  f(k, u, p, t+dt)
end


text:
include("test_examples_2d_part1.jl")
include("test_examples_2d_part2.jl")
include("test_examples_2d_part3.jl")

ext: jl


text:
module SingleCrystal

using LinearAlgebra
using JSON

# https://web.archive.org/web/20080324193801/http://cst-www.nrl.navy.mil/lattice/spcgrp/
# coordinate system -> standard cartesian
# primitive vectors -> describe the primitive cell
# basis vectors -> describe positions of atoms in the primitive cell
# translation of basis vectors by primitive vectors generate a crystal

# https://wiki.fysik.dtu.dk/ase/ase/spacegroup/spacegroup.html

"""
    get_chemical_info()

Returns instances `chemical_symbols`, `atomic_numbers` and `masses` for the setup of `Cell`.
"""
function get_chemical_info()
    chemical_symbols = [
        # 0
        "X",
        # 1
        "X", "He",
        # 2
        "Li", "Be", "B", "C", "N", "O", "F", "Ne",
        # 3
        "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar",
        # 4
        "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
        "Ga", "Ge", "As", "Se", "Br", "Kr",
        # 5
        "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd",
        "In", "Sn", "Sb", "Te", "I", "Xe",
        # 6
        "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy",
        "Ho", "Er", "Tm", "Yb", "Lu",
        "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi",
        "Po", "At", "Rn",
        # 7
        "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk",
        "Cf", "Es", "Fm", "Md", "No", "Lr",
        "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn", "Nh", "Fl", "Mc",
        "Lv", "Ts", "Og"
    ]

    atomic_numbers = Dict()
    for (Z, symbol) in enumerate(chemical_symbols)
        atomic_numbers[symbol] = Z
    end

    atomic_masses_iupac2016 = [
        1.0,  # X
        1.008,  # H [1.00784, 1.00811]
        4.002602,  # He
        6.94,  # Li [6.938, 6.997]
        9.0121831,  # Be
        10.81,  # B [10.806, 10.821]
        12.011,  # C [12.0096, 12.0116]
        14.007,  # N [14.00643, 14.00728]
        15.999,  # O [15.99903, 15.99977]
        18.998403163,  # F
        20.1797,  # Ne
        22.98976928,  # Na
        24.305,  # Mg [24.304, 24.307]
        26.9815385,  # Al
        28.085,  # Si [28.084, 28.086]
        30.973761998,  # P
        32.06,  # S [32.059, 32.076]
        35.45,  # Cl [35.446, 35.457]
        39.948,  # Ar
        39.0983,  # K
        40.078,  # Ca
        44.955908,  # Sc
        47.867,  # Ti
        50.9415,  # V
        51.9961,  # Cr
        54.938044,  # Mn
        55.845,  # Fe
        58.933194,  # Co
        58.6934,  # Ni
        63.546,  # Cu
        65.38,  # Zn
        69.723,  # Ga
        72.630,  # Ge
        74.921595,  # As
        78.971,  # Se
        79.904,  # Br [79.901, 79.907]
        83.798,  # Kr
        85.4678,  # Rb
        87.62,  # Sr
        88.90584,  # Y
        91.224,  # Zr
        92.90637,  # Nb
        95.95,  # Mo
        97.90721,  # 98Tc
        101.07,  # Ru
        102.90550,  # Rh
        106.42,  # Pd
        107.8682,  # Ag
        112.414,  # Cd
        114.818,  # In
        118.710,  # Sn
        121.760,  # Sb
        127.60,  # Te
        126.90447,  # I
        131.293,  # Xe
        132.90545196,  # Cs
        137.327,  # Ba
        138.90547,  # La
        140.116,  # Ce
        140.90766,  # Pr
        144.242,  # Nd
        144.91276,  # 145Pm
        150.36,  # Sm
        151.964,  # Eu
        157.25,  # Gd
        158.92535,  # Tb
        162.500,  # Dy
        164.93033,  # Ho
        167.259,  # Er
        168.93422,  # Tm
        173.054,  # Yb
        174.9668,  # Lu
        178.49,  # Hf
        180.94788,  # Ta
        183.84,  # W
        186.207,  # Re
        190.23,  # Os
        192.217,  # Ir
        195.084,  # Pt
        196.966569,  # Au
        200.592,  # Hg
        204.38,  # Tl [204.382, 204.385]
        207.2,  # Pb
        208.98040,  # Bi
        208.98243,  # 209Po
        209.98715,  # 210At
        222.01758,  # 222Rn
        223.01974,  # 223Fr
        226.02541,  # 226Ra
        227.02775,  # 227Ac
        232.0377,  # Th
        231.03588,  # Pa
        238.02891,  # U
        237.04817,  # 237Np
        244.06421,  # 244Pu
        243.06138,  # 243Am
        247.07035,  # 247Cm
        247.07031,  # 247Bk
        251.07959,  # 251Cf
        252.0830,  # 252Es
        257.09511,  # 257Fm
        258.09843,  # 258Md
        259.1010,  # 259No
        262.110,  # 262Lr
        267.122,  # 267Rf
        268.126,  # 268Db
        271.134,  # 271Sg
        270.133,  # 270Bh
        269.1338,  # 269Hs
        278.156,  # 278Mt
        281.165,  # 281Ds
        281.166,  # 281Rg
        285.177,  # 285Cn
        286.182,  # 286Nh
        289.190,  # 289Fl
        289.194,  # 289Mc
        293.204,  # 293Lv
        293.208,  # 293Ts
        294.214,  # 294Og
    ]

    masses = Dict(k => atomic_masses_iupac2016[atomic_numbers[k]] for k in chemical_symbols)

    return chemical_symbols, atomic_numbers, masses
end

function load_spgs(;file::String="$(@__DIR__)/spacegroup.json")
    return JSON.parsefile(file)
end

function load_refs(;file::String="$(@__DIR__)/ase-atoms.json")
    return JSON.parsefile(file)
end

"""
    parse_json_spacegroup(spgs,nr,setting)

Formats the info of a specific spacegroup, preparing its use in `get_symops`.
"""
function parse_json_spacegroup(spgs::Dict,nr::T0,setting::T1) where {T0<:Integer,T1<:Integer}
    spg = spgs["$(nr): $(setting)"]
    spg["subtrans"] = [Array{Float64}(reshape(v,3)) for v in spg["subtrans"]]
    spg["translations"] = [Array{Float64}(reshape(v,3)) for v in spg["translations"]]
    spg["rotations"] = [Array{Float64}(hcat(v...)) for v in spg["rotations"]]
    return spg
end

"""
    get_symops(spg)

Generates rotation and translation operations for the given spacegroup info (`spg` from `parse_json_spacegroup`).
"""
function get_symops(spg::Dict)
    #  `spg` is the output of `parse_json_spacegroup`
    parities = spg["centrosymmetric"] ? [1,-1] : [1]
    symops = []
    @assert length(spg["rotations"]) == length(spg["translations"])
    for (parity, trans_sub) in Iterators.product(parities, spg["subtrans"])
        for (rot, trans) in Iterators.zip(spg["rotations"], spg["translations"])
            push!(symops, (parity * rot,
                           (trans + trans_sub) .% 1))
        end
    end
    return symops
end

"""
    fold(x)

Folding `x` back into [0,1).
"""
function fold(x::T) where T<:Any
    return x < 0 ? x + 1 : x
end

"""
    get_equivalent_sites(basis, symops)

Computes equivalent sites of in `basis` using rotations and translations in `symops`.
"""
function get_equivalent_sites(basis::Array{Array{Float64,2},1}, symops)
    kinds, sites = [], []

    for (kind, pos) in enumerate(basis)
        for (rot, trans) in symops
            site = (transpose(pos * rot) + trans) .% 1
            site = fold.(site)
            isdifferent = !any([v ≈ site for v in sites])
            if ((length(sites) == 0) | isdifferent)
                push!(sites, site)
                append!(kinds, kind)
            end
        end
    end
    return sites, kinds
end

"""
    make_unit_vec(x)

Computes the normalized `x` vector.
"""
function make_unit_vec(x::Array{T,1}) where T<:Any
    return x / norm(x)
end

"""
    get_coords(a_direction, ab_normal)

Computes the coordinate system given `a_direction` (usually [1 0 0]) and the vector normal to `a` and `b` (`ab_normal`, usually [0 0 1]).
"""
function get_coords(a_direction::Array{Float64,1}, ab_normal::Array{Float64,1})
    @assert dot(a_direction, ab_normal) ≈ 0.
    _x = make_unit_vec(a_direction)
    z = make_unit_vec(ab_normal)

    x = _x - dot(_x, ab_normal) * z
    xyz = hcat(x, cross(z,x), z)
    return xyz
end

"""
    deg2rad(x)

Convert degrees to radians.
"""
function deg2rad(x::T) where T <: Real
    return x * π / 180.
end

"""
    get_cos(x)

Compute the cosine of x degrees.
"""
function get_cos(x::T) where T <: Real
    return x ≈ 90 ? 0 : cos(deg2rad(x))
end

"""
    get_cell_vectors(cellpar)

Computes 3x3 cell vectors (column-wise). cellpar = [a, b, c, α, β, γ]
"""
function get_cell_vectors(cellpar::Array{Float64,1})
    a, b, c, α, β, γ = cellpar
    cos_α = get_cos(α)
    cos_β = get_cos(β)
    cos_γ = get_cos(abs(γ))
    sin_γ = abs(γ) ≈ 90 ? sign(γ) : sin(deg2rad(γ))
    cos_α, cos_β, cos_γ, sin_γ

    cy = (cos_α - cos_β * cos_γ) / sin_γ
    abc = hcat([a; 0; 0], b*[cos_γ; sin_γ; 0], c*[cos_β; cy; √(1-cos_β*cos_β-cy*cy)])
    return abc
end

struct CoordinateSystem{T}
    x::Array{T}
    y::Array{T}
    z::Array{T}
    M::Matrix{T} # contains x,y,z as column vectors
end

"""
    CartesianCoords(T)

Initialization of the cartesian coordiantes using `CoordianteSystem`.
"""
CartesianCoords(T) = CoordinateSystem{T}([1; 0; 0], [0; 1; 0], [0; 0; 1], I(3))

"""
    PrimitiveVectors(cc)

Generates vectors for a primitive cell on the basis of a coordinate system, like the `CartesianCoords`.
"""
function PrimitiveVectors(cc::CoordinateSystem{T};
        A₁=[1; 0; 0], A₂=[0; 1; 0], A₃=[0; 0; 1]) where {T<:Real}
    # x = SVector{3}(sum([cc.x, cc.y, cc.z] .* A₁))
    x = cc.M * A₁
    y = cc.M * A₂
    z = cc.M * A₃
    M = hcat(x,y,z)
    return CoordinateSystem{T}(x, y, z, M)
end

"""
    Atom(; name, mass)

Generates an atom object for the creation of unit cells.
"""
struct Atom{T}
    name::String
    mass::T
end

Atom(name,mass) = Atom{Float64}(name,mass)
make_atom(name,mass) = Atom{typeof(mass)}(name,mass)

function Atom(;name::String="dummy",mass::T=42) where T<:Real
    return Atom{typeof(mass)}(name, mass)
end

struct Spacegroup
    nr::Int32 #spacegroup number
    setting::Int8 # spacegroup setting (1 or 2, inherited from ase.spacegroup)
    kinds::Array{Int16,1}
    sites::Array{Any,1}
end

Spacegroup() = Spacegroup(-1,-1,[-1],[0])

struct Cell{T}
    atoms::Array
    positions
    box::CoordinateSystem{T}
    edge_lengths::Array{T}
    # spacegroup related infos
    spacegroup::Spacegroup
end

Cell(atoms, positions, box, edge_lengths) = Cell{Float64}(atoms, positions, box, edge_lengths, Spacegroup())
Cell(atoms, positions, box, edge_lengths, spg) = Cell{Float64}(atoms, positions, box, edge_lengths, spg)

"""
    make_unitcell(basis, symbols, nr, setting, cellpar; spg_file=nothing, a_direction=[1.; 0.; 0.], ab_normal=[0.;0.;1.], atom_cb=make_atom)

Creates a `Cell` struct including atom information, atom coordinates, the cell vectors/box and the edge lengths of the box.
"""
function make_unitcell(basis, symbols, nr, setting, cellpar;
        a_direction=[1.; 0.; 0.], ab_normal=[0.; 0.; 1.],
        atom_cb::Function=make_atom)
    # get constants
    chemical_symbols, atomic_numbers, masses = get_chemical_info()

    # get stored symmetry operation components
    spgs = load_spgs()
    spg = parse_json_spacegroup(spgs, nr, setting)

    # prepare symmetry operations for application to `basis`
    symops = get_symops(spg)

    # generating equivalent sites
    sites, kinds = get_equivalent_sites(basis, symops)

    # computing the simulation cell/box
    xyz = get_coords(a_direction, ab_normal)
    abc = get_cell_vectors(cellpar)
    cell = abc * xyz

    # storing crystal properties in a `SingleCrystal.Cell` instance
    el2atom_map = Dict(el => atom_cb(el, masses[el]) for el in symbols)

    cc = CartesianCoords(Float64)
    box = PrimitiveVectors(cc, A₁=abc[:,1], A₂=abc[:,2], A₃=abc[:,3])

    _spg = Spacegroup(nr, setting, kinds, sites)

    crystal = Cell(
        [el2atom_map[symbols[v]] for v in kinds],
        [cell * v for v in  sites],
        box,
        [norm(abc[:,1]), norm(abc[:,2]), norm(abc[:,3])],
        _spg
    )
    return crystal
end

"""
    make_bcc_unitcell(symbol, a)

Convenience function wrapping `make_unitcell` to create a bcc unit cell.
"""
make_bcc_unitcell(symbol::String,a::T,atom_cb::Function=make_atom) where T<:Real = make_unitcell([[0. 0. 0.],], [symbol], 229, 1, [a,a,a,90.,90,90], atom_cb=atom_cb)

"""
    make_fcc_unitcell(element, a)

Convenience function wrapping `make_unitcell` to create a fcc unit cell.
"""
make_fcc_unitcell(symbol::String,a::T,atom_cb::Function=make_atom) where T<:Real = make_unitcell([[0. 0. 0.],], [symbol], 225, 1, [a,a,a,90.,90,90], atom_cb=atom_cb)

"""
    add_vacancies(c, ixs; random=false, n_vac=1)

Adds a vacancy / removes an atom to/from a `Cell` struct.
"""
function add_vacancies(c::Cell; ixs::Array{Int}=[1], random::Bool=false, n_vac::Int=1)
    # if `random=true` then `n_vac` are generated and `ixs` is ignored
    n_atoms = length(c.atoms)
    if random == true
        ixs = rand(1:n_atoms, n_vac)
    end
    atoms_vac = [c.atoms[i] for i in 1:n_atoms if !(i in ixs)]
    positions_vac = [c.positions[i] for i in 1:n_atoms if !(i in ixs)]
    type = typeof(positions_vac[1][1])
    return Cell(atoms_vac, positions_vac, c.box, c.edge_lengths, c.spacegroup)
end

"""
    make_supercell(c;nx=1,ny=1,nz=1)

Creates a supercell by replicating the unitcell `c` and translating it along the vectors in `cell.box` `nx` * `ny` * `nz` times.
"""
function make_supercell(c::Cell; nx::I=1, ny::I=1, nz::I=1) where {I <: Integer}
    @assert (nx > 0) & (ny > 0) & (nz > 0)
    atoms = []
    positions = []
    box = PrimitiveVectors(c.box, A₁=[nx; 0; 0], A₂=[0; ny; 0], A₃=[0; 0; nz])
    edge_lengths = [norm(box.x), norm(box.y), norm(box.z)]

    for i in 0:nx-1, j in 0:ny-1, k in 0:nz-1
        shift = c.box.M * [i, j, k] # [c.box.x * i + c.box.y * j + c.box.z * k]
        new_pos = [pos + shift for pos in c.positions]
        push!(positions, new_pos) # TODO:  ((3, 1), Size(3,)) of input arrays do not match
        push!(atoms, c.atoms)
    end
    supercell = Cell(vcat(atoms...), vcat(positions...), box, edge_lengths)
    return supercell
end

"""
    vector1D(c1,c2,box_size)

Computes the periodic 1D distance between `c1` and `c2` with `box_size` as the periodic boundary.
"""
function vector1D(c1::Real, c2::Real, box_size::Real)
    if c1 < c2
        return (c2 - c1) < (c1 - c2 + box_size) ? (c2 - c1) : (c2 - c1 - box_size)
    else
        return (c1 - c2) < (c2 - c1 + box_size) ? (c2 - c1) : (c2 - c1 + box_size)
    end
end

"""
    vector(c1,c2,box_size)

Convenience function wrapping `vector1D` for the case of a cube.
"""
vector(c1, c2, box_size::Real) = vector1D.(c1, c2, box_size) # all box vectors orthogonal and of equal length

"""
    vector(c1,c2,edge_lengths)

Convenience function wrapping `vector1D` for the case of a non-cubic box with still orthogonal box vectors.
"""
vector(c1, c2, edge_lengths::Array) = vector1D.(c1, c2, edge_lengths) # all box vectors orthogonal

function parse_json_crystal(json_ase_crystal::Dict{String,Any})
    ase_crystal = Dict()
    ase_crystal["spacegroup_kinds"] = Array{Int32}(json_ase_crystal["spacegroup_kinds"])
    ase_crystal["positions"] = [Array{Float64}(v) for v in json_ase_crystal["positions"]]
    ase_crystal["cell"] = hcat(json_ase_crystal["cell"]...)
    ase_crystal["symbols"] = Array{String}(json_ase_crystal["symbols"])
    return ase_crystal
end

function cell2dict(crystal::SingleCrystal.Cell)

    chemical_symbols, atomic_numbers, masses = get_chemical_info()

    d = Dict()
    d["numbers"] = [atomic_numbers[atom.name] for atom in crystal.atoms] #[atomic_numbers[symbols[v]] for v in kinds]
    d["positions"] = hcat(crystal.positions...)
    d["spacegroup_kinds"] = crystal.spacegroup.kinds
    d["cell"] = crystal.box.M
    d["pbc"] = [true,true,true]
    d["info"] = Dict("spacegroup"=>Dict("nr"=>crystal.spacegroup.nr, "setting"=>crystal.spacegroup.setting))
    return d
end

end


text:
#!/usr/bin/env julia
using Printf
using Distributed
using TimerOutputs
using CUDAdrv
using CUDAnative
using CuArrays

function ζ_gen_range(mstart,mend,e)
    s = 0.0
    for i=mstart:-1:mend
      s += (1.0/i)^e
    end
    s
end

function ζ_gen_thr_range(mstart,mend,e)
    s = 0
    cls  = 8  # cache line size in units of Float64 storage size
              # so we can avoid false sharing
    Nthr = Threads.nthreads()
    psum = zeros(Float64,Nthr*cls)
    ms   = zeros(Float64,Nthr*cls)
    me   = zeros(Float64,Nthr*cls)

    # preparing each thread, start, end, and partial sums
    for i=1:Nthr
        ms[i*cls]   = (N/Nthr)*i
        me[i*cls]   = 1 + (N/Nthr)*(i-1)
        psum[i*cls] = 0.0
    end

    # actual function we will execute per thread
    # psums will be updated per thread id.  This is why
    # the psums are spaced 1 cache line apart, versus adjacent indices
    function inner!(psum,ms,me,x)
        ndx     = Threads.threadid()
        for i=ms[ndx*cls]:-1:me[ndx*cls]
            psum[ndx*cls] += (1.0/i)^e
        end
	nothing
    end

    # Threaded inner loop, each thread has no dependence upon others
    Threads.@threads for t=1:Nthr
        inner!(psum,ms,me,2)
    end

    # now sum up the partial sums and return s
    s = sum(psum)
    s
end

function ζ_gen_gpu_range(N)

    Nthr = 512
    st   = 8

    psum_d = CuArrays.fill(0.0::Float64, Nthr*st)

    # actual function we will execute per thread
    # psums will be updated per thread id.  This is why
    # the psums are spaced 1 cache line apart, versus adjacent indices
    function innergpu!(N,ps)
        ndx     = threadIdx().x
        stride  = blockDim().x

        for i=ndx:stride:N
		#for i=N:-stride:ndx
	    	r = 1.0/i
            ps[ndx] += r*r
        end
        return nothing
    end

    # Threaded inner loop, each thread has no dependence upon others
    @cuda threads=Nthr  innergpu!(N,psum_d)

    # now sum up the partial sums and return s
    s = sum(psum_d)
    s
end

function ζ_2_pow2_range(mstart,mend)
    s = 0.0
    for i=mstart:-1:mend
      s += (1.0/i)^2
    end
    s
end

function ζ_2_sqr_2_range(mstart,mend)
    s = 0.0
    for i=mstart:-1:mend
      s += (1.0/i)*(1.0/i)
    end
    s
end

function ζ_gen_simd_range(mstart,mend,e)
    s = 0.0
    @simd for i=mstart:-1:mend
      s += (1.0/i)^e
    end
    s
end

function ζ_2_pow2_simd_range(mstart,mend)
    s = 0.0
    @simd for i=mstart:-1:mend
      s += (1.0/i)^2
    end
    s
end

function ζ_2_sqr_simd_range(mstart,mend)
    s = 0.0
    @simd for i=mstart:-1:mend
      s += (1.0/i)*(1.0/i)
    end
    s
end

to = TimerOutput()

N = 16000

# first bit of code to handle the compilation stage, otherwise you
# see this in the time output
s = ζ_gen_range(N,1,2) + ζ_2_pow2_range(N,1) + ζ_2_sqr_2_range(N,1) +
    ζ_gen_simd_range(N,1,2) + ζ_2_pow2_simd_range(N,1) +
    ζ_2_sqr_simd_range(N,1) + ζ_gen_thr_range(N,1,2) + ζ_gen_gpu_range(N)

s /= 8.0
@printf("s=%f\n",s)
s=0.0

@timeit to "ζ_gen_range"            s+=ζ_gen_range(N,1,2)
@timeit to "ζ_2_pow2_range"         s+=ζ_2_pow2_range(N,1)
@timeit to "ζ_2_sqr_2_range"        s+=ζ_2_sqr_2_range(N,1)
@timeit to "ζ_gen_simd_range"       s+=ζ_gen_simd_range(N,1,2)
@timeit to "ζ_2_pow2_simd_range"    s+=ζ_2_pow2_simd_range(N,1)
@timeit to "ζ_2_sqr_simd_range"     s+=ζ_2_sqr_simd_range(N,1)
@timeit to "ζ_gen_thr_range"        s+=ζ_gen_thr_range(N,1,2)
@timeit to "ζ_gen_gpu_range"        s+=ζ_gen_gpu_range(N)

s /= 8.0
@printf("s=%f\n",s)

show(to)
@printf("\n");

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
# This problem attempts to reproduce the results from the article
# A General Multiple-Relaxation-Time Boltzmann Collision Model
# by Xiaowen Shan and Hudong Chen

# struct LinearHydrodynamicsModes <: FluidFlowProblem
#     rho_0::Float64
#     u_max::Float64
#     ν::Float64
#     NX::Int64
#     NY::Int64
#     domain_size::Tuple{Float64, Float64}
# end

"""
Models thermal diffusion given by,

ρ = ρ_0 + ρ̃ sin(2pi y / Ly)
u = 0
θ = ρ_0 θ_0 / ρ

It is expected that the temperature θ decays at the rate of exp(- κ t)
"""
struct LinearizedThermalDiffusion <: FluidFlowProblem
    ρ_0::Float64
    θ_0::Float64

    ρ̃::Float64
    ũ::Float64

    ν::Float64
    κ::Float64

    domain_size::Tuple{Float64, Float64}

    u_max::Float64
    NX::Int64
    NY::Int64
end
function LinearizedThermalDiffusion(ν, κ, scale, NY = 4 * scale)
    u_max = 0.01 / scale

    return LinearizedThermalDiffusion(
        1.0,
        1.0,
        0.001,
        0.001,
        ν,
        κ,
        (2pi, 2pi),
        u_max,
        # 3,
        NY,
        NY,
    )
end

function delta_x(problem::LinearizedThermalDiffusion)
    return problem.domain_size[2] / problem.NY
end

function density(
    q::Quadrature,
    problem::LinearizedThermalDiffusion,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    ρ_0 = problem.ρ_0
    ρ̃ = problem.ρ̃

    return ρ_0 + ρ̃ * sin(y) * exp(-heat_diffusion(problem) * timestep)
end

function pressure(
    q::Quadrature,
    problem::LinearizedThermalDiffusion,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    ρ_0 = problem.ρ_0
    θ_0 = problem.θ_0

    # The constant pressure will enforce a temperature given by
    # θ = ρ_0 * θ_0 / ρ
    cs = q.speed_of_sound_squared
    return ρ_0 * θ_0
end

function velocity(
    problem::LinearizedThermalDiffusion,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    return [
        0.0
        0.0
    ]
end

"""
Models the lienarized transverse shear wave,

ρ = ρ_0
u = [ ũ sin(2pi y / Ly), 0 ]
θ = θ_0

It is expected that the temperature θ decays at the rate of exp(- κ t)
"""
struct LinearizedTransverseShearWave <: FluidFlowProblem
    ρ_0::Float64
    θ_0::Float64

    ρ̃::Float64
    ũ::Float64

    ν::Float64
    κ::Float64

    domain_size::Tuple{Float64, Float64}

    u_max::Float64
    NX::Int64
    NY::Int64
end
function LinearizedTransverseShearWave(ν, κ, scale, NY = 8 * scale)
    u_max = 0.01 / scale

    return LinearizedTransverseShearWave(
        1.0,
        1.0,
        0.001,
        0.001,
        ν,
        κ,
        (2pi, 2pi),
        u_max,
        NY,
        NY,
    )
end

function delta_x(problem::LinearizedTransverseShearWave)
    return problem.domain_size[2] / problem.NY
end

function density(
    q::Quadrature,
    problem::LinearizedTransverseShearWave,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    return problem.ρ_0
end

function pressure(
    q::Quadrature,
    problem::LinearizedTransverseShearWave,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    return problem.θ_0
end

function velocity(
    problem::LinearizedTransverseShearWave,
    x::Float64,
    y::Float64,
    timestep::Float64 = 0.0,
)
    return [
        problem.ũ * sin(y) * exp(-viscosity(problem) * timestep)
        0.0
    ]
end


text:
module SpringCollab2020NoDashRefillSpring

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/NoDashRefillSpring" NoDashRefillSpring(x::Integer, y::Integer, playerCanUse::Bool=true)
@mapdef Entity "SpringCollab2020/NoDashRefillSpringRight" NoDashRefillSpringRight(x::Integer, y::Integer)
@mapdef Entity "SpringCollab2020/NoDashRefillSpringLeft" NoDashRefillSpringLeft(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "No Dash Refill Spring (Up) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        NoDashRefillSpring
    ),
    "No Dash Refill Spring (Left) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        NoDashRefillSpringRight
    ),
    "No Dash Refill Spring (Right) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        NoDashRefillSpringLeft
    ),
)

function Ahorn.selection(entity::NoDashRefillSpring)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

function Ahorn.selection(entity::NoDashRefillSpringLeft)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

function Ahorn.selection(entity::NoDashRefillSpringRight)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

sprite = "objects/SpringCollab2020/noDashRefillSpring/00.png"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpring, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpringLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 24, 0, rot=pi / 2)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashRefillSpringRight, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, -8, 16, rot=-pi / 2)

end


text:
@testset "loss" begin
	include("softmaxLossTest.jl")
end


text:
# tests related to functional programming functions and styles

# map -- array.jl
@test isequal(map((x)->"$x"[end:end], [9:11]), ["9", "0", "1"])
# TODO: @test map!()
# map -- ranges.jl
@test isequal(map(i->sqrt(i), 1:5), [sqrt(i) for i in 1:5])
@test isequal(map(i->sqrt(i), 2:6), [sqrt(i) for i in 2:6])

# map on ranges should evaluate first value only once (#4453)
let io=IOBuffer(3)
    map(x->print(io,x), 1:2)
    @test takebuf_string(io)=="12"
end

# map over [] should return []
# issue #6719
@test isequal(typeof(map(x -> x, [])), Array{None,1})

# maps of tuples (formerly in test/core.jl) -- tuple.jl
@test map((x,y)->x+y,(1,2,3),(4,5,6)) == (5,7,9)
@test map((x,y)->x+y,
            (100001,100002,100003),
            (100004,100005,100006)) == (200005,200007,200009)

# maps of strings (character arrays) -- string.jl
@test map((c)->char(c+1), "abcDEF") == "bcdEFG"

# filter -- array.jl
@test isequal(filter(x->(x>1), [0 1 2 3 2 1 0]), [2, 3, 2])
# TODO: @test_throws isequal(filter(x->x+1, [0 1 2 3 2 1 0]), [2, 3, 2])
@test isequal(filter(x->(x>10), [0 1 2 3 2 1 0]), [])
@test isequal(filter((ss)->length(ss)==3, ["abcd", "efg", "hij", "klmn", "opq"]), ["efg", "hij", "opq"])

# zip and filter iterators
# issue #4718
@test collect(filter(x->x[1], zip([true, false, true, false],"abcd"))) == {(true,'a'),(true,'c')}

# enumerate (issue #6284)
let b = IOBuffer("1\n2\n3\n"), a = {}
    for (i,x) in enumerate(eachline(b))
        push!(a, (i,x))
    end
    @test a == {(1,"1\n"),(2,"2\n"),(3,"3\n")}
end

# zip eachline (issue #7369)
let zeb     = IOBuffer("1\n2\n3\n4\n5\n"),
    letters = ['a', 'b', 'c', 'd', 'e'],
    res     = {}
    for (number, letter) in zip(eachline(zeb), letters)
        push!(res, (int(strip(number)), letter))
    end
    @test res == {(1, 'a'), (2, 'b'), (3, 'c'), (4, 'd'), (5, 'e')}
end

max_stars_repo_name: jey/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
using BinaryLambdaCalculus
using Test

@testset "Encoding/decoding" begin
    for t in terms(0, 10)
        @test decode(encode(t)) == t
        @test decompress(compress(t)) == t
    end
end

@testset "Enumeration" begin
    function S(m::T, n::T) where T
        # explicit recursive formula for testing
        if n <= 1
            return 0
        else
            return T(m >= n - 1) + S(m + 1, n - 2) + sum(S(m, k) * S(m, n - 2 - k) for k = 0:(n-2))
        end
    end

    S₀ₙ = [0, 0, 0, 0, 1, 0, 1, 1, 2, 1, 6, 5, 13, 14, 37, 44, 101, 134, 298, 431]
    @test all(S₀ₙ[i+1] == S(0, i) for i = 0:19) # ensure the reference is correct...
    @test [grygiel_lescanne(0, n) for n = 0:19] == S₀ₙ

    for m = 1:3, n = 0:20
        @test grygiel_lescanne(m, n) == S(m, n)
        @test length(terms(m, n)) == S(m, n)
    end

    let m = 3, n = 100,
        s = (a, b) -> grygiel_lescanne(a, b, BigInt)
        @test s(m, n) == s(m + 1, n - 2) + sum(s(m, k) * s(m, n - 2 - k) for k = 0:(n-2))
    end
end


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

using LLVM
using LLVM.Interop

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

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

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

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

# make Julia compile and execute the function
push!(function_attributes(sum), EnumAttribute("alwaysinline"))
ptr = LLVM.ref(sum)
call_sum(x, y) = Base.llvmcall(ptr, Int32, Tuple{Int32, Int32}, x, y)
@show call_sum(x, y)


text:
######
##### Default Definitions
#####

const CachedTypes = Union{Object,CachedArray}

function setdirty!(x::CachedTypes, flag = true)
    meta = metadata(x)
    setdirty!(manager(x), meta, flag)
    setdirty!(meta, flag)
    return nothing
end

update!(x::CachedTypes) = update!(manager(x), metadata(x))
isdirty(x::CachedTypes) = isdirty(metadata(x))
id(x::CachedTypes) = getid(metadata(x))
pool(x::CachedTypes) = getpool(metadata(x))

function prefetch!(x::CachedTypes; kw...)
    _manager = manager(x)
    @spinlock alloc_lock(_manager) begin
        block = metadata(x)
        # Possibility that `maybe_cleanup!` cleans up this block.
        # Presumeably, this can't happen since we have `x` and `x` is holding onto the
        # object that owns this block ... so it should never have been garbage collected
        # in the first place.
        if maybe_cleanup!(_manager, getid(block)) == false
            prefetch!(block, getpolicy(_manager), _manager; kw...)
        end
    end
end

function evict!(x::CachedTypes; kw...)
    _manager = manager(x)
    @spinlock alloc_lock(_manager) begin
        maybe_cleanup!(_manager)
        evict!(x, getpolicy(_manager), _manager)
    end
    return nothing
end

function softevict!(x::CachedTypes)
    _manager = manager(x)
    @spinlock alloc_lock(_manager) softevict!(getpolicy(_manager), _manager, metadata(x))
end

#####
##### Annotation API
#####

function onobjects(f::F, x::Array{T}) where {F,T}
    isbitstype(T) || foreach(x -> onobjects(f, x), x)
    return nothing
end

onobjects(f::F, x::Object) where {F} = isnull(unsafe_pointer(x)) || f(x)
onobjects(f::F, x::CachedArray) where {F} = onobjects(f, x.object)
onobjects(f::F, x::Union{NamedTuple,Tuple}) where {F} = foreach(x -> onobjects(f, x), x)

"""
    onobjects(f, x)

Call function `f` on all [`CachedArrays.Objects`] that are reachable from `x`.
If `f` is a composite type, recursively calls `onobjects` on each member of `x`.
If `f` is some iterable like an `AbstractArrays` or `Tuple`, then `onobjects` is called on
each member of that array.

Behavior can be customized by extending `onobjects` for your custom types.
If you would like to opt out your type from `onobjects` altogether, the convenience macro
[`@blockobjects`] may be used.
"""
@generated function onobjects(f::F, x::T) where {F,T}
    (isbitstype(T) || iszero(fieldcount(T))) && return :()
    exprs = map(fieldnames(T)) do fieldname
        return :(onobjects(f, Core.getfield(x, $(QuoteNode(fieldname)))))
    end
    return quote
        $(exprs...)
    end
end

"""
    @blockobjects MyType

Restrict [`onobjects`](@ref) from recursing into the fields of `MyType`.
"""
macro blockobjects(type)
    return :(onobjects(f::F, _::$(esc(type))) where {F} = nothing)
end

#####
##### Recursively gather all objects
#####

function findobjects_builder!(q::Expr, tags::Vector{Any}, sym::Symbol, ::Type{T}) where {T}
    gf = GlobalRef(Core, :getfield)
    for f in Base.OneTo(fieldcount(T))
        TF = fieldtype(T, f)
        # Note: Objects are mutable structs.
        # Therefore, no `isbitstype` can contain a reference to an Object, allowing us
        # to skip all these types.
        skip = !Base.isconcretetype(TF) || Base.isbitstype(TF) || iszero(fieldcount(TF))
        skip && continue
        gfcall = Expr(:call, gf, sym, f)
        newsym = gensym(sym)
        # If this field type is an object, we're done.
        # Emit the expression to make the appropriate assignment to `newsym` and
        # record `newsym` as one of the symbols to use to build the final tuple.
        if TF <: Object
            push!(q.args, Expr(:(=), newsym, gfcall))
            push!(tags, newsym)
        else
            # Otherwise, there's the possiblity that THIS member is also a composite
            # struct that might contain an Object, so we need to recurse.
            # However, we don't want to needlessly bloat our generated function.
            #
            # SO, only append the sub-expressions if at least one Object is found during
            # this recursive call (i.e., `newtags` is not empty)
            newtags = []
            newq = :()
            findobjects_builder!(newq, newtags, newsym, TF)
            if !isempty(newtags)
                push!(q.args, Expr(:(=), newsym, gfcall))
                append!(q.args, newq.args)
                append!(tags, newtags)
            end
        end
    end
end

function findobjects_builder(::Type{T}) where {T}
    q = :()
    tags = []
    findobjects_builder!(q, tags, :x, T)
    return q, tags
end

"""
    findobjects(x) -> (Objects...,)
    findobjects_builder!(q::Expr, tags::Vector{Symbol}, sym::Symbol, ::Type{T})

Recursive compile-time function builder expression.
The goal of this function is to assist in building an expression that:

1. Builds up `getfield` calls on nested structs.
2. Records the LHS symbols of the `getfield` forest that correspond to `CachedArrays.Object`s
3. Removes misc. and unnecessary `getfield` calls to reduce the work required for the
    downstream compiler.

It achieves by reflecting on the type `T` and selectively expanding `q` and `tags`.

** WARNING **: A current limitations are:

    1. All recursively nested types should be concrete. If they aren't, they will be
       skipped, which may or may not be what you want.
    2. Arrays will not be recursed over.
"""
@generated function findobjects(x::T) where {T}
    body = Expr(:block)
    tags = []
    findobjects_builder!(body, tags, :x, T)
    tuple = :(($(tags...),))
    push!(body.args, tuple)
    return body
end

#####
##### noescape
#####

function noescape(
    manager::CacheManager, ::Val{Large}, f::F, args::Vararg{Any,N}; kw...
) where {Large,F,N}
    # Record the allocation ID in the manager before and after the call to `f`.
    # We the use `findobjects` on the returned result to gather any `Objects` that
    # might be present, extracting the `id` for these objects.
    start = readid(manager)
    result = f(args...; kw...)
    if Large
        saveids = Set{UInt}()
        onobjects(result) do o
            isnull(unsafe_pointer(o)) || push!(saveids, getid(metadata(o)))
        end
    else
        objects = findobjects(result)
        saveids = map(objects) do o
            # If an allocation happens to be for zero bytes, then CachedArrays
            # will return a null pointer.
            #
            # Thus, we need to perform a null check here in order to avoid a segfault in
            # these corner cases.
            return isnull(unsafe_pointer(o)) ? zero(UInt) : getid(metadata(o))
        end
    end
    stop = readid(manager)

    # For all allocations that were made between `start` and `stop`, perform
    # `unsafe_free` on those objects (if not present in the returned objects).
    @spinlock alloc_lock(manager) begin
        backedges = manager.map.dict
        for id in Iterators.filter(!in(saveids), start:stop)
            backedge = get(backedges, id, nothing)
            backedge === nothing && continue

            # Back edge exists - free the result
            unsafe_free(manager, backedge)
        end
    end
    return result
end

# So we can handle calls to our generated "__invoke__", the top level expression needs to
# either:
#
# 1. Just be a straight up function call.
# 2. Contain a call to `Base.invoke`.
#
# In the case of the latter, we only modify the outermost invoke.
function find_invoke(manager, annotations, expr::Expr)
    # Use a `Ref` for `success` to make the linter happy.
    processed_invoke = Ref(false)
    processed_tuple = Ref(false)
    newexpr = MacroTools.prewalk(expr) do subexpr
        processed_invoke[] && return subexpr
        if isa(subexpr, Expr) && subexpr.head == :tuple && processed_tuple[] == false
            processed_tuple[] = true
            return esc(subexpr)
        elseif isa(subexpr, Expr) && subexpr.head == :call && subexpr.args[1] == invoke
            args, kw = extract_kw(subexpr.args[2:end], 1)

            # Perform the proper nested escaping because dealing with nested macros
            # is SUPER annoying.
            #
            # In all honesty, maybe we should teach `@annotate` how to correctly deal with
            # calls to `@noescape` to simplify the whole implementation.
            args[1] = esc(args[1])
            kw = map(esc, kw)

            processed_invoke[] = true
            return :(noescape($manager, $(annotations...), invoke, $(args...); $(kw...)))
        end
        return subexpr
    end
    if processed_invoke[] == false || processed_tuple[] == false
        error("Could not find a suitable subexpression!")
    end
    return newexpr
end

function extract_kw(args, pos = 2)
    if args[pos] isa Expr && args[pos].head == :parameters
        kw = args[pos].args
        deleteat!(args, pos)
    else
        kw = Any[]
    end
    return args, kw
end

function noescape_impl(manager, params, fn)
    large = params[:large]
    annotations = (Val(large),)
    manager = esc(manager)
    # Simple case
    if fn.head == :call
        args, kw = extract_kw(fn.args)
        args = map(esc, args)
        kw = map(esc, kw)
        return :(noescape($manager, $(annotations...), $(args...); $(kw...)))
    else
        # More complicated case.
        return find_invoke(manager, annotations, fn)
    end
end

default_noescape_kw() = Dict(:large => false)

"""
    @noescape manager [args=false] [large=false] f(args...; kw...)

Inform the CacheManager `manager` that intermediate allocations made during the call
to `f` that are not present in the returned value are save to be reclaimed.

A trailing `tuple` argument `others` can provide a further list of types to inspect in order
to keep otherwise escaping types from being freed.
"""
macro noescape(manager, fn)
    return noescape_impl(manager, default_noescape_kw(), fn)
end

macro noescape(manager, args...)
    fn = last(args)

    # Default keyword arguments
    kw = default_noescape_kw()
    kw_type_map = Dict(k => typeof(v) for (k, v) in kw)
    for i in Base.OneTo(lastindex(args) - 1)
        arg = args[i]
        if arg.head != :(=)
            error("Expected keyword arguments to be in the form `a=b`!")
        end

        # Match implemented keywords
        name = arg.args[1]
        value = arg.args[2]
        expected_type = get(kw_type_map, name, nothing)
        expected_type === nothing && error("Unknown keyword argument $(name)!")
        if !isa(value, expected_type)
            msg = """
            Expected keyword argument $name to have type $expected_type.
            Instead, it has type $(typeof(value))!
            """
            error(msg)
        end
        kw[name] = value
    end
    # Process keywords
    return noescape_impl(manager, kw, fn)
end


text:
"""
Find agent in the model given the id
"""
function find_agent(id,model)
    first(filter(a->a.id == id,Set(allagents(model))))
end

function get_control_agent(asset_id::Int,model)
    return model.mapping_ctl_ntw[asset_id]
end

function get_controlled_assets(agent_id::Int,model)::Set{Int64}
    assets = filter(k->model.mapping_ctl_ntw[k] == agent_id,keys(model.mapping_ctl_ntw))
    #log_info("assets controlled by $(agent_id) are: $(length(assets))")
    return assets
end

function has_active_controlled_assets(agent::Agent,model)
    assets = get_controlled_assets(agent.id,model)

    sum_up = sum([ is_up(getindex(model,sne)) for sne in assets ])
    return sum_up > 0 ? true : false
end

function set_control_agent!(asset_id::Int, agent_id::Int, model)
    getindex(model,asset_id).controller_id = agent_id
    #TODO Consider removing this line below
    #To avoid putting info in model
    model.mapping_ctl_ntw[asset_id] = agent_id
end

"""
Given a SimNE id it returns its ntw node address.
"""
function get_address(eid::Int,g::AbstractGraph)::Int
    #res = filter(p->p[2] == eid,pairs(model.mapping_ntw_sne))
    #return !isempty(res) ? first(keys(res)) : -1
    return g[eid,:eid]
end

"""
Given a ntw node address it returns the corresponding SimNE id
"""
function get_eid(address::Int,model)::Int
    return model.mapping_ntw_sne[address]
end

function update_sne_address!(eid::Int,new_address::Int,model)
    #curr_address = get_address(eid,model)
    model.mapping_ntw_sne[new_address] = eid
end

"""
Update (ntw node) addresses of SimNE agents after removal of a given SimNE
"""
function update_addresses_removal!(dpn_id::Int,model)
    available_addr = get_address(dpn_id,model.ntw_graph)
    #log_info("Current length of g: $(nv(model.ntw_graph))")
    for addr::Int=available_addr:nv(model.ntw_graph)
        #log_info("Address $addr and its type: $(typeof(addr))")
        update_sne_address!(
            get_eid(addr+1,model),
            addr,
            model
            )
    end
    delete!(model.mapping_ntw_sne,length(model.ntw_graph)+1)
end

function circle_shape(h,k,r)
    θ = LinRange(0,2*π, 500)
    h .+ r*sin.(θ), k .+ r*cos.(θ)
end

function do_agent_step!(a::SimNE,model)
    #Process OF messages (packet data traffic)
    # log_info(model.ticks,a.id, "start step! $(get_state(a).up) ==> $(get_state(a).rul)")
    # log_info(model.ticks,a.id,12, "====> ports: $(get_port_edge_list(a))<====")
    is_up(a) && is_ready(a) ? in_packet_processing(a,model) : nothing
    clear_pending_query!(a,model)
    # log_info(model.ticks,a.id,"rqsted: $(a.requested_ctl)")
    # @debug("[$(model.ticks)]($(a.id)) end step")
    # if model.ticks == 100
    #     a.maintenance.deterioration_parameter = 0.7
    # end
    deteriorate!(a,model)
end

function do_agent_step!(a::Agent,model)
    #Schedule events regardless of state of ctl agent to make sure that controlled snes are brought back up
    #This is the first function to make sure that nodes that are down are recognised by others.
    do_events_step!(a,model)
    if is_up(a) && length(a.queue.data) > 0
        # for sprt in sne_print
        #     log_info(model.ticks,a.id," step!: {$(sprt.id)} $(get_state(sprt).flow_table) ===> all ports: $(get_port_edge_list(sprt)) ===> paths: $(a.paths)")
        # end

        ## Process OF Messages (SimNE to (sdn) control messages)
        is_ready(a) ? in_packet_processing(a,model) : log_info("queue of $(a.id) is empty")
        # Process inter-agent messages
        # log_info(model.ticks,a.id,"==> a.paths ==> $(a.paths)")
        do_receive_messages(a,model)

        # log_info(model.ticks,a.id,"ctl_agent step: $(neighbors(a.ntw_graph,4))")
    end
    # log_info(model.ticks,a.id,25,"pending msgs: $(length(a.
    # pending)) --> $(a.pending)")
    # log_info(model.ticks,a.id,25,"QUEUE --> $(a.queue.data)")

    # do_confidence_check!(a,model)

    # if !isempty(get_state(a).active_paths)
    #     log_info(model.ticks,a.id,"ap-->$(get_state(a).active_paths)")
    # end

end

"""
    Receives inter-agent messages
"""
function do_receive_messages(a::Agent,model)
    #TODO: check if another periodicity is required, rather than every tick
    #if !isempty(a.msgs_in) log_info("[$(model.ticks)]($(a.id)) in msgs: $(a.msgs_in)") end

    #senders = [ m.sid for m in a.msgs_in ]
    #log_info("[$(model.ticks)]($(a.id)) has $(length(a.msgs_in)) msgs to process from $senders" )

    if model.ctrl_model != GraphModel(1)
        for msg in a.msgs_in
            #log_info(msg)
            process_msg!(a,msg,model)
        end
    end
end

"""
Records benchmark for a query
"""
function record_benchmark!(bdir,run_label,aid,query_time,query,query_graph,query_paths,benchmark)
    if benchmark
        if !isdir(bdir)
        mkdir(bdir)
        end
        #benchmark block start
        b = @benchmark begin
            do_query($query_time,$query,$query_graph,$query_paths)
        end
        serialize( bdir * "/"* run_label *"_$(first(query))_$(last(query))_$(query_time)_$(aid)_bchmk.bin",b)
    end
end

## main functions

new_config(seed,ctl_model,ntw_topo,size,n_steps,drop_proportion,prob_random_walks,benchmark, animation,k,B,ctl_k,ctl_B,mnt_policy,ntw_services,mnt_wc_duration,mnt_bc_duration,mnt_wc_cost,mnt_bc_cost,traffic_dist_params,traffic_proportion,clear_cache_graph_freq,interval_tpt,pkt_size,pkt_per_tick,max_queue_ne,ofmsg_reattempt, max_cache_paths) =
    return ( seed = seed
            ,ctl_model=ctl_model
            ,ntw_topo = ntw_topo
            ,size=size
            ,n_steps=n_steps
            ,drop_proportion=drop_proportion
            ,prob_random_walks = prob_random_walks
            ,benchmark = benchmark
            ,animation = animation
            ,k=k
            ,B = B
            ,ctl_k=ctl_k
            ,ctl_B = ctl_B
            ,custom_topo = nothing
            ,ctl_custom_topo = nothing
            ,mnt_policy = mnt_policy
            ,ntw_services = ntw_services
            ,mnt_wc_duration = mnt_wc_duration
            ,mnt_bc_duration = mnt_bc_duration
            ,mnt_wc_cost = mnt_wc_cost
            ,mnt_bc_cost = mnt_bc_cost
            ,traffic_dist_params = traffic_dist_params
            ,traffic_packets = traffic_packets
            ,link_capacity = link_capacity
            ,clear_cache_graph_freq = clear_cache_graph_freq
            ,interval_tpt = interval_tpt
            ,pkt_size = pkt_size
            ,pkt_per_tick = pkt_per_tick
            ,max_queue_ne = max_queue_ne
            ,ofmsg_reattempt = ofmsg_reattempt
            ,max_cache_paths = max_cache_paths
            )

function get_dropping_nodes(drop_proportion)
    #TODO calcualte according to proportion
    return Dict(
       50 => [3]
        #55: Grid and Grid fails
        , 75 => [8]
        #, 115 => [11]
        # , 70 => [10]
        #50=>[3,6] #TODO simoultaneous drops only work in centralised control
    #,120=>[2]
    ) # drop time => drop node
end

"""
It returns the next event time for a given random number and rate of events
"""
function next_event_time(rn,λ)
    return -log(1.0-rn)/λ
end

"""
Clears cache of control agent
# TODO Opportunity to investigate ways to store relevant cache rather than clearing all
"""
function clear_cache!(a::Agent,model::ABM)
    if model.ticks - a.last_cache_cleared == model.clear_cache_graph_freq
        a.ntw_graph = a.base_ntw_graph
        a.last_cache_cleared = model.ticks
        a.previous_queries = Dict()
        a.paths = Dict()
    end

end

function clear_cache!(sne::SimNE,model::ABM)
    #placeholder
end

function to_string(s)
    sep = "; "
    return join([getfield(s,a) for a in fieldnames(typeof(s))],sep)
end

"""
Return SNEs that are up
"""
function get_live_snes(model)
    return [ sne.id for sne in allagents(model) if typeof(sne) == SimNE && get_state(sne).up  ]
end
"""
Return SNEs controlled by Agent a that are up
"""
function get_live_snes(a::Agent,model)
    controlled_snes = [ getindex(model,sid) for sid in get_controlled_assets(a.id,model)]
    return [ sne.id for sne in controlled_snes if get_state(sne).up  ]
end

# function init_logger(log_name)
#     loggers[log_name] = getlogger(log_name)
#     setlevel!(loggers[log_name], "info")
#     push!(loggers[log_name],getlogger(name="root"))
#     # push!(loggers[log_name], DefaultHandler(tempname(), DefaultFormatter("[{date} | {level} | {name}]: {msg}")))
#     return loggers[log_name]
# end


text:
using BinaryBuilder, Pkg

name = "SPIRV_LLVM_Translator"
repo = "https://github.com/KhronosGroup/SPIRV-LLVM-Translator.git"

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# Bash recipe for building across all platforms
script = raw"""
cd SPIRV-LLVM-Translator
install_license LICENSE.TXT

CMAKE_FLAGS=()

# Release build for best performance
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

# Install things into $prefix
CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})

# Explicitly use our cmake toolchain file and tell CMake we're cross-compiling
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_CROSSCOMPILING:BOOL=ON)

# Tell CMake where LLVM is
CMAKE_FLAGS+=(-DLLVM_DIR="${prefix}/lib/cmake/llvm")

# Build the library
CMAKE_FLAGS+=(-DBUILD_SHARED_LIBS=ON)

cmake -B build -S . -GNinja ${CMAKE_FLAGS[@]}
ninja -C build -j ${nproc} llvm-spirv install
install -Dm755 build/tools/llvm-spirv/llvm-spirv${exeext} -t ${bindir}
"""

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libLLVMSPIRVLib", "LLVMSPIRVLib"], :libLLVMSPIRV),
    ExecutableProduct("llvm-spirv", :llvm_spirv),
]

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
"""
    function update!(
        pcs::PotentialCalculationSetup{T};
        ::Nothing, # these two unused arguments are used such that the method
        ::Any;     # is similar to the GPU method for it.
        use_nthreads::Int = Base.Threads.nthreads(),
        depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
        is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
        only2d::Val{only_2d} = Val{false}()
)::Nothing where {T, depletion_handling_enabled, only_2d, _is_weighting_potential}

This function performs one iteration of the SOR. One iteration consists out of 4 steps:

    1) Iterate in parallel over all even points and update their potential.
    2) Apply the boundary conditions at the ends of the grid for all even points.
    3) Iterate in parallel over all odd points and update their potential.
    4) Apply the boundary conditions at the ends of the grid for all odd points.
"""
@inline function update!(
    pcs::PotentialCalculationSetup{T},
    ::Nothing, # these two unused arguments are used such that the method
    ::Any;     # is similar to the GPU method for it.
    use_nthreads::Int = Base.Threads.nthreads(),
    depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
    is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
    only2d::Val{only_2d} = Val{false}(),
)::Nothing where {T, depletion_handling_enabled, only_2d, _is_weighting_potential}
    outerloop!(pcs, use_nthreads, Val{true}(), depletion_handling, is_weighting_potential, only2d)
    apply_boundary_conditions!(pcs, Val{true}(), only2d)
    outerloop!(pcs, use_nthreads, Val{false}(), depletion_handling, is_weighting_potential, only2d)
    apply_boundary_conditions!(pcs, Val{false}(), only2d)
    nothing
end


text:
function hard_drop(cur_player::AbstractPlayer)
  cur_piece = cur_player.piece

  isnull(cur_piece) && return false

  did_move = hard_drop(cur_player, cur_piece, true)

  cur_player.clock.lock =
    Nullable{Base.Random.UUID}()

  glue_piece!(cur_player)

  did_move
end

function hard_drop(cur_player::AbstractPlayer, cur_piece::AbstractPiece, is_main_piece::Bool=false)
  cur_grid = cur_player.grid

  max_rows = cur_grid.rows
  max_rows += cur_grid.hidden_rows
  max_rows += 1 # for safe measure

  did_move = drop!(
    cur_player, cur_piece, cur_grid,
    -max_rows, is_main_piece, false
  )

  did_move
end


text:
using StingerGraphs
using Base.Test

function bfsutil(scale::Int64, edgefactor::Int64)
    graph = kronecker(scale, edgefactor)
    s = Stinger()
    for i in 1:size(graph, 2)
        insert_edge!(s, 0, graph[1, i], graph[2, i], 1, 1)
    end
    nv = get_nv(s)
    for src in 0:nv-1
        parents = bfs(s, src)
        bfstest(s, parents)
    end
    return s
end

function bfstest(s::Stinger, parents::Array{Int64, 1}; weight::Int64=1, etype::Int64=0)
    for i in 1:size(parents, 1)
        #Check if edge exists
        if parents[i]>-1
            @test edgeweight(s, parents[i], i-1, etype) != 0
        end
    end
end


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf
# Uses DFTK in place of LJ for the production stage as a proof-of-concept ab initio MD simulation
# Note that the choice of parameters is for demonstration purposes, and the results are non-physical

using Atomistic
using AtomsBase
using DFTK
using NBodySimulator
using Unitful
using UnitfulAtomic

setup_threading(n_blas = 4)

N = 8
element = :Ar
box_size = 1.5u"nm" # this number was chosen arbitrarily
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

eq_steps = 20000
eq_thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
eq_simulator = NBSimulator(Δt, eq_steps, thermostat = eq_thermostat)
potential = LennardJonesParameters(1.657e-21u"J", 0.34u"nm", 0.765u"nm")

eq_result = @time simulate(initial_system, eq_simulator, potential)

display(@time plot_temperature(eq_result, eq_simulator.steps ÷ 200))
display(@time plot_energy(eq_result, eq_simulator.steps ÷ 200))
display(@time plot_rdf(eq_result, potential.σ, 0.5))

ab_initio_steps = 200
ab_initio_simulator = NBSimulator(Δt, ab_initio_steps, t₀ = get_time_range(eq_result)[end])
dftk_potential = DFTKPotential(
    psp = ElementPsp(:Ar, psp = load_psp(list_psp(:Ar, functional = "lda")[1].identifier)),
    lattice = box_size * [[1.0 0 0]; [0 1.0 0]; [0 0 1.0]],
    Ecut = 5u"hartree", # very non-physical but fast for demonstration purposes
    kgrid = [1, 1, 1],
    damping = 0.7,
    mixing = LdosMixing()
)

ab_initio_result = @time simulate(get_system(eq_result), ab_initio_simulator, dftk_potential)

# Plotting on separate plots because the timespan is so much smaller than in the first phase
display(plot_temperature(ab_initio_result, 1))
display(plot_energy(ab_initio_result, 1))
display(@time plot_rdf(ab_initio_result, potential.σ))

write_nbs_animation(ab_initio_result, "artifacts/argon_ab_initio.gif")
write_ase_trajectory(ab_initio_result, dftk_potential.psp, dftk_potential.lattice, "artifacts/argon_ab_initio.traj")

;


text:
# this file contains some useful functions

"""
    iphase(n::Integer)
``(-1)^n``
"""
iphase(n::T) where {T <: Integer} = iseven(n) ? one(T) : -one(T)

"""
    is_same_parity(x::T, y::T) where {T <: Integer}
judge if two integers are same odd or same even
"""
is_same_parity(x::T, y::T) where {T <: Integer} = iseven(x ⊻ y)

"""
devide integer `n` with 2, and return how to show it,
for odd number, show `"n/2"`, else show the exact value.
"""
show_half(x::Integer) = iseven(x) ? "$(x >> 1)" : "$x/2"

"""
    check_jm(dj::T, dm::T) where {T <: Integer}
check if the m-quantum number if one of the components of the j-quantum number,
in other words, `m` and `j` has the same parity, and `abs(m) < j`
"""
check_jm(dj::T, dm::T) where {T <: Integer} = is_same_parity(dj, dm) & (abs(dm) <= dj)

"""
if `check_jm(j, m)` failed, show this message
"""
jm_mismatching_msg(dj::T, dm::T) where {T <: Integer} = begin
    "j = $(show_half(dj)), m = $(show_half(dm)) are not matched"
end

"""
    check_couple(dj1::T, dj2::T, dj3::T) where {T <: Integer}
check if three angular monentum number `j1, j2, j3` can couple
"""
check_couple(dj1::T, dj2::T, dj3::T) where {T <: Integer} = begin
    (dj1 >= 0) & (dj2 >= 0) & is_same_parity(dj1 + dj2, dj3) & (abs(dj1 - dj2) <= dj3 <= dj1 + dj2)
end

"""
if `check_couple(j1, j2, j3)` failed, show this message
"""
miss_couple_msg(dj1::T, dj2::T, dj3::T) where {T <: Integer} = begin
    "$(show_half(dj1)), $(show_half(dj2)) cannnot couple to $(show_half(dj3))"
end


text:
module Lyaponuv

export lyaponuv_k, lyaponuv, lyaponuv_exp, lyaponuv_next, match_pairs, attractor,
    follow_points, compute_norms, column_norms, get_next

function lyaponuv_k(time_series, J, m, ref)
    X = attractor(time_series, m, J)
    norms = compute_norms(X)
    pairs = match_pairs(norms)
    y = follow_points(pairs, norms, ref)
    return(norms,y)
end

function match_pairs(norms)
    M = size(norms)[1]
    pairs = Array(Int, M)
    for row in 1:M
        mn, idx = findmin(norms[row, :])
        pairs[row] = idx
    end
    return(pairs)
end

function attractor(time_series, m, J)
    N = length(time_series)
    M = N - (m - 1) * J
    X = Array(Float64, m, M)
    i = 1
    for i=1:M
        X[:,i] = time_series[i:J:(i+(m-1)*J)]
    end
    return(X)
end

function follow_points(pairs, norms, ref)
    y = Array(Float64, ref)
    M = size(norms)[1]
    for i=0:ref-1
        agg = 0
        count = 0
        for j=1:M
            jhat = pairs[j]+i
            jtrue = j+i

            if jhat <= M && jtrue <= M
                agg = agg + log(norms[jtrue, jhat])
               # agg = agg + log(vecnorm(X[:, jtrue] - X[:, jhat]))
                count = count + 1
            end
        end
        y[i+1] = agg/count # divide by delta-t also?
    end
    return(y)
end

function compute_norms(X)
    M = size(X)[2]
    norms = Array(Float64, M, M)
    for i=1:M
        norms[i,:] = column_norms(X, i)
    end
    return(norms)
end

function column_norms(X, i)
    M = size(X)[2]
    X_diff = X .- X[:, i]
    norm_vector = [vecnorm(X_diff[:, k]) for k=1:M]
    norm_vector[i] = 10^10
    return(norm_vector)
end

function lyaponuv_exp(series)
    nn = !isnan(series)
    A = ones(length(series), 2)
    A[:,1] = linspace(1, length(series), length(series))
    gradient = \(A, series)
    return(gradient[1])
end

function lyaponuv(time_series, J, m, ref)
    ts = lyaponuv_k(time_series, J, m, ref)[2]
    exponent = lyaponuv_exp(ts[isfinite(ts)])  ## only input those which are finite
    return(exponent)
end

function get_next(ts, m, M, norms, ref, J)

    attractor_array = attractor(ts, m, J)
    temp_norms = Array(Float64, M+1, M+1)
    temp_norms[1:M, 1:M] = norms
    col = column_norms(attractor_array, M+1)
    temp_norms[M+1, :] = col
    temp_norms[:, M+1] = col

    pairs=match_pairs(temp_norms)
    lyap_k_temp = follow_points(pairs, temp_norms, ref)
    return(lyaponuv_exp(lyap_k_temp))
end

function lyaponuv_next(time_series, J, m, ref, sample_size)
    ts_diff = time_series[2:end] - time_series[1:end-1]
    sigma = std(ts_diff)
    samples = randn(sample_size) * sigma + time_series[end]
    @time norms, lyap_k = lyaponuv_k(time_series, J, m, ref)
    true_exponent = lyaponuv_exp(lyap_k)
    exponents = Array(Float64, sample_size)
    M = size(norms)[1]
    tasks = Array(Future, sample_size)

    for i=1:sample_size
        s = samples[i]
        tasks[i] = @spawn get_next(vcat(time_series, s), m, M, norms, ref, J)
        #@time exponents[i] = get_next(vcat(time_series, s), m, M, norms, ref, J)
        #@printf("process: %d\n", i)
    end

    for i=1:sample_size
        exponents[i]=fetch(tasks[i])
    end

    diff = abs(exponents .-  true_exponent)
    val, idx = findmin(diff)
    println("Next Value:", samples[idx])
    return(samples[idx])
end

end

max_stars_repo_name: guypayeur/Generic-Pred

max_issues_repo_name: guypayeur/Generic-Pred

max_forks_repo_name: guypayeur/Generic-Pred


text:
using FastHistograms, StatsBase, BenchmarkTools, DataFrames, PrettyTables

function bench_noparallel(img1, img2)
    h_fw_a_np = create_fast_histogram(
        FastHistograms.FixedWidth(),
        FastHistograms.Arithmetic(),
        FastHistograms.NoParallelization(),
        [(0x00, 0xff, 16), (0x00, 0xff, 16)],
    )
    bench_fh_fw_a_np = @benchmarkable increment_bins!($h_fw_a_np, $img1, $img2)

    edge1 = [0.0, 0.1, 0.5, 0.9, 1.0]
    edge2 = [0.0, 0.5, 0.6, 1.0]
    h_vw_bs_np = create_fast_histogram(
        FastHistograms.VariableWidth(),
        FastHistograms.BinarySearch(),
        FastHistograms.NoParallelization(),
        [edge1, edge2],
    )
    bench_fh_vw_bs_np = @benchmarkable increment_bins!($h_vw_bs_np, $img1, $img2)

    h_fw_a_simd = create_fast_histogram(
        FastHistograms.FixedWidth(),
        FastHistograms.Arithmetic(),
        FastHistograms.SIMD(),
        [(0x00, 0xff, 16), (0x00, 0xff, 16)],
    )
    bench_fh_fw_a_simd = @benchmarkable increment_bins!($h_fw_a_simd, $img1, $img2)

    img1_vec = vec(img1)
    img2_vec = vec(img2)
    bench_sb_fw = @benchmarkable StatsBase.fit(
        StatsBase.Histogram,
        ($img1_vec, $img2_vec),
        (0x00:UInt8(16):0x10f, 0x00:UInt8(16):0x10f),
    )

    bench_sb_vw = @benchmarkable StatsBase.fit(StatsBase.Histogram, ($img1_vec, $img2_vec), ($edge1, $edge2))

    return bench_fh_fw_a_np, bench_fh_vw_bs_np, bench_fh_fw_a_simd, bench_sb_fw, bench_sb_vw
end

function bench_parallel(data)
    img_vec = vec(data)
    bench_sb_fw = @benchmarkable StatsBase.fit(
        StatsBase.Histogram,
        $img_vec,
        range(Float32(0.0); stop = Float32(512.0), length = 513),
    )

    h_fw_a_pt = create_fast_histogram(
        FastHistograms.FixedWidth(),
        FastHistograms.Arithmetic(),
        FastHistograms.PrivateThreads(),
        [(Float32(0.0), Float32(512.0), 512)],
    )
    bench_fh_fw_a_pt = @benchmarkable increment_bins!($h_fw_a_pt, $data)

    return bench_fh_fw_a_pt, bench_sb_fw
end

bench_df_row(trial) = [time(trial), gctime(trial), allocs(trial), memory(trial)]

function run_all_benchmarks(img_type, img_size)
    img1 = rand(img_type, img_size...)
    img2 = rand(img_type, img_size...)
    tune!.(bench_noparallel(img1, img2))
    results_noparallel = run.(bench_noparallel(img1, img2))

    df = DataFrame(
        "Rows" => ["Min Time (ns)", "GC Time (ns)", "Allocs (B)", "Memory (B)"],
        "FastHistograms:\nFixedWidth,\nArithmetic,\nNoParallelization" => bench_df_row(results_noparallel[1]),
        "FastHistograms:\nVariableWidth,\nBinarySearch,\nNoParallelization" => bench_df_row(results_noparallel[2]),
        "FastHistograms:\nFixedWidth,\nArithmetic,\nSIMD" => bench_df_row(results_noparallel[3]),
        "StatsBase:\nFixedWidth" => bench_df_row(results_noparallel[4]),
        "StatsBase:\nVariableWidth" => bench_df_row(results_noparallel[5]),
    )
    pretty_table(df, autowrap = true, linebreaks = true)
end

function run_all_benchmarks()
    types = [UInt8, Float32]
    sizes = [(40, 80), (256, 256)]
    for img_type in types
        for img_size in sizes
            println("Benchmarking type=$img_type, size=$img_size")
            run_all_benchmarks(img_type, img_size)
        end
    end

    ahtl_data = rand(Float32, 100000000) .* (512 - Float32(0.01))
    tune!.(bench_parallel(ahtl_data))
    results_parallel = run.(bench_parallel(ahtl_data))
    df = DataFrame(
        "Rows" => ["Min Time (ns)", "GC Time (ns)", "Allocs (B)", "Memory (B)"],
        "FastHistograms:\nFixedWidth,\nArithmetic,\nPrivateThreads" => bench_df_row(results_parallel[1]), # AHTL is 108249000 ns
        "StatsBase:\nFixedWidth" => bench_df_row(results_parallel[2]),
    )
    pretty_table(df, autowrap = true, linebreaks = true)
end

run_all_benchmarks()


text:
using ColPack
using MatrixMarket
using SparseArrays
using SparseDiffTools
using LinearAlgebra

const filename = "/scratch/jac_case9241pegase_17036.mtx"

orderings = Vector{String}()

push!(orderings, "NATURAL")
push!(orderings, "LARGEST_FIRST")
push!(orderings, "DYNAMIC_LARGEST_FIRST")
push!(orderings, "DISTANCE_TWO_LARGEST_FIRST")
push!(orderings, "SMALLEST_LAST")
push!(orderings, "DISTANCE_TWO_SMALLEST_LAST")
push!(orderings, "INCIDENCE_DEGREE")
push!(orderings, "DISTANCE_TWO_INCIDENCE_DEGREE")
push!(orderings, "RANDOM")

verbose = 1
A = MatrixMarket.mmread(filename)

@time sparsediff_coloring = SparseDiffTools.matrix_colors(A)
ncolor = size(unique(sparsediff_coloring),1)
@time colpack_coloring = ColPackColoring("/scratch/jac_case9241pegase_17036.mtx", ColPack.d1, "RANDOM", verbose)

for ordering in orderings
    colpack_coloring = ColPackColoring("/scratch/jac_case9241pegase_17036.mtx", ColPack.d1, ordering, verbose)
end
colors = get_colors(colpack_coloring)
@show length(unique(colors))


text:
using DashBootstrapComponents, DashHtmlComponents

tab1_content = dbc_card(
    dbc_cardbody([
        html_p("This is tab 1!", className="card-text"),
        dbc_button("Click here", color="success"),
    ]),
    className="mt-3",
);

tab2_content = dbc_card(
    dbc_cardbody([
        html_p("This is tab 2!", className="card-text"),
        dbc_button("Don't click here", color="danger"),
    ]),
    className="mt-3",
);

tabs = dbc_tabs([
    dbc_tab(tab1_content, label="Tab 1"),
    dbc_tab(tab2_content, label="Tab 2"),
    dbc_tab("This tab's content is never seen", label="Tab 3", disabled=true),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module logger

using Logging, LoggingExtras

# Add elapsed time to logger
function ElapsedTimeLoggerDecorator(logger)
	logStartTime = time()
	return TransformerLogger(logger) do log
		t = round(time() - logStartTime, digits=2)
		merge(log, (; message = "$(t) $(log.message)"))
	end
end

export ElapsedTimeLoggerDecorator

end

max_stars_repo_name: tfburns/BurnsHagaFukai

max_issues_repo_name: tfburns/BurnsHagaFukai

max_forks_repo_name: tfburns/BurnsHagaFukai


text:
# Ragel
# =====
#
# Utilities for the Ragel state machine compiler.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Ragel

export tryread!

using BufferedStreams
import Bio.IO: FileFormat, AbstractReader

# A type keeping track of a ragel-based parser's state.
type State{T<:BufferedInputStream}
    stream::T     # input stream
    cs::Int       # current DFA state of Ragel
    linenum::Int  # line number: parser is responsible for updating this
end

function State(initstate::Int, input::BufferedInputStream)
    return State(input, initstate, 1)
end

@inline function anchor!(state, p)
    state.stream.anchor = 1 + p
end

@inline function upanchor!(state)
    @assert state.stream.anchor != 0 "upanchor! called with no anchor set"
    anchor = state.stream.anchor
    state.stream.anchor = 0
    return anchor
end

# Parse an integer from an interval in the input buffer. This differs from
# `parse(Int, ...)` in that we don't have to copy or allocate anything, and it
# doesn't check that the characters are digits (we don't need to since this is
# already checked during parsing).
function parse_int64(buffer, firstpos, lastpos)
    if buffer[firstpos] == UInt8('-')
        sign = -1
        firstpos += 1
    else
        sign = +1
    end
    x = Int64(0)
    for i in firstpos:lastpos
        x = x * 10 + buffer[i] - UInt8('0')
    end
    return sign * x
end

# Actions
# -------

# Macros that help make common parsing tasks more succinct

macro anchor!()
    quote
        anchor!($(esc(:state)), $(esc(:p)))
    end
end

macro copy_from_anchor!(dest)
    quote
        firstpos = upanchor!($(esc(:state)))
        copy!($(esc(dest)), $(esc(:state)).stream.buffer, firstpos, $(esc(:p)))
    end
end

macro append_from_anchor!(dest)
    quote
        firstpos = upanchor!($(esc(:state)))
        append!($(esc(dest)), $(esc(:state)).stream.buffer, firstpos, $(esc(:p)))
    end
end

macro int64_from_anchor!()
    quote
        firstpos = upanchor!($(esc(:state)))
        parse_int64($(esc(:state)).stream.buffer, firstpos, $(esc(:p)))
    end
end

macro float64_from_anchor!()
    quote
        firstpos = upanchor!($(esc(:state))) - 1
        get(ccall(
            :jl_try_substrtod,
            Nullable{Float64},
            (Ptr{UInt8}, Csize_t, Csize_t),
            $(esc(:state)).stream.buffer, firstpos, $(esc(:p)) - firstpos
        ))
    end
end

macro ascii_from_anchor!()
    quote
        firstpos = upanchor!($(esc(:state)))
        n = $(esc(:p)) - firstpos + 1
        dst = Vector{UInt8}(n)
        copy!(dst, 1, $(esc(:state)).stream.buffer, firstpos, n)
        String(dst)
    end
end

macro load_from_anchor!(T)
    quote
        firstpos = upanchor!($(esc(:state)))
        unsafe_load(convert(Ptr{$(T)},
            pointer($(esc(:state)).stream.buffer, firstpos)))
    end
end

# return the current character
macro char()
    quote
        convert(Char, $(esc(:state)).stream.buffer[$(esc(:p))+1])
    end
end

# advance `p`, set the current state to the target state `ts` and escape from
# scanning input
macro yield(ts)
    esc(quote
        p += 1
        cs = $ts
        @goto yield
    end)
end

# read/read! function generators
# ------------------------------

# Define a read!/read function wrapping ragel-generated parser.  These macros
# handle some the dirty work of maintaining state, refilling buffers, etc.
macro generate_read!_function(machine_name, input_type, output_type, ragel_body)
    accept_state = Symbol(machine_name, "_first_final")
    error_state = Symbol(machine_name, "_error")

    return esc(quote
        function Base.read!(input::$(input_type), output::$(output_type))
            state = input.state

            # restore state variables used by Ragel (`p` is 0-based)
            cs = state.cs
            p = state.stream.position - 1
            pe = state.stream.available
            eof = -1
            data = state.stream.buffer

            if !isopen(state.stream)
                error("input stream is already closed")
            end

            if p ≥ pe
                # fill data buffer
                if BufferedStreams.fillbuffer!(state.stream) == 0
                    throw(EOFError())
                end
                p = state.stream.position - 1
                pe = state.stream.available
            end

            # run the parser until all input is consumed or a match is found
            while true
                $(ragel_body)

                state.cs = cs
                state.stream.position = p + 1

                if cs == $(error_state)
                    error("error parsing ", $(machine_name),
                          " input on line ", state.linenum)
                elseif p == eof  # exhausted all input data
                    if cs == $(accept_state)
                        throw(EOFError())
                    else
                        error("unexpected end of input while parsing ", $(machine_name))
                    end
                elseif p == pe  # exhausted filled input data
                    # refill data buffer
                    hits_eof = BufferedStreams.fillbuffer!(state.stream) == 0
                    p = state.stream.position - 1
                    pe = state.stream.available
                    if hits_eof
                        eof = pe
                    end
                else
                    # unreachable here
                    @assert false
                end
            end

            # `@yield` jumps here
            @label yield

            # save the current state and the reading position
            state.cs = cs
            state.stream.position = p + 1

            return output
        end
    end)
end

function Base.read(input::AbstractReader)
    return read!(input, eltype(input)())
end

"""
    tryread!(reader::AbstractReader, output)

Try to read the next element into `output` from `reader`.

The result is wrapped in `Nullable` and will be null if no entry is available.
"""
function tryread!(reader::AbstractReader, output)
    T = eltype(reader)
    try
        read!(reader, output)
        return Nullable{T}(output)
    catch ex
        if isa(ex, EOFError)
            return Nullable{T}()
        end
        rethrow()
    end
end

# Iterator
# --------

function Base.start(reader::AbstractReader)
    T = eltype(reader)
    nextone = T()
    if isnull(tryread!(reader, nextone))
        return Nullable{T}()
    else
        return Nullable{T}(nextone)
    end
end

Base.done(reader::AbstractReader, nextone) = isnull(nextone)

function Base.next(reader::AbstractReader, nextone)
    item = get(nextone)
    ret = copy(item)
    if isnull(tryread!(reader, item))
        return ret, Nullable{eltype(reader)}()
    else
        return ret, Nullable(item)
    end
end

end # module Ragel


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "NL2sol"
version = v"2.2.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/macd/NL2sol.f.git", "4bdb74e725cb9996095f7b46b9e99ff7fe8f3f8e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/NL2sol.f
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make
cd ..
mkdir -p ${libdir}
mv -f libnl2sol.${dlext} ${libdir}
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnl2sol", :libnl2sol)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil