text:
export extend

mutable struct NfToNfMor <: Map{AnticNumberField, AnticNumberField, HeckeMap, NfToNfMor}
  header::MapHeader{AnticNumberField, AnticNumberField}
  prim_img::nf_elem

  function NfToNfMor()
    z = new()
    z.header = MapHeader()
    return r
  end

  function NfToNfMor(K::AnticNumberField, L::AnticNumberField, y::nf_elem)
    z = new()
    z.prim_img = y

    function image(x::nf_elem)
      g = parent(K.pol)(x)
      return evaluate(g, y)
    end

    z.header = MapHeader(K, L, image)
    return z
  end
end

function Base.:(==)(f::NfToNfMor, g::NfToNfMor)
  if (domain(f) != domain(g)) || (codomain(f) != codomain(g))
    return false
  end

  return f.prim_img == g.prim_img
end

function *(f::NfToNfMor, g::NfToNfMor)
  domain(f) == codomain(g) || throw("Maps not compatible")

  a = gen(domain(g))
  y = f(g(a))

  return NfToNfMor(domain(g), codomain(f), y)
end

function show(io::IO, h::NfToNfMor)
  if domain(h) == codomain(h)
    println(io, "Automorphism of ", domain(h))
  else
    println(io, "Injection of ", domain(h), " into ", codomain(h))
  end
  println(io, "defined by ", gen(domain(h)), " -> ", h.prim_img)
end

function automorphisms(K::AnticNumberField)
  f = K.pol
  lr = roots(f, K)
  Aut = Hecke.NfToNfMor[]
  for r in lr
    push!(Aut, Hecke.NfToNfMor(K, K, r))
  end
  return Aut
end

hom(K::AnticNumberField, L::AnticNumberField, a::nf_elem) = NfToNfMor(K, L, a)


text:
import IntervalArithmetic: radius, mid, mag, mig, sup, inf
export radius, mid, mag, mig, sup, inf

for f in (:radius, :mid, :mag, :mig, :sup, :inf)
    @eval $f(a::NumberInterval) = $f(Interval(a))
end

ext: jl


