text:
## y'' = f(t, y, y')
## y(t₀) = y₀; y'(t₀) = y'₀
## kᵢ' = f(t₀+cᵢh, y₀+cᵢhy'₀+h²∑āᵢⱼk'ⱼ, y'₀+h∑aᵢⱼk'ⱼ)
## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
## y'₁ = y'₀ + h∑bᵢk'ᵢ

const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache,
                                           Nystrom4VelocityIndependentConstantCache,
                                           Nystrom5VelocityIndependentConstantCache,
                                           IRKN3ConstantCache,  IRKN4ConstantCache,
                                           DPRKN8ConstantCache, DPRKN12ConstantCache,
                                           ERKN4ConstantCache,  ERKN5ConstantCache}

function initialize!(integrator,cache::NystromCCDefaultInitialization)
 integrator.kshortsize = 2
 integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

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

const NystromDefaultInitialization = Union{Nystrom4Cache,
                                           Nystrom4VelocityIndependentCache,
                                           Nystrom5VelocityIndependentCache,
                                           IRKN3Cache,  IRKN4Cache,
                                           DPRKN8Cache, DPRKN12Cache,
                                           ERKN4Cache,  ERKN5Cache}

function initialize!(integrator,cache::NystromDefaultInitialization)
  @unpack fsalfirst,k = cache
  duprev,uprev = integrator.uprev.x

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

@muladd function perform_step!(integrator,cache::Nystrom4ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  k₁ = integrator.fsalfirst.x[1]
  halfdt = dt/2
  dtsq = dt^2
  eighth_dtsq = dtsq/8
  half_dtsq = dtsq/2
  ttmp = t+halfdt

  ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
  ku = uprev + halfdt*duprev + eighth_dtsq*k₁
  ## y'₁ = y'₀ + h∑bᵢk'ᵢ
  kdu = duprev + halfdt*k₁

  k₂ = f.f1(kdu,ku,p,ttmp)
  ku = uprev + halfdt*duprev + eighth_dtsq*k₁
  kdu = duprev + halfdt*k₂

  k₃ = f.f1(kdu,ku,p,ttmp)
  ku = uprev + dt*duprev + half_dtsq*k₃
  kdu = duprev + dt*k₃

  k₄ = f.f1(kdu,ku,p,t+dt)
  u = uprev + (dtsq/6)*(k₁ + k₂ + k₃) + dt*duprev
  du = duprev + (dt/6)*(k₁ + k₄ + 2*(k₂ + k₃))

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Nystrom4Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,fsalfirst,k₂,k₃,k₄,k = cache
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  k₁ = integrator.fsalfirst.x[1]
  halfdt = dt/2
  dtsq = dt^2
  eighth_dtsq = dtsq/8
  half_dtsq = dtsq/2
  ttmp = t+halfdt

  ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
  @.. broadcast=false ku = uprev + halfdt*duprev + eighth_dtsq*k₁
  ## y'₁ = y'₀ + h∑bᵢk'ᵢ
  @.. broadcast=false kdu = duprev + halfdt*k₁

  f.f1(k₂,kdu,ku,p,ttmp)
  @.. broadcast=false ku = uprev + halfdt*duprev + eighth_dtsq*k₁
  @.. broadcast=false kdu = duprev + halfdt*k₂

  f.f1(k₃,kdu,ku,p,ttmp)
  @.. broadcast=false ku = uprev + dt*duprev + half_dtsq*k₃
  @.. broadcast=false kdu = duprev + dt*k₃

  f.f1(k₄,kdu,ku,p,t+dt)
  @.. broadcast=false u = uprev + (dtsq/6)*(k₁ + k₂ + k₃) + dt*duprev
  @.. broadcast=false du = duprev + (dt/6)*(k₁ + k₄ + 2*(k₂ + k₃))

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

@muladd function perform_step!(integrator,cache::Nystrom4VelocityIndependentConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  k₁ = integrator.fsalfirst.x[1]
  halfdt = dt/2
  dtsq = dt^2
  eighth_dtsq = dtsq/8
  half_dtsq = dtsq/2
  ttmp = t+halfdt

  ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
  ku = uprev + halfdt*duprev + eighth_dtsq*k₁

  k₂ = f.f1(duprev,ku,p,ttmp)
  ku = uprev + dt*duprev + half_dtsq*k₂

  k₃ = f.f1(duprev,ku,p,t+dt)
  u = uprev + (dtsq/6)*(k₁ + 2*k₂) + dt*duprev
  du = duprev + (dt/6)*(k₁ + k₃ + 4*k₂)

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

@muladd function perform_step!(integrator,cache::Nystrom4VelocityIndependentCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,fsalfirst,k₂,k₃,k = cache
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  k₁ = integrator.fsalfirst.x[1]
  halfdt = dt/2
  dtsq = dt^2
  eighth_dtsq = dtsq/8
  half_dtsq = dtsq/2
  ttmp = t+halfdt

  ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ
  @.. broadcast=false ku = uprev + halfdt*duprev + eighth_dtsq*k₁

  f.f1(k₂,duprev,ku,p,ttmp)
  @.. broadcast=false ku = uprev + dt*duprev + half_dtsq*k₂

  f.f1(k₃,duprev,ku,p,t+dt)
  @.. broadcast=false u = uprev + (dtsq/6)*(k₁ + 2*k₂) + dt*duprev
  @.. broadcast=false du = duprev + (dt/6)*(k₁ + k₃ + 4*k₂)

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 3
  integrator.destats.nf2 += 1
end

@muladd function perform_step!(integrator,cache::IRKN3ConstantCache,repeat_step=false)
  @unpack t,dt,k,tprev,f,p = integrator
  duprev, uprev   = integrator.uprev.x
  duprev2,uprev2 = integrator.uprev2.x
  @unpack bconst1,bconst2,c1,a21,b1,b2,bbar1,bbar2 = cache
  k₁ = integrator.fsalfirst
  # if there's a discontinuity or the solver is in the first step
  if integrator.iter < 2 && !integrator.u_modified
    perform_step!(integrator,Nystrom4VelocityIndependentConstantCache())
    k = integrator.fsallast
    k1cache = ArrayParition((k.x[1],f.f1(duprev,uprev,p,t+c1*dt)))
    kdu = uprev + dt*(c1*duprev + dt*a21*k1cache.x[1])
    k₂.x[1] = f.f1(duprev,kdu,p,t+c1*dt)
    integrator.destats.nf += 2
  else
    kdu = uprev2 + dt*(c1*duprev2 + dt*a21*k1cache.x[1])
    ku  = uprev  + dt*(c1*duprev  + dt*a21*k1cache.x[2])

    k₂x1 = f.f1(duprev, ku, p, t+c1*dt)
    du = duprev + dt*(b1*k1cache.x[1] + bbar1*k1cache.x[1] + b2*(k₂x1-k₂.x[1]))
    u  = uprev + bconst1*dt*duprev + dt*(bconst2*duprev2 + dt*bbar2*(k₂x1-k₂.x[1]))

    integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
    integrator.destats.nf += 3
    integrator.destats.nf2 += 1
    copyto!(k₂.x[1],k₂.x[2])
    k1cache = ArrayPartition((k1cache.x[1],k.x[2]))
  end # end if
end

@muladd function perform_step!(integrator,cache::IRKN3Cache,repeat_step=false)
  @unpack t,dt,k,tprev,f,p = integrator
  du,u = integrator.u.x
  duprev, uprev   = integrator.uprev.x
  duprev2,uprev2 = integrator.uprev2.x
  uidx = eachindex(integrator.uprev.x[1])
  @unpack tmp,fsalfirst,k₂,k = cache
  @unpack bconst1,bconst2,c1,a21,b1,b2,bbar1,bbar2 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  k1cache = cache.tmp2
  k₁ = fsalfirst
  # if there's a discontinuity or the solver is in the first step
  if integrator.iter < 2 && !integrator.u_modified
    perform_step!(integrator,integrator.cache.onestep_cache)
    copyto!(k1cache.x[1], k.x[1])
    f.f1(k1cache.x[2],duprev,uprev,p,t+c1*dt)
    @.. broadcast=false kdu= uprev + dt*(c1*duprev + dt*a21*k1cache.x[2])
    f.f1(k₂.x[1],duprev,kdu,p,t+c1*dt)
    integrator.destats.nf += 2
  else
    @.. broadcast=false kdu = uprev2 + dt*(c1*duprev2 + dt*a21*k1cache.x[1])
    @.. broadcast=false ku  = uprev  + dt*(c1*duprev  + dt*a21*k1cache.x[2])

    f.f1(k₂.x[2], duprev, ku, p, t+c1*dt)
    @tight_loop_macros for i in uidx
      @inbounds u[i]  = uprev[i] + bconst1*dt*duprev[i] + dt*(bconst2*duprev2[i] + dt*bbar2*(k₂.x[2][i]-k₂.x[1][i]))
      @inbounds du[i] = duprev[i] + dt*(b1*k1cache.x[1][i] + bbar1*k1cache.x[2][i] + b2*(k₂.x[2][i]-k₂.x[1][i]))
    end
    f.f1(k.x[1],du,u,p,t+dt)
    f.f2(k.x[2],du,u,p,t+dt)
    integrator.destats.nf += 3
    integrator.destats.nf2 += 1
    copyto!(k₂.x[1],k₂.x[2])
    copyto!(k1cache.x[2],k1cache.x[1])
    copyto!(k1cache.x[1],k.x[1])
  end # end if
end

@muladd function perform_step!(integrator,cache::IRKN4Cache,repeat_step=false)
  @unpack t,dt,k,tprev,f,p = integrator
  du,u = integrator.u.x
  duprev, uprev   = integrator.uprev.x
  duprev2,uprev2 = integrator.uprev2.x
  uidx = eachindex(integrator.uprev.x[1])
  @unpack tmp,tmp2,fsalfirst,k₂,k₃,k = cache
  @unpack bconst1,bconst2,c1,c2,a21,a32,b1,b2,b3,bbar1,bbar2,bbar3 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  k1cache = integrator.cache.tmp2
  k₁ = fsalfirst
  # if there's a discontinuity or the solver is in the first step
  if integrator.iter < 2 && !integrator.u_modified
    perform_step!(integrator,integrator.cache.onestep_cache)
    copyto!(k1cache.x[1], k.x[1])
    f.f1(k1cache.x[2],duprev,uprev,p,t+c1*dt)
    @.. broadcast=false kdu= uprev + dt*(c1*duprev + dt*a21*k1cache.x[1])
    f.f1(k₂.x[1],duprev,kdu,p,t+c1*dt)
    @.. broadcast=false kdu= uprev + dt*(c2*duprev + dt*a32*k1cache.x[2])
    f.f1(k₃.x[1],duprev,kdu,p,t+c1*dt)
    integrator.destats.nf += 3
  else
    @.. broadcast=false ku  = uprev  + dt*(c1*duprev  + dt*a21*k1cache.x[1])
    @.. broadcast=false kdu = uprev2+ dt*(c1*duprev2 + dt*a21*k1cache.x[2])

    f.f1(k₂.x[2], duprev, ku, p, t+c1*dt)
    @.. broadcast=false ku  = uprev  + dt*(c2*duprev  + dt*a32*k₂.x[2])
    @.. broadcast=false kdu = uprev2 + dt*(c2*duprev2 + dt*a32*k₂.x[1])

    f.f1(k₃.x[2], duprev, ku, p, t+c2*dt)
    @tight_loop_macros for i in uidx
      @inbounds u[i]  = uprev[i] + dt*bconst1*duprev[i] + dt*(bconst2*duprev2[i] + dt*(bbar2*(k₂.x[2][i]-k₂.x[1][i]) + bbar3*(k₃.x[2][i]-k₃.x[1][i])))
      @inbounds du[i] = duprev[i] + dt*(b1*k1cache.x[1][i] + bbar1*k1cache.x[2][i] + b2*(k₂.x[2][i]-k₂.x[1][i]) + b3*(k₃.x[2][i]-k₃.x[1][i]))
    end
    f.f1(k.x[1],du,u,p,t+dt)
    f.f2(k.x[2],du,u,p,t+dt)
    integrator.destats.nf += 4
    integrator.destats.nf2 += 1
    copyto!(k₂.x[1],k₂.x[2])
    copyto!(k₃.x[1],k₃.x[2])
    copyto!(k1cache.x[2],k1cache.x[1])
    copyto!(k1cache.x[1],k.x[1])
  end # end if
end

@muladd function perform_step!(integrator,cache::Nystrom5VelocityIndependentConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, a21, a31, a32, a41, a42, a43, bbar1, bbar2, bbar3, b1, b2, b3, b4 = cache
  k₁ = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k₁)

  k₂ = f.f1(duprev,ku,p,t+c1*dt)
  ku = uprev + dt*(c2*duprev + dt*(a31*k₁ + a32*k₂))

  k₃ = f.f1(duprev,ku,p,t+c2*dt)
  ku = uprev + dt*(duprev + dt*(a41*k₁ + a42*k₂ + a43*k₃))

  k₄ = f.f1(duprev,ku,p,t+dt)
  u  = uprev + dt*(duprev + dt*(bbar1*k₁ + bbar2*k₂ + bbar3*k₃))
  du = duprev + dt*(b1*k₁ + b2*k₂ + b3*k₃ + b4*k₄)

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator,cache::Nystrom5VelocityIndependentCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  uidx = eachindex(integrator.uprev.x[1])
  @unpack tmp,fsalfirst,k₂,k₃,k₄,k = cache
  @unpack c1, c2, a21, a31, a32, a41, a42, a43, bbar1, bbar2, bbar3, b1, b2, b3, b4 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  k₁ = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k₁)

  f.f1(k₂,du,ku,p,t+c1*dt)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k₁ + a32*k₂))

  f.f1(k₃,du,ku,p,t+c2*dt)
  #@tight_loop_macros for i in uidx
  #  @inbounds ku[i] = uprev[i] + dt*(duprev[i] + dt*(a41*k₁[i] + a42*k₂[i] + a43*k₃[i]))
  #end
  @.. broadcast=false ku = uprev + dt*(duprev + dt*(a41*k₁ + a42*k₂ + a43*k₃))

  f.f1(k₄,duprev,ku,p,t+dt)
  #@tight_loop_macros for i in uidx
  #  @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(bbar1*k₁[i] + bbar2*k₂[i] + bbar3*k₃[i]))
  #  @inbounds du[i] = duprev[i] + dt*(b1*k₁[i] + b2*k₂[i] + b3*k₃[i] + b4*k₄[i])
  #end
  @.. broadcast=false u  = uprev + dt*(duprev + dt*(bbar1*k₁ + bbar2*k₂ + bbar3*k₃))
  @.. broadcast=false du = duprev + dt*(b1*k₁ + b2*k₂ + b3*k₃ + b4*k₄)
  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  return nothing
