text:
function initialize!(integrator, cache::Vern6ConstantCache)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 9) : (integrator.kshortsize = 12)
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:8
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast

  if !alg.lazy
    @inbounds for i in 10:12
      integrator.k[i] = zero(integrator.fsalfirst)
    end
  end
end

@muladd function perform_step!(integrator, cache::Vern6ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9 = cache.tab
  k1 = integrator.fsalfirst
  a = dt*a21
  k2 = f(uprev+a*k1, p, t + c1*dt)
  k3 = f(uprev+dt*(a31*k1+a32*k2), p, t + c2*dt)
  k4 = f(uprev+dt*(a41*k1       +a43*k3), p, t + c3*dt)
  k5 = f(uprev+dt*(a51*k1       +a53*k3+a54*k4), p, t + c4*dt)
  k6 = f(uprev+dt*(a61*k1       +a63*k3+a64*k4+a65*k5), p, t + c5*dt)
  k7 = f(uprev+dt*(a71*k1       +a73*k3+a74*k4+a75*k5+a76*k6), p, t + c6*dt)
  g8 =   uprev+dt*(a81*k1       +a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  k8 = f(g8, p, t+dt)
  u = uprev+dt*(a91*k1              +a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
  integrator.fsallast = f(u, p, t+dt); k9 = integrator.fsallast
  integrator.destats.nf += 8
  if typeof(integrator.alg) <: CompositeAlgorithm
    g9 = u
    ϱu = integrator.opts.internalnorm(k9 - k8,t)
    ϱd = integrator.opts.internalnorm(g9 - g8,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2;
  integrator.k[3]=k3; integrator.k[4]=k4;
  integrator.k[5]=k5; integrator.k[6]=k6;
  integrator.k[7]=k7; integrator.k[8]=k8;
  integrator.k[9]=k9

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab.extra
    k[10] = f(uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9]),p,t+c10*dt)
    k[11] = f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10]),p,t+c11*dt)
    k[12] = f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11]),p,t+c12*dt)
    integrator.destats.nf += 3
  end

  integrator.u = u
end

function initialize!(integrator, cache::Vern6Cache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 9) : (integrator.kshortsize = 12)
  integrator.fsalfirst = cache.k1 ; integrator.fsallast = cache.k9
  @unpack k = integrator
  resize!(k, integrator.kshortsize)
  k[1]=cache.k1; k[2]=cache.k2; k[3]=cache.k3;
  k[4]=cache.k4; k[5]=cache.k5; k[6]=cache.k6;
  k[7]=cache.k7; k[8]=cache.k8; k[9]=cache.k9 # Set the pointers

  if !alg.lazy
    k[10] = similar(cache.k1)
    k[11] = similar(cache.k1)
    k[12] = similar(cache.k1)
  end

  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::Vern6Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,utilde,tmp,atmp = cache
  a = dt*a21
  @.. tmp = uprev+a*k1
  f(k2, tmp, p, t + c1*dt)
  @.. tmp = uprev+dt*(a31*k1+a32*k2)
  f(k3, tmp, p, t + c2*dt)
  @.. tmp = uprev+dt*(a41*k1+a43*k3)
  f(k4, tmp, p, t + c3*dt)
  @.. tmp = uprev+dt*(a51*k1+a53*k3+a54*k4)
  f(k5, tmp, p, t + c4*dt)
  @.. tmp = uprev+dt*(a61*k1+a63*k3+a64*k4+a65*k5)
  f(k6, tmp, p, t + c5*dt)
  @.. tmp = uprev+dt*(a71*k1+a73*k3+a74*k4+a75*k5+a76*k6)
  f(k7, tmp, p, t + c6*dt)
  @.. tmp = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
  f(k8, tmp, p, t+dt)
  @.. u = uprev+dt*(a91*k1+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
  f(k9, u, p, t+dt)
  integrator.destats.nf += 8
  if integrator.alg isa CompositeAlgorithm
    g9 = u
    g8 = tmp
    @.. utilde = k9 - k8
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g9 - g8
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab.extra
    @unpack tmp = cache
    @.. tmp = uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9])
    f(k[10],tmp,p,t+c10*dt)
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10])
    f(k[11],tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11])
    integrator.destats.nf += 3
    f(k[12],tmp,p,t+c12*dt)
  end
  return nothing
end

function initialize!(integrator, cache::Vern7ConstantCache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 10) : (integrator.kshortsize = 16)
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  @inbounds for i in eachindex(integrator.k)
    integrator.k[i] = zero(integrator.uprev)./oneunit(integrator.t)
  end
end

@muladd function perform_step!(integrator, cache::Vern7ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,k,f,p = integrator
  @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107,b1,b4,b5,b6,b7,b8,b9,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9,btilde10 = cache.tab
  k1 = f(uprev, p, t)
  a = dt*a021
  k2 = f(uprev+a*k1, p, t + c2*dt)
  k3 = f(uprev+dt*(a031*k1+a032*k2), p, t + c3*dt)
  k4 = f(uprev+dt*(a041*k1       +a043*k3), p, t + c4*dt)
  k5 = f(uprev+dt*(a051*k1       +a053*k3+a054*k4), p, t + c5*dt)
  k6 = f(uprev+dt*(a061*k1       +a063*k3+a064*k4+a065*k5), p, t + c6*dt)
  k7 = f(uprev+dt*(a071*k1       +a073*k3+a074*k4+a075*k5+a076*k6), p, t + c7*dt)
  k8 = f(uprev+dt*(a081*k1       +a083*k3+a084*k4+a085*k5+a086*k6+a087*k7), p, t + c8*dt)
  g9 =   uprev+dt*(a091*k1          +a093*k3+a094*k4+a095*k5+a096*k6+a097*k7+a098*k8)
  g10=   uprev+dt*(a101*k1          +a103*k3+a104*k4+a105*k5+a106*k6+a107*k7)
  k9 = f(g9, p, t+dt)
  k10= f(g10, p, t+dt)
  integrator.destats.nf += 10
  u = uprev + dt*(b1*k1 + b4*k4 + b5*k5 + b6*k6 + b7*k7 + b8*k8 + b9*k9)
  if typeof(integrator.alg) <: CompositeAlgorithm
    ϱu = integrator.opts.internalnorm(k10 - k9,t)
    ϱd = integrator.opts.internalnorm(g10 - g9,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9 + btilde10*k10)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2;
  integrator.k[3]=k3; integrator.k[4]=k4;
  integrator.k[5]=k5; integrator.k[6]=k6;
  integrator.k[7]=k7; integrator.k[8]=k8;
  integrator.k[9]=k9; integrator.k[10]=k10
  integrator.u = u

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab.extra
    k[11] = f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]),p,t+c11*dt)
    k[12] = f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11]),p,t+c12*dt)
    k[13] = f(uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12]),p,t+c13*dt)
    k[14] = f(uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13]),p,t+c14*dt)
    k[15] = f(uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13]),p,t+c15*dt)
    k[16] = f(uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13]),p,t+c16*dt)
    integrator.destats.nf += 6
  end

end

function initialize!(integrator, cache::Vern7Cache)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10 = cache
  @unpack k = integrator
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 10) : (integrator.kshortsize = 16)
  resize!(k, integrator.kshortsize)
  k[1]=k1;k[2]=k2;k[3]=k3;k[4]=k4;k[5]=k5;k[6]=k6;k[7]=k7;k[8]=k8;k[9]=k9;k[10]=k10 # Setup pointers

  if !alg.lazy
    k[11] = similar(cache.k1)
    k[12] = similar(cache.k1)
    k[13] = similar(cache.k1)
    k[14] = similar(cache.k1)
    k[15] = similar(cache.k1)
    k[16] = similar(cache.k1)
  end
end

@muladd function perform_step!(integrator, cache::Vern7Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107,b1,b4,b5,b6,b7,b8,b9,btilde1,btilde4,btilde5,btilde6,btilde7,btilde8,btilde9,btilde10 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,tmp,atmp = cache
  f(k1, uprev, p, t)
  a = dt*a021
  @.. tmp = uprev+a*k1
  f(k2, tmp, p, t + c2*dt)
  @.. tmp = uprev+dt*(a031*k1+a032*k2)
  f(k3, tmp, p, t + c3*dt)
  @.. tmp = uprev+dt*(a041*k1+a043*k3)
  f(k4, tmp, p, t + c4*dt)
  @.. tmp = uprev+dt*(a051*k1+a053*k3+a054*k4)
  f(k5, tmp, p, t + c5*dt)
  @.. tmp = uprev+dt*(a061*k1+a063*k3+a064*k4+a065*k5)
  f(k6, tmp, p, t + c6*dt)
  @.. tmp = uprev+dt*(a071*k1+a073*k3+a074*k4+a075*k5+a076*k6)
  f(k7, tmp, p, t + c7*dt)
  @.. tmp = uprev+dt*(a081*k1+a083*k3+a084*k4+a085*k5+a086*k6+a087*k7)
  f(k8, tmp, p, t + c8*dt)
  @.. tmp = uprev+dt*(a091*k1+a093*k3+a094*k4+a095*k5+a096*k6+a097*k7+a098*k8)
  f(k9, tmp, p, t+dt)
  @.. tmp = uprev+dt*(a101*k1+a103*k3+a104*k4+a105*k5+a106*k6+a107*k7)
  f(k10, tmp, p, t+dt)
  @.. u = uprev + dt*(b1*k1 + b4*k4 + b5*k5 + b6*k6 + b7*k7 + b8*k8 + b9*k9)
  integrator.destats.nf += 10
  if integrator.alg isa CompositeAlgorithm
    g10 = u
    g9 = tmp
    @.. utilde = k10 - k9
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g10 - g9
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9 + btilde10*k10)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack tmp = cache
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab.extra
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9])
    f(k[11],tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11])
    f(k[12],tmp,p,t+c12*dt)
    @.. tmp = uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12])
    f(k[13],tmp,p,t+c13*dt)
    @.. tmp = uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13])
    f(k[14],tmp,p,t+c14*dt)
    tmp=  uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13])
    f(k[15],tmp,p,t+c15*dt)
    @.. tmp = uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13])
    f(k[16],tmp,p,t+c16*dt)
    integrator.destats.nf += 6
  end
  return nothing
