text:
using Distributions, DataFrames, DataFramesMeta, Gadfly

# @Input: a function f
# @Output: a function F whose input is a dataframe with the lengths and
# heights of a piecewise linear function and whose output is 1 if the
# dataframe crossed f and 0 otherwise.
function cross_func(f::Function)
  function F(df::DataFrame)
    x = [0.]
    y = [0.]
    append!(x,cumsum(df[:length]))
    append!(y,cumsum(df[:height]))
    return 1 - prod(y .>= f.(x))
  end
  return F
end

function cross_func(f::Function,y0::Real)
  function F(df::DataFrame)
    x = [0.]
    y = [y0]
    append!(x,cumsum(df[:length]))
    append!(y,y0 .+ cumsum(df[:height]))
    return 1 - prod(y .>= f.(x))
  end
  return F
end

# Tupple({ℓ_n},{ξ_n}) using the notation of the paper
struct ConvexMinorantStableMeander <: Distribution{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  d::Float64
  η::Float64
  δ::Float64
  γ::Float64
  κ::Float64
  Δ::Int64
  mAst::Int64
  ε::Float64
  ConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer,ε::Float64) =
    (β < -1 || β > 1 || (β == -1 && α <= 1) || α <= 0 || α > 2) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") : (0 >= γ || γ >= α) ?
    error("Parameters' requirements unmet: α>γ>0") : (0 >= δ || δ >= d || d >= 2/(α + β*(α <= 1 ? α : α-2))) ?
    error("Parameters' requirements unmet: 1/(αρ)>d>δ>0") : κ < 0 ?
    error("Parameters' requirements unmet: κ≥0") : Δ <= 0 ?
    error("Parameters' requirements unmet: Δ≥0") : mAst < 0 ?
    error("Parameters' requirements unmet: m*≥0") : ε <= 0 ?
    error("Parameters' requirements unmet: ε>0") : new(α, β, β*(α <= 1 ? 1 : (α-2)/α),(1+β*(α <= 1 ? 1 : (α-2)/α))/2 , d,
      etaF(d*(α + β*(α <= 1 ? α : α-2))/2)*(α + β*(α <= 1 ? α : α-2))/2 , δ, γ, κ + Int(ceil(max(2/(α+β*(α <= 1 ? α : α-2)),
      log(2)*2/(3*etaF(d*(α+β*(α <= 1 ? α : α-2))/2)*(α+β*(α <= 1 ? α : α-2)) ) ))), Δ, mAst, ε)
  ConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer) = ConvexMinorantStableMeander(α,β,d,δ,γ,κ,Δ,mAst,.5^16)
  ConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer) = ConvexMinorantStableMeander(α,β,d,δ,γ,κ,Δ,12)
  ConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real) = ConvexMinorantStableMeander(α,β,d,δ,γ,κ,
    Int(ceil(17*log(2)/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))))
  ConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real) = ConvexMinorantStableMeander(α,β,d,δ,γ,4)
  ConvexMinorantStableMeander(α::Real,β::Real) = ConvexMinorantStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α)
  ConvexMinorantStableMeander(α::Real,β::Real,Δ::Integer,ε::Real) = ConvexMinorantStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,4,Δ,12,ε)
  ConvexMinorantStableMeander(α::Real,β::Real,ε::Real) = ConvexMinorantStableMeander(α,β,
    Int(ceil(abs(log(ε/2))/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),ε)
end

import Distributions.params
function params(d::ConvexMinorantStableMeander)
  return (d.α,d.β,d.θ,d.ρ,
          d.d,d.η,d.δ,d.γ,d.κ,d.Δ,d.mAst,
          d.ε)
end

# Useful auxiliary function that reconstructs the faces from the sequences {S_n} and {U_n}
# @Input: Stability parameter and (almost) output of rand(ConvexMinorantStableMeander)
# @Output: DataFrame with heights, lengths, and slopes of the faces of the
# lower and upper convex minorants sorted by slopes and their distance
function faces(α::Real,D::Real,S::Array{Float64},U::Array{Float64})
  ℓ = cumprod(U)
  push!(ℓ,0.)
  for i = length(ℓ):-1:2
      ℓ[i] = ℓ[i-1]-ℓ[i]
  end
  ℓ[1] = 1 - ℓ[1]
  push!(ℓ,ℓ[end])

  H = vcat(S,[D,0.]) .* ℓ .^ (1/α)
  err = H[end-1]

  # Labels: 0 belongs to both estimates, 1 to upper estimate and -1 to lower estimate
  lab = Int.(zeros(length(ℓ)))
  lab[end-1] = 1
  lab[end] = -1
  S = H ./ ℓ

  df = DataFrame(length = ℓ, height = H, slope = S, label = lab)
  sort!(df, order(:slope, rev = false))

  return (df,err)
end

import Distributions.rand
# @Input: ConvexMinorantStableMeander
# @Output: (df,err)=(faces,error)
# Note: We follow the algorithms in the EpsStrongStable paper, with slight modifications
# The initial face's length should be replaced with 0 to get a lower estimate
function rand(d::ConvexMinorantStableMeander)
  ar = d.α*d.ρ
  et1 = 1 /(ar+d.η)
  ra = 1 /d.α
  rar = 1 /ar
  # Conditionally positive stable distribution
  dPos = PositiveStable(d.α,d.β)
  sep = d.d - d.δ
  e2 = 1 /(1 -exp(-sep))
  mom = mellin(dPos,d.γ)
  dega = d.δ*d.γ
  dega2 = -1 /(1 -exp(-dega))
  mAst = d.mAst + Int(ceil( max(0,log(mom)/dega) + rar ))
  # θ sequence
  U = Float64[]
  # Line 1
  x = Inf
  σ = 0
  # Line 2
  lagU = rand(Beta(d.ρ,1),d.Δ)
  lagS = rand(dPos,d.Δ)
  eps = d.ε / prod(lagU)^ra

  # Line 3 (first iteration)
  # Line 4 (Algorithm 3 with unconditional probabilities)
  S = rand(dPos,mAst)
  append!(S,chiStable1(dPos,d.δ,dega,dega2,mom,mAst-1))

  R = Float64[]
  C = [0.]
  F = Float64[]
  # Last value of C
  endC = 0.
  # First iteration
  # Line 5
  while (length(R) == σ) || (length(S) >= length(C))
    t = length(C)
    # Lines 6 and 7
    (C1,F1) = downRW(d.d,rar,d.η,et1,d.κ,x-endC)
    append!(F,F1)
    append!(C,endC .+ C1)
    endC = C[end]
    # Lines 8 and 9
    if reflectedProcess(d.d,d.η,et1,d.κ,x-endC)
      # Line 10
      x = endC + d.κ
      R1 = [maximum(C[t:end])]
      for i = (t-1):-1:(length(R)+1)
        push!(R1,max(R1[end],C[i]))
      end
      append!(R,R1[end:-1:1] .- C[(length(R)+1):t])
    else # Line 11
      # Lines 12 and 13
      (C1,F1) = upRW(d.d,d.η,et1,d.κ,x-endC)
      append!(F,F1)
      append!(C,endC .+ C1)
      endC = C[end]
    end # Line 14
  end # Line 15
  # Line 16
  while length(U) < length(S)
    push!(U,exp(-d.α*F[length(U)+1]))
  end
  # Lines 17, 18 and 19
  σ = 1
  D = exp(R[σ])*( e2*exp(-sep*(length(S)-σ-1)) + sum(S[σ:end] .*
    exp.(-d.d .* (0:(length(S)-σ))) .* (1 .- U[σ:end]) .^ ra) )
  eps /= U[σ]^ra
  push!(lagS,S[σ])
  push!(lagU,U[σ])
  # Line 20
  if D <= eps # Finished!
    # Line 21
    (df,err) = faces(d.α,D,lagS,lagU)
    return (df,err,σ)
  end # Line 22
  # Line 3 (second to last iterations)
  while true
    # Line 4 (Algorithm 3 with conditional probabilities)
    append!(S,chiStable2(dPos,d.δ,dega,dega2,mom,length(S)-σ-1))
    # Line 5
    while (length(R) == σ) || (length(S) >= length(C))
      t = length(C)-1
      # Lines 6 and 7
      (C1,F1) = downRW(d.d,rar,d.η,et1,d.κ,x-endC)
      append!(F,F1)
      append!(C,endC .+ C1)
      endC = C[end]
      # Lines 8 and 9
      if reflectedProcess(d.d,d.η,et1,d.κ,x-endC)
        # Line 10
        x = endC + d.κ
        R1 = [maximum(C[t:end])]
        for i = (t-1):-1:(length(R)+1)
          push!(R1,max(R1[end],C[i]))
        end
        append!(R,R1[end:-1:1] .- C[(length(R)+1):t])
      else # Line 11
        # Lines 12 and 13
        (C1,F1) = upRW(d.d,d.η,et1,d.κ,x-endC)
        append!(F,F1)
        append!(C,endC .+ C1)
        endC = C[end]
      end # Line 14
    end # Line 15
    # Line 16
    while length(U) < length(S)
      push!(U,exp(-d.α*F[length(U)+1]))
    end
    # Lines 17, 18 and 19
    σ += 1
    D = exp(R[σ])*( e2*exp(-sep*(length(S)-σ-1)) + sum(S[σ:end] .*
      exp.(-d.d .* (0:(length(S)-σ))) .* (1 .- U[σ:end]) .^ ra) )
    eps /= U[σ]^ra
    push!(lagS,S[σ])
    push!(lagU,U[σ])
    # Line 24
    if D <= eps # Finished!
      # Line 25
      (df,err) = faces(d.α,D,lagS,lagU)
      return (df,err,σ)
    end # Line 26
  end # Line 27
end

struct CrossConvexMinorantStableMeander <: Distribution{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  d::Float64
  η::Float64
  δ::Float64
  γ::Float64
  κ::Float64
  Δ::Int64
  mAst::Int64
  f::Function
  CrossConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer,f::Function) =
    (β < -1 || β > 1 || (β == -1 && α <= 1) || α <= 0 || α > 2) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") : (0 >= γ || γ >= α) ?
    error("Parameters' requirements unmet: α>γ>0") : (0 >= δ || δ >= d || d >= 2/(α + β*(α <= 1 ? α : α-2))) ?
    error("Parameters' requirements unmet: 1/(αρ)>d>δ>0") : κ < 0 ?
    error("Parameters' requirements unmet: κ≥0") : Δ <= 0 ?
    error("Parameters' requirements unmet: Δ≥0") : mAst < 0 ?
    error("Parameters' requirements unmet: m*≥0") : new(α, β, β*(α <= 1 ? 1 : (α-2)/α),(1+β*(α <= 1 ? 1 : (α-2)/α))/2 , d,
      etaF(d*(α + β*(α <= 1 ? α : α-2))/2)*(α + β*(α <= 1 ? α : α-2))/2 , δ, γ, κ + Int(ceil(max(2/(α+β*(α <= 1 ? α : α-2)),
      log(2)*2/(3*etaF(d*(α+β*(α <= 1 ? α : α-2))/2)*(α+β*(α <= 1 ? α : α-2)) ) ))), Δ, mAst, cross_func(f))
  CrossConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,f::Function) = CrossConvexMinorantStableMeander(α,β,d,δ,γ,κ,Δ,12,f)
  CrossConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,f::Function) = CrossConvexMinorantStableMeander(α,β,d,δ,γ,κ,
    Int(ceil(33*log(2)/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),f)
  CrossConvexMinorantStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,f::Function) = CrossConvexMinorantStableMeander(α,β,d,δ,γ,4,f)
  CrossConvexMinorantStableMeander(α::Real,β::Real,Δ::Integer,mAst::Integer,f::Function) = CrossConvexMinorantStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,4,Δ,mAst,f)
  CrossConvexMinorantStableMeander(α::Real,β::Real,Δ::Integer,f::Function) = CrossConvexMinorantStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,4,Δ,12,f)
  CrossConvexMinorantStableMeander(α::Real,β::Real,f::Function) = CrossConvexMinorantStableMeander(α,β,Int(ceil(33*log(2)/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),f)
end

# @Input: CrossConvexMinorantStableMeander
# @Output: The triple (B,df,err) where (df,err) is the output of the function 'faces' and
# B is a boolean indicating a crossing of the function f or not
# Note: We follow the algorithms in the EpsStrongStable paper, with slight modifications
# The initial face's length should be replaced with 0 to get a lower estimate
function rand(d::CrossConvexMinorantStableMeander)
  ar = d.α*d.ρ
  et1 = 1 /(ar+d.η)
  ra = 1 /d.α
  rar = 1 /ar
  # Conditionally positive stable distribution
  dPos = PositiveStable(d.α,d.β)
  sep = d.d - d.δ
  e2 = 1 /(1 -exp(-sep))
  mom = mellin(dPos,d.γ)
  dega = d.δ*d.γ
  dega2 = -1 /(1 -exp(-dega))
  mAst = d.mAst + Int(ceil( max(0,log(mom)/dega) + rar ))
  # θ sequence
  U = Float64[]
  # Line 1
  x = Inf
  σ = 0
  # Line 2
  lagU = rand(Beta(d.ρ,1),d.Δ)
  lagS = rand(dPos,d.Δ)

  # Line 3 (first iteration)
  # Line 4 (Algorithm 3 with unconditional probabilities)
  S = rand(dPos,mAst)
  append!(S,chiStable1(dPos,d.δ,dega,dega2,mom,mAst-1))

  R = Float64[]
  C = [0.]
  F = Float64[]
  # Last value of C
  endC = 0.
  # First iteration
  # Line 5
  while (length(R) == σ) || (length(S) >= length(C))
    t = length(C)
    # Lines 6 and 7
    (C1,F1) = downRW(d.d,rar,d.η,et1,d.κ,x-endC)
    append!(F,F1)
    append!(C,endC .+ C1)
    endC = C[end]
    # Lines 8 and 9
    if reflectedProcess(d.d,d.η,et1,d.κ,x-endC)
      # Line 10
      x = endC + d.κ
      R1 = [maximum(C[t:end])]
      for i = (t-1):-1:(length(R)+1)
        push!(R1,max(R1[end],C[i]))
      end
      append!(R,R1[end:-1:1] .- C[(length(R)+1):t])
    else # Line 11
      # Lines 12 and 13
      (C1,F1) = upRW(d.d,d.η,et1,d.κ,x-endC)
      append!(F,F1)
      append!(C,endC .+ C1)
      endC = C[end]
    end # Line 14
  end # Line 15
  # Line 16
  while length(U) < length(S)
    push!(U,exp(-d.α*F[length(U)+1]))
  end
  # Lines 17, 18 and 19
  σ = 1
  D = exp(R[σ])*( e2*exp(-sep*(length(S)-σ-1)) + sum(S[σ:end] .*
    exp.(-d.d .* (0:(length(S)-σ))) .* (1 .- U[σ:end]) .^ ra) )
  # Line 20
  # Line 21
  push!(lagS,S[1])
  push!(lagU,U[1])
  # Line 25
  (df,err) = faces(d.α,D,lagS,lagU)
  dfL = @linq df |>
      where(:label .< 1) |>
      select(:length,:height)
  dfU = @linq df |>
      where(:label .> -1) |>
      select(:length,:height)
  aux = d.f(dfL)
  if aux == d.f(dfU) # Crossing / non-crossing detection!
    return (aux,df,err)
  end # Line 22
  # Line 3 (second to last iterations)
  while true
    # Line 4 (Algorithm 3 with conditional probabilities)
    append!(S,chiStable2(dPos,d.δ,dega,dega2,mom,length(S)-σ-1))
    # Line 5
    while (length(R) == σ) || (length(S) >= length(C))
      t = length(C)-1
      # Lines 6 and 7
      (C1,F1) = downRW(d.d,rar,d.η,et1,d.κ,x-endC)
      append!(F,F1)
      append!(C,endC .+ C1)
      endC = C[end]
      # Lines 8 and 9
      if reflectedProcess(d.d,d.η,et1,d.κ,x-endC)
        # Line 10
        x = endC + d.κ
        R1 = [maximum(C[t:end])]
        for i = (t-1):-1:(length(R)+1)
          push!(R1,max(R1[end],C[i]))
        end
        append!(R,R1[end:-1:1] .- C[(length(R)+1):t])
      else # Line 11
        # Lines 12 and 13
        (C1,F1) = upRW(d.d,d.η,et1,d.κ,x-endC)
        append!(F,F1)
        append!(C,endC .+ C1)
        endC = C[end]
      end # Line 14
    end # Line 15
    # Line 16
    while length(U) < length(S)
      push!(U,exp(-d.α*F[length(U)+1]))
    end
    # Lines 17, 18 and 19
    σ += 1
    D = exp(R[σ])*( e2*exp(-sep*(length(S)-σ-1)) + sum(S[σ:end] .*
      exp.(-d.d .* (0:(length(S)-σ))) .* (1 .- U[σ:end]) .^ ra) )
    # Line 24
    push!(lagS,S[σ])
    push!(lagU,U[σ])
    # Line 25
    (df,err) = faces(d.α,D,lagS,lagU)
    dfL = @linq df |>
        where(:label .< 1) |>
        select(:length,:height)
    dfU = @linq df |>
        where(:label .> -1) |>
        select(:length,:height)
    aux = d.f(dfL)
    if aux == d.f(dfU)
      return (aux,df,err)
    end # Line 26
  end # Line 27
end

struct ConvexMinorantWeaklyStable <: Distribution{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  # Scale, drift and time horizon
  δ::Float64
  μ::Float64
  T::Float64
  # For the rest of the parameters, we have left/right versions
  d_l::Float64
  d_r::Float64
  η_l::Float64
  η_r::Float64
  δ_l::Float64
  δ_r::Float64
  γ_l::Float64
  γ_r::Float64
  κ_l::Float64
  κ_r::Float64
  Δ_l::Int64
  Δ_r::Int64
  mAst_l::Int64
  mAst_r::Int64
  ε::Float64
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,mAst_l::Integer,mAst_r::Integer,ε::Float64) =
    (β < -1 || β > 1 || α <= 0 || α > 2) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]") :
    (δ <= 0) ?
    error("Parameters' requirements unmet: δ>0") :
    (T <= 0) ?
    error("Parameters' requirements unmet: T>0") :
    (β == -1 && α <= 1) ? (0 >= γ_l || γ_l >= α) ?
      error("Parameters' requirements unmet: α>γ_l>0") : (0 >= δ_l || δ_l >= d_l || d_l >= 1/α) ?
      error("Parameters' requirements unmet: 1/α>d_l>δ_l>0") : κ_l < 0 ?
      error("Parameters' requirements unmet: κ_l≥0") : Δ_l <= 0 ?
      error("Parameters' requirements unmet: Δ_l≥0") : mAst_l < 0 ?
      error("Parameters' requirements unmet: m_l*≥0") : ε <= 0 ?
      error("Parameters' requirements unmet: ε>0") : new(α,-1.,-1.,0.,δ,μ,T,d_l,d_r,etaF(d_l*α)*α,0.,
        δ_l,δ_r,γ_l,γ_r,κ_l + Int(ceil(max(1/α,log(2)*2/(3*etaF(d_l*α)*α)))),κ_r,Δ_l,Δ_r,mAst_l,mAst_r,ε) :
    (β == 1 && α <= 1) ? (0 >= γ_r || γ_r >= α) ?
      error("Parameters' requirements unmet: α>γ_r>0") : (0 >= δ_r || δ_r >= d_r || d_r >= 1/α) ?
      error("Parameters' requirements unmet: 1/α>d_r>δ_r>0") : κ_r < 0 ?
      error("Parameters' requirements unmet: κ_r≥0") : Δ_r <= 0 ?
      error("Parameters' requirements unmet: Δ_r≥0") : mAst_r < 0 ?
      error("Parameters' requirements unmet: m_r*≥0") : ε <= 0 ?
      error("Parameters' requirements unmet: ε>0") : new(α,1.,1.,1.,δ,μ,T,d_l,d_r,0.,etaF(d_r*α)*α,
        δ_l,δ_r,γ_l,γ_r,κ_l,κ_r + Int(ceil(max(1/α,log(2)*2/(3*etaF(d_r*α)*α)))),Δ_l,Δ_r,mAst_l,mAst_r,ε) :
    (0 >= γ_l || γ_l >= α || 0 >= γ_r || γ_r >= α) ?
    error("Parameters' requirements unmet: α>γ_l>0 and α>γ_r>0") : (0 >= δ_l || δ_l >= d_l || 0 >= δ_r || δ_r >= d_r || d_l >= 2/(α - β*(α <= 1 ? α : α-2)) || d_r >= 2/(α + β*(α <= 1 ? α : α-2))) ?
    error("Parameters' requirements unmet: 1/(α(1-ρ))>d_l>δ_l>0 and 1/(αρ)>dr>δr>0") : (κ_l < 0 || κ_r < 0) ?
    error("Parameters' requirements unmet: κ_l,κ_r≥0") : (Δ_l <= 0 || Δ_r <= 0) ?
    error("Parameters' requirements unmet: Δ_l,Δr≥0") : (mAst_l < 0 || mAst_r < 0) ?
    error("Parameters' requirements unmet: m_l*,m_r*≥0") : ε <= 0 ?
    error("Parameters' requirements unmet: ε>0") : new(α, β, β*(α <= 1 ? 1 : (α-2)/α),(1+β*(α <= 1 ? 1 : (α-2)/α))/2 , δ, μ, T, d_l, d_r,
      etaF(d_l*(α - β*(α <= 1 ? α : α-2))/2)*(α - β*(α <= 1 ? α : α-2))/2, etaF(d_r*(α + β*(α <= 1 ? α : α-2))/2)*(α + β*(α <= 1 ? α : α-2))/2,
      δ_l, δ_r, γ_l, γ_r,
      κ_l + Int(ceil(max(2/(α-β*(α <= 1 ? α : α-2)),log(2)*2/(3*etaF(d_l*(α-β*(α <= 1 ? α : α-2))/2)*(α-β*(α <= 1 ? α : α-2)) ) ))),
      κ_r + Int(ceil(max(2/(α+β*(α <= 1 ? α : α-2)),log(2)*2/(3*etaF(d_r*(α+β*(α <= 1 ? α : α-2))/2)*(α+β*(α <= 1 ? α : α-2)) ) ))),
      Δ_l, Δ_r, mAst_l, mAst_r, ε)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,mAst_l::Integer,mAst_r::Integer) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r,.5^16)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,12,12)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,
    Int(ceil(17*log(2)/log(gamma(1+1/α+(1-β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1-β*(α <= 1 ? 1 : (α-2)/α))/2))))),
    Int(ceil(17*log(2)/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))))
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,4,4)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,Δ_l::Integer,Δ_r::Integer,ε::Real) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α,4,4,Δ_l,Δ_r,12,12,ε)
  ConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,ε::Real) = ConvexMinorantWeaklyStable(α,β,δ,μ,T,
    Int(ceil(abs(log(ε/2))/log(gamma(1+1/α+(1-β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1-β*(α <= 1 ? 1 : (α-2)/α))/2))))),
    Int(ceil(abs(log(ε/2))/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),ε)