text:
function AdamsIntegrator(arg0::JString, arg1::jint, arg2::jint, arg3::jdouble, arg4::jdouble, arg5::Vector{jdouble}, arg6::Vector{jdouble})
    return AdamsIntegrator((JString, jint, jint, jdouble, jdouble, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function AdamsIntegrator(arg0::JString, arg1::jint, arg2::jint, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble)
    return AdamsIntegrator((JString, jint, jint, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function integrate(obj::AdamsIntegrator, arg0::ExpandableODE, arg1::ODEState, arg2::jdouble)
    return jcall(obj, "integrate", ODEStateAndDerivative, (ExpandableODE, ODEState, jdouble), arg0, arg1, arg2)
end


text:
function run(s)
    # Your code here
    entries = split(s, "\n")
    entries = map(x -> parse(Int, x), entries)
    for el1 in entries
        for el2 in entries
            if el1 + el2 == 2020
                return el1 * el2
            end
        end
    end
    return 0
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
Base.ndims(::AbstractDomain) = error("Domain dimension is undefined")

struct DiscreteDomain <: AbstractDomain
    length::Int
end

struct Reals <: AbstractDomain end

struct Interval{T} <: AbstractDomain
  start::T
  stop::T
end

struct ProductDomain{T} <: AbstractDomain
    domains::T
end

Base.:*(x::AbstractDomain...) = ProductDomain((x...))

export Reals, DiscreteDomain, Interval, ProductDomain


text:
using DelaySSAToolkit
using Test, DiffEqJump, StaticArrays
#σ_off: Gon -> Goff
#σ_on: Goff -> Gon
#ρ_on: Gon -> Gon + N, triggers N => τ 0
#d: P -> 0
# 1. Gon, 2. Goff, 3. P
params = [0.1, 0.1, 10.0, 0.1, 10.0, 50.0]
σ_off, σ_on, ρ_on, d, τ, tf = params
rates = [σ_off, σ_on, ρ_on, d]
react_stoich = [[1 => 1], [2 => 1], [1 => 1], [3 => 1]]
net_stoich = [[1 => -1, 2 => 1], [1 => 1, 2 => -1], [3 => 1], [3 => -1]]
mass_action_jump = MassActionJump(rates, react_stoich, net_stoich; scale_rates=false)
jumpset = JumpSet((), (), nothing, mass_action_jump)
delay_trigger = Dict(3 => [1 => τ])
delay_complete = Dict(1 => [3 => -1])
delay_affect! = function (integrator, rng)
    i = rand(rng, 1:length(integrator.de_chan[1]))
    deleteat!(integrator.de_chan[1], i)
end
delay_interrupt = Dict(4 => delay_affect!)

delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
u0 = @SVector [0, 1, 0]
de_chan0 = [[]]
tspan = (0.0, tf)
dprob = DiscreteProblem(u0, tspan)
algs = [DelayRejection(), DelayDirect(), DelayMNRM(), DelayDirectCR()]
for alg in algs
    djprob = DelayJumpProblem(dprob, DelayRejection(), jumpset, delayjumpset, de_chan0, save_positions=(true, true), save_delay_channel=true)
    @info "Testing method $(alg)"
    sol = solve(djprob, SSAStepper())

    for i in eachindex(sol.u)
        @test sol.u[i][3] == length(sol.channel[i][1])
    end
end


text:
import DataFrames

using Dance.Router

function df(headers::Dict{String, String}) :: DataFrames.DataFrame
    return DataFrames.DataFrame(A = 1:4, B = ["M", "F", "F", "M"])
end

function dict_1(params_dict::Dict{Symbol, Union{Float64, Int64, String}}, headers::Dict{String, String}) :: Dict{Symbol, Float64}
    return Dict(:a => params_dict[:value])
end

function dict_2(params_dict::Dict{Symbol, Union{Float64, Int64, String}}, headers::Dict{String, String}) :: Dict{Symbol, Int64}
    return Dict(Symbol(params_dict[:key]) => params_dict[:value])
end

function hello(headers::Dict{String, String}) :: String
    return "Hello World"
end

route("/", hello; method=GET, endpoint=EP_HTML)
route_group(route_prefix="/dict", method=GET, endpoint=EP_HTML, [
    (path=r"/(?<value>\d*\.\d*)", action=dict_1)
    (path=r"/(?<key>\w+)/(?<value>\d{3})", action=dict_2)
])
route("/dataframe", df; method=GET, endpoint=EP_HTML)


text:
function initialize!(integrator, cache::ROCK2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  cache.max_stage = (integrator.alg.max_stages < 1 || integrator.alg.max_stages > 200) ? 200 : integrator.alg.max_stages
  cache.min_stage = (integrator.alg.min_stages > cache.max_stage) ? cache.max_stage : integrator.alg.min_stages
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ROCK2ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, fp1, fp2, recf = cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + abs(dt)*integrator.eigen_est)/0.811) + 1))
  mdeg = min(max(mdeg,cache.min_stage), cache.max_stage)
  cache.mdeg = max(mdeg, 3) - 2
  choosedeg!(cache)
  # recurrence
  # for the first stage
  tᵢ₋₁ = t + dt*recf[cache.start]
  tᵢ₋₂ = t + dt*recf[cache.start]
  tᵢ₋₃ = t
  uᵢ₋₂ = copy(uprev)
  uᵢ₋₁ = uprev + (dt*recf[cache.start])*fsalfirst
  cache.mdeg < 2 && ( u = uᵢ₋₁ )
  # for the second to the ms[cache.mdeg] th stages
  for i in 2:cache.mdeg
    μ, κ = recf[cache.start + (i - 2)*2 + 1], recf[cache.start + (i - 2)*2 + 2]
    ν = -1 - κ
    u = f(uᵢ₋₁, p, tᵢ₋₁)
    tᵢ₋₁ = dt*μ - ν*tᵢ₋₂ - κ*tᵢ₋₃
    u = (dt*μ)*u - ν*uᵢ₋₁ - κ*uᵢ₋₂
    i < cache.mdeg && (uᵢ₋₂ = uᵢ₋₁; uᵢ₋₁ = u)
    tᵢ₋₃ = tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
  end # end if
  # two-stage finishing procedure.
  δt₁ = dt*fp1[cache.deg_index]
  δt₂ = dt*fp2[cache.deg_index]
  uᵢ₋₂  = f(u, p, tᵢ₋₁)
  integrator.destats.nf += 1
  uᵢ₋₁ = u + δt₁*uᵢ₋₂
  tᵢ₋₁ += δt₁
  u = f(uᵢ₋₁, p, tᵢ₋₁)
  integrator.destats.nf += 1

  if integrator.opts.adaptive
    tmp = δt₂*(u - uᵢ₋₂)
    u = uᵢ₋₁ + δt₁*u + tmp
  else
    u = uᵢ₋₁ + δt₁*u + δt₂*(u - uᵢ₋₂)
  end
  # error estimate
  if integrator.opts.adaptive
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.u = u
end

function initialize!(integrator, cache::ROCK2Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  cache.constantcache.max_stage = (integrator.alg.max_stages < 1 || integrator.alg.max_stages > 200) ? 200 : integrator.alg.max_stages
  cache.constantcache.min_stage = (integrator.alg.min_stages > cache.constantcache.max_stage) ? cache.constantcache.max_stage : integrator.alg.min_stages

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

@muladd function perform_step!(integrator, cache::ROCK2Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack k, tmp, uᵢ₋₂, uᵢ₋₁, atmp = cache
  @unpack ms, fp1, fp2, recf = cache.constantcache
  ccache = cache.constantcache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((1.5 + abs(dt)*integrator.eigen_est)/0.811) + 1))
  mdeg = min(max(mdeg,ccache.min_stage), ccache.max_stage)
  ccache.mdeg = max(mdeg, 3) - 2
  choosedeg!(cache)
  # recurrence
  # for the first stage
  tᵢ₋₁ = t + dt*recf[ccache.start]
  tᵢ₋₂ = t + dt*recf[ccache.start]
  tᵢ₋₃ = t
  @.. uᵢ₋₂ = uprev
  @.. uᵢ₋₁ = uprev + (dt*recf[ccache.start])*fsalfirst
  ccache.mdeg < 2 && ( @.. u = uᵢ₋₁ )
  # for the second to the ms[ccache.mdeg] th stages
  for i in 2:ccache.mdeg
    μ, κ = recf[ccache.start + (i - 2)*2 + 1], recf[ccache.start + (i - 2)*2 + 2]
    ν = -1 - κ
    f(k, uᵢ₋₁, p, tᵢ₋₁)
    tᵢ₋₁ = dt*μ - ν*tᵢ₋₂ - κ*tᵢ₋₃
    @.. u = (dt*μ)*k - ν*uᵢ₋₁ - κ*uᵢ₋₂
    if i < ccache.mdeg
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = u
    end
    tᵢ₋₃ = tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
  end # end if
  # two-stage finishing procedure.
  δt₁ = dt*fp1[ccache.deg_index]
  δt₂ = dt*fp2[ccache.deg_index]
  f(k, u, p, tᵢ₋₁)
  integrator.destats.nf += 1
  @.. uᵢ₋₁ = u + δt₁*k
  if integrator.opts.adaptive
    @.. tmp = - δt₂*k
  else
    @.. u = -δt₂*k
  end
  c = DiffEqBase.value(sign(δt₁))*integrator.opts.internalnorm(δt₁,t)
  tᵢ₋₁ += c
  f(k, uᵢ₋₁, p, tᵢ₋₁)
  integrator.destats.nf += 1

  if integrator.opts.adaptive
    @.. tmp += δt₂*k
    @.. u = uᵢ₋₁ + δt₁*k + tmp
  else
    @.. u += uᵢ₋₁ + (δt₁ + δt₂)*k
  end

  # error estimate
  if integrator.opts.adaptive
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::ROCK4ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  cache.max_stage = (integrator.alg.max_stages < 1 || integrator.alg.max_stages > 152) ? 152 : integrator.alg.max_stages
  cache.min_stage = (integrator.alg.min_stages > cache.max_stage) ? cache.max_stage : integrator.alg.min_stages
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ROCK4ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, fpa, fpb, fpbe, recf = cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((3 + abs(dt)*integrator.eigen_est)/0.353) + 1))
  mdeg = min(max(mdeg,cache.min_stage), cache.max_stage)
  cache.mdeg = max(mdeg, 5) - 4
  choosedeg!(cache)
  # recurrence
  # for the first stage
  tᵢ₋₁ = t + dt*recf[cache.start]
  tᵢ₋₂ = t + dt*recf[cache.start]
  tᵢ₋₃ = t
  uᵢ₋₂ = copy(uprev)
  uᵢ₋₁ = uprev + (dt*recf[cache.start])*fsalfirst
  cache.mdeg < 2 && ( u = uᵢ₋₁ )
  # for the second to the cache.mdeg th stages
  for i in 2:cache.mdeg
    μ, κ = recf[cache.start + (i - 2)*2 + 1], recf[cache.start + (i - 2)*2 + 2]
    ν = - 1 - κ
    u = f(uᵢ₋₁, p,tᵢ₋₁)
    tᵢ₋₁ = dt*μ - ν*tᵢ₋₂ - κ*tᵢ₋₃
    u = (dt*μ)*u - ν*uᵢ₋₁ - κ*uᵢ₋₂
    i < cache.mdeg && (uᵢ₋₂ = uᵢ₋₁; uᵢ₋₁ = u)
    tᵢ₋₃ = tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
  end

  # These constants correspond to the Buther Tableau coefficients of explicit RK methods
  a₂₁ = dt*fpa[cache.deg_index][1]
  a₃₁ = dt*fpa[cache.deg_index][2]; a₃₂ = dt*fpa[cache.deg_index][3]
  a₄₁ = dt*fpa[cache.deg_index][4]; a₄₂ = dt*fpa[cache.deg_index][5]; a₄₃ = dt*fpa[cache.deg_index][6]
  B₁  = dt*fpb[cache.deg_index][1]; B₂  = dt*fpb[cache.deg_index][2]; B₃  = dt*fpb[cache.deg_index][3]; B₄  = dt*fpb[cache.deg_index][4]
  # coefficients of embedded method for error estimation
  B̂₁  = dt*(fpbe[cache.deg_index][1] - fpb[cache.deg_index][1])
  B̂₂  = dt*(fpbe[cache.deg_index][2] - fpb[cache.deg_index][2])
  B̂₃  = dt*(fpbe[cache.deg_index][3] - fpb[cache.deg_index][3])
  B̂₄  = dt*(fpbe[cache.deg_index][4] - fpb[cache.deg_index][4])
  B̂₅  = dt*fpbe[cache.deg_index][5]

  # 4-stage finishing procedure.
  # Stage-1
  uᵢ₋₁ = f(u, p, tᵢ₋₁)
  integrator.destats.nf += 1
  uᵢ₋₂ = u + a₃₁*uᵢ₋₁
  uᵢ₋₃ = u + a₄₁*uᵢ₋₁
  u    += B₁*uᵢ₋₁
  integrator.opts.adaptive && (tmp = B̂₁*uᵢ₋₁)
  uᵢ₋₁ = u + (a₂₁ - B₁)*uᵢ₋₁

  # Stage-2
  c₂ = a₂₁
  _c₂ = DiffEqBase.value(sign(c₂))*integrator.opts.internalnorm(c₂,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₂
  uᵢ₋₁ = f(uᵢ₋₁, p, tᵢ₋₂)
  integrator.destats.nf += 1
  uᵢ₋₂ += a₃₂*uᵢ₋₁
  uᵢ₋₃ += a₄₂*uᵢ₋₁
  u    += B₂*uᵢ₋₁
  integrator.opts.adaptive && (tmp += B̂₂*uᵢ₋₁)

  # Stage-3
  c₃ = a₃₁ + a₃₂
  _c₃ = DiffEqBase.value(sign(c₃))*integrator.opts.internalnorm(c₃,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₃
  uᵢ₋₂ = f(uᵢ₋₂, p, tᵢ₋₂)
  integrator.destats.nf += 1
  uᵢ₋₃ += a₄₃*uᵢ₋₂
  u    += B₃*uᵢ₋₂
  integrator.opts.adaptive && (tmp += B̂₃*uᵢ₋₂)

  #Stage-4
  c₄ = a₄₁ + a₄₂ + a₄₃
  _c₄ = DiffEqBase.value(sign(c₄))*integrator.opts.internalnorm(c₄,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₄
  uᵢ₋₃ = f(uᵢ₋₃, p, tᵢ₋₂)
  integrator.destats.nf += 1
  u    += B₄*uᵢ₋₃
  integrator.opts.adaptive && (tmp += B̂₄*uᵢ₋₃)

  uᵢ₋₁ = f(u, p, t + dt)
  integrator.destats.nf += 1

  #Error estimate (embedded method of order 3)
  if integrator.opts.adaptive
    tmp  += B̂₅*uᵢ₋₁
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = uᵢ₋₁
  integrator.u = u
end

function initialize!(integrator, cache::ROCK4Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  cache.constantcache.max_stage = (integrator.alg.max_stages < 1 || integrator.alg.max_stages > 152) ? 152 : integrator.alg.max_stages
  cache.constantcache.min_stage = (integrator.alg.min_stages > cache.constantcache.max_stage) ? cache.constantcache.max_stage : integrator.alg.min_stages

  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ROCK4Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack uᵢ₋₁, uᵢ₋₂, uᵢ₋₃, tmp, atmp, k= cache
  @unpack ms, fpa, fpb, fpbe, recf = cache.constantcache
  ccache = cache.constantcache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  mdeg = Int(floor(sqrt((3 + abs(dt)*integrator.eigen_est)/0.353) + 1))
  mdeg = min(max(mdeg,ccache.min_stage), ccache.max_stage)
  ccache.mdeg = max(mdeg, 5) - 4
  choosedeg!(cache)
  # recurrence
  # for the first stage
  tᵢ₋₁ = t + dt*recf[ccache.start]
  tᵢ₋₂ = t + dt*recf[ccache.start]
  tᵢ₋₃ = t
  @.. uᵢ₋₂ = uprev
  @.. uᵢ₋₁ = uprev + (dt*recf[ccache.start])*fsalfirst
  ccache.mdeg < 2 && ( @.. u = uᵢ₋₁ )
  # for the second to the ccache.mdeg th stages
  for i in 2:ccache.mdeg
    μ, κ = recf[ccache.start + (i - 2)*2 + 1], recf[ccache.start + (i - 2)*2 + 2]
    ν = -1 - κ
    f(k, uᵢ₋₁, p, tᵢ₋₁)
    tᵢ₋₁ = (dt*μ) - ν*tᵢ₋₂ - κ*tᵢ₋₃
    @.. u = (dt*μ)*k - ν*uᵢ₋₁ - κ*uᵢ₋₂
    if i < ccache.mdeg
      @.. uᵢ₋₂ = uᵢ₋₁
      @.. uᵢ₋₁ = u
    end
    tᵢ₋₃ = tᵢ₋₂
    tᵢ₋₂ = tᵢ₋₁
  end

  # These constants correspond to the Buther Tableau coefficients of explicit RK methods
  a₂₁ = dt*fpa[ccache.deg_index][1]
  a₃₁ = dt*fpa[ccache.deg_index][2]; a₃₂ = dt*fpa[ccache.deg_index][3]
  a₄₁ = dt*fpa[ccache.deg_index][4]; a₄₂ = dt*fpa[ccache.deg_index][5]; a₄₃ = dt*fpa[ccache.deg_index][6]
  B₁  = dt*fpb[ccache.deg_index][1]; B₂  = dt*fpb[ccache.deg_index][2]; B₃  = dt*fpb[ccache.deg_index][3]; B₄  = dt*fpb[ccache.deg_index][4]
  # coefficients of embedded method for error estimation
  B̂₁  = dt*(fpbe[ccache.deg_index][1] - fpb[ccache.deg_index][1])
  B̂₂  = dt*(fpbe[ccache.deg_index][2] - fpb[ccache.deg_index][2])
  B̂₃  = dt*(fpbe[ccache.deg_index][3] - fpb[ccache.deg_index][3])
  B̂₄  = dt*(fpbe[ccache.deg_index][4] - fpb[ccache.deg_index][4])
  B̂₅  = dt*fpbe[ccache.deg_index][5]

  # 4-stage finishing procedure.
  # Stage-1

  f(k, u, p, tᵢ₋₁)
  integrator.destats.nf += 1
  @.. uᵢ₋₂ = u + a₃₁*k
  @.. uᵢ₋₃ = u + a₄₁*k
  @.. uᵢ₋₁ = u + a₂₁*k
  @.. u    += B₁*k
  integrator.opts.adaptive && (@.. tmp = B̂₁*k)

  # Stage-2
  c₂ = a₂₁
  _c₂ = DiffEqBase.value(sign(c₂))*integrator.opts.internalnorm(c₂,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₂
  f(k, uᵢ₋₁, p, tᵢ₋₂)
  integrator.destats.nf += 1
  @.. uᵢ₋₂ += a₃₂*k
  @.. uᵢ₋₃ += a₄₂*k
  @.. u    += B₂*k
  integrator.opts.adaptive && (@.. tmp += B̂₂*k)

  # Stage-3
  c₃ = a₃₁ + a₃₂
  _c₃ = DiffEqBase.value(sign(c₃))*integrator.opts.internalnorm(c₃,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₃
  f(k, uᵢ₋₂, p, tᵢ₋₂)
  integrator.destats.nf += 1
  @.. uᵢ₋₃ += a₄₃*k
  @.. u    += B₃*k
  integrator.opts.adaptive && (@.. tmp += B̂₃*k)

  #Stage-4
  c₄ = a₄₁ + a₄₂ + a₄₃
  _c₄ = DiffEqBase.value(sign(c₄))*integrator.opts.internalnorm(c₄,t)
  tᵢ₋₂ = tᵢ₋₁ + _c₄
  f(k, uᵢ₋₃, p, tᵢ₋₂)
  integrator.destats.nf += 1
  @.. u    += B₄*k
  integrator.opts.adaptive && (tmp += B̂₄*k)

  f(k, u, p, t + dt)
  integrator.destats.nf += 1

  #Error estimate (embedded method of order 3)
  if integrator.opts.adaptive
    tmp  += B̂₅*k
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  @.. integrator.fsallast = k
  integrator.k[1] = integrator.fsalfirst
  integrator.destats.nf += 1
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::RKCConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::RKCConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  #maxm = max(2,Int(floor(sqrt(integrator.opts.internalnorm(integrator.opts.reltol,t)/(10*eps(integrator.opts.internalnorm(uprev,t)))))))
  maxm = 50
  mdeg = 1 + Int(floor(sqrt(1.54*abs(dt)*integrator.eigen_est + 1)))
  mdeg = (mdeg > maxm) ? maxm : mdeg

  w0 = 1 + 2/(13*(mdeg^2))
  temp1 = w0^2 - 1
  temp2 = sqrt(temp1)
  arg   = mdeg*log(w0 + temp2)
  w1    = (sinh(arg)*temp1) / (cosh(arg)*mdeg*temp2 - w0*sinh(arg))
  b1    = 1/((2*w0)^2)
  b2    = b1

  # stage-1
  gprev2 = copy(uprev)
  μs     = w1*b1
  gprev  = uprev + dt*μs*fsalfirst
  th2  = zero(eltype(u))
  th1  = μs
  z1   = w0
  z2   = one(eltype(u))
  dz1  = one(eltype(u))
  dz2  = zero(eltype(u))
  d2z1 = zero(eltype(u))
  d2z2 = zero(eltype(u))

  # stage 2 - mdeg
  for iter in 2:mdeg
    z   = 2*w0*z1 - z2
    dz  = 2*w0*dz1 - dz2 + 2*z1
    d2z = 2*w0*d2z1 - d2z2 + 4*dz1
    b   = d2z/(dz^2)
    νs  = 1 - z1*b1
    μ   = (2*w0*b)/b1
    ν   = - b/b2
    μs  = μ*w1/w0
    #using u as temporary storage
    u   = f(gprev, p, t + dt*th1)
    integrator.destats.nf += 1
    u   = μ*gprev + ν*gprev2  + (1 - μ - ν)*uprev + dt*μs*(u - νs*fsalfirst)
    th  = μ*th1 + ν*th2 + μs*(1 - νs)
    if (iter < mdeg)
      gprev2 = gprev
      gprev  = u
      th2  = th1
      th1  = th
      b2   = b1
      b1   = b
      z2   = z1
      z1   = z
      dz2  = dz1
      dz1  = dz
      d2z2 = d2z1
      d2z1 = d2z
    end
  end
  # error estimate
  if integrator.opts.adaptive
    tmp = 0.8*(uprev - u) + 0.4*dt*(fsalfirst + gprev)
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.u = u
end

function initialize!(integrator, cache::RKCCache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst  # done by pointers, no copying
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::RKCCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack k, tmp, gprev2, gprev, atmp = cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  #maxm = max(2,Int(floor(sqrt(integrator.opts.internalnorm(integrator.opts.reltol,t)/10eps(t)))))
  maxm = 50
  mdeg = 1 + Int(floor(sqrt(1.54*abs(dt)*integrator.eigen_est + 1)))
  mdeg = (mdeg > maxm) ? maxm : mdeg

  w0 = 1 + 2/(13*(mdeg^2))
  temp1 = w0^2 - 1
  temp2 = sqrt(temp1)
  arg   = mdeg*log(w0 + temp2)
  w1    = (sinh(arg)*temp1) / (cosh(arg)*mdeg*temp2 - w0*sinh(arg))
  b1    = 1/((2*w0)^2)
  b2    = b1

  # stage-1
  @.. gprev2 = uprev
  μs     = w1*b1
  @.. gprev  = uprev + dt*μs*fsalfirst
  th2  = zero(eltype(u))
  th1  = μs
  z1   = w0
  z2   = one(eltype(u))
  dz1  = one(eltype(u))
  dz2  = zero(eltype(u))
  d2z1 = zero(eltype(u))
  d2z2 = zero(eltype(u))

  # stage 2 - mdeg
  for iter in 2:mdeg
    z   = 2*w0*z1 - z2
    dz  = 2*w0*dz1 - dz2 + 2*z1
    d2z = 2*w0*d2z1 - d2z2 + 4*dz1
    b   = d2z/(dz^2)
    νs  = 1 - z1*b1
    μ   = (2*w0*b)/b1
    ν   = - b/b2
    μs  = μ*w1/w0
    f(k, gprev, p, t + dt*th1)
    integrator.destats.nf += 1
    @.. u   = μ*gprev + ν*gprev2  + (1 - μ - ν)*uprev + dt*μs*(k - νs*fsalfirst)
    th  = μ*th1 + ν*th2 + μs*(1 - νs)
    if (iter < mdeg)
      gprev2 = gprev
      gprev  = u
      th2  = th1
      th1  = th
      b2   = b1
      b1   = b
      z2   = z1
      z1   = z
      dz2  = dz1
      dz1  = dz
      d2z2 = d2z1
      d2z1 = d2z
    end
  end
  # error estimate
  if integrator.opts.adaptive
    @.. tmp = 0.8*(uprev - u) + 0.4*dt*(fsalfirst + gprev)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::IRKCConstantCache)
  @unpack uprev, p, t = integrator
  @unpack f1, f2 = integrator.f
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  cache.du₁ = f1(uprev,p,t)
  cache.du₂ = f2(uprev,p,t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.fsalfirst = cache.du₁ + cache.du₂

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::IRKCConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p,alg,fsalfirst = integrator
  @unpack minm,du₁,du₂,nlsolver = cache
  @unpack f1, f2 = integrator.f
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  # The the number of degree for Chebyshev polynomial
  #maxm = max(2,Int(floor(sqrt(integrator.opts.internalnorm(integrator.opts.reltol,t)/(10 *eps(integrator.opts.internalnorm(uprev,t)))))))
  maxm = 50
  mdeg = 1 + Int(floor(sqrt(1.54*abs(dt)*integrator.eigen_est + 1)))
  mdeg = min(maxm, max(minm, mdeg))

  ω₀    = 1 + 2/(13 * (mdeg^2))
  temp₁ = ω₀^2 - 1
  temp₂ = sqrt(temp₁)
  θ     = mdeg*log(ω₀ + temp₂)
  ω₁    = (sinh(θ)*temp₁)/(cosh(θ)*mdeg*temp₂ - ω₀*sinh(θ))
  Bⱼ₋₂  = 1/(4 * ω₀^2)
  Bⱼ₋₁  = 1/ω₀

  #stage-1
  f1ⱼ₋₂  = du₁
  gprev2 = copy(uprev)
  μs     = ω₁*Bⱼ₋₁
  μs₁    = μs

  # initial guess for implicit part
  # if alg.extrapolant == :linear
  #   nlsolver.z = dt*du₁
  # else # :constant
  #   nlsolver.z = zero(u)
  # end

  nlsolver.z = dt*du₁

  nlsolver.tmp = uprev + dt*μs₁*du₂
  nlsolver.γ   = μs₁
  nlsolver.c   = μs
  markfirststage!(nlsolver)
  z = nlsolve!(nlsolver, integrator, cache, false)
  # nlsolvefail(nlsolver) && return
  gprev = nlsolver.tmp + μs₁*z

  Cⱼ₋₂   = zero(eltype(u))
  Cⱼ₋₁   = μs
  Tⱼ₋₁   = ω₀
  Tⱼ₋₂   = one(eltype(u))
  Tⱼ₋₁′  = one(eltype(u))
  Tⱼ₋₂′  = zero(eltype(u))
  Tⱼ₋₁″  = zero(eltype(u))
  Tⱼ₋₂″  = zero(eltype(u))

  #stage- 2...mdeg
  for iter in 2:mdeg
    Tⱼ   = 2*ω₀*Tⱼ₋₁ - Tⱼ₋₂
    Tⱼ′  = 2*ω₀*Tⱼ₋₁′ + 2*Tⱼ₋₁ - Tⱼ₋₂′
    Tⱼ″  = 2*ω₀*Tⱼ₋₁″ + 4*Tⱼ₋₁′ - Tⱼ₋₂″
    Bⱼ   = Tⱼ″/(Tⱼ′^2)
    μ    = (2*ω₀*Bⱼ)/Bⱼ₋₁
    ν    = - Bⱼ/Bⱼ₋₂
    μs   = (μ*ω₁)/ω₀
    νs   = -(1 - Tⱼ₋₁*Bⱼ₋₁)*μs
    Cⱼ   = μ*Cⱼ₋₁ + ν*Cⱼ₋₂ + μs + νs

    f1ⱼ₋₁  = f1(gprev, p, t+Cⱼ₋₁*dt)
    f2ⱼ₋₁  = f2(gprev, p, t+Cⱼ₋₁*dt)
    integrator.destats.nf += 1
    integrator.destats.nf2 += 1
    nlsolver.tmp = (1-μ-ν)*uprev + μ*gprev + ν*gprev2 + dt*μs*f2ⱼ₋₁ + dt*νs*du₂ + (νs - (1 -μ-ν)*μs₁)*dt*du₁ - ν*μs₁*dt*f1ⱼ₋₂
    nlsolver.z   = dt*f1ⱼ₋₁
    nlsolver.c   = Cⱼ
    z = nlsolve!(nlsolver, integrator, cache, false)
    # ignoring newton method's convergence failure
    # nlsolvefail(nlsolver) && return
    u = nlsolver.tmp + μs₁*z
    if (iter < mdeg)
      f1ⱼ₋₂= f1ⱼ₋₁
      gprev2 = gprev
      gprev  = u
      Cⱼ₋₂   = Cⱼ₋₁
      Cⱼ₋₁   = Cⱼ
      Bⱼ₋₂   = Bⱼ₋₁
      Bⱼ₋₁   = Bⱼ
      Tⱼ₋₂   = Tⱼ₋₁
      Tⱼ₋₁   = Tⱼ
      Tⱼ₋₂′  = Tⱼ₋₁′
      Tⱼ₋₁′  = Tⱼ′
      Tⱼ₋₂″  = Tⱼ₋₁″
      Tⱼ₋₁″  = Tⱼ″
    end
  end

  cache.du₁ = f1(u, p, t+dt)
  cache.du₂ = f2(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  # error estimate
  if isnewton(nlsolver) && integrator.opts.adaptive
    update_W!(integrator, cache, dt, false)
    tmp = dt*(0.5*(cache.du₂ - du₂) + (0.5 - μs₁)*(cache.du₁ - du₁))
    tmp = _reshape(get_W(nlsolver) \ _vec(tmp), axes(tmp))
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  integrator.fsallast = cache.du₁ + cache.du₂
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::IRKCCache)
  @unpack uprev, p, t = integrator
  @unpack f1, f2 = integrator.f
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f1(cache.du₁, uprev, p, t)
  f2(cache.du₂, uprev, p, t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  @.. integrator.fsalfirst = cache.du₁ + cache.du₂
end

function perform_step!(integrator, cache::IRKCCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p,alg = integrator
  @unpack gprev,gprev2,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,du₁,du₂,atmp,nlsolver = cache
  @unpack tmp,z = nlsolver
  @unpack minm = cache.constantcache
  @unpack f1, f2 = integrator.f

  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)
  # The the number of degree for Chebyshev polynomial
  #maxm = max(2,int(floor(sqrt(integrator.opts.internalnorm(integrator.opts.reltol,t)/(10 *eps(integrator.opts.internalnorm(uprev,t)))))))
  maxm = 50
  mdeg = 1 + Int(floor(sqrt(1.54*abs(dt)*integrator.eigen_est + 1)))
  mdeg = (mdeg < minm) ? minm : mdeg
  mdeg = (mdeg >= maxm) ? maxm : mdeg

  ω₀    = 1 + 2/(13 *(mdeg^2))
  temp₁ = ω₀^2 - 1
  temp₂ = sqrt(temp₁)
  θ     = mdeg*log(ω₀ + temp₂)
  ω₁    = (sinh(θ)*temp₁)/(cosh(θ)*mdeg*temp₂ - ω₀*sinh(θ))
  Bⱼ₋₂  = 1 / (4 * ω₀^2)
  Bⱼ₋₁  = 1 / ω₀

  #stage-1
  f1ⱼ₋₂  = du₁
  @.. gprev2 = uprev
  μs     = ω₁*Bⱼ₋₁
  μs₁    = μs

  # initial guess
  # if alg.extrapolant == :linear
  #   @.. z = dt*du₁
  # else # :constant
  #   @.. z = zero(eltype(u))
  # end
  @.. nlsolver.z = dt*du₁

  @.. nlsolver.tmp = uprev + dt*μs₁*du₂
  nlsolver.γ   = μs₁
  nlsolver.c   = μs
  markfirststage!(nlsolver)
  z = nlsolve!(nlsolver, integrator, cache, false)
  # ignoring newton method's convergence failure
  # nlsolvefail(nlsolver) && return
  @.. gprev = nlsolver.tmp + μs₁*nlsolver.z

  Cⱼ₋₂   = zero(eltype(u))
  Cⱼ₋₁   = μs
  Tⱼ₋₁   = ω₀
  Tⱼ₋₂   = one(eltype(u))
  Tⱼ₋₁′  = one(eltype(u))
  Tⱼ₋₂′  = zero(eltype(u))
  Tⱼ₋₁″  = zero(eltype(u))
  Tⱼ₋₂″  = zero(eltype(u))

  #stage- 2...mdeg
  for iter in 2:mdeg
    Tⱼ   = 2*ω₀*Tⱼ₋₁ - Tⱼ₋₂
    Tⱼ′  = 2*ω₀*Tⱼ₋₁′ + 2*Tⱼ₋₁ - Tⱼ₋₂′
    Tⱼ″  = 2*ω₀*Tⱼ₋₁″ + 4*Tⱼ₋₁′ - Tⱼ₋₂″
    Bⱼ   = Tⱼ″/(Tⱼ′^2)
    μ    = (2*ω₀*Bⱼ)/Bⱼ₋₁
    ν    = - Bⱼ/Bⱼ₋₂
    μs   = (μ*ω₁)/ω₀
    νs   = -(1 - Tⱼ₋₁*Bⱼ₋₁)*μs
    Cⱼ   = μ*Cⱼ₋₁ + ν*Cⱼ₋₂ + μs + νs

    f1(f1ⱼ₋₁, gprev, p, t+Cⱼ₋₁*dt)
    f2(f2ⱼ₋₁, gprev, p, t+Cⱼ₋₁*dt)
    integrator.destats.nf += 1
    integrator.destats.nf2 += 1
    @.. nlsolver.tmp = (1-μ-ν)*uprev + μ*gprev + ν*gprev2 + dt*μs*f2ⱼ₋₁ + dt*νs*du₂ + (νs - (1-μ-ν)*μs₁)*dt*du₁ - ν*μs₁*dt*f1ⱼ₋₂
    @.. nlsolver.z   = dt*f1ⱼ₋₁
    nlsolver.c = Cⱼ

    z = nlsolve!(nlsolver, integrator, cache, false)
    # nlsolvefail(nlsolver) && return
    @.. u = nlsolver.tmp + μs₁*nlsolver.z
    if (iter < mdeg)
      @.. f1ⱼ₋₂  = f1ⱼ₋₁
      @.. gprev2 = gprev
      @.. gprev  = u
      Cⱼ₋₂   = Cⱼ₋₁
      Cⱼ₋₁   = Cⱼ
      Bⱼ₋₂   = Bⱼ₋₁
      Bⱼ₋₁   = Bⱼ
      Tⱼ₋₂   = Tⱼ₋₁
      Tⱼ₋₁   = Tⱼ
      Tⱼ₋₂′  = Tⱼ₋₁′
      Tⱼ₋₁′  = Tⱼ′
      Tⱼ₋₂″  = Tⱼ₋₁″
      Tⱼ₋₁″  = Tⱼ″
    end
  end

  @.. f1ⱼ₋₁ = du₁
  @.. f2ⱼ₋₁ = du₂
  f1(du₁, u, p, t+dt)
  f2(du₂, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  # error estimate
  if isnewton(nlsolver) && integrator.opts.adaptive
    update_W!(integrator, cache, dt, false)
    @.. gprev = dt*0.5*(du₂ - f2ⱼ₋₁) + dt*(0.5 - μs₁)*(du₁ - f1ⱼ₋₁)
    nlsolver.cache.linsolve(vec(tmp),get_W(nlsolver),vec(gprev),false)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  @.. integrator.fsallast = du₁ + du₂
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::ESERK4ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ESERK4ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, Cᵤ, Cₑ= cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est))+1)
  mdeg = (mdeg > 4000) ? 4000 : mdeg
  cache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = cache.mdeg
  start = cache.start
  internal_deg = cache.internal_deg
  α = 2.0/(mdeg^2)

  u = zero(uprev)
  tmp = zero(uprev)

  for i in 1:4
    hᵢ = dt/i
    tᵢ = t
    Sᵢ = zero(u)
    uᵢ₋₁ = uprev
    uᵢ₋₂ = zero(u)
    for j in 1:i
      r  = tᵢ
      Sᵢ = (cache.Bᵢ[start])*uᵢ₋₁
      for st in 1:mdeg
        k = f(uᵢ₋₁, p, r)
        integrator.destats.nf += 1

        if st%internal_deg == 1
          uᵢ = uᵢ₋₁ + α*hᵢ*k
        else
          uᵢ = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*hᵢ*k
        end
        q = convert(Int, floor(st/internal_deg))
        r = tᵢ + α*(st^2 + q*internal_deg^2)*hᵢ
        Sᵢ = Sᵢ + (cache.Bᵢ[start+st])*uᵢ
        if st < mdeg
          uᵢ₋₂ = uᵢ₋₁
          uᵢ₋₁ = uᵢ
        end
      end

      if j < i
        tᵢ = tᵢ + hᵢ
        uᵢ₋₁ = Sᵢ
      end
    end

    u = u + Cᵤ[i]*Sᵢ
    integrator.opts.adaptive && (tmp = tmp + Cₑ[i]*Sᵢ)
  end

  u = u/6
  if integrator.opts.adaptive
    tmp = tmp/6
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.u = u
end

function initialize!(integrator, cache::ESERK4Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ESERK4Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack uᵢ, uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, k = cache
  @unpack ms, Cᵤ, Cₑ = cache.constantcache
  ccache = cache.constantcache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est))+1)
  mdeg = (mdeg > 4000) ? 4000 : mdeg
  ccache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = ccache.mdeg
  start = ccache.start
  internal_deg = ccache.internal_deg
  α = 2.0/(mdeg^2)

  @.. u = zero(uprev)
  @.. tmp = zero(uprev)
  for i in 1:4
    hᵢ = dt/i
    tᵢ = t
    @.. Sᵢ = zero(u)
    @.. uᵢ₋₁ = uprev
    @.. uᵢ₋₂ = zero(u)
    for j in 1:i
      r  = tᵢ
      @.. Sᵢ = (cache.constantcache.Bᵢ[start])*uᵢ₋₁
      for st in 1:mdeg
        f(k, uᵢ₋₁, p, r)
        integrator.destats.nf += 1

        if st%internal_deg == 1
          @.. uᵢ = uᵢ₋₁ + α*hᵢ*k
        else
          @.. uᵢ = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*hᵢ*k
        end
        q = convert(Int, floor(st/internal_deg))
        r = tᵢ + α*(st^2 + q*internal_deg^2)*hᵢ
        @.. Sᵢ = Sᵢ + (cache.constantcache.Bᵢ[start+st])*uᵢ
        if st < mdeg
          @.. uᵢ₋₂ = uᵢ₋₁
          @.. uᵢ₋₁ = uᵢ
        end
      end

      if j < i
        tᵢ = tᵢ + hᵢ
        @.. uᵢ₋₁ = Sᵢ
      end
    end

    @.. u = u + Cᵤ[i]*Sᵢ
    integrator.opts.adaptive && (@.. tmp = tmp + Cₑ[i]*Sᵢ)
  end

  @.. u = u/6

  if integrator.opts.adaptive
    @.. tmp = tmp/6
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::ESERK5ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ESERK5ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, Cᵤ, Cₑ, Bᵢ= cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est/0.98))+1)
  mdeg = (mdeg > 2000) ? 2000 : mdeg
  cache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = cache.mdeg
  start = cache.start
  internal_deg = cache.internal_deg
  α = 100.0/(49.0*mdeg^2)

  u = zero(uprev)
  tmp = zero(uprev)
  for i in 1:5
    hᵢ = dt/i
    tᵢ = t
    Sᵢ = zero(u)
    uᵢ₋₁ = uprev
    uᵢ₋₂ = zero(u)
    for j in 1:i
      r  = tᵢ
      Sᵢ = (Bᵢ[start])*uᵢ₋₁
      for st in 1:mdeg
        k = f(uᵢ₋₁, p, r)
        integrator.destats.nf += 1

        if st%internal_deg == 1
          uᵢ = uᵢ₋₁ + α*hᵢ*k
        else
          uᵢ = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*hᵢ*k
        end
        q = convert(Int, floor(st/internal_deg))
        r = tᵢ + α*(st^2 + q*internal_deg^2)*hᵢ
        Sᵢ = Sᵢ + (Bᵢ[start+st])*uᵢ
        if st < mdeg
          uᵢ₋₂ = uᵢ₋₁
          uᵢ₋₁ = uᵢ
        end
      end

      if j < i
        tᵢ = tᵢ + hᵢ
        uᵢ₋₁ = Sᵢ
      end
    end

    u = u + Cᵤ[i]*Sᵢ
    integrator.opts.adaptive && (tmp = tmp + Cₑ[i]*Sᵢ)
  end

  u = u/24
  if integrator.opts.adaptive
    tmp = tmp/24
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.u = u
end

function initialize!(integrator, cache::ESERK5Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ESERK5Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack uᵢ, uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, k = cache
  @unpack ms, Cᵤ, Cₑ, Bᵢ = cache.constantcache
  ccache = cache.constantcache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est/0.98))+1)
  mdeg = (mdeg > 2000) ? 2000 : mdeg
  ccache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = ccache.mdeg
  start = ccache.start
  internal_deg = ccache.internal_deg
  α = 100.0/(49.0*mdeg^2)

  @.. u = zero(uprev)
  @.. tmp = zero(uprev)
  for i in 1:5
    hᵢ = dt/i
    tᵢ = t
    @.. Sᵢ = zero(u)
    @.. uᵢ₋₁ = uprev
    @.. uᵢ₋₂ = zero(u)
    for j in 1:i
      r  = tᵢ
      @.. Sᵢ = (Bᵢ[start])*uᵢ₋₁
      for st in 1:mdeg
        f(k, uᵢ₋₁, p, r)
        integrator.destats.nf += 1

        if st%internal_deg == 1
          @.. uᵢ = uᵢ₋₁ + α*hᵢ*k
        else
          @.. uᵢ = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*hᵢ*k
        end
        q = convert(Int, floor(st/internal_deg))
        r = tᵢ + α*(st^2 + q*internal_deg^2)*hᵢ
        @.. Sᵢ = Sᵢ + (Bᵢ[start+st])*uᵢ
        if st < mdeg
          @.. uᵢ₋₂ = uᵢ₋₁
          @.. uᵢ₋₁ = uᵢ
        end
      end

      if j < i
        tᵢ = tᵢ + hᵢ
        @.. uᵢ₋₁ = Sᵢ
      end
    end

    @.. u = u + Cᵤ[i]*Sᵢ
    integrator.opts.adaptive && (@.. tmp = tmp + Cₑ[i]*Sᵢ)
  end

  @.. u = u/24

  if integrator.opts.adaptive
    @.. tmp = tmp/24
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::SERK2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1
  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::SERK2ConstantCache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack ms, Bᵢ= cache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est/0.8))+1)
  mdeg = (mdeg > 250) ? 250 : mdeg
  cache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = cache.mdeg
  start = cache.start
  internal_deg = cache.internal_deg
  α = 1.0/(0.4*mdeg^2)

  uᵢ₋₁ = uprev
  uᵢ₋₂ = uprev
  Sᵢ   = Bᵢ[start]*uprev
  for i in 1:10
    k = f(uᵢ₋₁, p, t+(1+(i-1)*internal_deg^2)*α*dt)
    integrator.destats.nf += 1
    u    = uᵢ₋₁ + α*dt*k
    Sᵢ   = Sᵢ + Bᵢ[start + (i-1)*internal_deg + 1]*u
    uᵢ₋₂ = uᵢ₋₁
    uᵢ₋₁ = u
    for j in 2:internal_deg
      k = f(uᵢ₋₁, p, t+(j^2+(i-1)*internal_deg^2)*α*dt)
      integrator.destats.nf += 1
      u = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*dt*k
      Sᵢ= Sᵢ + Bᵢ[start+j+(i-1)*internal_deg]*u
      if j*i < mdeg
        uᵢ₋₂ = uᵢ₋₁
        uᵢ₋₁ = u
      end
    end
  end
  u = Sᵢ
  k = f(u,p,t+dt)
  integrator.destats.nf += 1

  if integrator.opts.adaptive
    tmp = u - uprev - dt*k
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = k
  integrator.u = u
