text:
abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end
################################################################################

# Shampine's Low-order Rosenbrocks

@cache mutable struct Rosenbrock23Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  weight::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

@cache mutable struct Rosenbrock32Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = zero(u) fsalfirst
  f₁ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  weight = similar(u, uEltypeNoUnits)
  tab = Rosenbrock23Tableau(constvalue(uBottomEltypeNoUnits))
  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)

  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2,Val(false))

  Rosenbrock23Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,weight,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = zero(u) fsalfirst
  f₁ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rosenbrock32Tableau(constvalue(uBottomEltypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,linsolve,jac_config,grad_config)
end

struct Rosenbrock23ConstantCache{T,TF,UF,JType,WType,F} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
  J::JType
  W::WType
  linsolve::F
end

function Rosenbrock23ConstantCache(T::Type,tf,uf,J,W,linsolve)
  tab = Rosenbrock23Tableau(T)
  Rosenbrock23ConstantCache(tab.c₃₂,tab.d,tf,uf,J,W,linsolve)
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rosenbrock23ConstantCache(constvalue(uBottomEltypeNoUnits),tf,uf,J,W,linsolve)
end

struct Rosenbrock32ConstantCache{T,TF,UF,JType,WType,F} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
  J::JType
  W::WType
  linsolve::F
end

function Rosenbrock32ConstantCache(T::Type,tf,uf,J,W,linsolve)
  tab = Rosenbrock32Tableau(T)
  Rosenbrock32ConstantCache(tab.c₃₂,tab.d,tf,uf,J,W,linsolve)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rosenbrock32ConstantCache(constvalue(uBottomEltypeNoUnits),tf,uf,J,W,linsolve)
end

################################################################################

### 3rd order specialized Rosenbrocks

struct Rosenbrock33ConstantCache{TF,UF,Tab,JType,WType,F} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
  J::JType
  W::WType
  linsolve::F
end

@cache mutable struct Rosenbrock33Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = ROS3PTableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock33Cache(u,uprev,du,du1,du2,k1,k2,k3,k4,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rosenbrock33ConstantCache(tf,uf,ROS3PTableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

@cache mutable struct Rosenbrock34Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas3Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock34Cache(u,uprev,du,du1,du2,k1,k2,k3,k4,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

struct Rosenbrock34ConstantCache{TF,UF,Tab,JType,WType,F} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
  J::JType
  W::WType
  linsolve::F
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rosenbrock34ConstantCache(tf,uf,Rodas3Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

################################################################################

### ROS34PW methods

@ROS34PW(:cache)

################################################################################

### ROS4 methods

@Rosenbrock4(:cache)
jac_cache(c::Rosenbrock4Cache) = (c.J,c.W)

###############################################################################

### Rodas methods

struct Rodas4ConstantCache{TF,UF,Tab,JType,WType,F} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
  J::JType
  W::WType
  linsolve::F
end

@cache mutable struct Rodas4Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas4Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rodas4ConstantCache(tf,uf,Rodas4Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas42Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rodas4ConstantCache(tf,uf,Rodas42Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas4PTableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rodas4ConstantCache(tf,uf,Rodas4PTableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

function alg_cache(alg::Rodas4P2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas4P2Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4P2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rodas4ConstantCache(tf,uf,Rodas4P2Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

################################################################################

### Rosenbrock5

struct Rosenbrock5ConstantCache{TF,UF,Tab,JType,WType,F} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
  J::JType
  W::WType
  linsolve::F
end

@cache mutable struct Rosenbrock5Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::rateType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  fsalfirst::rateType
  fsallast::rateType
  dT::rateType
  J::JType
  W::WType
  tmp::rateType
  atmp::uNoUnitsType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::rateType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = zero(rate_prototype)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
  tmp = zero(rate_prototype)
  atmp = similar(u, uEltypeNoUnits)
  tab = Rodas5Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

  tf = TimeGradientWrapper(f,uprev,p)
  uf = UJacobianWrapper(f,t,p)
  linsolve_tmp = zero(rate_prototype)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock5Cache(u,uprev,dense1,dense2,du,du1,du2,k1,k2,k3,k4,
                    k5,k6,k7,k8,
                    fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,
                    linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  tf = TimeDerivativeWrapper(f,u,p)
  uf = UDerivativeWrapper(f,t,p)
  J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
  linsolve = alg.linsolve(Val{:init},uf,u)
  Rosenbrock5ConstantCache(tf,uf,Rodas5Tableau(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
end

################################################################################

### RosenbrockW6S4O

@RosenbrockW6S4OS(:cache)


text:
using LinearAlgebra
using Statistics
using SparseArrays
using ProgressMeter
using Plots

function periodic_filter(index_vect::Vector{Int},dim_vect::Vector{Int})::Nothing
    """Indexes are updated periodically so that overflow and underflow situations are understood with respect to the tensors' dimensions"""
    for i in 1:length(index_vect)
        index_vect[i] = mod(index_vect[i]-1, dim_vect[i]) + 1
    end
    nothing
end

#Hash access function
function hash_function(index_vect::Vector{Int},dim_vect::Vector{Int},periodic_flag = true)::Int
    """Hash function to flatten tensors length(dim_vect)-tensors"""
    @assert length(index_vect) == length(dim_vect)
    my_index_vect = copy(index_vect)
    if periodic_flag
        periodic_filter(my_index_vect, dim_vect)
    end
    l = length(my_index_vect)
    s = 0
    for i in 1:l
        s = (my_index_vect[i] - 1) + s * dim_vect[i]
    end
    return s + 1 #Problem?
end

function hash_function(index_vect::Vector{Int},dim::Int)::Int
    """Hash function to flatten tensors length(dim_vect)-tensors, square case"""
    return hash_function(index_vect, dim * ones(Int,length(index_vect)))
end

function hash_function(i::Int,j::Int,dim_vect::Vector{Int})::Int
    """Hash function to flatten tensors length(dim_vect)-tensors, dimension 2 case"""
    return j + dim_vect[2]*i
end

#Terrain building, centered around 0
function vector_dimension(dim_vect::Vector{Int})::Int
    """Total dimension of tensor of a given dimension as a column vector."""
    dim = 1
    for i in dim_vect
        dim *= i
    end
    return dim
end

function multidimensional_indexes(dim_vect::Vector{Int})::Vector{Vector{Int}}
    """Supposed to get list of indexes possible given a dimensions' vector."""
    multi = []
    if length(dim_vect) == 1
        for i in 1:dim_vect[1]
            append!(multi,[[i]])
        end
        return multi
    else
        pred = multidimensional_indexes(dim_vect[1:(length(dim_vect)-1)])
        for i in 1:dim_vect[length(dim_vect)]
            for pred_ind in pred
                acc = copy(pred_ind)
                append!(acc,[i])
                append!(multi,[acc])
            end
        end
        return multi
    end
end

function terrain_builder(dim_vect::Vector{Int}, step_vect::Vector{Float64}, corner_point::Vector{Float64})::Array{Vector{Float64}}
    """Builds a terrain vector of a given step vector in the shape of cube which inferior corner is given."""
    terrain_vector = Array{Vector{Float64}}(undef, vector_dimension(dim_vect))
    possible_indexes = multidimensional_indexes(dim_vect)
    for index_vect in possible_indexes
        terrain_vector[hash_function(index_vect,dim_vect)] = index_vect .* step_vect + corner_point
    end
    return terrain_vector
end

#Matrix building
function first_order_differentiation_rigidity(rigidity::AbstractArray,index_vect::Vector{Int},dim_vect::Vector{Int},diff_direction::Int,step_in_direction::Float64, prefactor = 1.::Float64)::Nothing
    """Computes a matrix corresponding to a first order derivation in a periodical context. Side-effect action."""
    index_tilt = zeros(Int, length(index_vect))
    index_tilt[diff_direction] = 1

    index_vect_plus = index_vect + index_tilt
    index_vect_minus = index_vect - index_tilt

    i = hash_function(index_vect,dim_vect)
    j1 = hash_function(index_vect_plus,dim_vect)
    j2 = hash_function(index_vect_minus,dim_vect)

    rigidity[i,j1] += prefactor * 1/(2step_in_direction)
    rigidity[i,j2] += prefactor * (- 1/(2step_in_direction))

    nothing
end

function second_order_differentiation_rigidity(rigidity::AbstractArray,index_vect::Vector{Int},dim_vect::Vector{Int},diff_direction::Int,step_in_direction::Float64, prefactor = 1.::Float64)::Nothing
    """Computes a matrix corresponding to a second order derivation in a periodical context. Side-effect action."""
    index_tilt = zeros(Int, length(index_vect))
    index_tilt[diff_direction] = 1

    index_vect_plus = index_vect + index_tilt
    index_vect_minus = index_vect - index_tilt

    i = hash_function(index_vect,dim_vect)
    j1 = hash_function(index_vect_plus,dim_vect)
    j2 = hash_function(index_vect_minus,dim_vect)

    rigidity[i,j1] += prefactor * 1/(step_in_direction*step_in_direction)
    rigidity[i,j2] += prefactor * 1/(step_in_direction*step_in_direction)
    rigidity[i,i] += prefactor * -2/(step_in_direction*step_in_direction)

    nothing
end

#Matrix composition
function identity_matrix(matrix_dimension::Int,sparse_flag = true::Bool)::AbstractMatrix
    if sparse_flag
        identity = spzeros(Float64, matrix_dimension, matrix_dimension)
        for i in 1:matrix_dimension
            identity[matrix_dimension,matrix_dimension] = 1
        end
        return identity
    else
        return Matrix(1.0I, matrix_dimension, matrix_dimension)
    end
end

function total_regidity(terrain::Vector{Vector{Float64}}, step_vect::Vector{Float64}, D::Vector{Float64},C::Matrix{Float64},dim_vect::Vector{Int})::AbstractArray
    """Computes the total rigidity of the PDE problem"""
    tensor_dim = vector_dimension(dim_vect)
    rigidity = zeros(Float64, tensor_dim, tensor_dim)
    index_list = multidimensional_indexes(dim_vect)
    @showprogress "Assemblage de la rigidité" for index_vect in index_list
        current_position = C * terrain[hash_function(index_vect,dim_vect)]
        for i in 1:length(dim_vect)
            second_order_differentiation_rigidity(rigidity,index_vect, dim_vect, i, step_vect[i],D[i])
            first_order_differentiation_rigidity(rigidity, index_vect, dim_vect, i, step_vect[i],current_position[i])
        end
    end
    return rigidity + tr(C) * identity_matrix(tensor_dim,false)
end

#Gradient caluculations
function tensor_gradient(U::Vector{Float64},dim_vect::Vector{Int},step_vect::Vector{Float64})
    """Computes gradient through finite differences."""
    index_list = multidimensional_indexes(dim_vect)
    gradU = [zeros(length(dim_vect)) for k in 1:length(U)]
    for index_vect in index_list
        i = hash_function(index_vect,dim_vect)
        grad_vector = gradU[i]
        for diff_direction in 1:length(dim_vect)
            index_tilt = zeros(Int, length(index_vect))
            index_tilt[diff_direction] = 1

            index_vect_plus = index_vect + index_tilt
            index_vect_minus = index_vect - index_tilt

            j1 = hash_function(index_vect_plus,dim_vect)
            j2 = hash_function(index_vect_minus,dim_vect)

            grad_vector[diff_direction] = (U[j1] - U[j2])/(step_vect[diff_direction])
        end
        gradU[i] = grad_vector
    end
    return gradU
end

function tensor_gradient(U::Vector{Float64},dim_vect::Vector{Int},step::Float64)
    """Computes gradient through finite differences."""
    step_vect = step * ones(Float64,length(dim_vect))
    return tensor_gradient(U,dim_vect,step_vect)
end

#P - norms
function custom_norm_squared(vector::Vector{Float64}, matrix::Matrix{Float64})::Float64
    """Computes tXMX for given X and M."""
    return transpose(vector) * matrix * vector
end

#Derichlet Setup
function is_extremal(index_vect::Vector{Int},dim_vect::Vector{Int})::Bool
    """Checks if coordinate vector cooresponds to boundary point."""
    d = length(index_vect)
    @assert d == length(dim_vect)
    if d == 1
        return (index_vect[1] == 1)||(index_vect[1] == dim_vect[1])
    else
        return (index_vect[1] == 1)||(index_vect[1] == dim_vect[1]) || is_extremal(index_vect[2:d],dim_vect[2:d])
    end
end

function derichlet_reset(U::Vector{Float64},dim_vect::Vector{Int})::Nothing
    """Sets boundary points to 0 to enforce Derichlet conditions for the PDE. Side effect action."""
    index_list = multidimensional_indexes(dim_vect)
    for index_vect in index_list
        if is_extremal(index_vect,dim_vect)
            U[hash_function(index_vect,dim_vect)] = 0
        end
    end
    nothing
end

#show stuff
function show_2D_case(vector::Vector{Float64}, dim_vect::Vector{Int})
    """ Displays heatmap of function. """
    mat = reshape(vector,dim_vect[1],dim_vect[2])
    xgrid = collect(1:dim_vect[1])
    vgrid = collect(1:dim_vect[2])
    p = heatmap(xgrid,vgrid,mat)
end

#Initialisation
function terrain_parameters_initialisation(lattice_parameter::Vector{Float64},dim_vect::Vector{Int})
    """ Returns a terrain vector and the spacewise stepvector."""
    step_vect = lattice_parameter ./ dim_vect
    corner_point = - lattice_parameter /2
    return (terrain_builder(dim_vect,step_vect,corner_point), step_vect)
end

function terrain_parameters_initialisation(L::Float64, dim::Int, d::Int)
    """ Returns a terrain vector and the spacewise stepvector, in the case of a square terrain.
    Dim is the number of steps and d the number of physical dimensions of the PDE."""
    lattice_parameter = L* ones(Float64,d)
    dim_vect = dim * ones(Int,d)
    return terrain_parameters_initialisation(lattice_parameter,dim_vect)
end

function gaussian(entry_vector::Vector{Float64}, mean = 0.::Float64)::Float64
    mean_v = mean * ones(Float64,length(entry_vector))
    return exp(-(norm(entry_vector- mean_v)^2)/2)
end

function gaussian(entry_vector::Vector{Float64}, mean::Vector{Float64})::Float64
    return gaussian(entry_vector - mean)
end

#Special Handling
get_last_column(my_array) = my_array[:,length(my_array[1,:])]

maxa(array::AbstractArray) = reduce((x,y) -> max(x,y), array)

min(array::AbstractArray) = reduce((x,y) -> min(x,y), array)

function solve()
    n = 1
    d = 2n
    C = [0. 0. ; 0. 0.]
    D = [1.,1.]

    L = 15.
    dim = 100

    lattice_parameter = L* ones(Float64,d)
    dim_vect = dim * ones(Int,d)
    terrain, step_vect = terrain_parameters_initialisation(L, dim, d)
    println("Terrain done")

    U = gaussian.(terrain)
    rigidity = total_regidity(terrain, step_vect,D,C,dim_vect)
    println("Rigidity done")

    #######
    T = 1.
    J = 25
    time_step = T/J

    K = inv( Matrix(1.0I, vector_dimension(dim_vect), vector_dimension(dim_vect)) - time_step * rigidity )
    println("Rigidity inverted")

    U_tab = copy(U)
    @gif for j in 1:J
        if j%5 == 0
            @show j
        end
        show_2D_case(get_last_column(U_tab),dim_vect)
        U_new = K * get_last_column(U_tab)
        U_tab = [U_tab U_new]
    end
end

function debug_diff()
    n = 1
    d = 2n
    C = [0. 0. ; 0. 0.]
    D = [1.,1.]

    L = 15.
    dim = 3

    lattice_parameter = L* ones(Float64,d)
    dim_vect = dim * ones(Int,d)
    terrain, step_vect = terrain_parameters_initialisation(L, dim, d)
    println("Terrain done")

    U = gaussian.(terrain)
    rigidity = total_regidity(terrain, step_vect,D,C,dim_vect)
    println("Rigidity done")

    println("Sum of coefficients = $( sum(rigidity) )")
    println("Matrix is symmetric? $( issymmetric(rigidity) )")
end

ext: jl

max_stars_repo_name: RJ-McReady/MOPSI-Hypocoercivity

max_issues_repo_name: RJ-McReady/MOPSI-Hypocoercivity

max_forks_repo_name: RJ-McReady/MOPSI-Hypocoercivity


text:
# --- above === below ---
"
    This object manages the periodic mpo environments for an MPSMultiline
"
mutable struct PerMPOInfEnv{H<:MPOMultiline,V,S<:MPSMultiline} <: AbstractInfEnv
    opp :: H

    dependency :: S
    tol :: Float64
    maxiter :: Int

    lw :: PeriodicArray{V,2}
    rw :: PeriodicArray{V,2}

    lock :: ReentrantLock
end

function recalculate!(envs::PerMPOInfEnv,nstate::MPSMultiline)
    sameDspace = reduce((prev,i) -> prev && _lastspace(envs.lw[i...]) == _firstspace(nstate.CR[i...])',
        Iterators.product(1:size(nstate,1),1:size(nstate,2)),init=true);

    init = collect(zip(envs.lw[:,1],envs.rw[:,end]))
    if !sameDspace
        init = gen_init_fps(nstate,envs.opp,nstate)
    end

    (envs.lw,envs.rw) = mixed_fixpoints(nstate,envs.opp,nstate,init,tol = envs.tol, maxiter = envs.maxiter);
    envs.dependency = nstate;

    envs
end

environments(state::InfiniteMPS,opp::InfiniteMPO;kwargs...) = environments(convert(MPSMultiline,state),convert(MPOMultiline,opp);kwargs...);
function environments(state::MPSMultiline,mpo::MPOMultiline;tol = Defaults.tol,maxiter=Defaults.maxiter)
    (lw,rw) = mixed_fixpoints(state,mpo,state;tol = tol, maxiter = maxiter)

    PerMPOInfEnv(mpo,state,tol,maxiter,lw,rw,ReentrantLock())
end

mutable struct MixPerMPOInfEnv{H<:MPOMultiline,V,S<:MPSMultiline} <: AbstractInfEnv
    opp :: H

    above :: S
    dependency :: S

    tol :: Float64
    maxiter :: Int

    lw :: PeriodicArray{V,2}
    rw :: PeriodicArray{V,2}

    lock :: ReentrantLock
end

function recalculate!(envs::MixPerMPOInfEnv,nstate::MPSMultiline)
    sameDspace = reduce((prev,i) -> prev && _firstspace(envs.lw[i...]) == _firstspace(nstate.CR[i...]),
        Iterators.product(1:size(nstate,1),1:size(nstate,2)),init=true);

    init = collect(zip(envs.lw[:,1],envs.rw[:,end]))
    if !sameDspace
        init = gen_init_fps(envs.above,envs.opp,nstate)
    end

    (envs.lw,envs.rw) = mixed_fixpoints(envs.above,envs.opp,nstate,init,tol = envs.tol, maxiter = envs.maxiter);
    envs.dependency = nstate;

    envs
end

function environments(below::InfiniteMPS,toapprox::Tuple{<:InfiniteMPO,<:InfiniteMPS};kwargs...)
    (opp,above) = toapprox
    environments(convert(MPSMultiline,below),(convert(MPOMultiline,opp),convert(MPSMultiline,above));kwargs...);
end
function environments(below::MPSMultiline,toapprox::Tuple{<:MPOMultiline,<:MPSMultiline};tol = Defaults.tol,maxiter=Defaults.maxiter)
    (mpo,above) = toapprox;
    (lw,rw) = mixed_fixpoints(above,mpo,below;tol = tol, maxiter = maxiter)

    MixPerMPOInfEnv(mpo,above,below,tol,maxiter,lw,rw,ReentrantLock())
end

# --- utility functions ---

function gen_init_fps(above::MPSMultiline,mpo::MPOMultiline,below::MPSMultiline)
    T = eltype(above)

    map(1:size(mpo,1)) do cr
        L0::T = TensorMap(rand,eltype(T),space(below.AL[cr,1],1)*space(mpo[cr,1],1)',space(above.AL[cr,1],1))
        R0::T = TensorMap(rand,eltype(T),space(above.AR[cr,1],1)*space(mpo[cr,1],1),space(below.AR[cr,1],1))
        (L0,R0)
    end
end

function mixed_fixpoints(above::MPSMultiline,mpo::MPOMultiline,below::MPSMultiline,init = gen_init_fps(above,mpo,below);tol = Defaults.tol, maxiter = Defaults.maxiter)
    T = eltype(above);

    #sanity check
    (numrows,numcols) = size(above)
    @assert size(above) == size(mpo)
    @assert size(below) == size(mpo);

    lefties = PeriodicArray{T,2}(undef,numrows,numcols);
    righties = PeriodicArray{T,2}(undef,numrows,numcols);

    @sync for cr = 1:numrows
        @Threads.spawn begin
            (L0,R0) = $init[cr]

            (_,Ls,convhist) = eigsolve(x-> transfer_left(x,$mpo[cr,:],$above.AL[cr,:],$below.AL[cr+1,:]),L0,1,:LM,Arnoldi(tol = tol,maxiter=maxiter))
            convhist.converged < 1 && @info "left eigenvalue failed to converge $(convhist.normres)"
            (_,Rs,convhist) = eigsolve(x-> transfer_right(x,$mpo[cr,:],$above.AR[cr,:],$below.AR[cr+1,:]),R0,1,:LM,Arnoldi(tol = tol,maxiter=maxiter))
            convhist.converged < 1 && @info "right eigenvalue failed to converge $(convhist.normres)"

            $lefties[cr,1] = Ls[1]
            for loc in 2:numcols
                $lefties[cr,loc] = transfer_left($lefties[cr,loc-1],$mpo[cr,loc-1],$above.AL[cr,loc-1],$below.AL[cr+1,loc-1])
            end

            renormfact::eltype(T) = @tensor Ls[1][1,2,3]*above.CR[cr,0][3,4]*Rs[1][4,2,5]*conj(below.CR[cr+1,0][1,5])

            $righties[cr,end] = Rs[1]/sqrt(renormfact);
            $lefties[cr,1] /=sqrt(renormfact);

            for loc in numcols-1:-1:1
                $righties[cr,loc] = transfer_right($righties[cr,loc+1],$mpo[cr,loc+1],$above.AR[cr,loc+1],$below.AR[cr+1,loc+1])

                renormfact = @tensor lefties[cr,loc+1][1,2,3]*above.CR[cr,loc][3,4]*righties[cr,loc][4,2,5]*conj(below.CR[cr+1,loc][1,5])
                $righties[cr,loc]/=sqrt(renormfact)
                $lefties[cr,loc+1]/=sqrt(renormfact)
            end
        end
    end

    return (lefties,righties)
end


text:
# Based on the enumeration of all 7,462 Five-Card Poker Hand Equivalence Classes

using .HandCombinations

"""
    evaluate5(cards::Card...)
    evaluate5(::Card,::Card,::Card,::Card,::Card)
    evaluate5(::Tuple{Card,Card,Card,Card,Card})

This is PokerHandEvaluator.jl's _core_ method.

Returns the `rank` (an `Int`) from `1:7462`:

    (A♡,K♡,Q♡,J♡,10♡) -> 1
    ...
    (7♡,5♢,4♣,3♠,2♡) -> 7462

for 5-card hands _only_.
"""
function evaluate5 end

evaluate5(cards::Card...)::Int = evaluate5(cards)
# The product of 5 primes will be unique, and
# do not depend on the order, so we're using
# this to map card combinations to a single number
# for which a method, `evaluate_offsuit` and
# `evaluate_flush` are defined. First, we dispatch
# based on flush/non-flush:
function evaluate5(t::NTuple{N,Card})::Int where {N}
    @assert N == 5
    if suit(t[1]) == suit(t[2]) == suit(t[3]) == suit(t[4]) == suit(t[5])
        return hash_table_suited[prod(prime.(t))]
    else
        return hash_table_offsuit[prod(prime.(t))]
    end
end

const hash_table_suited = Dict{Int,Int}(hcat(
    map(enumerate(straight_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), i)       # Rows 1:10 (Straight flush)
    end...,
    map(enumerate(flush_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 323+i-1) # Rows 323:1599 (flush)
    end...,
))

const hash_table_offsuit = Dict{Int,Int}(hcat(
    map(enumerate(quad_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 11+i-1)   # Rows 11:166 (4 of a kind)
    end...,
    map(enumerate(full_house_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 167+i-1)  # Rows 167:322 (full house)
    end...,
    map(enumerate(straight_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 1600+i-1) # Rows 1600:1609 (off-suit straight)
    end...,
    map(enumerate(trip_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 1610+i-1) # Rows 1610:2467 (trips)
    end...,
    map(enumerate(two_pair_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 2468+i-1) # Rows 2468:3325 (two pair)
    end...,
    map(enumerate(pair_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 3326+i-1) # Rows 3326:6185 (pair)
    end...,
    map(enumerate(high_card_ranks())) do (i,card_ranks)
        Pair(prod(prime.(card_ranks)), 6186+i-1) # Rows 6186:7462 (high card)
    end...,
))


text:
using ThinkJulia

turtle = Turtle()
@svg begin
    forward(turtle, 100)
end

turtle = Turtle()
@svg begin
    forward(turtle, 100)
    turn(turtle, -90)
    forward(turtle, 100)
end

turtle = Turtle()
@svg begin
    for i in 1:4
        forward(turtle, 100)
        turn(turtle, -90)
    end
end

function polygon(t, n, len)
    angle = 360 / n
    for i in 1:n
        forward(t, len)
        turn(t, -angle)
    end
end

function circle(t, r)
    circumference = 2 * pi * r
    n = trunc(circumference / 3) + 3
    len = circumference / n
    polygon(t, n, len)
end

function polyline(t, n, len, angle)
    for i in 1:n
        forward(t, len)
        turn(t, -angle)
    end
end

function polygon(t, n, len)
    angle = 360 / n
    polyline(t, n, len, angle)
end

function arc(t, r, angle)
    arc_len = 2 * π * r * angle / 360
    n = trunc(arc_len / 3) + 1
    step_len = arc_len / n
    step_angle = angle / n
    polyline(t, n, step_len, step_angle)
end

function circle(t, r)
    arc(t, r, 360)
end


text:
using Test
using TypeClasses

# MonoidAlternative
# =================

@test neutral(Dict) == Dict()
@test neutral(Dict{String, Int}) == Dict{String, Int}()

# we can combine any Dict where elements can be combined
d1 = Dict(:a => "3", :b => "1")
d2 = Dict(:a => "5", :b => "9", :c => "15")

@test d1 ⊕ d2 == Dict(:a => "35", :b => "19", :c => "15")

@test d1 ⊘ d2 == Dict(:a => "3", :b => "1", :c => "15")

# FlipTypes
# =========

@test flip_types(Dict(:a => [1,2], :b => [3, 4])) == [
    Dict(:a => 1, :b => 3),
    Dict(:a => 1, :b => 4),
    Dict(:a => 2, :b => 3),
    Dict(:a => 2, :b => 4),
]

# the other way arround does not work, as Dict does not implement Functor/Applicative/Monad interfaces
# this is not a MethodError, but a special exception thrown intentionally by Base
@test_throws ErrorException flip_types([d1, d2])


text:
export Matsuno_fast,Matsuno

"""
	Matsuno_fast(param;kwargs)

Same as `Matsuno`, but faster.
"""
mutable struct Matsuno_fast <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info	:: String

    function Matsuno_fast(param::NamedTuple;
							mesh = Mesh(param),
							IL	    = false,
							dealias = false,
							label = "Matsuno" )

		# Set-up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used thereafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "Matsuno model.\n"
		if IL == true
			info *= "├─Steepness parameter ϵ=$ϵ (infinite depth case).\n"
		else
			info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
					scaling parameter ν=$nu.\n"
		end
		if dealias == true || dealias == 1
			info *= "└─Dealiasing with Orszag’s 3/2 rule. "
		else
			info *= "└─No dealiasing. "
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data
		x = mesh.x
		k = mesh.k
		if IL == true
			Tμ 	= -1im * sign.(k)
			G₀ 	= abs.(k)
		else
			Tμ 	= -1im* sign.(k) .* tanh.(sqrt(μ)*abs.(k))
			G₀ 	= sqrt(μ)*abs.(k).*tanh.(sqrt(μ)*abs.(k))
		end
		∂ₓ	=  1im * sqrt(μ)* k            # Differentiation
		if dealias == true || dealias == 1
			Π⅔    = abs.(k) .< (mesh.kmax-mesh.kmin)/3    # Dealiasing low-pass filter
		else
			Π⅔    = zero(k) .+ 1     		# No dealiasing (Π⅔=Id)
		end

        ζ = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₀ = zeros(Complex{Float64}, mesh.N)
        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

        Px  = plan_fft(ζ; flags = FFTW.MEASURE)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		    for i in eachindex(ζ)
		        ζ[i] = G₀[i] * U[i,1]
		    end

		    ldiv!(unew, Px, ζ )

		    for i in eachindex(ζ)
		        ζ[i] = ∂ₓ[i] * U[i,1]
		    end

		    ldiv!(I₁, Px, ζ)

		    unew  .*= I₁

		    mul!(I₁, Px, unew)

		    I₁  .*= ϵ .* Π⅔
		    I₁  .-= ζ

		    ldiv!(ζ, Px, view(U,:,1))
		    ldiv!(unew, Px, view(U,:,2))

		    I₂    .= ζ .* unew

		    mul!(I₃, Px, I₂)

		    I₃    .*= ∂ₓ

		    for i in eachindex(Tμ)
		        U[i,1]  = Tμ[i] * U[i,2]
		        I₀[i] = G₀[i] * U[i,2]
		    end

		    ldiv!(I₂, Px, I₀)

		    I₂    .*= ζ

		    mul!(ζ, Px, I₂)

		    ζ  .*= Tμ
		    I₃    .+= ζ
		    I₃    .*= ϵ .* Π⅔

		    for i in eachindex(I₃)
		        U[i,1] -= I₃[i]
		    end
			U[:,1] ./= sqrt(μ)/ν

		    I₃    .=  unew.^2

		    mul!(unew, Px, I₃)

		    unew  .*= ∂ₓ
		    unew  .*= ϵ/2/ν .* Π⅔
		    I₁    .-= unew

		    for i in eachindex(I₁)
		        U[i,2] =  I₁[i]/sqrt(μ)
		    end

		end

		# Build raw data from physical data.
		function mapto(data::InitialData)
			fftη = Π⅔ .* fft(data.η(x));
			fftv = Π⅔ .* fft(data.v(x));
			U = [fftη fftv-ϵ* Π⅔ .*fft(ifft(Tμ.*fftv).*ifft(∂ₓ.*fftη) )]
			return U
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U;n=10)
			∂ζ=ifft(∂ₓ.*U[:,1]);
			I₁.=view(U,:,2);I₂.=view(U,:,2);
			for j=1:n
				I₂.=I₁+ϵ*Π⅔ .* fft( ∂ζ .* ifft(Tμ.*I₂))
			end
			real(ifft(view(U,:,1))),real(ifft(I₂)),mesh.x
		end

		new(label, f!, mapto, mapfro, info )
    end
end

"""
	Matsuno(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the quadratic deep-water model proposed by [Matsuno](https://doi.org/10.1103/PhysRevLett.69.609).

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- optionally, `ν` the shallow/deep water scaling factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0`, or `μ=Inf`.
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- `IL`: Set the infinite-layer case if `IL=true` (or `μ=Inf`, or `ν=0`), in which case `ϵ` is the steepness parameter. Default is `false`.
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `dealias`: dealiasing with `1/3` Orlicz rule if `true` or no dealiasing if `false` (by default);
- `label`: a label for future references (default is `"Matsuno"`);

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `Matsuno.f!` to be called in explicit time-integration solvers;
2. a function `Matsuno.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `Matsuno.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct Matsuno <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info 	:: String

    function Matsuno(param::NamedTuple ;
						mesh = Mesh(param),
						IL	    = false,
						dealias = false,
						label = "Matsuno" )

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used thereafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "Matsuno model.\n"
		if IL == true
			info *= "├─Steepness parameter ϵ=$ϵ (infinite depth case).\n"
		else
			info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
					scaling parameter ν=$nu.\n"
		end
		if dealias == true || dealias == 1
			info *= "└─Dealiasing with Orszag’s 3/2 rule. "
		else
			info *= "└─No dealiasing. "
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data
		x = mesh.x
		k = mesh.k
		if IL == true
			Tμ 	= -1im * sign.(k)
			G₀ 	= abs.(k)
		else
			Tμ 	= -1im* sign.(k) .* tanh.(sqrt(μ)*abs.(k))
			G₀ 	= sqrt(μ)*abs.(k).*tanh.(sqrt(μ)*abs.(k))
		end
		∂ₓ	=  1im * sqrt(μ)* k            # Differentiation

		if dealias == true || dealias == 1
			Π⅔    = abs.(k) .< (mesh.kmax-mesh.kmin)/3 	# Dealiasing low-pass filter
		else
			Π⅔    = zero(k) .+ 1     		# No dealiasing (Π⅔=Id)
		end

        ζ = zeros(Complex{Float64}, mesh.N)
        unew = zeros(Complex{Float64}, mesh.N)

        I₁ = zeros(Complex{Float64}, mesh.N)
        I₂ = zeros(Complex{Float64}, mesh.N)
        I₃ = zeros(Complex{Float64}, mesh.N)

		# Evolution equations are ∂t U = f(U)
		function f!(U)

		   ζ .= ifft(U[:,1])
		   unew .= ifft(U[:,2])
		   I₃ .= fft(ifft(∂ₓ.*U[:,1]).*ifft(G₀.*U[:,1]))
		   I₁ .= Tμ.*U[:,2].-ϵ*Π⅔.*(Tμ.*fft(ζ.*ifft(G₀.*U[:,2])).+∂ₓ.*fft(ζ.*unew))
		   I₂ .= -ν*(∂ₓ.*U[:,1])+ϵ*Π⅔.*(ν*I₃-∂ₓ.*fft(unew.^2)/2)
		   #
		   U[:,1] .= I₁/sqrt(μ)/ν
		   U[:,2] .= I₂/sqrt(μ)/ν

		end

		# Build raw data from physical data.
		function mapto(data::InitialData)
			fftη = Π⅔ .* fft(data.η(x));
			fftv = Π⅔ .* fft(data.v(x));
			U = [fftη fftv-ϵ* Π⅔ .*fft(ifft(Tμ.*fftv).*ifft(∂ₓ.*fftη) )]
			return U
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U;n=10)
			∂ζ=ifft(∂ₓ.*U[:,1]);
			I₁.=U[:,2];I₂.=U[:,2];
			for j=1:n
				I₂.=I₁+ϵ*Π⅔ .* fft( ∂ζ .* ifft(Tμ.*I₂))
			end
			real(ifft(U[:,1])),real(ifft(I₂)),mesh.x
		end

		new(label, f!, mapto, mapfro, info )
    end
end


text:
"""
    abstract type Triangulation{Dc,Dp}

Abstract type representing an arbitrary tiling, tessellation,
or triangulation of a domain of parametric dimension `Dc` and
physical dimension `Dp`.

We define a triangulation from two basic ingredients:

- the cell-wise nodal coordinates of the cells in the triangulation, plus
- an interpolation of this cell-wise coordinates into the cells interior.

Note that this type represents general triangulations (not necessarily conforming),
which is the minimum geometrical information needed to perform cell-wise numerical integration.

The `Triangulation` interface is defined by overloading these methods:

- [`get_cell_coordinates(trian::Triangulation)`](@ref)
- [`get_reffes(trian::Triangulation)`](@ref)
- [`get_cell_type(trian::Triangulation)`](@ref)

Optional interface:

For triangulations living in a space of co-dimension 1, the following method can be defined:

- [`get_facet_normal(trian::Triangulation)`]

In some cases, concrete implementations want to override the default implementation of the following methods:

- [`restrict(f::AbstractArray, trian::Triangulation)`]
- [`get_cell_to_bgcell(f::AbstractArray, trian::Triangulation)`]

The (mandatory) `Triangulation` interface can be tested with

- [`test_triangulation`](@ref)

"""
abstract type Triangulation{Dc,Dp} <: GridapType end

"""
    get_cell_coordinates(trian::Triangulation) -> AbstractArray{Vector{<:Point{Dp}}}
"""
function get_cell_coordinates(trian::Triangulation)
  @abstractmethod
end

"""
    get_reffes(trian::Triangulation) -> Vector{LagrangianRefFE}
"""
function get_reffes(trian::Triangulation)
  @abstractmethod
end

"""
    get_cell_type(trian::Triangulation) -> AbstractVector{<:Integer}
"""
function get_cell_type(trian::Triangulation)
  @abstractmethod
end

"""
    get_facet_normal(trian::Triangulation)
"""
function get_facet_normal(trian::Triangulation)
  Dp = num_point_dims(trian)
  Dc = num_cell_dims(trian)
  if Dp == Dc + 1
    @abstractmethod
  else
    @unreachable "get_facet_normal does not make sense for this triangulation"
  end
end

"""
If the reference (resp. physical) space is the same for both triangulaitons.
"""
have_compatible_domains(a::Triangulation,b::Triangulation) = a===b

# Trait that signals if the triangulation is a sub-mesh of a background triangulation

abstract type TriangulationStyle end
struct BackgroundTriangulation <: TriangulationStyle end
struct SubTriangulation <: TriangulationStyle end

TriangulationStyle(::Type{<:Triangulation}) = BackgroundTriangulation()
TriangulationStyle(::T) where T<:Triangulation = TriangulationStyle(T)

# Default methods needed to be overloaded by Triangulations that are integration sub-meshes

"""
Returns the background triangulation.
"""
get_background_triangulation(trian::Triangulation) =
  get_background_triangulation(trian,TriangulationStyle(trian))
get_background_triangulation(trian::Triangulation,::BackgroundTriangulation) = trian
get_background_triangulation(trian::Triangulation,::SubTriangulation) = @abstractmethod

#"""
#    reindex(a::AbstractArray, trian::Triangulation)
#"""
#function reindex(a::AbstractArray,trian::Triangulation)
#  reindex(a,get_cell_to_bgcell(trian))
#end

"""
    get_cell_to_bgcell(trian::Triangulation)

Map from the indices in the sub-triangulation to the indices in the background triangulation
"""
get_cell_to_bgcell(trian::Triangulation) = get_cell_to_bgcell(trian,TriangulationStyle(trian))
get_cell_to_bgcell(trian::Triangulation,::BackgroundTriangulation) = IdentityVector(num_cells(trian))
get_cell_to_bgcell(trian::Triangulation,::SubTriangulation) = @abstractmethod
function get_cell_to_bgcell(trian::Triangulation,bgtrian::Triangulation)
  if have_compatible_domains(get_background_triangulation(trian),bgtrian)
    return get_cell_to_bgcell(trian)
  else
    @notimplemented
  end
end

function is_included(trian::Triangulation,bgtrian::Triangulation)
  if have_compatible_domains(get_background_triangulation(trian),bgtrian)
    true
  else
    false
  end
end

#"""
#    restrict(f::AbstractArray, trian::Triangulation)
#"""
#function restrict(f::AbstractArray,trian::Triangulation)
#  f
#end

"""
Return the cell-wise map that goes from the reference space of the sub-triangulation to
the reference space of the background triangulation
"""
get_cell_ref_map(trian::Triangulation) = get_cell_ref_map(trian,TriangulationStyle(trian))
get_cell_ref_map(trian::Triangulation,::BackgroundTriangulation) = Fill(GenericField(identity),num_cells(trian))
get_cell_ref_map(trian::Triangulation,::SubTriangulation) = @abstractmethod
function get_cell_ref_map(trian::Triangulation,bgtrian::Triangulation)
  if have_compatible_domains(get_background_triangulation(trian),bgtrian)
    return get_cell_ref_map(trian)
  else
    @notimplemented
  end
end
#"""
#Given an array aligned with the cells in the background triangulation, return another array
#aligned with the cells of the sub-triangulation. Do nothing if `trian` is already a
#background triangulation.
#"""
#change_cell_index(a::AbstractArray,trian::Triangulation) = change_cell_index(a,trian,TriangulationStyle(trian))
#function change_cell_index(a::AbstractArray,trian::Triangulation,::BackgroundTriangulation)
#  @assert length(a) == num_cells(trian)
#  a
#end
#function change_cell_index(a::AbstractArray,trian::Triangulation,::SubTriangulation)
#  bgtrian = get_background_triangulation(trian)
#  @assert length(a) == num_cells(bgtrian)
#  lazy_map(Reindex(a),get_cell_to_bgcell(trian))
#end
#
#"""
#Given an array (of arrays) of Field objects that is aligned with the sub-triangulation but
#the domain space of those fields corresponds to the reference space of the background triangulation,
#return another array (of arrays) of Field objects with domain in the sub-triangulation.
#Do nothing if `trian` is already a background triangulation.
#Note that one typically needs to call `change_cell_index` before calling `change_cell_domain`.
#"""
#change_cell_domain(a::AbstractArray,trian::Triangulation) = change_cell_domain(a,trian,TriangulationStyle(trian))
#function change_cell_domain(a::AbstractArray,trian::Triangulation,::BackgroundTriangulation)
#  @assert length(a) == num_cells(trian)
#  a
#end
#function change_cell_domain(a::AbstractArray,trian::Triangulation,::SubTriangulation)
#  @assert length(a) == num_cells(trian)
#  cell_ref_map = get_cell_ref_map(trian)
#  lazy_map(Broadcasting(∘),a,cell_map)
#end

"""
    struct SkeletonPair{L,R} <: GridapType
      plus::L
      minus::R
    end
"""
struct SkeletonPair{L,R} <: GridapType
  plus::L
  minus::R
end

function Base.getproperty(x::SkeletonPair, sym::Symbol)
  if sym == :⁺
    x.plus
  elseif sym == :⁻
    x.minus
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::SkeletonPair, private=false)
  (fieldnames(typeof(x))...,:⁺,:⁻)
end

"""
    test_triangulation(trian::Triangulation)
"""
function test_triangulation(trian::Triangulation{Dc,Dp}) where {Dc,Dp}
  @test num_cell_dims(trian) == Dc
  @test num_point_dims(trian) == Dp
  @test num_cell_dims(typeof(trian)) == Dc
  @test num_point_dims(typeof(trian)) == Dp
  cell_coords = get_cell_coordinates(trian)
  @test isa(cell_coords,AbstractArray{<:AbstractVector{<:Point}})
  reffes = get_reffes(trian)
  @test isa(reffes,AbstractVector{<:LagrangianRefFE{Dc}})
  cell_types = get_cell_type(trian)
  @test isa(cell_types,AbstractArray{<:Integer})
  ncells = num_cells(trian)
  @test ncells == length(cell_coords)
  @test ncells == length(cell_types)
  @test isa(TriangulationStyle(trian),TriangulationStyle)
  bgtrian = get_background_triangulation(trian)
  @test isa(bgtrian,Triangulation)
  cell_id = get_cell_to_bgcell(trian)
  @test isa(cell_id,AbstractArray) || isa(cell_id,SkeletonPair)
  cell_ref_map = get_cell_ref_map(trian)
  @test isa(cell_ref_map,AbstractArray) || isa(cell_ref_map,SkeletonPair)
end

# Some API

"""
    num_cells(trian::Triangulation) -> Int
"""
num_cells(trian::Triangulation) = length(get_cell_type(trian))

"""
    num_nodes(trian::Triangulation) -> Int
"""
num_nodes(trian::Triangulation) = length(get_node_coordinates(trian))

"""
    num_cell_dims(::Triangulation) -> Int
    num_cell_dims(::Type{<:Triangulation}) -> Int
"""
num_cell_dims(::Triangulation{Dc,Dp}) where {Dc,Dp} = Dc
num_cell_dims(::Type{<:Triangulation{Dc,Dp}}) where {Dc,Dp} = Dc

"""
    num_point_dims(::Triangulation) -> Int
    num_point_dims(::Type{<:Triangulation}) -> Int
"""
num_point_dims(::Triangulation{Dc,Dp}) where {Dc,Dp} = Dp
num_point_dims(::Type{<:Triangulation{Dc,Dp}}) where {Dc,Dp} = Dp

"""
    num_dims(::Triangulation) -> Int
    num_dims(::Type{<:Triangulation}) -> Int

Equivalent to `num_cell_dims`.
"""
num_dims(g::Triangulation{Dc}) where Dc = Dc
num_dims(::Type{<:Triangulation{Dc}}) where Dc = Dc

"""
    is_first_order(trian::Triangulation) -> Bool
"""
function is_first_order(trian::Triangulation)
  reffes = get_reffes(trian)
  all(map(is_first_order,reffes))
end

"""
    get_cell_reffe(trian::Triangulation) -> Vector{<:LagrangianRefFE}

It is not desirable to iterate over the resulting array
for large number of cells if the underlying reference FEs
are of different Julia type.
"""
function get_cell_reffe(trian::Triangulation)
  type_to_reffe = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  expand_cell_data(type_to_reffe,cell_to_type)
end

"""
"""
function get_cell_ref_coordinates(trian::Triangulation)
  type_to_reffe = get_reffes(trian)
  type_to_coords = map(get_node_coordinates,type_to_reffe)
  cell_to_type = get_cell_type(trian)
  expand_cell_data(type_to_coords,cell_to_type)
end

"""
    get_cell_shapefuns(trian::Triangulation) -> Vector{<:Field}
"""
function get_cell_shapefuns(trian::Triangulation)
  type_to_reffes = get_reffes(trian)
  cell_to_type = get_cell_type(trian)
  type_to_shapefuns = map(get_shapefuns, type_to_reffes)
  expand_cell_data(type_to_shapefuns,cell_to_type)
end

"""
    get_cell_map(trian::Triangulation) -> Vector{<:Field}
"""
function get_cell_map(trian::Triangulation)
  cell_to_coords = get_cell_coordinates(trian)
  cell_to_shapefuns = get_cell_shapefuns(trian)
  lazy_map(linear_combination,cell_to_coords,cell_to_shapefuns)
end

function get_node_coordinates(trian::Triangulation)
  @notimplemented
end

function get_cell_node_ids(trian::Triangulation)
  @notimplemented
end

function Quadrature(trian::Triangulation,args...;kwargs...)
  cell_ctype = get_cell_type(trian)
  ctype_polytope = map(get_polytope,get_reffes(trian))
  ctype_quad = map(p->Quadrature(p,args...;kwargs...),ctype_polytope)
  cell_quad = expand_cell_data(ctype_quad,cell_ctype)
end

#"""
#"""
#function CellField(object,trian::Triangulation)
#  cm = get_cell_map(trian)
#  convert_to_cell_field(object,cm)
#end
#

# Helpers for Triangulation

#"""
#    restrict(cf::CellField,trian::Triangulation)
#"""
#function restrict(cf::CellField,trian::Triangulation)
#  _cf = to_ref_space(cf)
#  a = get_array(_cf)
#  r = restrict(a,trian)
#  axs = reindex(get_cell_axes(cf),trian)
#  _restrict_cell_field(r,axs,MetaSizeStyle(cf),trian)
#end
#
#function _restrict_cell_field(r::AbstractArray,axs::AbstractArray,msize_style::Val,trian)
#  cm = get_cell_map(trian)
#  GenericCellField(r,cm,Val(true),axs,msize_style)
#end
#
#function _restrict_cell_field(r::SkeletonPair,axs::SkeletonPair,msize_style::Val,trian)
#  cm = get_cell_map(trian)
#  la = r.plus
#  ra = r.minus
#  l = GenericCellField(la,cm,Val(true),axs.plus,msize_style)
#  r = GenericCellField(ra,cm,Val(true),axs.minus,msize_style)
#  merge_cell_fields_at_skeleton(l,r)
#end
#
#"""
#    CellQuadrature(trian::Triangulation, degree::Integer)
#"""
#function CellQuadrature(trian::Triangulation, degree::Integer)
#  polytopes = map(get_polytope,get_reffes(trian))
#  cell_type = get_cell_type(trian)
#  CellQuadrature(degree,polytopes,cell_type)
#end
#
#"""
#    integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
#
#The `cell_field` is aligned with the cells in `trian`
#"""
#function integrate(cell_field,trian::Triangulation,quad::CellQuadrature)
#  cell_map = get_cell_map(trian)
#  integrate(cell_field,cell_map,quad)
#end
#
#function CellField(value::Number,trian::Triangulation,quad::CellQuadrature)
#  CellField(value,get_cell_map(trian),quad)
#end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

#

#

# In a bifurcation, there is a `parent` vessel (indicated as 1 in the figure)
# and two `daughter` vessels (2 and 3).
#
# The bifurcation is solved by imposing the conservation of mass and of
# __*static*__ pressure at the bifurcation node $b$. Three additional
# relations are obtained by extrapolating the three outgoing characteristics.
# The solution process is the same described for the
# [conjunction](conjuncitons.html) case.
#
# $U$ and $F$ vectors read
# $$
#   U = \{U_i\} = \left\{ \begin{array}{c}
#                           u_1 \\
#                           u_2 \\
#                           u_3 \\
#                           A_1^{1/4} \\
#                           A_2^{1/4} \\
#                           A_3^{1/4}
#                        \end{array} \right\} , \quad i =1, ..., 6,
# $$
# $$
#   F = \left\{f_i \right\} =
#   \begin{cases}
#     U_1 + 4k_1U_4 - W_1^* = 0, \\
#     U_2 - 4k_2U_5 - W_2^* = 0, \\
#     U_3 - 4k_3U_6 - W_3^* = 0, \\
#     U_1U_4^4 - U_2U_5^4 - U_3U_6^4 = 0, \\
#     \beta_1 \left(\tfrac{U_4^2}{A_{01}^{1/2}} -1 \right) -
#     \beta_2 \left(\tfrac{U_5^2}{A_{02}^{1/2}} -1 \right) = 0, \\
#     \beta_1 \left(\tfrac{U_4^2}{A_{01}^{1/2}} -1 \right) -
#     \beta_3 \left(\tfrac{U_6^2}{A_{03}^{1/2}} -1 \right) = 0, \\
#   \end{cases}
# $$
# and the Jacobian reads
# $$
#   J = \left[ \begin{array}{cccccc}
#         1 & 0 & 0 & 4k_1 & 0 & 0 \\
#         0 & 1 & 0 & 0 & -4k_2 & 0 \\
#         0 & 0 & 1 & 0 & 0 & -4k_3 \\
#         U_4^4 & -U_5^4 & -U_6^4 & 4U_1 U_4^3 & -4 U_2 U_5^3 & -4U_3 U_6^3 \\
#         0 & 0 & 0 & 2\beta_1 U_4/A_{01}^{1/2} &
#                    -2\beta_2 U_5/A_{02}^{1/2} & 0 \\
#         0 & 0 & 0 & 2\beta_1 U_4/A_{01}^{1/2} &
#                    -2\beta_2 U_5/A_{02}^{1/2} & 0 \\
#       \end{array} \right].
# $$
# The solution is performed by
# [`solveBifurcation`](bifurcations.html#solveBifurcation) and follows the
# same procedure of [`solveConjunction`](conjunctions.html#solveConjuntion)
# function. Refer to [conjunctions.jl](conjunctions.html) documentation for
# further informations.

# *function* __`solveBifurcation`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------------- --------------------------------------------------------
# `v1`                `::Vessel` data structure for the parent vessel.
#
# `v2`                `::Vessel` data structure for the first daughter vessel.
#
# `v3`                `::Vessel` data structure for the second daughter vessel.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# See [`solveConjunction`](conjunctions.html#solveConjunction) documentation
# for details.
# ----------------------------------------------------------------------------
#
function solveBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)

  #Unknowns vector
  U = [v1.u[end],
       v2.u[ 1 ],
       v3.u[ 1 ],
       sqrt(sqrt(v1.A[end])),
       sqrt(sqrt(v2.A[ 1 ])),
       sqrt(sqrt(v3.A[ 1 ]))]

  #Parameters vector
  k1 = sqrt(0.5*3*v1.gamma[end])
  k2 = sqrt(0.5*3*v2.gamma[ 1 ])
  k3 = sqrt(0.5*3*v3.gamma[ 1 ])
  k = [k1, k2, k3]

  W = calculateWstarBif(U, k)
  J = calculateJacobianBif(v1, v2, v3, U, k)
  F = calculateFofUBif(v1, v2, v3, U, k, W)

  #Newton-Raphson
  nr_toll_U = 1e-5
  nr_toll_F = 1e-5

  while true
    dU = J\(-F)
    # U_new = U + 0.01*dU
    U_new = U + dU

    if any(isnan(dot(F,F)))
      println(F)
      @printf "error at bifurcation with vessels %s, %s, and %s \n" v1.label v2.label v3.label
      break
    end

    u_ok = 0
    f_ok = 0
    for i in 1:length(dU)
      if abs(dU[i]) <= nr_toll_U || abs(F[i]) <= nr_toll_F
        u_ok += 1
        f_ok += 1
      end
    end

    if u_ok == length(dU) || f_ok == length(dU)
      U = U_new
      break
    else
      U = U_new
      W = calculateWstarBif(U, k)
      F = calculateFofUBif(v1, v2, v3, U, k, W)
    end
  end

  #Update vessel quantities
  v1.u[end] = U[1]
  v2.u[ 1 ] = U[2]
  v3.u[ 1 ] = U[3]

  v1.A[end] = U[4]*U[4]*U[4]*U[4]
  v2.A[ 1 ] = U[5]*U[5]*U[5]*U[5]
  v3.A[ 1 ] = U[6]*U[6]*U[6]*U[6]

  v1.Q[end] = v1.u[end]*v1.A[end]
  v2.Q[ 1 ] = v2.u[ 1 ]*v2.A[ 1 ]
  v3.Q[ 1 ] = v3.u[ 1 ]*v3.A[ 1 ]

  v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
  v2.P[ 1 ] = pressure(v2.A[ 1 ], v2.A0[ 1 ], v2.beta[ 1 ], v2.Pext)
  v3.P[ 1 ] = pressure(v3.A[ 1 ], v3.A0[ 1 ], v3.beta[ 1 ], v3.Pext)

  v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
  v2.c[ 1 ] = waveSpeed(v2.A[ 1 ], v2.gamma[ 1 ])
  v3.c[ 1 ] = waveSpeed(v3.A[ 1 ], v3.gamma[ 1 ])

end

# *function* __`calculateWstarBif`__ $\rightarrow$
# `W::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Riemann invariants are computed starting from unknown and parameters
# vectors as
# $$
#   W_1 = u_1 + 4c_1, \quad W_2 = u_2 - 4c_2, \quad W_3 = u_3 - 4c_3.
# $$
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `W`         `::Array` containing outgoing characteristics from the three
#             vessels.
# ----------------------------------------------------------------------------
#
function calculateWstarBif(U :: Array, k :: Array)

  W1 = U[1] + 4*k[1]*U[4]
  W2 = U[2] - 4*k[2]*U[5]
  W3 = U[3] - 4*k[3]*U[6]

  return [W1, W2, W3]
end

# *function* __`calculateFofUBif`__ $\rightarrow$ `F::Array`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `b`            `::Blood` data structure.
#
# `v1`           `::Vessel` parent vessel data structure.
#
# `v2`           `::Vessel` first daughter vessel data structure.
#
# `v3`           `::Vessel` second daughter vessel data structure.
#
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
#
# `W`            `::Array` outgoing characteristics array.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# $F$ array is computed by imposing the conservation of mass and static
# pressure at bifurcation node. $F$ reads
# $$
#   F = \left\{f_i \right\} =
#   \begin{cases}
#     U_1 + 4k_1U_4 - W_1^* = 0, \\
#     U_2 - 4k_2U_5 - W_2^* = 0, \\
#     U_3 - 4k_3U_6 - W_3^* = 0, \\
#     U_1U_4^4 - U_2U_5^4 - U_3U_6^4 = 0, \\
#     \beta_1 \left(\tfrac{U_4^2}{A_{01}^{1/2}} -1 \right) -
#     \beta_2 \left(\tfrac{U_5^2}{A_{02}^{1/2}} -1 \right) = 0, \\
#     \beta_1 \left(\tfrac{U_4^2}{A_{01}^{1/2}} -1 \right) -
#     \beta_3 \left(\tfrac{U_6^2}{A_{03}^{1/2}} -1 \right) = 0, \\
#   \end{cases}
# $$
# ----------------------------------------------------------------------------
# Total pressure conservation can be imposed by replacing `f5` and `f6` with
#
#     f5 = v1.beta*( U[4]^2 - sqrt(v1.A0) ) -
#        ( v2.beta*( U[5]^2 - sqrt(v2.A0) )
#
#     f6 = v1.beta*( U[4]^2 - sqrt(v1.A0) ) -
#        ( v3.beta*( U[6]^2 - sqrt(v3.A0) )
#
# respectively.
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `F`         `::Array` Newton's method relations.
# ----------------------------------------------------------------------------
#
function calculateFofUBif(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel,
                           U :: Array,   k :: Array,   W :: Array)

  f1 = U[1] + 4*k[1]*U[4] - W[1]

  f2 = U[2] - 4*k[2]*U[5] - W[2]

  f3 = U[3] - 4*k[3]*U[6] - W[3]

  f4 = U[1]*(U[4]*U[4]*U[4]*U[4]) - U[2]*(U[5]*U[5]*U[5]*U[5]) - U[3]*(U[6]*U[6]*U[6]*U[6])

  f5 = v1.beta[end]*(U[4]*U[4]/sqrt(v1.A0[end])  - 1  ) -
     ( v2.beta[ 1 ]*(U[5]*U[5]/sqrt(v2.A0[1]) - 1) )

  f6 = v1.beta[end]*(U[4]*U[4]/sqrt(v1.A0[end])  - 1  ) -
     ( v3.beta[ 1 ]*(U[6]*U[6]/sqrt(v3.A0[1]) - 1) )

  # f5 = v1.beta[end]*(U[4]*U[4]/sqrt(v1.A0[end])  - 1  ) -
  #    ( v2.beta[ 1 ]*(U[5]*U[5]/sqrt(v2.A0[1]) - 1) ) + 1060*0.5*(U[1]*U[1]-U[2]*U[2])
  #
  # f6 = v1.beta[end]*(U[4]*U[4]/sqrt(v1.A0[end])  - 1  ) -
  #    ( v3.beta[ 1 ]*(U[6]*U[6]/sqrt(v3.A0[1]) - 1) ) + 1060*0.5*(U[1]*U[1]-U[3]*U[3])

  return [f1, f2, f3, f4, f5, f6]
end

# *function* __`calculateJacobianBif`__ $\rightarrow$
# `W::Array{Float, 1}`
#
# ----------------------------------------------------------------------------
# Parameters:
# -------------- -------------------------------------------------------------
# `b`            `::Blood` data structure.
#
# `v1`           `::Vessel` parent vessel data structure.
#
# `v2`           `::Vessel` daughter vessel data structure.
#
# `U`            `::Array` junction unknown vector.
#
# `k`            `::Array` junction parameters vector.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The Jacobian is computed as
# $$
#   J = \left[ \begin{array}{cccccc}
#         1 & 0 & 0 & 4k_1 & 0 & 0 \\
#         0 & 1 & 0 & 0 & -4k_2 & 0 \\
#         0 & 0 & 1 & 0 & 0 & -4k_3 \\
#         U_4^4 & -U_5^4 & -U_6^4 & 4U_1 U_4^3 & -4 U_2 U_5^3 & -4U_3 U_6^3 \\
#         0 & 0 & 0 & 2\beta_1 U_4/A_{01}^{1/2} &
#                    -2\beta_2 U_5/A_{02}^{1/2} & 0 \\
#         0 & 0 & 0 & 2\beta_1 U_4/A_{01}^{1/2} &
#                    -2\beta_2 U_5/A_{02}^{1/2} & 0 \\
#       \end{array} \right].
# $$
# ----------------------------------------------------------------------------
#
# The conservation of total pressure is imposed by setting the following
# elements different than zero:
#
#     J[5,1] =  rho*U[1]
#     J[5,2] = -rho*U[2]
#     J[5,4] =  2*v1.beta*U[4]
#     J[5,5] = -2*v2.beta*U[5]
#
#     J[6,1] =  rho*U[1]
#     J[6,3] = -rho*U[3]
#     J[6,4] =  2*v1.beta*U[4]
#     J[6,6] = -2*v3.beta*U[6]
#
# ----------------------------------------------------------------------------
# Returns:
# ----------- ----------------------------------------------------------------
# `J`         `::Array{Float, 2}` Jacobian matrix.
# ----------------------------------------------------------------------------
#
function calculateJacobianBif(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel,
                               U :: Array,   k :: Array)

  J = eye(6)

  J[1,4] =  4*k[1]
  J[2,5] = -4*k[2]
  J[3,6] = -4*k[3]

  J[4,1] =  (U[4]*U[4]*U[4]*U[4])
  J[4,2] = -(U[5]*U[5]*U[5]*U[5])
  J[4,3] = -(U[6]*U[6]*U[6]*U[6])
  J[4,4] =  4*U[1]*(U[4]*U[4]*U[4])
  J[4,5] = -4*U[2]*(U[5]*U[5]*U[5])
  J[4,6] = -4*U[3]*(U[6]*U[6]*U[6])

  J[5,1] =  0.
  # J[5,1] = 1060*U[1]

  J[5,2] =  0.
  # J[5,2] = -1060*U[2]

  J[5,4] =  2*v1.beta[end]*U[4]/sqrt(v1.A0[end])
  J[5,5] = -2*v2.beta[ 1 ]*U[5]/sqrt(v2.A0[1])

  J[6,1] =  0.
  # J[6,1] = 1060*U[1]

  J[6,3] =  0.
  # J[6,3] = -1060U[3]

  J[6,4] =  2*v1.beta[end]*U[4]/sqrt(v1.A0[end])
  J[6,6] = -2*v3.beta[ 1 ]*U[6]/sqrt(v3.A0[1])

  return J
end


text:
@testset "Optimized 1Tree" begin
    module_path = dirname(pathof(TSPSolver))
    g = TSPSolver.simple_parse_tsp(joinpath(module_path, "../test/data/bier127.tsp"))
    root = TSPSolver.Root(g)

    tree, lb = TSPSolver.get_optimized_1tree(root.cost)
    @show lb
    @test lb < 118293.52381566973

    TSPSolver.fix_edge!(root, (2, 3))
    TSPSolver.fix_edge!(root, (2, 5))
    tree, lb = TSPSolver.get_optimized_1tree(root.cost)
    @test Edge(2,3) in tree || Edge(3,2) in tree
    @test Edge(2,5) in tree || Edge(5,2) in tree
end

@testset "Greedy" begin
    module_path = dirname(pathof(TSPSolver))
    g = TSPSolver.simple_parse_tsp(joinpath(module_path, "../test/data/bier127.tsp"))
    root = TSPSolver.Root(g)
    original_root = deepcopy(root)

    @time tour, ub = TSPSolver.greedy(root.g)
    @show ub
    @test length(tour) == nv(g)
    @test ub > 118293.52381566973

    # fix the first edges as it was in greedy anyway
    v1, v2 = tour[1], tour[2]
    extra_cost = TSPSolver.fix_edge!(root, (v1, v2))
    new_tour, new_ub = TSPSolver.greedy(root.g)
    @test new_ub + extra_cost ≈ ub
    @test new_tour == tour

    # fix a couple of edges and check that they were fixed correctly
    root = deepcopy(original_root)
    fixed_edges = zeros(Int, nv(g))
    fixed_edges[1] = 2
    fixed_edges[2] = 3
    fixed_edges[4] = 8
    extra_cost = 0.0
    extra_cost += TSPSolver.fix_edge!(root, (1,2))
    extra_cost += TSPSolver.fix_edge!(root, (2,3))
    extra_cost += TSPSolver.fix_edge!(root, (4,8))
    tour, ub = TSPSolver.greedy(root.g)
    ub += extra_cost
    c = 0
    for i in 1:nv(g)
        if fixed_edges[tour[i]] != 0
            c += 1
            @test tour[i+1] == fixed_edges[tour[i]]
        end
    end
    @test c == count(i->i!=0, fixed_edges)

    # fix a couple of edges and check that they were fixed correctly
    # also disallow some edges
    root = deepcopy(original_root)
    fixed_edges = zeros(Int, nv(root.g))
    fixed_edges[1] = 2
    fixed_edges[2] = 3
    fixed_edges[4] = 8
    extra_cost = 0.0
    extra_cost += TSPSolver.fix_edge!(root, (1,2))
    extra_cost += TSPSolver.fix_edge!(root, (2,3))
    extra_cost += TSPSolver.fix_edge!(root, (4,8))

    disallow_edges = Dict{Int, Set{Int}}()
    N = nv(g)
    disallow_edges[7] = Set([i for i in 1:N if !(i in [2,3,8,9])])
    for i in 1:N
        i in [2,3,8,9] && continue
        TSPSolver.disallow_edge!(root, (7, i))
    end

    tour, ub = TSPSolver.greedy(root.g)
    ub += extra_cost
    if tour === nothing
        @test isnan(ub)
    else
        c = 0
        for i in 1:N
            if fixed_edges[tour[i]] != 0
                c += 1
                @test tour[i+1] == fixed_edges[tour[i]]
            end

            if tour[i] == 7
                if i + 1 <= N
                    @test tour[i+1] == 9 # only 9 is not
                else
                    @test tour[1] == 9
                end
            end
        end
        @test c == count(i->i!=0, fixed_edges)
    end
end


text:
# QuEST_jl/src/base/state_init.jl
#

function cloneQureg(targetQureg ::QuEST_Types.Qureg, copyQureg ::QuEST_Types.Qureg) ::Nothing

    return ccall(:cloneQureg, Cvoid, (QuEST_Types.Qureg, QuEST_Types.Qureg), targetQureg, copyQureg)

end

function initBlankState(qureg ::QuEST_Types.Qureg) ::Nothing

    ccall(:initBlankState, Cvoid,
          (QuEST_Types.Qureg,),
          qureg)
end

function initClassicalState(qureg    ::QuEST_Types.Qureg,
                            stateIdx ::Integer)           ::Nothing

    ccall(:initClassicalState, Cvoid,
          (QuEST_Types.Qureg, Clonglong),
          qureg,              stateIdx)
end

function initPlusState(qureg ::QuEST_Types.Qureg) ::Nothing

    ccall(:initPlusState, Cvoid,
          (QuEST_Types.Qureg,),
          qureg)
end

function initPureState(qureg ::QuEST_Types.Qureg,
                       pure  ::QuEST_Types.Qureg) ::Nothing

    @assert qureg.numQubitsRepresented == pure.numQubitsRepresented
    @assert pure.isDensityMatrix == 0

    ccall(:initPureState, Cvoid,
          (QuEST_Types.Qureg,  QuEST_Types.Qureg),
          qureg,               pure)
end

function initStateFromAmps(qureg     ::QuEST_Types.Qureg,
                           amps_real ::Vector{Qreal},
                           amps_imag ::Vector{Qreal}      ) ::Nothing

    @assert length(amps_real) == qureg.numAmpsTotal
    @assert length(amps_imag) == qureg.numAmpsTotal

    ccall(:initStateFromAmps,
          Cvoid,
          (QuEST_Types.Qureg, Ptr{Qreal}, Ptr{Qreal}),
          qureg,              amps_real,  amps_imag)
end

function initZeroState(qureg ::QuEST_Types.Qureg) ::Nothing

    ccall(:initZeroState, Cvoid,
          (QuEST_Types.Qureg,),
          qureg)
end

function setAmps(qureg        ::QuEST_Types.Qureg,
                 startIdx     ::Integer,
                 amps_real    ::Vector{Qreal},
                 amps_imag    ::Vector{Qreal},
                 numAmps      ::Clonglong)   :: Nothing

    @assert numAmps             == length(amps_real)
    @assert numAmps             == length(amps_imag)
    @assert startIdx + numAmps  <= qureg.numAmpsTotal

    ccall(:setAmps, Cvoid,
          (QuEST_Types.Qureg, Clonglong, Ptr{Qreal}, Ptr{Qreal}, Clonglong),
          qureg,              startIdx,  amps_real,  amps_imag,  numAmps)
end

const __ℂ     = QuEST_Types.Complex
const __Qureg = QuEST_Types.Qureg

function setWeightedQureg(factor1    ::Complex{Qreal},
                          qureg1     ::QuEST_Types.Qureg,
                          factor2    ::Complex{Qreal},
                          qureg2     ::QuEST_Types.Qureg,
                          factorOut  ::Complex{Qreal},
                          quregOut   ::QuEST_Types.Qureg)  ::Nothing

    @assert Bool(qureg1.isDensityMatrix) == Bool(qureg2.isDensityMatrix) == Bool(quregOut.isDensityMatrix)
    @assert qureg1.numQubitsRepresented  == qureg2.numQubitsRepresented  == quregOut.numQubitsRepresented

    local factor1_quest   = __ℂ(real(factor1  ),imag(factor1  ))
    local factor2_quest   = __ℂ(real(factor2  ),imag(factor2  ))
    local factorOut_quest = __ℂ(real(factorOut),imag(factorOut))

    ccall(:setWeightedQureg, Cvoid,
          (__ℂ,          __Qureg,  __ℂ,           __Qureg, __ℂ,             __Qureg),
          factor1_quest, qureg1,   factor2_quest, qureg2,  factorOut_quest, quregOut)
end


text:
using Unicode

# Parser
# ======

# Missing
# -------

function fillcolumn!(col::Vector{Missing}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    for i in 1:nvals
        @inbounds col[end-nvals+i] = missing
    end
    return col
end

# Integer
# -------

function fillcolumn!(col::Vector{Int}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_integer(mem, start, length)
    end
    return col
end

function fillcolumn!(col::Vector{Union{Int,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(t)
            col[end-nvals+i] = parse_integer(mem, start, length)
        end
    end
    return col
end

const SAFE_INT_LENGTH = sizeof(string(typemax(Int))) - 1

@inline function parse_integer(mem::Memory, start::Int, length::Int)
    stop = start + length - 1
    if length > SAFE_INT_LENGTH
        # safe but slow fallback
        buf = IOBuffer()
        for i in start:stop
            write(buf, mem[i])
        end
        return parse(Int, String(take!(buf)))
    end
    i = start
    @inbounds b = mem[i]
    negative = false
    if b == UInt8('-')
        negative = true
        i += 1
    elseif b == UInt8('+')
        i += 1
    end
    n::Int = 0
    while i ≤ stop
        @inbounds b = mem[i]
        n = 10n + (b - UInt8('0'))
        i += 1
    end
    return negative ? -n : n
end

# Float
# -----

function fillcolumn!(col::Vector{Float64}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    buf = Vector{UInt8}(undef, 32)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_float!(buf, mem, start, length)
    end
    return col
end

function fillcolumn!(col::Vector{Union{Float64,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    buf = Vector{UInt8}(undef, 32)
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(tokens[c,i])
            col[end-nvals+i] = parse_float!(buf, mem, start, length)
        end
    end
    return col
end

@inline function parse_float!(buf::Vector{UInt8}, mem::Memory, start::Int, length::Int)
    if Base.length(buf) < length + 1
        resize!(buf, length + 1)
    end
    ccall(:memmove, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), pointer(buf), mem.ptr + start - 1, length)
    @inbounds buf[length+1] = 0x00  # terminate with NUL
    # TODO: The next line is the bottleneck of parsing lots of floats!
    hasvalue, value = ccall(:jl_try_substrtod, Tuple{Bool,Float64}, (Ptr{UInt8}, Csize_t, Csize_t), pointer(buf), 0, length)
    if !hasvalue
        throw(ReadError("failed to parse a floating-point number"))
    end
    return value
end

# Bool
# ----

function fillcolumn!(col::Vector{Bool}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_bool(mem, start, length)
    end
    return col
end

@inline function parse_bool(mem::Memory, start::Int, length::Int)
    c = mem[start]
    # No need to check all the bytes as the format is already validated.
    return (c == UInt8('f') || c == UInt8('F')) ? false : true
end

# String
# ------

const N_CACHED_STRINGS = 8

function fillcolumn!(col::Vector{String}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    cache = StringCache(N_CACHED_STRINGS)
    usecache = true
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        start, length = location(t)
        if kind(t) & QSTRING != 0
            s = qstring(mem, start, length, quot)
        elseif usecache
            s = allocate!(cache, mem.ptr + start - 1, length)
        else
            s = unsafe_string(mem.ptr + start - 1, length)
        end
        col[end-nvals+i] = s
        if usecache && i % 512 == 0 && 10 * cache.stats.hit < cache.stats.hit + cache.stats.miss
            # stop using cache because the cache hit rate is too low
            usecache = false
        end
    end
    if !usecache
        @debug "Cache is turned off"
    end
    @debug cache
    return col
end

function fillcolumn!(col::Vector{Union{String,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    cache = StringCache(N_CACHED_STRINGS)
    usecache = true
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(t)
            if kind(t) & QSTRING != 0
                s = qstring(mem, start, length, quot)
            elseif usecache
                s = allocate!(cache, mem.ptr + start - 1, length)
            else
                s = unsafe_string(mem.ptr + start - 1, length)
            end
            col[end-nvals+i] = s
            if usecache && i % 512 == 0 && 10 * cache.stats.hit < cache.stats.hit + cache.stats.miss
                # stop using cache because the cache hit rate is too low
                usecache = false
            end
        end
    end
    if !usecache
        @debug "Cache is turned off"
    end
    @debug cache
    return col
end

function qstring(mem::Memory, start::Int, length::Int, quot::UInt8)
    i = start
    skip = false
    buf = IOBuffer()
    while i < start + length
        if skip
            skip = false
        else
            x = mem[i]
            write(buf, x)
            skip = x == quot
        end
        i += 1
    end
    return String(take!(buf))
end

# Date
# ----

const DATE_REGEX = r"^\d{4}-\d{2}-\d{2}$"

function parse_date(s::String)
    return Date(s, dateformat"y-m-d")
end

function parse_date(col::Vector{String})
    out = Vector{Date}(undef, length(col))
    for i in 1:length(col)
        out[i] = parse_date(col[i])
    end
    return out
end

function parse_date(col::Vector{Union{String,Missing}})
    out = Vector{Union{String,Missing}}(undef, length(col))
    for i in 1:length(col)
        x = col[i]
        if x === missing
            out[i] = missing
        else
            out[i] = parse_date(x)
        end
    end
    return out
end

function is_date_like(col::Vector{<:Union{String,Missing}})
    i = 1
    n = 0
    while i ≤ length(col) && n < 3
        x = col[i]
        if x isa String
            if !occursin(DATE_REGEX, x)
                return false
            end
            n += 1
        end
        i += 1
    end
    return n > 0
end

# DateTime
# --------

const DATETIME_REGEX = r"^\d{4}-\d{2}-\d{2}[T ]\d{2}:\d{2}:\d{2}(:?\.\d+)?$"

function parse_datetime(s::String, hasT::Bool)
    if hasT
        return DateTime(s, dateformat"y-m-dTH:M:S.s")
    else
        return DateTime(s, dateformat"y-m-d H:M:S.s")
    end
end

function is_T_delimited(col::Vector{<:Union{String,Missing}})
    for x in col
        if !Base.ismissing(x)
            return occursin('T', x)
        end
    end
    return false
end

function is_datetime_like(col::Vector{<:Union{String,Missing}})
    # Check if the first three strings (if any) are datetime-like.
    i = 1
    n = 0
    while i ≤ length(col) && n < 3
        x = col[i]
        if x isa String
            if !occursin(DATETIME_REGEX, x)
                return false
            end
            n += 1
        end
        i += 1
    end
    return n > 0
end

function parse_datetime(col::Vector{String}, hasT::Bool)
    out = Vector{DateTime}(undef, length(col))
    for i in 1:length(col)
        out[i] = parse_datetime(col[i], hasT)
    end
    return out
end

function parse_datetime(col::Vector{Union{String,Missing}}, hasT::Bool)
    out = Vector{Union{DateTime,Missing}}(undef, length(col))
    for i in 1:length(col)
        x = col[i]
        if x === missing
            out[i] = missing
        else
            out[i] = parse_datetime(x, hasT)
        end
    end
    return out
end

# these are Julia keywords used by normalizename to prevent special keywords
# from becoming column names
const RESERVED = Set(["local", "global", "export", "let",
    "for", "struct", "while", "const", "continue", "import",
    "function", "if", "else", "try", "begin", "break", "catch",
    "return", "using", "baremodule", "macro", "finally",
    "module", "elseif", "end", "quote", "do"])

normalizename(name::Symbol) = name

# normalize name tries to make strings safe for use as Julia symbols and
# dataframe columns
function normalizename(name::String)
    uname = strip(Unicode.normalize(name))
    id = Base.isidentifier(uname) ? uname : map(c->Base.is_id_char(c) ? c : '_', uname)
    cleansed = string((isempty(id) || !Base.is_id_start_char(id[1]) || id in RESERVED) ? "_" : "", id)
    return Symbol(replace(cleansed, r"(_)\1+"=>"_"))
end


text:
seats = []
open("day05/day05.in") do io
    global map_str, passports
    for line in readlines(io)
        push!(seats, line)
    end
end

function seatID(s::String)::Int64
    id = 0
    for chi in 1:10
        id *= 2
        if s[chi] == 'B' || s[chi] == 'R'
            id += 1
        end
    end
    return id
end

seatID_max = 0
for seat_s in seats
    global seatID_max
    seatID_max = max(seatID_max, seatID(seat_s))
end
println(seatID_max)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code