end

# @Input: ConvexMinorantStableMeander
# @Output: (df,err_left,err_right)=(faces,error) where err_left and err_right
# are the errors produced by the pre- and post-minumum convex minorants (before tilting)
function rand(d::ConvexMinorantWeaklyStable)
  ra = 1/d.α
  ϵ = d.ε / (d.δ * d.T^ra)

  if d.α <= 1 && d.β == -1
    (df,err) = rand(ConvexMinorantStableMeander(d.α,1,d.d_l,d.δ_l,d.γ_l,d.κ_l,d.Δ_l,d.mAst_l,ϵ))
    sort!(df, order(:slope, rev = true))
    df[:length] = df[:length] .* d.T
    df[:height] = df[:height] .* (-d.δ * d.T^ra) .+ df[:length] .* d.μ
    df[:slope] = df[:slope] .* (-d.δ * d.T^(ra-1)) .+ d.μ
    return (df,err * d.δ * d.T^ra,0.)
  elseif d.α <= 1 && d.β == 1
    (df,err) = rand(ConvexMinorantStableMeander(d.α,1,d.d_r,d.δ_r,d.γ_r,d.κ_r,d.Δ_r,d.mAst_r,ϵ))
    df[:length] = df[:length] .* d.T
    df[:height] = df[:height] .* (d.δ * d.T^ra) .+ df[:length] .* d.μ
    df[:slope] = df[:slope] .* (d.δ * d.T^(ra-1)) .+ d.μ
    return (df,0,err * d.δ * d.T^ra)
  end
  τ = rand(Beta(1-d.ρ,d.ρ))
  T_l = d.T * τ
  T_r = d.T * (1-τ)

  ar_l, ar_r = d.α*(1-d.ρ), d.α*d.ρ
  et1_l, et1_r = 1/(ar_l+d.η_l), 1/(ar_r+d.η_r)
  rar_l, rar_r = 1/ar_l, 1/ar_r

  # Conditionally positive stable distribution
  dPos_l, dPos_r = PositiveStable(d.α,-d.β), PositiveStable(d.α,d.β)
  sep_l, sep_r = d.d_l - d.δ_l, d.d_r - d.δ_r
  e2_l, e2_r = 1 /(1 -exp(-sep_l)), 1 /(1 -exp(-sep_r))
  mom_l, mom_r = mellin(dPos_l,d.γ_l), mellin(dPos_r,d.γ_r)
  dega_l, dega_r = d.δ_l*d.γ_l, d.δ_r*d.γ_r
  dega2_l, dega2_r = -1 /(1 -exp(-dega_l)), -1 /(1 -exp(-dega_r))
  mAst_l, mAst_r = d.mAst_l + Int(ceil( max(0,log(mom_l)/dega_l) + rar_l )), d.mAst_r + Int(ceil( max(0,log(mom_r)/dega_r) + rar_r ))
  # θ sequence
  U_l, U_r = Float64[], Float64[]
  # Line 1
  x_l, x_r = Inf, Inf
  σ, σ = 0, 0
  # Line 2
  lagU_l, lagU_r = rand(Beta(1-d.ρ,1),d.Δ_l), rand(Beta(d.ρ,1),d.Δ_r)
  lagS_l, lagS_r = rand(dPos_l,d.Δ_l), rand(dPos_r,d.Δ_r)
  eps_l, eps_r = 1 / (τ * prod(lagU_l))^ra, 1 / ((1-τ) * prod(lagU_r))^ra

  # Line 3 (first iteration)
  # Line 4 (Algorithm 3 with unconditional probabilities)
  S_l, S_r = rand(dPos_l,mAst_l), rand(dPos_r,mAst_r)
  append!(S_l,chiStable1(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,mAst_l-1))
  append!(S_r,chiStable1(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,mAst_r-1))

  R_l, R_r = Float64[], Float64[]
  C_l, C_r = [0.], [0.]
  F_l, F_r = Float64[], Float64[]
  # Last value of C
  endC_l, endC_r = 0., 0.
  # First iteration
  # Line 5
  while (length(R_l) == σ) || (length(S_l) >= length(C_l))
    t_l = length(C_l)
    # Lines 6 and 7
    (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
    append!(F_l,F1_l)
    append!(C_l,endC_l .+ C1_l)
    endC_l = C_l[end]
    # Lines 8 and 9
    if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      # Line 10
      x_l = endC_l + d.κ_l
      R1_l = [maximum(C_l[t_l:end])]
      for i = (t_l-1):-1:(length(R_l)+1)
        push!(R1_l,max(R1_l[end],C_l[i]))
      end
      append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
    else # Line 11
      # Lines 12 and 13
      (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      append!(F_l,F1_l)
      append!(C_l,endC_l .+ C1_l)
      endC_l = C_l[end]
    end # Line 14
  end # Line 15 - left
  while (length(R_r) == σ) || (length(S_r) >= length(C_r))
    t_r = length(C_r)
    # Lines 6 and 7
    (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
    append!(F_r,F1_r)
    append!(C_r,endC_r .+ C1_r)
    endC_r = C_r[end]
    # Lines 8 and 9
    if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      # Line 10
      x_r = endC_r + d.κ_r
      R1_r = [maximum(C_r[t_r:end])]
      for i = (t_r-1):-1:(length(R_r)+1)
        push!(R1_r,max(R1_r[end],C_r[i]))
      end
      append!(R_r,R1_r[end:-1:1] .- C_r[(length(R_r)+1):t_r])
    else # Line 11
      # Lines 12 and 13
      (C1_r,F1_r) = upRW(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      append!(F_r,F1_r)
      append!(C_r,endC_r .+ C1_r)
      endC_r = C_r[end]
    end # Line 14
  end # Line 15 - right
  # Line 16
  while length(U_l) < length(S_l)
    push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
  end
  while length(U_r) < length(S_r)
    push!(U_r,exp(-d.α*F_r[length(U_r)+1]))
  end
  # Lines 17, 18 and 19
  σ, σ = 1, 1
  D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
      exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )
  D_r = exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r)-σ-1)) + sum(S_r[σ:end] .*
      exp.(-d.d_r .* (0:(length(S_r)-σ))) .* (1 .- U_r[σ:end]) .^ ra) )
  eps_l /= U_l[σ]^ra
  eps_r /= U_r[σ]^ra

  push!(lagS_l,S_l[σ])
  push!(lagS_r,S_r[σ])
  push!(lagU_l,U_l[σ])
  push!(lagU_r,U_r[σ])
  # Line 20
  if D_l/eps_l + D_r/eps_r <= ϵ # Finished!
    # Line 21
    (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
    (df_r,err_r) = faces(d.α,D_r,lagS_r,lagU_r)

    df_l = @linq df_l |>
         where(:label .> -1) |>
         select(:length,:height,:slope,:label)
    df_l[:label] = 0
    sort!(df_l, order(:slope, rev = true))
    df_l[:length] = df_l[:length] .* T_l
    df_l[:height] = df_l[:height] .* (-d.δ * T_l^ra) .+ df_l[:length] .* d.μ
    df_l[:slope] = df_l[:slope] .* (-d.δ * T_l^(ra-1)) .+ d.μ
    df_r[:length] = df_r[:length] .* T_r
    df_r[:height] = df_r[:height] .* (d.δ * T_r^ra) .+ df_r[:length] .* d.μ
    df_r[:slope] = df_r[:slope] .* (d.δ * T_r^(ra-1)) .+ d.μ

    append!(df_l,df_r)
    return (df_l,err_l * d.δ * T_l^ra,err_r * d.δ * T_r^ra)
  end # Line 22
  # Line 3 (second to last iterations)
  while true
    # Line 4 (Algorithm 3 with conditional probabilities)
    append!(S_l,chiStable2(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,length(S_l)-σ-1))
    append!(S_r,chiStable2(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,length(S_r)-σ-1))
    # Line 5
    while (length(R_l) == σ) || (length(S_l) >= length(C_l))
      t_l = length(C_l)-1
      # Lines 6 and 7
      (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      append!(F_l,F1_l)
      append!(C_l,endC_l .+ C1_l)
      endC_l = C_l[end]
      # Lines 8 and 9
      if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        # Line 10
        x_l = endC_l + d.κ_l
        R1_l = [maximum(C_l[t_l:end])]
        for i = (t_l-1):-1:(length(R_l)+1)
          push!(R1_l,max(R1_l[end],C_l[i]))
        end
        append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
      else # Line 11
        # Lines 12 and 13
        (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        append!(F_l,F1_l)
        append!(C_l,endC_l .+ C1_l)
        endC_l = C_l[end]
      end # Line 14
    end # Line 15
    while (length(R_r) == σ) || (length(S_r) >= length(C_r))
      t_r = length(C_r)-1
      # Lines 6 and 7
      (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      append!(F_r,F1_r)
      append!(C_r,endC_r .+ C1_r)
      endC_r = C_r[end]
      # Lines 8 and 9
      if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
        # Line 10
        x_r = endC_r + d.κ_r
        R1_r = [maximum(C_r[t_r:end])]
        for i = (t_r-1):-1:(length(R_r)+1)
          push!(R1_r,max(R1_r[end],C_r[i]))
        end
        append!(R_r,R1_r[end:-1:1] .- C_r[(length(R_r)+1):t_r])
      else # Line 11
        # Lines 12 and 13
        (C1_r,F1_r) = upRW(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
        append!(F_r,F1_r)
        append!(C_r,endC_r .+ C1_r)
        endC_r = C_r[end]
      end # Line 14
    end # Line 15
    # Line 16
    while length(U_l) < length(S_l)
      push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
    end
    while length(U_r) < length(S_r)
      push!(U_r,exp(-d.α*F_r[length(U_r)+1]))
    end
    # Lines 17, 18 and 19
    σ += 1
    σ += 1
    D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
        exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )
    D_r = exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r)-σ-1)) + sum(S_r[σ:end] .*
        exp.(-d.d_r .* (0:(length(S_r)-σ))) .* (1 .- U_r[σ:end]) .^ ra) )
    eps_l /= U_l[σ]^ra
    eps_r /= U_r[σ]^ra
    push!(lagS_l,S_l[σ])
    push!(lagS_r,S_r[σ])
    push!(lagU_l,U_l[σ])
    push!(lagU_r,U_r[σ])
    # Line 24
    if D_l/eps_l + D_r/eps_r <= ϵ # Finished!
      # Line 25
      (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
      (df_r,err_r) = faces(d.α,D_r,lagS_r,lagU_r)

      df_l = @linq df_l |>
           where(:label .> -1) |>
           select(:length,:height,:slope,:label)
      df_l[:label] = 0
      sort!(df_l, order(:slope, rev = true))
      df_l[:length] = df_l[:length] .* T_l
      df_l[:height] = df_l[:height] .* (-d.δ * T_l^ra) .+ df_l[:length] .* d.μ
      df_l[:slope] = df_l[:slope] .* (-d.δ * T_l^(ra-1)) .+ d.μ
      df_r[:length] = df_r[:length] .* T_r
      df_r[:height] = df_r[:height] .* (d.δ * T_r^ra) .+ df_r[:length] .* d.μ
      df_r[:slope] = df_r[:slope] .* (d.δ * T_r^(ra-1)) .+ d.μ

      append!(df_l,df_r)
      return (df_l,err_l * d.δ * T_l^ra,err_r * d.δ * T_r^ra)
    end # Line 26
  end # Line 27
end

struct CrossConvexMinorantWeaklyStable <: Distribution{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  # Scale, drift and time horizon
  δ::Float64
  μ::Float64
  T::Float64
  # For the rest of the parameters, we have left/right versions
  d_l::Float64
  d_r::Float64
  η_l::Float64
  η_r::Float64
  δ_l::Float64
  δ_r::Float64
  γ_l::Float64
  γ_r::Float64
  κ_l::Float64
  κ_r::Float64
  Δ_l::Int64
  Δ_r::Int64
  mAst_l::Int64
  mAst_r::Int64
  f::Function
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,mAst_l::Integer,mAst_r::Integer,f::Function) =
    (β < -1 || β > 1 || α <= 0 || α > 2) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]") :
    (δ <= 0) ?
    error("Parameters' requirements unmet: δ>0") :
    (T <= 0) ?
    error("Parameters' requirements unmet: T>0") :
    (β == -1 && α <= 1) ? (0 >= γ_l || γ_l >= α) ?
      error("Parameters' requirements unmet: α>γ_l>0") : (0 >= δ_l || δ_l >= d_l || d_l >= 1/α) ?
      error("Parameters' requirements unmet: 1/α>d_l>δ_l>0") : κ_l < 0 ?
      error("Parameters' requirements unmet: κ_l≥0") : Δ_l <= 0 ?
      error("Parameters' requirements unmet: Δ_l≥0") : mAst_l < 0 ?
      error("Parameters' requirements unmet: m_l*≥0") : new(α,-1.,-1.,0.,δ,μ,T,d_l,d_r,etaF(d_l*α)*α,0.,
        δ_l,δ_r,γ_l,γ_r,κ_l + Int(ceil(max(1/α,log(2)*2/(3*etaF(d_l*α)*α)))),κ_r,Δ_l,Δ_r,mAst_l,mAst_r,f) :
    (β == 1 && α <= 1) ? (0 >= γ_r || γ_r >= α) ?
      error("Parameters' requirements unmet: α>γ_r>0") : (0 >= δ_r || δ_r >= d_r || d_r >= 1/α) ?
      error("Parameters' requirements unmet: 1/α>d_r>δ_r>0") : κ_r < 0 ?
      error("Parameters' requirements unmet: κ_r≥0") : Δ_r <= 0 ?
      error("Parameters' requirements unmet: Δ_r≥0") : mAst_r < 0 ?
      error("Parameters' requirements unmet: m_r*≥0") : new(α,1.,1.,1.,δ,μ,T,d_l,d_r,0.,etaF(d_r*α)*α,
        δ_l,δ_r,γ_l,γ_r,κ_l,κ_r + Int(ceil(max(1/α,log(2)*2/(3*etaF(d_r*α)*α)))),Δ_l,Δ_r,mAst_l,mAst_r,f) :
    (0 >= γ_l || γ_l >= α || 0 >= γ_r || γ_r >= α) ?
    error("Parameters' requirements unmet: α>γ_l>0 and α>γ_r>0") : (0 >= δ_l || δ_l >= d_l || 0 >= δ_r || δ_r >= d_r || d_l >= 2/(α - β*(α <= 1 ? α : α-2)) || d_r >= 2/(α + β*(α <= 1 ? α : α-2))) ?
    error("Parameters' requirements unmet: 1/(α(1-ρ))>d_l>δ_l>0 and 1/(αρ)>dr>δr>0") : (κ_l < 0 || κ_r < 0) ?
    error("Parameters' requirements unmet: κ_l,κ_r≥0") : (Δ_l <= 0 || Δ_r <= 0) ?
    error("Parameters' requirements unmet: Δ_l,Δr≥0") : (mAst_l < 0 || mAst_r < 0) ?
    error("Parameters' requirements unmet: m_l*,m_r*≥0") : new(α, β, β*(α <= 1 ? 1 : (α-2)/α),(1+β*(α <= 1 ? 1 : (α-2)/α))/2, δ, μ, T, d_l, d_r,
      etaF(d_l*(α - β*(α <= 1 ? α : α-2))/2)*(α - β*(α <= 1 ? α : α-2))/2, etaF(d_r*(α + β*(α <= 1 ? α : α-2))/2)*(α + β*(α <= 1 ? α : α-2))/2,
      δ_l, δ_r, γ_l, γ_r,
      κ_l + Int(ceil(max(2/(α-β*(α <= 1 ? α : α-2)),log(2)*2/(3*etaF(d_l*(α-β*(α <= 1 ? α : α-2))/2)*(α-β*(α <= 1 ? α : α-2)) ) ))),
      κ_r + Int(ceil(max(2/(α+β*(α <= 1 ? α : α-2)),log(2)*2/(3*etaF(d_r*(α+β*(α <= 1 ? α : α-2))/2)*(α+β*(α <= 1 ? α : α-2)) ) ))),
      Δ_l, Δ_r, mAst_l, mAst_r, f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,12,12,f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,
    Int(ceil(17*log(2)/log(gamma(1+1/α+(1-β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1-β*(α <= 1 ? 1 : (α-2)/α))/2))))),
    Int(ceil(17*log(2)/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,d_l,d_r,δ_l,δ_r,γ_l,γ_r,4,4,f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α,f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,Δ_l::Integer,Δ_r::Integer,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α,4,4,Δ_l,Δ_r,12,12,f)
  CrossConvexMinorantWeaklyStable(α::Real,β::Real,δ::Real,μ::Real,T::Real,f::Function) = CrossConvexMinorantWeaklyStable(α,β,δ,μ,T,
    Int(ceil(abs(log(.5^16))/log(gamma(1+1/α+(1-β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1-β*(α <= 1 ? 1 : (α-2)/α))/2))))),
    Int(ceil(abs(log(.5^16))/log(gamma(1+1/α+(1+β*(α <= 1 ? 1 : (α-2)/α))/2)/(gamma(1+1/α)*gamma(1+(1+β*(α <= 1 ? 1 : (α-2)/α))/2))))),f)