end

function initialize!(integrator, cache::Vern8ConstantCache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 13) : (integrator.kshortsize = 21)
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  @inbounds for i in eachindex(integrator.k)
    integrator.k[i] = zero(integrator.uprev)./oneunit(integrator.t)
  end
end

@muladd function perform_step!(integrator, cache::Vern8ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13 = cache.tab
  k1 = f(uprev, p, t)
  a = dt*a0201
  k2 = f(uprev+a*k1, p, t + c2*dt)
  k3 = f(uprev+dt*(a0301*k1+a0302*k2), p, t + c3*dt)
  k4 = f(uprev+dt*(a0401*k1       +a0403*k3), p, t + c4*dt)
  k5 = f(uprev+dt*(a0501*k1       +a0503*k3+a0504*k4), p, t + c5*dt)
  k6 = f(uprev+dt*(a0601*k1                +a0604*k4+a0605*k5), p, t + c6*dt)
  k7 = f(uprev+dt*(a0701*k1                +a0704*k4+a0705*k5+a0706*k6), p, t + c7*dt)
  k8 = f(uprev+dt*(a0801*k1                +a0804*k4+a0805*k5+a0806*k6+a0807*k7), p, t + c8*dt)
  k9 = f(uprev+dt*(a0901*k1                +a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8), p, t + c9*dt)
  k10= f(uprev+dt*(a1001*k1                +a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9), p, t + c10*dt)
  k11= f(uprev+dt*(a1101*k1                +a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10), p, t + c11*dt)
  g12=  uprev+dt*(a1201*k1                +a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
  g13=  uprev+dt*(a1301*k1                +a1304*k4+a1305*k5+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10)
  k12= f(g12, p, t+dt)
  k13= f(g13, p, t+dt)
  integrator.destats.nf += 13
  u = uprev + dt*(b1*k1 + b6*k6 + b7*k7 + b8*k8 + b9*k9 + b10*k10 + b11*k11 + b12*k12)
  if typeof(integrator.alg) <: CompositeAlgorithm
    ϱu = integrator.opts.internalnorm(k13 - k12,t)
    ϱd = integrator.opts.internalnorm(g13 - g12,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9 + btilde10*k10 + btilde11*k11 + btilde12*k12 + btilde13*k13)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1]=k1; integrator.k[2]=k2;
  integrator.k[3]=k3; integrator.k[4]=k4;
  integrator.k[5]=k5; integrator.k[6]=k6;
  integrator.k[7]=k7; integrator.k[8]=k8;
  integrator.k[9]=k9; integrator.k[10]=k10;
  integrator.k[11]=k11; integrator.k[12]=k12;
  integrator.k[13]=k13
  integrator.u = u

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab.extra
    k[14] = f(uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12]),p,t+c14*dt)
    k[15] = f(uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14]),p,t+c15*dt)
    k[16] = f(uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15]),p,t+c16*dt)
    k[17] = f(uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16]),p,t+c17*dt)
    k[18] = f(uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17]),p,t+c18*dt)
    k[19] = f(uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17]),p,t+c19*dt)
    k[20] = f(uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17]),p,t+c20*dt)
    k[21] = f(uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17]),p,t+c21*dt)
    integrator.destats.nf += 8
  end

end

function initialize!(integrator, cache::Vern8Cache)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13 = cache
  @unpack k = integrator
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 13) : (integrator.kshortsize = 21)
  resize!(k, integrator.kshortsize)
  k[1]=k1;k[2]=k2;k[3]=k3;k[4]=k4;k[5]=k5;k[6]=k6;k[7]=k7;k[8]=k8;k[9]=k9;k[10]=k10;k[11]=k11;k[12]=k12;k[13]=k13 # Setup pointers

  if !alg.lazy
    for i in 14:21
      k[i] = similar(cache.k1)
    end
  end
end

@muladd function perform_step!(integrator, cache::Vern8Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310,b1,b6,b7,b8,b9,b10,b11,b12,btilde1,btilde6,btilde7,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,tmp,atmp = cache
  f(k1, uprev, p, t)
  a = dt*a0201
  @.. tmp = uprev+a*k1
  f(k2, tmp, p, t + c2*dt)
  @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
  f(k3, tmp, p, t + c3*dt)
  @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
  f(k4, tmp, p, t + c4*dt)
  @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
  f(k5, tmp, p, t + c5*dt)
  @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
  f(k6, tmp, p, t + c6*dt)
  @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
  f(k7, tmp, p, t + c7*dt)
  @.. tmp = uprev+dt*(a0801*k1+a0804*k4+a0805*k5+a0806*k6+a0807*k7)
  f(k8, tmp, p, t + c8*dt)
  @.. tmp = uprev+dt*(a0901*k1+a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)
  f(k9, tmp, p, t + c9*dt)
  @.. tmp = uprev+dt*(a1001*k1+a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
  f(k10, tmp, p, t + c10*dt)
  @.. tmp = uprev+dt*(a1101*k1+a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
  f(k11, tmp, p, t + c11*dt)
  @.. tmp = uprev+dt*(a1201*k1+a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
  f(k12, tmp, p, t+dt)
  @.. u = uprev+dt*(a1301*k1+a1304*k4+a1305*k5+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10)
  f(k13, u, p, t+dt)
  integrator.destats.nf += 13
  if integrator.alg isa CompositeAlgorithm
    g13 = u
    g12 = tmp
    @.. utilde = k13 - k12
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g13 - g12
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  @.. u = uprev + dt*(b1*k1 + b6*k6 + b7*k7 + b8*k8 + b9*k9 + b10*k10 + b11*k11 + b12*k12)
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde6*k6 + btilde7*k7 + btilde8*k8 + btilde9*k9 + btilde10*k10 + btilde11*k11 + btilde12*k12 + btilde13*k13)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab.extra
    @unpack tmp = cache
    @.. tmp = uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12])
    f(k[14],tmp,p,t+c14*dt)
    @.. tmp = uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14])
    f(k[15],tmp,p,t+c15*dt)
    @.. tmp = uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15])
    f(k[16],tmp,p,t+c16*dt)
    @.. tmp = uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16])
    f(k[17],tmp,p,t+c17*dt)
    @.. tmp = uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17])
    f(k[18],tmp,p,t+c18*dt)
    @.. tmp = uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17])
    f(k[19],tmp,p,t+c19*dt)
    @.. tmp = uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17])
    f(k[20],tmp,p,t+c20*dt)
    @.. tmp = uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17])
    integrator.destats.nf += 8
    f(k[21],tmp,p,t+c21*dt)
  end
  return nothing
end

function initialize!(integrator, cache::Vern9ConstantCache)
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 10) : (integrator.kshortsize = 20)
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

  # Avoid undefined entries if k is an array of arrays
  @inbounds for i in eachindex(integrator.k)
    integrator.k[i] = zero(integrator.uprev)./oneunit(integrator.t)
  end
end

@muladd function perform_step!(integrator, cache::Vern9ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613,b1,b8,b9,b10,b11,b12,b13,b14,b15,btilde1,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13,btilde14,btilde15,btilde16 = cache.tab
  k1 = f(uprev, p, t)
  a = dt*a0201
  k2 = f(uprev+a*k1, p, t + c1*dt)
  k3 = f(uprev+dt*(a0301*k1+a0302*k2), p, t + c2*dt)
  k4 = f(uprev+dt*(a0401*k1       +a0403*k3), p, t + c3*dt)
  k5 = f(uprev+dt*(a0501*k1       +a0503*k3+a0504*k4), p, t + c4*dt)
  k6 = f(uprev+dt*(a0601*k1                +a0604*k4+a0605*k5), p, t + c5*dt)
  k7 = f(uprev+dt*(a0701*k1                +a0704*k4+a0705*k5+a0706*k6), p, t + c6*dt)
  k8 = f(uprev+dt*(a0801*k1                                  +a0806*k6+a0807*k7), p, t + c7*dt)
  k9 = f(uprev+dt*(a0901*k1                                  +a0906*k6+a0907*k7+a0908*k8), p, t + c8*dt)
  k10 =f(uprev+dt*(a1001*k1                                  +a1006*k6+a1007*k7+a1008*k8+a1009*k9), p, t + c9*dt)
  k11= f(uprev+dt*(a1101*k1                                  +a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10), p, t + c10*dt)
  k12= f(uprev+dt*(a1201*k1                                  +a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11), p, t + c11*dt)
  k13= f(uprev+dt*(a1301*k1                                  +a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10+a1311*k11+a1312*k12), p, t + c12*dt)
  k14= f(uprev+dt*(a1401*k1                                  +a1406*k6+a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13), p, t + c13*dt)
  g15=   uprev+dt*(a1501*k1                                  +a1506*k6+a1507*k7+a1508*k8+a1509*k9+a1510*k10+a1511*k11+a1512*k12+a1513*k13+a1514*k14)
  g16=   uprev+dt*(a1601*k1                                  +a1606*k6+a1607*k7+a1608*k8+a1609*k9+a1610*k10+a1611*k11+a1612*k12+a1613*k13)
  k15= f(g15, p, t+dt)
  k16= f(g16, p, t+dt)
  integrator.destats.nf += 16
  u = uprev + dt*(b1*k1+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12+b13*k13+b14*k14+b15*k15)
  if typeof(integrator.alg) <: CompositeAlgorithm
    ϱu = integrator.opts.internalnorm(k16 - k15,t)
    ϱd = integrator.opts.internalnorm(g16 - g15,t)
    integrator.eigen_est = ϱu/ϱd
  end
  if integrator.opts.adaptive
    utilde = dt*(btilde1*k1 + btilde8*k8 + btilde9*k9 + btilde10*k10 + btilde11*k11 + btilde12*k12 + btilde13*k13 + btilde14*k14 + btilde15*k15 + btilde16*k16)
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  # k2, k3,k4,k5,k6,k7 are not used in the code (not even in interpolations), we dont need their pointers.
  # So we mapped k[2] (from integrator) with k8 (from cache), k[3] with k9 and so on.
  integrator.k[1]=k1; integrator.k[2]=k8;
  integrator.k[3]=k9; integrator.k[4]=k10;
  integrator.k[5]=k11; integrator.k[6]=k12;
  integrator.k[7]=k13; integrator.k[8]=k14;
  integrator.k[9]=k15; integrator.k[10]=k16
  integrator.u = u

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab.extra
    k[11] = f(uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9]),p,t+c17*dt)
    k[12] = f(uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11]),p,t+c18*dt)
    k[13] = f(uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12]),p,t+c19*dt)
    k[14] = f(uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13]),p,t+c20*dt)
    k[15] = f(uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14]),p,t+c21*dt)
    k[16] = f(uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15]),p,t+c22*dt)
    k[17] = f(uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15]),p,t+c23*dt)
    k[18] = f(uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15]),p,t+c24*dt)
    k[19] = f(uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15]),p,t+c25*dt)
    k[20] = f(uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15]),p,t+c26*dt)
    integrator.destats.nf += 10
  end