end

function initialize!(integrator, cache::DPRKN6ConstantCache)
  duprev,uprev = integrator.uprev.x
  integrator.kshortsize = 3
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)

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

  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast

end

@muladd function perform_step!(integrator,cache::DPRKN6ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, c3, c4, c5, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a63, a64, a65, b1, b3, b4, b5, bp1, bp3, bp4, bp5, bp6, btilde1, btilde2, btilde3, btilde4, btilde5, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6 = cache
  k1 = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k1)

  k2 = f.f1(duprev,ku,p,t+dt*c1)
  ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  k3 = f.f1(duprev,ku,p,t+dt*c2)
  ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  k4 = f.f1(duprev,ku,p,t+dt*c3)
  ku = uprev + dt*(c4*duprev + dt*(a51*k1 + a52*k2 + a53*k3 + a54*k4))

  k5 = f.f1(duprev,ku,p,t+dt*c4)
  ku = uprev + dt*(c5*duprev + dt*(a61*k1          + a63*k3 + a64*k4 + a65*k5)) # no a62

  k6 = f.f1(duprev,ku,p,t+dt*c5)
  u  = uprev + dt*(duprev + dt*(b1 *k1 + b3 *k3 + b4 *k4 + b5 *k5)) # b1 -- b5, no b2
  du = duprev                + dt*(bp1*k1 + bp3*k3 + bp4*k4 + bp5*k5 + bp6*k6) # bp1 -- bp6, no bp2

  #=
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(bhat1*k1.x[2][i] + bhat2*k2.x[2][i] + bhat3*k3.x[2][i]))
    @inbounds du[i] = duprev[i]+ dt*(bphat1*k1.x[2][i] + bphat3*k3.x[2][i] + bphat4*k4.x[2][i] + bphat5*k5.x[2][i] + bphat6*k6.x[2][i])
  end
  =#

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.k[1] = ArrayPartition(integrator.fsalfirst.x[1],k2)
  integrator.k[2] = ArrayPartition(k3, k4)
  integrator.k[3] = ArrayPartition(k5, k6)
  integrator.destats.nf += 6
  integrator.destats.nf2 += 1

  if integrator.opts.adaptive
    dtsq = dt^2
    uhat  = dtsq*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5)
    duhat = dt*(bptilde1*k1 + bptilde3*k3 + bptilde4*k4 + bptilde5*k5 + bptilde6*k6)
    utilde = ArrayPartition((duhat,uhat))
    atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