end

# @Input: CrossConvexMinorantStableMeander
# @Output: The triple (B,df,err_left,err_right) where
# B is a boolean indicating a crossing of the function f or not
function rand(d::CrossConvexMinorantWeaklyStable)
  ra = 1/d.α

  if d.α <= 1 && d.β == -1
    sc_l = d.δ * d.T^ra
    sc1_l = d.δ * d.T^(ra-1)

    F = cross_func(d.f)

    ar_l = d.α*(1-d.ρ)
    et1_l = 1/(ar_l+d.η_l)
    rar_l = 1/ar_l

    # Conditionally positive stable distribution
    dPos_l = PositiveStable(d.α,1)
    sep_l = d.d_l - d.δ_l
    e2_l = 1 /(1 -exp(-sep_l))
    mom_l = mellin(dPos_l,d.γ_l)
    dega_l = d.δ_l*d.γ_l
    dega2_l = -1 /(1 -exp(-dega_l))
    mAst_l = d.mAst_l + Int(ceil( max(0,log(mom_l)/dega_l) + rar_l ))
    # θ sequence
    U_l = Float64[]
    # Line 1
    x_l = Inf
    σ = 0
    # Line 2
    lagU_l = rand(Uniform(),d.Δ_l)
    lagS_l = rand(dPos_l,d.Δ_l)

    # Line 3 (first iteration)
    # Line 4 (Algorithm 3 with unconditional probabilities)
    S_l = rand(dPos_l,mAst_l)
    append!(S_l,chiStable1(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,mAst_l-1))

    R_l = Float64[]
    C_l = [0.]
    F_l = Float64[]
    # Last value of C
    endC_l = 0.
    # First iteration
    # Line 5
    while (length(R_l) == σ) || (length(S_l) >= length(C_l))
      t_l = length(C_l)
      # Lines 6 and 7
      (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      append!(F_l,F1_l)
      append!(C_l,endC_l .+ C1_l)
      endC_l = C_l[end]
      # Lines 8 and 9
      if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        # Line 10
        x_l = endC_l + d.κ_l
        R1_l = [maximum(C_l[t_l:end])]
        for i = (t_l-1):-1:(length(R_l)+1)
          push!(R1_l,max(R1_l[end],C_l[i]))
        end
        append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
      else # Line 11
        # Lines 12 and 13
        (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        append!(F_l,F1_l)
        append!(C_l,endC_l .+ C1_l)
        endC_l = C_l[end]
      end # Line 14
    end # Line 15 - left
    # Line 16
    while length(U_l) < length(S_l)
      push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
    end
    # Lines 17, 18 and 19
    σ = 1
    D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
        exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )

    push!(lagS_l,S_l[σ])
    push!(lagU_l,U_l[σ])
    # Line 20
    # Construct the dataframes
    (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
    df_l = @linq df_l |>
         where(:label .> -1) |>
         select(:length,:height)
    # Rescale and add the drift
    sort!(df_l, order(:slope, rev = true))
    df_l[:length] = df_l[:length] .* d.T
    df_l[:height] = df_l[:height]  .* (-sc_l) .+ df_l[:length] .* d.μ
    df_l[:slope] = df_l[:slope] .* (-sc1_l) .+ d.μ
    # Construct (the dataframes with the faces of) the sandwiching convex functions
    err_l *= sc_l
    # Check for crossing:
    aux = F(df_l)
    G = cross_func(d.f,-err_l)
    if aux == G(df_l) # Crossing / non-crossing detection!
      return (aux,df_l,err_l,0.)
    end # Line 22
    # Line 3 (second to last iterations)
    while true
      # Line 4 (Algorithm 3 with conditional probabilities)
      append!(S_l,chiStable2(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,length(S_l)-σ-1))
      # Line 5
      while (length(R_l) == σ) || (length(S_l) >= length(C_l))
        t_l = length(C_l)-1
        # Lines 6 and 7
        (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        append!(F_l,F1_l)
        append!(C_l,endC_l .+ C1_l)
        endC_l = C_l[end]
        # Lines 8 and 9
        if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
          # Line 10
          x_l = endC_l + d.κ_l
          R1_l = [maximum(C_l[t_l:end])]
          for i = (t_l-1):-1:(length(R_l)+1)
            push!(R1_l,max(R1_l[end],C_l[i]))
          end
          append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
        else # Line 11
          # Lines 12 and 13
          (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
          append!(F_l,F1_l)
          append!(C_l,endC_l .+ C1_l)
          endC_l = C_l[end]
        end # Line 14
      end # Line 15
      # Line 16
      while length(U_l) < length(S_l)
        push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
      end
      # Lines 17, 18 and 19
      σ += 1
      D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
          exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )
      push!(lagS_l,S_l[σ])
      push!(lagU_l,U_l[σ])
      # Line 24
      # Construct the dataframes
      (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
      # Rescale and add the drift
      sort!(df_l, order(:slope, rev = true))
      df_l[:length] = df_l[:length] .* d.T
      df_l[:height] = df_l[:height]  .* (-sc_l) .+ df_l[:length] .* d.μ
      df_l[:slope] = df_l[:slope] .* (-sc1_l) .+ d.μ
      # Construct (the dataframes with the faces of) the sandwiching convex functions
      df = @linq df_l |>
          where(:label .> -1) |>
          select(:length,:height)
      err_l *= sc_l
      # Check for crossing:
      aux = F(df)
      G = cross_func(d.f,-err_l)
      if aux == G(df) # Crossing / non-crossing detection!
        return (aux,df_l,err_l,0.)
      end # Line 26
    end # Line 27
  elseif d.α <= 1 && d.β == 1
    (aux,df,err) = rand(CrossConvexMinorantStableMeander(d.α,1,d.d_r,d.δ_r,d.γ_r,d.κ_r,d.Δ_r,d.mAst_r,d.f))
    df[:length] = df[:length] .* d.T
    df[:height] = df[:height] .* (d.δ * d.T^ra) .+ df[:length] .* d.μ
    df[:slope] = df[:slope] .* (d.δ * d.T^(ra-1)) .+ d.μ
    return (aux,df,0,err * d.δ * d.T^ra)
  end
  τ = rand(Beta(1-d.ρ,d.ρ))
  T_l = d.T * τ
  T_r = d.T * (1-τ)

  sc_l = d.δ * T_l^ra
  sc_r = d.δ * T_r^ra
  sc1_l = d.δ * T_l^(ra-1)
  sc1_r = d.δ * T_r^(ra-1)

  F = cross_func(d.f)

  ar_l, ar_r = d.α*(1-d.ρ), d.α*d.ρ
  et1_l, et1_r = 1/(ar_l+d.η_l), 1/(ar_r+d.η_r)
  rar_l, rar_r = 1/ar_l, 1/ar_r

  # Conditionally positive stable distribution
  dPos_l, dPos_r = PositiveStable(d.α,-d.β), PositiveStable(d.α,d.β)
  sep_l, sep_r = d.d_l - d.δ_l, d.d_r - d.δ_r
  e2_l, e2_r = 1 /(1 -exp(-sep_l)), 1 /(1 -exp(-sep_r))
  mom_l, mom_r = mellin(dPos_l,d.γ_l), mellin(dPos_r,d.γ_r)
  dega_l, dega_r = d.δ_l*d.γ_l, d.δ_r*d.γ_r
  dega2_l, dega2_r = -1 /(1 -exp(-dega_l)), -1 /(1 -exp(-dega_r))
  mAst_l, mAst_r = d.mAst_l + Int(ceil( max(0,log(mom_l)/dega_l) + rar_l )), d.mAst_r + Int(ceil( max(0,log(mom_r)/dega_r) + rar_r ))
  # θ sequence
  U_l, U_r = Float64[], Float64[]
  # Line 1
  x_l, x_r = Inf, Inf
  σ = 0
  # Line 2
  lagU_l, lagU_r = rand(Beta(1-d.ρ,1),d.Δ_l), rand(Beta(d.ρ,1),d.Δ_r)
  lagS_l, lagS_r = rand(dPos_l,d.Δ_l), rand(dPos_r,d.Δ_r)

  # Line 3 (first iteration)
  # Line 4 (Algorithm 3 with unconditional probabilities)
  S_l, S_r = rand(dPos_l,mAst_l), rand(dPos_r,mAst_r)
  append!(S_l,chiStable1(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,mAst_l-1))
  append!(S_r,chiStable1(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,mAst_r-1))

  R_l, R_r = Float64[], Float64[]
  C_l, C_r = [0.], [0.]
  F_l, F_r = Float64[], Float64[]
  # Last value of C
  endC_l, endC_r = 0., 0.
  # First iteration
  # Line 5
  while (length(R_l) == σ) || (length(S_l) >= length(C_l))
    t_l = length(C_l)
    # Lines 6 and 7
    (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
    append!(F_l,F1_l)
    append!(C_l,endC_l .+ C1_l)
    endC_l = C_l[end]
    # Lines 8 and 9
    if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      # Line 10
      x_l = endC_l + d.κ_l
      R1_l = [maximum(C_l[t_l:end])]
      for i = (t_l-1):-1:(length(R_l)+1)
        push!(R1_l,max(R1_l[end],C_l[i]))
      end
      append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
    else # Line 11
      # Lines 12 and 13
      (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      append!(F_l,F1_l)
      append!(C_l,endC_l .+ C1_l)
      endC_l = C_l[end]
    end # Line 14
  end # Line 15 - left
  while (length(R_r) == σ) || (length(S_r) >= length(C_r))
    t_r = length(C_r)
    # Lines 6 and 7
    (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
    append!(F_r,F1_r)
    append!(C_r,endC_r .+ C1_r)
    endC_r = C_r[end]
    # Lines 8 and 9
    if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      # Line 10
      x_r = endC_r + d.κ_r
      R1_r = [maximum(C_r[t_r:end])]
      for i = (t_r-1):-1:(length(R_r)+1)
        push!(R1_r,max(R1_r[end],C_r[i]))
      end
      append!(R_r,R1_r[end:-1:1] .- C_r[(length(R_r)+1):t_r])
    else # Line 11
      # Lines 12 and 13
      (C1_r,F1_r) = upRW(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      append!(F_r,F1_r)
      append!(C_r,endC_r .+ C1_r)
      endC_r = C_r[end]
    end # Line 14
  end # Line 15 - right
  # Line 16
  while length(U_l) < length(S_l)
    push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
  end
  while length(U_r) < length(S_r)
    push!(U_r,exp(-d.α*F_r[length(U_r)+1]))
  end
  # Lines 17, 18 and 19
  σ = 1
  D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
      exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )
  D_r = exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r)-σ-1)) + sum(S_r[σ:end] .*
      exp.(-d.d_r .* (0:(length(S_r)-σ))) .* (1 .- U_r[σ:end]) .^ ra) )

  push!(lagS_l,S_l[σ])
  push!(lagS_r,S_r[σ])
  push!(lagU_l,U_l[σ])
  push!(lagU_r,U_r[σ])
  # Line 20
  # Construct the dataframes
  (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
  (df_r,err_r) = faces(d.α,D_r,lagS_r,lagU_r)

  # Rescale and add the drift
  df_l = @linq df_l |>
       where(:label .> -1) |>
       select(:length,:height,:slope,:label)
  df_l[:label] = 0
  sort!(df_l, order(:slope, rev = true))
  df_l[:length] = df_l[:length] .* T_l
  df_l[:height] = df_l[:height]  .* (-sc_l) .+ df_l[:length] .* d.μ
  df_l[:slope] = df_l[:slope] .* (-sc1_l) .+ d.μ
  df_r[:length] = df_r[:length] .* T_r
  df_r[:height] = df_r[:height] .* sc_r .+ df_r[:length] .* d.μ
  df_r[:slope] = df_r[:slope] .* sc1_r .+ d.μ
  df_l1 = deepcopy(df_l)
  df_r1 = deepcopy(df_r)

  # Construct (the dataframes with the faces of) the sandwiching convex functions
  dfL_l = @linq df_l1 |>
       select(:length,:height)
  dfU_l = copy(dfL_l)
  dfL_r = @linq df_r1 |>
      where(:label .< 1) |>
      select(:length,:height)
  dfU_r = @linq df_r1 |>
      where(:label .> -1) |>
      select(:length,:height)
  append!(dfL_l,dfL_r)
  append!(dfU_l,dfU_r)
  err_l *= sc_l

  # Check for crossing:
  aux = F(dfU_l)
  G = cross_func(d.f,-err_l)
  if aux == G(dfL_l) # Crossing / non-crossing detection!
    append!(df_l,df_r)
    return (aux,df_l,err_l,err_r * sc_r,σ)
  end # Line 22
  # Line 3 (second to last iterations)
  while true
    # Line 4 (Algorithm 3 with conditional probabilities)
    append!(S_l,chiStable2(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,length(S_l)-σ-1))
    append!(S_r,chiStable2(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,length(S_r)-σ-1))
    # Line 5
    while (length(R_l) == σ) || (length(S_l) >= length(C_l))
      t_l = length(C_l)-1
      # Lines 6 and 7
      (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
      append!(F_l,F1_l)
      append!(C_l,endC_l .+ C1_l)
      endC_l = C_l[end]
      # Lines 8 and 9
      if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        # Line 10
        x_l = endC_l + d.κ_l
        R1_l = [maximum(C_l[t_l:end])]
        for i = (t_l-1):-1:(length(R_l)+1)
          push!(R1_l,max(R1_l[end],C_l[i]))
        end
        append!(R_l,R1_l[end:-1:1] .- C_l[(length(R_l)+1):t_l])
      else # Line 11
        # Lines 12 and 13
        (C1_l,F1_l) = upRW(d.d_l,d.η_l,et1_l,d.κ_l,x_l-endC_l)
        append!(F_l,F1_l)
        append!(C_l,endC_l .+ C1_l)
        endC_l = C_l[end]
      end # Line 14
    end # Line 15
    while (length(R_r) == σ) || (length(S_r) >= length(C_r))
      t_r = length(C_r)-1
      # Lines 6 and 7
      (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
      append!(F_r,F1_r)
      append!(C_r,endC_r .+ C1_r)
      endC_r = C_r[end]
      # Lines 8 and 9
      if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
        # Line 10
        x_r = endC_r + d.κ_r
        R1_r = [maximum(C_r[t_r:end])]
        for i = (t_r-1):-1:(length(R_r)+1)
          push!(R1_r,max(R1_r[end],C_r[i]))
        end
        append!(R_r,R1_r[end:-1:1] .- C_r[(length(R_r)+1):t_r])
      else # Line 11
        # Lines 12 and 13
        (C1_r,F1_r) = upRW(d.d_r,d.η_r,et1_r,d.κ_r,x_r-endC_r)
        append!(F_r,F1_r)
        append!(C_r,endC_r .+ C1_r)
        endC_r = C_r[end]
      end # Line 14
    end # Line 15
    # Line 16
    while length(U_l) < length(S_l)
      push!(U_l,exp(-d.α*F_l[length(U_l)+1]))
    end
    while length(U_r) < length(S_r)
      push!(U_r,exp(-d.α*F_r[length(U_r)+1]))
    end
    # Lines 17, 18 and 19
    σ += 1
    D_l = exp(R_l[σ])*( e2_l*exp(-sep_l*(length(S_l)-σ-1)) + sum(S_l[σ:end] .*
        exp.(-d.d_l .* (0:(length(S_l)-σ))) .* (1 .- U_l[σ:end]) .^ ra) )
    D_r = exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r)-σ-1)) + sum(S_r[σ:end] .*
        exp.(-d.d_r .* (0:(length(S_r)-σ))) .* (1 .- U_r[σ:end]) .^ ra) )
    push!(lagS_l,S_l[σ])
    push!(lagS_r,S_r[σ])
    push!(lagU_l,U_l[σ])
    push!(lagU_r,U_r[σ])
    # Line 24
    # Construct the dataframes
    (df_l,err_l) = faces(d.α,D_l,lagS_l,lagU_l)
    (df_r,err_r) = faces(d.α,D_r,lagS_r,lagU_r)

    # Rescale and add the drift
    df_l = @linq df_l |>
         where(:label .> -1) |>
         select(:length,:height,:slope,:label)
    df_l[:label] = 0
    sort!(df_l, order(:slope, rev = true))
    df_l[:length] = df_l[:length] .* T_l
    df_l[:height] = df_l[:height]  .* (-sc_l) .+ df_l[:length] .* d.μ
    df_l[:slope] = df_l[:slope] .* (-sc1_l) .+ d.μ
    df_r[:length] = df_r[:length] .* T_r
    df_r[:height] = df_r[:height] .* sc_r .+ df_r[:length] .* d.μ
    df_r[:slope] = df_r[:slope] .* sc1_r .+ d.μ
    df_l1 = deepcopy(df_l)
    df_r1 = deepcopy(df_r)

    # Construct (the dataframes with the faces of) the sandwiching convex functions
    dfL_l = @linq df_l1 |>
         select(:length,:height)
    dfU_l = copy(dfL_l)
    dfL_r = @linq df_r1 |>
        where(:label .< 1) |>
        select(:length,:height)
    dfU_r = @linq df_r1 |>
        where(:label .> -1) |>
        select(:length,:height)
    append!(dfL_l,dfL_r)
    append!(dfU_l,dfU_r)
    err_l *= sc_l

    # Check for crossing:
    aux = F(dfU_l)
    G = cross_func(d.f,-err_l)
    if aux == G(dfL_l) # Crossing / non-crossing detection!
      append!(df_l,df_r)
      return (aux,df_l,err_l,err_r * sc_r,σ)
    end # Line 26
  end # Line 27
end

export ConvexMinorantStableMeander, CrossConvexMinorantStableMeander, ConvexMinorantWeaklyStable, CrossConvexMinorantWeaklyStable, rand, params


text:
@testset "Test Second Order Derivative" begin
    @test isapprox(biderivative(x->x^2, 1, 0.00001), 2; atol=1e-5)
    @test isapprox(biderivative(x->1/sin(x), 1, 0.000001), 1/2*(3+cos(2))*csc(1)^3; atol=1e-5)
end

@testset "Test Third Order Derivative" begin
    @test isapprox(triderivative(sin, 1, 0.00001), -cos(1); atol=1e-5)
    @test isapprox(triderivative(x->x^3, 1, 0.00001), 6; atol=1e-5)
    @test isapprox(triderivative(x->x^4, 1, 0.0000001), 24; atol=1e-5)
end

@testset "Test Second Order Derivative" begin
    @test isapprox(quadderivative(sin, 1, 0.000001), sin(1); atol=1e-5)
end

ext: jl


text:
# B.C. Goodwin
# Oscillatory behavior in enzymatic control processes
# https://doi.org/10.1016/0065-2571(65)90067-1

using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
    x2'(t) = alpha * x1(t) - beta * x2(t),
    x3'(t) = gama * x2(t) - delta * x3(t),
    x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
    y(t) = x1(t)
)

@time println(assess_identifiability(ode))


text:
module BSplines

export splineBasisFuns,
       cubicSplineBasisFuns,
       quadraticSplineBasisFuns,
       linearSplineBasisFuns,
       splineBFE,
       cubicSplineBFE,
       quadraticSplineBFE,
       linearSplineBFE

include("code.jl")
##include("Fourier.jl")
##include("Wavelets.jl")

end  ## module


text:
__precompile__()

module BAPOMCP

#=
Current constraints:
- action space discrete
- action space same for all states, histories
- no built-in support for history-dependent rollouts (this could be added though)
- initial n and initial v are 0
=#

using POMDPs
using Parameters
using POMDPToolbox
using ParticleFilters
using CPUTime
using Colors
using Distributions

import POMDPs: action, solve, updater, requirements_info
import POMDPToolbox: action_info

using MCTS
import MCTS: convert_estimator, estimate_value, node_tag, tooltip_tag, default_action

using D3Trees

export
    POMCPSolver,
    POMCPPlanner,

    action,
    solve,
    updater,

    NoDecision,
    AllSamplesTerminal,
    ExceptionRethrow,
    ReportWhenUsed,
    default_action,

    BeliefNode,
    AbstractPOMCPSolver,

    PORollout,
    FORollout,
    RolloutEstimator,
    FOValue,

    D3Tree,
    node_tag,
    tooltip_tag

abstract type AbstractPOMCPSolver <: Solver end

"""
    POMCPSolver(#=keyword arguments=#)

Partially Observable Monte Carlo Planning Solver.

## Keyword Arguments

- `max_depth::Int`
    Rollouts and tree expension will stop when this depth is reached.
    default: `20`

- `c::Float64`
    UCB exploration constant - specifies how much the solver should explore.
    default: `1.0`

- `tree_queries::Int`
    Number of iterations during each action() call.
    default: `1000`

- `max_time::Float64`
    Maximum time for planning in each action() call.
    default: `Inf`

- `tree_in_info::Bool`
    If `true`, returns the tree in the info dict when action_info is called.
    default: `false`

- `estimate_value::Any`
    Function, object, or number used to estimate the value at the leaf nodes.
    default: `RolloutEstimator(RandomSolver(rng))`
    - If this is a function `f`, `f(pomdp, s, h::BeliefNode, steps)` will be called to estimate the value.
    - If this is an object `o`, `estimate_value(o, pomdp, s, h::BeliefNode, steps)` will be called.
    - If this is a number, the value will be set to that number
    Note: In many cases, the simplest way to estimate the value is to do a rollout on the fully observable MDP with a policy that is a function of the state. To do this, use `FORollout(policy)`.

- `default_action::Any`
    Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
    default: `ExceptionRethrow()`
    - If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
    - If this is a Policy `p`, `action(p, belief)` will be called.
    - If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and if this method is not implemented, `a` will be returned directly.

- `rng::AbstractRNG`
    Random number generator.
    default: `Base.GLOBAL_RNG`
"""
@with_kw mutable struct POMCPSolver <: AbstractPOMCPSolver
    max_depth::Int          = 20
    c::Float64              = 1.0
    tree_queries::Int       = 1000
    max_time::Float64       = Inf
    tree_in_info::Bool      = false
    default_action::Any     = ExceptionRethrow()
    rng::AbstractRNG        = Base.GLOBAL_RNG
    estimate_value::Any     = RolloutEstimator(RandomSolver(rng))
end

struct POMCPTree{A,O}
    # for each observation node
    total_n::Vector{Int}            #total number of times visited this node
    children::Vector{Vector{Int}}   #children of node
    o_labels::Vector{O}             #Observation label searched for this node
    o_lookup::Dict{Tuple{Int, O}, Int}  #Look-up action node index + observation --> observation node index

    # for each action node
    n::Vector{Int}                  # number of times visited this node
    v::Vector{Float64}              # Value of action node
    a_labels::Vector{A}             # Action label selected for this node
end

function POMCPTree(pomdp::POMDP, sz::Int=1000)
    acts = collect(iterator(actions(pomdp)))
    A = action_type(pomdp)
    O = obs_type(pomdp)
    sz = min(100_000, sz)
    return POMCPTree{A,O}(sizehint!(Int[0], sz),
                          sizehint!(Vector{Int}[collect(1:length(acts))], sz),
                          sizehint!(Array{O}(1), sz),

                          sizehint!(Dict{Tuple{Int,O},Int}(), sz),

                          sizehint!(zeros(Int, length(acts)), sz),
                          sizehint!(zeros(Float64, length(acts)), sz),
                          sizehint!(acts, sz)
                         )
end

"""DEFINE TREE
"""
#Returns index of observation node
function insert_obs_node!(t::POMCPTree, pomdp::POMDP, ha::Int, o)
    push!(t.total_n, 0)
    push!(t.children, sizehint!(Int[], n_actions(pomdp)))
    push!(t.o_labels, o)
    hao = length(t.total_n)
    t.o_lookup[(ha, o)] = hao
    for a in iterator(actions(pomdp))
        n = insert_action_node!(t, hao, a)
        push!(t.children[hao], n)
    end
    return hao
end

#Returns index of action node
function insert_action_node!(t::POMCPTree, h::Int, a)
    push!(t.n, 0)
    push!(t.v, 0.0)
    push!(t.a_labels, a)
    return length(t.n)
end

abstract type BeliefNode <: AbstractStateNode end

struct POMCPObsNode{A,O} <: BeliefNode
    tree::POMCPTree{A,O}
    node::Int
end

"""DEFINE POLICY (i.e. Planner)
"""
mutable struct POMCPPlanner{P, SE, RNG} <: Policy
    solver::POMCPSolver
    problem::P
    solved_estimator::SE
    rng::RNG
    _best_node_mem::Vector{Int}
    _tree::Nullable
end

function POMCPPlanner(solver::POMCPSolver, pomdp::POMDP)
    #Create a solved estimator to be used for MC Rollouts
    #Uses rollout estimator random choice
    se = convert_estimator(solver.estimate_value, solver, pomdp)
    return POMCPPlanner(solver, pomdp, se, solver.rng, Int[], Nullable())
end

Base.srand(p::POMCPPlanner, seed) = srand(p.rng, seed)

function updater(p::POMCPPlanner)
    P = typeof(p.problem)
    S = state_type(P) #<<-----------------------------!!!!!!!!!!!!!!!!!!
    A = action_type(P)
    O = obs_type(P)
    if !@implemented ParticleFilters.obs_weight(::P, ::S, ::A, ::S, ::O)
        return UnweightedParticleFilter(p.problem, p.solver.tree_queries, rng=p.rng)
    end
    return SIRParticleFilter(p.problem, p.solver.tree_queries, rng=p.rng)
end

# TODO (maybe): implement this for history-dependent policies
#=
immutable AOHistory
    tree::POMCPTree
    tail::Int
end

length
getindex
=#

include("solver.jl")

include("exceptions.jl")
include("rollout.jl")
include("visualization.jl")
include("requirements_info.jl")

end # module


text:
"""
  Runs the energy system model.
"""

module energisystemprojekt

using JuMP, AxisArrays, Gurobi, UnPack, PlotlyJS, CSV, DataFrames

export runmodel

include("input_energisystemprojekt.jl")
include("model_energisystemprojekt.jl")

function runmodel()
    # TODO: kolla plotly!

    input = read_input()

    model = buildmodel(input)

    @unpack m, InstalledCapacity, Electricity, Emission, GeneratedElectricity, TransmissionCapacity = model
    @unpack REGION, PLANT, HOUR, numregions, numplants = input

    println("\nSolving model...")

    status = optimize!(m)

    if termination_status(m) == MOI.OPTIMAL
        println("\nSolve status: Optimal")
    elseif termination_status(m) == MOI.TIME_LIMIT && has_values(m)
        println("\nSolve status: Reached the time-limit")
    else
        error("The model was not solved correctly.")
    end

    # Calculates installed capacity in each region
    RegionalCapacity = AxisArray(zeros(numregions), REGION)
    for r in REGION
       RegionalCapacity[r] = sum(value.(InstalledCapacity[r, p]) for p in PLANT)
    end
    # Installed capacity values
    Capacity = AxisArray(zeros(numregions, numplants), REGION, PLANT)
    for r in REGION
      for p in PLANT
       Capacity[r, p] = value.(InstalledCapacity[r, p])
      end
    end
    # Calculates annual produktion in each region
    AnnualProduction = AxisArray(zeros(numregions), REGION)
    for r in REGION
       AnnualProduction[r] = sum(value.(GeneratedElectricity[r, h]) for h in HOUR)
    end
    # Production in Germany hour 147-651
    GermanyGenerators = AxisArray(zeros(651-147+1), 147:651)
    for h in 147:651
       GermanyGenerators[h-146] = value.(GeneratedElectricity[:DE, h])
    end
    Combination = [:DeSe, :SeDk, :DkDe]
    TransCap = AxisArray(zeros(numregions), Combination) # "Germany-Sweden", "Sweden-Denmark", "Denmark-Germany"
    TransCap[Combination[1]] = value.(TransmissionCapacity[:DE, :SE])
    TransCap[Combination[2]] = value.(TransmissionCapacity[:SE, :DK])
    TransCap[Combination[3]] = value.(TransmissionCapacity[:DK, :DE])

    Cost_result = objective_value(m)/1000000 # €->M€
    Capacity_result = value.(InstalledCapacity)
    Emission_result = value.(Emission) # Mton CO2

    println("Cost (M€): ", Cost_result)
    # println("Capacity (MW): ", Capacity_result)
    println("CO2 Emission (Mton): ", Emission_result)
    println("Annual production (MWh): ", AnnualProduction)
    println("Regional Capacity (MW): ", RegionalCapacity)
    println("Transmission Capacity (MW): ", TransCap)
    # println("Production in Germany (MW): ", GermanyGenerators)

    # Plots:
    # dfInstalledCapacities = DataFrame(A=RegionalCapacity, B=[:Germany, :Sweden, :Denmark])
    # plot_InstalledCapacities() = plot(dfInstalledCapacities, x=dfInstalledCapacities.B, y = dfInstalledCapacities.A, kind = "bar", Layout(title="Installed Capacities (MW)"))

    df = DataFrame(Wind=Capacity[:, :Wind], PV=Capacity[:, :PV], Gas=Capacity[:, :Gas],
      Hydro=Capacity[:, :Hydro], Region=[:Germany, :Sweden, :Denmark], Batteries=Capacity[:, :Batteries], Transmission=Capacity[:, :Transmission], Nuclear=Capacity[:, :Nuclear])
    dfCapacity = stack(df)#, Not([:H]))
    println(dfCapacity)
    plot_Capacity() = plot(dfCapacity, x=:Region, y=:value, color=:variable, kind="bar", Layout(title="Installed Capacities (MW)", barmode="stack"))
    display(plot_Capacity())

    # dfGeneratedElectricity = DataFrame(A=AnnualProduction, B=[:Germany, :Sweden, :Denmark])
    # plot_GeneratedElectricity() = plot(dfGeneratedElectricity, x=dfGeneratedElectricity.B, y = dfGeneratedElectricity.A, kind = "bar", Layout(title="Annual Production (MWh)"))
    #
    # dfGermanyGenerators = DataFrame(A=GermanyGenerators, B=147:651)
    # plot_GermanyGenerators() = plot(dfGermanyGenerators, x=dfGermanyGenerators.B, y = dfGermanyGenerators.A, kind="scatter", mode="lines", Layout(title="Production in Germany hour 147-651 (MW)"))
    #
    # dfTransmissionCapacity = DataFrame(A=TransCap, B=["Germany-Sweden", "Sweden-Denmark", "Denmark-Germany"])
    # plot_TransmissionCapacity() = plot(dfTransmissionCapacity, x=dfTransmissionCapacity.B, y = dfTransmissionCapacity.A, kind = "bar", Layout(title="Transmission Capacity (MW)"))

    # p = [plot_Capacity(); plot_GeneratedElectricity(); plot_GermanyGenerators(); plot_TransmissionCapacity()]
    # relayout!(p, title_text="Exercise 3 plots")
    # display(p)

    # df = dataset(DataFrame, "medals")
    # long_df = stack(df, Not([:nation]), variable_name="medal", value_name="count")
    #
    # plot(long_df, kind="bar", x=:nation, y=:count, color=:medal, Layout(title="Long-Form Input", barmode="stack"))

    nothing

    #Exercise:      1,     2      3      4
    # Total cost: 37238, 66305, 49121, 43619  [M€]
    # Emissions:   139,   13.9,  13.9   13.9 [Mton]

end #runmodel

end # module


text:
abstract type Chromosome end

mutable struct RealChromosome <: Chromosome
    genes::Array{Float64,1}
    cost::Float64
end

function Base.:*(c1::RealChromosome, num::Float64)::RealChromosome
    RealChromosome(c1.genes * num, Inf64)
end

function Base.:*(num::Float64, c1::RealChromosome)::RealChromosome
    RealChromosome(c1.genes * num, Inf64)
end

function Base.:+(c1::RealChromosome, c2::RealChromosome)::RealChromosome
    RealChromosome(c1.genes .+ c2.genes, Inf64)
end

function Base.:-(c1::RealChromosome, c2::RealChromosome)::RealChromosome
    RealChromosome(c1.genes .- c2.genes, Inf64)
end

function Base.isless(c1::RealChromosome, c2::RealChromosome)::Bool
    return c1.cost < c2.cost
end

function LinearCrossover(c1::RealChromosome, c2::RealChromosome)::Tuple{RealChromosome,RealChromosome,RealChromosome}
    offspring1 = 0.5 * c1 + 0.5 * c2
    offspring2 = 1.5 * c1 - 0.5 * c2
    offspring3 = 1.5 * c2 - 0.5 * c1
    return (offspring1, offspring2, offspring3)
end

function ArithmeticCrossOver(c1::RealChromosome, c2::RealChromosome)::Tuple{RealChromosome,RealChromosome}
    alpha = rand()
    return alpha * c1 + (1.0 - alpha) * c2
end

function Mutate(c::RealChromosome, prob::Float64)::RealChromosome
    genes = copy(c.genes)
    @inbounds for i in 1:length(genes)
        if rand(1)[1] < prob
            genes[i] += randn(1)[1]
        end
    end
    return RealChromosome(genes, Inf64)
end

function TournamentSelection(pop::Array{RealChromosome,1})::Tuple{RealChromosome,RealChromosome}
    n = length(pop)
    indices = sample(collect(1:n), 4, replace=false)
    lucky1 = nothing
    lucky2 = nothing
    if pop[indices[1]].cost < pop[indices[2]].cost
        lucky1 = pop[indices[1]]
    else
        lucky1 = pop[indices[2]]
    end
    if pop[indices[3]].cost < pop[indices[4]].cost
        lucky2 = pop[indices[3]]
    else
        lucky2 = pop[indices[4]]
    end
    return (lucky1, lucky2)
end

function createPopulation(popsize::Int,
        chsize::Int,
        mins::Array{Float64,1},
        maxs::Array{Float64,1})::Array{RealChromosome,1}
    pop = Array{RealChromosome,1}(undef, popsize)
    @inbounds for i in 1:popsize
        c = RealChromosome(mins .+ rand(chsize) .* (maxs - mins), Inf64)
        pop[i] = c
    end
    return pop
end

# Change the number of threads before starting Julia
# for using multi-core fitness evaluation.
# Example:
# $ export JULIA_NUM_THREADS=4
function Evaluate(pop::Array{RealChromosome,1}, fcost::Function)::Array{RealChromosome,1}
    Base.Threads.@threads for i in eachindex(pop)
        @inbounds pop[i].cost = fcost(pop[i].genes)
    end
    return pop
end

function Generation(pop::Array{RealChromosome,1},
                    fcost::Function,
                    elitism::Int,
                    pcross::Float64,
                    pmutate::Float64)::Array{RealChromosome,1}
    popsize = length(pop)
    newpop = []
    pop = sort(Evaluate(pop, fcost))
    @inbounds for i in 1:elitism
        push!(newpop, pop[i])
    end
    @inbounds while length(newpop) != popsize
        parent1, parent2 = TournamentSelection(pop)
        winner1 = nothing
        winner2 = nothing
        if rand(1)[1] < pcross
            offspring1, offspring2, offspring3 = LinearCrossover(parent1, parent2)
            offpop = sort(Evaluate([offspring1, offspring2, offspring3], fcost))
            winner1 = Mutate(offpop[1], pmutate)
            winner2 = Mutate(offpop[2], pmutate)
            Evaluate([winner1, winner2], fcost)
        else
            winner1, winner2 = parent1, parent2
        end
        if (length(newpop) < popsize)
            push!(newpop, winner1)
        end
        if (length(newpop) < popsize)
            push!(newpop, winner2)
        end
    end
    return sort(newpop)
end

function ga(
            popsize::Int,
            chsize::Int,
            fcost::Function,
            mins::Array{Float64,1},
            maxs::Array{Float64,1},
            pcross::Float64,
            pmutate::Float64,
            elitisim::Int,
            iterations::Int
            )::Array{RealChromosome,1}
    pop = createPopulation(popsize, chsize, mins, maxs)
    for i in 1:iterations
        pop = Generation(pop, fcost, elitisim, pcross, pmutate)
    end
    return pop
end

max_stars_repo_name: angeris/LinRegOutliers

max_issues_repo_name: angeris/LinRegOutliers

max_forks_repo_name: angeris/LinRegOutliers


text:
# Autogenerated wrapper script for startin_jll for x86_64-apple-darwin
export libstartin

JLLWrappers.@generate_wrapper_header("startin")
JLLWrappers.@declare_library_product(libstartin, "@rpath/libstartin.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libstartin,
        "lib/libstartin.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()