end

function initialize!(integrator, cache::Vern9Cache)
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16 = cache
  @unpack k = integrator
  alg = unwrap_alg(integrator, false)
  alg.lazy ? (integrator.kshortsize = 10) : (integrator.kshortsize = 20)
  resize!(k, integrator.kshortsize)
  # k2, k3,k4,k5,k6,k7 are not used in the code (not even in interpolations), we dont need their pointers.
  # So we mapped k[2] (from integrator) with k8 (from cache), k[3] with k9 and so on.
  k[1]=k1;k[2]=k8;k[3]=k9;k[4]=k10;k[5]=k11;k[6]=k12;k[7]=k13;k[8]=k14;k[9]=k15;k[10]=k16 # Setup pointers

  if !alg.lazy
    for i in 11:20
      k[i] = similar(cache.k1)
    end
  end
end

@muladd function perform_step!(integrator, cache::Vern9Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  uidx = eachindex(integrator.uprev)
  @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613,b1,b8,b9,b10,b11,b12,b13,b14,b15,btilde1,btilde8,btilde9,btilde10,btilde11,btilde12,btilde13,btilde14,btilde15,btilde16 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,utilde,tmp,atmp = cache
  f(k1, uprev, p, t)
  a = dt*a0201
  @.. tmp = uprev+a*k1
  f(k2, tmp, p, t + c1*dt)
  @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
  f(k3, tmp, p, t + c2*dt)
  @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
  f(k4, tmp, p, t + c3*dt)
  @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
  f(k5, tmp, p, t + c4*dt)
  @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
  f(k6, tmp, p, t + c5*dt)
  @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
  f(k7, tmp, p, t + c6*dt)
  @.. tmp = uprev+dt*(a0801*k1+a0806*k6+a0807*k7)
  f(k8, tmp, p, t + c7*dt)
  @.. tmp = uprev+dt*(a0901*k1+a0906*k6+a0907*k7+a0908*k8)
  f(k9, tmp, p, t + c8*dt)
  @.. tmp = uprev+dt*(a1001*k1+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
  f(k10, tmp, p, t + c9*dt)
  @.. tmp = uprev+dt*(a1101*k1+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
  f(k11, tmp, p, t + c10*dt)
  @.. tmp = uprev+dt*(a1201*k1+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
  f(k12, tmp, p, t + c11*dt)
  @.. tmp = uprev+dt*(a1301*k1+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10+a1311*k11+a1312*k12)
  f(k13, tmp, p, t + c12*dt)
  @.. tmp = uprev+dt*(a1401*k1+a1406*k6+a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13)
  f(k14, tmp, p, t + c13*dt)
  @.. tmp = uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8+a1509*k9+a1510*k10+a1511*k11+a1512*k12+a1513*k13+a1514*k14)
  f(k15, tmp, p, t+dt)
  @.. u = uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9+a1610*k10+a1611*k11+a1612*k12+a1613*k13)
  f(k16, u, p, t+dt)
  integrator.destats.nf += 16
  if integrator.alg isa CompositeAlgorithm
    g16 = u
    g15 = tmp
    @.. utilde = k16 - k15
    ϱu = integrator.opts.internalnorm(utilde,t)
    @.. utilde = g16 - g15
    ϱd = integrator.opts.internalnorm(utilde,t)
    integrator.eigen_est = ϱu/ϱd
  end
  @.. u = uprev + dt*(b1*k1+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12+b13*k13+b14*k14+b15*k15)
  if integrator.opts.adaptive
    @.. utilde = dt*(btilde1*k1 + btilde8*k8 + btilde9*k9 + btilde10*k10 + btilde11*k11 + btilde12*k12 + btilde13*k13 + btilde14*k14 + btilde15*k15 + btilde16*k16)
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  alg = unwrap_alg(integrator, false)
  if !alg.lazy && (integrator.opts.adaptive == false || integrator.EEst <= 1.0)
    k = integrator.k
    @unpack tmp = cache
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab.extra
    @.. tmp = uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9])
    f(k[11],tmp,p,t+c17*dt)
    @.. tmp = uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11])
    f(k[12],tmp,p,t+c18*dt)
    @.. tmp = uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12])
    f(k[13],tmp,p,t+c19*dt)
    @.. tmp = uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13])
    f(k[14],tmp,p,t+c20*dt)
    @.. tmp = uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14])
    f(k[15],tmp,p,t+c21*dt)
    @.. tmp = uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15])
    f(k[16],tmp,p,t+c22*dt)
    @.. tmp = uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15])
    f(k[17],tmp,p,t+c23*dt)
    @.. tmp  = uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15])
    f(k[18],tmp,p,t+c24*dt)
    @.. tmp = uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15])
    f(k[19],tmp,p,t+c25*dt)
    @.. tmp = uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15])
    integrator.destats.nf += 10
    f(k[20],tmp,p,t+c26*dt)
  end
  return nothing
end


text:
const urlApi = "https://api.eveonline.com/"

api_to_filename(url::AbstractString) = replace(rsplit(url, '/'; limit = 2, keep = false)[2], ".aspx", "")
api_cache_name(url::AbstractString) = joinpath(cacheDir, api_to_filename(url))

const timeFormat = "y-m-d H:M:S"

function xml_read(fileName::AbstractString)
	local xmlStr, xdoc
	try
		xmlStr = readstring(fileName)
	catch e
		if isa(e, SystemError)
			return nothing
		else
			rethrow()
		end
	end
	try
		xdoc = LightXML.parse_string(xmlStr)
	catch e
		if isa(e, LightXML.XMLParseError)
			return nothing
		else
			rethrow()
		end
	end
	xroot = LightXML.root(xdoc)
	return xroot
end

function xml_write(fileName::AbstractString, data)
	dirName = dirname(fileName)
	if !isdir(dirName)
		mkpath(dirName)
	end
	open(fileName, "w+") do s
		write(s, data)
	end
end

function check_time(fileName::AbstractString, xroot::LightXML.XMLElement)
	serverTime = LightXML.content(LightXML.find_element(xroot, "currentTime"))
	serverCached = LightXML.content(LightXML.find_element(xroot, "cachedUntil"))
	serverTimeF = Dates.datetime2unix(Dates.DateTime(serverTime, timeFormat))
	serverCachedF = Dates.datetime2unix(Dates.DateTime(serverCached, timeFormat))
	writeTimeF = mtime(fileName)
	validTimeF = writeTimeF + (serverCachedF - serverTimeF)
	if validTimeF <= time()
		return nothing
	end
	return xroot
end

xml_read_timed(fileName::AbstractString) = check_time(fileName, xml_read_timed(fileName))

function xml_join_attribute_values(xElem::LightXML.XMLElement, keyAttribs::Array)
	values = map(k->LightXML.attribute(xElem, k), keyAttribs)
	if any(v->v==nothing, values)
		return nothing
	end
	return values
end

function xml_find_by_attribute(xElem::LightXML.XMLElement, key::Array, keyAttribs::Array)
	matches = LightXML.XMLElement[]
	for e in LightXML.child_elements(xElem)
		attribVals = xml_join_attribute_values(e, keyAttribs)
		if attribVals == key
			push!(matches, e)
		end
	end
	return matches
end

function add_dataframe_row!(df::DataFrames.DataFrame, keyAttribs::Array, values::Array)
	attrSyms = map(Symbol, keyAttribs)
	newCols = setdiff(attrSyms, names(df))
	for c in newCols
		df[c] = DataFrames.DataArray(AbstractString, size(df, 1))
	end
	row = DataFrames.DataArray(AbstractString, size(df, 2))
	attrIndices = indexin(attrSyms, names(df))
	for i in 1:length(attrIndices)
		row[attrIndices[i]] = values[i]
	end
	push!(df, row)
end