end

function initialize!(integrator, cache::SERK2Cache)
  integrator.kshortsize = 2
  resize!(integrator.k, integrator.kshortsize)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::SERK2Cache, repeat_step=false)
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  @unpack uᵢ₋₁, uᵢ₋₂, Sᵢ, tmp, atmp, k = cache
  @unpack ms, Bᵢ = cache.constantcache
  ccache = cache.constantcache
  alg = unwrap_alg(integrator, true)
  alg.eigen_est === nothing ? maxeig!(integrator, cache) : alg.eigen_est(integrator)

  mdeg = Int(floor(sqrt(abs(dt)*integrator.eigen_est/0.8))+1)
  mdeg = (mdeg > 250) ? 250 : mdeg
  ccache.mdeg = mdeg
  choosedeg_SERK!(integrator,cache)
  mdeg = ccache.mdeg
  start = ccache.start
  internal_deg = ccache.internal_deg
  α = 1.0/(0.4*mdeg^2)

  @.. uᵢ₋₁ = uprev
  @.. uᵢ₋₂ = uprev
  @.. Sᵢ   = Bᵢ[start]*uprev
  for i in 1:10
    f(k, uᵢ₋₁, p, t+(1+(i-1)*internal_deg^2)*α*dt)
    integrator.destats.nf += 1
    @.. u    = uᵢ₋₁ + α*dt*k
    @.. Sᵢ   = Sᵢ + Bᵢ[start + (i-1)*internal_deg + 1]*u
    @.. uᵢ₋₂ = uᵢ₋₁
    @.. uᵢ₋₁ = u
    for j in 2:internal_deg
      f(k, uᵢ₋₂, p, t+(j^2+(i-1)*internal_deg^2)*α*dt)
      integrator.destats.nf += 1
      @.. u = 2*uᵢ₋₁ - uᵢ₋₂ + 2*α*dt*k
      @.. Sᵢ= Sᵢ + Bᵢ[start+j+(i-1)*internal_deg]*u
      if j < mdeg
        @.. uᵢ₋₂ = uᵢ₋₁
        @.. uᵢ₋₁ = u
      end
    end
  end
  @.. u = Sᵢ
  f(k, u, p, t+dt)
  integrator.destats.nf += 1

  if integrator.opts.adaptive
    @.. tmp = u - uprev - dt*k
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast = k
  integrator.u = u
