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

const CachesInpHamilton = Union{Symplectic2Cache,Symplectic3Cache,
Symplectic4Cache,Symplectic45Cache,Symplectic5Cache,
Symplectic6Cache,Symplectic62Cache,
McAte8Cache,KahanLi8Cache,SofSpa10Cache,}
const CachesInpNewton = Union{VelocityVerletCache,}

const CachesNipHamilton = Union{Symplectic2ConstantCache,Symplectic3ConstantCache,
Symplectic4ConstantCache,Symplectic45ConstantCache,Symplectic5ConstantCache,
Symplectic6ConstantCache,Symplectic62ConstantCache,
McAte8ConstantCache,KahanLi8ConstantCache,SofSpa10ConstantCache,}
const CachesNipNewton = Union{VelocityVerletConstantCache,}

# some of the algorithms are designed only for the case
# f.f2(p, q, pa, t) = p which is the Newton/Lagrange equations
# If called with different functions (which are possible in the Hamiltonian case)
# an exception is thrown to avoid silently calculate wrong results.
verify_f2(::C, f, p, q, pa, t) where C<:CachesNipHamilton = f(p, q, pa, t)
verify_f2(::C, f, res, p, q, pa, t) where C<:CachesInpHamilton = f(res, p, q, pa, t)

function verify_f2(::C, f, p, q, pa, t) where C<:CachesNipNewton
    res = f(p, q, pa, t)
    res == p ? p : throwex()
end
function verify_f2(::C, f, res, p, q, pa, t) where C<:CachesInpNewton
    f(res, p, q, pa, t)
    res == p ? res : throwex()
end
throwex() = throw(ArgumentError("This algorithm is invalid if f2(p, q, t) != p"))

function initialize!(integrator,cache::C) where C<:Union{CachesInpHamilton,CachesInpNewton}
  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)
  verify_f2(cache, 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{CachesNipHamilton,CachesNipNewton}
  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 = verify_f2(cache, integrator.f.f2, duprev, uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.fsalfirst = ArrayPartition((kdu,ku))
  integrator.k[2] = integrator.fsalfirst
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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  kdu = f.f1(du,u,p,tnew)
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 3
  integrator.destats.nf2 += 2

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  integrator.k[1] = integrator.k[2]
  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
  kuprev = integrator.k[2].x[2]
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f1(kdu,du,u,p,tnew)
  f.f2(ku,du,u,p,tnew)
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  kdu = f.f1(du,u,p,tnew)
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 4
  integrator.destats.nf2 += 3
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f1(kdu,du,u,p,tnew)
  f.f2(ku,du,u,p,tnew)
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  kdu = f.f1(du,u,p,tnew)
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 5
  integrator.destats.nf2 += 4
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f1(kdu,du,u,p,tnew)
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 5
  integrator.destats.nf2 += 4
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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)
  if typeof(integrator.alg) <: McAte42
    du = du + dt*a5*kdu
    kdu = f.f1(du,u,p,tnew)
    integrator.destats.nf += 1
  end
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 5
  integrator.destats.nf2 += 5
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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)
  if typeof(integrator.alg) <: McAte42
    @.. du = du + dt*a5*kdu
    f.f1(kdu,du,u,p,tnew)
    integrator.destats.nf += 1
  end
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 5
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  kdu = f.f1(du,u,p,tnew)
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 7
  integrator.destats.nf2 += 6
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  du,u = integrator.u.x
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f1(kdu,du,u,p,tnew)
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 7
  integrator.destats.nf2 += 6
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 8
  integrator.destats.nf2 += 8
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 8
  integrator.destats.nf2 += 8
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 10
  integrator.destats.nf2 += 10
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 10
  integrator.destats.nf2 += 10
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 16
  integrator.destats.nf2 += 16
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 16
  integrator.destats.nf2 += 16
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  ku = f.f2(du,u,p,tnew)

  integrator.destats.nf += 18
  integrator.destats.nf2 += 18
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f2(ku,du,u,p,tnew)

  integrator.destats.nf += 18
  integrator.destats.nf2 += 18
  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],ku)
  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
  kuprev = integrator.k[2].x[2]
  # update position
  u = uprev + dt*b1*kuprev
  # 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
  ku = f.f2(du,u,p,tnew)
  integrator.destats.nf += 36
  integrator.destats.nf2 += 36
  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((kdu,ku))
  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
  kuprev = integrator.k[2].x[2]
  kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  # update position
  @.. u = uprev + dt*b1*kuprev
  # 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
  f.f2(ku,du,u,p,tnew)
  integrator.destats.nf += 36
  integrator.destats.nf2 += 36
  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],ku)
  copyto!(integrator.k[2].x[1],kdu)
end


text:
module PProfAllocsTest

import PProf
import Profile
using ProtoBuf

using Test

const out = tempname()

@testset "basic profiling" begin
    Profile.Allocs.clear()
    Profile.Allocs.@profile sample_rate=1.0 begin
        # Profile compilation
        @eval foo(x,y) = x * y + y / x
        @eval foo(2, 3)
    end

    # Write the profile
    outf = PProf.Allocs.pprof(out=out, web=false)

    # Read the exported profile
    prof = open(io->readproto(io, PProf.perftools.profiles.Profile()), outf, "r")

    # Verify that we exported stack trace samples:
    @test length(prof.sample) > 0
    # Verify that we exported frame information
    @test length(prof.location) > 0
    @test length(prof._function) > 0
    @test length(prof.sample_type) >= 2  # allocs and size

end

end  # module PProfAllocsTest

ext: jl


text:
@test inf(1.) == inf(Float64) == Inf
@test inf(Float16(1.)) == inf(Float16) == Inf16
@test_throws InexactError inf(Int)

for x = 0.:0.1:10.
    @test log1mexp(x) ≈ log(1 - exp(-x))
end

for x = -100.:0.
    @test log1pexp(x) ≈ log(1 + exp(x)) atol=1e-14