function xml_find(xElem::LightXML.XMLElement, keyArr::Array, keyIndex::Int, results::DataFrames.DataFrame)
	keyAttr = LightXML.attribute(xElem, "key")
	keyAttribs = keyAttr != nothing? split(keyAttr, ","; keep = false) : ["name"]
	key = keyArr[keyIndex]
	attrKey = split(key, ":")
	matchAttribs = length(attrKey) > 1
	if matchAttribs
		if !isempty(attrKey[1])
			if attrKey[1] == "*"
				keyAttribs = map(LightXML.name, LightXML.attributes(xElem))
			else
				keyAttribs = split(attrKey[1], ",")
			end
		end
		key = attrKey[2]
	end
	if keyIndex == length(keyArr) && isempty(key)
		attr = xml_join_attribute_values(xElem, keyAttribs)
		if attr != nothing
			add_dataframe_row!(results, keyAttribs, attr)
			return
		end
	end
	if key == "*"
		children = LightXML.child_elements(xElem)
	elseif matchAttribs
		children = xml_find_by_attribute(xElem, split(key, ","), keyAttribs)
	else
		children = LightXML.get_elements_by_tagname(xElem, key)
	end
	for child in children
		if keyIndex == length(keyArr)
			add_dataframe_row!(results, ["_content"], [LightXML.content(child)])
		else
			xml_find(child, keyArr, keyIndex+1, results)
		end
	end
end

function convert_column_type(arr::DataFrames.DataArray)
	types       = [Int, Float64]
	minInd = 1
	for e in arr
		if DataFrames.isna(e)
			continue
		end
		i = minInd
		while i <= length(types)
			try
				parse(types[i], e)
				break
			catch
				minInd = i+1
				if minInd > length(types)
					return arr
				end
			end
			i += 1
		end
	end
	newArr = DataFrames.DataArray(types[minInd], length(arr))
	for i = 1:length(arr)
		if !DataFrames.isna(arr[i])
			newArr[i] = parse(types[minInd], arr[i])
		end
	end
	return newArr
end

# Matched paths are slash-delimited. Elements are of the kind [attribNames][:][attribValuesOrTag]
# attribNames and attribValuesOrTag are comma-delimited lists of attribute names and values to match in the parent's children.
# If attribNames is omitted, then the value of the parent's key attribute is used for the names, or if the parent has no key attribute,
# attribNames defaults to "name".
# attribValuesOrTag can be empty only for the last element, in which case a comma-separated list of the values of attributes in attribNames is returned.
# If there isn't a colon in the element, then attribValuesOrTag is a tag name to match.
# If the last element is a tag, then the contents of the tagged node are returned.
# If attribValuesOrTag="*", then all child nodes are processed for matches. If attribNames="*", then all current node's attributes are selected
# (in the order they appear in the node).

function xml_find(xElem::LightXML.XMLElement, path::AbstractString)
	results = DataFrames.DataFrame()
	xml_find(xElem, split(path, "/"; keep = false), 1, results)
	for colName in names(results)
		results[colName] = convert_column_type(results[colName])
	end
	return results
end

global userAgent = "EveInvent Julia APIv2 library, greyalex2003@yahoo.com"
set_api_user_agent(agent::AbstractString) = (global userAgent = agent)

function get_api_url(url::AbstractString; params::Dict = Dict(), auth::Dict = Dict(), ignoreCache::Bool = false)
	fileName = api_cache_name(url)
	data = ignoreCache? nothing : xml_read(fileName)
	if data == nothing
		params = merge(params, auth)
		dataStr = Requests.format_query_str(params)
		headers = Dict("Content-Type"=>"application/x-www-form-urlencoded",
				   "Content-Length"=>string(length(dataStr)),
				   "User-Agent"=>userAgent)
		info("POST to $url with data: $dataStr")
		resp = Requests.post(url; data=dataStr, headers=headers)
		if resp.status != 200
			info("Response status: $(resp.status)")
			return nothing
		end
		xml_write(fileName, resp.data)
		xdoc = LightXML.parse_string(String(resp.data))
		data = LightXML.root(xdoc)
	end
	return data
end

get_api(endpoint::AbstractString; params::Dict = Dict(), auth::Dict = Dict(), ignoreCache::Bool = false) = get_api_url(urlApi * endpoint * ".xml.aspx"; params=params, auth=auth, ignoreCache=ignoreCache)

ext: jl


text:
module  TestTypesDatasets

using PopGenCore
using DataFrames
using PooledArrays
using Test

cats = @nancycats;
sharks = @gulfsharks;

@testset "Data Types and Datasets" begin
    @testset "Dataset dimensions" begin
        @test size(cats.sampleinfo) == (237,3)
        @test size(cats.genodata) == (2133,4)
        @test size(sharks.sampleinfo) == (212,5)
        @test size(sharks.genodata) == (468308,4)
    end

    @testset "Dataset column names" begin
        meta_names_sorted = ["latitude","longitude","name","ploidy","population"]
        loci_names_sorted = ["genotype","locus","name","population"]
        @test sort(names(cats.sampleinfo)) == meta_names_sorted[3:end]
        @test sort(names(cats.genodata)) == loci_names_sorted
        @test sort(names(sharks.sampleinfo)) == meta_names_sorted
        @test sort(names(sharks.genodata)) == loci_names_sorted
    end

    @testset "Nancycats column types" begin
        @test eltype(cats.sampleinfo.name) <: AbstractString
        @test eltype(cats.sampleinfo.population) <: AbstractString
        @test typeof(cats.sampleinfo.ploidy) == Vector{Int8}
        @test_throws ArgumentError typeof(cats.sampleinfo.latitude)
        @test_throws ArgumentError typeof(cats.sampleinfo.longitude)
        @test typeof(cats.genodata.name) <: PooledArray
        @test eltype(cats.genodata.name) <: AbstractString
        @test typeof(cats.genodata.population) <: PooledArray
        @test eltype(cats.genodata.population) <: AbstractString
        @test typeof(cats.genodata.locus) <: PooledArray
        @test eltype(cats.genodata.locus) <: AbstractString
        @test typeof(cats.genodata.genotype) <: GenoArray
        @test eltype(cats.genodata.genotype) <: Union{Missing, Genotype}
    end

    @testset "Gulfsharks column types" begin
        @test eltype(sharks.sampleinfo.name) <: AbstractString
        @test eltype(sharks.sampleinfo.population) <: AbstractString
        @test eltype(sharks.sampleinfo.ploidy) == Int8
        @test eltype(sharks.sampleinfo.latitude) ==  Float64
        @test eltype(sharks.sampleinfo.longitude) == Float64
        @test typeof(sharks.genodata.name) <: PooledArray
        @test eltype(sharks.genodata.name) <: AbstractString
        @test typeof(sharks.genodata.population) <: PooledArray
        @test eltype(sharks.genodata.population) <: AbstractString
        @test typeof(sharks.genodata.locus) <: PooledArray
        @test eltype(sharks.genodata.locus) <: AbstractString
        @test typeof(sharks.genodata.genotype) <: GenoArray
        @test typeof(sharks.genodata.genotype) <: PooledArray
        @test eltype(sharks.genodata.genotype) <: Union{Missing, Genotype}
    end
end

end  # module


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ThingType <: AbstractThingType
    label::Label
    is_root::Bool
end

# Porting note: Java client calls into RequestBuilder but it really
# has nothing to to with requests... I think it's probably better
# migrating the function here.
function proto(t::AbstractThingType)
    Proto._Type(
        label = t.label.name,
        encoding = encoding(t)
    )
    # return ThingTypeRequestBuilder.proto_thing_type(t.label, encoding(t))
end

# Contract: all subtypes of AbstractThingType should have these two fields
is_root(t::AbstractThingType) = t.is_root
label(t::AbstractThingType) = t.label

# ------------------------------------------------------------------------
# Remote functions
# ------------------------------------------------------------------------

function set_supertype(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = ThingTypeRequestBuilder.set_supertype_req(r.concept.label)
    res = execute(r.transaction, req)
end

function get_supertype(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = TypeRequestBuilder.get_supertype_req(r.concept.label)
    res = execute(r.transaction, req)
    typ = res.type_res.type_get_supertype_res._type
    return instantiate(typ)
end

function get_supertypes(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = TypeRequestBuilder.get_supertypes_req(r.concept.label)
    res = execute(r.transaction, req)
    typs = res.type_res_part.type_get_supertypes_res_part.types
    return instantiate.(typs)
end

function get_subtypes(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = TypeRequestBuilder.get_subtypes_req(r.concept.label)
    res = execute(r.transaction, req)
    typs = res.type_res_part.type_get_subtypes_res_part.types
    return instantiate.(typs)
end

function get_instances(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = ThingTypeRequestBuilder.get_instances_req(r.concept.label)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_instances_res_part.things for r in res)))
end

function set_abstract(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = ThingTypeRequestBuilder.set_abstract_req(r.concept.label)
    execute(r.transaction, req)
end

function unset_abstract(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = ThingTypeRequestBuilder.unset_abstract_req(r.concept.label)
    execute(r.transaction, req)
end

function set_plays(
    r::RemoteConcept{C,T},
    role_type::AbstractRoleType,
    overridden_role_type::Optional{AbstractRoleType} = nothing
) where {C <: AbstractThingType, T <: AbstractCoreTransaction}
    req = ThingTypeRequestBuilder.set_plays_req(
        r.concept.label,
        proto(role_type),
        overridden_role_type
    )
    execute(r.transaction, req)
end

function set_owns(
    r::RemoteConcept{C,T},
    attribute_type::AbstractAttributeType,
    is_key::Bool = false
) where {C <: AbstractThingType, T <: AbstractCoreTransaction}
    req = ThingTypeRequestBuilder.set_owns_req(
        r.concept.label,
        is_key,
        proto(attribute_type)
    )
    execute(r.transaction, req)
end

function get_owns(
    r::RemoteConcept{C,T},
    value_type::Optional{EnumType} = nothing,
    keys_only::Bool = false
) where {C <: AbstractThingType, T <: AbstractCoreTransaction}
    req = ThingTypeRequestBuilder.get_owns_req(r.concept.label, value_type, keys_only)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_owns_res_part.attribute_types for r in res)))
end

function get_plays(r::RemoteConcept{C,T}) where {
    C <: AbstractThingType, T <: AbstractCoreTransaction
}
    req = ThingTypeRequestBuilder.get_plays_req(r.concept.label)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_plays_res_part.roles for r in res)))