function initialize!(integrator, cache::DPRKN6Cache)
  @unpack fsalfirst,k = cache
  duprev,uprev = integrator.uprev.x

  integrator.fsalfirst = fsalfirst
  integrator.fsallast = k
  integrator.kshortsize = 3
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = ArrayPartition(cache.fsalfirst.x[1],cache.k2)
  integrator.k[2] = ArrayPartition(cache.k3, cache.k4)
  integrator.k[3] = ArrayPartition(cache.k5, cache.k6)
  integrator.f.f1(integrator.fsallast.x[1],duprev,uprev,integrator.p,integrator.t)
  integrator.f.f2(integrator.fsallast.x[2],duprev,uprev,integrator.p,integrator.t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
end

@muladd function perform_step!(integrator,cache::DPRKN6Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k5,k6,k,utilde = cache
  @unpack c1, c2, c3, c4, c5, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a63, a64, a65, b1, b3, b4, b5, bp1, bp3, bp4, bp5, bp6, btilde1, btilde2, btilde3, btilde4, btilde5, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  uidx = eachindex(integrator.uprev.x[2])
  k1 = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k1)

  f.f1(k2,du,ku,p,t+dt*c1)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  f.f1(k3,du,ku,p,t+dt*c2)
  @.. broadcast=false ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  f.f1(k4,du,ku,p,t+dt*c3)
  @.. broadcast=false ku = uprev + dt*(c4*duprev + dt*(a51*k1 + a52*k2 + a53*k3 + a54*k4))

  f.f1(k5,du,ku,p,t+dt*c4)
  @.. broadcast=false ku = uprev + dt*(c5*duprev + dt*(a61*k1 + a63*k3 + a64*k4 + a65*k5)) # no a62

  f.f1(k6,du,ku,p,t+dt*c5)

  @.. broadcast=false u  = uprev  + dt*(duprev + dt*(b1 *k1 + b3 *k3 + b4 *k4 + b5 *k5)) # b1 -- b5, no b2
  @.. broadcast=false du = duprev + dt*(bp1*k1 + bp3*k3 + bp4*k4 + bp5*k5 + bp6*k6) # bp1 -- bp6, no bp2

  #=
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(bhat1*k1.x[2][i] + bhat2*k2.x[2][i] + bhat3*k3.x[2][i]))
    @inbounds du[i] = duprev[i]+ dt*(bphat1*k1.x[2][i] + bphat3*k3.x[2][i] + bphat4*k4.x[2][i] + bphat5*k5.x[2][i] + bphat6*k6.x[2][i])
  end
  =#

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 6
  integrator.destats.nf2 += 1
  if integrator.opts.adaptive
    duhat, uhat = utilde.x
    dtsq = dt^2
    @.. broadcast=false uhat  = dtsq*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5)
    @.. broadcast=false duhat = dt*(bptilde1*k1 + bptilde3*k3 + bptilde4*k4 + bptilde5*k5 + bptilde6*k6)
    calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::DPRKN8ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, c3, c4, c5, c6, c7, c8, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, a81, a82, a83, a84, a85, a86, a87, a91, a93, a94, a95, a96, a97, b1, b3, b4, b5, b6, b7, bp1, bp3, bp4, bp5, bp6, bp7, bp8, btilde1, btilde3, btilde4, btilde5, btilde6, btilde7, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, bptilde8, bptilde9 = cache
  k1 = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k1)

  k2 = f.f1(duprev,ku,p,t+dt*c1)
  ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  k3 = f.f1(duprev,ku,p,t+dt*c2)
  ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  k4 = f.f1(duprev,ku,p,t+dt*c3)
  ku = uprev + dt*(c4*duprev + dt*(a51*k1 + a52*k2 + a53*k3 + a54*k4))

  k5 = f.f1(duprev,ku,p,t+dt*c4)
  ku = uprev + dt*(c5*duprev + dt*(a61*k1 + a62*k2 + a63*k3 + a64*k4 + a65*k5))

  k6 = f.f1(duprev,ku,p,t+dt*c5)
  ku = uprev + dt*(c6*duprev + dt*(a71*k1 + a72*k2 + a73*k3 + a74*k4 + a75*k5 + a76*k6))

  k7 = f.f1(duprev,ku,p,t+dt*c6)
  ku = uprev + dt*(c7*duprev + dt*(a81*k1 + a82*k2 + a83*k3 + a84*k4 + a85*k5 + a86*k6 + a87*k7))

  k8 = f.f1(duprev,ku,p,t+dt*c7)
  ku = uprev + dt*(c8*duprev + dt*(a91*k1 + a93*k3 + a94*k4 + a95*k5 + a96*k6 + a97*k7)) # no a92 & a98

  k9 = f.f1(duprev,ku,p,t+dt*c8)
  u  = uprev + dt*(duprev + dt*(b1 *k1 + b3 *k3 + b4 *k4 + b5 *k5 + b6 *k6 + b7 *k7)) # b1 -- b7, no b2
  du = duprev+ dt*(bp1*k1 + bp3*k3 + bp4*k4 + bp5*k5 + bp6*k6 + bp7*k7 + bp8*k8) # bp1 -- bp8, no bp2

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 9
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast

  if integrator.opts.adaptive
    dtsq = dt^2
    uhat  = dtsq*(btilde1*k1 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7)
    duhat = dt*(bptilde1*k1 + bptilde3*k3 + bptilde4*k4 + bptilde5*k5 + bptilde6*k6 + bptilde7*k7 + bptilde8*k8 + bptilde9*k9)
    utilde = ArrayPartition((duhat,uhat))
    atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::DPRKN8Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k5,k6,k7,k8,k9,k,utilde = cache
  @unpack c1, c2, c3, c4, c5, c6, c7, c8, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, a81, a82, a83, a84, a85, a86, a87, a91, a93, a94, a95, a96, a97, b1, b3, b4, b5, b6, b7, bp1, bp3, bp4, bp5, bp6, bp7, bp8, btilde1, btilde3, btilde4, btilde5, btilde6, btilde7, bptilde1, bptilde3, bptilde4, bptilde5, bptilde6, bptilde7, bptilde8, bptilde9 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  uidx = eachindex(integrator.uprev.x[2])
  k1 = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k1)

  f.f1(k2,duprev,ku,p,t+dt*c1)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  f.f1(k3,duprev,ku,p,t+dt*c2)
  @.. broadcast=false ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  f.f1(k4,duprev,ku,p,t+dt*c3)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c4*duprev[i] + dt*(a51*k1[i] + a52*k2[i] + a53*k3[i] + a54*k4[i]))
  end

  f.f1(k5,duprev,ku,p,t+dt*c4)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c5*duprev[i] + dt*(a61*k1[i] + a62*k2[i] + a63*k3[i] + a64*k4[i] + a65*k5[i]))
  end

  f.f1(k6,duprev,ku,p,t+dt*c5)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c6*duprev[i] + dt*(a71*k1[i] + a72*k2[i] + a73*k3[i] + a74*k4[i] + a75*k5[i] + a76*k6[i]))
  end

  f.f1(k7,duprev,ku,p,t+dt*c6)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c7*duprev[i] + dt*(a81*k1[i] + a82*k2[i] + a83*k3[i] + a84*k4[i] + a85*k5[i] + a86*k6[i] + a87*k7[i]))
  end

  f.f1(k8,duprev,ku,p,t+dt*c7)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c8*duprev[i] + dt*(a91*k1[i] + a93*k3[i] + a94*k4[i] + a95*k5[i] + a96*k6[i] + a97*k7[i])) # no a92 & a98
  end

  f.f1(k9,duprev,ku,p,t+dt*c8)
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b3 *k3[i] + b4 *k4[i] + b5 *k5[i] + b6 *k6[i] + b7 *k7[i])) # b1 -- b7, no b2
    @inbounds du[i] = duprev[i]+ dt*(bp1*k1[i] + bp3*k3[i] + bp4*k4[i] + bp5*k5[i] + bp6*k6[i] + bp7*k7[i] + bp8*k8[i]) # bp1 -- bp8, no bp2
  end

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 9
  integrator.destats.nf2 += 1
  if integrator.opts.adaptive
    duhat, uhat = utilde.x
    dtsq = dt^2
    @tight_loop_macros for i in uidx
      @inbounds uhat[i]  = dtsq*(btilde1*k1[i] + btilde3*k3[i] + btilde4*k4[i] + btilde5*k5[i] + btilde6*k6[i] + btilde7*k7[i])
      @inbounds duhat[i] = dt*(bptilde1*k1[i] + bptilde3*k3[i] + bptilde4*k4[i] + bptilde5*k5[i] + bptilde6*k6[i] + bptilde7*k7[i] + bptilde8*k8[i] + bptilde9*k9[i])
    end
    calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::DPRKN12ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14, c15, c16, a21, a31, a32, a41, a42, a43, a51, a53, a54, a61, a63, a64, a65, a71, a73, a74, a75, a76, a81, a84, a85, a86, a87, a91, a93, a94, a95, a96, a97, a98, a101, a103, a104, a105, a106, a107, a108, a109, a111, a113, a114, a115, a116, a117, a118, a119, a1110, a121, a123, a124, a125, a126, a127, a128, a129, a1210, a1211, a131, a133, a134, a135, a136, a137, a138, a139, a1310, a1311, a1312, a141, a143, a144, a145, a146, a147, a148, a149, a1410, a1411, a1412, a1413, a151, a153, a154, a155, a156, a157, a158, a159, a1510, a1511, a1512, a1513, a1514, a161, a163, a164, a165, a166, a167, a168, a169, a1610, a1611, a1612, a1613, a1614, a1615, a171, a173, a174, a175, a176, a177, a178, a179, a1710, a1711, a1712, a1713, a1714, a1715, b1, b7, b8, b9, b10, b11, b12, b13, b14, b15, bp1, bp7, bp8, bp9, bp10, bp11, bp12, bp13, bp14, bp15, bp16, bp17, btilde1, btilde7, btilde8, btilde9, btilde10, btilde11, btilde12, btilde13, btilde14, btilde15, bptilde1, bptilde7, bptilde8, bptilde9, bptilde10, bptilde11, bptilde12, bptilde13, bptilde14, bptilde15, bptilde16, bptilde17 = cache
  k1 = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k1)

  k2 = f.f1(duprev,ku,p,t+dt*c1)
  ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  k3 = f.f1(duprev,ku,p,t+dt*c2)
  ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  k4 = f.f1(duprev,ku,p,t+dt*c3)
  ku = uprev + dt*(c4*duprev + dt*(a51*k1 + a53*k3 + a54*k4)) # no a52

  k5 = f.f1(duprev,ku,p,t+dt*c4)
  ku = uprev + dt*(c5*duprev + dt*(a61*k1 + a63*k3 + a64*k4 + a65*k5)) # no a62

  k6 = f.f1(duprev,ku,p,t+dt*c5)
  ku = uprev + dt*(c6*duprev + dt*(a71*k1 + a73*k3 + a74*k4 + a75*k5 + a76*k6)) # no a72

  k7 = f.f1(duprev,ku,p,t+dt*c6)
  ku = uprev + dt*(c7*duprev + dt*(a81*k1 + a84*k4 + a85*k5 + a86*k6 + a87*k7)) # no a82, a83

  k8 = f.f1(duprev,ku,p,t+dt*c7)
  ku = uprev + dt*(c8*duprev + dt*(a91*k1 + a93*k3 + a94*k4 + a95*k5 + a96*k6 + a97*k7 + a98*k8)) # no a92

  k9 = f.f1(duprev,ku,p,t+dt*c8)
  ku = uprev + dt*(c9*duprev + dt*(a101*k1 + a103*k3 + a104*k4 + a105*k5 + a106*k6 + a107*k7 + a108*k8 + a109*k9)) # no a102

  k10 = f.f1(duprev,ku,p,t+dt*c9)
  ku = uprev + dt*(c10*duprev + dt*(a111*k1 + a113*k3 + a114*k4 + a115*k5 + a116*k6 + a117*k7 + a118*k8 + a119*k9 + a1110*k10)) # no a112

  k11 = f.f1(duprev,ku,p,t+dt*c10)
  ku = uprev + dt*(c11*duprev + dt*(a121*k1 + a123*k3 + a124*k4 + a125*k5 + a126*k6 + a127*k7 + a128*k8 + a129*k9 + a1210*k10 + a1211*k11)) # no a122

  k12 = f.f1(duprev,ku,p,t+dt*c11)
  ku = uprev + dt*(c12*duprev + dt*(a131*k1 + a133*k3 + a134*k4 + a135*k5 + a136*k6 + a137*k7 + a138*k8 + a139*k9 + a1310*k10 + a1311*k11 + a1312*k12)) # no a132

  k13 = f.f1(duprev,ku,p,t+dt*c12)
  ku = uprev + dt*(c13*duprev + dt*(a141*k1 + a143*k3 + a144*k4 + a145*k5 + a146*k6 + a147*k7 + a148*k8 + a149*k9 + a1410*k10 + a1411*k11 + a1412*k12 + a1413*k13)) # no a142

  k14 = f.f1(duprev,ku,p,t+dt*c13)
  ku = uprev + dt*(c14*duprev + dt*(a151*k1 + a153*k3 + a154*k4 + a155*k5 + a156*k6 + a157*k7 + a158*k8 + a159*k9 + a1510*k10 + a1511*k11 + a1512*k12 + a1513*k13 + a1514*k14)) # no a152

  k15 = f.f1(duprev,ku,p,t+dt*c14)
  ku = uprev + dt*(c15*duprev + dt*(a161*k1 + a163*k3 + a164*k4 + a165*k5 + a166*k6 + a167*k7 + a168*k8 + a169*k9 + a1610*k10 + a1611*k11 + a1612*k12 + a1613*k13 + a1614*k14 + a1615*k15)) # no a162

  k16 = f.f1(duprev,ku,p,t+dt*c15)
  ku = uprev + dt*(c16*duprev + dt*(a171*k1 + a173*k3 + a174*k4 + a175*k5 + a176*k6 + a177*k7 + a178*k8 + a179*k9 + a1710*k10 + a1711*k11 + a1712*k12 + a1713*k13 + a1714*k14 + a1715*k15)) # no a172, a1716

  k17 = f.f1(duprev,ku,p,t+dt*c16)
  u  = uprev + dt*(duprev + dt*(b1 *k1 + b7 *k7+ b8 *k8+ b9 *k9+ b10 *k10+ b11 *k11+ b12 *k12+ b13 *k13+ b14 *k14+ b15 *k15)) # b1 & b7 -- b15
  du = duprev+ dt*(bp1*k1 + bp7*k7+ bp8*k8+ bp9*k9+ bp10*k10+ bp11*k11+ bp12*k12+ bp13*k13+ bp14*k14+ bp15*k15+ bp16*k16+ bp17*k17) # bp1 & bp7 -- bp17

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 17
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  if integrator.opts.adaptive
    dtsq = dt^2
    uhat  = dtsq*(btilde1*k1 + btilde7*k7 + btilde8*k8 + btilde9*k9 + btilde10*k10 + btilde11*k11 + btilde12*k12 + btilde13*k13 + btilde14*k14 + btilde15*k15) # btilde1 & btilde7 -- btilde15
    duhat = dt*(bptilde1*k1 + bptilde7*k7 + bptilde8*k8 + bptilde9*k9 + bptilde10*k10 + bptilde11*k11 + bptilde12*k12 + bptilde13*k13 + bptilde14*k14 + bptilde15*k15 + bptilde16*k16 + bptilde17*k17) # bptilde1 & bptilde7 -- bptilde17
    utilde = ArrayPartition((duhat,uhat))
    atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::DPRKN12Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k,utilde = cache
  @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14, c15, c16, a21, a31, a32, a41, a42, a43, a51, a53, a54, a61, a63, a64, a65, a71, a73, a74, a75, a76, a81, a84, a85, a86, a87, a91, a93, a94, a95, a96, a97, a98, a101, a103, a104, a105, a106, a107, a108, a109, a111, a113, a114, a115, a116, a117, a118, a119, a1110, a121, a123, a124, a125, a126, a127, a128, a129, a1210, a1211, a131, a133, a134, a135, a136, a137, a138, a139, a1310, a1311, a1312, a141, a143, a144, a145, a146, a147, a148, a149, a1410, a1411, a1412, a1413, a151, a153, a154, a155, a156, a157, a158, a159, a1510, a1511, a1512, a1513, a1514, a161, a163, a164, a165, a166, a167, a168, a169, a1610, a1611, a1612, a1613, a1614, a1615, a171, a173, a174, a175, a176, a177, a178, a179, a1710, a1711, a1712, a1713, a1714, a1715, b1, b7, b8, b9, b10, b11, b12, b13, b14, b15, bp1, bp7, bp8, bp9, bp10, bp11, bp12, bp13, bp14, bp15, bp16, bp17, btilde1, btilde7, btilde8, btilde9, btilde10, btilde11, btilde12, btilde13, btilde14, btilde15, bptilde1, bptilde7, bptilde8, bptilde9, bptilde10, bptilde11, bptilde12, bptilde13, bptilde14, bptilde15, bptilde16, bptilde17 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  uidx = eachindex(integrator.uprev.x[2])
  k1 = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k1)

  f.f1(k2,duprev,ku,p,t+dt*c1)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  f.f1(k3,duprev,ku,p,t+dt*c2)
  @.. broadcast=false ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  f.f1(k4,duprev,ku,p,t+dt*c3)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c4*duprev[i] + dt*(a51*k1[i] + a53*k3[i] + a54*k4[i])) # no a52
  end

  f.f1(k5,duprev,ku,p,t+dt*c4)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c5*duprev[i] + dt*(a61*k1[i] + a63*k3[i] + a64*k4[i] + a65*k5[i])) # no a62
  end

  f.f1(k6,duprev,ku,p,t+dt*c5)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c6*duprev[i] + dt*(a71*k1[i] + a73*k3[i] + a74*k4[i] + a75*k5[i] + a76*k6[i])) # no a72
  end

  f.f1(k7,duprev,ku,p,t+dt*c6)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c7*duprev[i] + dt*(a81*k1[i] + a84*k4[i] + a85*k5[i] + a86*k6[i] + a87*k7[i])) # no a82, a83
  end

  f.f1(k8,duprev,ku,p,t+dt*c7)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c8*duprev[i] + dt*(a91*k1[i] + a93*k3[i] + a94*k4[i] + a95*k5[i] + a96*k6[i] + a97*k7[i] + a98*k8[i])) # no a92
  end

  f.f1(k9,duprev,ku,p,t+dt*c8)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c9*duprev[i] + dt*(a101*k1[i] + a103*k3[i] + a104*k4[i] + a105*k5[i] + a106*k6[i] + a107*k7[i] + a108*k8[i] + a109*k9[i])) # no a102
  end

  f.f1(k10,duprev,ku,p,t+dt*c9)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c10*duprev[i] + dt*(a111*k1[i] + a113*k3[i] + a114*k4[i] + a115*k5[i] + a116*k6[i] + a117*k7[i] + a118*k8[i] + a119*k9[i] + a1110*k10[i])) # no a112
  end

  f.f1(k11,duprev,ku,p,t+dt*c10)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c11*duprev[i] + dt*(a121*k1[i] + a123*k3[i] + a124*k4[i] + a125*k5[i] + a126*k6[i] + a127*k7[i] + a128*k8[i] + a129*k9[i] + a1210*k10[i] + a1211*k11[i])) # no a122
  end

  f.f1(k12,duprev,ku,p,t+dt*c11)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c12*duprev[i] + dt*(a131*k1[i] + a133*k3[i] + a134*k4[i] + a135*k5[i] + a136*k6[i] + a137*k7[i] + a138*k8[i] + a139*k9[i] + a1310*k10[i] + a1311*k11[i] + a1312*k12[i])) # no a132
  end

  f.f1(k13,duprev,ku,p,t+dt*c12)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c13*duprev[i] + dt*(a141*k1[i] + a143*k3[i] + a144*k4[i] + a145*k5[i] + a146*k6[i] + a147*k7[i] + a148*k8[i] + a149*k9[i] + a1410*k10[i] + a1411*k11[i] + a1412*k12[i] + a1413*k13[i])) # no a142
  end

  f.f1(k14,duprev,ku,p,t+dt*c13)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c14*duprev[i] + dt*(a151*k1[i] + a153*k3[i] + a154*k4[i] + a155*k5[i] + a156*k6[i] + a157*k7[i] + a158*k8[i] + a159*k9[i] + a1510*k10[i] + a1511*k11[i] + a1512*k12[i] + a1513*k13[i] + a1514*k14[i])) # no a152
  end

  f.f1(k15,duprev,ku,p,t+dt*c14)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c15*duprev[i] + dt*(a161*k1[i] + a163*k3[i] + a164*k4[i] + a165*k5[i] + a166*k6[i] + a167*k7[i] + a168*k8[i] + a169*k9[i] + a1610*k10[i] + a1611*k11[i] + a1612*k12[i] + a1613*k13[i] + a1614*k14[i] + a1615*k15[i])) # no a162
  end

  f.f1(k16,duprev,ku,p,t+dt*c15)
  @tight_loop_macros for i in uidx
    @inbounds ku[i] = uprev[i] + dt*(c16*duprev[i] + dt*(a171*k1[i] + a173*k3[i] + a174*k4[i] + a175*k5[i] + a176*k6[i] + a177*k7[i] + a178*k8[i] + a179*k9[i] + a1710*k10[i] + a1711*k11[i] + a1712*k12[i] + a1713*k13[i] + a1714*k14[i] + a1715*k15[i])) # no a172, a1716
  end

  f.f1(k17,duprev,ku,p,t+dt*c16)
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b7 *k7[i]+ b8 *k8[i]+ b9 *k9[i]+ b10 *k10[i]+ b11 *k11[i]+ b12 *k12[i]+ b13 *k13[i]+ b14 *k14[i]+ b15 *k15[i])) # b1 & b7 -- b15
    @inbounds du[i] = duprev[i]+ dt*(bp1*k1[i] + bp7*k7[i]+ bp8*k8[i]+ bp9*k9[i]+ bp10*k10[i]+ bp11*k11[i]+ bp12*k12[i]+ bp13*k13[i]+ bp14*k14[i]+ bp15*k15[i]+ bp16*k16[i]+ bp17*k17[i]) # bp1 & bp7 -- bp17
  end

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 17
  integrator.destats.nf2 += 1
  if integrator.opts.adaptive
    duhat, uhat = utilde.x
    dtsq = dt^2
    @tight_loop_macros for i in uidx
      @inbounds uhat[i]  = dtsq*(btilde1*k1[i] + btilde7*k7[i] + btilde8*k8[i] + btilde9*k9[i] + btilde10*k10[i] + btilde11*k11[i] + btilde12*k12[i] + btilde13*k13[i] + btilde14*k14[i] + btilde15*k15[i]) # btilde1 & btilde7 -- btilde15
      @inbounds duhat[i] = dt*(bptilde1*k1[i] + bptilde7*k7[i] + bptilde8*k8[i] + bptilde9*k9[i] + bptilde10*k10[i] + bptilde11*k11[i] + bptilde12*k12[i] + bptilde13*k13[i] + bptilde14*k14[i] + bptilde15*k15[i] + bptilde16*k16[i] + bptilde17*k17[i]) # bptilde1 & bptilde7 -- bptilde17
    end
    calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::ERKN4ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4, bptilde1, bptilde2, bptilde3, bptilde4 = cache
  k1 = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k1)

  k2 = f.f1(duprev,ku,p,t+dt*c1)
  ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  k3 = f.f1(duprev,ku,p,t+dt*c2)
  ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  k4 = f.f1(duprev,ku,p,t+dt*c3)
  u  = uprev + dt*(duprev + dt*(b1 *k1 + b2 *k2 + b3 *k3 + b4 *k4))
  du = duprev                + dt*(bp1*k1 + bp2*k2 + bp3*k3 + bp4*k4)

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast

  if integrator.opts.adaptive
    dtsq = dt^2
    uhat  = dtsq*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
    duhat = dt*(bptilde1*k1 + bptilde2*k2 + bptilde3*k3 + bptilde4*k4)
    utilde = ArrayPartition((duhat,uhat))
    atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::ERKN4Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k,utilde = cache
  @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4, bptilde1, bptilde2, bptilde3, bptilde4 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  uidx = eachindex(integrator.uprev.x[2])
  k1 = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k1)

  f.f1(k2,duprev,ku,p,t+dt*c1)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  f.f1(k3,duprev,ku,p,t+dt*c2)
  @.. broadcast=false ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  f.f1(k4,duprev,ku,p,t+dt*c3)
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b2 *k2[i] + b3 *k3[i] + b4 *k4[i]))
    @inbounds du[i] = duprev[i]                + dt*(bp1*k1[i] + bp2*k2[i] + bp3*k3[i] + bp4*k4[i])
  end

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  if integrator.opts.adaptive
    duhat, uhat = utilde.x
    dtsq = dt^2
    @tight_loop_macros for i in uidx
      @inbounds uhat[i]  = dtsq*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i])
      @inbounds duhat[i] = dt*(bptilde1*k1[i] + bptilde2*k2[i] + bptilde3*k3[i] + bptilde4*k4[i])
    end
    calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::ERKN5ConstantCache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  duprev,uprev = integrator.uprev.x
  @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4 = cache
  k1 = integrator.fsalfirst.x[1]

  ku = uprev + dt*(c1*duprev + dt*a21*k1)

  k2 = f.f1(duprev,ku,p,t+dt*c1)
  ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  k3 = f.f1(duprev,ku,p,t+dt*c2)
  ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  k4 = f.f1(duprev,ku,p,t+dt*c3)
  u  = uprev + dt*(duprev + dt*(b1 *k1 + b2 *k2 + b3 *k3 + b4 *k4))
  du = duprev                + dt*(bp1*k1 + bp2*k2 + bp3*k3 + bp4*k4)

  integrator.u = ArrayPartition((du,u))
  integrator.fsallast = ArrayPartition((f.f1(du,u,p,t+dt),f.f2(du,u,p,t+dt)))
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  if integrator.opts.adaptive
    dtsq = dt^2
    uhat = dtsq*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4)
    atmp = calculate_residuals(uhat, integrator.uprev.x[2], integrator.u.x[2], integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
end

@muladd function perform_step!(integrator,cache::ERKN5Cache,repeat_step=false)
  @unpack t,dt,f,p = integrator
  du,u = integrator.u.x
  duprev,uprev = integrator.uprev.x
  @unpack tmp,atmp,fsalfirst,k2,k3,k4,k,utilde = cache
  @unpack c1, c2, c3, a21, a31, a32, a41, a42, a43, b1, b2, b3, b4, bp1, bp2, bp3, bp4, btilde1, btilde2, btilde3, btilde4 = cache.tab
  kdu,ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2]
  uidx = eachindex(integrator.uprev.x[2])
  k1 = integrator.fsalfirst.x[1]

  @.. broadcast=false ku = uprev + dt*(c1*duprev + dt*a21*k1)

  f.f1(k2,duprev,ku,p,t+dt*c1)
  @.. broadcast=false ku = uprev + dt*(c2*duprev + dt*(a31*k1 + a32*k2))

  f.f1(k3,duprev,ku,p,t+dt*c2)
  @.. broadcast=false ku = uprev + dt*(c3*duprev + dt*(a41*k1 + a42*k2 + a43*k3))

  f.f1(k4,duprev,ku,p,t+dt*c3)
  @tight_loop_macros for i in uidx
    @inbounds u[i]  = uprev[i] + dt*(duprev[i] + dt*(b1 *k1[i] + b2 *k2[i] + b3 *k3[i] + b4 *k4[i]))
    @inbounds du[i] = duprev[i]                + dt*(bp1*k1[i] + bp2*k2[i] + bp3*k3[i] + bp4*k4[i])
  end

  f.f1(k.x[1],du,u,p,t+dt)
  f.f2(k.x[2],du,u,p,t+dt)
  integrator.destats.nf += 4
  integrator.destats.nf2 += 1
  if integrator.opts.adaptive
    duhat, uhat = utilde.x
    dtsq = dt^2
    @tight_loop_macros for i in uidx
      @inbounds uhat[i]  = dtsq*(btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i])
    end
    calculate_residuals!(atmp.x[2], uhat, integrator.uprev.x[2], integrator.u.x[2], integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp.x[2],t)
  end
