text:
@muladd function perform_step!(integrator,cache::SROCK1ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator

  maxeig!(integrator, cache)
  cache.mdeg = Int(floor(sqrt(2*dt*integrator.opts.internalnorm(integrator.eigen_est,t))+1)) # this is the spectral radius estimate to choose optimal stage
  choose_deg!(integrator,cache)

  mdeg = cache.mdeg
  η  = cache.optimal_η
  ω₀ = 1.0 + (η/(mdeg^2))
  ωSq = (ω₀^2) - 1.0
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  if alg_interpretation(integrator.alg) == :Stratonovich
    α  = cosh(mdeg*cosh_inv)/(2*ω₀*cosh((mdeg-1)*cosh_inv))
    γ  = 1/(2*α)
    β  = -γ
  end

  uᵢ₋₂ = copy(uprev)
  k = integrator.f(uprev,p,t)
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t+dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t
  gₘ₋₁ = zero(k)
  gₘ₋₂ = zero(k)

  #stage 1
  uᵢ₋₁ = uprev + (dt*ω₁/ω₀)*k

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*(Tᵢ₋₁/Tᵢ)
    ν = 2*ω₀*(Tᵢ₋₁/Tᵢ)
    κ = (-Tᵢ₋₂/Tᵢ)
    k = integrator.f(uᵢ₋₁,p,tᵢ₋₁)

    u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
    if (i > mdeg - 2) && alg_interpretation(integrator.alg) == :Stratonovich
      if i == mdeg - 1
        gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
        if typeof(W.dW) <: Number || !is_diagonal_noise(integrator.sol.prob)
          u += α*(gₘ₋₂*W.dW)
        else
          u .+= α .*  gₘ₋₂ .* W.dW
        end
      else
        gₘ₋₁ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
        if typeof(W.dW) <: Number || !is_diagonal_noise(integrator.sol.prob)
          u += β*(gₘ₋₂*W.dW) + γ*(gₘ₋₁*W.dW)
        else
          u .+= (β .* gₘ₋₂ .+ γ .* gₘ₋₁) .* W.dW
        end
      end
    elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
      if typeof(W.dW) <: Number
        gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
        uᵢ₋₂ = uᵢ₋₁ + sqrt(dt)*gₘ₋₂
        gₘ₋₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
        u += gₘ₋₂*W.dW + 1/(2.0*sqrt(dt))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - dt)
      elseif is_diagonal_noise(integrator.sol.prob)
        gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
        uᵢ₋₂ .= uᵢ₋₁ .+ sqrt(dt) .* gₘ₋₂
        gₘ₋₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
        u .+= gₘ₋₂ .* W.dW .+ (1/(2.0*sqrt(dt))) .* (gₘ₋₁ .- gₘ₋₂) .* (W.dW .^ 2 .- dt)
      else
        gₘ₋₂ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
        u += gₘ₋₂*W.dW
      end
    end

    if i < mdeg
      tᵢ = μ*dt + ν*tᵢ₋₁ + κ*tᵢ₋₂
      uᵢ₋₂ = uᵢ₋₁
      uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCK1Cache,f=integrator.f)
  @unpack uᵢ₋₁,uᵢ₋₂,k, gₘ₋₁, gₘ₋₂ = cache
  @unpack t,dt,uprev,u,W,p = integrator
  ccache = cache.constantcache
  maxeig!(integrator, cache)
  ccache.mdeg = Int(floor(sqrt(2*dt*integrator.opts.internalnorm(integrator.eigen_est,t))+1))   # this is the spectral radius estimate to choose optimal stage
  choose_deg!(integrator,cache)

  mdeg = ccache.mdeg
  η  = ccache.optimal_η
  ω₀ = 1 + η/(mdeg^2)
  ωSq = ω₀^2 - 1
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  if alg_interpretation(integrator.alg) == :Stratonovich
    α  = cosh(mdeg*cosh_inv)/(2*ω₀*cosh((mdeg-1)*cosh_inv))
    γ  = 1/(2*α)
    β  = -γ
  end

  @.. uᵢ₋₂ = uprev
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t + dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t

  #stage 1
  #this take advantage of the fact that cache.k === cache.fsalfirst
  #and this has already been done i maxeig!  i.e. integrator.f(fsalfirst, uprev, p, t)
  # integrator.f(k,uprev,p,t)
  @.. uᵢ₋₁ = uprev + (dt*ω₁/ω₀)*k

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*Tᵢ₋₁/Tᵢ
    ν = 2*ω₀*Tᵢ₋₁/Tᵢ
    κ = - Tᵢ₋₂/Tᵢ
    integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
    @.. u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
    if (i > mdeg - 2) && alg_interpretation(integrator.alg) == :Stratonovich
      if i == mdeg - 1
        integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
        if typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob)
          @.. u += α*gₘ₋₂*W.dW
        else
          mul!(k,gₘ₋₂,W.dW)
          @.. u += α*k
        end
      else
        integrator.g(gₘ₋₁,uᵢ₋₁,p,tᵢ₋₁)
        if typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob)
          @.. u += (β*gₘ₋₂ + γ*gₘ₋₁)*W.dW
        else
          mul!(k,gₘ₋₂,W.dW)
          @.. u += β*k
          mul!(k,gₘ₋₁,W.dW)
          @.. u += γ*k
        end
      end
    elseif (i == mdeg) && alg_interpretation(integrator.alg) == :Ito
      if typeof(W.dW) <: Number || is_diagonal_noise(integrator.sol.prob)
        integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
        @.. uᵢ₋₂ = uᵢ₋₁ + sqrt(dt)*gₘ₋₂
        integrator.g(gₘ₋₁,uᵢ₋₂,p,tᵢ₋₁)
        @.. u += gₘ₋₂*W.dW + 1/(2.0*sqrt(dt))*(gₘ₋₁ - gₘ₋₂)*(W.dW^2 - dt)
      else
        integrator.g(gₘ₋₂,uᵢ₋₁,p,tᵢ₋₁)
        mul!(uᵢ₋₁,gₘ₋₂,W.dW)
        @.. u += uᵢ₋₁
      end
    end

    if i < mdeg
      tᵢ = dt*μ + ν*tᵢ₋₁ + κ*tᵢ₋₂
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCK2ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, recf2, mα, mσ, mτ = cache

  gen_prob = !((is_diagonal_noise(integrator.sol.prob)) || (typeof(W.dW) <: Number) || (length(W.dW) == 1))
  gen_prob && (vec_χ = 2 .* floor.(false .* W.dW .+ 0.5 .+ rand(eltype(W.dW),length(W.dW))) .- 1.0)

  maxeig!(integrator, cache)
  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  cache.mdeg = max(3,min(cache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = cache.mdeg
  start     = cache.start
  deg_index = cache.deg_index
  α = mα[deg_index]
  σ = (1.0-α)*0.5 + α*mσ[deg_index]
  τ = 0.5*((1.0-α)^2) + 2*α*(1.0-α)*mσ[deg_index] + (α^2.0)*(mσ[deg_index]*(mσ[deg_index]+mτ[deg_index]))

  sqrt_dt   = sqrt(dt)

  μ = recf[start]  # here κ = 0
  tᵢ = t + α*dt*μ
  tᵢ₋₁ = tᵢ
  tᵢ₋₂ = t

  # stage 1
  uᵢ₋₂ = uprev
  uᵢ = integrator.f(uprev,p,t)
  uᵢ₋₁ = uprev + α*dt*μ*uᵢ

  # stages 2 upto s-2
  for i in 2:mdeg
    μ, κ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    ν    = 1.0 + κ
    uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
    uᵢ   = α*dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
    uᵢ₋₂ = uᵢ₋₁
    uᵢ₋₁ = uᵢ
    tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end

  #2 stage-finishing procedure
  #stage s-1
  μ, κ = recf2[(deg_index-1)*4 + 1], recf2[(deg_index-1)*4 + 2]
  ν    = 1.0 + κ
  uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)

  tₓ   = tᵢ₋₁ + 2*dt*τ                    # So that we don't have to calculate f(uₛ₋₂) again
  uₓ   = uᵢ₋₁ + 2*dt*τ*uᵢ                 # uₓ and tₓ represent u_star
  u    = uᵢ₋₁ + (2*σ - 0.5)*dt*uᵢ

  uᵢ   = α*dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
  tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
  tᵢ₋₂ = tᵢ₋₁
  tᵢ₋₁ = tᵢ
  uᵢ₋₂ = uᵢ₋₁
  uᵢ₋₁ = uᵢ

  #stage s
  μ, κ = recf2[(deg_index-1)*4 + 3], recf2[(deg_index-1)*4 + 4]
  ν    = 1.0 + κ
  uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
  uᵢ   = α*dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
  tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂

  # Now uᵢ₋₂ = uₛ₋₂, uᵢ₋₁ = uₛ₋₁, uᵢ = uₛ
  # Similarly tᵢ₋₂ = tₛ₋₂, tᵢ₋₁ = tₛ₋₁, tᵢ = tₛ

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    Gₛ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
    u  += Gₛ .* W.dW
    Gₛ = integrator.g(uᵢ,p,tᵢ)
    uₓ += Gₛ .* W.dW

    uₓ = integrator.f(uₓ,p,tₓ)
    u  += (1//2*dt)*uₓ
    uₓ  = 0.5 .* Gₛ .* (W.dW .^ 2 .- dt)
    uᵢ₋₂ = uᵢ + uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u   += (1//2)*Gₛ₁
    uᵢ₋₂ = uᵢ - uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u   -= (1//2)*Gₛ₁

    uₓ  = Gₛ*sqrt_dt
    uᵢ₋₂ = uᵢ + uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u += 1//4 .* W.dW .* (Gₛ₁ .- 2 .* Gₛ)
    uᵢ₋₂ = uᵢ - uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u  += 1//4 .* W.dW .* Gₛ₁
  else
    Gₛ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
    u += Gₛ*W.dW

    Gₛ = integrator.g(uᵢ,p,tᵢ)
    uₓ += Gₛ*W.dW

    uₓ = integrator.f(uₓ,p,tₓ)
    u  += (1//2)*dt*uₓ
    for i in 1:length(W.dW)
      WikJ = W.dW[i]; WikJ2 = vec_χ[i]
      WikRange = 1//2 .* (W.dW .* WikJ .- (1:length(W.dW) .== i) .* dt) #.- (1:length(W.dW) .> i) .* dt .* vec_χ .+ (1:length(W.dW) .< i) .* dt .* WikJ2)
      uₓ = Gₛ*WikRange
      WikRange = 0.5 .* (1:length(W.dW) .== i)
      uᵢ₋₂ = uᵢ + uₓ
      Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
      u   += (Gₛ₁*WikRange)
      uᵢ₋₂ = uᵢ - uₓ
      Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
      u   -= (Gₛ₁*WikRange)
    end

    uₓ = sqrt_dt*(Gₛ*vec_χ)
    uᵢ₋₂ = uᵢ + uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u += 1//4*(Gₛ₁*W.dW) - 1//2*(Gₛ*W.dW)

    uᵢ₋₂ = uᵢ - uₓ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ)
    u += 1//4*(Gₛ₁*W.dW)
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCK2Cache,f=integrator.f)
  @unpack uᵢ, uₓ, uᵢ₋₁, uᵢ₋₂, k, Gₛ, Gₛ₁, vec_χ, WikRange = cache

  @unpack t,dt,uprev,u,W,p = integrator

  @unpack recf, recf2, mα, mσ, mτ = cache.constantcache
  ccache = cache.constantcache
  gen_prob = !((is_diagonal_noise(integrator.sol.prob)) || (typeof(W.dW) <: Number) || (length(W.dW) == 1))

  maxeig!(integrator, cache)
  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  ccache.mdeg = max(3,min(ccache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = ccache.mdeg
  start     = ccache.start
  deg_index = ccache.deg_index
  α = mα[deg_index]
  σ = (1.0-α)*0.5 + α*mσ[deg_index]
  τ = 0.5*((1.0-α)^2) + 2*α*(1.0-α)*mσ[deg_index] + (α^2.0)*(mσ[deg_index]*(mσ[deg_index]+mτ[deg_index]))

  sqrt_dt   = sqrt(dt)
  if gen_prob
    vec_χ .= 0.5 .+ rand(eltype(W.dW),length(W.dW))
    @.. vec_χ = 2*floor(vec_χ) - 1.0
  end

  μ = recf[start]  # here κ = 0
  tᵢ = t + α*dt*μ
  tᵢ₋₁ = tᵢ
  tᵢ₋₂ = t

  # stage 1
  @.. uᵢ₋₂ = uprev
  #this take advantage of the fact that cache.k === cache.fsalfirst
  #and this has already been done i maxeig!  i.e. integrator.f(fsalfirst, uprev, p, t)
  # integrator.f(k,uprev,p,t)
  @.. uᵢ₋₁ = uprev + α*dt*μ*k

  # stages 2 upto s-2
  for i in 2:mdeg
    μ, κ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    ν    = 1.0 + κ
    integrator.f(k,uᵢ₋₁,p,t)
    @.. uᵢ   = α*dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
    @.. uᵢ₋₂ = uᵢ₋₁
    @.. uᵢ₋₁ = uᵢ
    tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end

  #2 stage-finishing procedure
  #stage s-1
  μ, κ = recf2[(deg_index-1)*4 + 1], recf2[(deg_index-1)*4 + 2]
  ν    = 1.0 + κ
  integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)

  tₓ   = tᵢ₋₁ + 2*dt*τ                    # So that we don't have to calculate f(uₛ₋₂) again
  @.. uₓ   = uᵢ₋₁ + 2*dt*τ*k                 # uₓ and tₓ represent u_star
  @.. u    = uᵢ₋₁ + (2*σ - 0.5)*dt*k

  @.. uᵢ   = α*dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
  tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
  tᵢ₋₂ = tᵢ₋₁
  tᵢ₋₁ = tᵢ
  @.. uᵢ₋₂ = uᵢ₋₁
  @.. uᵢ₋₁ = uᵢ

  #stage s
  μ, κ = recf2[(deg_index-1)*4 + 3], recf2[(deg_index-1)*4 + 4]
  ν    = 1.0 + κ
  integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
  @.. uᵢ   = α*dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
  tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂

  # Now uᵢ₋₂ = uₛ₋₂, uᵢ₋₁ = uₛ₋₁, uᵢ = uₛ
  # Similarly tᵢ₋₂ = tₛ₋₂, tᵢ₋₁ = tₛ₋₁, tᵢ = tₛ

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    integrator.g(Gₛ,uᵢ₋₁,p,tᵢ₋₁)
    @.. u  += Gₛ*W.dW
    integrator.g(Gₛ,uᵢ,p,tᵢ)
    @.. uₓ += Gₛ*W.dW
    integrator.f(k,uₓ,p,tₓ)
    @.. u  += (1//2)*dt*k

    @.. uₓ  = Gₛ*((W.dW^2 - dt)/2)
    @.. uᵢ₋₂ = uᵢ + uₓ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
    @.. u   += (1//2)*Gₛ₁
    @.. uᵢ₋₂ = uᵢ - uₓ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
    @.. u   -= (1//2)*Gₛ₁

    @.. uₓ = sqrt_dt*Gₛ
    @.. uᵢ₋₂ = uᵢ + uₓ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
    @.. u += (1//4*W.dW)*(Gₛ₁ - 2*Gₛ)
    @.. uᵢ₋₂ = uᵢ - uₓ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
    @.. u += (1//4*W.dW)*Gₛ₁
  else
      integrator.g(Gₛ,uᵢ₋₁,p,tᵢ₋₁)
      mul!(uᵢ₋₂,Gₛ,W.dW)
      @.. u += uᵢ₋₂

      integrator.g(Gₛ,uᵢ,p,tᵢ)
      mul!(uᵢ₋₁,Gₛ,W.dW)
      @.. uₓ += uᵢ₋₁

      integrator.f(k,uₓ,p,tₓ)
      @.. u  += (1//2)*dt*k

      for i in 1:length(W.dW)
        WikJ = W.dW[i]; WikJ2 = vec_χ[i]
        WikRange .= 1//2 .* (W.dW .* WikJ .- (1:length(W.dW) .== i) .* dt )#.+ (1:length(W.dW) .< i) .* dt .* WikJ2 .- (1:length(W.dW) .> i) .* dt .* vec_χ)
        mul!(uₓ,Gₛ,WikRange)
        @.. uᵢ₋₂ = uᵢ + uₓ
        WikRange .= 1//2 .* (1:length(W.dW) .== i)
        integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
        mul!(uᵢ₋₂,Gₛ₁,WikRange)
        @.. u   += uᵢ₋₂
        @.. uᵢ₋₂ = uᵢ - uₓ
        integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
        mul!(uᵢ₋₂,Gₛ₁,WikRange)
        @.. u   -= uᵢ₋₂
      end
      @.. WikRange = vec_χ*sqrt_dt
      mul!(uₓ,Gₛ,WikRange)

      @.. uᵢ₋₂ = uᵢ + uₓ
      integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
      mul!(uᵢ₋₂,Gₛ₁,W.dW)
      # mul!(uᵢ₋₁,Gₛ,W.dW) This is already been calculated
      @.. u += 1//4*(uᵢ₋₂ - 2*uᵢ₋₁)

      @.. uᵢ₋₂ = uᵢ - uₓ
      integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ)
      mul!(uᵢ₋₂,Gₛ₁,W.dW)
      @.. u += 1//4*uᵢ₋₂
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCKEMConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator

  maxeig!(integrator, cache)
  if integrator.alg.strong_order_1
    cache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
  else
    cache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
  end
  cache.mdeg = max(3,min(cache.mdeg,200))
  choose_deg!(integrator,cache)

  mdeg = cache.mdeg
  η  = cache.optimal_η
  ω₀ = one(eltype(u)) + (η/(mdeg^2))
  ωSq = ω₀^2 - one(eltype(u))
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  uᵢ₋₂ = copy(uprev)
  k = integrator.f(uprev,p,t)
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t+dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t

  #stage 1
  uᵢ₋₁ = uprev + (dt*ω₁/ω₀)*k

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*(Tᵢ₋₁/Tᵢ)
    ν = 2*ω₀*(Tᵢ₋₁/Tᵢ)
    κ = (-Tᵢ₋₂/Tᵢ)
    k = integrator.f(uᵢ₋₁,p,tᵢ₋₁)

    u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
    tᵢ = μ*dt + ν*tᵢ₋₁ + κ*tᵢ₋₂

    if i < mdeg
      uᵢ₋₂ = uᵢ₋₁
      uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end

  Gₛ = integrator.g(u,p,tᵢ)
  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    u += Gₛ .* W.dW
  else
    u += Gₛ*W.dW
  end

  if integrator.alg.strong_order_1
    if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || (is_diagonal_noise(integrator.sol.prob))
      uᵢ₋₂  = 1//2 .* Gₛ .* (W.dW .^ 2 .- dt)
      tmp = u + uᵢ₋₂
      Gₛ  = integrator.g(tmp,p,tᵢ)
      uᵢ₋₁ = 0.5*Gₛ
      tmp = u - uᵢ₋₂
      Gₛ  = integrator.g(tmp,p,tᵢ)
      uᵢ₋₁ -= 0.5*Gₛ
      u += uᵢ₋₁
    else
      for i in 1:length(W.dW)
        WikJ = W.dW[i]
        WikRange = 1//2 .* (W.dW .* WikJ - (1:length(W.dW) .== i) .* dt)
        uᵢ₋₂ = Gₛ*WikRange
        WikRange = 1//2 .* (1:length(W.dW) .== i)
        tmp = u + uᵢ₋₂
        Gₛ₁ = integrator.g(tmp,p,tᵢ)
        uᵢ₋₁ = Gₛ₁*WikRange
        tmp = u - uᵢ₋₂
        Gₛ₁ = integrator.g(tmp,p,tᵢ)
        uᵢ₋₁ -= Gₛ₁*WikRange
        u += uᵢ₋₁
      end
    end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCKEMCache,f=integrator.f)
  @unpack uᵢ₋₁,uᵢ₋₂,tmp,k,Gₛ,Gₛ₁,WikRange = cache
  @unpack t,dt,uprev,u,W,p = integrator
  ccache = cache.constantcache
  maxeig!(integrator, cache)
  if integrator.alg.strong_order_1
    ccache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.19)+1))
  else
    ccache.mdeg = Int(floor(sqrt(dt*integrator.opts.internalnorm(integrator.eigen_est,t)/0.33)+1))
  end
  ccache.mdeg = max(3,min(ccache.mdeg,200))
  choose_deg!(integrator,cache)

  mdeg = ccache.mdeg
  η  = ccache.optimal_η
  ω₀ = one(eltype(u)) + η/(mdeg^2)
  ωSq = ω₀^2 - one(eltype(u))
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  @.. uᵢ₋₂ = uprev
  integrator.f(k,uprev,p,t)
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t + dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t

  #stage 1
  @.. uᵢ₋₁ = uprev + (dt*ω₁/ω₀)*k

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*Tᵢ₋₁/Tᵢ
    ν = 2*ω₀*Tᵢ₋₁/Tᵢ
    κ = - Tᵢ₋₂/Tᵢ
    integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
    @.. u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
    tᵢ = dt*μ + ν*tᵢ₋₁ + κ*tᵢ₋₂

    if i < mdeg
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end

  integrator.g(Gₛ,u,p,tᵢ)
  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    @.. u += Gₛ*W.dW
  else
    mul!(uᵢ₋₁,Gₛ,W.dW)
    u += uᵢ₋₁
  end

  if integrator.alg.strong_order_1
    if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || (is_diagonal_noise(integrator.sol.prob))
      @.. uᵢ₋₂  = 1//2*Gₛ*(W.dW^2 - dt)
      @.. tmp = u + uᵢ₋₂
      integrator.g(Gₛ,tmp,p,tᵢ)
      @.. uᵢ₋₁ = 0.5*Gₛ
      @.. tmp = u - uᵢ₋₂
      integrator.g(Gₛ,tmp,p,tᵢ)
      @.. uᵢ₋₁ -= 0.5*Gₛ
      u += uᵢ₋₁
    else
      for i in 1:length(W.dW)
        WikJ = W.dW[i]
        WikRange .= 1//2 .* (WikJ .* W.dW .- (1:length(W.dW) .== i) .* dt)
        mul!(uᵢ₋₂,Gₛ,WikRange)
        WikRange .= 1//2 .* (1:length(W.dW) .== i)
        @.. tmp = u + uᵢ₋₂
        integrator.g(Gₛ₁,tmp,p,tᵢ)
        mul!(uᵢ₋₁,Gₛ₁,WikRange)
        @.. tmp = u - uᵢ₋₂
        integrator.g(Gₛ₁,tmp,p,tᵢ)
        @.. u += uᵢ₋₁
        mul!(uᵢ₋₁,Gₛ₁,WikRange)
        @.. u -= uᵢ₋₁
      end
    end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SKSROCKConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator

  maxeig!(integrator, cache)
  η = convert(typeof(t),0.05)
  mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
  mdeg = max(3,min(mdeg,200))

  ω₀ = 1.0 + (η/(mdeg^2))
  ωSq = (ω₀^2) - 1.0
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  μ, ν , κ = ω₁/ω₀, mdeg*ω₁/2, mdeg*ω₁/ω₀
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t+dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t

  #stage 1
  Gₛ = integrator.g(uprev,p,t)
  if (typeof(W.dW) <: Number) || !is_diagonal_noise(integrator.sol.prob)
    u = Gₛ*W.dW
  else
    u = Gₛ .* W.dW
  end

  if integrator.alg.post_processing
    uᵢ₋₁ = uprev + ν*u
    uᵢ₋₂ = integrator.f(uprev,p,t)
    uᵢ₋₁ = integrator.f(uᵢ₋₁,p,t)
    uᵢ₋₁ = uprev + (μ*dt)*uᵢ₋₁ + κ*u + cache.mα[mdeg-1]*dt*(uᵢ₋₁ - 2*uᵢ₋₂)
    uᵢ₋₂ = uprev - ν*u
    uᵢ₋₂ = integrator.f(uᵢ₋₂,p,t)
    uᵢ₋₁ += (cache.mα[mdeg-1]*dt)*uᵢ₋₂
  else
    uᵢ₋₁ = uprev + ν*u
    uᵢ₋₂ = integrator.f(uᵢ₋₁,p,t)
    uᵢ₋₁ = uprev + (μ*dt)*uᵢ₋₂ + κ*u
  end

  uᵢ₋₂ = uprev

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*(Tᵢ₋₁/Tᵢ)
    ν = 2*ω₀*(Tᵢ₋₁/Tᵢ)
    κ = (-Tᵢ₋₂/Tᵢ)
    u = integrator.f(uᵢ₋₁,p,tᵢ₋₁)

    u = dt*μ*u + ν*uᵢ₋₁ + κ*uᵢ₋₂
    tᵢ = μ*dt + ν*tᵢ₋₁ + κ*tᵢ₋₂

    if i < mdeg
      uᵢ₋₂ = uᵢ₋₁
      uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end
  if integrator.alg.post_processing && (t+dt >= integrator.sol.prob.tspan[2])
    Gₛ = integrator.g(u,p,tᵢ)
    if (typeof(W.dW) <: Number) || is_diagonal_noise(integrator.sol.prob)
      uᵢ₋₁ = Gₛ
    else
      WikRange = 1 .* (1:length(W.dW) .== 1)
      uᵢ₋₁ = Gₛ*WikRange
    end
    winc = rand()*6
    if winc < 1.0
      u -= (sqrt(3*dt)*ccache.mc[mdeg-1])*uᵢ₋₁
    elseif winc < 2
      u += (sqrt(3*dt)*ccache.mc[mdeg-1])*uᵢ₋₁
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SKSROCKCache,f=integrator.f)
  @unpack uᵢ₋₁,uᵢ₋₂,k, Gₛ, WikRange = cache
  @unpack t,dt,uprev,u,W,p = integrator

  ccache = cache.constantcache
  maxeig!(integrator, cache)
  η = convert(typeof(t),0.05)
  mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t) + 1.5)/(2-η*4/3))+1))
  mdeg = max(3,min(mdeg,200))

  ω₀ = 1.0 + (η/(mdeg^2))
  ωSq = (ω₀^2) - 1.0
  Sqrt_ω = sqrt(ωSq)
  cosh_inv = log(ω₀ + Sqrt_ω)             # arcosh(ω₀)
  ω₁ = (Sqrt_ω*cosh(mdeg*cosh_inv))/(mdeg*sinh(mdeg*cosh_inv))

  μ, ν , κ = ω₁/ω₀, mdeg*ω₁/2, mdeg*ω₁/ω₀
  Tᵢ₋₂ = one(eltype(u))
  Tᵢ₋₁ = convert(eltype(u),ω₀)
  Tᵢ   = Tᵢ₋₁
  tᵢ₋₁ = t+dt*(ω₁/ω₀)
  tᵢ   = tᵢ₋₁
  tᵢ₋₂ = t

  #stage 1
  integrator.g(Gₛ,uprev,p,t)
  if (typeof(W.dW) <: Number) || is_diagonal_noise(integrator.sol.prob)
    @.. u = Gₛ*W.dW
  else
    mul!(u,Gₛ,W.dW)
  end

  if integrator.alg.post_processing
    @.. uᵢ₋₂ = uprev + ν*u
    integrator.f(k,uᵢ₋₂,p,t)
    @.. uᵢ₋₁ = uprev + (μ*dt)*k + κ*u + (ccache.mα[mdeg-1]*dt)*k
    integrator.f(k,uprev,p,t)
    @.. uᵢ₋₁ -= (ccache.mα[mdeg-1]*dt*2)*k
    @.. uᵢ₋₂ = uprev - ν*u
    integrator.f(k,uᵢ₋₂,p,t)
    @.. uᵢ₋₁ += (ccache.mα[mdeg-1]*dt)*k
  else
    @.. uᵢ₋₁ = uprev + ν*u
    integrator.f(k,uᵢ₋₁,p,t)
    @.. uᵢ₋₁ = uprev + (μ*dt)*k + κ*u
  end

  @.. uᵢ₋₂ = uprev

  for i in 2:mdeg
    Tᵢ = 2*ω₀*Tᵢ₋₁ - Tᵢ₋₂
    μ = 2*ω₁*(Tᵢ₋₁/Tᵢ)
    ν = 2*ω₀*(Tᵢ₋₁/Tᵢ)
    κ = (-Tᵢ₋₂/Tᵢ)
    integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)

    @.. u = dt*μ*k + ν*uᵢ₋₁ + κ*uᵢ₋₂
    tᵢ = μ*dt + ν*tᵢ₋₁ + κ*tᵢ₋₂

    if i < mdeg
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = u
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
      Tᵢ₋₂ = Tᵢ₋₁
      Tᵢ₋₁ = Tᵢ
    end
  end

  if integrator.alg.post_processing && (t+dt >= integrator.sol.prob.tspan[2])
    integrator.g(Gₛ,u,p,tᵢ)
    # println(Gₛ/length(W.dW))
    if (typeof(W.dW) <: Number) || is_diagonal_noise(integrator.sol.prob)
      @.. uᵢ₋₁ = Gₛ
    else
      WikRange .= 1 .* (1:length(W.dW) .== 1)
      mul!(uᵢ₋₁,Gₛ,WikRange)
    end
    winc = rand()*6
    if winc < 1.0
      @.. u -= (sqrt(3*dt)*ccache.mc[mdeg-1])*uᵢ₋₁
    elseif winc < 2
      @.. u += (sqrt(3*dt)*ccache.mc[mdeg-1])*uᵢ₋₁
    end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::TangXiaoSROCK2ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, recf2, mα, mσ, mτ, mn̂, c1, c2 = cache

  n̂ = mn̂[integrator.alg.version_num]

  maxeig!(integrator, cache)
  (integrator.alg.version_num <= 2) && (cache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
  (integrator.alg.version_num > 2) && (cache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))

  cache.mdeg = max(4,min(cache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = cache.mdeg
  start     = cache.start
  deg_index = cache.deg_index
  α = mα[integrator.alg.version_num]
  σ = (1.0-α)*0.5 + α*mσ[deg_index]
  τ = 0.5*((1.0-α)^2) + 2*α*(1.0-α)*mσ[deg_index] + (α^2.0)*(mσ[deg_index]*(mσ[deg_index]+mτ[deg_index]))

  η₁ = (rand() < 0.5) ? -1 : 1
  η₂ = (rand() < 0.5) ? -1 : 1
  sqrt_dt   = sqrt(dt)

  Û₁ = zero(u)
  Û₂ = zero(u)
  t̂₁ = t̂₂ = zero(t)
  tᵢ =  tᵢ₋₁ = tᵢ₋₂ = tₓ = t
  uᵢ₋₂ = uprev

  for i in 0:mdeg+1
    if i == 1
      μ = recf[start]
      tᵢ = tᵢ₋₁ = t + α*dt*μ

      uᵢ = integrator.f(uprev,p,t)
      uᵢ₋₁ = uprev + α*dt*μ*uᵢ
    elseif i > 1 && i <= mdeg
      μ, ν, κ = recf[start + 2*(i-2) + 1], 1.0 + recf[start + 2*(i-2) + 2], recf[start + 2*(i-2) + 2]

      uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
      uᵢ   = α*dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
      tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    elseif i == mdeg + 1
      μ, ν, κ = recf2[(deg_index-1)*4 + 1], 1 + recf2[(deg_index-1)*4 + 2], recf2[(deg_index-1)*4 + 2]

      uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)

      tₓ   = tᵢ₋₁ + 2*τ*dt
      uₓ   = uᵢ₋₁ + (2*τ*dt)*uᵢ
      u    = uᵢ₋₁ + (2*σ - 1//2)*dt*uᵢ

      uᵢ   = α*dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
      tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    end

    j = i - mdeg - 1 + n̂
    if j > 0
      j += cache.start_mcs - 1
      if i == 0
        Û₁ += c1[j]*uprev
        t̂₁ += c1[j]*t
        Û₂ += c2[j]*uprev
        t̂₂ += c2[j]*t
      elseif i  == 1
        Û₁ += c1[j]*uᵢ₋₁
        t̂₁ += c1[j]*tᵢ₋₁
        Û₂ += c2[j]*uᵢ₋₁
        t̂₂ += c2[j]*tᵢ₋₁
      else
        Û₁ += c1[j]*uᵢ
        t̂₁ += c1[j]*tᵢ
        Û₂ += c2[j]*uᵢ
        t̂₂ += c2[j]*tᵢ
      end
    end
    if i > 1 && i < mdeg + 1
      uᵢ₋₂ = uᵢ₋₁
      uᵢ₋₁ = uᵢ
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
    end
  end

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1)
    Gₛ = integrator.g(Û₁,p,t̂₁)
    uₓ += Gₛ*W.dW

    uₓ = integrator.f(uₓ,p,tₓ)
    u  += (1//2*dt)*uₓ + Gₛ*((W.dW^2 - dt)/(η₁*sqrt_dt) - W.dW)
    Û₁ -= (η₁*sqrt_dt/2)*Gₛ
    Û₂ += (η₁*sqrt_dt/2)*Gₛ

    Gₛ = integrator.g(Û₂,p,t̂₂)
    u += Gₛ*W.dW

    Gₛ = integrator.g(Û₁,p,t̂₁)
    u += Gₛ*(W.dW - (W.dW^2 - dt)/(η₁*sqrt_dt))
  elseif is_diagonal_noise(integrator.sol.prob)

    Gₛ = integrator.g(Û₁,p,t̂₁)
    uᵢ₋₁ = Gₛ .* W.dW

    Û₁ -= (1//2*η₁*sqrt_dt)*Gₛ
    Û₂ += (1//2*η₁*sqrt_dt)*Gₛ

    uₓ += uᵢ₋₁
    uₓ = integrator.f(uₓ,p,tₓ)
    u  += (1//2)*dt*uₓ

    u .+= Gₛ .* ((W.dW .^ 2 .- dt) ./ (η₁*sqrt_dt) .- W.dW)

    Gₛ = integrator.g(Û₂,p,t̂₂)
    u .+= Gₛ .* W.dW

    Gₛ = integrator.g(Û₁,p,t̂₁)
    u .-= Gₛ .* ((W.dW .^ 2 .- dt) ./ (η₁*sqrt_dt) .- W.dW)
  else
      Gₛ = integrator.g(Û₁,p,t̂₁)

      for i in 1:length(W.dW)
        (i == 1) && (uᵢ₋₁ = @view(Gₛ[:,i])*W.dW[i])
        (i != 1) && (uᵢ₋₁ += @view(Gₛ[:,i])*W.dW[i])
      end

      uₓ += uᵢ₋₁
      uₓ = integrator.f(uₓ,p,tₓ)

      u  += (1//2*dt)*uₓ - uᵢ₋₁

      for i in 1:length(W.dW)
        uᵢ₋₁ = Û₁ - (1//2*η₁*sqrt_dt)*@view(Gₛ[:,i])
        Gₛ₁ = integrator.g(uᵢ₋₁,p,t̂₁)
        u += @view(Gₛ₁[:,i])*W.dW[i] + (@view(Gₛ[:,i]) - @view(Gₛ₁[:,i]))*((W.dW[i]^2 - dt)/(η₁*sqrt_dt))
      end

      for i in 1:length(W.dW)
        for j in 1:length(W.dW)
          (i > j) && (WikJ = (1//2)*(1+η₂)*W.dW[j])
          (i < j) && (WikJ = (1//2)*(1-η₂)*W.dW[j])
          (i == j) && (WikJ = (1//2)*(η₁*sqrt_dt))

          uᵢ₋₁ += @view(Gₛ[:,j])*WikJ
        end
        Gₛ₁ = integrator.g(uᵢ₋₁,p,t̂₂)
        u += @view(Gₛ₁[:,i])*W.dW[i]
      end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::TangXiaoSROCK2Cache,f=integrator.f)
  @unpack uᵢ, uₓ, uᵢ₋₁, uᵢ₋₂, Û₁, Û₂, k, Gₛ, Gₛ₁ = cache
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, recf2, mα, mσ, mτ, mn̂, c1, c2 = cache.constantcache

  n̂ = mn̂[integrator.alg.version_num]
  ccache = cache.constantcache

  maxeig!(integrator, cache)
  (integrator.alg.version_num <= 2) && (ccache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1)))
  (integrator.alg.version_num > 2) && (ccache.mdeg = Int(floor(sqrt((dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.611)+1)))
  ccache.mdeg = max(4,min(ccache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = ccache.mdeg
  start     = ccache.start
  deg_index = ccache.deg_index
  α = convert(eltype(u),1.33)
  σ = (1.0-α)*0.5 + α*mσ[deg_index]
  τ = 0.5*((1.0-α)^2) + 2*α*(1.0-α)*mσ[deg_index] + (α^2.0)*(mσ[deg_index]*(mσ[deg_index]+mτ[deg_index]))

  η₁ = (rand() < 0.5) ? -1 : 1
  η₂ = (rand() < 0.5) ? -1 : 1
  sqrt_dt   = sqrt(dt)

  @.. Û₁ = zero(eltype(u))
  @.. Û₂ = zero(eltype(u))
  t̂₁ = t̂₂ = tₓ = zero(t)
  tᵢ =  tᵢ₋₁ = tᵢ₋₂ = t

  for i in 0:mdeg+1
    if i == 1
      μ = recf[start]
      tᵢ = tᵢ₋₁ = t + α*dt*μ

      @.. uᵢ₋₂ = uprev
      integrator.f(k,uprev,p,t)
      @.. uᵢ₋₁ = uprev + α*dt*μ*k
    elseif i > 1 && i <= mdeg
      μ, ν, κ = recf[start + 2*(i-2) + 1], 1.0 + recf[start + 2*(i-2) + 2], recf[start + 2*(i-2) + 2]

      integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
      @.. uᵢ   = α*dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
      tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    elseif i == mdeg + 1
      μ, ν, κ = recf2[(deg_index-1)*4 + 1], 1.0 + recf2[(deg_index-1)*4 + 2], recf2[(deg_index-1)*4 + 2]

      integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)

      tₓ   = tᵢ₋₁ + 2*τ*dt
      @.. uₓ   = uᵢ₋₁ + (2*τ*dt)*k
      @.. u    = uᵢ₋₁ + (2*σ - 1//2)*dt*k

      @.. uᵢ   = α*dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
      tᵢ   = α*dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    end

    j = i - mdeg - 1 + n̂
    if j > 0
      j += ccache.start_mcs - 1
      if i == 0
        @.. Û₁ += c1[j]*uprev
        t̂₁ += c1[j]*t
        @.. Û₂ += c2[j]*uprev
        t̂₂ += c2[j]*t
      elseif i  == 1
        @.. Û₁ += c1[j]*uᵢ₋₁
        t̂₁ += c1[j]*tᵢ₋₁
        @.. Û₂ += c2[j]*uᵢ₋₁
        t̂₂ += c2[j]*tᵢ₋₁
      else
        @.. Û₁ += c1[j]*uᵢ
        t̂₁ += c1[j]*tᵢ
        @.. Û₂ += c2[j]*uᵢ
        t̂₂ += c2[j]*tᵢ
      end
    end

    if i > 1 && i < mdeg + 1
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = uᵢ
      tᵢ₋₂ = tᵢ₋₁
      tᵢ₋₁ = tᵢ
    end
  end

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    integrator.g(Gₛ,Û₁,p,t̂₁)
    @.. uₓ += Gₛ*W.dW

    integrator.f(k,uₓ,p,tₓ)
    @.. u  += (1//2*dt)*k + Gₛ*((W.dW^2 - dt)/(η₁*sqrt_dt) - W.dW)
    @.. Û₁ -= (η₁*sqrt_dt/2)*Gₛ
    @.. Û₂ += (η₁*sqrt_dt/2)*Gₛ

    integrator.g(Gₛ,Û₂,p,t̂₂)
    @.. u += Gₛ*W.dW

    integrator.g(Gₛ,Û₁,p,t̂₁)
    @.. u += Gₛ*(W.dW - (W.dW^2 - dt)/(η₁*sqrt_dt))
  else
      integrator.g(Gₛ,Û₁,p,t̂₁)

      for i in 1:length(W.dW)
        (i == 1) && (@.. uᵢ₋₁ = @view(Gₛ[:,i])*W.dW[i])
        (i > 1) && (@.. uᵢ₋₁ += @view(Gₛ[:,i])*W.dW[i])
      end

      @.. uₓ += uᵢ₋₁
      integrator.f(k,uₓ,p,tₓ)

      @.. u  += (1//2*dt)*k - uᵢ₋₁

      for i in 1:length(W.dW)
        @.. uᵢ₋₁ = Û₁ - (1//2*η₁*sqrt_dt)*@view(Gₛ[:,i])
        integrator.g(Gₛ₁,uᵢ₋₁,p,t̂₁)
        @.. u += @view(Gₛ₁[:,i])*W.dW[i] + (@view(Gₛ[:,i]) - @view(Gₛ₁[:,i]))*((W.dW[i]^2 - dt)/(η₁*sqrt_dt))
      end

      for i in 1:length(W.dW)
        for j in 1:length(W.dW)
          (i > j) && (WikJ = (1//2)*(1+η₂)*W.dW[j])
          (i < j) && (WikJ = (1//2)*(1-η₂)*W.dW[j])
          (i == j) && (WikJ = (1//2)*(η₁*sqrt_dt))

          @.. uᵢ₋₁ += @view(Gₛ[:,j])*WikJ
        end
        @.. uᵢ₋₁ += Û₂
        integrator.g(Gₛ₁,uᵢ₋₁,p,t̂₂)
        @.. u += @view(Gₛ₁[:,i])*W.dW[i]
      end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::KomBurSROCK2ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, mσ, mτ, mδ = cache

  gen_prob = !((is_diagonal_noise(integrator.sol.prob)) || (typeof(W.dW) <: Number) || (length(W.dW) == 1))

  maxeig!(integrator, cache)
  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  cache.mdeg = max(6,min(cache.mdeg,200))-2
  choose_deg!(integrator,cache)

  # here mdeg == s in the paper
  mdeg      = cache.mdeg + 2
  start     = cache.start
  deg_index = cache.deg_index
  σ = mσ[deg_index]
  τ = mτ[deg_index]

  sqrt_dt   = sqrt(dt)
  (gen_prob) && (vec_χ = 2 .* floor.( 0.5 .+ false .* W.dW .+ rand(length(W.dW))) .- 1.0)

  tᵢ₋₂ = t
  uᵢ₋₂ = uprev
  μ = recf[start]
  tᵢ = tᵢ₋₁ = t + dt*μ
  u = integrator.f(uprev,p,t)
  uᵢ₋₁ = uprev + dt*μ*u

  for i in 2:mdeg - 4
    μ, θ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    u  = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
    u  = dt*μ*u + (1 + θ)*uᵢ₋₁ - θ*uᵢ₋₂
    tᵢ = dt*μ + (1 + θ)*tᵢ₋₁ - θ*tᵢ₋₂

    uᵢ₋₂ = uᵢ₋₁
    uᵢ₋₁ = u
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end

  μₛ₋₃ = dt*recf[start + 2*(mdeg - 5) + 1]
  θₛ₋₃ = recf[start + 2*(mdeg - 5) + 2]
  μₛ₋₂ = dt*recf[start + 2*(mdeg - 4) + 1]
  θₛ₋₂ = recf[start + 2*(mdeg - 4) + 2]

  δ₁ = dt*mδ[(deg_index - 1)*8 + 1]; δ₂ = dt*mδ[(deg_index - 1)*8 + 2]; δ₃ = dt*mδ[(deg_index - 1)*8 + 3]
  δ₄ = mδ[(deg_index - 1)*8 + 4]; δ₅ = mδ[(deg_index - 1)*8 + 5]; δ₆ = mδ[(deg_index - 1)*8 + 6]
  δ₇ = mδ[(deg_index - 1)*8 + 7]; δ₈ = mδ[(deg_index - 1)*8 + 8]

  C₁ = (1 + θₛ₋₂)*μₛ₋₃

  u = uᵢ₋₁ + θₛ₋₃*(1 + θₛ₋₂)*(uᵢ₋₁ - uᵢ₋₂)
  uᵢ₋₁ += θₛ₋₃*(uᵢ₋₁ - uᵢ₋₂)
  uᵢ₋₂ = u
  ttmp = tᵢ₋₁ + θₛ₋₃*(1 + θₛ₋₂)*(tᵢ₋₁ - tᵢ₋₂)
  tᵢ₋₁ += θₛ₋₃*(tᵢ₋₁ - tᵢ₋₂)
  tᵢ₋₂ = ttmp

  if typeof(W.dW) <: Number || length(W.dW) == 1 || is_diagonal_noise(integrator.sol.prob)
    # stage s-3
    yₛ₋₃ = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    Xₛ₋₃ = integrator.g(utmp,p,ttmp)
    u += C₁*yₛ₋₃ + 1//8 .* W.dW .* Xₛ₋₃

    #stage s-2
    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₄ .* W.dW .* Xₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    yₛ₋₂ = integrator.f(utmp,p,ttmp)

    utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3 .* W.dW .* Xₛ₋₃
    ttmp = tᵢ₋₂ + C₁
    Xₛ₋₂ = integrator.g(utmp,p,ttmp)
    u += μₛ₋₂*yₛ₋₂ + 3//8 .* W.dW .* Xₛ₋₂

    #stage s-1
    utmp = uᵢ₋₂ + μₛ₋₂*yₛ₋₂ + C₁*yₛ₋₃ + δ₅ .* W.dW .* Xₛ₋₃ + δ₄ .* W.dW .* Xₛ₋₂
    ttmp = tᵢ₋₂ + μₛ₋₂ + C₁
    yₛ₋₁ = integrator.f(utmp,p,ttmp)

    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ + 1//12 .* W.dW .* Xₛ₋₃ + 1//4 .* W.dW .* Xₛ₋₂
    ttmp = tᵢ₋₁ + μₛ₋₃ + δ₁
    Xₛ₋₁ = integrator.g(utmp,p,ttmp)
    u += (σ - τ)*dt*yₛ₋₁ + 3//8 .* W.dW .* Xₛ₋₁

    #stage s
    utmp = uᵢ₋₂ + C₁*yₛ₋₃ + μₛ₋₂*yₛ₋₂ + σ*dt*yₛ₋₁ + δ₆ .* W.dW .* Xₛ₋₃ + δ₇ .* W.dW .* Xₛ₋₂ + δ₈ .* W.dW .* Xₛ₋₁
    ttmp = tᵢ₋₂ + C₁ + μₛ₋₂ + σ*dt
    utmp = integrator.f(utmp,p,ttmp)
    u += (σ + τ)*dt*utmp

    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ - 5//4 .* W.dW .* Xₛ₋₃ + 1//4 .* W.dW .* Xₛ₋₂ + 2 .* W.dW .* Xₛ₋₁
    ttmp = tᵢ₋₁ + μₛ₋₃ + δ₂ + δ₃
    Xₛ₋₁ = integrator.g(utmp,p,ttmp)
    u += 1//8 .* W.dW .* Xₛ₋₁
  else
    # stage s-3
    yₛ₋₃ = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    Xₛ₋₃ = integrator.g(utmp,p,ttmp)
    SXₛ₋₃ = Xₛ₋₃*W.dW

    u += C₁*yₛ₋₃ + 1//8*SXₛ₋₃

    #stage s-2
    utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₄*SXₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    yₛ₋₂ = integrator.f(utmp,p,ttmp)
    Xₛ₋₂ = zero(Xₛ₋₃)

    for i in 1:length(W.dW)
      WikRange = W.dW .* (1:length(W.dW) .== i)
      # utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3*@view(Xₛ₋₃[:,i])*W.dW[i]
      utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3*(Xₛ₋₃*WikRange)
      ttmp = tᵢ₋₂ + C₁
      Gₛ = integrator.g(utmp,p,ttmp)
      WikRange = 1 .* (1:length(W.dW) .== i)
      # @view(Xₛ₋₂[:,i]) .=  @view(Gₛ[:,i])
      Xₛ₋₂ .+=  Gₛ .* WikRange
    end
    SXₛ₋₂ = Xₛ₋₂*W.dW
    u += μₛ₋₂*yₛ₋₂ + 3//8*SXₛ₋₂

    #stage s-1
    utmp = uᵢ₋₂ + μₛ₋₂*yₛ₋₂ + C₁*yₛ₋₃ + δ₅*SXₛ₋₃ + δ₄*SXₛ₋₂
    ttmp = tᵢ₋₂ + μₛ₋₂ + C₁
    yₛ₋₁ = integrator.f(utmp,p,ttmp)

    Xₛ₋₁ = zero(Xₛ₋₂)
    for i in 1:length(W.dW)
      WikRange = W.dW .* (1:length(W.dW) .== i)
      # utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ - 1//6*W.dW[i]*@view(Xₛ₋₃[:,i]) - 1//2*W.dW[i]*@view(Xₛ₋₂[:,i]) + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂ - 1//6*(Xₛ₋₃*WikRange) - 1//2*(Xₛ₋₂*WikRange)
      ttmp = tᵢ₋₁ + μₛ₋₃ + δ₁
      Gₛ = integrator.g(utmp,p,ttmp)
      WikRange = 1 .* (1:length(W.dW) .== i)
      #@view(Xₛ₋₁[:,i]) .= @view(Gₛ[:,i])
      Xₛ₋₁ .+= Gₛ .* WikRange
    end
    SXₛ₋₁ = Xₛ₋₁*W.dW
    u += (σ - τ)*dt*yₛ₋₁ + 3//8*SXₛ₋₁

    #stage s
    utmp = uᵢ₋₂ + C₁*yₛ₋₃ + μₛ₋₂*yₛ₋₂ + σ*dt*yₛ₋₁ + δ₆*SXₛ₋₃ + δ₇*SXₛ₋₂ + δ₈*SXₛ₋₁
    ttmp = tᵢ₋₂ + C₁ + μₛ₋₂ + σ*dt
    utmp = integrator.f(utmp,p,ttmp)
    u += (σ + τ)*dt*utmp

    SXₛ₋₁ = zero(uprev)
    for i in 1:length(W.dW)
      WikRange = W.dW .* (1:length(W.dW) .== i)
      # utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ - 3//2*W.dW[i]*@view(Xₛ₋₃[:,i]) - 1//2*W.dW[i]*@view(Xₛ₋₂[:,i]) + 2*W.dW[i]*@view(Xₛ₋₁[:,i]) + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂ - 3//2*(Xₛ₋₃*WikRange) - 1//2*(Xₛ₋₂*WikRange) + 2*(Xₛ₋₁*WikRange)
      ttmp = tᵢ₋₁ + μₛ₋₃ + δ₂ + δ₃
      Gₛ = integrator.g(utmp,p,ttmp)
      # SXₛ₋₁ += W.dW[i]*@view(Gₛ[:,i])
      SXₛ₋₁ += Gₛ*WikRange
    end
    u += 1//8*SXₛ₋₁

    for i in 1:length(W.dW)
      SXₛ₋₁ = zero(uprev)
      for j in 1:length(W.dW)
        WikRange = 1 .* (1:length(W.dW) .== j)
        if j != i
          # SXₛ₋₁ += (i > j ? -vec_χ[i]*W.dW[j] : vec_χ[j]*W.dW[i] )*@view(Xₛ₋₃[:,j])
          SXₛ₋₁ += (i > j ? -vec_χ[i]*W.dW[j] : vec_χ[j]*W.dW[i] )*(Xₛ₋₃*WikRange)
        end
      end
      ttmp = tᵢ₋₁ + μₛ₋₃
      utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ - 1//4*SXₛ₋₁
      Gₛ = integrator.g(utmp,p,ttmp)
      utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + 1//4*SXₛ₋₁
      Xₛ₋₁ = integrator.g(utmp,p,ttmp)
      sqrt_dt *= (length(W.dW) - 1)
      WikRange = 1 .* (1:length(W.dW) .== i)
      # u += sqrt_dt*(@view(Xₛ₋₁[:,i]) - @view(Gₛ[:,i]))
      u += sqrt_dt*(Xₛ₋₁*WikRange - Gₛ*WikRange)
    end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::KomBurSROCK2Cache,f=integrator.f)
  @unpack utmp, uᵢ₋₁, uᵢ₋₂, k, yₛ₋₁, yₛ₋₂, yₛ₋₃, SXₛ₋₁, SXₛ₋₂, SXₛ₋₃, Gₛ, Xₛ₋₁, Xₛ₋₂, Xₛ₋₃, vec_χ = cache
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, mσ, mτ, mδ = cache.constantcache

  ccache = cache.constantcache
  gen_prob = !((is_diagonal_noise(integrator.sol.prob)) || (typeof(W.dW) <: Number) || (length(W.dW) == 1))

  maxeig!(integrator, cache)
  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  ccache.mdeg = max(6,min(ccache.mdeg,200))-2
  choose_deg!(integrator,cache)

  # here mdeg == s in the paper
  mdeg      = ccache.mdeg + 2
  start     = ccache.start
  deg_index = ccache.deg_index
  σ = mσ[deg_index]
  τ = mτ[deg_index]

  sqrt_dt   = sqrt(dt)
  (gen_prob) && (vec_χ .= 2 .* floor.(0.5 .+ false .* vec_χ .+ rand(length(vec_χ))) .- 1.0)

  tᵢ₋₂ = t
  @.. uᵢ₋₂ = uprev
  μ = recf[start]
  tᵢ = tᵢ₋₁ = t + dt*μ
  integrator.f(k,uprev,p,t)
  @.. uᵢ₋₁ = uprev + dt*μ*k

  for i in 2:mdeg-4
    μ, θ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
    @.. u  = dt*μ*k + (1 + θ)*uᵢ₋₁ - θ*uᵢ₋₂
    tᵢ = dt*μ + (1 + θ)*tᵢ₋₁ - θ*tᵢ₋₂

    @.. uᵢ₋₂ = uᵢ₋₁
    @.. uᵢ₋₁ = u
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end
  μₛ₋₃ = dt*recf[start + 2*(mdeg - 5) + 1]
  θₛ₋₃ = recf[start + 2*(mdeg - 5) + 2]
  μₛ₋₂ = dt*recf[start + 2*(mdeg - 4) + 1]
  θₛ₋₂ = recf[start + 2*(mdeg - 4) + 2]

  δ₁ = dt*mδ[(deg_index - 1)*8 + 1]; δ₂ = dt*mδ[(deg_index - 1)*8 + 2]; δ₃ = dt*mδ[(deg_index - 1)*8 + 3]
  δ₄ = mδ[(deg_index - 1)*8 + 4]; δ₅ = mδ[(deg_index - 1)*8 + 5]; δ₆ = mδ[(deg_index - 1)*8 + 6]
  δ₇ = mδ[(deg_index - 1)*8 + 7]; δ₈ = mδ[(deg_index - 1)*8 + 8]

  C₁ = (1 + θₛ₋₂)*μₛ₋₃

  @.. u = uᵢ₋₁ + θₛ₋₃*(1 + θₛ₋₂)*(uᵢ₋₁ - uᵢ₋₂)
  @.. uᵢ₋₁ += θₛ₋₃*(uᵢ₋₁ - uᵢ₋₂)
  @.. uᵢ₋₂ = u
  ttmp = tᵢ₋₁ + θₛ₋₃*(1 + θₛ₋₂)*(tᵢ₋₁ - tᵢ₋₂)
  tᵢ₋₁ += θₛ₋₃*(tᵢ₋₁ - tᵢ₋₂)
  tᵢ₋₂ = ttmp

  if typeof(W.dW) <: Number || length(W.dW) == 1 || is_diagonal_noise(integrator.sol.prob)
    # stage s-3
    integrator.f(yₛ₋₃,uᵢ₋₁,p,tᵢ₋₁)
    @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    integrator.g(Xₛ₋₃,utmp,p,ttmp)
    @.. u += C₁*yₛ₋₃ + 1//8*W.dW*Xₛ₋₃

    #stage s-2
    @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₄*W.dW*Xₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    integrator.f(yₛ₋₂,utmp,p,ttmp)

    @.. utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3*W.dW*Xₛ₋₃
    ttmp = tᵢ₋₂ + C₁
    integrator.g(Xₛ₋₂,utmp,p,ttmp)
    @.. u += μₛ₋₂*yₛ₋₂ + 3//8*W.dW*Xₛ₋₂

    #stage s-1
    @.. utmp = uᵢ₋₂ + μₛ₋₂*yₛ₋₂ + C₁*yₛ₋₃ + δ₅*W.dW*Xₛ₋₃ + δ₄*W.dW*Xₛ₋₂
    ttmp = tᵢ₋₂ + μₛ₋₂ + C₁
    integrator.f(yₛ₋₁,utmp,p,ttmp)

    @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ + 1//12*W.dW*Xₛ₋₃ + 1//4*W.dW*Xₛ₋₂
    ttmp = tᵢ₋₁ + μₛ₋₃ + δ₁
    integrator.g(Xₛ₋₁,utmp,p,ttmp)
    @.. u += (σ - τ)*dt*yₛ₋₁ + 3//8*W.dW*Xₛ₋₁

    #stage s
    @.. utmp = uᵢ₋₂ + C₁*yₛ₋₃ + μₛ₋₂*yₛ₋₂ + σ*dt*yₛ₋₁ + δ₆*W.dW*Xₛ₋₃ + δ₇*W.dW*Xₛ₋₂ + δ₈*W.dW*Xₛ₋₁
    ttmp = tᵢ₋₂ + C₁ + μₛ₋₂ + σ*dt

    @.. uᵢ₋₁ = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ - 5//4*W.dW*Xₛ₋₃ + 1//4*W.dW*Xₛ₋₂ + 2*W.dW*Xₛ₋₁
    tᵢ₋₁ = tᵢ₋₁ + μₛ₋₃ + δ₂ + δ₃

    integrator.f(yₛ₋₁,utmp,p,ttmp)
    integrator.g(Xₛ₋₁,uᵢ₋₁,p,tᵢ₋₁)
    @.. u += (σ + τ)*dt*yₛ₋₁ + 1//8*W.dW*Xₛ₋₁
  else
    # stage s-3
    integrator.f(yₛ₋₃,uᵢ₋₁,p,tᵢ₋₁)
    @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    integrator.g(Xₛ₋₃,utmp,p,ttmp)
    mul!(SXₛ₋₃,Xₛ₋₃,W.dW)
    @.. u += C₁*yₛ₋₃ + 1//8*SXₛ₋₃

    #stage s-2
    @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₄*SXₛ₋₃
    ttmp = tᵢ₋₁ + μₛ₋₃
    integrator.f(yₛ₋₂,utmp,p,ttmp)
    for i in 1:length(W.dW)
      WikRange .= 1 .* (1:length(W.dW) .== i) .* W.dW
      # @.. utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3*@view(Xₛ₋₃[:,i])*W.dW[i]
      mul!(SXₛ₋₂,Xₛ₋₃,WikRange)
      @.. utmp = uᵢ₋₂ + C₁*yₛ₋₃ + 2//3*SXₛ₋₂
      ttmp = tᵢ₋₂ + C₁
      integrator.g(Gₛ,utmp,p,ttmp)
      WikRange .= 1 .* (1:length(W.dW) .== i)
      # @.. @view(Xₛ₋₂[:,i]) =  @view(Gₛ[:,i])
      @.. Xₛ₋₂ =  Gₛ*W.dW
    end
    mul!(SXₛ₋₂,Xₛ₋₂,W.dW)
    @.. u += μₛ₋₂*yₛ₋₂ + 3//8*SXₛ₋₂

    #stage s-1
    @.. utmp = uᵢ₋₂ + μₛ₋₂*yₛ₋₂ + C₁*yₛ₋₃ + δ₅*SXₛ₋₃ + δ₄*SXₛ₋₂
    ttmp = tᵢ₋₂ + μₛ₋₂ + C₁
    integrator.f(yₛ₋₁,utmp,p,ttmp)
    for i in 1:length(W.dW)
      WikRange .= 1 .* (1:length(W.dW) .== i) .* W.dW
      # @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ - 1//6*W.dW[i]*@view(Xₛ₋₃[:,i]) - 1//2*W.dW[i]*@view(Xₛ₋₂[:,i]) + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      @.. utmp = uᵢ₋₁ + μₛ₋₃*yₛ₋₃ + δ₁*yₛ₋₂ + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      mul!(SXₛ₋₁,Xₛ₋₃,WikRange)
      @.. utmp += 1//6*SXₛ₋₁
      mul!(SXₛ₋₁,Xₛ₋₂,WikRange)
      @.. utmp += 1//2*SXₛ₋₁
      ttmp = tᵢ₋₁ + μₛ₋₃ + δ₁
      integrator.g(Gₛ,utmp,p,ttmp)
      WikRange .= 1 .* (1:length(W.dW) .== i)
      # @.. @view(Xₛ₋₁[:,i]) = @view(Gₛ[:,i])
      @.. Xₛ₋₁ = Gₛ*WikRange
    end
    mul!(SXₛ₋₁,Xₛ₋₁,W.dW)
    @.. u += (σ - τ)*dt*yₛ₋₁ + 3//8*SXₛ₋₁

    #stage s
    @.. utmp = uᵢ₋₂ + C₁*yₛ₋₃ + μₛ₋₂*yₛ₋₂ + σ*dt*yₛ₋₁ + δ₆*SXₛ₋₃ + δ₇*SXₛ₋₂ + δ₈*SXₛ₋₁
    ttmp = tᵢ₋₂ + C₁ + μₛ₋₂ + σ*dt
    # utmp = integrator.f(utmp,p,ttmp)
    # u += (τ/σ)*dt*utmp
    #
    # SXₛ₋₁ = zero(uprev)
    # for i in 1:length(W.dW)
    #   utmp = uᵢ₋₁ + θₛ₋₃*(uᵢ₋₁ - uᵢ₋₂) + μₛ₋₃*yₛ₋₃ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ - 3//2*W.dW[i]*@view(Xₛ₋₃[:,i]) -
    #           1//2*W.dW[i]*@view(Xₛ₋₂[:,i]) + 2*W.dW[i]*@view(Xₛ₋₁[:,i]) + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
    #   ttmp = tᵢ₋₁ + θₛ₋₃*(tᵢ₋₁ - tᵢ₋₂) + μₛ₋₃ + δ₂ + δ₃
    #   Gₛ = integrator.g(utmp,p,ttmp)
    #   SXₛ₋₁ += W.dW[i]*@view(Gₛ[:,i])
    # end
    # u += 1//8*SXₛ₋₁

    # memory optimisation
    @.. uᵢ₋₁ += μₛ₋₃*yₛ₋₃
    tᵢ₋₁ += μₛ₋₃
    #now we have uᵢ₋₂ and yₛ₋₃ free
    integrator.f(yₛ₋₃,utmp,p,ttmp)
    @.. SXₛ₋₁ = zero(uprev)
    for i in 1:length(W.dW)
      WikRange .= 1 .* (1:length(W.dW) .== i) .* W.dW
      # @.. uᵢ₋₂ = uᵢ₋₁ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ - 3//2*W.dW[i]*@view(Xₛ₋₃[:,i]) - 1//2*W.dW[i]*@view(Xₛ₋₂[:,i]) + 2*W.dW[i]*@view(Xₛ₋₁[:,i]) + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      @.. uᵢ₋₂ = uᵢ₋₁ + δ₂*yₛ₋₂ + δ₃*yₛ₋₁ + 1//4*SXₛ₋₃ + 3//4*SXₛ₋₂
      mul!(utmp,Xₛ₋₃,WikRange)
      @.. uᵢ₋₂ -= 3//2*utmp
      mul!(utmp,Xₛ₋₂,WikRange)
      @.. uᵢ₋₂ -= 1//2*utmp
      mul!(utmp,Xₛ₋₁,WikRange)
      @.. uᵢ₋₂ += 2*utmp

      tᵢ₋₂ = tᵢ₋₁ + δ₂ + δ₃
      integrator.g(Gₛ,uᵢ₋₂,p,tᵢ₋₂)
      # @.. SXₛ₋₁ += W.dW[i]*@view(Gₛ[:,i])
      mul!(utmp,Gₛ,WikRange)
      @.. SXₛ₋₁ += utmp
    end
    @.. u += (σ + τ)*dt*yₛ₋₃ + 1//8*SXₛ₋₁

    for i in 1:length(W.dW)
      @.. SXₛ₋₁ = zero(uprev)
      for j in 1:length(W.dW)
        WikRange .= 1 .* (1:length(W.dW) .== j)
        mul!(utmp,Xₛ₋₃,WikRange)
        if j != i
          # (i > j) && (@.. SXₛ₋₁ += -vec_χ[i]*W.dW[j]*@view(Xₛ₋₃[:,j]))
          # (i < j) && (@.. SXₛ₋₁ += vec_χ[j]*W.dW[i]*@view(Xₛ₋₃[:,j]))
          (i > j) && (@.. SXₛ₋₁ += -vec_χ[i]*W.dW[j]*utmp)
          (i < j) && (@.. SXₛ₋₁ += vec_χ[j]*W.dW[i]*utmp)
        end
      end
      ttmp = tᵢ₋₁
      @.. utmp = uᵢ₋₁ - 1//4*SXₛ₋₁
      integrator.g(Gₛ,utmp,p,ttmp)
      @.. utmp = uᵢ₋₁ + 1//4*SXₛ₋₁
      integrator.g(Xₛ₋₁,utmp,p,ttmp)
      sqrt_dt *= (length(W.dW) - 1)
      # @.. u += sqrt_dt*(@view(Xₛ₋₁[:,i]) - @view(Gₛ[:,i]))
      WikRange .= 1 .* (1:length(W.dW) .== i)
      mul!(utmp,Xₛ₋₁,WikRange)
      @.. u += sqrt_dt*utmp
      mul!(utmp,Gₛ,WikRange)
      @.. u -= sqrt_dt*utmp
    end
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCKC2ConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack recf, mσ, mτ = cache

  maxeig!(integrator, cache)
  cache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  cache.mdeg = max(3,min(cache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = cache.mdeg
  start     = cache.start
  deg_index = cache.deg_index
  σ = mσ[deg_index]
  τ = mτ[deg_index]

  sqrt_dt   = sqrt(dt)

  μ = recf[start]  # here κ = 0
  tᵢ = t + dt*μ
  tᵢ₋₁ = tᵢ
  tᵢ₋₂ = t

  # stage 1
  uᵢ₋₂ = uprev
  uᵢ = integrator.f(uprev,p,t)
  uᵢ₋₁ = uprev + dt*μ*uᵢ

  # stages 2 upto s-2
  for i in 2:mdeg
    μ, κ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    ν    = 1.0 + κ
    uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
    uᵢ   = dt*μ*uᵢ + ν*uᵢ₋₁ - κ*uᵢ₋₂
    uᵢ₋₂ = uᵢ₋₁
    uᵢ₋₁ = uᵢ
    tᵢ   = dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end

  #2 stage-finishing procedure
  #stage s-1
  uᵢ   = integrator.f(uᵢ₋₁,p,tᵢ₋₁)
  uᵢ₋₂ = uᵢ₋₁ + dt*σ*uᵢ
  u = uᵢ₋₁ + dt*(σ - τ)*uᵢ
  tᵢ₋₂ = tᵢ₋₁ + dt*σ

  #stage s
  uᵢ₋₂ = integrator.f(uᵢ₋₂,p,tᵢ₋₂)
  u += dt*(σ + τ)*uᵢ₋₂

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    Gₛ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
    u += Gₛ .* W.dW

    uᵢ₋₂ = uᵢ₋₁ .+ ((W.dW .^ 2 .- dt) ./ (2 .* sqrt_dt)) .* Gₛ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
    u += 1//2 .* sqrt_dt .* Gₛ₁

    uᵢ₋₂ = uᵢ₋₁ .- ((W.dW .^ 2 .- dt) ./ (2 .* sqrt_dt)) .* Gₛ
    Gₛ₁ = integrator.g(uᵢ₋₂,p,tᵢ₋₁)
    u -= 1//2 .* sqrt_dt .* Gₛ₁

  else
    Gₛ = integrator.g(uᵢ₋₁,p,tᵢ₋₁)
    uᵢ₋₂ = Gₛ*W.dW
    u += uᵢ₋₂

    for i in 1:length(W.dW)
      WikRange = 1 .* (1:length(W.dW) .== i)
      uᵢ = Gₛ*WikRange
      tmp = uᵢ₋₁ - 1//2*sqrt_dt*uᵢ + 1//2*(W.dW[i]/sqrt_dt)*uᵢ₋₂
      Gₛ₁ = integrator.g(tmp,p,tᵢ₋₁)
      u += 1//2*sqrt_dt*(Gₛ₁*WikRange)
      tmp = uᵢ₋₁ + 1//2*sqrt_dt*uᵢ - 1//2*(W.dW[i]/sqrt_dt)*uᵢ₋₂
      Gₛ₁ = integrator.g(tmp,p,tᵢ₋₁)
      u -= 1//2*sqrt_dt*(Gₛ₁*WikRange)
    end
  end
  integrator.u = u
end

@muladd function perform_step!(integrator,cache::SROCKC2Cache,f=integrator.f)
  @unpack uᵢ, tmp, uᵢ₋₁, uᵢ₋₂, k, Gₛ, Gₛ₁, WikRange = cache
  @unpack t,dt,uprev,u,W,p = integrator

  @unpack recf, mσ, mτ = cache.constantcache
  ccache = cache.constantcache

  maxeig!(integrator, cache)
  ccache.mdeg = Int(floor(sqrt((2*dt*integrator.opts.internalnorm(integrator.eigen_est,t)+1.5)/0.811)+1))
  ccache.mdeg = max(3,min(ccache.mdeg,200))-2
  choose_deg!(integrator,cache)

  mdeg      = ccache.mdeg
  start     = ccache.start
  deg_index = ccache.deg_index
  σ = mσ[deg_index]
  τ = mτ[deg_index]

  sqrt_dt   = sqrt(dt)
  μ = recf[start]  # here κ = 0
  tᵢ = t + dt*μ
  tᵢ₋₁ = tᵢ
  tᵢ₋₂ = t

  # stage 1
  @.. uᵢ₋₂ = uprev
  integrator.f(k,uprev,p,t)
  @.. uᵢ₋₁ = uprev + dt*μ*k

  # stages 2 upto s-2
  for i in 2:mdeg
    μ, κ = recf[start + 2*(i-2) + 1], recf[start + 2*(i-2) + 2]
    ν    = 1.0 + κ
    integrator.f(k,uᵢ₋₁,p,t)
    @.. uᵢ   = dt*μ*k + ν*uᵢ₋₁ - κ*uᵢ₋₂
    @.. uᵢ₋₂ = uᵢ₋₁
    @.. uᵢ₋₁ = uᵢ
    tᵢ   = dt*μ + ν*tᵢ₋₁ - κ*tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
    tᵢ₋₁ = tᵢ
  end

  #2 stage-finishing procedure
  #stage s-1
  integrator.f(k,uᵢ₋₁,p,tᵢ₋₁)
  @.. uᵢ₋₂ = uᵢ₋₁ + dt*σ*k
  @.. u = uᵢ₋₁ + dt*(σ - τ)*k
  tᵢ₋₂ = tᵢ₋₁ + dt*σ

  #stage s
  integrator.f(k,uᵢ₋₂,p,tᵢ₋₂)
  @.. u += dt*(σ + τ)*k

  if (typeof(W.dW) <: Number) || (length(W.dW) == 1) || is_diagonal_noise(integrator.sol.prob)
    integrator.g(Gₛ,uᵢ₋₁,p,tᵢ₋₁)
    @.. u += Gₛ*W.dW

    @.. uᵢ₋₂ = uᵢ₋₁ + ((W.dW^2 - dt)/(2*sqrt_dt))*Gₛ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ₋₁)
    @.. u += 1//2*sqrt_dt*Gₛ₁

    @.. uᵢ₋₂ = uᵢ₋₁ - ((W.dW^2 - dt)/(2*sqrt_dt))*Gₛ
    integrator.g(Gₛ₁,uᵢ₋₂,p,tᵢ₋₁)
    @.. u -= 1//2*sqrt_dt*Gₛ₁
  else
    integrator.g(Gₛ,uᵢ₋₁,p,tᵢ₋₁)
    mul!(uᵢ₋₂,Gₛ,W.dW)
    @.. u += uᵢ₋₂

    for i in 1:length(W.dW)
      WikRange .= 1 .* (1:length(W.dW) .== i)
      mul!(uᵢ,Gₛ,WikRange)

      @.. tmp = uᵢ₋₁ - 1//2*sqrt_dt*uᵢ + 1//2*(W.dW[i]/sqrt_dt)*uᵢ₋₂
      integrator.g(Gₛ₁,tmp,p,tᵢ₋₁)
      mul!(tmp,Gₛ₁,WikRange)
      @.. u += 1//2*sqrt_dt*tmp

      @.. tmp = uᵢ₋₁ + 1//2*sqrt_dt*uᵢ - 1//2*(W.dW[i]/sqrt_dt)*uᵢ₋₂
      integrator.g(Gₛ₁,tmp,p,tᵢ₋₁)
      mul!(tmp,Gₛ₁,WikRange)
      @.. u -= 1//2*sqrt_dt*tmp
    end
  end

  integrator.u = u
end


text:
using Test
using SentinelMissings

x  = [0.1,0.2,0.3,-9999.0]
xs = as_sentinel(x,-9999.0)
@test eltype(xs)<:SentinelMissings.SentinelMissing
@test isequal.(copy(xs),[0.1,0.2,0.3,missing]) |> all
@test typeof(copy(xs)) == Vector{Union{Missing,Float64}}
@test isequal.(xs[1:3], [0.1,0.2,0.3]) |> all
@test isapprox.((xs .+ 0.1)[1:3], [0.2,0.3,0.4]) |> all
xs[3] = missing
@test x[3] == -9999.0

ext: jl


text:
export minL2NormSPG

function minL2NormSPG(A::Tx, b, rank1::Int, sigma, iter)where {ETx<:Number, Tx<:AbstractMatrix{ETx}}
xLS_jl = zeros(Complex{Float64}, rank1)
if length(b) != 0
    xinit   = A'*b;
    tau     = 0.0;
    #sigma   = sigmafact*norm(b[:],2);
    # Choose options for GenSPGL
        opts = spgOptions(  optTol = Float32(1e-5),
                        bpTol = Float32(1e-5),
                        decTol = Float32(1e-4),
                        project = l2Norm_project,
                        primal_norm = l2Norm_primaldual,
                        dual_norm = l2Norm_primaldual,
                        proxy = true,
                        ignorePErr = true,
                        iterations = iter,
                        funCompositeR = funCompR2,
                        verbosity = 1)
    # Create Params Dict
        afunT(x) = reshape(x,rank1,1);
        params = Dict("nr"=> rank1,
                    "Ind"=> vec(b) .== 0,
                    "numr"=> rank1,
                    "numc"=> 1,
                    "mode"=> 1,
                    "ls"=> 1,
                    "logical"=> 0,
                    "funPenalty"=> funLS)
    time1 = @elapsed xLS_jl, r, g, info = spgl1(A, vec(b), x = vec(xinit), tau = tau, sigma = sigma,  options = opts, params = params);
    println("LTime = ",string(time1))
    time1 = @elapsed xLS_jl, r, g, info = spgl1(A, vec(b), x = vec(xinit), tau = tau, sigma = sigma,  options = opts, params = params);
    println("LTime2 = ",string(time1))
end
xLS_jl = convert(Array{Float32,1},real(xLS_jl)) + (convert(Array{Float32,1},imag(xLS_jl)))im;
return xLS_jl
end


text:
using SpineInterface
using JuMP
using Clp

db_url = ARGS[1]

using_spinedb(db_url)

optimizer = optimizer_with_attributes(Clp.Optimizer)
m = Model(optimizer)

# Get time indices
time = collect(demand(node=first(node())).indexes)
node__time = [(node=n, time=t) for n in node() for t in time]
node__unit__time = [(node=n, unit=u, time=t) for (n, u) in node__unit() for t in time]
node__node__time = [(node1=n1, node2=n2, time=t) for (n1, n2) in node__node() for t in time]

@variable(m, flow[[(node=n, unit=u, time=t) for (n, u, t) in node__unit__time]] >= 0)
#@variable(m, transfer[[(node1=n1, node2=n2, time=t) for (n1, n2, t) in node__node__time]])

@constraint(m, unit_capacity[(n, u, t) in node__unit__time],
            flow[(node=n, unit=u, time=t)] <= capacity(unit=u))

@constraint(m, nodal_balance[(n, t) in node__time],
            sum(flow[(node=n, unit=u, time=t)]
                for u in node__unit(node=n))
            == demand(node=n, time=t))

@objective(m, Min,
           sum(op_cost(unit=u) * flow[(node=n, unit=u, time=t)]
               for (n, u, t) in node__unit__time))

optimize!(m)

# Write results
output_db_url = ARGS[2]

parameters_1D = Dict(
    "objective" =>
    Dict(
        (model=:MyJuliaModel,) => objective_value(m)
    )
)
write_parameters(parameters_1D, output_db_url)

parameters_3D = Dict(
    "flow" =>
    Dict(
        (node=n, unit=u) => SpineInterface.Map(time, [value(flow[(node=n, unit=u, time=t)]) for t in time])
        for (n, u) in node__unit()
    )
)
write_parameters(parameters_3D, output_db_url)

parameters_2D = Dict(
    "price" =>
    Dict(
        (node=n,) => SpineInterface.Map(time, [dual(nodal_balance[(node=n, time=t)]) for t in time])
        for n in node()
    )
)
write_parameters(parameters_2D, output_db_url)


text:
import Plots.plot

function plot(env::Envelope, variable::Symbol = env.output_for_envelope_; points::Bool = true, legend::Symbol = :best)
    func_names = env.dd_function_[!,:name]
    intervals = size(env.dd_interval_)[1]
    plt = plot(; xlabel = string(env.x_name_), ylabel = string(variable))
    for i in 1:intervals
        int_func = env.dd_interval_[i, :func]
        int_ends = (env.dd_interval_[i, :interval_start] , env.dd_interval_[i, :interval_end])
        for j in 1:length(func_names)
            func = func_names[j]
            style = func == int_func ? :solid : :dash
            plot_options = (label = string(func), color = j, linestyle = style)
            plt = plot(env.splines_[variable][func], int_ends;derivs = false, plot_options = plot_options, plt = plt)
        end
    end
    if points
        plt = plot!(plt, env.dd_evals_[!,env.x_name_], env.dd_evals_[!,variable], seriestype=:scatter; color = :black, label = "points", legend=:none)
    end
    plt = plot!(plt; legend=legend)
    return plt
end


text:
"""
    activated()

Returns true if the Pretend framework is enabled. The default setting is OFF so
that the mocking code is optimized aways by the compiler.
"""
activated() = false

"""
    activate()

Activate the Pretend framework. Normally, it is called at the beginning of a test script.
"""
activate() = @eval(activated() = true)

"""
    deactivate()

Deactivate the Pretend framework. It is not commonly used because it is already the
default setting.
"""
deactivate() = @eval(activated() = false)


text:
using Test
using Bijections

b = Bijection(3,"Hello")
@test b[3] == "Hello"
@test b("Hello") == 3
@test haskey(b, 3)
@test !haskey(b, 4)
@test !haskey(b, "Hello")

b[2] = "Bye"
@test b[2] == "Bye"
@test b("Bye") == 2
@test inverse(b,"Bye") == 2

@test domain(b) == Set([2,3])
@test image(b) == Set(["Bye", "Hello"])

@test length(b) == 2
@test !isempty(b)

bb = inv(b)
@test bb["Bye"] == 2
bb["xx"] = 4
@test length(bb) != length(b)

bb = active_inv(b)
@test bb["Bye"] == 2
b[0] = "Ciao"
@test bb["Ciao"] == 0
@test inv(bb) == b

# iteration test
dom_list = [x for (x,y) in b]
@test Set(dom_list) == domain(b)

# conversion to a Dict
d = Dict(b)
@test d[0] == b[0]


text:
## functions to set/read the different parameters for the extraction
##
## return the default parameters
##
function set_default_params()::meta
    def = meta()

    ## Directories
    def.rootdir =  ENV["GAIA_ROOT"]
    # def.wdir    = @sprintf("%s/e2e_products", def.rootdir)
    # def.votdir  = @sprintf("%s/votable", def.wdir)
    # def.plotdir = @sprintf("%s/plots-select" , def.wdir)
    # def.ocdir   = @sprintf("%s/oc", def.wdir)
    def.wdir= "."
    def.votdir= "."
    def.plotdir= "."
    def.ocdir= "."
    def.votname = "test.vot"
    def.prefile= "ocres"

    ## placeholder for the uuid
    def.uuid= UUID("cfc395e8-590f-11e8-1f13-43a2532b2fa8")

    ## for standalone
    ##
    def.optim           = "no"             ## yes|no for weighting/dbscan optimization
    # if optim no
    def.w3d             = 7.0               ## w3d weighting
    def.wvel            = 8.0               ## wvel weighting
    def.whrd            = 2.5               ## whrd weighting
    def.eps             = 2.1               ## epsilon DBSCAN
    def.mcl             = 18                ## min_cluster DBSCAN
    def.mnei            = 7                 ## min_neighbor DBSCAN
    def.maxdist         = 2100.0            ## maximum distance to filter stars in pc
    def.mindist         = 0.0               ## minimum distance to filter stars in pc

    ## MCMC
    ##
    def.minQc    = 2.7
    def.minQn    = 40
    def.maxQn    = 5000
    def.forcedminstars = 30
    def.mingoodsolution = 10        # minimum good solution for Q in check_qminqstar
    def.niterqminq = 500            # check_qminqstar..

    def.nburnout  = 500             # burn-in iterations
    def.nchain    = 5000            # MCMC iterations==number of CHAINS
    def.maxiter   = 50000           # maximum iteration

    ### Cycle parameters
    ###
    def.cyclemax            = 3         ## maximum of cycles
    def.minstarselection    = 50        ## minimum of stars to select solution in a cycle...
    def.minstarstop         = 50        ## condition to stop cycling
    def.minchainreached     = 100       ## minimum chain to analyze solution
    def.qcminstop           = 1.5       ## more condition on Qc to stop cycling after the first pass
    def.wratiominstop       = 0.2       ## minimum ratio btwn w3d and wvel (otherwise not an OC)

    ## prior settings
    ##
    def.epsmean   = 2.5             # epsilon mean
    def.epsdisp   = 1.5             # epsilon dispersion
    def.min_nei   = 10              # min_nei mean
    def.min_cl    = 15              # min_cl mean
    def.ncoredisp = 10              # dispersion for min_nei and min_cl
    def.w3dmean   = 6.0
    def.w3ddisp   = 4.0
    def.wvelmean  = 6.0
    def.wveldisp  = 4.0
    def.whrdmean  = 2.0
    def.whrddisp  = 1.5

    ## random walk parameters
    ##
    def.eps_rw        = 0.4             # eps dispersion
    def.minnei_rw     = 4               # min_nei dispersion
    def.mincl_rw      = 4               # min_cl dispersion
    def.w3d_rw        = 0.8             # w3d dispersion
    def.wvel_rw       = 0.8             # wvel dispersion
    def.whrd_rw       = 0.8             # whrd dispersion

    ## Metrics & parameters find_clusters
    def.clustermax    = 1000
    def.aperture2d    = 1.5
    def.maxaperture2d = 15
    def.aperture3d    = 3
    def.maxaperture3d = 20
    def.aperturev     = 3.0
    def.maxaperturev  = 20
    def.nboot         = 50          # bootstrap number

    def.labels        = "Qc"        # method to get labels: Qc|Qn

    return(def)
end

## extracting meta from .ext file
function read_params(file, verbose=true)
    par= set_default_params()
    println("## All parameters set to default...")
    t= readdlm(file,comments=true, comment_char='#')

    s= size(t)
    for i in 1:s[1]
        if verbose println(t[i,:]) end
        if t[i,2] == "="
            set_param!(par,t[i,1],t[i,3])
        else
            p= t[i,1]
            println("### parsing parameter file could be wrong near $p ...")
        end
    end

    println("## Parameters read from $file")

    return(par)
end

## for parsing .ext file
function set_param!(def, parstr,value)
    if parstr == "rootdir" def.rootdir= value end
    if parstr == "wdir"    def.wdir= value end
    if parstr == "votdir"  def.votdir= value end
    if parstr == "plotdir" def.plotdir= value end
    if parstr == "ocdir"   def.ocdir= value end
    if parstr == "votname" def.votname= value end
    if parstr == "prefile" def.prefile= value end
    if parstr == "uuid" def.uuid= value end

    if parstr == "optim" def.optim= value end
    if parstr == "w3d" def.w3d= value end
    if parstr == "wvel" def.wvel= value end
    if parstr == "whrd" def.whrd= value end
    if parstr == "eps" def.eps=  value end
    if parstr == "mcl" def.mcl= value end
    if parstr == "mnei" def.mnei= value end
    if parstr == "maxdist" def.maxdist= value end
    if parstr == "mindist" def.mindist= value end

    if parstr == "minQc"           def.minQc= value end
    if parstr == "minQn"           def.minQn= value end
    if parstr == "maxQn"           def.maxQn= value end
    if parstr == "forcedminstars"  def.forcedminstars= value end
    if parstr == "mingoodsolution" def.mingoodsolution= value end
    if parstr == "niterqminq"      def.niterqminq= value end

    if parstr == "nburnout" def.nburnout= value end
    if parstr == "nchain"   def.nchain= value end
    if parstr == "maxiter"  def.maxiter= value end

    if parstr == "cyclemax" def.cyclemax= value  end
    if parstr == "minstarselection" def.minstarselection= value  end
    if parstr == "minstarstop" def.minstarstop= value  end
    if parstr == "minchainreached" def.minchainreached= value  end
    if parstr == "qcminstop" def.qcminstop= value  end
    if parstr == "wratiominstop" def.wratiominstop= value  end

    if parstr == "epsmean"   def.epsmean= value end
    if parstr == "epsdisp"   def.epsdisp= value end
    if parstr == "min_nei"   def.min_cl= value end
    if parstr == "min_cl"    def.min_cl= value end
    if parstr == "ncoredisp" def.ncoredisp= value end
    if parstr == "w3dmean"   def.w3dmean= value end
    if parstr == "w3ddisp"   def.w3ddisp= value end
    if parstr == "wvelmean"  def.wvelmean= value end
    if parstr == "wveldisp"  def.wveldisp= value end
    if parstr == "whrdmean"  def.whrdmean= value end
    if parstr == "whrddisp"  def.whrddisp= value end

    if parstr == "eps_rw"  def.eps_rw = value end
    if parstr == "minnei_rw"  def.minnei_rw= value end
    if parstr == "mincl_rw"  def.mincl_rw= value end
    if parstr == "w3d_rw"  def.w3d_rw= value end
    if parstr == "wvel_rw"  def.wvel_rw= value end
    if parstr == "whrd_rw"  def.whrd_rw= value end

    ## Metrics find_clusters
    if parstr == "clustermax"     def.clustermax= value end
    if parstr == "aperture2d"     def.aperture2d= value end
    if parstr == "maxaperture2d"  def.maxaperture2d= value end
    if parstr == "aperture3d"     def.aperture3d= value end
    if parstr == "maxaperture3d"  def.maxaperture3d= value end
    if parstr == "aperturev"      def.aperturev= value end
    if parstr == "maxaperturev"   def.maxaperturev= value end
    if parstr == "nboot"          def.nboot= value end

    if parstr == "labels"  def.labels= value end

end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


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

using Markdown
using InteractiveUtils

# ╔═╡ 92ae949c-2a00-11eb-0f11-6570ad25144b
using Revise, DifferentialEquations

# ╔═╡ 897e68d2-2a3f-11eb-32e4-0b56dec432e8
using NBodyLeapFrog

# ╔═╡ 4aaf7264-2a45-11eb-1a8e-7976e5bbc634
using StaticArrays

# ╔═╡ 4f659e96-2a45-11eb-03ab-01dbf5e82eb4
tosvec(x::V3) = SVector(x.x, x.y, x.z)

# ╔═╡ 7e0214da-2a45-11eb-26a5-9bb0c6b1518a
NBodyLeapFrog.acceleration(x::SVector, y::SVector, m, G) = NBodyLeapFrog.acceleration(V3(x...), V3(y...), m, G)

# ╔═╡ e0db6922-2a0d-11eb-040f-8ba9cd6d6f5b
function solar_problem(r0, v0; dt, n, params=nothing)
	function dp(v::AbstractVector{T}, x, params, t) where T
		@show x
		res = zeros(T, length(x))
		for i=1:length(res)
			for j=1:length(res)
				if i!= j
					res[i] += NBodyLeapFrog.acceleration(x[i], x[j], params.m[j], params.G)
				end
			end
		end
		return res
	end
	function dq(v, x, params, t)
		@show v
		v
	end
	ODEProblem(DynamicalODEFunction{false}(dp,dq), ArrayPartition(v0,r0), (0.0, n*dt), params)
end

# ╔═╡ e7a72072-2a49-11eb-3bec-3766113cf185
ArrayPartition

# ╔═╡ bf0d3bc4-2a3f-11eb-2c09-c7e1e709235d
v, r, m = let
	n = 100
    planets = Bodies.chunit_day2year.(Bodies.set)
    getfield.(planets, :v),
    getfield.(planets, :r),
    getfield.(planets, :m)
end

# ╔═╡ a4179630-2a48-11eb-0738-53f92fca7ef6
prob = solar_problem(r, v, dt=0.01, n=4000, params=(m=m, G=NBodyLeapFrog.G_year_AU))

# ╔═╡ d6dc5db2-2a48-11eb-2e5a-2f3574592e7e
DiffEqBase.recursive_unitless_bottom_eltype(::Type{V3{T}}) where T = T#V3{T}

# ╔═╡ 0a5d0150-2a49-11eb-3a04-cbe528382566
function DiffEqBase.recursive_unitless_eltype(::Type{V3{T}}) where T
	V3{T}
end

# ╔═╡ b47297c8-2a48-11eb-0f80-b352b8100a64
DiffEqBase.recursive_unitless_bottom_eltype(r |> typeof)

# ╔═╡ 6b99eb4a-2a49-11eb-2a9a-b5cf52a26ded
begin
	Base.one(::Type{V3{T}}) where T = one(T)
	Base.one(::V3{T}) where T = one(T)
	Base.oneunit(x::V3) = one(x)
	Base.:(/)(x::V3, y::Real) = V3(x.x/y, x.y/y, x.z/y)
	Base.broadcastable(x::V3) = Ref(x)
	Base.isnan(x::V3) = isnan(x.x) || isnan(x.y) || isnan(x.z)
end

# ╔═╡ 2a4f0bce-2a45-11eb-19e8-41512b3e66be
 #solve(prob)

# ╔═╡ 462eb43c-2a01-11eb-1726-593e84125e48
alg = McAte2()

# ╔═╡ 2f7fb8e6-2a01-11eb-229c-210763e47bca
integrator = DiffEqBase.__init(prob, alg; dt=0.01, tstops=40.0)

# ╔═╡ 32f10746-2a01-11eb-01ec-6b90be753181
 solve!(integrator)

# ╔═╡ 372f56c6-2a01-11eb-0826-432eefb6198e
integrator.sol

# ╔═╡ Cell order:
# ╠═92ae949c-2a00-11eb-0f11-6570ad25144b
# ╠═897e68d2-2a3f-11eb-32e4-0b56dec432e8
# ╠═4aaf7264-2a45-11eb-1a8e-7976e5bbc634
# ╠═4f659e96-2a45-11eb-03ab-01dbf5e82eb4
# ╠═7e0214da-2a45-11eb-26a5-9bb0c6b1518a
# ╠═e0db6922-2a0d-11eb-040f-8ba9cd6d6f5b
# ╠═e7a72072-2a49-11eb-3bec-3766113cf185
# ╠═bf0d3bc4-2a3f-11eb-2c09-c7e1e709235d
# ╠═a4179630-2a48-11eb-0738-53f92fca7ef6
# ╠═d6dc5db2-2a48-11eb-2e5a-2f3574592e7e
# ╠═0a5d0150-2a49-11eb-3a04-cbe528382566
# ╠═b47297c8-2a48-11eb-0f80-b352b8100a64
# ╠═6b99eb4a-2a49-11eb-2a9a-b5cf52a26ded
# ╠═2a4f0bce-2a45-11eb-19e8-41512b3e66be
# ╠═462eb43c-2a01-11eb-1726-593e84125e48
# ╠═2f7fb8e6-2a01-11eb-229c-210763e47bca
# ╠═32f10746-2a01-11eb-01ec-6b90be753181
# ╠═372f56c6-2a01-11eb-0826-432eefb6198e


text:
###############################################################################
#subtraction
@universal function sub(a::Ubound, b::Unum)
  lb = resolve_lower(a.lower - b)
  hb = resolve_upper(a.upper - b)

  (is_ulp(lb) && is_ulp(hb)) ? resolve_as_utype!(lb, hb) : B(lb, hb)
end

@universal function sub(a::Unum, b::Ubound)
  lb = resolve_lower(a - b.upper)
  hb = resolve_upper(a - b.lower)

  (is_ulp(lb) && is_ulp(hb)) ? resolve_as_utype!(lb, hb) : B(lb, hb)
end

@universal function sub(a::Ubound, b::Ubound)
  lb = resolve_lower(a.lower - b.upper)
  hb = resolve_upper(a.upper - b.lower)

  is_sss(lb) && is_sss(hb) && (@signof(lb) == @signof(hb)) && return sss(U, @signof(lb))
  is_mmr(lb) && is_mmr(hb) && (@signof(lb) == @signof(hb)) && return mmr(U, @signof(lb))

  (is_ulp(lb) && is_ulp(hb)) ? resolve_as_utype!(lb, hb) : B(lb, hb)
end

@universal function additiveinverse!(a::Ubound)
  hb = additiveinverse!(a.lower)
  a.lower = additiveinverse!(a.upper)
  a.upper = hb
  return a
end

#unary subtraction creates a new unum and flips it.
@universal function -(x::Ubound)
  additiveinverse!(copy(x))
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Pipe
using LinearAlgebra

function calc_densities(matrix::Array{Float64,2})

    resolution = (200, 200)

    densities = zeros(Float64, resolution...)

    cartesians = performBinning(matrix, resolution)

    map(cidx -> densities[cidx] += 1, cartesians)

    return (cartesians, densities)
end

function logtransform(matrix)
    @pipe matrix |> _ .+ 1 |> log10.(_)
end

# utility broadcasting function
function replaceZeroIndexes(idxs, newval=0)
    idxs[idxs .== 0] .= 1

    idxs
end

# make binning
function performBinning(matrix, resolution, columns=[1,2])
    mins = minimum(matrix[:,columns], dims=1)
    maxs = maximum(matrix[:,columns], dims=1)

    interval = maxs - mins

    cartesians = matrix |>
                (m -> m[:,columns]) |> # get the first two columns
                (tc -> tc .- mins) |> # subtract min values for each column
                (tc -> tc ./ interval) |> # divide with interval to get relative range from 0 to 1
                (tc ->  tc .* hcat(resolution...))  |> # multiply with resolution to get bin number
                (tc -> trunc.(Int64, tc)) |> # floor to get bins as integers
                (tc -> replaceZeroIndexes(tc, 1)) |> # handle edge case for bins on position zero
                tc -> [CartesianIndex(tc[i,:]...) for i in 1:size(tc,1)] # get list of indexes to facilitate broadcasting
end

# run smoothing for one dimension
function smooth1D(Y, lambda)

    m, n = size(Y)

    E = Matrix{Int}(I, m, m)
    D1 = diff(E, dims=1)
    D2 = diff(D1, dims=1)
    P = lambda.^2 .* D2' * D2 + lambda .* 2 .* D1' * D1;

    Z = (E + P) \ Y

end

function getDensity(densities, cartesians)
    [densities[gd] for gd in cartesians]
end

"""
Create a density scatter2D plot with smoothing
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.105.459&rep=rep1&type=pdf

...
# Arguments
- `x:: Array{Float64,1}` : x axis values
- `y:: Array{Float64,1}` : y axis values
- `lambda:: Int64=20` : input for the smoothing function
...
"""
function smoothe(x::Array{Float64,1}, y::Array{Float64,1}, lambda=20)::Tuple{Array{Float64,1},Array{Int64,1}}

    cartesians, densities = calc_densities(hcat(x, y))

    nbins = size(densities)

    zvalues = @pipe densities |>
        logtransform(_) |>
        smooth1D(_, nbins[2] / lambda) |>
        _' |>
        smooth1D(_, nbins[1] / lambda) |>
        _' |>
        getDensity(_, cartesians)

    uqindexes = uniqueIndexes(cartesians)

    return (zvalues,uqindexes)

end

function uniqueIndexes(cartesians::Array{CartesianIndex{2},1})::Array{Int64,1}
    ctindexed  = [(i,ct) for (i,ct) in enumerate(cartesians)]
    uqindexes = unique(x->x[2],ctindexed) |> uqCtindexes -> getfield.(uqCtindexes,1)

    uqindexes
end


text:
function TDMParser()
    return TDMParser(())
end

function characters(obj::DefaultHandler, arg0::Vector{jchar}, arg1::jint, arg2::jint)
    return jcall(obj, "characters", void, (Vector{jchar}, jint, jint), arg0, arg1, arg2)
end

function declaration(obj::ContentHandler, arg0::JString, arg1::JString, arg2::JString)
    return jcall(obj, "declaration", void, (JString, JString, JString), arg0, arg1, arg2)
end

function end_document(obj::DefaultHandler)
    return jcall(obj, "endDocument", void, ())
end

function end_element(obj::DefaultHandler, arg0::JString, arg1::JString, arg2::JString)
    return jcall(obj, "endElement", void, (JString, JString, JString), arg0, arg1, arg2)
end

function end_prefix_mapping(obj::DefaultHandler, arg0::JString)
    return jcall(obj, "endPrefixMapping", void, (JString,), arg0)
end

function error(obj::DefaultHandler, arg0::SAXParseException)
    return jcall(obj, "error", void, (SAXParseException,), arg0)
end

function fatal_error(obj::DefaultHandler, arg0::SAXParseException)
    return jcall(obj, "fatalError", void, (SAXParseException,), arg0)
end

function get_conventions(obj::TDMParser)
    return jcall(obj, "getConventions", IERSConventions, ())
end

function get_data_context(obj::TDMParser)
    return jcall(obj, "getDataContext", DataContext, ())
end

function get_file_format(obj::TDMParser)
    return jcall(obj, "getFileFormat", TDMParser_TDMFileFormat, ())
end

function get_mission_reference_date(obj::TDMParser)
    return jcall(obj, "getMissionReferenceDate", AbsoluteDate, ())
end

function ignorable_whitespace(obj::DefaultHandler, arg0::Vector{jchar}, arg1::jint, arg2::jint)
    return jcall(obj, "ignorableWhitespace", void, (Vector{jchar}, jint, jint), arg0, arg1, arg2)
end

function is_simple_eop(obj::TDMParser)
    return jcall(obj, "isSimpleEOP", jboolean, ())
end

function notation_decl(obj::DefaultHandler, arg0::JString, arg1::JString, arg2::JString)
    return jcall(obj, "notationDecl", void, (JString, JString, JString), arg0, arg1, arg2)
end

function parse(obj::TDMParser, arg0::InputStream)
    return jcall(obj, "parse", TDMFile, (InputStream,), arg0)
end

function parse(obj::TDMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", TDMFile, (InputStream, JString), arg0, arg1)
end

function parse(obj::TDMParser, arg0::JString)
    return jcall(obj, "parse", TDMFile, (JString,), arg0)
end

function parse_key_value(obj::TDMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parseKeyValue", TDMFile, (InputStream, JString), arg0, arg1)
end

function parse_xml(obj::TDMParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parseXml", TDMFile, (InputStream, JString), arg0, arg1)
end

function processing_instruction(obj::DefaultHandler, arg0::JString, arg1::JString)
    return jcall(obj, "processingInstruction", void, (JString, JString), arg0, arg1)
end

function resolve_entity(obj::DefaultHandler, arg0::JString, arg1::JString)
    return jcall(obj, "resolveEntity", InputSource, (JString, JString), arg0, arg1)
end

function set_document_locator(obj::DefaultHandler, arg0::Locator)
    return jcall(obj, "setDocumentLocator", void, (Locator,), arg0)
end

function skipped_entity(obj::DefaultHandler, arg0::JString)
    return jcall(obj, "skippedEntity", void, (JString,), arg0)
end

function start_document(obj::DefaultHandler)
    return jcall(obj, "startDocument", void, ())
end

function start_element(obj::DefaultHandler, arg0::JString, arg1::JString, arg2::JString, arg3::Attributes)
    return jcall(obj, "startElement", void, (JString, JString, JString, Attributes), arg0, arg1, arg2, arg3)
end

function start_prefix_mapping(obj::DefaultHandler, arg0::JString, arg1::JString)
    return jcall(obj, "startPrefixMapping", void, (JString, JString), arg0, arg1)
end

function unparsed_entity_decl(obj::DefaultHandler, arg0::JString, arg1::JString, arg2::JString, arg3::JString)
    return jcall(obj, "unparsedEntityDecl", void, (JString, JString, JString, JString), arg0, arg1, arg2, arg3)
end

function warning(obj::DefaultHandler, arg0::SAXParseException)
    return jcall(obj, "warning", void, (SAXParseException,), arg0)
end

function with_conventions(obj::TDMParser, arg0::IERSConventions)
    return jcall(obj, "withConventions", TDMParser, (IERSConventions,), arg0)
end

function with_data_context(obj::TDMParser, arg0::DataContext)
    return jcall(obj, "withDataContext", TDMParser, (DataContext,), arg0)
end

function with_file_format(obj::TDMParser, arg0::TDMParser_TDMFileFormat)
    return jcall(obj, "withFileFormat", TDMParser, (TDMParser_TDMFileFormat,), arg0)
end

function with_mission_reference_date(obj::TDMParser, arg0::AbsoluteDate)
    return jcall(obj, "withMissionReferenceDate", TDMParser, (AbsoluteDate,), arg0)
end

function with_simple_eop(obj::TDMParser, arg0::jboolean)
    return jcall(obj, "withSimpleEOP", TDMParser, (jboolean,), arg0)
end


text:
using JTetra
using Base.Test
include("derivative-test.jl")

println("==== Test derivatives of mass matrices ====")

# read tetrahedral mesh
mesh = getTetraMesh("points.txt","tetras.txt")

# scalar and tensor cell properties for mass matrices
a = exp(randn(mesh.nc))
A = zeros(mesh.nc,6)
for i = 1:mesh.nc
  # generate symmetric positive definite 3x3 matrix
  Ai  = randn(3,3)
  u,V = eig(Ai + Ai')
  B   = V' * diagm(exp(u)) * V
  A[i,1] = B[1,1]
  A[i,2] = B[2,2]
  A[i,3] = B[3,3]
  A[i,4] = B[1,2]
  A[i,5] = B[1,3]
  A[i,6] = B[2,3]
end

# random edge vector
u = rand(mesh.ne)

# derivative test
f(m)    = getEdgeMassMatrix(mesh, m) * u
df(v,m) = getdEdgeMassMatrix(mesh, m, u) * vec(v)

println("Edge mass matrix: isotropic")
checkDerivativeIsotropic,error,order = checkDerivativeMax(f,df,a)
@test checkDerivativeIsotropic

println("Edge mass matrix: anisotropic")
checkDerivativeAnisotropic,error,order = checkDerivativeMax(f,df,A)
@test checkDerivativeAnisotropic

print("Passed!\n")

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
function assert_nonnegative_diagonal(cov)
    if !all(diag(cov) .>= 0)
        @error "Non-positive variances" cov diag(cov)
        error("The provided covariance has non-positive entries on the diagonal!")
    end
end


text:
module Utils_cons

export eᵢ,
       idx2var,
       idx2var_arr,
       idx2var_dic,
       idx2varxx̄ᵀtyȳᵀ,
       get_Gᴿ

"""The standard basis vector eᵢ in dimension n"""
eᵢ(n::Int,i::Int) = [j==i for j in 1:n]

"""Takes an array of exponents α's and gives array of same shape L(xᵅ)
Input: Array B of multi-indices: α
Output: L of x of array of indices: L(x.(B))
 """

idx2var(var,index) = isempty(index) ?  0 : map(α -> var[α],index)
idx2var(var,index,coef) = sum(coef .* idx2var(var, index))
idx2var(var,index,coef,g) = idx2var(var,g .+ index,coef)

make_gen_aff(arr) = reshape([arr...],size(arr)...)
idx2var_arr(var,index,coef) = make_gen_aff(map((x,y)-> Utils_cons.idx2var(var,x,y),index,coef))
idx2var_arr(var,index,coef,g) = make_gen_aff(map((x,y)-> Utils_cons.idx2var(var,x,y,g),index,coef))

function idx2var_dic(var,index,coef)
    kl = [keys(coef)...]
    f(b) = make_gen_aff(idx2var_arr(var, index[b],coef[b]))
    return Dict(zip(kl,map(b->f(b),kl)))
end
function idx2var_dic(var,index,coef,g)
    kl = [keys(coef)...]
    f(b) = make_gen_aff(idx2var_arr(var, index[b],coef[b],g))
    return Dict(zip(kl,map(b->f(b),kl)))
end

## Complex side
"""Returns Lᴿ(ℜe(xx̄ᵀ⊗yȳᵀ)), Lᴿ(ℑm(xx̄ᵀ⊗yȳᵀ)) """
function idx2varxx̄ᵀtyȳᵀ(Lx,xx̄ᵀtyȳᵀ)
    T =  Dict()
    for r in ["real", "imag"]
        sleut = [keys(xx̄ᵀtyȳᵀ[r])...]
        T[r] =   idx2var(Lx,xx̄ᵀtyȳᵀ[r][popfirst!(sleut)])
        for k in sleut
            (k[1] == '+') ? T[r] += idx2var(Lx,xx̄ᵀtyȳᵀ[r][k]) : T[r] -= idx2var(Lx,xx̄ᵀtyȳᵀ[r][k])
        end
    end
    return T
end

""" ?"""
function get_Gᴿ(xx̄ᵀ_tens_yȳᵀ) #???????
    Rxx̄ᵀ_tens_yȳᵀ = xx̄ᵀ_tens_yȳᵀ["real"]
    Ixx̄ᵀ_tens_yȳᵀ = xx̄ᵀ_tens_yȳᵀ["imag"]
    rkeys = [keys(Rxx̄ᵀ_tens_yȳᵀ)...]
    ikeys = [keys(Ixx̄ᵀ_tens_yȳᵀ)...]

    Gℝ = Dict()
    sgn_mat = Dict()
    for k in 1:8
        Re =    xx̄ᵀ_tens_yȳᵀ["real"][rkeys[k]]
        Im =    xx̄ᵀ_tens_yȳᵀ["imag"][ikeys[k]]
        Gℝ[k] = vcat(hcat(Re,Im), hcat(Im,Re))

        rsgn = rkeys[k][1] == '+' ? 1 : -1
        isgn = ikeys[k][1] == '+' ? 1 : -1
        sgn_mat[k] = [rsgn -isgn; isgn rsgn]
    end
    return Gℝ,sgn_mat
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
module TestAppend

include("preamble.jl")

@testset begin
    @test append!!([0.0], [1.0]) == [0.0, 1.0]
    @test append!!([0], [1.0]) == [0.0, 1.0]
    @test append!!([0], (1.0,)) == [0.0, 1.0]
    @test append!!((0,), (1,)) === (0, 1)
    @test append!!((a=0,), pairs((b=1,))) === (a=0, b=1)
    @test append!!(ImmutableDict(:a=>1), ImmutableDict(:b=>2)) ==
        ImmutableDict(ImmutableDict(:a=>1), :b=>2)
    @test append!!("a", "b") === "ab"
    @test append!!(SVector(0), [1])::Vector == [0, 1]
    @test append!!([0], SVector(1))::Vector == [0, 1]
end

end  # module


text:
using .APIUtils

using CEnum

using Printf

# low-level wrappers
const CUdeviceptr = CuPtr{Cvoid}
const CUarray = CuArrayPtr{Cvoid}
include("libcuda_common.jl")
include("error.jl")
include("libcuda.jl")
include("libcuda_aliases.jl")

# high-level wrappers
include("types.jl")
include("version.jl")
include("devices.jl")
include("context.jl")
include(joinpath("context", "primary.jl"))
include("stream.jl")
include("memory.jl")
include("module.jl")
include("events.jl")
include("execution.jl")
include("profile.jl")
include("occupancy.jl")


text:
# Copied from ModellingToolkit.

"""
    @named x = f(...)

For variable assiments like the above, this is equivalent to:
`x = f(..., name=:x)`.

This shorthand copied from ModellingToolkit is a handy way to
cut down on typing and make sure that names of objects are kept
in sync with variable names.
"""
macro named(expr)
    name, call = split_assign(expr)
    if Meta.isexpr(name, :ref)
        name, idxs = name.args
        check_name(name)
        esc(_named_idxs(name, idxs, :($(gensym()) -> $call)))
    else
        check_name(name)
        esc(:($name = $(_named(name, call))))
    end
end
export @named

macro named(name::Symbol, idxs, call)
    esc(_named_idxs(name, idxs, call))
end

function _named(name, call, runtime=false)
    has_kw = false
    call isa Expr || throw(Meta.ParseError("The rhs must be an Expr. Got $call."))
    if length(call.args) >= 2 && call.args[2] isa Expr
        # canonicalize to use `:parameters`
        if call.args[2].head === :kw
            call.args[2] = Expr(:parameters, Expr(:kw, call.args[2].args...))
            has_kw = true
        elseif call.args[2].head === :parameters
            has_kw = true
        end
    end

    if !has_kw
        param = Expr(:parameters)
        if length(call.args) == 1
            push!(call.args, param)
        else
            insert!(call.args, 2, param)
        end
    end

    kws = call.args[2].args

    if !any(kw->(kw isa Symbol ? kw : kw.args[1]) == :name, kws) # don't overwrite `name` kwarg
        pushfirst!(kws, Expr(:kw, :name, runtime ? name : Meta.quot(name)))
    end
    call
end

function _named_idxs(name::Symbol, idxs, call)
    if call.head !== :->
        throw(ArgumentError("Not an anonymous function"))
    end
    if !isa(call.args[1], Symbol)
        throw(ArgumentError("not a single-argument anonymous function"))
    end
    sym, ex = call.args
    ex = Base.Cartesian.poplinenum(ex)
    ex = _named(:(Symbol($(Meta.quot(name)), :_, $sym)), ex, true)
    ex = Base.Cartesian.poplinenum(ex)
    :($name = $map($sym->$ex, $idxs))
end

check_name(name) = name isa Symbol || throw(Meta.ParseError("The lhs must be a symbol (a) or a ref (a[1:10]). Got $name."))

function split_assign(expr)
    if !(expr isa Expr && expr.head === :(=) && expr.args[2].head === :call)
        throw(ArgumentError("expression should be of the form `sys = foo(a, b)`"))
    end
    name, call = expr.args
end


text:
layers(::Type{WorldClim{Weather}}) = (:tmin, :tmax, :prec)

"""
    getraster(T::Type{WorldClim{Weather}}, [layer::Union{Symbol,Tuple}]; date) => Vector{String}
    getraster(T::Type{WorldClim{Weather}}, layer::Symbol, date)

Download WorldClim weather data, choosing `layer` from: `$(layers(WorldClim{Weather}))`.

Without a layer argument, all layers will be downloaded, and a tuple of paths is returned.
If the data is already downloaded the path will be returned.
"""
function getraster(T::Type{WorldClim{Weather}}, layer::Symbol; date)
    getraster(T::Type{WorldClim{Weather}}, layer::Symbol, date)
end
function getraster(T::Type{WorldClim{Weather}}, layer::Symbol, date::Tuple)
    getraster(T, layer, _date_sequence(date, Month(1)))
end
function getraster(T::Type{WorldClim{Weather}}, layer::Symbol, date::Dates.TimeType)
    decadestart = Date.(1960:10:2020)
    for i in eachindex(decadestart[1:end-1])
        # At least one date is in the decade
        date >= decadestart[i] && date < decadestart[i+1] || continue
        zip_path = zippath(T, layer; decade=decadestart[i])
        _maybe_download(zipurl(T, layer; decade=decadestart[i]), zip_path)
        zf = ZipFile.Reader(zip_path)
        raster_path = rasterpath(T, layer; date=date)
        mkpath(dirname(raster_path))
        if !isfile(raster_path)
            raster_name = rastername(T, layer, date)
            println("Writing $(raster_path)...")
            write(raster_path, read(_zipfile_to_read(raster_name, zf)))
        end
        push!(raster_paths, raster_path)
        close(zf)
        return raster_path
    end
    error("Date $date not between 1960 and 2020")
end
function getraster(T::Type{WorldClim{Weather}}, layer::Symbol, dates::AbstractArray)
    getraster.(T, layer, dates)
end

const WEATHER_DECADES = Dict(Date(1960) => "1960-1969",
                             Date(1970) => "1970-1979",
                             Date(1980) => "1980-1989",
                             Date(1990) => "1990-1999",
                             Date(2000) => "2000-2009",
                             Date(2010) => "2010-2018")

rastername(T::Type{<:WorldClim{Weather}}, layer; date) =
    joinpath("wc2.1_2.5m_$(layer)_$(_date2string(T, date)).tif")
zipname(T::Type{<:WorldClim{Weather}}, layer; decade) =
    "wc2.1_2.5m_$(layer)_$(WEATHER_DECADES[decade]).zip"
zipurl(T::Type{<:WorldClim{Weather}}, layer; decade) =
    joinpath(WORLDCLIM_URI, "hist", zipname(T, layer; decade))
zippath(T::Type{<:WorldClim{Weather}}, layer; decade) =
    joinpath(rasterpath(T), "zips", zipname(T, layer; decade))

# Utility methods

_dateformat(::Type{<:WorldClim}) = DateFormat("yyyy-mm")

_filename2date(T::Type{<:WorldClim}, fn::AbstractString) =
    _string2date(T, basename(fn)[findfirst(r"\d\d\d\d-\d\d", basename(fn))])


text:
using Revise
using PseudoArcLengthContinuation, SparseArrays, LinearAlgebra, Plots, Setfield
const PALC = PseudoArcLengthContinuation

normInf = x -> norm(x, Inf)

k = 2
F = (x, p) -> (@. p + x - x^(k+1)/(k+1))
Jac_m = (x, p) -> diagm(0 => 1  .- x.^k)

opt_newton0 = PALC.NewtonPar(tol = 1e-11, verbose = true)
	out0, hist, flag = @time PALC.newton(
		x -> F(x, 1.),
		[0.8],
		opt_newton0)

opts = PALC.ContinuationPar(dsmax = 0.1, dsmin = 1e-3, ds = -0.001, maxSteps = 130, pMin = -3., pMax = 3., saveSolEveryNsteps = 0, computeEigenValues = true, detectBifurcation = true, newtonOptions = NewtonPar(verbose = false))

println("\n"*"#"^120)
	br0, u1, _ = @time PALC.continuation(F,Jac_m,[0.8],1.,opts;verbosity=0,printSolution = (x,p) -> x[1]) #130 => 73 points
	plotBranch(br0);title!("")

br0.branch[:,:]'
br0.n_unstable
br0.n_imag
br0.sol
br0.eig


text:
#=
# 1D Poisson, Dirichlet bc
# CG, Linear element
# Simplest test possible
=#

### If the Femshop package has already been added, use this line #########
using Femshop # Note: to add the package, first do: ]add "https://github.com/paralab/femshop.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Femshop)
#     include("../Femshop.jl");
#     using .Femshop
# end
##########################################################################

init_femshop("poisson1d");

# Optionally generate a log
useLog("poisson1dlog")

# Set up the configuration (order doesn't matter)
domain(1)                      # dimension
functionSpace(order=3)         # basis function polynomial order

# Specify the problem (mesh comes first)
mesh(LINEMESH, elsperdim=20)   # build uniform LINEMESH with 20 elements

u = variable("u")              # make a scalar variable with symbol u
testSymbol("v")                # sets the symbol for a test function

boundary(u, 1, DIRICHLET, 0)  # boundary condition for BID 1 is Dirichlet with value 0

# Write the weak form
coefficient("f", "-100*pi*pi*sin(10*pi*x)*sin(pi*x) - pi*pi*sin(10*pi*x)*sin(pi*x) + 20*pi*pi*cos(10*pi*x)*cos(pi*x)")
weakForm(u, "-grad(u)*grad(v) - f*v")

solve(u);

# exact solution is sin(10*pi*x)*sin(pi*x)
# check error
maxerr = 0;
exact(x) = sin(10*pi*x)*sin(pi*x);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    err = abs(u.values[i] - exact(x));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# solution is stored in the variable's "values"
# using Plots
# pyplot();
# display(plot(Femshop.grid_data.allnodes[:], u.values[:], markershape=:circle, legend=false))

# Dump things to the log if desired
log_dump_config();
log_dump_prob();

finalize_femshop() # Finish writing and close any files


text:
@testset "Simple usage of HDF5 in $src with $scheduling_config" for scheduling_config in [
    "default scheduling",
    "parallelism encouraged",
    "parallelism and batches encouraged",
],
src in ["Internet", "S3"]
    use_job_for_testing(scheduling_config_name = scheduling_config) do
        use_data()

        for _ in 1:2
            x = if src == "S3"
                read_hdf5(joinpath("s3://", get_cluster_s3_bucket_name(), "fillval.h5/DS1"))
            else
                read_hdf5(joinpath("https://github.com/banyan-team/banyan-julia/raw/v0.1.1/BanyanArrays/test/res", "fillval.h5/DS1"))
            end

            # Test basic case of reading from remote file
            x_length_collect = length(x)
            @test x_length_collect == 600000
            x_size_collect = size(x)
            @test x_size_collect == (1000, 600)
            x_sum_collect = collect(sum(x))
            @test x_sum_collect == 32100000
            x_sum_collect = collect(sum(x))
            @test x_sum_collect == 32100000
        end
    end
end

@testset "Reading/writing 2D arrays with HDF5 in $src with $scheduling_config" for scheduling_config in
                                                                                [
        "default scheduling",
        "parallelism encouraged",
        "parallelism and batches encouraged",
    ],
    src in ["Internet", "S3"]

    use_job_for_testing(scheduling_config_name = scheduling_config) do
        use_data(src)

        # Determine where to read from

        # TODO: Make this more general by creating S3 bucket and uploading
        # file from test/res for testing
        # TODO: Use version of `pt_lib_info.json` with replication actually removed
        path = if src == "Internet"
            # TODO: Test Internet
            # "https://github.com/banyan-team/banyan-julia/blob/v0.1.1/BanyanArrays/test/res/fillval.h5?raw=true",
            "https://github.com/banyan-team/banyan-julia/raw/v0.1.1/BanyanArrays/test/res/fillval.h5"
        else
            # The file is produced in S3 using:
            # wget https://support.hdfgroup.org/ftp/HDF5/examples/files/exbyapi/h5ex_d_fillval.h5
            # AWS_DEFAULT_PROFILE=banyan-testing aws s3 cp \
            # h5ex_d_fillval.h5 \
            # s3://banyan-cluster-data-pumpkincluster0-3e15290827c0c584/fillval.h5
            # Then, repeat the dataset by (100, 100)
            joinpath("s3://", get_cluster_s3_bucket_name(), "fillval.h5")
        end

        # Perform computation

        x = read_hdf5(joinpath(path, "DS1"))
        @show collect(length(x))
        x = map(e -> e * 10, x)

        steps = if startswith(path, "s3://")
            1:3
        else
            # It's fine to ignore steps 2 and 3 b/c there's no reason I
            # can think of for why reading from Internet and then writing
            # to S3 would not work just because we read from the Internet
            1:1
        end
        for step in steps
            if step == 2
                # TODO: Maybe delete these datasets/files before writing
                # to them to ensure that we don't have issues where we need
                # to fsync the data into S3 before being able to read data.
                # Test writing to and reading from dataset in group in
                # group in same file
                copied_path = joinpath(path, "copies", "DS2")
                write_hdf5(x, copied_path)
                x = read_hdf5(copied_path)
            elseif step == 3
                # Test writing to different file and then reading from it
                copied_path = path[1:end-3] * "_copy.h5/DS1"
                write_hdf5(x, copied_path)
                x = read_hdf5(copied_path)
            end

            # Test basic case of reading from remote file
            x_length_collect = length(x)
            @test x_length_collect == 600000
            x_size_collect = size(x)
            @test x_size_collect == (1000, 600)
            x_sum_collect = collect(sum(x)) # here?
            @test x_sum_collect == 321000000 # incorrect
            x_sum_collect = collect(sum(x))
            @test x_sum_collect == 321000000 # incorrect
            x_minimum_collect = collect(minimum(x))
            @test x_minimum_collect == -60
            x_maximum_collect = collect(maximum(x))
            @test x_maximum_collect == 990
            x_length_collect = length(x)
            @test x_length_collect == 600000
            x_size_collect = size(x)
            @test x_size_collect == (1000, 600)
        end

        # TODO: Re-enable this test once we ensure that we can write out small
        # enough datasets without unnecessary batching. This test involves
        # multiple reductions on the same dataset but the scheduler doesn't
        # think the dataset can stay in memory between evaluations so it tries
        # to persist it to disk. And this fails because we don't support
        # writing HDF5 datasets with strings yet.
        # run_with_job("Reading/writing string arrays with HDF5") do job
        #     for path in [
        #         # "https://github.com/banyan-team/banyan-julia/blob/v0.1.1/BanyanArrays/test/res/vlstring.h5?raw=true",
        #         "https://github.com/banyan-team/banyan-julia/raw/v0.1.1/BanyanArrays/test/res/vlstring.h5",
        #         "s3://banyan-cluster-data-pumpkincluster0-3e15290827c0c584/vlstring.h5",
        #     ]
        #         x = read_hdf5(joinpath(path, "DS1"))
        #         x = map(identity, x)

        #         # TODO: Use all 3 steps so that we can test out writing once
        #         # we get writing strings to work
        #         steps = 1:1
        #         # steps = if startswith(path, "s3://")
        #         #     1:3
        #         # else
        #         #     1:1
        #         # end
        #         for step in steps
        #             if step == 2
        #                 # Test writing to and reading from dataset in group in
        #                 # group in same file
        #                 copied_path = joinpath(path, "copies", "DS2")
        #                 write_hdf5(x, copied_path)
        #                 x = read_hdf5(copied_path)
        #             elseif step == 3
        #                 # Test writing to different file and then reading from it
        #                 copied_path = path[1:end-3] * "_copy.h5/DS1"
        #                 write_hdf5(x, copied_path)
        #                 x = read_hdf5(copied_path)
        #             end

        #             # Test basic case of reading from remote file
        #             x_length_collect = length(x)
        #             @test x_length_collect == 400
        #             x_size_collect = size(x)
        #             @test x_size_collect == (400,)
        #             x_minimum_collect = collect(minimum(x))
        #             @test x_minimum_collect == "Parting"
        #             x_maximum_collect = collect(maximum(x))
        #             @test x_maximum_collect == "sweet"
        #             x_length_collect = length(x)
        #             @test x_length_collect == 400
        #             x_size_collect = size(x)
        #             @test x_size_collect == (400,)
        #         end
        #     end
        # end
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# First generate some random data and test functions in Base on it
const NVALS = 1000

input = Dict(
    t=>[ (randindomain(t, NVALS, domain),) for (_, _, domain) in base_unary_complex ]
        for t in (ComplexF32, ComplexF64)
)

fns = [x[1:2] for x in base_unary_complex]

# output = Dict(
#     t=>[ fns[i](input[t][i]...) for i = 1:length(fns) ]
#         for t in (ComplexF32, ComplexF64)
# )

@testset "Definitions and Comparison with Base for Complex" begin

  for t in (ComplexF32, ComplexF64), i = 1:length(fns)

    base_fn = eval(:($(fns[i][1]).$(fns[i][2])))
    vml_fn = eval(:(IntelVectorMath.$(fns[i][2])))

    Test.@test which(vml_fn, typeof(input[t][i])).module == IntelVectorMath

    # Test.test_approx_eq(output[t][i], fn(input[t][i]...), "Base $t $fn", "IntelVectorMath $t $fn")
    Test.@test vml_fn(input[t][i]...) ≈ base_fn.(input[t][i]...)

  end

end


text:
module ErgodicRadioLocalization

using FEBOL
using ErgodicControl
using ErgodicControlPlots
using Reel

export ErgodicPolicyR2
export ErgodicPolicySE2
export vis2

export NoisyBearing

include("r2.jl")
include("se2.jl")
include("bearing.jl")
include("df2.jl")
include("vis2.jl")

end # module


text:
module Example

using Twinkle, Twinkle.FlexUI
using Colors
using StaticArrays
using GeometryBasics
using Images, ImageSegmentation

println("Start [$(splitext(basename(@__FILE__))[1])]")

app = App()
prop!(app, :title, "Twinkle Example - Images")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

#---------------------------------------------------------------
# Define Variables
#---------------------------------------------------------------

filename = addVariable!(app, Variable(name="filename", type="string",value="C:/t/TestImages/horse.jpg"))
size = addVariable!(app, Variable(name="size", type="int",value=10))
details = addVariable!(app, Variable(name="details", type="string",value=10))

src_image = addVariable!(app, Variable(name="src_image", type="image", value=""))
dst_image = addVariable!(app, Variable(name="dst_image", type="image", value=""))

#---------------------------------------------------------------
# Define Controls
#---------------------------------------------------------------

ui = VContainer(
    H1Label("Segmantation Example"),
    HContainer(
        Button(
            text="Load a file",
            fileTypes=".png,.jpg",
            variable="filename",
        ),
        Label(text="Current file [\$(filename)]"),
    ),
    Container(
        direction = "row warp",
        align = "space-between center",
        children = [
            Slider(
                text="Segments Size",
                trailing_text="[\$()]",
                min=2,
                max=100,
                value=8,
                variable="size"
            ),
            Label(text="\$(details)"),
        ],
    ),
    H1Label("Image Viewer allowing Pan (left-drag) and Zoom (wheel)"),
    PanZoom(
        style="width: 100%; height: 60vh;  border: 1px solid black;",
        content = Container(
            direction = "row",
            children = [
                Card(
                    title="Source Image",
                    # style=card_style,
                    content=VContainer(
                        Image(
                            source="\$(src_image)",
                        ),
                    ),
                ),
                Card(
                    title="Segmented Image",
                    # style=card_style,
                    content=VContainer(
                        Image(
                            source="\$(dst_image)",
                        ),
                    ),
                ),
            ]
        ),
    ),

)
FlexUI.controls!(app, ui)

#---------------------------------------------------------------
# load the image at the start and on name change
#---------------------------------------------------------------
src_img = nothing
function loadImage()
    @info "Loading Image [$(filename[])]"
    global src_img = load(filename[])
end

on(filename) do val
    loadImage()
end

#---------------------------------------------------------------
# the render function
#---------------------------------------------------------------
function processImages()
    global src_img

    segments = felzenszwalb(src_img, size[])
    dst_img = map(i->segment_mean(segments,i), labels_map(segments))
    src_image["png"] = src_img
    dst_image["png"] = dst_img
    details[] = "$segments"
end

function render()
    processImages();
end
renderFunction!(app, render)

#---------------------------------------------------------------
# Run the application
#---------------------------------------------------------------
loadImage()
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module