end


text:
module BayesTree

  using
    Distributions,
    StatsBase

  export
   #type
   Decision,
   DecisionTree,
   DecisionNode,
   DecisionBranch,
   RegressionLeaf,
   ClassificationLeaf,
   BartProbability,
   BartOptions,
   BartStateParameters,
   BartTree,
   BartState,
   BartAdditiveTree,
   Bart,
   BartLeaf,
   Tree,
   Node,
   Branch,
   Leaf,

   #methods
   leaves,
   valid_node,
   bart_options,
   depth,
   branches,
   grand_branches,
   not_grand_branches,
   length,
   parent,
   node_nonterminal,
   log_tree_prior,
   node_loglikelihood,
   fit,
   predict,
   node_grow!,
   node_prune!,
   node_grow_prune!,
   change_decision_rule!,
   swap_decision_rule!,
   normalize,
   denormalize,
   initialize_bart_state,
   model_selection

   include("tree.jl")
   include("decision_tree.jl")
   include("bart_new.jl")

end # module

max_stars_repo_name: mathcg/BayesTree

max_issues_repo_name: mathcg/BayesTree

max_forks_repo_name: mathcg/BayesTree


text:
@testset "Graph Transformer Layer" begin
    num_nodes = 4
    num_edges = 10
    hidden_dim = 16
    num_heads = 2
    g = SimpleDiGraph(num_nodes)
    add_edge!(g, 1, 1)
    add_edge!(g, 2, 2)
    add_edge!(g, 3, 3)
    add_edge!(g, 4, 4)
    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 3)
    add_edge!(g, 2, 4)
    add_edge!(g, 3, 4)
    node_features = rand(hidden_dim, num_nodes)    # [feature_dim, num_num_nodes]
    edge_features = rand(hidden_dim, num_edges)   # [feature_dim, num_num_edges]
    graph = FeaturedGraph(g, nf=node_features, ef=edge_features, directed=:directed)

    layer = GraphTransformerLayer(hidden_dim, 2)
    output_graph = layer(graph)
    @test size(output_graph.nf) == (hidden_dim, num_nodes)
    @test size(output_graph.ef) == (hidden_dim, num_edges)

    # Test the gradients
    y = rand(hidden_dim, num_nodes + num_edges)
    grads = gradient(params(layer)) do
        new_graph = layer(graph)
        logits = cat(new_graph.nf, new_graph.ef, dims=2)
        return Flux.logitcrossentropy(logits, y)
    end
    for param in params(layer)
        @test size(param) == size(grads[param])
    end