end

for x = 0.:100.
    @test log1pexp(x) ≈ x + log(exp(-x) + 1) atol=1e-14
end

for x = 0.:10., y = 0.:10.
    @test logadd(log(x), log(y)) ≈ log(x + y)
    @test logsub(log(x), log(y)) ≈ log(abs(x - y)) atol=1e-14
end

for x = 0.:10., y = 0.:10., z = 0.:10.
    @test logadd(log(x), log(y), log(z)) ≈ log(x + y + z)
end

for rep = 1 : 100
    x = 10rand(10)
    lx = log.(x)
    if logadd(lx...) ≈ log(sum(x))
        @test true
        @test logadd(lx[1],lx[2],lx[3],lx[4],lx[5],lx[6],lx[7],lx[8],lx[9],lx[10]) ≈ log(sum(x))
    else
        println("logadd reduce test failed, with x=$x")
        @test false
    end
end


text:
@testset "Residue" begin
    residue = Residue("ALA", 4)
    @test name(residue) == "ALA"
    @test id(residue) == 4

    residue = Residue("GUA")
    @test name(residue) == "GUA"

    remove_chemfiles_warning() do
        @test_throws ChemfilesError id(residue)
    end

    @test size(residue) == 0
    add_atom!(residue, 0);
    add_atom!(residue, 56);
    add_atom!(residue, 30);
    @test size(residue) == 3

    add_atom!(residue, 56)
    size(residue) == 3

    @test contains(residue, 56) == true

    @test atoms(residue) == UInt64[0, 30, 56]

    topology = Topology()
    @test count_residues(topology) == 0

    remove_chemfiles_warning() do
        @test_throws ChemfilesError Residue(topology, 3)
    end

    add_residue!(topology, residue)
    residue = Residue(topology, 0)
    @test size(residue) == 3

    copy = deepcopy(residue)
    @test size(copy) == 3

    add_atom!(copy, 42)
    @test size(copy) == 4
    @test size(residue) == 3
end


text:
"""
    DOCFREE.jl
Here you are, you found an free package :-)
"""
module DocFree

export plusTwo

"""
    plusTwo(x)
Sum the numeric "2" to whatever it receives as input
# Arguments
* `x`: The amount to which we want to add 2
# Notes
* Notes can go here
# Examples
```julia
julia> five = plusTwo(3)
5
```
"""
plusTwo(x) = return x+2

end # module


text:
module test_julia_indices_linearindices

using Test

T = NTuple{N,UnitRange{Int}} where N
@test T{1}((1:2,)) == (1:2,)
@test T{2}((1:2, 3:4)) == (1:2, 3:4)

@test LinearIndices{1,Tuple{UnitRange{Int}}}((1:2,)) == LinearIndices((1:2,))
@test LinearIndices{2,Tuple{UnitRange{Int},UnitRange{Int}}}((1:2, 3:4)) == LinearIndices((1:2, 3:4))

end # module test_julia_indices

module test_julia_indices_end

using Test

A = [1,2,3]
@test A[end:-1:1] == [3,2,1]

end # module test_julia_indices_end


text:
using FileIO
using Test

struct MimeSaveTestType
end

@testset "FileIO" begin
    include("query.jl")
    include("loadsave.jl")
    include("error_handling.jl")
    include("test_mimesave.jl")
end


text:
CI = MOI.ConstraintIndex
SAF = MOI.ScalarAffineFunction

"""
    hull_reformulation!(constr::ConstraintRef{<:AbstractModel, CI{SAF{T},V}}, bin_var, args...) where {T,V}

Apply the hull reformulation to a linear constraint.

    hull_reformulation!(constr::ConstraintRef, bin_var, eps, i, j, k)

Apply the hull reformulation to a nonlinear constraint (includes quadratic) at index k of constraint j in disjunct i.

    hull_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, eps, i, j, k)

Call the hull reformulation on a constraint at index k of constraint j in disjunct i.
"""
function hull_reformulation!(constr::ConstraintRef{<:AbstractModel, CI{SAF{T},V}}, bin_var, args...) where {T,V}
    #check constraint type
    i = args[2] #get disjunct index
    bin_var_ref = constr.model[bin_var][i]
    #replace each variable with its disaggregated version
    for var_ref in constr.model[:gdp_variable_refs]
        #get disaggregated variable reference
        var_name_i = name_disaggregated_variable(var_ref, bin_var, i)
        var_i_ref = variable_by_name(constr.model, var_name_i)
        #check var_ref is present in the constraint
        coeff = normalized_coefficient(constr, var_ref)
        iszero(coeff) && continue #if not present, skip
        #swap variable for disaggregated variable
        set_normalized_coefficient(constr, var_ref, 0) #remove original variable
        set_normalized_coefficient(constr, var_i_ref, coeff) #add disaggregated variable
    end
    #multiply RHS constant by binary variable
    rhs = normalized_rhs(constr) #get rhs
    set_normalized_rhs(constr, 0) #set rhs to 0
    set_normalized_coefficient(constr, bin_var_ref, -rhs) #add binary variable (same as multiplying rhs constant by binary variable)