end


text:
function latexprint(io::IO, ct::CoefTable; tableenvir::AbstractString = "tabular", alignment::AbstractString="", caption::AbstractString="", centering::Bool=true)

    cols = ct.cols; rownms = ct.rownms; colnms = ct.colnms;
    nc = length(cols)
    nr = length(cols[1])
    if length(rownms) == 0
        rownms = [lpad("[$i]",floor(Integer, log10(nr))+3) for i in 1:nr]
    end
    rnwidth = max(4,maximum([length(nm) for nm in rownms]) + 1)
    rownms = [rpad(nm,rnwidth) for nm in rownms]
    widths = [length(cn)::Int for cn in colnms]
    str = AbstractString[isa(cols[j][i], AbstractString) ? cols[j][i] :
        sprint(showcompact,cols[j][i]) for i in 1:nr, j in 1:nc]
    for j in 1:nc
        for i in 1:nr
            lij = length(str[i,j])
            if lij > widths[j]
                widths[j] = lij
            end
        end
    end
    widths .+= 1
    @printf(io, "\\begin{%s}{%s}\n", tableenvir, alignment)
    length(caption)>0 && @printf(io, "\n\\caption{%s}", caption)
    (centering && length(caption)>0) && println(io,"\n\\centering")
    @printf(io, "\n%s\n", "\\hline")

    println(io, " " ^ rnwidth * "&" *
            join([lpad(string(colnms[i]), widths[i]) for i = 1:nc], "&") * "\\\\")
    for i = 1:nr
        println(io, rownms[i] * "&" *
                join([lpad(string(str[i, j]), widths[j]) for j = 1:nc], "&") * "\\\\")
    end
    @printf(io, "\n%s", "\\hline")
    @printf(io, "\n\\end{%s}", tableenvir)
end

latexprint(ct; kwargs...) = latexprint(STDIN, ct; kwargs...)

function latexprint(io::IO, m::MixtureModel; rounding::Int = 3, surrounding::Bool=false)
    p = round.(probs(m), rounding)
    μ = Float64[round(m.components[i].μ, rounding) for i in eachindex(p)]
    σ = Float64[round(m.components[i].σ, rounding) for i in eachindex(p)]
    surrounding && println(io, "\\[")
    println(io, join( [string(p[i]) * " Normal(" * string(μ[i]) * ", " * string(σ[i]) * "^2)" for i in eachindex(p)], " + "))
    surrounding && println(io, "\\]")
end


text:
using Pkg
pkg"activate .."
push!(LOAD_PATH, "../src/")

using Documenter, GeneralizedMonteCarlo

makedocs(
    sitename = "Generalized Monte Carlo",
    modules = [GeneralizedMonteCarlo],
    authors = "Scott Greenberg",
    pages = [
        "index.md",
    ]
)


text:
function alignment_objective(X, gmmx::AbstractGMM, gmmy::AbstractGMM, rot=nothing, trl=nothing, pσ=nothing, pϕ=nothing)
    return -overlap(AffineMap(X...)(gmmx), gmmy, pσ, pϕ)
end

# alignment objective for rigid rotation (i.e. the first stage of TIV-GOGMA)
function rot_alignment_objective(X, gmmx::AbstractGMM, gmmy::AbstractGMM, rot=nothing, trl=nothing, pσ=nothing, pϕ=nothing)
    if isnothing(trl)
        zro = zero(promote_type(numbertype(gmmx), numbertype(gmmy)))
        trl = (zro, zro, zro)
    end
    return alignment_objective((X..., trl...), gmmx, gmmy, pσ, pϕ)
end

# alignment objective for translation (i.e. the second stage of TIV-GOGMA)
function trl_alignment_objective(X, gmmx::AbstractGMM, gmmy::AbstractGMM, rot=nothing, trl=nothing, pσ=nothing, pϕ=nothing)
    if isnothing(rot)
        zro = zero(promote_type(numbertype(gmmx), numbertype(gmmy)))
        rot = (zro, zro, zro)
    end
    return alignment_objective((rot..., X...), gmmx, gmmy, pσ, pϕ)
end

"""
    obj, pos = local_align(gmmx, gmmy, block, pσ=nothing, pϕ=nothing; objfun=alignment_objective, rot=nothing, trl=nothing, rtol=1e-9, maxevals=100)

Performs local alignment within the specified `block` using L-BFGS to minimize objective function `objfun` for the provided GMMs, `gmmx` and `gmmy`.
"""
function local_align(gmmx::AbstractGMM, gmmy::AbstractGMM, block=fullBlock(gmmx,gmmy), pσ=nothing, pϕ=nothing;
                     objfun=alignment_objective, rot=nothing, trl=nothing, maxevals=100)
    # prepare pairwise widths and weights
    if isnothing(pσ) || isnothing(pϕ)
        pσ, pϕ = pairwise_consts(gmmx, gmmy)
    end

    # set optimization bounds
    lower = [r[1] for r in block.ranges]
    upper = [r[2] for r in block.ranges]

    # set initial guess at the center of the block
    initial_X = [x for x in block.center]

    # local optimization within the block
    f(X) = objfun(X, gmmx, gmmy, rot, trl, pσ, pϕ)
    res = optimize(f, lower, upper, initial_X, Fminbox(LBFGS()), Optim.Options(f_calls_limit=maxevals); autodiff = :forward)
    return Optim.minimum(res), tuple(Optim.minimizer(res)...)
end

local_align(gmmx::AbstractGMM, gmmys::Union{AbstractVector{AbstractSingleGMM},AbstractVector{AbstractMultiGMM}}, block,  pσ=nothing, pϕ=nothing; kwargs...) =
    local_align(gmmx, combine(gmmys), block, pσ, pϕ; kwargs...)


text:
# Closest Correlation matrix test matrix

# Original problem has the following format:
# min_X   1/2 ||X-C||^2
# s.t.    Xii = 1
#         X ⪴ 0

using COSMO, Test, LinearAlgebra, SparseArrays, Random

# This test is precision agnostic and will be run with TestFloat precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

# this function creates a matrix A that slices out the diagonal entries Xii of a vectorized square matrix x=vec(X)
function create_diagonal_extractor(n, type::Type{T} = Float64) where {T <: AbstractFloat}
  A = spzeros(type, n, n^2)
  A[1,1 ] = 1
  for iii = 2:n-1
    col = (iii - 1) * (n + 1)
    A[iii, col + 1] = 1
  end
  A[n, n^2] = 1
  return A
end

function diagonal_equal_one(A, atol)
  for iii = 1:size(A, 1)
    (abs(A[iii, iii] - 1) > atol) && return false
  end
  return true
end
@testset "SDP - Closest Correlation Matrix" begin

for T in UnitTestFloats
  if T != BigFloat
    rng = Random.MersenneTwister(12345)
    xMin = -one(T)
    xMax = one(T)
    n = 50
    C = xMin .+ randn(rng, T, n, n)*(xMax - xMin)
    c = vec(C)

    n2 = n^2
    m = n + n2

    P = spdiagm(0 => ones(T, n2))
    q = -vec(C)
    r = T(0.5) * dot(c, c)

    A1 = create_diagonal_extractor(n, T)
    b1 = -ones(T, n)
    cs1 = COSMO.Constraint(A1, b1, COSMO.ZeroSet)

    A2 = spdiagm(0 => ones(T, n2))
    b2 = zeros(T, n2)
    cs2 = COSMO.Constraint(A2, b2, COSMO.PsdCone)
    constraints = [cs1; cs2]

    model = COSMO.Model{T}()
    assemble!(model, P, q, constraints, settings = COSMO.Settings{T}(eps_abs = T(1e-4), eps_rel = T(1e-4)))

    res = COSMO.optimize!(model)

    Xsol = reshape(res.x, n, n)

    @testset "SDP - Closest Correlation Matrix (T = $(T))" begin
      @test res.status == :Solved
      @test diagonal_equal_one(Xsol, 1e-5)
      @test minimum(real.(eigen(Xsol).values)) > -1e-3
    end
  end
end
end
nothing


text:
"""
    get_classes(divblock)

Return the classe(s) of a div block. E.g. `@@c1,c2` will return `"c1 c2"` so
that it can be injected in a `<div class="..."`.
"""
function get_classes(b::Block)::String
    return replace(subs(b.open.ss, 3:lastindex(b.open.ss)), "," => " ")
end

"""
    prepare_text(blocks)

For a text block, replace the remaining tokens for special characters.
"""
function prepare_text(b::Block; tohtml=true)::String
    c = b.ss
    # if there's no tokens over that content, return
    isempty(b.inner_tokens) && return String(c)
    # otherwise inject as appropriate
    parent = parent_string(c)
    io = IOBuffer()
    head = from(c)
    for t in b.inner_tokens
        t.name in MD_IGNORE && continue
        write(io, subs(parent, head, prev_index(t)))
        write(io, insert(t; tohtml))
        head = next_index(t)
    end
    write(io, subs(parent, head, to(c)))
    return String(take!(io))
end