end

max_stars_repo_name: DwaraknathT/Flux-Graph-Transformer

max_issues_repo_name: DwaraknathT/Flux-Graph-Transformer

max_forks_repo_name: DwaraknathT/Flux-Graph-Transformer


text:
struct DSTU2 <: FHIRVersion
end

struct STU3 <: FHIRVersion
end

struct UnknownFHIRVersion <: FHIRVersion
end

struct OtherFHIRVersion{V} <: FHIRVersion
end


text:
# In Base/intfuncs.jl, x^p returns zero(x) when p == 0
# Since one(PolyVar) and one(Monomial) do not return
# a PolyVar and a Monomial, this results in type instability
# Defining the specific methods solve this problem and also make
# them a lot faster
Base.:(^)(x::PolyVar{C}, i::Int) where {C} = Monomial{C}([x], [i])
Base.:(^)(x::Monomial{true}, i::Int) = Monomial{true}(x.vars, i*x.z)

myminivect(x::T, y::T) where {T} = [x, y]
function myminivect(x::S, y::T) where {S,T}
    U = promote_type(S, T)
    [U(x), U(y)]
end

Base.:(+)(x::DMonomialLike, y::DMonomialLike) = Term(x) + Term(y)
Base.:(-)(x::DMonomialLike, y::DMonomialLike) = Term(x) - Term(y)

