text:
# http://www.chimica.unipd.it/antonino.polimeno/pubblica/downloads/JChemPhys_101_4062.pdf

function initialize!(integrator,cache::SymplecticEulerConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  # Do the calculation pre
  # So that way FSAL interpolation
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu = integrator.f.f1(duprev,uprev,integrator.p,integrator.t)
  kuprev = integrator.f.f2(duprev,uprev,integrator.p,integrator.t)
  @muladd du = duprev + integrator.dt*kdu
  ku = integrator.f.f2(du,uprev,integrator.p,integrator.t)
  integrator.destats.nf2 += 1
  integrator.destats.nf += 2
  integrator.fsalfirst = ArrayPartition((kdu,kuprev))
  integrator.fsallast = ArrayPartition((zero(kdu),ku))
end

@muladd function perform_step!(integrator,cache::SymplecticEulerConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  kuprev = integrator.fsalfirst.x[2]
  u = uprev + dt*kuprev
  # Now actually compute the step
  # Do it at the end for interpolations!
  kdu = f.f1(duprev,u,p,t)
  du = duprev + dt*kdu
  ku = f.f2(du,u,p,t)
  integrator.destats.nf2 += 1
  integrator.destats.nf += 1

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function initialize!(integrator,cache::SymplecticEulerCache)
  integrator.kshortsize = 2
  @unpack k,fsalfirst = cache
  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  # Do the calculation pre
  # So that way FSAL interpolation
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kuprev = integrator.fsalfirst.x[2]
  kdu,ku = integrator.fsallast.x
  integrator.f.f1(kdu,duprev,uprev,integrator.p,integrator.t)
  integrator.f.f2(kuprev,duprev,uprev,integrator.p,integrator.t)
  @muladd @. du = duprev + integrator.dt*kdu
  integrator.f.f2(ku,du,uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 2
end

@muladd function perform_step!(integrator,cache::SymplecticEulerCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kuprev = integrator.fsalfirst.x[2]
  kdu,ku = integrator.fsallast.x
  @. u = uprev + dt*kuprev
  # Now actually compute the step
  # Do it at the end for interpolations!
  integrator.destats.nf2 += 1
  integrator.destats.nf += 1
  f.f1(kdu,duprev,u,p,t)
  @. du = duprev + dt*kdu
  f.f2(ku,du,u,p,t)
end

function initialize!(integrator,cache::C) where
    {C<:Union{VelocityVerletCache,Symplectic2Cache,Symplectic3Cache,Symplectic4Cache,
              Symplectic45Cache,Symplectic5Cache,Symplectic6Cache,Symplectic62Cache,
              McAte8Cache,KahanLi8Cache,SofSpa10Cache}}
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k

  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast

  duprev,uprev = integrator.uprev.x
  integrator.f.f1(integrator.k[2].x[1],duprev,uprev,integrator.p,integrator.t)
  integrator.f.f2(integrator.k[2].x[2],duprev,uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
end

function initialize!(integrator,cache::C) where
    {C<:Union{VelocityVerletConstantCache,Symplectic2ConstantCache,
              Symplectic3ConstantCache,Symplectic4ConstantCache,
              Symplectic45ConstantCache,Symplectic5ConstantCache,
              Symplectic6ConstantCache,Symplectic62ConstantCache,
              McAte8ConstantCache,KahanLi8ConstantCache,SofSpa10ConstantCache}}
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  duprev,uprev = integrator.uprev.x
  kdu  = integrator.f.f1(duprev,uprev,integrator.p,integrator.t)
  ku = integrator.f.f2(duprev,uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.fsalfirst = ArrayPartition((kdu,ku))
end

@muladd function perform_step!(integrator,cache::VelocityVerletConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  # x(t+Δt) = x(t) + v(t)*Δt + 1/2*a(t)*Δt^2
  ku = f.f1(duprev,uprev,p,t)
  dtsq = dt^2
  half = cache.half
  u = uprev + dt*duprev + dtsq*(half*ku)
  kdu = f.f1(duprev,u,p,t+dt)
  integrator.destats.nf += 2
  # v(t+Δt) = v(t) + 1/2*(a(t)+a(t+Δt))*Δt
  du = duprev + dt*(half*ku + half*kdu)

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::VelocityVerletCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # x(t+Δt) = x(t) + v(t)*Δt + 1/2*a(t)*Δt^2
  f.f1(ku,duprev,uprev,p,t)
  dtsq = dt^2
  half = cache.half
  @. u = uprev + dt*duprev + dtsq*(half*ku)
  f.f1(kdu,duprev,u,p,t+dt)
  integrator.destats.nf += 2
  # v(t+Δt) = v(t) + 1/2*(a(t)+a(t+Δt))*Δt
  @. du = duprev + dt*(half*ku + half*kdu)
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[1],kdu)
  copyto!(integrator.k[2].x[2],du)
end

@muladd function perform_step!(integrator,cache::Symplectic2ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,b1,b2 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu
  integrator.destats.nf += 2
  integrator.destats.nf2 += 1

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic2Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,b1,b2 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu
  integrator.destats.nf += 2
  integrator.destats.nf2 += 1
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic3ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,b1,b2,b3 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  integrator.destats.nf += 3
  integrator.destats.nf2 += 2
  du = du + dt*a3*kdu
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic3Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,b1,b2,b3 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu
  integrator.destats.nf += 3
  integrator.destats.nf2 += 2
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic4ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,b1,b2,b3,b4 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu
  integrator.destats.nf += 4
  integrator.destats.nf2 += 3
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic4Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,b1,b2,b3,b4 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu
  integrator.destats.nf += 4
  integrator.destats.nf2 += 3
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic45ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,b1,b2,b3,b4,b5 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  integrator.destats.nf += 5
  integrator.destats.nf2 += 4
  if typeof(integrator.alg) <: McAte42
    du = du + dt*a5*kdu
  end
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic45Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,b1,b2,b3,b4,b5 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  integrator.destats.nf += 5
  integrator.destats.nf2 += 4
  if typeof(integrator.alg) <: McAte42
    @. du = du + dt*a5*kdu
  end
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic5ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + t+a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu
  integrator.destats.nf += 6
  integrator.destats.nf2 += 5
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic5Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + t+a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu
  integrator.destats.nf += 6
  integrator.destats.nf2 += 5
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic6ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b7*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b8*ku

  kdu = f.f1(du,u,p,tnew)
  # @. du = du + dt*a8*kdu
  integrator.destats.nf += 8
  integrator.destats.nf2 += 7
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic6Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b7*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b8*ku

  f.f1(kdu,du,u,p,tnew)
  # @. du = du + dt*a8*kdu
  integrator.destats.nf += 8
  integrator.destats.nf2 += 7
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::Symplectic62ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b7*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b8*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b9*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b10*ku

  kdu = f.f1(du,u,p,tnew)
  # @. du = du + dt*a10*kdu
  integrator.destats.nf += 10
  integrator.destats.nf2 += 9
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Symplectic62Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b7*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b8*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b9*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b10*ku

  f.f1(kdu,du,u,p,tnew)
  # @. du = du + dt*a10*kdu
  integrator.destats.nf += 10
  integrator.destats.nf2 += 9
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::McAte8ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b7*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b8*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b9*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b10*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b11*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b12*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b13*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b14*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b15*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b16*ku

  kdu = f.f1(du,u,p,tnew)
  # @. du = du + dt*a16*kdu
  integrator.destats.nf += 16
  integrator.destats.nf2 += 15
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::McAte8Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b7*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b8*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b9*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b10*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b11*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b12*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b13*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b14*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b15*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b16*ku

  f.f1(kdu,du,u,p,tnew)
  # @. du = du + dt*a16*kdu
  integrator.destats.nf += 16
  integrator.destats.nf2 += 15
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::KahanLi8ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b7*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b8*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b9*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b10*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b11*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b12*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b13*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b14*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b15*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b16*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a16*kdu

  tnew = tnew + a16*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b17*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a17*kdu

  tnew = tnew + a17*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b18*ku

  kdu = f.f1(du,u,p,tnew)
  # @. du = du + dt*a18*kdu
  integrator.destats.nf += 18
  integrator.destats.nf2 += 17
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::KahanLi8Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b7*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b8*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b9*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b10*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b11*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b12*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b13*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b14*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b15*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b16*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a16*kdu

  tnew = tnew + a16*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b17*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a17*kdu

  tnew = tnew + a17*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b18*ku

  f.f1(kdu,du,u,p,tnew)
  # @. du = du + dt*a18*kdu
  integrator.destats.nf += 18
  integrator.destats.nf2 += 17
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end

@muladd function perform_step!(integrator,cache::SofSpa10ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
          a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,
          a35,a36,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,
          b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,
          b35,b36 = cache
  duprev,uprev = integrator.uprev.x
  # update position
  u = uprev + dt*b1*duprev
  # update velocity
  kdu = f.f1(duprev,u,p,integrator.t)
  du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b2*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b3*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b4*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b5*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b6*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b7*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b8*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b9*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b10*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b11*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b12*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b13*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b14*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b15*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b16*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a16*kdu

  tnew = tnew + a16*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b17*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a17*kdu

  tnew = tnew + a17*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b18*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a18*kdu

  tnew = tnew + a18*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b19*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a19*kdu

  tnew = tnew + a19*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b20*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a20*kdu

  tnew = tnew + a20*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b21*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a21*kdu

  tnew = tnew + a21*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b22*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a22*kdu

  tnew = tnew + a22*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b23*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a23*kdu

  tnew = tnew + a23*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b24*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a24*kdu

  tnew = tnew + a24*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b25*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a25*kdu

  tnew = tnew + a25*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b26*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a26*kdu

  tnew = tnew + a26*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b27*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a27*kdu

  tnew = tnew + a27*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b28*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a28*kdu

  tnew = tnew + a28*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b29*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a29*kdu

  tnew = tnew + a29*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b30*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a30*kdu

  tnew = tnew + a30*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b31*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a31*kdu

  tnew = tnew + a31*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b32*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a32*kdu

  tnew = tnew + a32*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b33*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a33*kdu

  tnew = tnew + a33*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b34*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a34*kdu

  tnew = tnew + a34*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b35*ku

  kdu = f.f1(du,u,p,tnew)
  du = du + dt*a35*kdu

  tnew = tnew + a35*dt
  ku = f.f2(du,u,p,tnew)
  u = u + dt*b36*ku

  kdu = f.f1(du,u,p,tnew)
  # @. du = du + dt*a30*kdu
  integrator.destats.nf += 36
  integrator.destats.nf2 += 35
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,du))
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::SofSpa10Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,
          a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,
          a35,a36,
          b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,
          b19,b20,b21,b22,b23,b24,b25,b26,b27,b28,b29,b30,b31,b32,b33,b34,
          b35,b36 = cache.tab
  duprev,uprev = integrator.uprev.x
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @. u = uprev + dt*b1*duprev
  # update velocity
  f.f1(kdu,duprev,u,p,integrator.t)
  @. du = duprev + dt*a1*kdu
  # update position & velocity
  tnew = t+a1*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b2*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a2*kdu

  # update position & velocity
  tnew = tnew + a2*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b3*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a3*kdu

  # update position & velocity
  tnew = tnew + a3*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b4*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a4*kdu

  # update position & velocity
  tnew = tnew + a4*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b5*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a5*kdu

  tnew = tnew + a5*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b6*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a6*kdu

  tnew = tnew + a6*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b7*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a7*kdu

  tnew = tnew + a7*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b8*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a8*kdu

  tnew = tnew + a8*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b9*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a9*kdu

  tnew = tnew + a9*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b10*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a10*kdu

  tnew = tnew + a10*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b11*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a11*kdu

  tnew = tnew + a11*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b12*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a12*kdu

  tnew = tnew + a12*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b13*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a13*kdu

  tnew = tnew + a13*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b14*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a14*kdu

  tnew = tnew + a14*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b15*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a15*kdu

  tnew = tnew + a15*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b16*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a16*kdu

  tnew = tnew + a16*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b17*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a17*kdu

  tnew = tnew + a17*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b18*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a18*kdu

  tnew = tnew + a18*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b19*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a19*kdu

  tnew = tnew + a19*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b20*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a20*kdu

  tnew = tnew + a20*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b21*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a21*kdu

  tnew = tnew + a21*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b22*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a22*kdu

  tnew = tnew + a22*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b23*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a23*kdu

  tnew = tnew + a23*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b24*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a24*kdu

  tnew = tnew + a24*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b25*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a25*kdu

  tnew = tnew + a25*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b26*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a26*kdu

  tnew = tnew + a26*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b27*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a27*kdu

  tnew = tnew + a27*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b28*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a28*kdu

  tnew = tnew + a28*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b29*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a29*kdu

  tnew = tnew + a29*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b30*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a30*kdu

  tnew = tnew + a30*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b31*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a31*kdu

  tnew = tnew + a31*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b32*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a32*kdu

  tnew = tnew + a32*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b33*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a33*kdu

  tnew = tnew + a33*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b34*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a34*kdu

  tnew = tnew + a34*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b35*ku

  f.f1(kdu,du,u,p,tnew)
  @. du = du + dt*a35*kdu

  tnew = tnew + a35*dt
  f.f2(ku,du,u,p,tnew)
  @. u = u + dt*b36*ku

  f.f1(kdu,du,u,p,tnew)
  # @. du = du + dt*a30*kdu
  integrator.destats.nf += 36
  integrator.destats.nf2 += 35
  copyto!(integrator.k[1].x[1],integrator.k[2].x[1])
  copyto!(integrator.k[1].x[2],integrator.k[2].x[2])
  copyto!(integrator.k[2].x[2],du)
  copyto!(integrator.k[2].x[1],kdu)
end


text:
  # Structured solver interface
# ========================== #
"""
    StructuredModel(stochasticprogram::StochasticProgram, solver::AbstractStructuredSolver)

Return an instance of `AbstractStructuredModel` based on `stochasticprogram` and the given `solver`.

See also: [`optimize_structured!`](@ref), [`fill_solution!`](@ref)
"""
function StructuredModel(stochasticprogram::StochasticProgram, solver::AbstractStructuredSolver)
    throw(MethodError(StructuredModel, (stochasticprogram, solver)))
end
"""
    internal_solver(solver::AbstractStructuredSolver)

Return an `AbstractMathProgSolver`, if available, from `solver`.
"""
function internal_solver(solver::AbstractStructuredSolver)
    throw(MethodError(optimsolver, solver))
end
"""
    optimize_structured!(structuredmodel::AbstractStructuredModel)

Optimize the `AbstractStructuredModel`, which also optimizes the `stochasticprogram` it was instantiated from.

See also: [`fill_solution!`](@ref)
"""
function optimize_structured!(structuredmodel::AbstractStructuredModel)
    throw(MethodError(optimize_structured!, structuredmodel))
end
"""
    fill_solution!(stochasticprogram::StochasticProgram, structuredmodel::AbstractStructuredModel)

Fill in the optimal solution in `stochasticprogram` after a call to `optimize_structured!`. Should fill in the first stage result and second stage results for each available scenario.

See also: [`optimize_structured!`](@ref)
"""
function fill_solution!(stochasticprogram::StochasticProgram, structuredmodel::AbstractStructuredModel)
    throw(MethodError(fill_solution!, stochasticprogram, structuredmodel))
end
"""
    solverstr(solver::AbstractStructuredSolver)

Optionally, return a string identifier of `AbstractStructuredSolver`.
"""
function solverstr(::AbstractStructuredSolver)
    return "Unnamed structured solver"
end
# Sample-based solver interface
# ========================== #
"""
    SampledModel(stochasticmodel::StochasticModel, sampler::AbstractSampler, solver::AbstractSampledSolver)

Return an instance of `AbstractSampledModel` based on `stochasticmodel`, `sampler` and the given `solver`.

See also: [`optimize_sampled!`](@ref), [`stochastic_solution`](@ref)
"""
function SampledModel(stochasticmodel::StochasticModel, solver::AbstractSampledSolver)
    throw(MethodError(StructuredModel, (stochasticprogram, solver)))
end
"""
    optimize_sampled!(sampledmodel::AbstractSampledmodel, sampler::AbstractSampler, confidence::AbstractFloat)

Optimize the `AbstractSampledModel` to the given `confidence` level, using `sampler` to generate scenarios. This should approximately optimize the `stochasticmodel` the sampled model was instantiated from.

See also: [`stochastic_solution`](@ref)
"""
function optimize_sampled!(sampledmodel::AbstractSampledModel, sampler::AbstractSampler, confidence::AbstractFloat)
    throw(MethodError(optimize_sampled!, sampledmodel, sampler))
end
"""
    stochastic_solution(sampledmodel::AbstractSampledmodel)

Generate a `StochasticSolution` from `sampledmodel` after a call to `optimize_sampled!`. The solution should include an approximately optimal first-stage decision, an an approximate optimal value and a confidence interval around the true optimum of the original stochastic model.

See also: [`optimize_sampled!`](@ref), [`StochasticSolution`](@ref)
"""
function stochastic_solution(sampledmodel::AbstractSampledModel)
    throw(MethodError(stochastic_solution, sampledmodel))
end
"""
    internal_solver(solver::AbstractSampledSolver)

Return an `AbstractMathProgSolver`, if available, from `solver`.
"""
function internal_solver(solver::AbstractSampledSolver)
    throw(MethodError(optimsolver, solver))
end
"""
    solverstr(solver::AbstractStructuredSolver)

Optionally, return a string identifier of `AbstractStructuredSolver`.
"""
function solverstr(::AbstractSampledSolver)
    return "Unnamed sample-based solver"
end
# ========================== #

ext: jl


text:
module MRIsim

#IMPORT PACKAGES
using LinearAlgebra: Matrix
import Base.*, Base.+, Base.-, Base./, Base.vcat, Base.size,
       Base.copy, Base.Threads.@spawn, Base.Threads.@threads,
       Base.angle, Base.abs, Base.getproperty, Base.one
using Random, LinearAlgebra, FFTW, Images, Printf, MAT, PlotlyJS,
       ProgressMeter, CUDA, Parameters, Interpolations

global γ = 42.5e6; #Hz/T gyromagnetic constant for H1

#CORE
include("Grad.jl")
include("RF.jl")
include("DAC.jl")
include("Sequence.jl")
include("Phantom.jl")
include("Magnetization.jl")
include("Recon.jl")
include("Simulator.jl")
include("Display.jl")

#UNDER DEVELOPMENT
include("DiffusionModel.jl")

#Main
export Grad, delay, dur, DAC, RF, Sequence, Phantom, Mag
#RF-related
export Rx, Ry, Rz, Q, Un
#Secondary
export PulseDesigner
# Display
export plot_seq, plot_grads_moments, plot_ksapce_trajectory
# Simulator
export simulate

#GUI
using Blink, Interact, AssetRegistry, JLD2, FileIO
!Blink.AtomShell.isinstalled() && Blink.AtomShell.install()
include("SpinLab.jl")

export SpinLab

end


text:
### mcvar and mcse stand for Monte Carlo variance and Monte Carlo error respectively

## Monte Carlo variance assuming IID samples

mcvar_iid(v::AbstractArray) = var(v)/length(v)

mcvar_iid(v::AbstractArray, region) = mapslices(mcvar_iid, v, region)

mcvar_iid(s::VariableNState{Univariate}) = mcvar_iid(s.value)

mcvar_iid(s::VariableNState{Multivariate}, i::Integer) = mcvar_iid(s.value[i, :])

mcvar_iid(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size) = eltype(s)[mcvar_iid(s, i) for i in r]

## Monte Carlo standard error assuming IID samples

mcse_iid(v::AbstractArray) = sqrt(mcvar_iid(v))

mcse_iid(v::AbstractArray, region) = mapslices(mcse_iid, v, region)

mcse_iid(s::VariableNState{Univariate}) = mcse_iid(s.value)

mcse_iid(s::VariableNState{Multivariate}, i::Integer) = mcse_iid(s.value[i, :])

mcse_iid(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size) = eltype(s)[mcse_iid(s, i) for i in r]

## Monte Carlo variance using batch means
## Reference:
## Flegal J.M and Jones G.L
## Batch Means and Spectral Variance Estimators in Markov chain Monte Carlo
## Annals of Statistics, 2010, 38 (2), pp 1034-1070

function mcvar_bm{T}(v::AbstractArray{T}; batchlen::Integer=100)
  nbatches = div(length(v), batchlen)
  @assert nbatches > 1 "Choose batch size such that the number of batches is greather than one"
  nbsamples = nbatches*batchlen
  batchmeans = T[mean(v[((j-1)*batchlen+1):(j*batchlen)]) for j = 1:nbatches]
  return batchlen*var(batchmeans)/nbsamples
end

mcvar_bm(v::AbstractArray, region; batchlen::Integer=100) = mapslices(x -> mcvar_bm(x, batchlen=batchlen), v, region)

mcvar_bm(s::VariableNState{Univariate}; batchlen::Integer=100) = mcvar_bm(s.value, batchlen=batchlen)

mcvar_bm(s::VariableNState{Multivariate}, i::Integer; batchlen::Integer=100) = mcvar_bm(s.value[i, :], batchlen=batchlen)

mcvar_bm(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; batchlen::Integer=100) =
  eltype(s)[mcvar_bm(s, i, batchlen=batchlen) for i in r]

## Monte Carlo standard error using batch means

mcse_bm(v::AbstractArray; batchlen::Integer=100) = sqrt(mcvar_bm(v, batchlen=batchlen))

mcse_bm(v::AbstractArray, region; batchlen::Integer=100) = mapslices(x -> mcse_bm(x, batchlen=batchlen), v, region)

mcse_bm(s::VariableNState{Univariate}; batchlen::Integer=100) = mcse_bm(s.value, batchlen=batchlen)

mcse_bm(s::VariableNState{Multivariate}, i::Integer; batchlen::Integer=100) = mcse_bm(s.value[i, :], batchlen=batchlen)

mcse_bm(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; batchlen::Integer=100) =
  eltype(s)[mcse_bm(s, i, batchlen=batchlen) for i in r]

## Initial monotone sequence estimator (IMSE) of Monte Carlo variance
## Reference:
## Geyer C.J
## Practical Markov Chain Monte Carlo
## Statistical Science, 1992, 7 (4), pp 473-483

function mcvar_imse{T}(v::AbstractVector{T}; maxlag::Integer=length(v)-1)
  k = convert(Integer, floor((maxlag-1)/2))
  m = k+1

  # Preallocate memory
  g = Array(T, k+1)

  # Calculate empirical autocovariance
  acv = autocov(v, 0:maxlag)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Create the monotone sequence of g
  if m > 1
    for j = 2:m
      if g[j] > g[j-1]
        g[j] = g[j-1]
      end
    end
  end

  # Calculate the initial monotone sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(v)
end

mcvar_imse(v::AbstractArray; maxlag::Integer=length(v)-1) = mcvar_imse(vec(v), maxlag=maxlag)

mcvar_imse(v::AbstractArray, region; maxlag::Integer=length(v)-1) = mapslices(x -> mcvar_imse(x, maxlag=maxlag), v, region)

mcvar_imse(v::AbstractArray, region) = mapslices(x -> mcvar_imse(x, maxlag=length(x)-1), v, region)

mcvar_imse(s::VariableNState{Univariate}; maxlag::Integer=s.n-1) = mcvar_imse(s.value, maxlag=maxlag)

mcvar_imse(s::VariableNState{Multivariate}, i::Integer; maxlag::Integer=s.n-1) = mcvar_imse(s.value[i, :], maxlag=maxlag)

mcvar_imse(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; maxlag::Integer=s.n-1) =
  eltype(s)[mcvar_imse(s, i, maxlag=maxlag) for i in r]

## Initial monotone sequence estimator (IMSE) of Monte Carlo standard error

mcse_imse(v::AbstractArray; maxlag::Integer=length(v)-1) = sqrt(mcvar_imse(v, maxlag=maxlag))

mcse_imse(v::AbstractArray, region; maxlag::Integer=length(v)-1) = mapslices(x -> mcse_imse(x, maxlag=maxlag), v, region)

mcse_imse(v::AbstractArray, region) = mapslices(x -> mcse_imse(x, maxlag=length(x)-1), v, region)

mcse_imse(s::VariableNState{Univariate}; maxlag::Integer=s.n-1) = mcse_imse(s.value, maxlag=maxlag)

mcse_imse(s::VariableNState{Multivariate}, i::Integer; maxlag::Integer=s.n-1) = mcse_imse(s.value[i, :], maxlag=maxlag)

mcse_imse(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; maxlag::Integer=s.n-1) =
  eltype(s)[mcse_imse(s, i, maxlag=maxlag) for i in r]

## Initial positive sequence estimator (IPSE) of Monte Carlo variance
## Reference:
## Geyer C.J
## Practical Markov Chain Monte Carlo
## Statistical Science, 1992, 7 (4), pp 473-483

function mcvar_ipse{T}(v::AbstractVector{T}; maxlag::Integer=length(v)-1)
  k = convert(Integer, floor((maxlag-1)/2))

  # Preallocate memory
  g = Array(T, k+1)
  m = k+1

  # Calculate empirical autocovariance
  acv = autocov(v, 0:maxlag)

  # Calculate \hat{G}_{n, m} from the autocovariances, see pp. 477 in Geyer
  for j = 0:k
    g[j+1] = acv[2*j+1]+acv[2*j+2]
    if g[j+1] <= 0
      m = j
      break
    end
  end

  # Calculate the initial positive sequence estimator
  return (-acv[1]+2*sum(g[1:m]))/length(v)
end

mcvar_ipse(v::AbstractArray; maxlag::Integer=length(v)-1) = mcvar_ipse(vec(v), maxlag=maxlag)

mcvar_ipse(v::AbstractArray, region; maxlag::Integer=length(v)-1) = mapslices(x -> mcvar_ipse(x, maxlag=maxlag), v, region)

mcvar_ipse(v::AbstractArray, region) = mapslices(x -> mcvar_ipse(x, maxlag=length(x)-1), v, region)

mcvar_ipse(s::VariableNState{Univariate}; maxlag::Integer=s.n-1) = mcvar_ipse(s.value, maxlag=maxlag)

mcvar_ipse(s::VariableNState{Multivariate}, i::Integer; maxlag::Integer=s.n-1) = mcvar_ipse(s.value[i, :], maxlag=maxlag)

mcvar_ipse(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; maxlag::Integer=s.n-1) =
  eltype(s)[mcvar_ipse(s, i, maxlag=maxlag) for i in r]

## Initial positive sequence estimator (IPSE) of Monte Carlo standard error

mcse_ipse(v::AbstractArray; maxlag::Integer=length(v)-1) = sqrt(mcvar_ipse(v, maxlag=maxlag))

mcse_ipse(v::AbstractArray, region; maxlag::Integer=length(v)-1) = mapslices(x -> mcse_ipse(x, maxlag=maxlag), v, region)

mcse_ipse(v::AbstractArray, region) = mapslices(x -> mcse_ipse(x, maxlag=length(x)-1), v, region)

mcse_ipse(s::VariableNState{Univariate}; maxlag::Integer=s.n-1) = mcse_ipse(s.value, maxlag=maxlag)

mcse_ipse(s::VariableNState{Multivariate}, i::Integer; maxlag::Integer=s.n-1) = mcse_ipse(s.value[i, :], maxlag=maxlag)

mcse_ipse(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; maxlag::Integer=s.n-1) =
  eltype(s)[mcse_ipse(s, i, maxlag=maxlag) for i in r]

## Wrapper function for estimating Monte Carlo variance using various approaches

mcvar_types = (:iid, :bm, :imse, :ipse)
mcvar_functions = Dict{Symbol, Function}(zip(mcvar_types, (mcvar_iid, mcvar_bm, mcvar_imse, mcvar_ipse)))
mcse_functions = Dict{Symbol, Function}(zip(mcvar_types, (mcse_iid, mcse_bm, mcse_imse, mcse_ipse)))

mcvar(v::AbstractArray; vtype::Symbol=:imse, args...) = mcvar_functions[vtype](v; args...)

mcvar(v::AbstractArray, region; vtype::Symbol=:imse, args...) = mcvar_functions[vtype](v, region; args...)

mcvar(s::VariableNState{Univariate}; vtype::Symbol=:imse, args...) = mcvar_functions[vtype](s; args...)

mcvar(s::VariableNState{Multivariate}, i::Integer; vtype::Symbol=:imse, args...) = mcvar_functions[vtype](s, i; args...)

mcvar(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; vtype::Symbol=:imse, args...) =
  mcvar_functions[vtype](s, r; args...)

## Wrapper function for estimating Monte Carlo standard error using various approaches

mcse(v::AbstractArray; vtype::Symbol=:imse, args...) = mcse_functions[vtype](v; args...)

mcse(v::AbstractArray, region; vtype::Symbol=:imse, args...) = mcse_functions[vtype](v, region; args...)

mcse(s::VariableNState{Univariate}; vtype::Symbol=:imse, args...) = mcse_functions[vtype](s; args...)

mcse(s::VariableNState{Multivariate}, i::Integer; vtype::Symbol=:imse, args...) = mcse_functions[vtype](s, i; args...)

mcse(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size; vtype::Symbol=:imse, args...) =
  mcse_functions[vtype](s, r; args...)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    composite(dh::DrillHole; interval=1.0, zone=nothing, mode=:equalcomp,
    mincomp=0.5, gap=0.001)

Composite a drill hole object considering the given parameters. Outputs a
new composited `DrillHole` object.

## Parameters:

* `dh`       - desurveyed drill hole
* `interval` - composite length
* `zone`     - zone column name; if considered, intervals composited together
  must have the same zone value
* `mode`     - method for compositing (see below the options available)
* `mincomp`  - minimum composite lenght; smaller intervals are discarded.
* `gap`      - two intervals are not composited together if the spacing between
  them exceeds the `gap` value

## Methods:

* `:equalcomp` - seeks to create composites with the exact `interval` length;
  borders are discarded if have length below `mincomp`. Max composite length = `interval`
* `:nodiscard` - composite lengths are defined seeking to include all possible
  intervals with length above `mincomp`. Max composite length = 1.5*`interval`
"""
function composite(dhf::DrillHole; interval::Number=1.0, zone=nothing,
  mode=:equalcomp, mincomp=0.5, gap=0.001)
  # initial assertions
  @assert mode ∈ [:equalcomp, :nodiscard] "invalid method"
  @assert mincomp <= interval "mincomp must be <= interval"

  # copy drill hole and get column names
  dh   = copy(dhf.table)
  pars = dhf.pars
  bh, from, to = pars.holeid, pars.from, pars.to

  # get group of intervals to composite within
  gps = []; c = 1
  for i in 1:size(dh,1)
    i == 1 && push!(gps,c)
    i == 1 && continue

    # check gaps, holes and zones to create/separate different groups
    t1 = (dh[i-1,to] - dh[i,from]) > gap
    t2 = dh[i-1,bh] != dh[i,bh]
    t3 = isnothing(zone) ? false : !isequal(dh[i-1,zone], dh[i,zone])
    (t1 || t2 || t3) && (c += 1)
    push!(gps,c)
  end

  # get numeric columns to composite and retain main columns
  maincols = Symbol.([bh,from,to,:LENGTH])
  zone != nothing && push!(maincols,Symbol(zone))
  num     = eltype.(eachcol(dh)) .<: Union{Missing, Number}
  numcols = setdiff(Symbol.(names(dh))[num],maincols)
  numcols = setdiff(numcols,[:X,:Y,:Z])
  select!(dh,union(maincols,numcols))

  # output table and auxiliary grouping variable
  comps = DataFrame()
  dh[!,"_GP_"] = gps

  # loop all grouped intervals
  for grp in groupby(dh,:_GP_)
    # get group length and ignore it if below mincomp; create composited group table
    len = sum(grp[!,:LENGTH])
  	len < mincomp && continue
    tab = DataFrame()

    # get composited intervals for :nodiscard compositing mode
    if mode == :nodiscard
      div = len/interval
      if div <= 1 || isinteger(div)
        nbint = ceil(Int,div)
      else
        prenbint = [floor(Int,div), ceil(Int,div)]
        closest  = argmin([abs(interval-len/x) for x in prenbint])
        nbint    = prenbint[closest]
        closest == 1 && len/nbint < mincomp && (nbint = prenbint[2])
      end

      intlen  = nbint > 1 ? len/nbint : len
      minfrom = minimum(grp[!,:FROM])
      maxto   = maximum(grp[!,:TO])

      tab[!,from]     = collect(minfrom:intlen:maxto)[1:nbint]
      tab[!,to]       = tab[!,from] .+ intlen
      tab[!,bh]      .= grp[1,bh]
      tab[!,:LENGTH] .= intlen
      zone != nothing && (tab[!,zone] .= grp[1,zone])

    # get composited intervals for :equalcomp compositing mode
    elseif mode == :equalcomp
      div   = len/interval
      nbint = floor(Int,div)
      last  = (len-nbint*interval)
      last >= mincomp && (nbint += 1)
      last <  mincomp && (last = interval)

      minfrom = minimum(grp[!,:FROM])
      maxto   = maximum(grp[!,:TO])

      tab[!,from]    = collect(minfrom:interval:maxto)[1:nbint]
      tab[!,to]      = tab[!,from] .+ interval
      tab[end,to]   -= (interval-last)
      tab[!,bh]     .= grp[1,bh]
      tab[!,:LENGTH] = tab[!,to] - tab[!,from]
      zone != nothing && (tab[!,zone] .= grp[1,zone])
    end

    # initialize numeric columns that will be composited
    for col in numcols
      tab[!,col] .= -9999.9999
    end
    allowmissing!(tab, numcols)

    # averaging numeric values for each composited interval
    wcols = union(Symbol.([from,to]), numcols)
    for i in 1:size(tab,1)
      f, t = tab[i,from], tab[i,to]
      vals = grp[(grp[!,to] .> f) .& (grp[!,from] .< t), wcols]
      size(vals,1) == 0 && (tab[i,col] = missing)
      size(vals,1) == 0 && continue
      vals[1,from] = f
      vals[end,to] = t
      leng = vals[:,to]-vals[:,from]
      for col in numcols
        ivalid = findall(!ismissing,vals[!,col])
        n = length(ivalid)
        wgt = n >= 1 ? weights(leng[ivalid]) : missing
        n == 0 && (tab[i,col] = missing)
        n == 1 && (tab[i,col] = vals[1,col])
        n >  1 && (tab[i,col] = mean(vals[ivalid,col], wgt))
      end
    end

    # add composited intervals to the output comps table
    comps = vcat(comps, tab, cols=:union)
  end

  # add coordinates to the composites and return a drill hole object
  fillxyz!(comps, dhf.trace, pars)
  DrillHole(comps, dhf.trace, pars, dhf.warns)
end


text:
using SeaPearl: BasicMetrics, AbstractMetrics, AbstractTakeObjective, TakeObjective, ValueSelection, LearnedHeuristic

include("dataframes.jl")
include("plots.jl")

export storedata
export plotNodeVisited
export plotScoreVariation
export plotRewardVariation


text:
struct EnvelopeMatrix{T} <: AbstractMatrix{Tuple{T, T}}
    H::Int
    m::Int
    n::Int
    tree::Vector{Tuple{T, T}}
end

Base.size(arg::EnvelopeMatrix) = (arg.n + 1, arg.n + 1)

function rowenvelope(A::SparseMatrixCSC{Tv, Ti}) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)
        H = cllog2(n)
        tree = undefs(Tuple{Ti, Ti}, (1 << (H + 1)) - 1)
        for j = 1 : n
            if A.colptr[j] < A.colptr[j + 1]
                tree[(1 << H) - 1 + j] = (A.rowval[A.colptr[j]], A.rowval[A.colptr[j + 1] - 1])
            else
                tree[(1 << H) - 1 + j] = (m + 1, 0)
            end
        end
        for j = (1 << H) + n : (1 << (H + 1)) - 1
            tree[j] = (m + 1, 0)
        end

        for j = (1 << (H + 1)) - 2:-2:1
            (left_lo, left_hi) = tree[j]
            (right_lo, right_hi) = tree[j + 1]
            tree[j >> 1] = (min(left_lo, right_lo), max(left_hi, right_hi))
        end
        return EnvelopeMatrix{Ti}(H, m, n, tree)
    end
end

function Base.getindex(arg::EnvelopeMatrix{T}, j, j′) where {T}
    @inbounds begin
        H = arg.H
        m = arg.m
        n = arg.n
        tree = arg.tree

        (res_lo, res_hi) = (T(m + 1), T(0))

        j = ((1 << H) - 1) + j
        j′ = ((1 << H) - 1) + (j′ - 1)

        while j <= j′
            (left_lo, left_hi) = tree[j]
            (right_lo, right_hi) = tree[j′]
            (res_lo, res_hi) = (min(res_lo, left_lo, right_lo), max(res_hi, left_hi, right_hi))
            j = (j + 1) >> 1
            j′ = (j′ - 1) >> 1
        end

        return (res_lo, res_hi)
    end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module HeadEye

using ..OpticSim
using ..OpticSim.Geometry
using ..OpticSim.Repeat
using ..OpticSim.Emitters
using ..OpticSim.Repeat

using Colors
using StaticArrays, Statistics
using LinearAlgebra
import Makie, GeometryBasics, FileIO, MeshIO

include("Misc.jl")
include("Arrangement.jl")
include("HeadEyeClasses.jl")
include("Examples.jl")

end # module HeadEye
export HeadEye


text:
using RealSense

include("rs-utils.jl")
include("gl-utils.jl")

# initialize GLFW window
glfwWidth = 1280
glfwHeight = 720
window = create_window_context(glfwWidth, glfwHeight)

# prepare drawing via OpenGL
window_frame = GLfloat[-1.0,-1.0,   1.0,-1.0,   1.0, 1.0,
                        1.0, 1.0,  -1.0, 1.0,  -1.0,-1.0]
vertex_texcoords = GLfloat[0.0,1.0,  1.0,1.0,  1.0,0.0,
                           1.0,0.0,  0.0,0.0,  0.0,1.0]
# create Vertex Buffer Objects
vbo = Ref{GLuint}(0)
glGenBuffers(1, vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
glBufferData(GL_ARRAY_BUFFER, sizeof(window_frame), window_frame, GL_STATIC_DRAW)
texcoordsVBO = Ref{GLuint}(0)
glGenBuffers(1, texcoordsVBO)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_texcoords), vertex_texcoords, GL_STATIC_DRAW)
# create Vertex Array Objects
vao = Ref{GLuint}(0)
glGenVertexArrays(1, vao)
glBindVertexArray(vao[])
glBindBuffer(GL_ARRAY_BUFFER, vao[])
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)
# create shader program
vertSource= joinpath(@__DIR__, "stream.vert")
fragSource = joinpath(@__DIR__, "stream.frag")
shaderProgramID = createprogram(vertSource, fragSource)
# generate texture handle
tex = Ref{GLuint}(0)
glGenTextures(1, tex)

# initialize RealSense camera
err = Ref{Ptr{rs2_error}}(C_NULL)
# create a context
ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)

device = get_a_realsense_device(ctx)
print_device_information(device)

sensor = get_a_sensor_from_a_device(device)

profile = choose_a_streaming_profile(sensor)
sleep(1)

function on_frame_callback(f::Ptr{rs2_frame}, u::Ptr{Void})
    rs2_enqueue_frame(f, u)
end
on_frame_ptr = cfunction(on_frame_callback, Void, (Ptr{rs2_frame},Ptr{Void}))

rs2_open(sensor, profile, err)
checkerror(err)

frame_queue = rs2_create_frame_queue(1, err)
rs2_start(sensor, on_frame_ptr, frame_queue, err)
checkerror(err)

# OpenGL rendering loop
while !GLFW.WindowShouldClose(window)
    updatefps!(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, glfwWidth, glfwHeight)
    glUseProgram(shaderProgramID)

    frame = rs2_wait_for_frame(frame_queue, 5000, err)
    checkerror(err)
    data = Ptr{UInt8}(rs2_get_frame_data(frame, err))
    checkerror(err)
    width = rs2_get_frame_width(frame, err)
    height = rs2_get_frame_height(frame, err)
    glBindTexture(GL_TEXTURE_2D, tex[])
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, data)
    mipmap()
    glBindVertexArray(vao[])
    glDrawArrays(GL_TRIANGLES, 0, 6)
    rs2_release_frame(frame)

    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end

rs2_stop(sensor, err)
checkerror(err)

rs2_close(sensor, err)
checkerror(err)

# release rs2 resources
rs2_delete_frame_queue(frame_queue)
rs2_delete_sensor(sensor)
rs2_delete_device(device)
rs2_delete_context(ctx)

# release OpenGL resources
GLFW.DestroyWindow(window)


text:
using SubstitutionTilings
using SubstitutionTilings.Penrose
import SubstitutionTilings.Penrose: ψ
using Test

using Luxor

L = Penrose.Qζ

width = 1000*4
height = 1000*3
sc = 80
@draw begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 12, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(0.1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    origin()
    setline(10)
    draw(first_tile, sc, "black", :stroke)
end width height
@draw begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 12, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(0.1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    setline(10)
    origin()
    draw(first_tile, sc, "black", :stroke)
end width height# "canonical1.png"
@png begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 12, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(0.1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 11, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(2)

    for tile in tiling
        origin()
        draw(tile, sc*1.618, "black", :stroke)
    end
    setline(10)
    draw(first_tile, sc, "black", :stroke)
end width height "canonical2.png"
@png begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 12, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(0.1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 10, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(2)

    for tile in tiling
        origin()
        draw(tile, sc*1.618^2, "black", :stroke)
    end
    setline(10)
    draw(first_tile, sc, "black", :stroke)
end width height "canonical3.png"

width = 800
height = 600
@png begin
    sc = 80
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    quadrants = Tiler(width, height, 2, 2, margin=5)

    for (pos, n) in quadrants
        if n % 2 == 0
            first_tile = hkite(0, false, L(0))
        else
            first_tile = hdart(0, false, L(0))
        end
        tiling = substitute(penrose(), [first_tile], div(n-1,2),)
        for tile in tiling
            origin()
            translate(pos)
            scale(sc)
            sethue(colors[Penrose.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
        end
    end
end width height
@png begin
    sc = 120
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    quadrants = Tiler(width, height, 2, 3, margin=5)

    for (pos, n) in quadrants
        tiling = substitute(penrose(), [first_tile], n,)
        for tile in tiling
            origin()
            translate(pos)
            scale(sc/1.618^n)
            sethue(colors[Penrose.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
        end
    end
end width height
@png begin
    sc = 120
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hdart(0, false, L(0))
    quadrants = Tiler(width, height, 2, 3, margin=5)

    for (pos, n) in quadrants
        tiling = substitute(penrose(), [first_tile], n,)
        for tile in tiling
            origin()
            translate(pos)
            scale(sc/1.618^n)
            sethue(colors[Penrose.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
        end
    end
end width height
@png begin
    sc = 80
    setline(1)
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    quadrants = Tiler(width, height, 1, 2, margin=5)

    for (pos, n) in quadrants
        if n % 2 == 0
            first_tile = hkite(0, false, L(0))
        else
            first_tile = hdart(0, false, L(0))
        end
        tiling = Penrose.force([first_tile])
        for tile in tiling
            origin()
            translate(pos)
            scale(sc)
            sethue(colors[Penrose.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
        end
        origin()
        translate(pos)
        scale(sc)
        sethue("black")
        transform(embed_aff(first_tile[1]))
        draw(first_tile[2], :stroke)
    end
end width height
@png begin
    sc = 18
    setline(1)
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    quadrants = Tiler(width, height, 1, 2, margin=80)

    for (pos, n) in quadrants
        if n % 2 == 0
            first_tile = hkite(0, false, L(0))
        else
            first_tile = hdart(0, false, L(0))
        end
        tiling = unique(substitute(forced_penrose(), ([first_tile]), 4))
        for tile in tiling
            origin()
            translate(pos)
            scale(sc)
            sethue(colors[Penrose.color(tile)])
            transform(embed_aff(tile[1]))
            draw(tile[2], :fill)
        end
        origin()
        translate(pos)
        scale(sc*1.618^4)
        sethue("black")
        transform(embed_aff(first_tile[1]))
        draw(first_tile[2], :stroke)
    end
end width height

width = 1000
height = 1000
sc = 80
@draw begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 12, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    origin()
    setline(5)
    draw(first_tile, sc, "black", :stroke)
end width height
@draw begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    tiling = substitute(penrose(), [first_tile], 11, Penrose.in_bounds, (w=width/sc, h=height/sc))
    #tiling  = pentagon
    println(typeof(tiling))
    setline(1)

    for tile in substitute(penrose(), tiling, 1)
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    origin()
    setline(5)
    draw(collect(tiling)[80], sc*1.618, "black", :stroke)
    origin()
    setline(1)
    draw(first_tile, sc, "black", :stroke)
end width height
@draw begin
    colors = ["#DD93FC", "#E7977A", "#9B70AF", "#A0644F",]
    first_tile = hkite(0, false, L(0))
    preimg = substitute(penrose(), [first_tile], 7)
    tiling = substitute(penrose(), preimg, 1)
    g = inv(collect(tiling)[49][1])
    tiling = g*tiling
    preimg = dilate(ϕ, g)*preimg
    preimg_tile = collect(preimg)[600]
    setline(1)

    for tile in tiling
        origin()
        draw(tile, sc, colors[Penrose.color(tile)], :fill)
        draw(tile, sc, "black", :stroke)
    end
    origin()
    setline(1)
    draw(first_tile, sc, "black", :stroke)
    origin()
    setline(5)
    draw(preimg_tile, sc*1.618, "black", :stroke)
end width height


text:
module MTDS

using LinearAlgebra, StatsBase, Distributions
using Flux   # must be Flux@0.90
using BSON, YAML, ArgCheck, Dates, Formatting

export create_mtlds, MTLDS_variational, create_mtpd, MTPD_variational, create_mtgru,
       MTGRU_variational

include("util.jl")
include("modelutils.jl")
using .util
using .modelutils

module model

    using ..LinearAlgebra
    using ..Flux, ..ArgCheck, ..StatsBase, ..YAML
    using ..Flux: Tracker, gate
    using ..Flux.Tracker: istracked, @grad
    using ..Distributions, Formatting

    import ..save, ..load!
    import ..modelutils: randn_repar, posterior_sample, LookupTable, PartialSplit, MultiDense, mlp,
        get_final_layer_dim, get_first_layer_dim, BRNNenc, save, load!
    using ..util
    import ..util: unsqueeze, get_strtype_wo_params, MaskedArray, stacked

    export load!, forward, encode, reconstruct, elbo, elbo_w_kl,
           MTLDS_variational, create_mtlds,
           MTPD_variational, create_mtpd,
           MTGRU_variational, create_mtgru

    # Supertype
    include("super.jl")

    # Subtypes
    include("core-mtlds.jl")
    include("subtype-lds.jl")
    include("subtype-mtpd.jl")
    include("objective.jl")
    include("core-mtrnn.jl")
    include("subtype-rnn.jl")
end

using .model

end # module


text:
using CUDA.CUDNN:
    cudnnTensorDescriptor,
    cudnnCreateTensorDescriptor,
    cudnnFilterDescriptor,
    cudnnDataType,
    cudnnDataType_t,
    CUDNN_TENSOR_NCHW,
    CUDNN_STATUS_SUCCESS,
    @retry_reclaim

@testset "cudnn/tensor" begin
    x = CUDA.rand(1,1,1,2)

    TD = cudnnTensorDescriptor
    FD = cudnnFilterDescriptor
    DT = cudnnDataType

    @test TD(x) isa TD
    @test TD(CUDNN_TENSOR_NCHW, DT(eltype(x)), Cint(ndims(x)), Cint[reverse(size(x))...]) isa TD
    td = TD(x)
    @test TD(td.ptr) isa TD
    @test Base.unsafe_convert(Ptr, TD(td.ptr)) isa Ptr

    @test FD(x) isa FD
    @test FD(DT(eltype(x)),CUDNN_TENSOR_NCHW,Cint(ndims(x)),Cint[reverse(size(x))...]) isa FD
    fd = FD(x)
    @test FD(fd.ptr) isa FD
    @test Base.unsafe_convert(Ptr, FD(fd.ptr)) isa Ptr

    @test DT(Float32) isa cudnnDataType_t

    @test (@retry_reclaim(x->(x!==CUDNN_STATUS_SUCCESS),cudnnCreateTensorDescriptor(Ref{Ptr{Cvoid}}(C_NULL)))) isa Nothing
end


text:
using ImageFiltering, OffsetArrays, ComputationalResources
using Test

import Base: *, +

struct WrappedFloat
    x::Float64
end
Base.zero(w::WrappedFloat) = WrappedFloat(0.0)
Base.zero(::Type{WrappedFloat}) = WrappedFloat(0.0)
*(w::WrappedFloat, x::Float64) = WrappedFloat(w.x*x)
+(w1::WrappedFloat, w2::WrappedFloat) = WrappedFloat(w1.x+w2.x)

@testset "1d" begin
    img = 1:8
    # Exercise all the different ways to call imfilter
    kern = centered([1/3,1/3,1/3])
    imgf = imfilter(img, kern)
    r = CPU1(Algorithm.FIR())
    @test imfilter(Float64, img, kern) == imgf
    @test imfilter(Float64, img, (kern,)) == imgf
    @test imfilter(Float64, img, kern, "replicate") == imgf
    @test imfilter(Float64, img, (kern,), "replicate") == imgf
    @test_throws ArgumentError imfilter(Float64, img, (kern,), "inner")
    @test_throws ArgumentError imfilter(Float64, img, (kern,), "nonsense")
    @test imfilter(Float64, img, kern, "replicate", Algorithm.FIR()) == imgf
    @test imfilter(Float64, img, (kern,), "replicate", Algorithm.FIR()) == imgf
    @test imfilter(r, img, kern) == imgf
    @test imfilter(r, img, kern) == imgf
    @test imfilter(r, Float64, img, kern) == imgf
    @test imfilter(r, Float64, img, (kern,)) == imgf
    @test imfilter(r, img, kern, "replicate") == imgf
    @test imfilter(r, Float64, img, kern, "replicate") == imgf
    @test imfilter(r, img, (kern,), "replicate") == imgf
    @test imfilter(r, Float64, img, (kern,), "replicate") == imgf
    @test_throws MethodError imfilter(r, img, (kern,), "replicate", Algorithm.FIR())
    @test_throws MethodError imfilter(r, Float64, img, (kern,), "replicate", Algorithm.FIR())
    out = similar(imgf)
    @test imfilter!(out, img, kern) == imgf
    @test imfilter!(out, img, (kern,)) == imgf
    @test imfilter!(out, img, (kern,), "replicate") == imgf
    @test imfilter!(out, img, (kern,), "replicate", Algorithm.FIR()) == imgf
    @test imfilter!(r, out, img, kern) == imgf
    @test imfilter!(r, out, img, (kern,)) == imgf
    @test imfilter!(r, out, img, kern, "replicate") == imgf
    @test imfilter!(r, out, img, (kern,), "replicate") == imgf
    @test_throws MethodError imfilter!(r, out, img, (kern,), "replicate", Algorithm.FIR())

    # Element-type widening (issue #17)
    v = fill(0xff, 10)
    kern = centered(fill(0xff, 3))
    @info "Two warnings are expected"
    @test_throws InexactError imfilter(v, kern)
    vout = imfilter(UInt32, v, kern)
    @test eltype(vout) == UInt32
    @test all(x->x==0x0002fa03, vout)

    # Cascades don't result in out-of-bounds values
    k1 = centered([0.25, 0.5, 0.25])
    k2 = OffsetArray([0.5, 0.5], 1:2)
    casc = imfilter(img, (k1, k2, k1))
    A0 = padarray(img, Pad(:replicate, (2,), (4,)))
    A1 = imfilter(A0, k1, Inner())
    @test A1 ≈ OffsetArray([1.0,1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75,8.0,8.0,8.0], 0:11)
    A2 = imfilter(A1, k2, Inner())
    @test A2 ≈ OffsetArray([1.625,2.5,3.5,4.5,5.5,6.5,7.375,7.875,8.0,8.0], 0:9)
    A3 = imfilter(A2, k1, Inner())
    @test casc ≈ A3
    @test size(casc) == size(img)
    # copy! kernels, presence/order doesn't matter
    kc = centered([1])
    @test ImageFiltering.iscopy(kc)
    @test imfilter(img, (k1,)) ≈ [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75]
    @test imfilter(img, (kc, k1)) ≈ [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75]
    @test imfilter(img, (k1, kc)) ≈ [1.25,2.0,3.0,4.0,5.0,6.0,7.0,7.75]

    # FFT without padding
    img = collect(1:8)
    img[1] = img[8] = 0
    out = imfilter!(CPU1(Algorithm.FFT()), similar(img, Float64), img, kern, NoPad())
    @test out ≈ imfilter(img, kern)

    # Inputs that have non-1 axes
    img = OffsetArray(zeros(11), -5:5)
    img[0] = 1
    for border in ("replicate", "circular", "symmetric", "reflect",
                   Fill(zero(eltype(img))), Inner(1), NA())
        imgf = imfilter(img, centered([0.25, 0.5, 0.25]), border)
        @test imgf[-1] == imgf[1] == 0.25
        @test imgf[0] == 0.5
        inds = axes(imgf,1)
        @test all(x->x==0, imgf[first(inds):-2]) && all(x->x==0, imgf[2:last(inds)])
    end

    # Issue #110
    img = reinterpret(WrappedFloat, rand(128))
    kern = centered(rand(31))
    imgf = imfilter(img, kern)
    @test imgf == imfilter(img, kern, ImageFiltering.FIR())
end

@testset "2d widening" begin
    # issue #17
    img = fill(typemax(Int16), 10, 10)
    kern = centered(Int16[1 2 2 2 1])
    @test_throws InexactError imfilter(img, kern)
    ret = imfilter(Int32, img, kern)
    @test eltype(ret) == Int32
    @test all(x->x==262136, ret)
end

@testset "3d" begin
    img = trues(10,10,10)
    kernel = centered(trues(3,3,3)/27)
    for border in ("replicate", "circular", "symmetric", "reflect", Fill(true))
        for alg in (Algorithm.FIR(), Algorithm.FIRTiled(), Algorithm.FFT())
            @test imfilter(img, kernel, border) ≈ img
        end
    end
    # Fill(0)
    target = convert(Array{Float64}, img)
    for i in (1,10)
        target[:,:,i] .= target[:,i,:] .= target[i,:,:] .= 2/3
    end
    for i in (1,10)
        for j in (1,10)
            target[:,i,j] .= target[i,:,j] .= target[i,j,:] .= (2/3)^2
        end
    end
    for i in (1,10)
        for j in (1,10)
            for k in (1,10)
                target[i,j,k] = (2/3)^3
            end
        end
    end
    @test imfilter(img, kernel, Fill(0)) ≈ target
    @test imfilter(img, kernel, Fill(0), Algorithm.FFT()) ≈ target
    # Inner
    @test imfilter(img, kernel, Inner()) ≈ OffsetArray(trues(8,8,8), 2:9, 2:9, 2:9)
end

nothing


text:
include("solution_file_1a.jl")
include("solution_file_2a.jl")
include("solution_file_1b.jl")
include("solution_file_2b.jl")


text:
module TestBifurcations
using Test

RUNINFO = []

TEST_GROUPS = Dict{String, Vector{String}}(
    # ENV["CI_GROUP"] => test files
    "0" => [
        "test_utils.jl",
        "test_normal_form.jl",
        "test_smoke.jl",
        "test_jacobian.jl",
        "test_residual_jacobian.jl",
        "test_calcium.jl",
        "test_predator_prey.jl",
        "test_bazykin_85.jl",
        "test_bautin.jl",
        "test_examples.jl",
        "test_vs_svector.jl",
        "test_duffing_van_der_pol.jl",
        "test_hopf_to_lc.jl",
    ],
    "1" => [
        "test_reparametrization.jl",
        "test_reparametrized_bautin.jl",
    ],
    "2" => [
        "test_fold_lc.jl",
        "test_morris_lecar.jl",
        "test_diffeq_jacobian.jl",
    ],
)
TEST_GROUPS["all"] = vcat(last.(sort(collect(TEST_GROUPS), by=first))...)

@testset "$file" for file in TEST_GROUPS[get(ENV, "CI_GROUP", "all")]
    _, time, bytes, gctime, memallocs =  @timed @time include(file)
    push!(RUNINFO, Dict(
        :file => file,
        :time => time,
        :bytes => bytes,
        :gctime => gctime,
        :memallocs => memallocs,
    ))
end

function print_runinfo(io = stdout, runinfo = RUNINFO)
    io = IOContext(io, :compact => true)  # reduce number of digits shown
    titles = ["File", "Time (sec)", "%GC", "GiB"]
    table = Array{Any}(undef, (length(runinfo) + 1, length(titles)))
    table[1, :] .= titles
    for (i, rec) in enumerate(runinfo)
        table[i + 1, 1] = rec[:file]
        table[i + 1, 2] = rec[:time]
        table[i + 1, 3] = rec[:gctime] / rec[:time] * 100
        table[i + 1, 4] = rec[:bytes] / 2^30
    end
    Base.print_matrix(io, table)
    println(io)
end

try
    println()
    print_runinfo()
    println()
catch err
    @warn "Error while printing RUNINFO" exception=err
end

end  # module


text:
using MyExamplePackage
using Test

@testset "MyExamplePackage.jl" begin
    # 2x + 3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end


text:
module mesh_examples

using FinEtools
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionCircle, CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_solid, space_aspectratio
using PlotlyJS
using JSON

function curve_mesh()
    L = 42
    xyz = [0 0 0;
    0 L/4 L*1/4;
    L/4 L/4 L*2/4;
    L/4 0 L*3/4;
    0 0 L]
    nL = 20

    cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
    fens, fes = frame_member(xyz, nL, cs)
    plots = cat(plot_nodes(fens),
        plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4),
        dims = 1)
    # push!(plots, plot_nodes(fens))
    render(plots; aspectratio = [1.0 1.0 4.0])
    true
end # curve_mesh

function line_mesh_solid()
    L = 40.2
    xyz = [0 0 0;
    0 0 L]
    nL = 1

    cs = CrossSectionCircle(s -> 2.5, s -> [0.0, 1.0, 0.0])
    fens, fes = frame_member(xyz, nL, cs)
    plots = cat(plot_nodes(fens),
        plot_solid(fens, fes);
        dims = 1)
    # push!(plots, plot_nodes(fens))
    render(plots; aspectratio = [1.0 1.0 1.0])
    true
end # curve_mesh

function curve_mesh_solid()
    L = 42
    xyz = [0 0 0;
    0 L/4 L*1/4;
    L/4 L/4 L*2/4;
    L/4 0 L*3/4;
    0 0 L]
    nL = 20

    cs = CrossSectionCircle(s -> 0.5, s -> [0.0, 0.0, 1.0])
    fens, fes = frame_member(xyz, nL, cs)
    plots = cat(plot_nodes(fens),
        plot_solid(fens, fes);
        dims = 1)
    # push!(plots, plot_nodes(fens))
    render(plots; aspectratio = space_aspectratio(fens.xyz))
    true
end # curve_mesh

function argyr_l_frame()
    # Parameters:
    E=71240.0;#MPa
    nu=0.31;# Poisson ratio
    rho=5e-9;
    b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

    ##
    # Choose the mass formulation:
    mass_type=2;

    ##
    # Reference frequencies
    reffs = [11.2732, 30.5269]
    neigvs = 2;

    # Select the number of elements per half the length.
    xyz =
    n=8;
    members = []
    push!(members, frame_member([0 0 L; L 0 L], n, cs))
    push!(members, frame_member([L 0 L; L 0 0], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    plots = cat(plot_space_box([[0 -L/2 0]; [L L/2 L]]),
        plot_nodes(fens),
        plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4), dims = 1)
    render(plots; aspectratio = [1.0 1.0 1.0])
    true
end # argyr_l_frame

function argyr_l_frame_movable()
    # Parameters:
    E=71240.0;#MPa
    nu=0.31;# Poisson ratio
    rho=5e-9;
    b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

    ##
    # Choose the mass formulation:
    mass_type=2;

    ##
    # Reference frequencies
    reffs = [11.2732, 30.5269]
    neigvs = 2;

    # Select the number of elements per half the length.
    xyz =
    n=1;
    members = []
    push!(members, frame_member([0 0 L; L 0 L], n, cs))
    push!(members, frame_member([L 0 L; L 0 0], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    plots = cat(plot_space_box([[0 -L/2 0]; [L L/2 L]]),
        plot_nodes(fens),
        plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4), dims = 1)
    p = render(plots; aspectratio = [1.0 1.0 1.0])
    @show p.plot.data
    true
end # argyr_l_frame

function argyr_l_frame_solid()
    # Parameters:
    E=71240.0;#MPa
    nu=0.31;# Poisson ratio
    rho=5e-9;
    b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

    ##
    # Choose the mass formulation:
    mass_type=2;

    ##
    # Reference frequencies
    reffs = [11.2732, 30.5269]
    neigvs = 2;

    # Select the number of elements per half the length.
    xyz =
    n=8;
    members = []
    push!(members, frame_member([0 0 L; L 0 L], n, cs))
    push!(members, frame_member([L 0 L; L 0 0], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    plots = cat(plot_space_box([[0 -L/2 0]; [L L/2 L]]),
        plot_nodes(fens),
        plot_solid(fens, fes; facecolor = "rgb(155, 0, 0)"), dims = 1)
    render(plots; aspectratio = [1.0 1.0 1.0])
    true
end # argyr_l_frame

function curve_mesh_change_view()
    L = 42
    xyz = [0 0 0;
    0 L/4 L*1/4;
    L/4 L/4 L*2/4;
    L/4 0 L*3/4;
    0 0 L]
    nL = 20

    cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
    fens, fes = frame_member(xyz, nL, cs)
    plots = cat(plot_nodes(fens),
        plot_midline(fens, fes; color = "rgb(155, 155, 255)", lwidth = 4),
        dims = 1)
    # push!(plots, plot_nodes(fens))
    pl = render(plots; aspectratio = [1.0 1.0 4.0])
    savejson(pl, "plot.json")
    true
end # curve_mesh

function allrun()
    println("#####################################################")
    println("# curve_mesh ")
    curve_mesh()
    println("#####################################################")
    println("# curve_mesh_solid ")
    curve_mesh_solid()
    println("#####################################################")
    println("# argyr_l_frame ")
    argyr_l_frame()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
@with_kw mutable struct TextStyle
    font ::Option{String}  = ∅ #
    hei  ::Option{Float64} = ∅ #
    color::Option{Color}   = ∅ #
end

@with_kw mutable struct LineStyle
    lstyle::Option{Int}     = ∅ #
    lwidth::Option{Float64} = ∅ #
    smooth::Option{Bool}    = ∅ #
    color ::Option{Color}   = ∅ #
end

@with_kw mutable struct MarkerStyle
    marker::Option{String} = ∅ #
    msize::Option{Float64} = ∅ #
    color::Option{Color}   = ∅ #
end

@with_kw mutable struct BarStyle
    color::Option{Color}    = ∅ #
    fill::Colorant          = c"white" #
#   pattern::Option{String}   =  .... see page 39 of manual, test first
end

@with_kw mutable struct FillStyle
    fill::Colorant = c"cornflowerblue"
end

@with_kw mutable struct BoxplotStyle
    # box and whisker styling
    bwidth::Float64   = 0.6 # width of the box
    wwidth::Float64   = 0.3 # width of the whiskers
    wrlength::Float64 = 1.5 # whisker length is wrlength * IQR, if INF will be min-max
    blstyle::LineStyle = LineStyle(lstyle=1, lwidth=0, color=c"black")
    # median
    mlstyle::LineStyle = LineStyle(lstyle=1, lwidth=0, color=c"seagreen")
    # mean
    mshow::Bool          = true
    mmstyle::MarkerStyle = MarkerStyle(marker="fdiamond", msize=.4, color=c"dodgerblue")
    # outliers
    oshow::Bool          = true
    omstyle::MarkerStyle = MarkerStyle(marker="+", msize=.5, color=c"tomato")
end

"""
    str(m::MarkerStyle)

Internal function to help in the specific case where a line with markers of different
color than the line is required. In that case a subroutine has to be written to help
GLE, see [`add_sub_marker!`](@ref).
"""
str(m::MarkerStyle) = "mk_$(m.marker)_$(col2str(m.color; str=true))"


text:
using GeoArrays
using ArchGDAL

infile = "N:/DATA/3h metrology data/Data_forcing_01dy_010deg/nc/lrad_ITPCAS-CMFD_V0106_B-01_01dy_010deg_198301-198312.nc"
r = GeoArrays.read(infile)

dataset = ArchGDAL.unsafe_read(infile)
A = ArchGDAL.read(dataset)


text:
@testset "Local Level With Cycle Model" begin
    rj_temp = CSV.read(StateSpaceModels.RJ_TEMPERATURE, DataFrame).Values

    @test has_fit_methods(LocalLevelCycle)
    model = LocalLevelCycle(rj_temp)
    fit!(model)

    # forecasting
    forec = forecast(model, 52)
    @test monotone_forecast_variance(forec)

    smoother = kalman_smoother(model)
    alpha = get_smoothed_state(smoother)
    @test maximum(alpha[:, 1]) >= 296
    @test minimum(alpha[:, 1]) <= 296
    @test maximum(alpha[:, 2]) >= 7.5
    @test minimum(alpha[:, 2]) <= -7.5
end