"""
    insert(token)

For tokens representing special characters, insert the relevant string.
"""
function insert(t::Token; tohtml=true)::String
    stname = String(t.name)
    s = String(t.ss)  # safe default
    if t.name == :CHAR_HTML_ENTITY
        s = String(t.ss)
    elseif startswith(stname, "CHAR_") && tohtml # CHAR_*
        id = stname[6:end]
        s = "&#$(id);"
    elseif t.name == :CAND_EMOJI
        # check if it's a valid emoji
        s = get(emoji_symbols, "\\$(t.ss)", s)
    end
    tohtml || (s = replace(s, "&" => "\\&"))
    return s
end


text:
include("meshing/meshing.jl")
include("init_model.jl")
include("init_receiver.jl")
include("init_source.jl")


text:
using Observers
using DataFrames
using Test
using JLD2

@testset "Observer" begin
  # Series for π/4
  f(k) = (-1)^(k+1)/(2k-1)

  function my_iterative_function(niter; observer!, observe_step)
    π_approx = 0.0
    for n in 1:niter
      π_approx += f(n)
      if iszero(n % observe_step)
        update!(observer!; π_approx = 4π_approx, iteration = n)
      end
    end
    return 4π_approx
  end

  # Measure the relative error from π at each iteration
  err_from_π(; π_approx, kwargs...) = abs(π - π_approx) / π

  # Record which iteration we are at
  iteration(; iteration, kwargs...) = iteration

  obs = Observer(["Error" => err_from_π, "Iteration" => iteration])

  @test length(obs) == 2
  @test results(obs, "Error") == []
  @test results(obs, "Iteration") == []

  niter = 10000
  observe_step = 1000
  π_approx = my_iterative_function(niter; observer! = obs, observe_step = observe_step)

  for res in obs
    @test length(last(last(res))) == niter ÷ observe_step
  end

  save("outputdata.jld2", results(obs))
  obs_load = load("outputdata.jld2")

  for (k,v) in obs_load
    @test v ≈ last(obs[k])
  end

  obs = Observer(["Error" => err_from_π, "Iteration" => iteration])
  obs["nofunction"] = nothing

  niter = 10000
  observe_step = 1000
  π_approx = my_iterative_function(niter; observer! = obs, observe_step = observe_step)

  @test results(obs,"nofunction") == []

  #
  # List syntax
  #

  obs = Observer("Error" => err_from_π, "Iteration" => iteration)

  @test length(obs) == 2
  @test results(obs, "Error") == []
  @test results(obs, "Iteration") == []

  niter = 10000
  observe_step = 1000
  π_approx = my_iterative_function(niter; observer! = obs, observe_step = observe_step)

  for res in obs
    @test length(last(last(res))) == niter ÷ observe_step
  end

  #
  # Function list syntax
  #

  obs = Observer(err_from_π, iteration)

  @test length(obs) == 2
  @test results(obs, err_from_π) == []
  @test results(obs, "err_from_π") == []
  @test results(obs, iteration) == []
  @test results(obs, "iteration") == []

  niter = 10000
  observe_step = 1000
  π_approx = my_iterative_function(niter; observer! = obs, observe_step = observe_step)

  for res in obs
    @test length(last(last(res))) == niter ÷ observe_step
  end

  f1(x::Int) = x^2
  f2(x::Int, y::Float64) = x + y
  f3(_,_,t::Tuple) = first(t)
  f4(x::Int; a::Float64) = x * a
  f5(x::Int; a::Float64, b::Float64) = x * a + b

  function my_other_iterative_function(; observer!)
    k = 2
    x = k
    y = k * √2
    t = (x+2*y,0,0)
    a = y^2
    b = 3.0
    update!(observer!, x, y, t; a = a, b = b)
  end

  obs = Observer(["f1" => f1, "f2" => f2, "f3" => f3, "f4" => f4, "f5" => f5])

  my_other_iterative_function(; observer! = obs)

  @test results(obs,"f1")[1] ≈ 4.0
  @test results(obs,"f2")[1] ≈ 2.0 + 2*√2
  @test results(obs,"f3")[1] ≈ 2.0 + 4*√2
  @test results(obs,"f4")[1] ≈ 16.0
  @test results(obs,"f5")[1] ≈ 19.0
end

@testset "Observer constructed from functions" begin
  # Series for π/4
  f(k) = (-1)^(k+1)/(2k-1)

  function my_iterative_function(niter; observer!, observe_step)
    π_approx = 0.0
    for n in 1:niter
      π_approx += f(n)
      if iszero(n % observe_step)
        update!(observer!; π_approx = 4π_approx, iteration = n)
      end
    end
    return 4π_approx
  end

  # Measure the relative error from π at each iteration
  err_from_π(; π_approx, kwargs...) = abs(π - π_approx) / π

  # Record which iteration we are at
  iteration(; iteration, kwargs...) = iteration
  obs = Observer([err_from_π, iteration])

  @test length(obs) == 2
  @test results(obs, "err_from_π") == []
  @test results(obs, "iteration") == []

  niter = 10000
  observe_step = 1000
  π_approx = my_iterative_function(niter; observer! = obs, observe_step = observe_step)

  @test length(results(obs, "err_from_π")) == niter ÷ observe_step
  @test length(results(obs, "iteration")) == niter ÷ observe_step
  @test length(results(obs, err_from_π)) == niter ÷ observe_step
  @test length(results(obs, iteration)) == niter ÷ observe_step
  f1 = err_from_π
  f2 = iteration
  @test length(results(obs, f1)) == niter ÷ observe_step
  @test length(results(obs, f2)) == niter ÷ observe_step
  @test_throws KeyError results(obs, "f1")
  @test_throws KeyError results(obs, "f2")
end

@testset "save only last value" begin
  f(x::Int, y::Float64) = x + y
  g(x::Int) = x

  function Observers.update!(::typeof(g), results, result)
    empty!(results)
    push!(results, result)
  end

  function my_yet_another_iterative_function(niter::Int; observer!)
    for k in 1:niter
      x = k
      y = k * √2
      update!(observer!, x, y)
    end
  end

  obs = Observer(["g" => g, "f" => f])

  my_yet_another_iterative_function(100; observer! = obs)

  @test length( results(obs,"g")) == 1
  @test length( results(obs,"f")) == 100
end

@testset "empty" begin
  f(x) = 2x
  function iterative(niter; observer!)
    for k in 1:niter
      update!(observer!, k)
    end
  end
  obs0 = Observer(["f" => f])

  obs1 = copy(obs0)
  @test obs0 == obs1
  iterative(10; observer! = obs1)
  @test obs1 ≠ obs0
  empty_results!(obs1)
  @test obs1 == obs0

  iterative(10; observer! = obs1)
  obs1 = empty_results(obs0)
  @test obs1 == obs0
end

@testset "Test element types of Array" begin
  f(k, x::Int, y::Float64) = x + y
  g(k, x) = k < 20 ? 0 : exp(im * x)

  function iterative_function(niter::Int; observer!)
    for k in 1:niter
      x = k
      y = x * √2
      update!(observer!, k, x, y)
    end
  end

  obs = Observer(["f" => f, "g" => g])
  iterative_function(100; observer! = obs)
  res = results(obs)
  @test length(res["f"]) == 100
  @test length(res["g"]) == 100
  @test res["f"] isa Vector{Float64}
  @test res["g"] isa Vector{ComplexF64}

  obs = Observer(["f" => f, "g" => g])
  iterative_function(10; observer! = obs)
  res = results(obs)
  @test length(res["f"]) == 10
  @test length(res["g"]) == 10
  @test res["f"] isa Vector{Float64}
  @test res["g"] isa Vector{Int}
end

@testset "Conversion to DataFrames" begin
  f(k, x::Int, y::Float64) = x + y
  g(k, x) = k < 20 ? 0 : exp(im * x)
  function iterative_function(niter::Int; observer!)
    for k in 1:niter
      x = k
      y = x * √2
      update!(observer!, k, x, y)
    end
  end
  obs = Observer(["f" => f, "g" => g])
  iterative_function(100; observer! = obs)
  res = results(obs)
  df = DataFrame(results(obs))
  @test df.f == res["f"]
  @test df.g == res["g"]
  @test length(df.f) == 100
  @test length(df.g) == 100
  @test df.f isa Vector{Float64}
  @test df.g isa Vector{ComplexF64}
end


text:
function general_BlockDiagonal_tests(rng, blocks)
    d = block_diagonal(blocks)
    Ps, Qs = size.(blocks, 1), size.(blocks, 2)

    @testset "general" begin
        @test blocksizes(d, 1) == Ps
        @test blocksizes(d, 2) == Qs

        @test view(d, Block(1, 1)) == blocks[1]
        @test view(d, Block(2, 2)) == blocks[2]
        @test view(d, Block(1, 2)) == zeros(Ps[1], Qs[2])
        @test view(d, Block(2, 1)) == zeros(Ps[2], Qs[1])

        @test d[Block(1, 1)] == view(d, Block(1, 1))
    end
end

function BlockDiagonal_mul_tests(rng, blocks)
    D, Ps = block_diagonal(blocks), size.(blocks, 1)
    Dmat = Matrix(D)

    U = UpperTriangular(D)

    xs, ys = [randn(rng, P) for P in Ps], [randn(rng, P) for P in Ps]
    y, x = mortar(ys), mortar(xs)

    # Matrix-Vector product
    @test mul!(y, D, x) ≈ Dmat * Vector(x)
    @test mul!(y, D, x) == D * x
    @test mul!(y, U, x) ≈ Matrix(U) * Vector(x)
    @test mul!(y, U, x) == U * x

    Qs = [3, 4]
    X = mortar([randn(rng, P, Q) for P in Ps, Q in Qs])
    Y = mortar([randn(rng, P, Q) for P in Ps, Q in Qs])

    # Matrix-Matrix product
    @test mul!(Y, D, X) ≈ Dmat * X
    @test mul!(Y, U, X) ≈ Matrix(U) * Matrix(X)
    @test mul!(Y, D, X) == D * X
    @test mul!(Y, U, X) == U * X