end


text:
@inline function add_fpfp_db(x::T, y::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(add_fpfp_dd(x, y))
end

@inline function sub_fpfp_db(x::T, y::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(sub_fpfp_dd(x, y))
end

@inline function mul_fpfp_db(x::T, y::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(mul_fpfp_dd(x, y))
end

@inline function dvi_fpfp_db(x::T, y::T) where {T<:IEEEFloat}
    return DoubleFloat{T}(dvi_fpfp_dd(x, y))
end

ext: jl


text:
# error analysis for ISDF in BSE

# %%
# loading packages
using BenchmarkTools, JLD2, FileIO, LinearAlgebra, FFTW, Statistics, Arpack
using BSE_k_ISDF

BLAS.set_num_threads(1)
FFTW.set_num_threads(1)

# general problem setup
isdir("1d") || mkdir("1d")

l = 1.5
V_sp = r -> 20 * cos(4π / l * r) + 0.2 * sin(2π / l * r)

V_func = (r_1, r_2) -> 1 /  sqrt((r_1 - r_2)^2 + 0.01)
W_func = (r_1, r_2, l, L) -> 0.0625 * (3.0 + sin(2π / l * r_1)) * (3.0 + cos(4π / l * r_2)) *
    exp(-abs(BSE_k_ISDF.supercell_difference(r_1, r_2, L))^2 / (2 * (4 * l)^2)) * V_func(BSE_k_ISDF.supercell_difference(r_1, r_2, L), 0.)

# fixed parameters
N_r = 128

N_core = 0
N_v = 4
N_c = 5
N_k = 256

# set up problem
sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_r, N_k);
prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func);

u_v, u_c = prob.u_v, prob.u_c

V_grid = V_func.(prob.prob.r_super, 0)
W_grid = W_func.(prob.prob.r_super, 0, l, N_k * l)
save("1d/example_$(N_r)_$(N_v)_$(N_c)_$(N_k).jld2", "r_super", prob.prob.r_super, "V_grid", V_grid, "W_grid", W_grid, "E_v", prob.E_v, "E_c", prob.E_c, "k_bz", prob.prob.k_bz)

# variable parameters
N_μ_vec = 1:40

N_sub = 30
F_vv = BSE_k_ISDF.qrcp(u_v, N_sub);
F_cc = BSE_k_ISDF.qrcp(u_c, N_sub);
F_vc = BSE_k_ISDF.qrcp(u_v, u_c, N_sub^2);

# for plotting
rc = BSE_k_ISDF.lattice_matrix(prob) * BSE_k_ISDF.r_lattice(prob)
r_μ_vv = rc[:, F_vv.p];
r_μ_cc = rc[:, F_cc.p];
r_μ_vc = rc[:, F_vc.p];
save("1d/real_space_grid_$(N_sub).jld2", "rc", rc, "r_μ_vv", r_μ_vv, "r_μ_cc", r_μ_cc, "r_μ_vc", r_μ_vc)

errors_Z_vv = []
errors_Z_cc = []
errors_Z_vc = []

@time for N_μ_vv in N_μ_vec
    r_μ_vv_indices = F_vv.p[1:N_μ_vv]
    ζ_vv = BSE_k_ISDF.assemble_ζ(u_v, r_μ_vv_indices)
    # error_Z_vv = BSE_k_ISDF.isdf_error_estimate(u_v, ζ_vv, r_μ_vv_indices, N_k_samples)
    error_Z_vv = BSE_k_ISDF.isdf_error(u_v, ζ_vv, r_μ_vv_indices)

    push!(errors_Z_vv, error_Z_vv)
end
@time for N_μ_cc in N_μ_vec
    r_μ_cc_indices = F_cc.p[1:N_μ_cc]
    ζ_cc = BSE_k_ISDF.assemble_ζ(u_c, r_μ_cc_indices)
    # error_Z_cc = BSE_k_ISDF.isdf_error_estimate(u_c, ζ_cc, r_μ_cc_indices, N_k_samples)
    error_Z_cc = BSE_k_ISDF.isdf_error(u_c, ζ_cc, r_μ_cc_indices)

    push!(errors_Z_cc, error_Z_cc)
end
@time for N_μ_vc in N_μ_vec
    r_μ_vc_indices = F_vc.p[1:N_μ_vc]
    ζ_vc = BSE_k_ISDF.assemble_ζ(u_v, u_c, r_μ_vc_indices)
    # error_Z_vc = BSE_k_ISDF.isdf_error_estimate(u_v, u_c, ζ_vc, r_μ_vc_indices, N_k_samples)
    error_Z_vc = BSE_k_ISDF.isdf_error(u_v, u_c, ζ_vc, r_μ_vc_indices)

    push!(errors_Z_vc, error_Z_vc)
end

save("1d/errors_Z.jld2", "N_μ_vec", N_μ_vec, "errors_Z_vv", errors_Z_vv, "errors_Z_cc", errors_Z_cc, "errors_Z_vc", errors_Z_vc)

# %% set up reference Hamiltonian

# fixed parameters
N_unit = 128

N_core = 0
N_v = 4
N_c = 5
N_k = 256

# set up problem
sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_unit, N_k);
prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func);

# compute reference Hamiltonian
t_V_entry = @benchmark BSE_k_ISDF.V_entry($prob.v_hat, 1, 1, 15, 1, 1, 60, $prob.u_v, $prob.u_c, $l, $(BSE_k_ISDF.size_r(prob)), $(BSE_k_ISDF.size_k(prob)))
w_hat, q_2bz_ind, q_2bz_shift = BSE_k_ISDF.compute_w_hat(prob);
ikkp2iq = BSE_k_ISDF.ikkp2iq_matrix(BSE_k_ISDF.k_lattice(prob), BSE_k_ISDF.q_lattice(prob));
t_W_entry = @benchmark BSE_k_ISDF.W_entry($w_hat, 1, 1, 15, 1, 1, 60, $prob.u_v, $prob.u_c, $l, $(BSE_k_ISDF.size_r(prob)), $ikkp2iq, $q_2bz_ind, $q_2bz_shift)

println("estimated time to assemble H_exact is ", (mean(t_V_entry).time + mean(t_W_entry).time) * (N_v * N_c * N_k)^2 * 1e-9, " seconds")
sleep(1) # to make sure the estimated time is printed

D = BSE_k_ISDF.setup_D(prob);
V_exact = BSE_k_ISDF.assemble_exact_V(prob)
W_exact = BSE_k_ISDF.assemble_exact_W(prob)
H_exact = D + 2 * V_exact - W_exact

@time F = eigen(H_exact)