end
function hull_reformulation!(constr::ConstraintRef, bin_var, eps, i, j, k)
    eps = get_reform_param(eps, i, j, k)
    #create symbolic variables (using Symbolics.jl)
    sym_vars = Dict(
        symbolic_variable(var_ref) => symbolic_variable(name_disaggregated_variable(var_ref, bin_var, i))
        for var_ref in constr.model[:gdp_variable_refs]
    )
    ϵ = eps #epsilon parameter for perspective function (See Furman, Sawaya, Grossmann [2020] perspecive function)
    bin_var_sym = Symbol("$bin_var[$i]")
    λ = Num(Symbolics.Sym{Float64}(bin_var_sym))
    FSG1 = Num(Symbolics.Sym{Float64}(gensym())) #this will become: [(1-ϵ)⋅λ + ϵ] (See Furman, Sawaya, Grossmann [2020] perspecive function)
    FSG2 = Num(Symbolics.Sym{Float64}(gensym())) #this will become: ϵ⋅(1-λ) (See Furman, Sawaya, Grossmann [2020] perspecive function)

    #parse constr
    op, lhs, rhs = parse_constraint(constr)
    replace_Symvars!(lhs, constr.model) #convert JuMP variables into Symbolic variables
    gx = eval(lhs) #convert the LHS of the constraint into a Symbolic expression
    #use symbolic substitution to obtain the following expression:
    #[(1-ϵ)⋅λ + ϵ]⋅g(v/[(1-ϵ)⋅λ + ϵ]) - ϵ⋅g(0)⋅(1-λ) <= 0
    #first term
    g1 = FSG1*substitute(gx, Dict(var => var_i/FSG1 for (var,var_i) in sym_vars))
    #second term
    g0 = substitute(gx, Dict(var => 0 for var in keys(sym_vars)))
    @assert !isinf(g0.val) "Hull reformulation has failed for non-linear constraint $constr: $gx is not defined at 0. Perspective function is undetermined."
    g2 = FSG2*g0
    #create perspective function and simplify
    pers_func = simplify(g1 - g2, expand = true)
    #replace FSG expressions & simplify
    pers_func = substitute(pers_func, Dict(FSG1 => (1-ϵ)*λ+ϵ,
                                           FSG2 => ϵ*(1-λ)))
    pers_func = simplify(pers_func)
    replace_constraint(constr, pers_func, op, rhs)
end
hull_reformulation!(constr::AbstractArray{<:ConstraintRef}, bin_var, eps, i, j, k) =
    hull_reformulation!(constr[k], bin_var, eps, i, j, k)

"""
    disaggregate_variables(m::Model, disj, bin_var)

Disaggregate all variables in the model and tag them with the disjunction name.
"""
function disaggregate_variables(m::Model, disj, bin_var)
    #check that variables are bounded
    var_refs = m[:gdp_variable_refs]
    @assert all((has_upper_bound.(var_refs) .&& has_lower_bound.(var_refs)) .|| is_binary.(var_refs)) "All variables must be bounded to perform the Hull reformulation."
    #reformulate variables
    obj_dict = object_dictionary(m)
    bounds_dict = :variable_bounds_dict in keys(obj_dict) ? obj_dict[:variable_bounds_dict] : Dict() #NOTE: should pass as an keyword argument
    for var_name in m[:gdp_variable_names]
        var = m[var_name]
        #define UB and LB
        LB, UB = get_bounds(var, bounds_dict)
        #disaggregate variable and add bounding constraints
        for i in eachindex(disj)
            var_name_i_str = "$(var_name)_$bin_var$i"
            var_name_i = Symbol(var_name_i_str)
            #create disaggregated variable
            m[var_name_i] = add_disaggregated_variable(m, var, LB, UB, var_name_i_str)
            #apply bounding constraints on disaggregated variable
            var_i_lb = "$(var_name_i)_lb"
            var_i_ub = "$(var_name_i)_ub"
            m[Symbol(var_i_lb)] = @constraint(m, LB * m[bin_var][i] .- m[var_name_i] .<= 0, base_name = var_i_lb)
            m[Symbol(var_i_ub)] = @constraint(m, m[var_name_i] .- UB * m[bin_var][i] .<= 0, base_name = var_i_ub)
        end
    end
end

"""
    add_disaggregated_variable(m::Model, var::VariableRef, LB, UB, base_name)

Disaggreagate a variable with lower bound `LB`, upper bound `UB`, and name `base_name`.

    add_disaggregated_variable(m::Model, var::AbstractArray{VariableRef}, LB, UB, base_name)

Disaggregate a variable block stored in an Array or DenseAxisArray.

    add_disaggregated_variable(m::Model, var::Containers.SparseAxisArray, LB, UB, base_name)

Disaggregate a variable block stored in a SparseAxisArray.
"""
function add_disaggregated_variable(m::Model, var::VariableRef, LB, UB, base_name)
    @variable(
        m,
        lower_bound = min(LB,0),
        upper_bound = max(UB,0),
        binary = is_binary(var),
        integer = is_integer(var),
        base_name = base_name
    )
end
function add_disaggregated_variable(m::Model, var::AbstractArray{VariableRef}, LB, UB, base_name)
    idxs = Iterators.product(axes(var)...)
    var_i_array = [
        add_disaggregated_variable(m, var[idx...], LB[idx...], UB[idx...], "$base_name[$(join(idx,","))]")
        for idx in idxs
    ]
    return containerize(var, var_i_array)
end
function add_disaggregated_variable(m::Model, var::Containers.SparseAxisArray, LB, UB, base_name)
    idxs = keys(var.data)
    var_i_dict = Dict(
        idx => add_disaggregated_variable(m, var[idx], LB[idx], UB[idx], "$base_name[$(join(idx,","))]")
        for idx in idxs
    )
    return Containers.SparseAxisArray(var_i_dict)
end
containerize(var::Array, arr) = arr
containerize(var::Containers.DenseAxisArray, arr) = Containers.DenseAxisArray(arr, axes(var)...)

"""
    sum_disaggregated_variables(m::Model, disj, bin_var)

Add constraint that global variable is equal to the sum of the disaggregated copies.
"""
function sum_disaggregated_variables(m::Model, disj, bin_var)
    for var in m[:gdp_variable_refs]
        dis_vars = []
        for i in eachindex(disj)
            var_name_i = name_disaggregated_variable(var, bin_var, i)
            var_i = variable_by_name(m, var_name_i)
            push!(dis_vars, var_i)
        end
        aggr_con = "$(var)_$(bin_var)_aggregation"
        m[Symbol(aggr_con)] = @constraint(m, var == sum(dis_vars), base_name = aggr_con)
    end