_getindex(p::Polynomial, i) = p[i]
_getindex(t::Term, i) = t
function plusorminus(p::TermPoly{C, S}, q::TermPoly{C, T}, op) where {C, S, T}
    varsvec = [_vars(p), _vars(q)]
    allvars, maps = mergevars(varsvec)
    nvars = length(allvars)
    U = Base.promote_op(op, S, T)
    a = Vector{U}()
    Z = Vector{Vector{Int}}()
    i = j = 1
    while i <= nterms(p) || j <= nterms(q)
        z = zeros(Int, nvars)
        if j > nterms(q) || (i <= nterms(p) && _getindex(p, i).x > _getindex(q, j).x)
            t = _getindex(p, i)
            z[maps[1]] = t.x.z
            α = convert(U, t.α)
            i += 1
        elseif i > nterms(p) || _getindex(q, j).x > _getindex(p, i).x
            t = _getindex(q, j)
            z[maps[2]] = t.x.z
            α = convert(U, op(t.α))
            j += 1
        else
            t = _getindex(p, i)
            z[maps[1]] = t.x.z
            s = _getindex(q, j)
            α = op(t.α, s.α)
            i += 1
            j += 1
        end
        push!(a, α)
        push!(Z, z)
    end

    Polynomial(a, MonomialVector{C}(allvars, Z))
end

Base.:(+)(x::TermPoly{C}, y::TermPoly{C}) where C = plusorminus(x, y, +)
Base.:(-)(x::TermPoly{C}, y::TermPoly{C}) where C = plusorminus(x, y, -)
Base.:(+)(x::TermPoly{C}, y::Union{Monomial,PolyVar}) where C = x + Term{C}(y)
Base.:(+)(x::Union{Monomial,PolyVar}, y::TermPoly{C}) where C = Term{C}(x) + y

Base.:(-)(x::TermPoly{T}, y::DMonomialLike) where T = x - Term{T}(y)
Base.:(-)(x::DMonomialLike, y::TermPoly{T}) where T = Term{T}(x) - y

Base.:(-)(p::Polynomial) = Polynomial(-p.a, p.x)

include("mult.jl")


text:
using JLD
using Plots
using JSON

# local
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/gradient_utils.jl")
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/functions_julia.jl")
# path="/home/genis/cluster_home/model_fitting_julia/"

#portatil
#
path="/home/emma/fit_behaviour_julia/"
# include(path*"/scripts/functions/functions_julia.jl")
# include(path*"/scripts/functions/utils.jl")
# include(path*"/scripts/functions/analysis.jl")
# include(path*"/scripts/functions/parameter_functions.jl")
# include(path*"/scripts/functions/simulations.jl")

SUBJECTS=[1,2,4,5]

#model_fit="urgency_x0_biasd_st3_dxdt42"
model_fit="StLinearDw_dx"

##### PARAM simulations ######3

tau=1.0
Tframe=0.5
Nframes=10
T=Tframe*Nframes

dt=tau/100.0
NT=T/dt
#############################

#plotly(ticks=:native)
#fig,axes=subplots(4,4,figsize=(20/2.54,20/2.54))

#colors_subject=["#cf4396","#59b94c","#a258c6","#a8b546","#5a6dc5",
#"#d7a044","#ab8ad0","#5a8437","#d43f4b","#60c395","#d47baf",
#"#358461","#c8683a","#4faad1","#8c7537","#b75d6c"]

for isubject in 1:length(SUBJECTS)

    subject=SUBJECTS[isubject]

    path_save=path*"figures/"
    filename_save=path_save*"check_PR_fit_PR_sim"*model_fit*"_subject"*string(subject)*".jld"
    data=JLD.load(filename_save)
    PR_sim=data["PR_sim"]
    PR_fit=data["PR_fit"]
    #axes[isubject].plot([0,1],[0,1],"k--")
    #axes[isubject].plot(PR_sim,PR_fit[1:length(PR_sim)],".",color=colors_subject[isubject])
    #axes[isubject].set_title(string(subject) )
    #if isubject==4
        #axes[isubject].set_xlabel("PR simulations")
    #end
    #if isubject==1 || isubject==4
        #axes[isubject].set_ylabel("PR Brunton")
    #end
end

#fig.tight_layout()
#fig.show()
#path_save=path*"figures/"
#fig.savefig(path_save*"PR_sim_PR_fit_model"*model_fit*".png")

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using WhereTraits
using Test

# Test standard dispatch
# ======================
  @macroexpand @traits g(a) = a

@testset "standard dispatch" begin
  # s === standard
  @traits_test fs(a::A) where {A} = A
  @test fs(1) == Int
  @test fs("hi") == String
  @test_throws MethodError fs(1, 2)

  @traits_test fs2(a::Vector{A}) where {Integer<:A<:Number} = sum(a)
  @test fs2(convert(Vector{Integer}, [1,2,3])) == 6
  @test_throws MethodError fs2([1,2,3])
end

@traits fb2(a::Vector{A}) where {A, Base.isconcretetype(A)} = length(a)

# Test boolean dispatch
# =====================

@testset "bool dispatch" begin
  # b === bool
  @traits_test fb(a::Vector) where {!isempty(a)} = a[1]
  @traits_test fb(a::Vector) where {isempty(a)} = nothing
  @traits_show_implementation fb

  @test isnothing(fb([]))
  @test fb([43]) == 43
  @test_throws MethodError fb((1,2,3))
  @traits_test fb2(a::Vector{A}) where {A, Base.isconcretetype(A)} = length(a)
  @traits_test fb2(a::Vector{A}) where {A, !Base.isconcretetype(A)} = length(a) + 1
  @test fb2([]) == 1  # Any is not concrete, length == 0
  @test fb2([3,5,.6]) == 3  # Float is concrete, length == 3
  @test_throws MethodError fb2((1,2,3))

  @traits_test fb3(a) where {a} = 1
  @traits_test fb3(a) where {!a} = 2

  @traits_show_implementation fb3

  @test fb3(true) == 1
  @test fb3(false) == 2
  @test_throws MethodError fb3(true, false)
end

# Test :: dispatch
# ================

@testset ":: dispatch" begin
  #cc === colon colon
  @traits_test fcc(a::A) where {A, Base.IteratorSize(A)::Base.HasLength} = length(a)
  @traits_test fcc(a::A) where {A, Base.IteratorSize(A)::Base.HasShape} = size(a)
  @traits_test fcc(a::A) where {A, Base.IteratorSize(A)::Base.SizeUnknown} = nothing
  @traits_show_implementation fcc

  @test fcc(1:4) == (4,)
  @test fcc("hallo") == 5
  @test_throws MethodError fcc(Base.Iterators.repeated(2))  # IteratorSize == Infinite

  @traits_test fcc2(a) where {eltype(a)::Type{Int}} = 5 + sum(a)
  @traits_test fcc2(a) where {eltype(a)::Type{String}} = "[$(join(a, ","))]"
  @test fcc2([1,2,3,4]) == 15
  @test fcc2(["a", "b"]) == "[a,b]"
  @test_throws MethodError fcc2([2.0, 5.0])

  @traits_test fcc3(a::Val{T}) where {T, T::Int} = T + 4
  @traits_test fcc3(a::Val{T}) where {T, T::Symbol} = T
  @test fcc3(Val(1)) == 5
  @test fcc3(Val(:hi)) == :hi
  @test_throws MethodError fcc3(Val(true))
end

# Test <: dispatch
# ================

@testset "<: dispatch" begin
  #sc === smaller colon
  @traits_test fsc(a::A) where {A, eltype(A)<:Int} = 5 + sum(a)
  @traits_test fsc(a::A) where {A, eltype(A)<:String} = "[$(join(a, ","))]"
  @test fsc([1,2,3,4]) == 15
  @test fsc(["a", "b"]) == "[a,b]"
  @test_throws MethodError fsc([2.0, 5.0])

  @traits_test fsc2(a) where {eltype(a)<:Int} = 5 + sum(a)
  @traits_test fsc2(a) where {eltype(a)<:String} = "[$(join(a, ","))]"
  @test fsc2([1,2,3,4]) == 15
  @test fsc2(["a", "b"]) == "[a,b]"
  @traits_show_implementation fsc2
  @test_throws MethodError fsc2([2.0, 5.0])

  @traits_test fsc2(a) where {eltype(a)<:Int} = 5 + sum(a)
  @traits_test fsc2(a) where {eltype(a)<:String} = "[$(join(a, ","))]"
  @test fsc2([1,2,3,4]) == 15
  @test fsc2(["a", "b"]) == "[a,b]"
  @traits_show_implementation fsc2
  @test_throws MethodError fsc2([2.0, 5.0])
end

# Test default
# ============

@testset "default args" begin
  @traits_test fda(a, b = 3) = a + b
  @traits_show_implementation fda
  @test fda(1) == 4

  # defaults will be preserved because they have been translated to extra methods
  # which are currently not tracked by @traits_test due to complexity
  # this is the same behaviour as for standard julia functions
  @traits_test fda(a, b) = a + b + 9
  @test fda(1) == 13
  @test_throws MethodError fda()

  # but they can be overwritten as usual
  @traits_test fda(a, b=1) = a + b + 99
  @test fda(1) == 101

  # or more drastically by an overwrite
  @traits_test fda(a) = 3 + a
  @test fda(1) == 4

  # also reassigning default values will overwrite main clause
  @traits_test fda(a, b=50) = a + b + 999
  @test fda(1) == 1050

  # check that also traits get lowered correctly
  @traits_test fda2(a::A, b::B=50) where {A<:Number, B<:Number, eltype(A) == Int, A == B, eltype(b) <: Number} = a + b + 999
  @traits_show_implementation fda2
  @test fda2(1) == 1050  # should not throw "UndefVarError: B not defined"
  @test_throws MethodError fda2(1.0)
  @test_throws MethodError fda2(1.0, 4)
end

# Test kwargs
# ===========

@testset "default kwargs" begin
  # kwargs are not used for standard julia dispatch, and neither for traits dispatch
  # they can be seen as a list of extra configurations for a concrete dispatch
  @traits_test fdkw(a; b = 3) = a + b
  @traits_show_implementation fdkw
  @test fdkw(1) == 4

  @traits_test fdkw(a; b=1) = a + b + 99
  @test fdkw(1) == 101

  # or more drastically by an overwrite
  @traits_test fdkw(a) = 3 + a
  @test fdkw(1) == 4

  # overwriting works
  @traits_test fdkw(a; b=50) = a + b + 999
  @test fdkw(1) == 1050

  # or more drastically by an overwrite
  @traits_test fdkw(a) = 3 + a
  @test fdkw(1) == 4

  # overwriting again works too
  @traits_test fdkw(a; b=50) = a + b + 999
  @test fdkw(1) == 1050
end

# Test overwritability
# ====================

@testset "overwritability" begin
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.HasLength} = length(a)
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.HasShape} = shape(a)
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.SizeUnknown} = 3
  state1 = @traits_get_state ow
  expr1 = @traits_show_implementation ow

  # @test length(expr1.args) == 4
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.SizeUnknown} = 8
  state2 = @traits_get_state ow
  expr2 = @traits_show_implementation ow

  @test length(state1) == length(state2)
  for (store1, store2) in zip(state1, state2)
    @test length(store1.definitions.inners) == length(store2.definitions.inners)
  end

  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.HasLength} = length(a)
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.HasShape} = size(a)
  @traits_test ow(a::A) where {A, Base.IteratorSize(A)::Base.SizeUnknown} = nothing
  state3 = @traits_get_state ow
  expr3 = @traits_show_implementation ow
  @test length(state1) == length(state3)
  for (store1, store3) in zip(state1, state3)
    @test length(store1.definitions.inners) == length(store3.definitions.inners)
  end
end

# Test varargs and varkwargs
# ==========================

@testset "Vararg varkwarg" begin
  # super default case
  @traits_test f(args...; kwargs...) = (args, kwargs)
  @traits_show_implementation f

  a, k = f(1,2,3,4, a=4, b=5)
  @test a == (1,2,3,4)
  @test k[:a] == 4
  @test k[:b] == 5
  @test sort(collect(keys(k))) == [:a, :b]
end

# Test complex type dispatch
# ==========================

# Tuple{Tuple{Tuple{Vector{Some}}, Vararg}}

# TODO @testset still interacts even with @traits_test ... unbelievable...
# @testset "complex dispatch" begin
  @traits_test g(a::Tuple{Number, Tuple{Int, Vector}}) = 1
  @traits_show_implementation g
  names(Main, all=true)
  @traits_test function g(a::Tuple{Number, Tuple{Int, T}}) where {T <: Vector, eltype(T)::Type{Int}}
    2
  end
  @traits_show_implementation g

  @test g((1.0, (2, [:a, :b]))) == 1
  @test g((1.0, (2, [1, 2]))) == 2
# end

# TODO add case with Union type
# TODO add case for diagonal dispatch `f(a::T, b::T) where T`
# TODO add case for arg without name

# TODO add case for function in submodule like `@traits_test mysubmodule.myfunc(a) = 4`
# TODO add case for overwriting same target function from different modules like `@traits Target.myfunc(...)` from
# both module A and module B

# TODO test that `@traits f(a, b)` indeed inserts two where variables T1 and T2

# TODO test that traits are normalized correctly
# TODO we currently do not normalize the traits function names
# TODO e.g. using both `Base.IteratorSize(a)` and `IteratorSize(a)` should result in the same trait, but don't do currently


text:
"""

**`Helium.writehe`** -*Function*.

    `Helium.writehe(matrix, source;
                    colNames::Array{String,1} = [""],
                    rowNames::Array{String,1} = [""],
                    supplement::Array{String,2} = ["" ""])` => `file`

Write a matrix (2D array) to a binary .he file, given as an IO argument or
String/FilePaths.jl type representing the file name to write. The binary files
Helium format includes a header that describes the size and the data type of
the matrix, as well as the endianness of the file.
The first 56 bytes contains:
[number of rows, number of columns, data type, endianness,
 if column names exists, if row names exists, number of column of supplement]

"""
function writehe(mat, heFile::String;
                 colNames::Array{String,1} = [""],
                 rowNames::Array{String,1} = [""],
                 supplement::Array{String,2} = ["" ""])

    # Check name ending by .he
    if heFile[end-2:end] != ".he"
       heFile = string(heFile, ".he")
    end

    # Check path
    try
      realpath(heFile[1:end-length(basename(heFile))])
    catch e
      error("Path incorrect, check directory!!!")
    end

    # Check for row and column names
    hascolnames = colNames[1] != ""
    hasrownames = rowNames[1] != ""
    hassupplement = supplement != ["" ""]

    if hassupplement
      numcolsupp = size(supplement)[2]
    else
      numcolsupp = 0
    end

    # Create a dictionary to map size, data type and endianness in integer
    dictCtrl =  Dict("Float64"=>Int64(0xf64),
                     "Int64"=>Int64(0xe64),
                     "Float32"=>Int64(0xf32),
                     "Int32"=>Int64(0xe32),
                     "Bool"=>Int64(0xb00),
                     "Little-endian"=>Int64(0x04030201),
                     "Big-endian"=>Int64(0x01020304));

    dataType = string(typeof(mat[1]))

    # Check endianness
    if ENDIAN_BOM == 0x04030201
        endianness = "Little-endian"
    else
        endianness = "Big-endian"
    end

    # Create the header of the binary file
    matCtrl = [size(mat)[1]; size(mat)[2];
               dictCtrl[dataType]; dictCtrl[endianness];
               hascolnames;hasrownames; numcolsupp]
    matCtrl = convert(Array{Int64}, matCtrl)

    # Write HEADER
    open(heFile, "w") do file
             write(file, matCtrl)
    end

    # Write DATA
    open(heFile, "a") do file
             write(file, mat)
    end

    # Write COL NAMES
    if hascolnames
       open(heFile, "a") do file
                write(file, string(join(colNames, ","),"\n" ))
       end
    end

    # Write ROW NAMES
    if hasrownames
       open(heFile, "a") do file
                write(file, string(join(rowNames, ","),"\n" ))
       end
    end

    # Write SUPPLEMENT
    if hassupplement

       global i = 1

       open(heFile, "a") do file
          while i <= size(supplement)[2]
                 write(file, string(join(supplement[:,i], ","),"\n" ))
                 global i += 1
          end
       end
    end

end