save("1d/H_exact_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "H_exact", H_exact, "ev", F.values, "ef", F.vectors)

# %% copute error in H for different errors in M (it might not be necessary to compute this)

N_unit = 128

N_core = 0
N_v = 4
N_c = 5
N_k = 256

H_exact = load("1d/H_exact_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "H_exact")

N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc = load("1d/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

N_μ_vv_vec = []
N_μ_cc_vec = []
N_μ_vc_vec = []
errors_H = []

Z_tol_vec = [0.8, 0.4, 0.2, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
for Z_tol in Z_tol_vec
    N_μ_vv = N_μ_vec[findfirst(errors_Z_vv .<= Z_tol)]
    N_μ_cc = N_μ_vec[findfirst(errors_Z_cc .<= Z_tol)]
    N_μ_vc = N_μ_vec[findfirst(errors_Z_vc .<= Z_tol)]

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)

    H = BSE_k_ISDF.setup_H(prob, isdf)

    @time H_full = Matrix(H)

    error_H = norm(H_full - H_exact) / norm(H_exact)

    push!(N_μ_vv_vec, N_μ_vv)
    push!(N_μ_cc_vec, N_μ_cc)
    push!(N_μ_vc_vec, N_μ_vc)
    push!(errors_H, error_H)
end

save("1d/errors_H_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "Z_tol_vec", Z_tol_vec, "N_μ_vv_vec", N_μ_vv_vec, "N_μ_cc_vec", N_μ_cc_vec, "N_μ_vc_vec", N_μ_vc_vec, "errors_H", errors_H)

# %% compute reference absorption spectrum

# fixed parameters
N_unit = 128

N_core = 0
N_v = 4
N_c = 5
N_k = 256

# broadening
σ = 0.1
g = (ω, σ) -> 1 / π * σ / (ω^2 + σ^2)
# energy range
E_min, E_max = 0.0, 20.0
Erange = E_min:0.01:E_max
# parameter for lanczos
N_iter = 200

# set up problem
sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_unit, N_k)
prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func)

H_exact, ev, ef = load("1d/H_exact_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "H_exact", "ev", "ef")

d = BSE_k_ISDF.optical_absorption_vector(prob, 1)
weights = [abs2(ef[:, n]' * d) for n in 1:size(H_exact, 1)]
optical_absorption = BSE_k_ISDF.lanczos_optical_absorption(real.(ev), weights, ω -> g(ω, σ), Erange, 8 * π^2 / (l * N_k))
optical_absorption_lanc = BSE_k_ISDF.lanczos_optical_absorption(H_exact, normalize(d), N_iter, ω -> g(ω, σ), Erange, norm(d)^2 * 8 * π^2 / (l * N_k))

save("1d/optical_absorption_ref_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "Erange", Erange, "optical_absorption", optical_absorption)
save("1d/optical_absorption_lanc_ref_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(N_iter).jld2", "Erange", Erange, "optical_absorption_lanc", optical_absorption_lanc)

# %% compute absorption spectra for different k

# fixed parameters
N_unit = 128

N_core = 0
N_v = 4
N_c = 5

N_μ_cc = 22
N_μ_vv = 16
N_μ_vc = 21

# broadening
σ = 0.1
g = (ω, σ) -> 1 / π * σ / (ω^2 + σ^2)
# energy range
E_min, E_max = 0.0, 20.0
Erange = E_min:0.01:E_max
# parameter for lanczos
N_iter = 200

#variable parameters
N_k_vec = 2 .^(4:10)

for N_k in N_k_vec
    sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_unit, N_k)
    prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func)

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)

    H = BSE_k_ISDF.setup_H(prob, isdf)

    optical_absorption_lanc = BSE_k_ISDF.lanczos_optical_absorption(prob, isdf, 1, N_iter, ω -> g(ω, σ), Erange)
    ev, ef = eigs(H, which=:SR, nev = 1, maxiter=1000)

    # save results
    save("1d/optical_absorption_lanczos_$(N_unit)_$(N_k)_$(N_iter).jld2", "Erange", Erange, "optical_absorption_lanc", optical_absorption_lanc)

    save("1d/eigs_$(N_unit)_$(N_k).jld2", "ev", ev, "ef", ef)
end

# %% compute absorption spectra for different N_μ

# fixed parameters
N_unit = 128

N_core = 0
N_v = 4
N_c = 5
N_k = 256

sp_prob = BSE_k_ISDF.SPProblem1D(V_sp, l, N_unit, N_k)
prob = BSE_k_ISDF.BSEProblem1D(sp_prob, N_core, N_v, N_c, V_func, W_func)

N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc = load("1d/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

# broadening
σ = 0.1
g = (ω, σ) -> 1 / π * σ / (ω^2 + σ^2)
# energy range
E_min, E_max = 0.0, 20.0
Erange = E_min:0.01:E_max
# parameter for lanczos
N_iter = 200

#variable parameters
Z_tol_vec = [0.8, 0.4, 0.2, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
for Z_tol in Z_tol_vec
    N_μ_vv = N_μ_vec[findfirst(errors_Z_vv .<= Z_tol)]
    N_μ_cc = N_μ_vec[findfirst(errors_Z_cc .<= Z_tol)]
    N_μ_vc = N_μ_vec[findfirst(errors_Z_vc .<= Z_tol)]
#     @show N_μ_vv, N_μ_cc, N_μ_vc
# end

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)

    H = BSE_k_ISDF.setup_H(prob, isdf)

    optical_absorption_lanc = BSE_k_ISDF.lanczos_optical_absorption(prob, isdf, 1, N_iter, ω -> g(ω, σ), Erange)

    ev, ef = eigs(H, which=:SR, nev = 1, maxiter=1000)

    # save results
    save("1d/optical_absorption_lanczos_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(N_iter)_$(Z_tol).jld2", "Erange", Erange, "optical_absorption_lanc", optical_absorption_lanc)
    save("1d/eigs_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(Z_tol).jld2", "ev", ev, "ef", ef)
end

## compute errors

N_unit = 128
N_v = 4
N_c = 5
N_k = 256
N_iter = 200

Z_tol_vec = [0.8, 0.4, 0.2, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6]
evs = [load("1d/eigs_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(Z_tol).jld2", "ev")[1] for Z_tol in Z_tol_vec]
optical_absorption_lanc = [load("1d/optical_absorption_lanczos_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(N_iter)_$(Z_tol).jld2", "optical_absorption_lanc") for Z_tol in Z_tol_vec]

ev_ref = minimum(real.(load("1d/H_exact_$(N_unit)_$(N_v)_$(N_c)_$(N_k).jld2", "ev")))
Erange, optical_absorption_ref = load("1d/optical_absorption_lanc_ref_$(N_unit)_$(N_v)_$(N_c)_$(N_k)_$(N_iter).jld2", "Erange", "optical_absorption_lanc")

errors_optical_absorption = []
errors_ground_state_energy = []
for i in 1:length(Z_tol_vec)
    push!(errors_optical_absorption, norm(optical_absorption_lanc[i] - optical_absorption_ref, 1) / norm(optical_absorption_ref, 1))
    push!(errors_ground_state_energy, abs(evs[i] - ev_ref))
end

save("1d/errors_spectrum.jld2", "Z_tol_vec", Z_tol_vec,  "errors_optical_absorption", errors_optical_absorption, "errors_ground_state_energy", errors_ground_state_energy)


text:
using Test
using DiffEqPhysics, ForwardDiff, OrdinaryDiffEq
using StaticArrays, LinearAlgebra, Random

test_solve(prob...) = mapreduce(p->solve(p, Tsit5(), dt=1//2).u, ==, prob)

p0, q0   = rand(2)
H(dθ, θ, p) = dθ / 2 - 9.8 * cos(θ)
dp(dθ, θ, p, t) = - 9.8 * sin(θ)
dq(dθ, θ, p, t) = 0.5
acc      = (v, x, p, t) -> ForwardDiff.derivative(x -> -H(v[1], x, p), x[1])
vel      = (v, x, p, t) -> ForwardDiff.derivative(v -> H(v, x[1], p), v[1])
prob_1   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))

@testset "prob1($h)" for h in (H, (dp, dq))
    prob1    = HamiltonianProblem(h, p0, q0, (0., 10.))
    @test test_solve(prob1, prob_1)
end

p0 = @SVector rand(2)
q0 = @SVector rand(2)
H4(dθ, θ, p, t) = dθ[1] / 2 + dθ[2] / 2 - 9.8 * cos(θ[1]) - 9.8 * cos(θ[2]) + t
dq(dθ, θ, p, t) = @SVector [0.5, 0.5]
dp(dθ, θ, p, t) = @SVector [-9.8 * sin(θ[1]), -9.8 * sin(θ[2])]
acc      = (v, x, p, t) -> ForwardDiff.gradient(x -> -H4(v, x, p, t), x)
vel      = (v, x, p, t) -> ForwardDiff.gradient(v -> H4(v, x, p, t), v)
prob_2   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))

@testset "prob2($h)" for h in (H4, (dp, dq))
    prob2    = HamiltonianProblem(h, p0, q0, (0., 10.))
    @test test_solve(prob2, prob_2)
end

H(dθ, θ, p) = (dθ / 2 - 9.8 * cos.(θ))[1]
p0, q0 = [rand(5) for i in 1:2]
dq(dx, dΘ, θ, p, t) = (dx .= [0.5, 0, 0, 0, 0])
dp(dv, dθ, θ, p, t) = (dv .= [(-9.8 * sin.(θ))[1], 0, 0, 0, 0])
acc = (dv, v, x, p, t) -> ForwardDiff.gradient!(dv, x -> -H(v, x, p), x)
vel = (dx, v, x, p, t) -> ForwardDiff.gradient!(dx, v -> H(v, x, p), v)
prob_3 = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))

@testset "prob3($h)" for h in (H, (dp, dq))
    prob3  = HamiltonianProblem(h, p0, q0, (0., 10.))
    @test test_solve(prob3, prob_3)
end

H4(dθ, θ, p, t) = (dθ / 2 - 9.8 * cos.(θ))[1] + t
p0, q0 = [rand(5) for i in 1:2]
dq(dx, dΘ, θ, p, t) = (dx .= [0.5, 0, 0, 0, 0])
dp(dv, dθ, θ, p, t) = (dv .= [(-9.8 * sin.(θ))[1], 0, 0, 0, 0])
acc = (dv, v, x, p, t) -> ForwardDiff.gradient!(dv, x -> -H4(v, x, p, t), x)
vel = (dx, v, x, p, t) -> ForwardDiff.gradient!(dx, v -> H4(v, x, p, t), v)
prob_4 = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))

@testset "prob4($h)" for h in (H4, (dp, dq))
    prob4  = HamiltonianProblem(h, p0, q0, (0., 10.))
    @test test_solve(prob4, prob_4)
end


text:
using Graphics: @mustimplement

export Amplifier, AmplifierMode, AMP_VOLTAGE_MODE, AMP_CURRENT_MODE, AmplifierVoltageMode,
       AMP_HIGH_VOLTAGE_MODE, AMP_LOW_VOLTAGE_MODE, getAmplifiers, getAmplifier, state,
       turnOn, turnOff, mode, voltageMode, matchingNetwork, temperature, toCurrentMode,
       toVoltageMode, toLowVoltageMode, toHighVoltageMode

@enum AmplifierMode begin
  AMP_VOLTAGE_MODE
  AMP_CURRENT_MODE
end

function convert(::Type{AmplifierMode}, x::String)
  if lowercase(x) == "voltage"
    return AMP_VOLTAGE_MODE
  elseif  lowercase(x) == "current"
    return AMP_CURRENT_MODE
  else
    throw(ScannerConfigurationError("The given amplifier mode `$x` for is not valid. Please use `voltage` or `current`."))
  end
end

@enum AmplifierVoltageMode begin
  AMP_HIGH_VOLTAGE_MODE
  AMP_LOW_VOLTAGE_MODE
end

function convert(::Type{AmplifierVoltageMode}, x::String)
  if lowercase(x) == "low voltage"
    return AMP_LOW_VOLTAGE_MODE
  elseif  lowercase(x) == "high voltage"
    return AMP_HIGH_VOLTAGE_MODE
  else
    throw(ScannerConfigurationError("The given amplifier mode `$x` for is not valid. Please use `low voltage` or `high voltage`."))
  end
end

abstract type Amplifier <: Device end

Base.close(amp::Amplifier) = nothing

@mustimplement state(amp::Amplifier)
@mustimplement turnOn(amp::Amplifier)
@mustimplement turnOff(amp::Amplifier)
@mustimplement mode(amp::Amplifier)::AmplifierMode
@mustimplement mode(amp::Amplifier, mode::AmplifierMode)
@mustimplement voltageMode(amp::Amplifier)::AmplifierVoltageMode
@mustimplement voltageMode(amp::Amplifier, mode::AmplifierVoltageMode)
@mustimplement matchingNetwork(amp::Amplifier)::Integer
@mustimplement matchingNetwork(amp::Amplifier, network::Integer)
@mustimplement temperature(amp::Amplifier)::typeof(1.0u"°C")

getAmplifiers(scanner::MPIScanner) = getDevices(scanner, Amplifier)
function getAmplifier(scanner::MPIScanner)
  amplifiers = getGaussMeters(scanner)
  if length(amplifiers) > 1
    error("The scanner has more than one amplifier device. Therefore, a single amplifier cannot be retrieved unambiguously.")
  else
    return amplifiers[1]
  end
end

"""
Sets the amplifier to current mode.
"""
toCurrentMode(amp::Amplifier) = mode(amp, AMP_CURRENT_MODE)

"""
Sets the amplifier to voltage mode.
"""
toVoltageMode(amp::Amplifier) = mode(amp, AMP_VOLTAGE_MODE)

"""
Sets the amplifier to low voltage mode.
"""
toLowVoltageMode(amp::Amplifier) = voltageMode(amp, AMP_LOW_VOLTAGE_MODE)

"""
Sets the amplifier to high voltage mode.
"""
toHighVoltageMode(amp::Amplifier) = voltageMode(amp, AMP_HIGH_VOLTAGE_MODE)

include("SimulatedAmplifier.jl")
include("HubertAmplifier.jl")


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

# Register a Prc to a clock. Check its id and change it apropriately.
function _register!(clk::Clock, p::Prc)
    id = p.id
    while haskey(clk.processes, id)
        if isa(id, Float64)
            id = nextfloat(id)
        elseif isa(id, Int)
            id += 1
        elseif isa(id, String)
            s = split(id, "#")
            id = length(s) > 1 ? chop(id)*string(s[end][end]+1) : id*"#1"
        else
            throw(ArgumentError("process id $id is duplicate, cannot convert!"))
        end
    end
    clk.processes[id] = p
    p.id = id
end

# Put a [`Prc`](@ref) in a loop which can be broken by a `ClockException`.
# - `p::Prc`:
# - `cycles=Inf`: determine, how often the loop should be run.
function _loop(p::Prc, cycles::T) where {T<:Number}
    threadid() > 1 && (p.clk = _clock(pclock(p.clk)))
    _register!(p.clk, p)
    while cycles > 0
        try
            p.f(p.clk, p.arg...; p.kw...)
        catch exc
            if isa(exc, PrcException)
                exc.event == Stop() && break
            end
            rethrow(exc)
        end
        cycles -= 1
    end
    p.clk.processes = delete!(p.clk.processes, p.id)
end

# startup a `Prc` as a task in a loop.
function _startup!(c::C, p::Prc, cycles::T, cid::Int, spawn::Bool) where {C<:AbstractClock,T<:Number}

    function startit()
        p.task = @task _loop(p, cycles)
        yield(p.task)
        return p.task
    end

    t = Task(nothing)
    cid = _cid(c, cid, spawn)
    if cid == c.id
        t = startit()
    else
        @threads for i in 1:nthreads()
            i == cid && (t = startit())
        end
    end
    return t