end


text:
module MIMEs

export mime_from_extension, mime_from_path, extension_from_mime, charset_from_mime, compressible_from_mime, contenttype_from_mime

const _mimedb, _ext2mime, _mime2ext = include(joinpath(@__DIR__, "..", "mimedb", "mimedb.jlon"))

"""
```julia
mime_from_extension(query::AbstractString[, default::T=nothing])::Union{MIME,T}
```

# Examples:
```julia
mime_from_extension(".json") == MIME"application/json"()
mime_from_extension("html") == MIME"text/html"()
mime_from_extension("asdfff") == nothing
mime_from_extension("asdfff", MIME"text/plain"()) == MIME"text/plain"()
```
"""
function mime_from_extension(query::AbstractString, default=nothing)
    m = get(_ext2mime, lowercase(lstrip(query, '.')), nothing)
    m === nothing ? default : MIME(m)
end

"""
```julia
mime_from_path(path::AbstractString[, default::T=nothing])::Union{MIME,T}
```

Return the MIME type of the file at `path`, based on the file extension.

# Examples:
```julia
mime_from_path("hello.json") == MIME"application/json"()
mime_from_path("/home/fons/wow.html") == MIME"text/html"()
mime_from_path("/home/fons/wow") == nothing
mime_from_path("/home/fons/wow", MIME"text/plain"()) == MIME"text/plain"()
```
"""
function mime_from_path(path::AbstractString, default=nothing)
    mime_from_extension(splitext(path)[2], default)
end

"""
```julia
extension_from_mime(mime::MIME[, default::T=""])::Union{String,T}
```

Return the most common file extension used for files of the given MIME type.

# Examples:
```julia
extension_from_mime(MIME"application/json"()) == ".json"
extension_from_mime(MIME"text/html"()) == ".html"
extension_from_mime(MIME"text/blablablaa"()) == ""
extension_from_mime(MIME"text/blablablaa"(), ".bin") == ".bin"
```
"""
function extension_from_mime(mime::MIME, default="")
    exs = get(_mime2ext, string(mime), nothing)
    if exs === nothing || isempty(exs)
        default
    else
        "." * first(exs)
    end
end

"""
```julia
compressible_from_mime(mime::MIME)::Bool
```

Whether a file of the given MIME type can/should be gzipped.

# Examples:
```julia
compressible_from_mime(MIME"text/html"()) == true
compressible_from_mime(MIME"image/png"()) == false
compressible_from_mime(MIME"text/blablablaa"()) == false
```
"""
function compressible_from_mime(mime::MIME)
    c = get(_mimedb, string(mime), nothing)
    if c === nothing
        false
    else
        get(c, "compressible", false)::Bool
    end
end

"""
```julia
charset_from_mime(mime::MIME[, default::T=nothing])::Union{String,T}
```

The default charset associated with this type, if any. If not known, text MIMEs default to "UTF-8". Possible values are: $(Set([get(x, "charset", nothing) for x in (_mimedb |> values)]) |> collect |> string).

# Examples:
```julia
charset_from_mime(MIME"application/json"()) == "UTF-8"
charset_from_mime(MIME"application/x-bogus"()) == nothing
charset_from_mime(MIME"text/blablablaa"()) == "UTF-8" # because it is a `text/` mime
charset_from_mime(MIME"text/blablablaa"(), "LATIN-1") == "LATIN-1"
```
"""
function charset_from_mime(mime::MIME)
    fallback() = istextmime(mime) ? "UTF-8" : nothing
    c = get(_mimedb, string(mime), nothing)
    if c === nothing
        fallback()
    else
        get(fallback, c, "charset")::Union{String,Nothing}
    end
end

function charset_from_mime(mime::MIME, default)
    c = get(_mimedb, string(mime), nothing)
    if c === nothing
        default
    else
        get(c, "charset", default)
    end
end

"""
```julia
contenttype_from_mime(mime::MIME)::String
```

Turn a MIME into a Content-Type header value, which might include the `charset` parameter.

# Examples:
```julia
contenttype_from_mime(MIME"application/json"()) == "application/json; charset=utf-8"
contenttype_from_mime(MIME"application/x-bogus"()) == "application/x-bogus"
contenttype_from_mime(mime_from_extension(".png", MIME"application/octet-stream"())) == "image/png"
```

# See also:
[`charset_from_mime`](@ref)
"""
contenttype_from_mime(mime::MIME) = let c = charset_from_mime(mime)
    c === nothing ? string(mime) : "$(string(mime)); charset=$(lowercase(c))"
end

end

max_stars_repo_name: fonsp/MIMEs

max_issues_repo_name: fonsp/MIMEs

max_forks_repo_name: fonsp/MIMEs


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: dataexchange
using AWS.Compat
using AWS.UUIDs

"""
    CancelJob()

This operation cancels a job. Jobs can be cancelled only when they are in the WAITING state.

# Required Parameters
- `JobId`: The unique identifier for a job.

"""
cancel_job(JobId; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/jobs/$(JobId)"; aws_config=aws_config)
cancel_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/jobs/$(JobId)", args; aws_config=aws_config)

"""
    CreateDataSet()

This operation creates a data set.

# Required Parameters
- `AssetType`: The type of file your data is stored in. Currently, the supported asset type is S3_SNAPSHOT.
- `Description`: A description for the data set. This value can be up to 16,348 characters long.
- `Name`: The name of the data set.

# Optional Parameters
- `Tags`: A data set tag is an optional label that you can assign to a data set when you create it. Each tag consists of a key and an optional value, both of which you define. When you use tagging, you can also use tag-based access control in IAM policies to control access to these data sets and revisions.
"""
create_data_set(AssetType, Description, Name; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets", Dict{String, Any}("AssetType"=>AssetType, "Description"=>Description, "Name"=>Name); aws_config=aws_config)
create_data_set(AssetType, Description, Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AssetType"=>AssetType, "Description"=>Description, "Name"=>Name), args)); aws_config=aws_config)

"""
    CreateJob()

This operation creates a job.

# Required Parameters
- `Details`: The details for the CreateJob request.
- `Type`: The type of job to be created.

"""
create_job(Details, Type; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/jobs", Dict{String, Any}("Details"=>Details, "Type"=>Type); aws_config=aws_config)
create_job(Details, Type, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/jobs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Details"=>Details, "Type"=>Type), args)); aws_config=aws_config)

"""
    CreateRevision()

This operation creates a revision for a data set.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.

# Optional Parameters
- `Comment`: An optional comment about the revision.
- `Tags`: A revision tag is an optional label that you can assign to a revision when you create it. Each tag consists of a key and an optional value, both of which you define. When you use tagging, you can also use tag-based access control in IAM policies to control access to these data sets and revisions.
"""
create_revision(DataSetId; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets/$(DataSetId)/revisions"; aws_config=aws_config)
create_revision(DataSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets/$(DataSetId)/revisions", args; aws_config=aws_config)

"""
    DeleteAsset()

This operation deletes an asset.

# Required Parameters
- `AssetId`: The unique identifier for an asset.
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

"""
delete_asset(AssetId, DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)"; aws_config=aws_config)
delete_asset(AssetId, DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", args; aws_config=aws_config)

"""
    DeleteDataSet()

This operation deletes a data set.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.

"""
delete_data_set(DataSetId; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
delete_data_set(DataSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)", args; aws_config=aws_config)

"""
    DeleteRevision()

This operation deletes a revision.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

"""
delete_revision(DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
delete_revision(DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", args; aws_config=aws_config)

"""
    GetAsset()

This operation returns information about an asset.

# Required Parameters
- `AssetId`: The unique identifier for an asset.
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

"""
get_asset(AssetId, DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)"; aws_config=aws_config)
get_asset(AssetId, DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", args; aws_config=aws_config)

"""
    GetDataSet()

This operation returns information about a data set.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.

"""
get_data_set(DataSetId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
get_data_set(DataSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)", args; aws_config=aws_config)

"""
    GetJob()

This operation returns information about a job.

# Required Parameters
- `JobId`: The unique identifier for a job.

"""
get_job(JobId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs/$(JobId)"; aws_config=aws_config)
get_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs/$(JobId)", args; aws_config=aws_config)

"""
    GetRevision()

This operation returns information about a revision.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

"""
get_revision(DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
get_revision(DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", args; aws_config=aws_config)

"""
    ListDataSetRevisions()

This operation lists a data set's revisions sorted by CreatedAt in descending order.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.

# Optional Parameters
- `maxResults`: The maximum number of results returned by a single call.
- `nextToken`: The token value retrieved from a previous call to access the next page of results.
"""
list_data_set_revisions(DataSetId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions"; aws_config=aws_config)
list_data_set_revisions(DataSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions", args; aws_config=aws_config)

"""
    ListDataSets()

This operation lists your data sets. When listing by origin OWNED, results are sorted by CreatedAt in descending order. When listing by origin ENTITLED, there is no order and the maxResults parameter is ignored.

# Optional Parameters
- `maxResults`: The maximum number of results returned by a single call.
- `nextToken`: The token value retrieved from a previous call to access the next page of results.
- `origin`: A property that defines the data set as OWNED by the account (for providers) or ENTITLED to the account (for subscribers).
"""
list_data_sets(; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets"; aws_config=aws_config)
list_data_sets(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets", args; aws_config=aws_config)

"""
    ListJobs()

This operation lists your jobs sorted by CreatedAt in descending order.

# Optional Parameters
- `dataSetId`: The unique identifier for a data set.
- `maxResults`: The maximum number of results returned by a single call.
- `nextToken`: The token value retrieved from a previous call to access the next page of results.
- `revisionId`: The unique identifier for a revision.
"""
list_jobs(; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs"; aws_config=aws_config)
list_jobs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs", args; aws_config=aws_config)

"""
    ListRevisionAssets()

This operation lists a revision's assets sorted alphabetically in descending order.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

# Optional Parameters
- `maxResults`: The maximum number of results returned by a single call.
- `nextToken`: The token value retrieved from a previous call to access the next page of results.
"""
list_revision_assets(DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets"; aws_config=aws_config)
list_revision_assets(DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets", args; aws_config=aws_config)

"""
    ListTagsForResource()

This operation lists the tags on the resource.

# Required Parameters
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.

"""
list_tags_for_resource(resource_arn; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("GET", "/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    StartJob()

This operation starts a job.

# Required Parameters
- `JobId`: The unique identifier for a job.

"""
start_job(JobId; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/jobs/$(JobId)"; aws_config=aws_config)
start_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/jobs/$(JobId)", args; aws_config=aws_config)

"""
    TagResource()

This operation tags a resource.

# Required Parameters
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tags`: A label that consists of a customer-defined key and an optional value.

"""
tag_resource(resource_arn, tags; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

This operation removes one or more tags from a resource.

# Required Parameters
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tagKeys`: The key tags.

"""
untag_resource(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateAsset()

This operation updates an asset.

# Required Parameters
- `AssetId`: The unique identifier for an asset.
- `DataSetId`: The unique identifier for a data set.
- `Name`: The name of the asset. When importing from Amazon S3, the S3 object key is used as the asset name. When exporting to Amazon S3, the asset name is used as default target S3 object key.
- `RevisionId`: The unique identifier for a revision.

"""
update_asset(AssetId, DataSetId, Name, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
update_asset(AssetId, DataSetId, Name, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), args)); aws_config=aws_config)

"""
    UpdateDataSet()

This operation updates a data set.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.

# Optional Parameters
- `Description`: The description for the data set.
- `Name`: The name of the data set.
"""
update_data_set(DataSetId; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
update_data_set(DataSetId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)", args; aws_config=aws_config)

"""
    UpdateRevision()

This operation updates a revision.

# Required Parameters
- `DataSetId`: The unique identifier for a data set.
- `RevisionId`: The unique identifier for a revision.

# Optional Parameters
- `Comment`: An optional comment about the revision.
- `Finalized`: Finalizing a revision tells AWS Data Exchange that your changes to the assets in the revision are complete. After it's in this read-only state, you can publish the revision to your products.
"""
update_revision(DataSetId, RevisionId; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
update_revision(DataSetId, RevisionId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", args; aws_config=aws_config)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgApisMetaV1StatusDetails <: SwaggerModel
    causes::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1StatusCause} } # spec name: causes
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    retryAfterSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: retryAfterSeconds
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid

    function IoK8sApimachineryPkgApisMetaV1StatusDetails(;causes=nothing, group=nothing, kind=nothing, name=nothing, retryAfterSeconds=nothing, uid=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("causes"), causes)
        setfield!(o, Symbol("causes"), causes)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("retryAfterSeconds"), retryAfterSeconds)
        setfield!(o, Symbol("retryAfterSeconds"), retryAfterSeconds)
        validate_property(IoK8sApimachineryPkgApisMetaV1StatusDetails, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1StatusDetails

const _property_map_IoK8sApimachineryPkgApisMetaV1StatusDetails = Dict{Symbol,Symbol}(Symbol("causes")=>Symbol("causes"), Symbol("group")=>Symbol("group"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"), Symbol("retryAfterSeconds")=>Symbol("retryAfterSeconds"), Symbol("uid")=>Symbol("uid"))
const _property_types_IoK8sApimachineryPkgApisMetaV1StatusDetails = Dict{Symbol,String}(Symbol("causes")=>"Vector{IoK8sApimachineryPkgApisMetaV1StatusCause}", Symbol("group")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String", Symbol("retryAfterSeconds")=>"Int32", Symbol("uid")=>"String")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1StatusDetails }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1StatusDetails))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1StatusDetails }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1StatusDetails[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1StatusDetails }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1StatusDetails[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1StatusDetails)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1StatusDetails }, name::Symbol, val)
end


text:
# Author: Lukasz Bratos

include("diff.jl")

f(x) = exp(x)
g(x) = x^2.0 * sin(x)

for i in [5, 10, 15]
    Diff.rysujNnfx(f, 0.0, 1.0, i, "plot5a$i")
    Diff.rysujNnfx(g, -1.0, 1.0, i, "plot5b$i")
end


text:
using Test

using QuantumControlBase
using QuantumControlBase.TestUtils

@testset "dummy optimization" begin

    println("")
    problem = dummy_control_problem()
    result = optimize(problem; method=:dummymethod)
    @test result.converged

end


text:
using Literate, Documenter, Superfluids

Literate.markdown("src/index.jl", "./src", documenter = true)

makedocs(modules = [Superfluids], sitename = "Superfluids.jl")


text:
@inline function setDandΔs!(mechanism::Mechanism, diagonal::DiagonalEntry, body::Body)
    diagonal.D = ∂dyn∂vel(body, mechanism.Δt)
    diagonal.Δs = dynamics(mechanism, body)
    return
end

@inline function extendDandΔs!(mechanism::Mechanism, diagonal::DiagonalEntry, body::Body, ineqc::InequalityConstraint)
    diagonal.D += schurD(ineqc, body, mechanism.Δt) # + SMatrix{6,6,Float64,36}(1e-5*I)
    diagonal.Δs += schurf(mechanism, ineqc, body)
    return
end

@inline function setDandΔs!(mechanism::Mechanism, diagonal::DiagonalEntry{T,N}, eqc::EqualityConstraint) where {T,N}
    diagonal.D = @SMatrix zeros(T, N, N)
    # μ = 1e-5
    # diagonal.D = SMatrix{N,N,T,N*N}(μ*I) # TODO Positiv because of weird system? fix generally
    diagonal.Δs = g(mechanism, eqc)
    return
end

@inline function setLU!(mechanism::Mechanism, offdiagonal::OffDiagonalEntry, bodyid::Int64, eqc::EqualityConstraint)
    offdiagonal.L = -∂g∂pos(mechanism, eqc, bodyid)'
    offdiagonal.U = ∂g∂vel(mechanism, eqc, bodyid)
    return
end

@inline function setLU!(mechanism::Mechanism, offdiagonal::OffDiagonalEntry, eqc::EqualityConstraint, bodyid::Int64)
    offdiagonal.L = ∂g∂vel(mechanism, eqc, bodyid)
    offdiagonal.U = -∂g∂pos(mechanism, eqc, bodyid)'
    return
end

@inline function setLU!(offdiagonal::OffDiagonalEntry{T,N1,N2}) where {T,N1,N2}
    offdiagonal.L = @SMatrix zeros(T, N2, N1)
    offdiagonal.U = offdiagonal.L'
    return
end

@inline function updateLU1!(offdiagonal::OffDiagonalEntry, diagonal::DiagonalEntry, maintograndchild::OffDiagonalEntry, childtograndchild::OffDiagonalEntry)
    D = diagonal.D
    offdiagonal.L -= maintograndchild.L * D * childtograndchild.U
    offdiagonal.U -= childtograndchild.L * D * maintograndchild.U
    return
end

@inline function updateLU2!(offdiagonal::OffDiagonalEntry, diagonal::DiagonalEntry)
    Dinv = diagonal.Dinv
    offdiagonal.L = offdiagonal.L * Dinv
    offdiagonal.U = Dinv * offdiagonal.U
    return
end

@inline function updateD!(diagonal::DiagonalEntry, childdiagonal::DiagonalEntry, fillin::OffDiagonalEntry)
    diagonal.D -= fillin.L * childdiagonal.D * fillin.U
    return
end

function invertD!(diagonal::DiagonalEntry)
    diagonal.Dinv = inv(diagonal.D)
    return
end

@inline function LSol!(diagonal::DiagonalEntry, child::DiagonalEntry, fillin::OffDiagonalEntry)
    diagonal.Δs -= fillin.L * child.Δs
    return
end

function DSol!(diagonal::DiagonalEntry)
    diagonal.Δs = diagonal.Dinv * diagonal.Δs
    return
end

@inline function USol!(diagonal::DiagonalEntry, parent::DiagonalEntry, fillin::OffDiagonalEntry)
    diagonal.Δs -= fillin.U * parent.Δs
    return
end

function factor!(graph::Graph, ldu::SparseLDU)
    for id in graph.dfslist
        sucs = successors(graph, id)
        for childid in sucs
            offdiagonal = getentry(ldu, (id, childid))
            for grandchildid in sucs
                grandchildid == childid && break
                if hasdirectchild(graph, childid, grandchildid)
                    updateLU1!(offdiagonal, getentry(ldu, grandchildid), getentry(ldu, (id, grandchildid)), getentry(ldu, (childid, grandchildid)))
                end
            end
            updateLU2!(offdiagonal, getentry(ldu, childid))
        end

        diagonal = getentry(ldu, id)

        for childid in successors(graph, id)
            updateD!(diagonal, getentry(ldu, childid), getentry(ldu, (id, childid)))
        end
        invertD!(diagonal)
    end
end

function solve!(mechanism)
    ldu = mechanism.ldu
    graph = mechanism.graph
    dfslist = graph.dfslist

    for id in dfslist
        diagonal = getentry(ldu, id)

        for childid in successors(graph, id)
            LSol!(diagonal, getentry(ldu, childid), getentry(ldu, (id, childid)))
        end
    end

    for id in graph.rdfslist
        diagonal = getentry(ldu, id)

        DSol!(diagonal)

        for parentid in predecessors(graph, id)
            USol!(diagonal, getentry(ldu, parentid), getentry(ldu, (parentid, id)))
        end

        for childid in ineqchildren(graph, id)
            eliminatedSol!(mechanism, getineqentry(ldu, childid), diagonal, getbody(mechanism, id), getineqconstraint(mechanism, childid))
        end
    end
end

@inline function s0tos1!(component::Component)
    component.s1 = component.s0
    return
end

@inline function s1tos0!(component::Component)
    component.s0 = component.s1
    return
end

@inline function s0tos1!(ineqc::InequalityConstraint)
    ineqc.s1 = ineqc.s0
    ineqc.γ1 = ineqc.γ0
    return
end

@inline function s1tos0!(ineqc::InequalityConstraint)
    ineqc.s0 = ineqc.s1
    ineqc.γ0 = ineqc.γ1
    return
end

@inline function normΔs(component::Component)
    d = component.s1 - component.s0
    return dot(d, d)
end

@inline function normΔs(ineqc::InequalityConstraint)
    d1 = ineqc.s1 - ineqc.s0
    d2 = ineqc.γ1 - ineqc.γ0
    return dot(d1, d1) + dot(d2, d2)
end

function eliminatedSol!(mechanism::Mechanism, ineqentry::InequalityEntry, diagonal::DiagonalEntry, body::Body, ineqc::InequalityConstraint)
    Δt = mechanism.Δt
    μ = mechanism.μ
    No = 2

    φ = g(mechanism, ineqc)

    Nx = ∂g∂pos(mechanism, ineqc, body)
    Nv = ∂g∂vel(mechanism, ineqc, body)

    γ1 = ineqc.γ1
    s1 = ineqc.s1

    Δv = diagonal.Δs
    ineqentry.Δγ = γ1 ./ s1 .* φ - μ ./ s1 - γ1 ./ s1 .* (Nv * Δv)
    ineqentry.Δs = s1 .- μ ./ γ1 - s1 ./ γ1 .* ineqentry.Δγ

    return
end


text:
function powerconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where T <: PowerSystems.Storage

    pstin = m[:pstin]
    pstout = m[:pstout]
    name_index = m[:pstin].axes[1]
    time_index = m[:pstin].axes[2]

    (length(pstin.axes[2]) != time_periods) ? @error("Length of time dimension inconsistent") : true
    (length(pstout.axes[2]) != time_periods) ? @error("Length of time dimension inconsistent") : true

    pmax_in = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(pstin))), name_index, time_index)
    pmax_out= JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(pstout))), name_index, time_index)
    pmin_in = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(pstin))), name_index, time_index)
    pmin_out= JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(pstout))), name_index, time_index)

    (pstin.axes[1] !== pstout.axes[1]) ? @warn("Input/Output variables indexes are inconsistent") : true

    for t in pstin.axes[2], (ix, name) in enumerate(pstin.axes[1])
        if name == devices[ix].name
            pmin_in[name, t] = @constraint(m, pstin[name, t] <= devices[ix].inputactivepowerlimits.min)
            pmin_out[name, t] = @constraint(m, pstout[name, t] <= devices[ix].outputactivepowerlimits.min)
            pmax_in[name, t] = @constraint(m, pstin[name, t] <= devices[ix].inputactivepowerlimits.max)
            pmax_out[name, t] = @constraint(m, pstout[name, t] <= devices[ix].outputactivepowerlimits.max)
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :pmax_in, pmax_in)
    JuMP.register_object(m, :pmax_out, pmax_out)
    JuMP.register_object(m, :pmin_in, pmin_in)
    JuMP.register_object(m, :pmin_out, pmin_out)

    return m
end

function energybookkeeping(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64; ini_cond = 0.0) where T <: PowerSystems.GenericBattery

    pstin = m[:pstin]
    pstout = m[:pstout]
    ebt = m[:ebt]
    name_index = m[:ebt].axes[1]
    time_index = m[:ebt].axes[2]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent in E_bt") : true
    (pstin.axes[1] !== time_index) ? @warn("Input/Output and Battery Energy variables indexes are inconsistent") : true

    bookkeep_bt = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(ebt))), name_index, time_index)

    # TODO: Add Initial SOC for storage for sequential simulation
    for t1 = time_index[1], (ix,name) in enumerate(name_index)
        if name == devices[ix].name
            bookkeep_bt[name,t1] = @constraint(m,ebt[name,t1] == devices[ix].energy -  pstout[name,t1]/devices[ix].efficiency.out + pstin[name,t1]*devices[ix].efficiency.in)
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    for t in time_index[2:end], (ix,name) in enumerate(name_index)
        if name == devices[ix].name
            bookkeep_bt[name,t] = @constraint(m,ebt[name,t] == ebt[name,t-1] -  pstout[name,t]/devices[ix].efficiency.out + pstin[name,t]*devices[ix].efficiency.in)
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :book_keep, bookkeep_bt)

    return m

end

function energyconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where T <: PowerSystems.GenericBattery

    ebt = m[:ebt]
    name_index = m[:ebt].axes[1]
    time_index = m[:ebt].axes[2]

    (length(ebt.axes[2]) != time_periods) ? @error("Length of time dimension inconsistent") : true

    energylimit_bt = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(ebt))), name_index, time_index)

    for t in time_index, (ix,name) in enumerate(name_index)
        if name == devices[ix].name
            energylimit_bt[name,t] = @constraint(m,ebt[name,t] <= devices[ix].capacity.max)
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :energystoragelimit, energylimit_bt)

    return m
end


text:
using Printf, Distributions

struct Grid
    cells::BitArray{2}
    hwall::BitArray{2}
    vwall::BitArray{2}
end
function Grid(p::AbstractFloat, m::Integer=10, n::Integer=10)
    cells = fill(false, m, n)
    hwall = rand(Bernoulli(p), m + 1, n)
    vwall = rand(Bernoulli(p), m, n + 1)
    vwall[:, 1] = true
    vwall[:, end] = true
    return Grid(cells, hwall, vwall)
end

function Base.show(io::IO, g::Grid)
    H = (" .", " _")
    V = (":", "|")
    C = (" ", "#")
    ind = findfirst(g.cells[end, :] .& .!g.hwall[end, :])
    percolated = !iszero(ind)
    println(io, "$(size(g.cells, 1))×$(size(g.cells, 2)) $(percolated ? "Percolated" : "Not percolated") grid")
    for r in 1:size(g.cells, 1)
        println(io, "    ", join(H[w+1] for w in g.hwall[r, :]))
        println(io, " $(r % 10)) ", join(V[w+1] * C[c+1] for (w, c) in zip(g.vwall[r, :], g.cells[r, :])))
    end
    println(io, "    ", join(H[w+1] for w in g.hwall[end, :]))
    if percolated
        println(io, " !)  ", "  " ^ (ind - 1), '#')
    end
end

function floodfill!(m::Integer, n::Integer, cells::AbstractMatrix{<:Integer},
                    hwall::AbstractMatrix{<:Integer}, vwall::AbstractMatrix{<:Integer})
    # fill cells
    cells[m, n] = true
    percolated = false
    # bottom
    if m < size(cells, 1) && !hwall[m+1, n] && !cells[m+1, n]
        percolated = percolated || floodfill!(m + 1, n, cells, hwall, vwall)
    # The Bottom
    elseif m == size(cells, 1) && !hwall[m+1, n]
        return true
    end
    # left
    if n > 1 && !vwall[m, n] && !cells[m, n-1]
        percolated = percolated || floodfill!(m, n - 1, cells, hwall, vwall)
    end
    # right
    if n < size(cells, 2) && !vwall[m, n+1] && !cells[m, n+1]
        percolated = percolated || floodfill!(m, n + 1, cells, hwall, vwall)
    end
    # top
    if m > 1 && !hwall[m, n] && !cells[m-1, n]
        percolated = percolated || floodfill!(m - 1, n, cells, hwall, vwall)
    end
    return percolated
end
function pourontop!(g::Grid)
    m, n = 1, 1
    percolated = false
    while !percolated && n ≤ size(g.cells, 2)
        percolated = !g.hwall[m, n] && floodfill!(m, n, g.cells, g.hwall, g.vwall)
        n += 1
    end
    return percolated
end

function main(probs, nrep::Integer=1000)
    sampleprinted = false
    pcount = zeros(Int, size(probs))
    for (i, p) in enumerate(probs), _ in 1:nrep
        g = Grid(p)
        percolated = pourontop!(g)
        if percolated
            pcount[i] += 1
            if !sampleprinted
                println(g)
                sampleprinted = true
            end
        end
    end
    return pcount ./ nrep
end

probs = collect(10:-1:0) ./ 10
percprobs = main(probs)

println("Fraction of 1000 tries that percolate through:")
for (pr, pp) in zip(probs, percprobs)
    @printf("\tp = %.3f ⇒ freq. = %5.3f\n", pr, pp)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export predictLayer

type predictLayer <:Layer
  bottom ::Array{Layer}
  top    ::Array{Layer}
  topRange ::Vector

  rᵢ     ::AFArray{Float32}
  Predict      ::Function # predict function
  ∂Predict     ::Function # derivative of activation function
  rₒ     ::AFArray{Float32}

  ∂rᵢ_L   ::AFArray{Float32}
  ∂rₒ_L   ::AFArray{Float32}

  lock   ::Bool
end

predictLayer(inDim::Int,outDim::Int)=
  predictLayer(Layer[],
               Layer[],
               Vector{AFArray}(),

               randn(AFArray{Float32},inDim,batchsize),
               softmax,  # id,
               ∂softmax, #∂id,
               randn(AFArray{Float32},outDim,batchsize),

               randn(AFArray{Float32},inDim,batchsize),
               randn(AFArray{Float32},outDim,batchsize),

               false)