end

function BlockDiagonal_chol_tests(rng, blocks)

    D, Ps = block_diagonal(blocks), size.(blocks, 1)
    Dmat = Matrix(D)

    C, Cmat = cholesky(D), cholesky(Dmat)

    @test C.U ≈ Cmat.U
    @test logdet(C) ≈ logdet(Cmat)

    Csym = cholesky(Symmetric(D))
    @test C.U ≈ Csym.U

    # Test backprop for accessing `U`.
    U_diag, back_diag = Zygote.pullback(D->cholesky(D).U, D)
    U_dens, back_dens = Zygote.pullback(D->cholesky(D).U, Matrix(D))

    @test U_diag ≈ U_dens

    Ū = block_diagonal([randn(rng, P, P) for P in Ps])
    D̄_diag = first(back_diag(Ū))
    D̄_dens = first(back_dens(Matrix(Ū)))
    @test Matrix(D̄_diag) ≈ D̄_dens
    @test D̄_diag isa BlockDiagonal

    # Test backprop for logdet of a Cholesky.
    l_diag, l_back_diag = Zygote.pullback(D->logdet(cholesky(D)), D)
    l_dens, l_back_dens = Zygote.pullback(D->logdet(cholesky(D)), Matrix(D))

    @test l_diag ≈ l_dens

    l̄ = randn(rng)
    D̄_diag = first(l_back_diag(l̄))
    D̄_dens = first(l_back_dens(l̄))
    @test Matrix(D̄_diag) ≈ D̄_dens
    @test D̄_diag isa BlockDiagonal
end

function BlockDiagonal_add_tests(rng, blks; grad=true)

    D = block_diagonal(blks)
    Dmat = Matrix(D)
    A = randn(rng, size(D))

    A_copy = copy(A)
    C = A_copy + D
    @test A_copy == A
    @test C == A + Dmat

    if grad == true
        @assert length(blks) == 2
        adjoint_test(
            (A, b1, b2)->A + block_diagonal([b1, b2]),
            randn(rng, size(A)), A, blks[1], blks[2],
        )
    end
end

@timedtestset "BlockDiagonal" begin
    @timedtestset "Matrix" begin
        rng, Ps, Qs = MersenneTwister(123456), [2, 3], [4, 5]
        vs = [randn(rng, Ps[1], Qs[1]), randn(rng, Ps[2], Qs[2])]
        general_BlockDiagonal_tests(rng, vs)

        As = [randn(rng, Ps[n], Ps[n]) for n in eachindex(Ps)]
        blks = [As[n] * As[n]' + I for n in eachindex(As)]
        BlockDiagonal_mul_tests(rng, blks)
        BlockDiagonal_mul_tests(rng, UpperTriangular.(blks))
        BlockDiagonal_mul_tests(rng, Hermitian.(blks))
        BlockDiagonal_mul_tests(rng, Symmetric.(blks))
        BlockDiagonal_chol_tests(rng, blks)
        BlockDiagonal_add_tests(rng, blks; grad=false)
    end
    @timedtestset "Diagonal{T, <:Vector{T}}" begin
        rng, Ps = MersenneTwister(123456), [2, 3]
        vs = [Diagonal(randn(rng, Ps[n])) for n in eachindex(Ps)]
        general_BlockDiagonal_tests(rng, vs)

        blocks = [Diagonal(ones(P) + exp.(randn(rng, P))) for P in Ps]
        BlockDiagonal_add_tests(rng, blocks; grad=false)
        @timedtestset "cholesky" begin
            x, ȳ = randn(rng, sum(Ps)), randn(rng, sum(Ps))
            adjoint_test((X, blks)->cholesky(block_diagonal(blks)).U \ X, ȳ, x, blocks)

            X, Ȳ = randn(rng, sum(Ps), 7), randn(rng, sum(Ps), 7)
            adjoint_test((X, blks)->cholesky(block_diagonal(blks)).U \ X, Ȳ, X, blocks)
            adjoint_test(blks->logdet(cholesky(block_diagonal(blks))), randn(rng), blocks)
        end
    end
    @timedtestset "Negation" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6, 7]
        A = block_diagonal([randn(rng, P, P) for P in Ps])

        @test Matrix(-A) == -Matrix(A)
        @test -A isa BlockDiagonal

        Y_diag, back_diag = Zygote.pullback(-, A)
        Y_dens, back_dens = Zygote.pullback(-, Matrix(A))

        Ȳ = block_diagonal([randn(rng, P, P) for P in Ps])
        @test Y_diag == -A
        @test Matrix(first(back_diag(Ȳ))) == first(back_dens(Matrix(Ȳ)))
        @test first(back_diag(Ȳ)) isa BlockDiagonal
    end
    @timedtestset "adjoint" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])

        @test Matrix(A') == Matrix(A)'
        @test A' isa BlockDiagonal

        Y, back = Zygote.pullback(adjoint, A)
        Y_dens, back_dens = Zygote.pullback(adjoint, Matrix(A))
        Ȳ = block_diagonal([randn(rng, P, P) for P in Ps])
        @test Y == A'
        @test Matrix(first(back(Ȳ))) == first(back_dens(Matrix(Ȳ)))
        @test first(back(Ȳ)) isa BlockDiagonal
    end
    @timedtestset "transpose" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])

        @test Matrix(transpose(A)) == transpose(Matrix(A))
        @test transpose(A) isa BlockDiagonal

        Y, back = Zygote.pullback(transpose, A)
        Y_dens, back_dens = Zygote.pullback(transpose, Matrix(A))
        Ȳ = block_diagonal([randn(rng, P, P) for P in Ps])
        @test Y == transpose(A)
        @test Matrix(first(back(Ȳ))) == first(back_dens(Matrix(Ȳ)))
        @test first(back(Ȳ)) isa BlockDiagonal
    end
    @timedtestset "UpperTriangular" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])

        @test Matrix(UpperTriangular(A)) == UpperTriangular(Matrix(A))
        @test UpperTriangular(A) isa BlockDiagonal

        B_diag, back_diag = Zygote.pullback(UpperTriangular, A)
        B_dens, back_dens = Zygote.pullback(UpperTriangular, Matrix(A))
        @test Matrix(B_diag) == B_dens

        B̄ = block_diagonal([randn(rng, P, P) for P in Ps])
        Ā_diag, Ā_dens = first(back_diag(B̄)), first(back_dens(Matrix(B̄)))
        @test Ā_diag == Ā_dens
        @test Ā_diag isa BlockDiagonal
    end
    @timedtestset "Symmetric" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        S = Symmetric(A)
        @test S == Symmetric(Matrix(A))
        @test S isa BlockDiagonal

        S_diag, back_diag = Zygote.pullback(Symmetric, A)
        S_dens, back_dens = Zygote.pullback(Symmetric, Matrix(A))
        @test S_diag ≈ S_dens
    end
    @timedtestset "BlockDiagonal * BlockDiagonal" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        B = block_diagonal([randn(rng, P, P) for P in Ps])

        @test Matrix(A * B) ≈ Matrix(A) * Matrix(B)
        @test A * B isa BlockDiagonal

        Y_diag, back_diag = Zygote.pullback(*, A, B)
        Y_dens, back_dens = Zygote.pullback(*, Matrix(A), Matrix(B))

        Ȳ = block_diagonal([randn(rng, P, P) for P in Ps])
        @test Y_diag == A * B

        Ā_diag, B̄_diag = back_diag(Ȳ)
        Ā_dens, B̄_dens = back_dens(Matrix(Ȳ))

        @test Matrix(Ā_diag) ≈ Ā_dens
        @test Matrix(B̄_diag) ≈ B̄_dens

        @test Ā_diag isa BlockDiagonal
        @test B̄_diag isa BlockDiagonal
    end
    @timedtestset "BlockDiagonal * Matrix" begin
        rng, Ps, Q = MersenneTwister(123456), [4, 5, 6], 11
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        B = randn(rng, sum(Ps), Q)
        @test Matrix(A * B) ≈ Matrix(A) * B
        @test Matrix(A * collect(B')') ≈ Matrix(A) * B

        Y_diag, back_diag = Zygote.pullback(*, A, B)
        Y_dens, back_dens = Zygote.pullback(*, Matrix(A), B)
        @test Y_diag ≈ Y_dens

        Ȳ = randn(rng, sum(Ps), Q)
        Ā_diag, B̄_diag = back_diag(Ȳ)
        Ā_dens, B̄_dens = back_dens(Ȳ)
        @test Matrix(Ā_diag) ≈ Ā_dens
        @test B̄_diag ≈ B̄_dens
        @test_broken Ā_diag isa BlockDiagonal
    end
    @timedtestset "BlockDiagonal * Vector" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        x = randn(rng, sum(Ps))
        @test Vector(A * x) ≈ Matrix(A) * x

        Y_diag, back_diag = Zygote.pullback(*, A, x)
        Y_dens, back_dens = Zygote.pullback(*, Matrix(A), x)
        @test Y_diag ≈ Y_dens

        ȳ = randn(rng, sum(Ps))
        Ā_diag, x̄_diag = back_diag(ȳ)
        Ā_dens, x̄_dens = back_dens(ȳ)
        @test Matrix(Ā_diag) ≈ Ā_dens
        @test x̄_diag ≈ x̄_dens
        @test_broken Ā_diag isa BlockDiagonal
    end
    @timedtestset "ldiv(BlockDiagonal, Matrix)" begin
        rng, Ps, Q = MersenneTwister(123456), [4, 5, 6], 11
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        B = randn(rng, sum(Ps), Q)
        @test Matrix(A \ B) ≈ Matrix(A) \ B

        Y_diag, back_diag = Zygote.pullback(\, A, B)
        Y_dens, back_dens = Zygote.pullback(\, Matrix(A), B)
        @test Y_diag ≈ Y_dens

        Ȳ = randn(rng, sum(Ps), Q)
        Ā_diag, B̄_diag = back_diag(Ȳ)
        Ā_dens, B̄_dens = back_dens(Ȳ)
        @test_broken Matrix(Ā_diag) ≈ Ā_dens # we're not checking the right bits of the matrix here
        @test B̄_diag ≈ B̄_dens
        @test Ā_diag isa BlockDiagonal
        @test blocksizes(Ā_diag, 1) == blocksizes(A, 1)
        @test blocksizes(Ā_diag, 2) == blocksizes(A, 2)
    end
    @timedtestset "ldiv(BlockDiagonal, Vector)" begin
        rng, Ps = MersenneTwister(123456), [4, 5, 6]
        A = block_diagonal([randn(rng, P, P) for P in Ps])
        B = randn(rng, sum(Ps))
        @test Vector(A \ B) ≈ Matrix(A) \ B

        Y_diag, back_diag = Zygote.pullback(\, A, B)
        Y_dens, back_dens = Zygote.pullback(\, Matrix(A), B)
        @test Y_diag ≈ Y_dens

        Ȳ = randn(rng, sum(Ps))
        Ā_diag, B̄_diag = back_diag(Ȳ)
        Ā_dens, B̄_dens = back_dens(Ȳ)
        @test_broken Matrix(Ā_diag) ≈ Ā_dens # we're not checking the right bits of the matrix here
        @test B̄_diag ≈ B̄_dens
        @test Ā_diag isa BlockDiagonal
        @test blocksizes(Ā_diag, 1) == blocksizes(A, 1)
        @test blocksizes(Ā_diag, 2) == blocksizes(A, 2)
    end
end


text:
using Revise
using SortingLab
using Base.Test, BenchmarkTools

N = 100_000_000
K = 100

# test fsortperm
a = rand(1:Int(N/K), N);
ca = copy(a);
@code_warntype fsortperm(a)
@benchmark fsortperm($a) samples = 5 seconds  = 20
# @code_warntype SortingLab.fsortperm_int_range_counting_sort(a, Int(N/K), 1)

# @code_warntype SortingLab.fsortperm_msd_hybrid(a, Int(N/K), 1, 10)
# @code_warntype SortingLab.fsortperm_msd_hybrid(a, Int(N/K), 1, 10, 1024)
# @code_warntype SortingLab.fsortperm_int_range_lsd(a, Int(N/K), 1, 10)

@time ab = fsortperm(a);
@test ca == a
@test issorted(a[ab])

@time aba = SortingLab.fsortperm_int_range_lsd(a, Int(N/K), 1, 10)

a2 = rand(1:2048, N);
@time a2b = SortingLab.fsortperm(a2);
issorted(a2[a2b])

@time a2b = SortingLab.fsortperm(a2, rev = true);
issorted(a2[a2b], rev = true)

@benchmark fsortperm($a2) samples = 5 seconds  = 20

a2 = rand(1:500, N);
@time a2b = SortingLab.fsortperm(a2);
@benchmark fsortperm($a2) samples = 5 seconds  = 20

@time SortingLab.fsortperm_int_range_counting_sort(a2, 2048, 1)

@benchmark SortingLab.fsortperm_int_range_counting_sort($a2, 2048, 1)

# @time SortingLab.fsortperm_msd_hybrid(a, Int(N/K), 1, 10)
# @time SortingLab.fsortperm_msd_hybrid(a, Int(N/K), 1, 10, 1024)
# @time SortingLab.fsortperm_int_range_lsd(a, Int(N/K), 1, 10)

@benchmark SortingLab.fsortperm_msd_hybrid($a, Int(N/K), 1, 10) samples = 5 seconds = 15
@benchmark SortingLab.fsortperm_msd_hybrid($a, Int(N/K), 1, 10, 1024) samples = 5 seconds = 15
@benchmark SortingLab.fsortperm_msd_hybrid($a, Int(N/K), 1, 10, 2048) samples = 5 seconds = 15
@benchmark SortingLab.fsortperm_int_range_lsd($a, Int(N/K), 1, 10) samples = 5 seconds = 15

@benchmark fsortperm_msd_hybrid($a, Int(N/K), 1, 10) seconds = 3600
# BenchmarkTools.Trial:
#   memory estimate:  3.74 GiB
#   allocs estimate:  12711
#   --------------
#   minimum time:     2.069 s (8.46% GC)
#   mean time:        2.334 s (17.08% GC)
#   --------------
#   evals/sample:     1
@benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 1024) seconds = 3600
@benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 2048) seconds = 3600
@benchmark fsortperm_int_range_lsd($a, 1_000_000, 1, 10) seconds = 3600

@benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10) samples = 5 seconds = 15
@benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 1024) samples = 5 seconds = 15
@benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 2048) samples = 5 seconds = 15
@benchmark fsortperm_int_range_lsd($a, 1_000_000, 1, 10) samples = 5 seconds = 15
# BenchmarkTools.Trial:
#   memory estimate:  3.74 GiB
#   allocs estimate:  12711
#   --------------
#   minimum time:     2.448 s (10.43% GC)
#   median time:      2.530 s (15.98% GC)
#   mean time:        2.532 s (15.69% GC)
#   maximum time:     2.668 s (18.73% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

# Main> @benchmark fsortperm_int_range_lsd($a, 1_000_000, 1, 10) samples = 5 seconds = 15
# BenchmarkTools.Trial:
#   memory estimate:  2.24 GiB
#   allocs estimate:  13
#   --------------
#   minimum time:     2.256 s (0.49% GC)
#   median time:      2.399 s (9.61% GC)
#   mean time:        2.394 s (9.05% GC)
#   maximum time:     2.530 s (13.07% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

# Main> @benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 1024) samples = 5 seconds = 15
# BenchmarkTools.Trial:
#   memory estimate:  4.85 GiB
#   allocs estimate:  17599
#   --------------
#   minimum time:     2.595 s (8.32% GC)
#   median time:      2.795 s (14.35% GC)
#   mean time:        2.770 s (14.73% GC)
#   maximum time:     2.904 s (17.48% GC)
#   --------------
#   samples:          5
#   evals/sample:     1

# Main> @benchmark fsortperm_msd_hybrid($a, 1_000_000, 1, 10, 2048) samples = 5 seconds = 15
# BenchmarkTools.Trial:
#   memory estimate:  4.85 GiB
#   allocs estimate:  17599
#   --------------
#   minimum time:     2.592 s (7.99% GC)
#   median time:      2.779 s (15.41% GC)
#   mean time:        2.762 s (15.46% GC)
#   maximum time:     2.902 s (20.23% GC)
#   --------------
#   samples:          5
#   evals/sample:     1


text:
##
cd(@__DIR__())
codes = Dict(
    '.' => 1,
    'L' => 2,
    '#' => 3,
)
function check_cartesian(c, rows, cols)
    i, j = c.I
    1 <= i <= rows && 1 <= j <= cols && return c
    missing
end
directions = [
    CartesianIndex(-1,-1),
    CartesianIndex(-1,0),
    CartesianIndex(-1,1),
    CartesianIndex(0,-1),
    CartesianIndex(0,1),
    CartesianIndex(1,-1),
    CartesianIndex(1,0),
    CartesianIndex(1,1),
]
emptyseat(v) = v == codes['L']
occupiedseat(v) = v == codes['#']
isseat(v) = emptyseat(v) || occupiedseat(v)
function find_seat(c, d, grid)
    rows, cols = size(grid)
    cur = c + d
    while !ismissing(check_cartesian(cur, rows, cols)) && !isseat(grid[cur])
        cur += d
    end
    check_cartesian(cur, rows, cols)
end
function get_neighbours(c, grid)
    [find_seat(c, d, grid) for d ∈ directions]
end
function compute_neighbours(grid::AbstractMatrix{Int})
    rows, cols = size(grid)
    neigh = Array{Union{Missing,CartesianIndex{2}}}(undef, 8, rows, cols)
    for c ∈ CartesianIndices(grid)
        neigh[:, c] .= get_neighbours(c, grid)
    end
    neigh
end
function rules(v, ns)
    emptyseat(v) && return count(==(codes['#']), ns) == 0 ? codes['#'] : v
    occupiedseat(v) && return count(==(codes['#']), ns) >= 5 ? codes['L'] : v
    return v
end
flip((x,y)) = (y,x)
flip(c::CartesianIndex{2}) = CartesianIndex(flip(c.I))
function update_cells!(grid, old, neighbours, rules)
    copy!(old, grid)
    for c ∈ CartesianIndices(grid)
        ns = [old[n] for n ∈ skipmissing(view(neighbours, :, c))]
        grid[c] = rules(old[c], ns)
    end
end
function show_grid(grid)
    invcodes = ['.', 'L', '#']
    rows, cols = size(grid)
    for j ∈ 1:cols
        for i ∈ 1:rows
            print(invcodes[grid[i, j]])
        end
        print('\n')
    end
end

##
seats = reduce(hcat, map(eachline("input.txt")) do line
    map(c -> codes[c], collect(line))
end);
neighbours = compute_neighbours(seats);

##
grid = copy(seats)
old = fill!(similar(grid), 0);

##
while grid != old
    update_cells!(grid, old, neighbours, rules)
end

##
count(==(codes['#']), grid)