end

"""
    process!([clk], prc, cycles; <keyword arguments>)

Register a [`Prc`](@ref) to a clock, start an asynchronous task
executing the process function in a loop and return the `id` it
was registered with. It can then be found under `clk.processes[id]`.

# Arguments
- `c<:AbstractClock`: if not provided, the process runs under `𝐶`,
- `prc::Prc`: it contains a function and its arguments,
- `cycles<:Number=Inf`: number of loop cycles the process should run,

# Keyword arguments
- `cid::Int=clk.id`: if cid ≠ clk.id, assign the event to the parallel clock
    with id == cid. This overrides `spawn`,
- `spawn::Bool=false`: if true, the process may be scheduled on another thread
    in parallel and registered to the thread specific clock.
"""
function process!(c::C, p::Prc, cycles::T=Inf;
                  cid::Int=c.id, spawn::Bool=false) where {C<:AbstractClock,T<:Number}
    p.clk = c
    t = _startup!(c, p, cycles, cid, spawn)
    t.state === :failed && return t
    p.id
end
process!(p::Prc, cycles::T=Inf; kwargs...) where {T<:Number} = process!(𝐶, p, cycles; kwargs...)

# wakeup a process waiting for a `Condition`
# two yields for giving back control first and then enabling a
# take! on a channel (heuristic solution of the clock overrun problem)
_wakeup(c::Condition) = (notify(c); yield(); yield())

"""
```
delay!(clk, Δt)
delay!(clk, T, t)
```
Delay (suspend) a process for a time interval `Δt` on the clock `clk`.

# Arguments
- `clk::Clock`,
- `Δt`: time interval, `Number` or `Distribution`,
- `T::Timing`: only `until` is accepted,
- `t`: time delay, `Number` or `Distribution`.
"""
function delay!(clk::Clock, Δt::N) where N<:Number
    c = Condition()
    event!(clk, ()->_wakeup(c), after, Δt)
    wait(c)
end
function delay!(clk::Clock, T::Timing, t::N) where N<:Number
    @assert T == until "bad Timing $T for delay!"
    if t > clk.time
        c = Condition()
        event!(clk, ()->_wakeup(c), t)
        wait(c)
    else
        now!(clk, fun(println, stderr, "warning: delay until $t ≤ τ=$(tau(clk))"))
    end
end
delay!(clk::Clock, Δt::X) where X<:Distribution = delay!(clk, rand(Δt))
delay!(clk::Clock, T::Timing, t::X) where X<:Distribution = delay!(clk, T, rand(t))

"""
    wait!(clk, cond)

Delay (suspend) a process on a clock clk until a condition has become true.

# Arguments
- `clk::Clock`,
- `cond<:Action`: a condition, true if all expressions or functions therein return true.
"""
function wait!(clk::Clock, cond::A) where {A<:Action}
    if all(_evaluate(cond))   # all conditions met
        return         # return immediately
    else
        c = Condition()
        event!(clk, ()->_wakeup(c), cond)
        wait(c)
    end
end

"""
    interrupt!(p::Prc, ev, value)

Interrupt a [`Prc`](@ref) by throwing a [`PrcException`](@ref) to it.
"""
function interrupt!(p::Prc, ev, value)
    schedule(p.task, PrcException(ev, value), error=true)
    yield()
end

"Stop a Prc"
stop!(p::Prc, value=nothing) = interrupt!(p, Stop(), value)

"""
    now!(clk::Clock, ex::A) where {A<:Action}

Transfer an IO-operation `ex` to `clk` or if it is a parallel clock
to the master clock (on thread 1). The clock executes it before
proceeding to the next time step.
"""
function now!(clk::Clock, ex::A) where {A<:Action}
    cid = clk.ac isa Ref ? 1 : threadid()
    event!(clk, ex, clk.time, cid=cid)
end

"""
```
print(clk::Clock, [io::IO], x, xs...)
```
Create a [`now!`](@ref) event to a busy clock `clk` to `print(x, xs...)`
to `io` (or to `stdout`).

If the clock is not busy, `x` and `xs...` are printed as usual.
`print(clk)` still prints `repr(clk)`.
"""
function Base.print(clk::Clock, io::IO, x, xs...)
    if clk.state == Busy()
        now!(clk, fun(print, io, x))
        for x_ in xs
            now!(clk, fun(print, io, x_))
        end
    else
        print(io, x, xs...)
    end
end
Base.print(clk::Clock, x, xs...) = print(clk, stdout, x, xs...)

"""
```
println(clk::Clock, [io::IO], xs...)
```
Create a [`now!`](@ref) event to a busy clock `clk` to `println(x, xs...)`
to `io` (or to `stdout`).

If the clock is not busy, `x` and `xs...` are printed as usual.
`println(clk)` still prints `repr(clk)`.
"""
Base.println(clk::Clock, io::IO, x, xs...) = print(clk, io, x, xs..., '\n')
Base.println(clk::Clock, x, xs...) = print(clk, x, xs..., '\n')


text:
## karhunen_loeve.jl : approximate Gaussian random field generator using a Karhunen-Lo\`eve decomposition

"""
    KarhunenLoeve(n)

Returns a [`GaussianRandomFieldGenerator`](@ref) using a Karhunen-Loève (KL) expansion with `n` terms.

# Optional Arguments for `GaussianRandomField`
- `quad::QuadratureRule`: quadrature rule used for the integral equation (see [`QuadratureRule`](@ref)), default is `EOLE()`.
- `nq::Integer`: number of quadrature points in each dimension, where we require `nq^d > n`. Default is `nq = ceil(n^(1/d))`.
- `eigensolver::EigenSolver`: which method to use for the eigenvalue decomposition (see [`AbstractEigenSolver`](@ref)). The default is `EigenSolver()`.

# Examples
```jldoctest label1
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, KarhunenLoeve(300), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 300 terms

julia> plot_eigenvalues(grf) # plot the eigenvalue decay
[...]

julia> plot_eigenfunction(grf, 4) # plots the fourth eigenfunction
[...]

```
If more terms `n` are used in the expansion, the approximation becomes better.
```jldoctest label1; setup=:(import Random; Random.seed!(12345))
julia> for n in [1, 2, 5, 10, 20, 50, 100, 200, 500, 1000]
           grf = GaussianRandomField(cov, KarhunenLoeve(n), pts, pts)
           println(rel_error(grf))
       end
0.6983828486813854
0.454941868632304
0.23231277904920067
0.10079295313241687
0.026470201665282467
0.009784266729696567
0.003565488318168386
0.0010719081249264129
0.00019809766995382283
4.085273649512278e-5

```
!!! note
    Techniqually, the KL expansion is computed using the Nystrom method. For nonstructured grids, we use a bounding box approach. Try using the `Spectral` method if this is not what you want.

!!! warning
    To avoid an *end effect* in the eigenvalue decay, choose `nq^d ≥ 5n`.

```jldoctest label1
julia> grf = GaussianRandomField(cov, KarhunenLoeve(300), pts, pts, quad=GaussLegendre())
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 300 terms

julia> grf = GaussianRandomField(cov, KarhunenLoeve(300), pts, pts, nq=40)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 300 terms

```
See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`CirculantEmbedding`](@ref)
"""
struct KarhunenLoeve{n} <: GaussianRandomFieldGenerator
    function KarhunenLoeve{n}() where n
        n > 0 || throw(DomainError(n, "in KarhunenLoeve{n}, number of terms n must be positive"))
        new{n}()
    end
end

KarhunenLoeve(n::Integer) = KarhunenLoeve{n}()

function _GaussianRandomField(mean, cov::CovarianceFunction{d}, method::KarhunenLoeve{n},
                              pts...;
                              nq::Integer=ceil(typeof(n), n^(1/d)),
                              quad::QuadratureRule=EOLE(),
                              eigensolver::AbstractEigenSolver=EigsSolver()) where {d,n}
    # check if number of terms and number of quadrature points are compatible
    nq = nq > 0 ? ntuple(i -> nq, d) : length.(pts)
    prod_nq = prod(nq)
    prod_nq < n && throw(ArgumentError("too many terms requested, increase nq or lower n"))
    # adjustment for ARPACK error when looking for ALL eigenvalues
    eigensolver isa EigsSolver && prod_nq == n && (nq = nq .+ 1)

    # determine bounding box when irregular mesh
    is_irr = pts[1] isa AbstractMatrix
    a = is_irr ? Tuple(minimum(pts[1],dims=2)) : minimum.(pts)
    b = is_irr ? Tuple(maximum(pts[1],dims=2)) : maximum.(pts)

    # compute quadrature nodes and weights
    struc = get_nodes_and_weights.(nq, a, b, Ref(quad))
    nodes = first.(struc)
    weights = last.(struc)

    # eigenvalue problem
    C = apply(cov, nodes...)
    W = d == 1 ? Diagonal(weights...) : Diagonal(kron(weights...))
    W .= sqrt.(W)
    B = Symmetric(W * C * W) # should be symmetric and positive semi-definite
    isposdef(B) || @warn "equivalent eigenvalue problem is not SPD, results may be wrong or inaccurate"

    # solve
    eigenval, eigenfunc = compute(B, n, eigensolver)

    # compute eigenfunctions in nodes
    K = apply(cov.cov, pts, nodes)
    Λ = Diagonal(1 ./ eigenval)
    eigenfunc = K * W * eigenfunc * Λ

    m = findfirst(x -> x < 0, eigenval)
    if m != nothing
        m -= 1
        @warn "$(length(eigenval) - m) negative eigenvalues ≥ $(eigenval[end]) detected,
            Gaussian random field will be approximated (ignoring all negative eigenvalues)"
    else
        m = n
    end

    resize!(eigenval, m)
    eigenfunc = eigenfunc[:, 1:m]

    for id=1:m
    if(eigenfunc[1,id]<0)
        eigenfunc[:,id]=eigenfunc[:,id].*-1
    end
    end

    # store eigenvalues and eigenfunctions
    eigenval .= sqrt.(eigenval) # note: store sqrt of eigenval for more efficient sampling
    data = SpectralData(eigenval, eigenfunc)

    GaussianRandomField{KarhunenLoeve{length(eigenval)},typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean,cov,pts,data)
end

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{KarhunenLoeve{n}}) where n = n

# relative error in the KL approximation
"""
    rel_error(grf)

Returns the relative error in the Karhunen-Loève approximation of the random field, computed as

``1 - \\displaystyle\\frac{\\sum \\theta_j^2}{\\sigma^2 \\int_D \\mathrm{d}x}``.

Only useful for fields defined on a rectangular domain.

# Examples
```jldoctest; setup=:(import Random; Random.seed!(12345))
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, KarhunenLoeve(300), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a KL expansion with 300 terms

julia> rel_error(grf)
0.00046070730242930846

```
"""
function rel_error(grf::GaussianRandomField{<:KarhunenLoeve})
    Leb = prod(grf.pts) do point
        a, b = extrema(point)
        b - a
    end
    1 - sum(abs2, grf.data.eigenval) / (std(grf.cov)^2 * Leb)
end

# sample function for both Spectral() and KarhunenLoeve(n) type
function _sample(grf::GaussianRandomField{<:Union{Spectral,KarhunenLoeve}}, xi)
    grf.mean + std(grf.cov) * reshape((grf.data.eigenfunc * Diagonal(grf.data.eigenval)) * xi, size(grf.mean))
end

Base.show(io::IO, ::KarhunenLoeve{n}) where n = print(io, "KL expansion with $(n) terms")


text:
using IMinuit
using Test
# using DataFrames, CSV

# @testset "IMinuit.jl" begin
    # Write your tests here.
# end

@testset "minimizing simple functions: " begin
    f(x) = x[1]^2 + (x[2]-1)^2 + (x[3]-2)^4
    f1(x, y, z) = x^2 + (y-1)^2 + (z-2)^4

    m = Minuit(f, [1,1,0])
    m1 = Minuit(f1, x = 1, y = 1, z = 0)

    migrad(m)
    migrad(m1)

    @test typeof(m) == ArrayFit
    @test typeof(m1) == Fit

    # the second parameter should be about 1
    @test get(m.values, 1) ≈ 1
    @test get(m1.values, 1) ≈ 1

    @test length(get_contours_all(m, f, npts = 3)) == 19
    @test length(get_contours_all(m1, f1, npts = 3)) == 19
end
#
# @testset "a simple but realistic fit: " begin
#     datadf = DataFrame!(CSV.File("../docs/testdata.csv"))
#     data = Data(datadf)
#
#     M = 3.686; mπ = 0.14; mJ = 3.097;
#
#     λ(x, y, z) = x^2 + y^2 + z^2 - 2x*y - 2y*z - 2z*x
#
#     # a simple function that will be used to fit the data: QCD multipole expansion model for ψ'→J/ψπ⁺π⁻
#     # The important ππ FSI effect is not taken into account
#     # bg is just for introducing a third parameter
#     function dist(w, N, c, bg)
#         if (w ≤ 2mπ || w ≥ M-mJ)
#             res = 0.0
#         else
#             q1 = sqrt(λ(w^2, mπ^2, mπ^2))/(2w)
#             q2 = sqrt(λ(M^2, w^2, mJ^2))/(2M)
#             res = N * q1 * q2 * (w^2 - c*mπ^2)^2 + bg
#         end
#         return res * 1e6
#     end
#
#     # dist(w, par) = dist(w, par...);
#
#     χsq(par) = chisq(dist, data, par)
#     χsq1(N, c, bg) = chisq(dist, data, (N, c, bg));
#
#     gradf(par) = gradient(χsq, par)
#     fit = Minuit(χsq, [1, 2, 0], error = 0.1*ones(3), grad = gradf)
#     fit.strategy = 0;
#
#     fit1 = Minuit(χsq1, N = 1, c = 2, bg = 0, error_N = 0.1, error_c = 0.1, error_bg = 0.1)
#     fit1.strategy = 0;
#
#     migrad(fit)
#
#     migrad(fit1)
#
#     @test isapprox(get(fit.values, 0), 2.61, rtol = 1e-2)
#     @test isapprox(get(fit1.values, 0), 2.61, rtol = 1e-2)
#
# end


text:
"""
    ThresholdError <: Exception

Is thrown when a Threshold limit is exceed.

# Fields
* limit::Float64 - the threshold limit.
* value::Float64 - the missing data ratio identified
"""
struct ThresholdError <: Exception
    limit::Float64
    value::Float64
end

function Base.showerror(io::IO, err::ThresholdError)
    println(
        io,
        "ThresholdError: Missing data limit exceeded $(err.limit) ($(err.value))",
    )
end

"""
    Threshold(; limit=0.1)

Assert that the ratio of missing values in the provided dataset does not exceed to
specified limit.

# Keyword Arguments
* `limit::Real`: Allowed proportion of missing values (should be between 0.0 and 1.0).
"""
struct Threshold <: Validator
    limit::Float64
    weights::Union{AbstractWeights, Nothing}
end

Threshold(; limit=0.1, weights=nothing) = Threshold(limit, weights)

function _validate(data::AbstractArray{Union{T, Missing}}, t::Threshold) where T
    mratio = count(ismissing, data) / length(data)
    mratio > t.limit && throw(ThresholdError(t.limit, mratio))
    return data
end

"""
    WeightedThreshold(; limit, weights)

Assert that the weighted proportion missing values in the provided dataset does not exceed
to specified limit. The weighed proportion is calculated as
`sum(weights[ismissing.(data)]) / sum(weights)`

# Keyword Arguments
* `limit::Real`: Allowed proportion of missing values (should be between 0.0 and 1.0).
* `weights::AbstractWeights`: A set of statistical weights to use when evaluating the importance
  of each observation.
"""
struct WeightedThreshold{W <: AbstractArray{<:Real}} <: Validator
    limit::Float64
    weights::W
end

WeightedThreshold(; limit, weights) = WeightedThreshold(limit, weights)

function _validate(data::AbstractArray{Union{T, Missing}}, wt::WeightedThreshold) where T
    if size(data) != size(wt.weights)
        throw(DimensionMismatch(string(
            "Input has dimensions $(size(data)), but thresholds weights ",
            "has dimensions $(size(wt.weights))"
        )))
    end

    val = sum(wt.weights[ismissing.(data)]) / sum(wt.weights)
    val > wt.limit && throw(ThresholdError(wt.limit, val))

    return data
end


text:
using LLVM
using LLVM.Compat
using CUDAnative

##
# TODO:
# - Do we need to support address-spaces
# - Alignment
##

abstract type AbstractSemaphore end

const SUCCESS = Int32(0)
const ERROR_TIMEOUT = Int32(11)
const ERROR_INVALID = Int32(22)

# must match verbs_exp enum
const GDS_WAIT_COND = Int32
const GDS_WAIT_COND_GEQ = Int32(0)
const GDS_WAIT_COND_EQ  = Int32(1)
const GDS_WAIT_COND_AND = Int32(2)
const GDS_WAIT_COND_NOR = Int32(3)

# Not currently needed
# mutable struct Semaphore{T} <: AbstractSemaphore
#     sem::T
#     value::T
# end
#
# const Sem32 = Semaphore{UInt32}
# const Sem64 = Semaphore{UInt64}
#
# @inline function Base.getindex(sem::ISemaphore{T}) where T
#     # volatile load from field
# end
#
# @inline function Base.setindex(sem::ISemaphore{T}, val::T) where T
#     # volatile store to field
# end

struct ISemaphore{T} <: AbstractSemaphore
    sem::Ptr{T}
    value::T
    ISemaphore() = new(C_NULL, 0)
end
const ISem32 = ISemaphore{UInt32}
const ISem64 = ISemaphore{UInt64}

@generated function Base.getindex(sem::ISemaphore{T}) where T
    eltyp = convert(LLVMType, T)
    T_ptr = convert(LLVMType, Ptr{T})
    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr]
    llvm_f, _ create_function(eltyp, param_types)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        ld = load!(builder, ptr)
        LLVM.API.LLVMSetVolatile(ld, LLVM.True)
        ret!(builder, ld)
    end

    call_function(llvm_f, T, Tuple{Ptr{T}}, :(sem.sem,))
end

@generated function Base.setindex!(sem::ISemaphore{T}, val::T) where T
    eltyp = convert(LLVMType, T)
    T_ptr = convert(LLVMType, Ptr{T})
    T_actual_ptr = LLVM.PointerType(eltyp)

    # create a function
    param_types = [T_ptr, eltyp]
    llvm_f, _ create_function(LLVM.VoidType(JuliaContext()), param_types)

    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = inttoptr!(builder, parameters(llvm_f)[1], T_actual_ptr)
        val = parameters(llvm_f)[2]
        st = store!(builder, val, ptr)
        LLVM.API.LLVMSetVolatile(st, LLVM.True)
        ret!(builder)
    end

    call_function(llvm_f, Cvoid, Tuple{Ptr{T}, T}, :(sem.sem, val))
end

@inline function release(sem::Semaphore)
    @cuassert(C_NULL != sem[])
    sem[] = sem.value
end

@inline function wait(sem::Semaphore, cond)
    ret = SUCCESS
    if cond == GDS_WAIT_COND_EQ
        ret = wait(sem, ==)
    elseif cond == GDS_WAIT_COND_GEQ
        ret = wait(sem, >=)
    elseif cond == GDS_WAIT_COND_AND
        # TODO: this is just (a&b) in C
        ret = wait(sem, (a, b) -> (a & b) == 1)
    elseif cond == GDS_WAIT_COND_NOR
        ret = wait(sem, (a, b) -> ~(a | b) != 0)
    else
        ret = ERROR_INVALID
    end
    return ret
end

@inline function wait(sem::Semaphore, cond::F) where F<:Function
    timeout = clock(UInt64) + UInt64(1) << 32 # this is marked volatile
    while (true)
        if cond(sem[], sem.value)
            return SUCCESS
        end
        threadfence_block()
        if clock(UInt64) >= timeout
            return ERROR_TIMEOUT
        end
    end
end


text:
module TestDefaults

using Mimi
using Test

import Mimi: model_params

@defcomp A begin
    p1 = Parameter(default = 1)
    p2 = Parameter{Symbol}()
end

m = Model()
set_dimension!(m, :time, 1:10)
add_comp!(m, A)

add_shared_param!(m, :p2, :hello)
connect_param!(m, :A, :p2, :p2)

# So far only :p2 is in the model definition's dictionary
@test :p2 in keys(model_params(m))
@test length(model_params(m)) == 2

run(m)

# :p1's value is it's default
@test m[:A, :p1] == 1

# This errors because p1 is unshared
@test_throws ErrorException update_param!(m, :p1, 10)
update_param!(m, :A, :p1, 10)

# :p1 still not in the dictionary because unshared
@test !(:p1 in keys(model_params(m)))

# now add it as a shared parameter
add_shared_param!(m, :model_p1, 20)
connect_param!(m, :A, :p1, :model_p1)

# Now there is a :model_p1 in the model definition's dictionary but not :p1
@test !(:p1 in keys(model_params(m)))
@test :model_p1 in keys(model_params(m))

run(m)
@test m[:A, :p1] == 20

# Now we can use update_param! but only for the model parameter name and exclusively as a shared parameter
@test_throws ErrorException update_param!(m, :p1, 11)
update_param!(m, :model_p1, 30)
run(m)
@test m[:A, :p1] == 30

# convert explicitly back to being unshared
@test_throws ErrorException update_param!(m, :A, :p1, 40)
disconnect_param!(m, :A, :p1)
update_param!(m, :A, :p1, 40)
run(m)
@test m[:A, :p1] == 40
@test Mimi.model_param(m, :model_p1).value == 30

# try setting defaults for Strings or Symbols
@defcomp A begin
    p1 = Parameter{String}(default = "hello")
    p2 = Parameter{Symbol}(default = :goodbye)
end

m = Model()
set_dimension!(m, :time, 1:10)
add_comp!(m, A)
run(m)

@test m[:A, :p1] == "hello"
@test m[:A, :p2] == :goodbye

end


text:
function sieve(limit::Integer)
    limit <= 0 && error("Invalid limit")
    nums = fill(true, limit)
    nums[1] = false
    primes = []
    while (i = findfirst(nums)) > 0
        push!(primes, i)
        for j = find(nums)
            if j % i == 0
                nums[j] = false
            end
        end
    end
    primes
end


text:
include("_lines.jl")
include("_plots.jl")

include("validate.jl")


text:
# Univariate

const FillVectorOfUnivariate{
    S <: ValueSupport,
    T <: UnivariateDistribution{S},
    Tdists <: Fill{T, 1},
} = VectorOfUnivariate{S, T, Tdists}

function filldist(dist::UnivariateDistribution, N::Int)
    return product_distribution(Fill(dist, N))
end
filldist(d::Normal, N::Int) = TuringMvNormal(fill(d.μ, N), d.σ)

function Distributions.logpdf(
    dist::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return _logpdf(dist, x)
end
function Distributions.logpdf(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return _logpdf(dist, x)
end
ZygoteRules.@adjoint function Distributions.logpdf(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return ZygoteRules.pullback(_logpdf, dist, x)
end

function _logpdf(
    dist::FillVectorOfUnivariate,
    x::AbstractVector{<:Real},
)
    return _flat_logpdf(dist.v.value, x)
end
function _logpdf(
    dist::FillVectorOfUnivariate,
    x::AbstractMatrix{<:Real},
)
    return _flat_logpdf_mat(dist.v.value, x)
end

function _flat_logpdf(dist, x)
    if toflatten(dist)
        f, args = flatten(dist)
        return sum(f.(args..., x))
    else
        return sum(map(x) do x
            logpdf(dist, x)
        end)
    end
end
function _flat_logpdf_mat(dist, x)
    if toflatten(dist)
        f, args = flatten(dist)
        return vec(sum(f.(args..., x), dims = 1))
    else
        temp = map(x -> logpdf(dist, x), x)
        return vec(sum(temp, dims = 1))
    end
end

function Distributions.rand(rng::Random.AbstractRNG, d::FillVectorOfUnivariate)
    return rand(rng, d.v.value, length(d))
end
function Distributions.rand(rng::Random.AbstractRNG, d::FillVectorOfUnivariate, n::Int)
    return rand(rng, d.v.value, length(d), n)
end

const FillMatrixOfUnivariate{
    S <: ValueSupport,
    T <: UnivariateDistribution{S},
    Tdists <: Fill{T, 2},
} = MatrixOfUnivariate{S, T, Tdists}

function filldist(dist::UnivariateDistribution, N1::Integer, N2::Integer)
    return MatrixOfUnivariate(Fill(dist, N1, N2))
end
function Distributions.logpdf(dist::FillMatrixOfUnivariate, x::AbstractMatrix{<:Real})
    return _flat_logpdf(dist.dists.value, x)
end
function Distributions.rand(rng::Random.AbstractRNG, dist::FillMatrixOfUnivariate)
    return rand(rng, dist.dists.value, length.(dist.dists.axes)...,)
end

# Multivariate

const FillVectorOfMultivariate{
    S <: ValueSupport,
    T <: MultivariateDistribution{S},
    Tdists <: Fill{T, 1},
} = VectorOfMultivariate{S, T, Tdists}

function filldist(dist::MultivariateDistribution, N::Int)
    return VectorOfMultivariate(Fill(dist, N))
end
function Distributions.logpdf(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return _logpdf(dist, x)
end
function _logpdf(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return sum(logpdf(dist.dists.value, x))
end
ZygoteRules.@adjoint function Distributions.logpdf(
    dist::FillVectorOfMultivariate,
    x::AbstractMatrix{<:Real},
)
    return ZygoteRules.pullback(_logpdf, dist, x)
end
function Distributions.rand(rng::Random.AbstractRNG, dist::FillVectorOfMultivariate)
    return rand(rng, dist.dists.value, length.(dist.dists.axes)...,)
end


text:
# Test 2D modeling
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

# Set parallel if specified
if nw > 1 && nworkers() < nw
   addprocs(nw-nworkers() + 1; exeflags=["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
end

@everywhere using JOLI
@everywhere using JUDI, LinearAlgebra, Test, Distributed

### Model
model, model0, dm = setup_model(tti, viscoacoustic, nlayer; n=(101, 101), d=(10., 10.))
q, srcGeometry, recGeometry, f0 = setup_geom(model; nsrc=2, tn=500f0)
dt = srcGeometry.dt[1]

# Modeling operators
println("Generic modeling and misc test with ", nlayer, " layers and tti: ", tti)

ftol = 1f-5
######################## WITH DENSITY ############################################

cases = [(true, true), (true, false), (false, true), (false, false)]

@testset "Generic tests with limit_m = $(limit_m)  and save to disk = $(disk)" for (limit_m, disk)=cases
	@timeit TIMEROUTPUT "Modeling  limit_m=$(limit_m), save_to_disk=$(disk)" begin
		# Options structures
		opt = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0, f0=f0,
						file_path=pwd(),	# path to files
						file_name="shot_record")	# saves files as file_name_xsrc_ysrc.segy

		opt0 = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0, f0=f0,
						file_path=pwd(),	# path to files
						file_name="smooth_shot_record")	# saves files as file_name_xsrc_ysrc.segy

		optJ = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0, f0=f0,
						file_path=pwd(),	# path to files
						file_name="linearized_shot_record")	# saves files as file_name_xsrc_ysrc.segy

		# Setup operators
		Pr = judiProjection(recGeometry)
		F = judiModeling(model; options=opt)
		F0 = judiModeling(model0; options=opt0)
		Ps = judiProjection(srcGeometry)

		# Combined operator Pr*F*adjoint(Ps)
		Ffull = judiModeling(model, srcGeometry, recGeometry)
		J = judiJacobian(Pr*F0*adjoint(Ps),q; options=optJ) # equivalent to J = judiJacobian(Ffull,q)

		# Nonlinear modeling
		d1 = Pr*F*adjoint(Ps)*q	# equivalent to d = Ffull*q
		dfull = Ffull*q
		@test isapprox(get_data(d1), dfull, rtol=ftol)

		qad = (Ps*adjoint(F))*adjoint(Pr)*d1
		qfull = adjoint(Ffull)*d1
		@test isapprox(qad, qfull, rtol=ftol)

		# fwi objective function
		f, g = fwi_objective(model0, q, d1; options=opt)
		f, g = fwi_objective(model0, getindex(q,1), getindex(d1,1); options=opt)

		# Indexing (per source)
		for inds=[2, [1, 2]]
			dsub = getindex(dfull, inds)
			qsub = getindex(q, inds)
			Fsub = getindex(F, inds)
			Jsub = getindex(J, inds)
			Ffullsub = getindex(Ffull, inds)
			Pssub = getindex(Ps, inds)
			Prsub = getindex(Pr, inds)
			ds1 = Ffullsub*qsub
			ds2 = Prsub * Fsub * adjoint(Pssub) *qsub
			@test isapprox(ds1, get_data(ds2), rtol=ftol)
			@test isapprox(ds1, dsub, rtol=ftol)
			@test isapprox(get_data(ds2), dsub, rtol=ftol)
		end

		# vcat, norms, dot
		dcat = [d1, d1]
		@test isapprox(norm(d1)^2, .5f0*norm(dcat)^2)
		@test isapprox(dot(d1, d1), norm(d1)^2)
	end
end

############################# Full wavefield ############################################

@testset "Basic judiWavefield modeling tests" begin
	@timeit TIMEROUTPUT "Wavefield modeling" begin
		opt = Options(dt_comp=dt, f0=f0)
		F = judiModeling(model; options=opt)
		Fa = adjoint(F)
		Ps = judiProjection(srcGeometry)
		Pr = judiProjection(recGeometry)

		# Return wavefields
		u = F * adjoint(Ps) * q

		# Adjoint from data
		dobs = Pr*F*u
		v = Fa*adjoint(Pr)*dobs

		a = dot(u, v)
		b = dot(dobs, dobs)
		@printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", b, a, (a-b)/(a+b))
		@test isapprox(a/(a+b), b/(a+b), atol=1f-5, rtol=0)

		# Forward from data
		qa = Ps*Fa*v

		a = dot(u, v)
		b = dot(q, qa)
		@printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", b, a, (a-b)/(a+b))
		@test isapprox(a/(a+b), b/(a+b), atol=1f-5, rtol=0)

		# Wavefields as source + return wavefields
		u2 = F*u
		v2 = Fa*v

		a = dot(u2, v)
		b = dot(v2, u)
		@printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", a, b, (a-b)/(a+b))
		@test isapprox(a/(a+b), b/(a+b), atol=1f-4, rtol=0)
	end
end


text:
function load_buffer(handle::Vector{UInt8}, size::Int)

    # get key as UInt32, then get string value instead of len 1 array
    shmkey = reinterpret(UInt32, handle)[1]

    # use <sys/ipc.h> from C standard library to get shared memory id
    # validate that shmget returns a valid id
    shmid = ccall((:shmget, "libc"), Cint, (Cuint, Int32, Int32), shmkey, size, 0)
    shmid == -1 ? error("Invalid shared memory key: $shmkey") : nothing

    # with shmid, get shared memory start address
    ptr = ccall((:shmat, "libc"), Ptr{Nothing}, (Cint, Ptr{Nothing}, Cint), shmid, C_NULL, 0)

    # makes a zero-copy reference to memory, true gives ownership to julia
    # validate that memory no longer needs to be released using MapD methods
    return unsafe_wrap(Array, convert(Ptr{UInt8}, ptr), size)

end


text:
using MinimumViableTables
using Test

using MinimumViableTables: _issubset, _issetequal, _headsubset, _makevectors, _add, _all,
                           clean, searchsortedlastless, searchsortedfirstgreater,
                           _valuetype, _cat_types, _intersect, _setdiff, _find

# TODO Make an @inferred-like macro to test inference with constant propagation

include("interval.jl")

include("util.jl")
include("project.jl")
include("rename.jl")
include("index.jl")
include("table.jl")
include("rowoperations.jl")
include("columnoperations.jl")
include("makeindex.jl")
include("predicates.jl")

include("custom_arrays.jl")
include("producttable.jl")


text:
################
# Julien  Pascal
################

#################################
# PATHS
# change "path_main" if necessary
home = ENV["HOME"]
path_main = "/home/julien/Final-Project-Julien-Pascal"
cd(path_main)

# Tests:
include("test/runtests.jl")

# Check for the existence of crucial files:
# A. Check that the empirical moments have been calculated
# Otherwise has to launch the python code "Calculate_Moments.py"
if isfile("tables/table_moments.csv") == false
	println("tables/table_moments.csv is missing")
	error("Execute Calculate_Moments.py to create it")
end

# B. Check that the parameters of the model have been estimated
# Otherwise estimate them by the method of simulated moments
if isfile("tables/estimated_parameters.csv") == false
	println("tables/estimated_parameters.csv is missing")
	println("launching the simulated method of moments procedure")
	println("this may take a while")

	include("src/estimate_parameters.jl")

	using estimate_parameters
	#= first argument: minimization algorithm to use
		- "NLopt": use Controlled Random Search (CRS) with local mutation
		- "DE": Differential Evolution
		- "SA": Simulated Annealing
	# second argument: 1-step or 2-step procedure (using the optimal weighting matrix)
		- "1Step": use the Identity matrix as the weighting matrix
		- "2Step": use the optimal weighting matrix
	# by default: "NLopt", "1Step"
	=#
	MSM("LOCAL", "1Step")
end

# C. Are the value functions for the wages there?
# If not, calculate them by value function iteration
if isfile("surpluses/W_high.jld") == false  || isfile("surpluses/W_high_star.jld") == false  || isfile("surpluses/W_low_star.jld") == false  || isfile("surpluses/W_low.jld") == false
	println("Wages surpluses are missing")
	println("Calculating them by value function iteration")
	println("this may take a while")

	include("src/calculate_wages.jl")
	using calculate_wages
	Wages() #2380 seconds to run
end

################################################
# D. If the file are there, simulate one economy
################################################
simulate = "false" #set to "true" if want to simulate one economy
if simulate == "true"
	include("src/simulation_economy.jl")
	using simulation_economy

	#In two steps:
	#Simulate the model and store the results:
	model_ouput = simulation_economy.execute_simulation()

	#Analyse the economy:
	analysis_output = simulation_economy.analyse_economy(model_ouput, 10)
end

#Can be done in one step:
#Simulate_and_Analyse_Economy()

######################################
# E. Test what is driving inequalities
######################################
bargaining = "true" #set to "true" if want to test the bargaining power theory
if bargaining == "true"
	include("src/Bargaining_power_test.jl")
	using Bargaining_power_test

	Bargaining_power_test.create_and_analyse_bargaining_power()
end


text:
include("pci.jl")
include("pmu.jl")
include("imc.jl")
include("cha.jl")


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ 19befc40-7302-11ec-2216-a533c6a003d4
# hideall
let
    docs_dir = dirname(dirname(@__DIR__))
    pkg_dir = dirname(docs_dir)

    using Pkg: Pkg
    Pkg.activate(docs_dir)
    Pkg.develop(; path=pkg_dir)
    Pkg.instantiate()

    # Putting the include here to avoid Pluto getting confused about cell order.
    include(joinpath(docs_dir, "src", "tutorials_utils.jl"))
end;

# ╔═╡ 544e6e15-dce2-408b-aed2-b7aa7c16b669
using CSV

# ╔═╡ c2008c37-756f-43e1-beb1-9471ab054229
using DataFrames

# ╔═╡ a9b80331-2ca3-4c1e-beec-d3f7c9227e58
using TuringGLM

# ╔═╡ a8873af4-f4eb-4285-aabe-67abf9656570
md"""
For our example on **Negative Binomial Regression**, let's use a famous dataset called `roaches` (Gelman & Hill, 2007), which is data on the efficacy of a
pest management system at reducing the number of roaches in urban apartments.
It has 262 observations and the following variables:

* `y` -- number of roaches caught.
* `roach1` -- pretreatment number of roaches.
* `treatment` -- binary/dummy (0 or 1) for treatment indicator.
* `senior` -- binary/dummy (0 or 1) for only elderly residents in building.
* `exposure2` -- number of days for which the roach traps were used
"""

# ╔═╡ 5f29861c-cd0e-4f33-bf4e-3cb53faa1e7a
url = "https://github.com/TuringLang/TuringGLM.jl/raw/main/data/roaches.csv";

# ╔═╡ ac9084ba-9287-44df-843e-5ba33fde2055
roaches = CSV.read(download(url), DataFrame)

# ╔═╡ e715eda7-bcbf-4292-b95f-003816d47342
md"""
Using `y` as dependent variable and `roach1`, `treatment`, and `senior` as independent variables:
"""

# ╔═╡ fcec7550-4902-4907-a3a0-ef34b0f0cc21
fm = @formula(y ~ roach1 + treatment + senior)

# ╔═╡ 124aeb01-7661-4402-a2fa-77d2771b686c
md"""
We instantiate our model with `turing_model` passing a keyword argument `model=NegativeBinomial` to indicate that the model is a negative binomial regression:
"""

# ╔═╡ 7ea3a50d-3b0d-4cfd-b311-806d6ae59c1a
model = turing_model(fm, roaches; model=NegativeBinomial);

# ╔═╡ 597ce5d7-df3e-44e3-a154-e06e64894854
chn = sample(model, NUTS(), 2_000);

# ╔═╡ e8aec499-f48f-4e7f-93a6-ba38881ce147
# hide
plot_chains(chn)

# ╔═╡ abd3b2b0-af16-4895-9391-72cab5fda507
md"""
## References

Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge university press.
"""

# ╔═╡ Cell order:
# ╠═19befc40-7302-11ec-2216-a533c6a003d4
# ╠═a8873af4-f4eb-4285-aabe-67abf9656570
# ╠═544e6e15-dce2-408b-aed2-b7aa7c16b669
# ╠═c2008c37-756f-43e1-beb1-9471ab054229
# ╠═5f29861c-cd0e-4f33-bf4e-3cb53faa1e7a
# ╠═ac9084ba-9287-44df-843e-5ba33fde2055
# ╠═a9b80331-2ca3-4c1e-beec-d3f7c9227e58
# ╠═e715eda7-bcbf-4292-b95f-003816d47342
# ╠═fcec7550-4902-4907-a3a0-ef34b0f0cc21
# ╠═124aeb01-7661-4402-a2fa-77d2771b686c
# ╠═7ea3a50d-3b0d-4cfd-b311-806d6ae59c1a
# ╠═597ce5d7-df3e-44e3-a154-e06e64894854
# ╠═e8aec499-f48f-4e7f-93a6-ba38881ce147
# ╠═abd3b2b0-af16-4895-9391-72cab5fda507


text:
@testset "HW6" begin
    using DMUStudent
    using DMUStudent.HW6
    using POMDPs
    using POMDPTesting
    using QMDP
    using POMDPSimulators
    using POMDPPolicies
    import POMDPModelTools: render, Uniform
    using Compose
    using SARSOP

    small = LaserTagPOMDP(size=(4,3), n_obstacles=3)

    @test has_consistent_distributions(small)

    @test simulate(RolloutSimulator(max_steps=100), small, RandomPolicy(small)) isa Float64

    solve(QMDPSolver(), small)

    sp = LTState([3, 2], [2, 3])
    r = render(small, (o=[1,1,1,1], bp=Uniform(states(small)), sp=sp))
    filename = "/tmp/lasertag.svg"
    draw(SVG(filename), r)

    r = render(small, (done=true,))

    @test HW6.evaluate(RandomPolicy(LaserTagPOMDP())).score < 0.0
    @test HW6.evaluate(RandomPolicy(LaserTagPOMDP()), "zachary.sunberg@colorado.edu").score < 0.0

    tiny = LaserTagPOMDP(size=(2,2), n_obstacles=0)

    solve(SARSOPSolver(), tiny)
end


text:
using Documenter
using CoexistenceHoles

push!(LOAD_PATH,"../src/")

makedocs(
    sitename="CoexistenceHoles.jl",
    modules = [CoexistenceHoles],
    pages = Any[
        "Home" => "index.md"
        "Overview" => "overview.md"
        "User Guide" => [
            "examples.md",
            "tutorial.md",
            "glossary.md",
            "R_use.md"
            ]
    ],
)

deploydocs(
    repo = "github.com/SyntheticDynamics/CoexistenceHoles.jl.git"
)

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes