text:
using Distributions

struct StableMeander <: ContinuousUnivariateDistribution
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  StableMeander(α,β) = ( β<-1 || β>1 || (β==-1 && α<=1) || α<=0 || α>2 ) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") :
    β==1 && α<=1 ? Stable(α,1) :
    new(α,β,β*(α<=1 ? 1 : (α-2)/α),(1+β*(α<=1 ? 1 : (α-2)/α))/2)
end

import Distributions.params
function params(d::StableMeander)
  return (d.α,d.β,d.θ,d.ρ)
end

import Distributions.minimum
function minimum(d::StableMeander)
  return 0
end

import Distributions.maximum
function maximum(d::StableMeander)
  return Inf
end

import Distributions.insupport
function insupport(d::StableMeander,x::Real)
  return x>=0
end

import Distributions.mean
function mean(d::StableMeander)
  return pi/(sin(pi/d.α) * gamma(d.ρ + 1/d.α) * gamma(1 - d.ρ))
end

struct PreciseStableMeander <: Sampleable{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  d::Float64
  η::Float64
  δ::Float64
  γ::Float64
  κ::Float64
  Δ::Int64
  mAst::Int64
  ε::Float64
  PreciseStableMeander(α::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, ε)
  PreciseStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer) = PreciseStableMeander(α,β,d,δ,γ,κ,Δ,mAst,.5^32)
  PreciseStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer) = PreciseStableMeander(α,β,d,δ,γ,κ,Δ,12)
  PreciseStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real) = PreciseStableMeander(α,β,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))))))
  PreciseStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real) = PreciseStableMeander(α,β,d,δ,γ,4)
  PreciseStableMeander(α::Real,β::Real) = PreciseStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α)
  PreciseStableMeander(α::Real,β::Real,Δ::Integer,ε::Real) = PreciseStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,4,Δ,12,ε)
  PreciseStableMeander(α::Real,β::Real,ε::Real) = PreciseStableMeander(α,β,
    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

function precise_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer,ε::Real)
  return PreciseStableMeander(d.α,d.β,dr,δ,γ,κ,Δ,mAst,ε)
end

function precise_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer)
  return PreciseStableMeander(d.α,d.β,dr,δ,γ,κ,Δ,mAst)
end

function precise_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,Δ::Integer)
  return PreciseStableMeander(d.α,d.β,dr,δ,γ,κ,Δ)
end

function precise_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real)
  return PreciseStableMeander(d.α,d.β,dr,δ,γ,κ)
end

function precise_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real)
  return PreciseStableMeander(d.α,d.β,dr,δ,γ)
end

function precise_sampler(d::StableMeander,Δ::Integer,ε::Real)
  return PreciseStableMeander(d.α,d.β,Δ,ε)
end

function precise_sampler(d::StableMeander,ε::Real)
  return PreciseStableMeander(d.α,d.β,ε)
end

function precise_sampler(d::StableMeander)
  return PreciseStableMeander(d.α,d.β)
end

import Base.length
function length(d::PreciseStableMeander)
  return 3
end

using LambertW
# @Input: normalised drift d'=αρd∈(0,1)
# @Returns: normalised η'=η/(αρ)>0
function etaF(d::Real)
  return -lambertw(-d*exp(-d),-1)/d-1
end

# @Input: drift d, η, 1/(αρ+η), and a bound M we want to know if we exceed
# @Returns: a sample of 1{R_0 <= M}
function reflectedProcess(d::Real,et::Real,et1::Real,M::Real)
  T = rand(Exponential())/et
  if T <= M
    return true
  end
  Sn = [0.]
  while Sn[end] < T
    push!(Sn,Sn[end]-rand(Exponential())*et1+d)
  end
  return maximum(Sn[1:(end-1)]) <= M
end

# @Input: drift d, η, 1/(αρ+η), a bound M that it should not exceed,
# and another bound M1 we want to know if we exceed
# @Returns: a sample of 1{R_0 <= M1} given R_0<M
function reflectedProcess(d::Real,et::Real,et1::Real,M::Real,M1::Real)
  m = exp(-M-d)
  Sn = [0.]
  while true
    T = -log(rand(Uniform(m,1)))/et # Exponential conditioned on <= M+d, since T>M+d always yields a rejection
    if T <= M1
      return true
    end
    while Sn[end] < T
      push!(Sn,Sn[end]-rand(Exponential())*et1+d)
    end
    max = maximum(Sn[1:(end-1)])
    if max <= M
      return max <= M1
    else
      Sn = [0.]
    end
  end
end

# @Input: drift d, η, 1/(αρ+η), κ and an upper bound x
# @Returns: A simulated run of the random walk up to the hitting time of the barrier
# -2κ conditioned on never exceeding x
function downRW(d::Real,rar::Real,et::Real,et1::Real,κ::Real,x::Real)
  Sn = [0.]
  INCn = [0.]
  if isfinite(x)
    while true
      while Sn[end] > -2 *κ
        push!(INCn,rand(Exponential())*rar)
        push!(Sn,Sn[end]-INCn[end]+d)
      end
      if maximum(Sn) < x && reflectedProcess(d,et,et1,x-Sn[end])
        return (Sn[2:end],INCn[2:end])
      else
        Sn = [0.]
        INCn = [0.]
      end
    end
  else
    while Sn[end] > -2 *κ
      push!(INCn,rand(Exponential())*rar)
      push!(Sn,Sn[end]-INCn[end]+d)
    end
    return (Sn[2:end],INCn[2:end])
  end
end

# @Input: drift d, η, 1/(αρ+η), κ, an upper bound x that will be attained
# @Returns: A simulated run of the random walk up to the hitting time of the barrier
# κ conditioned on never exceeding x
function upRW(d::Real,et::Real,et1::Real,κ::Real,x::Real)
  Sn = [0.]
  INCn = [0.]
  if x == Inf
    while true
      while Sn[end] < κ
        push!(INCn,rand(Exponential())*et1)
        push!(Sn,Sn[end]-INCn[end]+d)
      end
      if rand(Uniform()) <= exp(-et*Sn[end])
        return (Sn[2:end],INCn[2:end])
      else
        Sn = [0.]
        INCn = [0.]
      end
    end
  else
    while true
      while Sn[end] < κ
        push!(INCn,rand(Exponential())*et1)
        push!(Sn,Sn[end]-INCn[end]+d)
      end
      if rand(Uniform()) <= exp(-et*Sn[end]) && reflectedProcess(d,et,et1,x-Sn[end]) # The condition maximum(Sn) < x1 is trivially satisfied by construction of κ
        return (Sn[2:end],INCn[2:end])
      else
        Sn = [0.]
        INCn = [0.]
      end
    end
  end
end

# Blanchet-Sigman algorithm
# @Input: drift d, η, 1/(αρ+η), κ and an upper bound x, last value of random walk s,
# @Returns: A simulated run of the random walk up to the first time we know an upper bound for the walk,
# the increments of the RW, the new upper bound, and the last time at which we may compute D_n
function BSAlgorithm(d::Real,et::Real,et1::Real,κ::Real,x::Real,s::Real)
  Sn = [Float64(s)]
  INCn = [0.]
  x1 = x
  t = 0

  while true
    # Step 1
    (Sn1,INCn1) = downRW(d,et,et1,κ,x1-Sn[end])
    # Step 2
    append!(Sn,Sn[end] .+ Sn1)
    append!(INCn,INCn1)
    t = t+length(Sn1)
    # Step 3 and 4
    if reflectedProcess(d,et,et1,x1-Sn[end],κ)
      x1 = Sn[end]+κ
      t = t-length(Sn1)
      break
    else
      (Sn1,INCn1) = upRW(d,et,et1,κ,x1-Sn[end])
      append!(Sn,Sn[end] .+ Sn1)
      append!(INCn,INCn1)
      t = t+length(Sn1)
    end
  end

  return (Sn[2:end],INCn[2:end],x1,t)
end

# Algorithm 3 - 1 (unconditional probabilities)
# @Input: S^+(α,ρ), δ, dega = δ*γ, dega2 = -1/(1-exp(-dega)), mom = E(S_1^γ), shift = m-k >= m*-1
function chiStable1(d::PositiveStable,δ::Real,dega::Real,dega2::Real,mom::Real,shift::Integer)
  U = rand(Uniform())
  n = shift
  S = Float64[]
  while true
    n = n+1
    aux = exp(δ*n)
    p = cdf(d,aux) # 1-pn
    aux1 = exp(-(n+1)*dega)*mom
    q = p*exp(dega2*aux1/(1-aux1)) # qn
    if U > p
      push!(S,cond_rand(d,aux,Inf))
      U = (U-p)/(1-p)
    elseif U < q
      push!(S,cond_rand(d,0,aux))
      return S
    else
      push!(S,cond_rand(d,0,aux))
      U = U/p
    end
  end
end

# Algorithm 3 - 2 (conditional probabilities)
# @Input: S^+(α,ρ), δ, dega = δ*γ, dega2 = -1/(1-exp(-dega)), mom = E(S_1^γ)
function chiStable2(d::PositiveStable,δ::Real,dega::Real,dega2::Real,mom::Real,shift::Integer)
  U = rand(Uniform())
  n = shift
  S = Float64[]
  while true
    n = n+1
    aux = exp(δ*n)
    aux0 = exp(δ*(n+1))
    p = cdf(d,aux)/cdf(d,aux0) # 1-pn
    aux1 = exp(-(n+1)*dega)*mom
    q = p*exp(dega2*aux1/(1 -aux1)) # qn
    if U > p
      push!(S,cond_rand(d,aux,aux0))
      U = (U-p)/(1 -p)
    elseif U < q
      push!(S,cond_rand(d,0,aux))
      return S
    else
      push!(S,cond_rand(d,0,aux))
      U = U/p
    end
  end
end

import Distributions.rand

# Algorithm 9
# @Input: PreciseStableMeander
# @Output: A random sample from a the law d, σ
function rand(d::PreciseStableMeander)
  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
  # Line 20
  if D <= eps # Finished!
    # Line 21
    D = (prod(lagU)*U[1])^ra * D
    X = (1-U[1])^ra * S[1]
    for i = d.Δ:-1:1
      X = lagU[i]^ra * X + (1-lagU[i])^ra * lagS[i]
    end
    return (X,X+D,σ)
  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
    # Line 24
    if D <= eps # Finished!
      D = (prod(U[1:σ])*prod(lagU))^ra * D
      X = 0.
      for i = σ:-1:1
        X = U[i]^ra * X + (1-U[i])^ra * S[i]
      end
      for i = d.Δ:-1:1
        X = lagU[i]^ra * X + (1-lagU[i])^ra * lagS[i]
      end
      # Line 25
      return (X,X+D,σ)
    end # Line 26
  end # Line 27
end

# @Input: StableMeander
# @Output: A random sample (lower estimate) from the law d
function rand(d::StableMeander)
  return rand(precise_sampler(d))[1]
end

function rand(d::StableMeander,N::Integer)
  return [rand(d) for i =1:N]
end

function length(d::PreciseStableMeander)
  return 3
end

# Distribution created to sample ε-strongly and localise the function according
# to the partition induced by a discrete function f.
struct LocStableMeander <: Distribution{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  d::Float64
  η::Float64
  δ::Float64
  γ::Float64
  κ::Float64
  Δ::Int64
  mAst::Int64
  f::Function
  LocStableMeander(α::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, f)
  LocStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,f::Function) = LocStableMeander(α,β,d,δ,γ,κ,Δ,12,f)
  LocStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,κ::Real,f::Function) = LocStableMeander(α,β,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)
  LocStableMeander(α::Real,β::Real,d::Real,δ::Real,γ::Real,f::Function) = LocStableMeander(α,β,d,δ,γ,4,f)
  LocStableMeander(α::Real,β::Real,Δ::Integer,f::Function) = LocStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,4,Δ,f)
  LocStableMeander(α::Real,β::Real,f::Function) = LocStableMeander(α,β,(2/3)*2/(α+β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,f)
end

function local_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,mAst::Integer,f::Function)
  return LocStableMeander(d.α,d.β,dr,δ,γ,κ,Δ,mAst,f)
end

function local_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,Δ::Integer,f::Function)
  return LocStableMeander(d.α,d.β,dr,δ,γ,κ,Δ,f)
end

function local_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,κ::Real,f::Function)
  return LocStableMeander(d.α,d.β,dr,δ,γ,κ,f)
end

function local_sampler(d::StableMeander,dr::Real,δ::Real,γ::Real,f::Function)
  return LocStableMeander(d.α,d.β,dr,δ,γ,f)
end

function local_sampler(d::StableMeander,Δ::Integer,f::Function)
  return LocStableMeander(d.α,d.β,Δ,f)
end

function local_sampler(d::StableMeander,f::Function)
  return LocStableMeander(d.α,d.β,f)
end

# Algorithm 9
# @Input: LocStableMeander
# @Output: A random sample from (f(StableMeander),UE,LE),
# where UE and LE are the upper and lower estimates of StableMeander.
function rand(d::LocStableMeander)
  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 = 1 / 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
  # Line 20
  # Line 21
  X = (1-U[1])^ra * S[1]
  for i = d.Δ:-1:1
    X = lagU[i]^ra * X + (1-lagU[i])^ra * lagS[i]
  end
  aux = d.f(X)
  if aux == d.f(X+D/eps)
    return (aux,X,X+D/eps,σ)
  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) )
    X += (1-U[σ])^ra * S[σ] / eps
    eps /= U[σ]^ra
    # Line 24
    aux = d.f(X)
    if aux == d.f(X+D/eps)
      return (aux,X,X+D/eps,σ)
    end # Line 26
  end # Line 27
end

function length(d::LocStableMeander)
  return 4
end

# Finite-dimensional stable meander

struct MvStableMeander <: ContinuousMultivariateDistribution
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  t::Array{Float64,1}
  MvStableMeander(α,β,t) = ( β<-1 || β>1 || (β==-1 && α<=1) || α<=0 || α>2 ) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") :
    !( 0 < t[1] && prod([ t[i-1] < t[i] for i = 2:length(t)]) ) ? error("Parameters' requirements unmet: 0 < t[1] < ... < t[end]") :
    new(α,β,β*(α<=1 ? 1 : (α-2)/α),(1+β*(α<=1 ? 1 : (α-2)/α))/2,t)
end

import Distributions.params
function params(d::MvStableMeander)
  return (d.α,d.β,d.θ,d.ρ,d.t)
end

import Distributions.insupport
function insupport(d::MvStableMeander,x::Array{Float64,1})
  return prod(x .>= 0)
end

import Base.length
function length(d::MvStableMeander)
  return length(d.t)
end

import Distributions.rand
function rand(d::MvStableMeander)
  if d.α <= 1 && d.β == 1
    X = rand(Stable(d.α,1),length(d.t))
    X[1] *= d.t[1]^(1/d.α)
    for i=2:length(d.t)
      X[i] *= (d.t[i]-d.t[i-1])^(1/d.α)
    end
    return cumsum(X)
  else
    return rand(precise_sampler(d))[1]
  end
end

# Does not really support α≤1 and β=1 (produces an exact MvStableMeander)
struct PreciseMvStableMeander <: Sampleable{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  t::Array{Float64,1}
  # 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
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},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,ε::Real) =
    ( β<-1 || β>1 || (β==-1 && α<=1) || α<=0 || α>2 ) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") :
    (β == 1 && α <= 1) ? MvStableMeander(α,1,t) :
    (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") : !( 0 < t[1] && prod([ t[i-1] < t[i] for i = 2:length(t)]) ) ?
    error("Parameters' requirements unmet: 0 < t[1] < ... < t[end]") :
    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, ε)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},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) = PreciseMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r,.5^16)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer) = PreciseMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,12,12)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real) = PreciseMvStableMeander(α,β,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))))))
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real) = PreciseMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,4,4)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1}) = PreciseMvStableMeander(α,β,t,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},Δ_l::Integer,Δ_r::Integer,ε::Real) = PreciseMvStableMeander(α,β,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,ε)
  PreciseMvStableMeander(α::Real,β::Real,t::Array{Real,1},ε::Real) = PreciseMvStableMeander(α,β,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

function precise_sampler(d::MvStableMeander,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,ε::Real)
  return PreciseMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r,ε)
end

function precise_sampler(d::MvStableMeander,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)
  return PreciseMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r)
end

function precise_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer)
  return PreciseMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r)
end

function precise_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real)
  return PreciseMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r)
end

function precise_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real)
  return PreciseMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r)
end

function precise_sampler(d::MvStableMeander,Δ_l::Integer,Δ_r::Integer,ε::Real)
  return PreciseMvStableMeander(d.α,d.β,d.t,Δ_l,Δ_r,ε)
end

function precise_sampler(d::MvStableMeander,ε::Real)
  return PreciseMvStableMeander(d.α,d.β,d.t,ε)
end

function precise_sampler(d::MvStableMeander)
  return PreciseMvStableMeander(d.α,d.β,d.t)
end

import Base.length
function length(d::PreciseMvStableMeander)
  return Int(2*length(d.t) + 1)
end

function rand(d::PreciseMvStableMeander)
  ra = 1 /d.α

  if length(d.t) == 1
    return rand(PreciseStableMeander(d.α,d.β,d.d_r,δ_r,γ_r,κ_r,Δ_r,mAst_r,ε)) * t[1]^ra
  end

  m_r = length(d.t)
  m_l = Int(m_r-1)

  while true
    τ = rand(Beta(1-d.ρ,d.ρ),m_l)
    T_l = [d.t[i]-d.t[i-1] for i = 2:length(d.t)] .* τ
    T_r = [dt[1]]
    append!(T_r, [d.t[i]-d.t[i-1] for i = 2:length(d.t)] .* (1 .- τ))

    T_l = T_l .^ ra
    T_r = T_r .^ ra

    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[] for i = 1:m_l], [Float64[] for i = 1:m_r]
    # Line 1
    x_l, x_r = [Inf for i = 1:m_l], [Inf for i = 1:m_r]
    σ = 0
    # Line 2
    lagU_l, lagU_r = [rand(Beta(1-d.ρ,1),d.Δ_l) for i = 1:m_l], [rand(Beta(d.ρ,1),d.Δ_r) for i = 1:m_r]
    lagS_l, lagS_r = [rand(dPos_l,d.Δ_l) for i = 1:m_l], [rand(dPos_r,d.Δ_r) for i = 1:m_r]
    prod_l, prod_r = [prod(lagU_l[i]) for i = 1:m_l], [prod(lagU_r[i]) for i = 1:m_r]
    # Line 3 (first iteration)
    # Line 4 (Algorithm 3 with unconditional probabilities)
    S_l, S_r = [rand(dPos_l,mAst_l) for i = 1:m_l], [rand(dPos_r,mAst_r) for i = 1:m_r]
    for i = 1:m_l
      append!(S_l[i],chiStable1(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,mAst_l-1))
    end
    for i = 1:m_r
      append!(S_r[i],chiStable1(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,mAst_r-1))
    end

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

    for i = 1:m_l
      push!(lagS_l[i],S_l[i][σ])
      push!(lagU_l[i],U_l[i][σ])
    end
    for i = 1:m_r
      push!(lagS_r[i],S_r[i][σ])
      push!(lagU_r[i],U_r[i][σ])
    end
    # Line 20
    # Line 21
    L_l = [(prod_l[i]*U_l[i][1])^ra for i = 1:m_l]
    L_r = [(prod_r[i]*U_r[i][1])^ra for i = 1:m_r]

    X_l = [(1-U_l[i][1])^ra * S_l[i][1] for i = 1:m_l]
    X_r = [(1-U_r[i][1])^ra * S_r[i][1] for i = 1:m_r]

    for i = 1:m_l
      for j = d.Δ_l:-1:1
        X_l[i] = lagU_l[i][j]^ra * X_l[i] + (1-lagU_l[i][j])^ra * lagS_l[i][j]
      end
    end
    for i = 1:m_r
      for j = d.Δ_r:-1:1
        X_r[i] = lagU_r[i][j]^ra * X_r[i] + (1-lagU_r[i][j])^ra * lagS_r[i][j]
      end
    end

    CX_lu = [0.]
    append!(CX_lu, cumsum(X_l .* T_l))
    CX_ld = [0.]
    append!(CX_ld, cumsum((X_l .+ D_l .* L_l) .* T_l))
    CX_ru = cumsum((X_r .+ D_r .* L_r) .* T_r)
    CX_rd = cumsum(X_r .* T_r)
    X_u = CX_ru .- CX_lu
    X_d = CX_rd .- CX_ld

    if !(prod(X_u .> 0))
      continue
    elseif X_u[end] - X_d[end] < d.ε && prod(X_d .> 0)# Precision reached!
      return (X_d,X_u,σ)
    end # Line 22
    # Line 3 (second to last iterations)
    while true
      # Line 4 (Algorithm 3 with conditional probabilities)
      for i = 1:m_l
        append!(S_l[i],chiStable2(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,length(S_l[i])-σ-1))
      end
      for i = 1:m_r
        append!(S_r[i],chiStable2(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,length(S_r[i])-σ-1))
      end
      # Line 5
      for i = 1:m_l
        while (length(R_l[i]) == σ) || (length(S_l[i]) >= length(C_l[i]))
          t_l = length(C_l[i])-1
          # Lines 6 and 7
          (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l[i]-endC_l[i])
          append!(F_l[i],F1_l)
          append!(C_l[i],endC_l[i] .+ C1_l)
          endC_l[i] = C_l[i][end]
          # Lines 8 and 9
          if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l[i]-endC_l[i])
            # Line 10
            x_l[i] = endC_l[i] + d.κ_l
            R1_l = [maximum(C_l[i][t_l:end])]
            for j = (t_l-1):-1:(length(R_l[i])+1)
              push!(R1_l,max(R1_l[end],C_l[i][j]))
            end
            append!(R_l[i],R1_l[end:-1:1] .- C_l[i][(length(R_l[i])+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[i]-endC_l[i])
            append!(F_l[i],F1_l)
            append!(C_l[i],endC_l .+ C1_l)
            endC_l[i] = C_l[i][end]
          end # Line 14
        end # Line 15
      end
      for i = 1:m_r
        while (length(R_r[i]) == σ) || (length(S_r[i]) >= length(C_r[i]))
          t_r = length(C_r[i])-1
          # Lines 6 and 7
          (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r[i]-endC_r[i])
          append!(F_r[i],F1_r)
          append!(C_r[i],endC_r[i] .+ C1_r)
          endC_r[i] = C_r[i][end]
          # Lines 8 and 9
          if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r[i]-endC_r[i])
            # Line 10
            x_r[i] = endC_r[i] + d.κ_r
            R1_r = [maximum(C_r[t_r:end])]
            for j = (t_r-1):-1:(length(R_r[i])+1)
              push!(R1_r,max(R1_r[end],C_r[i][j]))
            end
            append!(R_r[i],R1_r[end:-1:1] .- C_r[i][(length(R_r[i])+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[i]-endC_r[i])
            append!(F_r[i],F1_r)
            append!(C_r[i],endC_r[i] .+ C1_r)
            endC_r[i] = C_r[i][end]
          end # Line 14
        end # Line 15
      end
      # Line 16
      for i = 1:m_l
        while length(U_l[i]) < length(S_l[i])
          push!(U_l[i],exp(-d.α*F_l[i][length(U_l[i])+1]))
        end
      end
      for i = 1:m_r
        while length(U_r[i]) < length(S_r[i])
          push!(U_r[i],exp(-d.α*F_r[i][length(U_r[i])+1]))
        end
      end
      # Lines 17, 18 and 19
      σ += 1
      D_l = [exp(R_l[i][σ])*( e2_l*exp(-sep_l*(length(S_l[i])-σ-1)) + sum(S_l[i][σ:end] .*
          exp.(-d.d_l .* (0:(length(S_l[i])-σ))) .* (1 .- U_l[i][σ:end]) .^ ra) ) for i = 1:m_l]
      D_r = [exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r[i])-σ-1)) + sum(S_r[i][σ:end] .*
          exp.(-d.d_r .* (0:(length(S_r[i])-σ))) .* (1 .- U_r[i][σ:end]) .^ ra) ) for i = 1:m_r]
      for i = 1:m_l
        push!(lagS_l[i],S_l[i][σ])
        push!(lagU_l[i],U_l[i][σ])
      end
      for i = 1:m_r
        push!(lagS_r[i],S_r[i][σ])
        push!(lagU_r[i],U_r[i][σ])
      end
      # Line 24
      X_l .+= [L_l[i]*(1-U_l[i][σ])^ra * S_l[i][σ] for i = 1:m_l]
      X_r .+= [L_r[i]*(1-U_r[i][σ])^ra * S_r[i][σ] for i = 1:m_r]

      L_l .*= [U_l[i][σ]^ra for i = 1:m_l]
      L_r .*= [U_r[i][σ]^ra for i = 1:m_r]

      CX_lu = [0.]
      append!(CX_lu, cumsum(X_l .* T_l))
      CX_ld = [0.]
      append!(CX_ld, cumsum((X_l .+ D_l .* L_l) .* T_l))
      CX_ru = cumsum((X_r .+ D_r .* L_r) .* T_r)
      CX_rd = cumsum(X_r .* T_r)
      X_u = CX_ru .- CX_lu
      X_d = CX_rd .- CX_ld

      if !(prod(X_u .> 0))
        continue
      elseif X_u[end] - X_d[end] < d.ε && prod(X_d .> 0)# Precision reached!
        return (X_d,X_u,σ)
      end # Line 26
    end # Line 27
  end
end

struct LocMvStableMeander <: Sampleable{Multivariate,Continuous}
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  t::Array{Float64,1}
  # 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
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},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 || (β==-1 && α<=1) || α<=0 || α>2 ) ?
    error("Parameters' requirements unmet: (α,β)∈(0,2]×[-1,1]-(0,1]×{-1}") :
    #(β == 1 && α <= 1) ? MvStableMeander(α,1,t) :
    (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 < t[1] && prod([ t[i-1] < t[i] for i = 2:length(t)]) ) ?
    error("Parameters' requirements unmet: 0 < t[1] < ... < t[end]") :
    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)
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},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) = LocMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r,.5^16,f)
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,f::Function) = LocMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,12,12,f)
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,f::Function) = LocMvStableMeander(α,β,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))))))
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,f::Function) = LocMvStableMeander(α,β,t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,4,4,f)
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},f::Function) = LocMvStableMeander(α,β,t,
    (2/3)*2/(α-β*(α <= 1 ? α : α-2)),(2/3)*2/(α+β*(α <= 1 ? α : α-2)),
    (1/3)*2/(α-β*(α <= 1 ? α : α-2)),(1/3)*2/(α+β*(α <= 1 ? α : α-2)),.95*α,.95*α,f)
  LocMvStableMeander(α::Real,β::Real,t::Array{Real,1},Δ_l::Integer,Δ_r::Integer,f::Function) = LocMvStableMeander(α,β,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)
end

function local_sampler(d::MvStableMeander,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)
  return LocMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,mAst_l,mAst_r,f)
end

function local_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,Δ_l::Integer,Δ_r::Integer,f::Function)
  return LocMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,Δ_l,Δ_r,f)
end

function local_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,κ_l::Real,κ_r::Real,f::Function)
  return LocMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,κ_l,κ_r,f)
end

function local_sampler(d::MvStableMeander,d_l::Real,d_r::Real,δ_l::Real,δ_r::Real,γ_l::Real,γ_r::Real,f::Function)
  return LocMvStableMeander(d.α,d.β,d.t,d_l,d_r,δ_l,δ_r,γ_l,γ_r,f)
end

function local_sampler(d::MvStableMeander,Δ_l::Integer,Δ_r::Integer,f::Function)
  return LocMvStableMeander(d.α,d.β,d.t,Δ_l,Δ_r,f)
end

function local_sampler(d::MvStableMeander,f::Function)
  return LocMvStableMeander(d.α,d.β,d.t,f)
end

import Base.length
function length(d::LocMvStableMeander)
  return Int(2*length(d.t) + 2)
end

function rand(d::LocMvStableMeander)
  ra = 1 /d.α

  if d.β == 1 && d.α <= 1
    return f(rand(MvStableMeander(d.α,1,d.t)))
  end

  if length(d.t) == 1
    g(x) = d.f(x* t[1]^ra)
    return rand(LocStableMeander(d.α,d.β,d.d_r,δ_r,γ_r,κ_r,Δ_r,mAst_r,g))
  end

  m_r = length(d.t)
  m_l = Int(m_r-1)

  while true
    τ = rand(Beta(1-d.ρ,d.ρ),m_l)
    T_l = [d.t[i]-d.t[i-1] for i = 2:length(d.t)] .* τ
    T_r = [dt[1]]
    append!(T_r, [d.t[i]-d.t[i-1] for i = 2:length(d.t)] .* (1 .- τ))

    T_l = T_l .^ ra
    T_r = T_r .^ ra

    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[] for i = 1:m_l], [Float64[] for i = 1:m_r]
    # Line 1
    x_l, x_r = [Inf for i = 1:m_l], [Inf for i = 1:m_r]
    σ = 0
    # Line 2
    lagU_l, lagU_r = [rand(Beta(1-d.ρ,1),d.Δ_l) for i = 1:m_l], [rand(Beta(d.ρ,1),d.Δ_r) for i = 1:m_r]
    lagS_l, lagS_r = [rand(dPos_l,d.Δ_l) for i = 1:m_l], [rand(dPos_r,d.Δ_r) for i = 1:m_r]
    prod_l, prod_r = [prod(lagU_l[i]) for i = 1:m_l], [prod(lagU_r[i]) for i = 1:m_r]
    # Line 3 (first iteration)
    # Line 4 (Algorithm 3 with unconditional probabilities)
    S_l, S_r = [rand(dPos_l,mAst_l) for i = 1:m_l], [rand(dPos_r,mAst_r) for i = 1:m_r]
    for i = 1:m_l
      append!(S_l[i],chiStable1(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,mAst_l-1))
    end
    for i = 1:m_r
      append!(S_r[i],chiStable1(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,mAst_r-1))
    end

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

    for i = 1:m_l
      push!(lagS_l[i],S_l[i][σ])
      push!(lagU_l[i],U_l[i][σ])
    end
    for i = 1:m_r
      push!(lagS_r[i],S_r[i][σ])
      push!(lagU_r[i],U_r[i][σ])
    end
    # Line 20
    # Line 21
    L_l = [(prod_l[i]*U_l[i][1])^ra for i = 1:m_l]
    L_r = [(prod_r[i]*U_r[i][1])^ra for i = 1:m_r]

    X_l = [(1-U_l[i][1])^ra * S_l[i][1] for i = 1:m_l]
    X_r = [(1-U_r[i][1])^ra * S_r[i][1] for i = 1:m_r]

    for i = 1:m_l
      for j = d.Δ_l:-1:1
        X_l[i] = lagU_l[i][j]^ra * X_l[i] + (1-lagU_l[i][j])^ra * lagS_l[i][j]
      end
    end
    for i = 1:m_r
      for j = d.Δ_r:-1:1
        X_r[i] = lagU_r[i][j]^ra * X_r[i] + (1-lagU_r[i][j])^ra * lagS_r[i][j]
      end
    end

    CX_lu = [0.]
    append!(CX_lu, cumsum(X_l .* T_l))
    CX_ld = [0.]
    append!(CX_ld, cumsum((X_l .+ D_l .* L_l) .* T_l))
    CX_ru = cumsum((X_r .+ D_r .* L_r) .* T_r)
    CX_rd = cumsum(X_r .* T_r)
    X_u = CX_ru .- CX_lu
    X_d = CX_rd .- CX_ld

    if !(prod(X_u .> 0))
      continue
    end
    aux = d.f(X_u)
    if aux == d.f(X_d) && prod(X_d .> 0)# Precision reached!
      return (aux,X_d,X_u,σ)
    end # Line 22
    # Line 3 (second to last iterations)
    while true
      # Line 4 (Algorithm 3 with conditional probabilities)
      for i = 1:m_l
        append!(S_l[i],chiStable2(dPos_l,d.δ_l,dega_l,dega2_l,mom_l,length(S_l[i])-σ-1))
      end
      for i = 1:m_r
        append!(S_r[i],chiStable2(dPos_r,d.δ_r,dega_r,dega2_r,mom_r,length(S_r[i])-σ-1))
      end
      # Line 5
      for i = 1:m_l
        while (length(R_l[i]) == σ) || (length(S_l[i]) >= length(C_l[i]))
          t_l = length(C_l[i])-1
          # Lines 6 and 7
          (C1_l,F1_l) = downRW(d.d_l,rar_l,d.η_l,et1_l,d.κ_l,x_l[i]-endC_l[i])
          append!(F_l[i],F1_l)
          append!(C_l[i],endC_l[i] .+ C1_l)
          endC_l[i] = C_l[i][end]
          # Lines 8 and 9
          if reflectedProcess(d.d_l,d.η_l,et1_l,d.κ_l,x_l[i]-endC_l[i])
            # Line 10
            x_l[i] = endC_l[i] + d.κ_l
            R1_l = [maximum(C_l[i][t_l:end])]
            for j = (t_l-1):-1:(length(R_l[i])+1)
              push!(R1_l,max(R1_l[end],C_l[i][j]))
            end
            append!(R_l[i],R1_l[end:-1:1] .- C_l[i][(length(R_l[i])+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[i]-endC_l[i])
            append!(F_l[i],F1_l)
            append!(C_l[i],endC_l .+ C1_l)
            endC_l[i] = C_l[i][end]
          end # Line 14
        end # Line 15
      end
      for i = 1:m_r
        while (length(R_r[i]) == σ) || (length(S_r[i]) >= length(C_r[i]))
          t_r = length(C_r[i])-1
          # Lines 6 and 7
          (C1_r,F1_r) = downRW(d.d_r,rar_r,d.η_r,et1_r,d.κ_r,x_r[i]-endC_r[i])
          append!(F_r[i],F1_r)
          append!(C_r[i],endC_r[i] .+ C1_r)
          endC_r[i] = C_r[i][end]
          # Lines 8 and 9
          if reflectedProcess(d.d_r,d.η_r,et1_r,d.κ_r,x_r[i]-endC_r[i])
            # Line 10
            x_r[i] = endC_r[i] + d.κ_r
            R1_r = [maximum(C_r[t_r:end])]
            for j = (t_r-1):-1:(length(R_r[i])+1)
              push!(R1_r,max(R1_r[end],C_r[i][j]))
            end
            append!(R_r[i],R1_r[end:-1:1] .- C_r[i][(length(R_r[i])+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[i]-endC_r[i])
            append!(F_r[i],F1_r)
            append!(C_r[i],endC_r[i] .+ C1_r)
            endC_r[i] = C_r[i][end]
          end # Line 14
        end # Line 15
      end
      # Line 16
      for i = 1:m_l
        while length(U_l[i]) < length(S_l[i])
          push!(U_l[i],exp(-d.α*F_l[i][length(U_l[i])+1]))
        end
      end
      for i = 1:m_r
        while length(U_r[i]) < length(S_r[i])
          push!(U_r[i],exp(-d.α*F_r[i][length(U_r[i])+1]))
        end
      end
      # Lines 17, 18 and 19
      σ += 1
      D_l = [exp(R_l[i][σ])*( e2_l*exp(-sep_l*(length(S_l[i])-σ-1)) + sum(S_l[i][σ:end] .*
          exp.(-d.d_l .* (0:(length(S_l[i])-σ))) .* (1 .- U_l[i][σ:end]) .^ ra) ) for i = 1:m_l]
      D_r = [exp(R_r[σ])*( e2_r*exp(-sep_r*(length(S_r[i])-σ-1)) + sum(S_r[i][σ:end] .*
          exp.(-d.d_r .* (0:(length(S_r[i])-σ))) .* (1 .- U_r[i][σ:end]) .^ ra) ) for i = 1:m_r]
      for i = 1:m_l
        push!(lagS_l[i],S_l[i][σ])
        push!(lagU_l[i],U_l[i][σ])
      end
      for i = 1:m_r
        push!(lagS_r[i],S_r[i][σ])
        push!(lagU_r[i],U_r[i][σ])
      end
      # Line 24
      X_l .+= [L_l[i]*(1-U_l[i][σ])^ra * S_l[i][σ] for i = 1:m_l]
      X_r .+= [L_r[i]*(1-U_r[i][σ])^ra * S_r[i][σ] for i = 1:m_r]

      L_l .*= [U_l[i][σ]^ra for i = 1:m_l]
      L_r .*= [U_r[i][σ]^ra for i = 1:m_r]

      CX_lu = [0.]
      append!(CX_lu, cumsum(X_l .* T_l))
      CX_ld = [0.]
      append!(CX_ld, cumsum((X_l .+ D_l .* L_l) .* T_l))
      CX_ru = cumsum((X_r .+ D_r .* L_r) .* T_r)
      CX_rd = cumsum(X_r .* T_r)
      X_u = CX_ru .- CX_lu
      X_d = CX_rd .- CX_ld

      if !(prod(X_u .> 0))
        continue
      end
      aux = d.f(X_u)
      if aux == d.f(X_d) && prod(X_d .> 0)# Precision reached!
        return (aux,X_d,X_u,σ)
      end # Line 26
    end # Line 27
  end
end

export StableMeander, LocStableMeander, PreciseStableMeander, MvStableMeander, LocMvStableMeander, PreciseMvStableMeander, rand, precise_sampler, local_sampler, minimum, maximum, insupport, mean, params, length


text:
@testset "static assignment to/from array" begin
    inner = StaticAssignment((a=1., b=2.),NamedTuple())
    outer = StaticAssignment((c=3.,), (d=inner, e=inner))

    arr = to_array(outer, Float64)
    @test to_array(outer, Float64) == Float64[3.0, 1.0, 2.0, 1.0, 2.0]

    trie = from_array(outer, Float64[1, 2, 3, 4, 5])
    @test trie[:c] == 1.0
    @test trie[:d => :a] == 2.0
    @test trie[:d => :b] == 3.0
    @test trie[:e => :a] == 4.0
    @test trie[:e => :b] == 5.0
    @test length(get_internal_nodes(trie)) == 2
    @test length(get_leaf_nodes(trie)) == 1
    inner1 = get_internal_node(trie, :d)
    @test length(get_leaf_nodes(inner1)) == 2
    @test length(get_internal_nodes(inner1)) == 0
    inner2 = get_internal_node(trie, :e)
    @test length(get_leaf_nodes(inner2)) == 2
    @test length(get_internal_nodes(inner2)) == 0
end

@testset "dynamic assignment to/from array" begin
    outer = DynamicAssignment()
    set_leaf_node!(outer, :c, 3.)
    inner = DynamicAssignment()
    set_leaf_node!(inner, :a, 1.)
    set_leaf_node!(inner, :b, 2.)
    set_internal_node!(outer, :d, inner)
    set_internal_node!(outer, :e, inner)

    arr = to_array(outer, Float64)
    @test to_array(outer, Float64) == Float64[3.0, 1.0, 2.0, 1.0, 2.0]

    trie = from_array(outer, Float64[1, 2, 3, 4, 5])
    @test trie[:c] == 1.0
    @test trie[:d => :a] == 2.0
    @test trie[:d => :b] == 3.0
    @test trie[:e => :a] == 4.0
    @test trie[:e => :b] == 5.0
    @test length(get_internal_nodes(trie)) == 2
    @test length(get_leaf_nodes(trie)) == 1
    inner1 = get_internal_node(trie, :d)
    @test length(get_leaf_nodes(inner1)) == 2
    @test length(get_internal_nodes(inner1)) == 0
    inner2 = get_internal_node(trie, :e)
    @test length(get_leaf_nodes(inner2)) == 2
    @test length(get_internal_nodes(inner2)) == 0
end

@testset "internal node vector assignment to/from array" begin
    inner = DynamicAssignment()
    set_leaf_node!(inner, :a, 1.)
    set_leaf_node!(inner, :b, 2.)
    outer = vectorize_internal([inner, inner, inner])

    arr = to_array(outer, Float64)
    @test to_array(outer, Float64) == Float64[1, 2, 1, 2, 1, 2]

    trie = from_array(outer, Float64[1, 2, 3, 4, 5, 6])
    @test trie[1 => :a] == 1.0
    @test trie[1 => :b] == 2.0
    @test trie[2 => :a] == 3.0
    @test trie[2 => :b] == 4.0
    @test trie[3 => :a] == 5.0
    @test trie[3 => :b] == 6.0
    @test length(get_internal_nodes(trie)) == 3
end

@testset "dynamic assignment merge" begin
    inner = DynamicAssignment()
    set_leaf_node!(inner, :x, 1)
    trie1 = DynamicAssignment()
    set_leaf_node!(trie1, :a, 1.)
    set_leaf_node!(trie1, :b, 2.)
    set_internal_node!(trie1, :c, inner)
    set_internal_node!(trie1, :shared, inner)
    trie2 = DynamicAssignment()
    set_leaf_node!(trie2, :d, 3.)
    set_internal_node!(trie2, :e, inner)
    set_internal_node!(trie2, :f, inner)
    inner2 = DynamicAssignment()
    set_leaf_node!(inner2, :y, 4.)
    set_internal_node!(trie2, :shared, inner2)
    trie = merge(trie1, trie2)
    @test trie[:a] == 1.
    @test trie[:b] == 2.
    @test trie[:d] == 3.
    @test trie[:c => :x] == 1
    @test trie[:e => :x] == 1
    @test trie[:f => :x] == 1
    @test trie[:shared => :x] == 1
    @test trie[:shared => :y] == 4.
    @test length(get_internal_nodes(trie)) == 4
    @test length(get_leaf_nodes(trie)) == 3
end

@testset "static assignment merge" begin
    inner = DynamicAssignment()
    set_leaf_node!(inner, :x, 1)
    inner2 = DynamicAssignment()
    set_leaf_node!(inner2, :y, 4.)
    trie1 = StaticAssignment((a=1., b=2.), (c=inner, shared=inner))
    trie2 = StaticAssignment((d=3.,), (e=inner, f=inner, shared=inner2))
    trie = merge(trie1, trie2)
    @test trie[:a] == 1.
    @test trie[:b] == 2.
    @test trie[:d] == 3.
    @test trie[:c => :x] == 1
    @test trie[:e => :x] == 1
    @test trie[:f => :x] == 1
    @test trie[:shared => :x] == 1
    @test trie[:shared => :y] == 4.
    @test length(get_internal_nodes(trie)) == 4
    @test length(get_leaf_nodes(trie)) == 3
end

@testset "dynamic assignment overwrite" begin

    # overwrite leaf node with leaf node
    assignment = DynamicAssignment()
    assignment[:x] = 1
    assignment[:x] = 2
    @test assignment[:x] == 2
    @test !has_internal_node(assignment, :x)

    # overwrite leaf node with an (empty) internal node
    assignment = DynamicAssignment()
    assignment[:x] = 1
    node = DynamicAssignment(); node[:y] = 2
    set_internal_node!(assignment, :x, node)
    @test !has_leaf_node(assignment, :x)
    @test has_internal_node(assignment, :x)

    # overwrite internal node node with a leaf node
    assignment = DynamicAssignment()
    assignment[:x => :y] = 1
    assignment[:x] = 2
    @test !has_internal_node(assignment, :x)
    @test assignment[:x] == 2

    # overwrite internal node node with internal node
    assignment = DynamicAssignment()
    assignment[:x => :y] = 1
    node = DynamicAssignment(); node[:z] = 2
    set_internal_node!(assignment, :x,  node)
    @test has_internal_node(assignment, :x)
    @test !has_leaf_node(assignment, :x => :y)
    @test assignment[:x => :z] == 2

    # illegal set address under leaf node
    assignment = DynamicAssignment()
    assignment[:x] = 1
    threw = false
    try
        assignment[:x => :y] = 2
    catch
        threw = true
    end
    @test threw

    # illegal set address under leaf node
    assignment = DynamicAssignment()
    assignment[:x] = 1
    node = DynamicAssignment(); node[:z] = 2
    threw = false
    try
        set_internal_node!(assignment, :x => :y, node)
    catch
        threw = true
    end
    @test threw
end

ext: jl

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
const TOKEN_PREFIX = '/'

macro j_str(token)
    Pointer(token)
end

"""
    _unescape_jpath(raw::String)

Transform escaped characters in JPaths back to their original value.
https://tools.ietf.org/html/rfc6901
"""
function _unescape_jpath(raw::AbstractString)
    m = match(r"%([0-9A-F]{2})", raw)
    if m !== nothing
        for c in m.captures
            raw = replace(raw, "%$(c)" => Char(parse(UInt8, "0x$(c)")))
        end
    end
    return raw
end

function _last_element_to_type!(jk)
    if !occursin("::", jk[end])
        return Any
    end
    x = split(jk[end], "::")
    jk[end] = String(x[1])
    if x[2] == "string"
        return String
    elseif x[2] == "number"
        return Union{Int, Float64}
    elseif x[2] == "object"
        return OrderedCollections.OrderedDict{String, Any}
    elseif x[2] == "array"
        return Vector{Any}
    elseif x[2] == "boolean"
        return Bool
    elseif x[2] == "null"
        return Missing
    else
        error(
            "You specified a type that JSON doesn't recognize! Instead of " *
            "`::$(x[2])`, you must use one of `::string`, `::number`, " *
            "`::object`, `::array`, `::boolean`, or `::null`."
        )
    end
end

"""
    Pointer(token::AbstractString; shift_index::Bool = false)

A JSON Pointer is a Unicode string containing a sequence of zero or more
reference tokens, each prefixed by a '/' (%x2F) character.

Follows IETF JavaScript Object Notation (JSON) Pointer https://tools.ietf.org/html/rfc6901.

## Arguments

- `shift_index`: shift given index by 1 for compatibility with original JSONPointer.

## Non-standard extensions

- Index numbers starts from `1` instead of `0`

- User can declare type with '::T' notation at the end. For example
  `/foo::string`. The type `T` must be one of the six types supported by JSON:
  * `::string`
  * `::number`
  * `::object`
  * `::array`
  * `::boolean`
  * `::null`

## Examples

    Pointer("/a")
    Pointer("/a/3")
    Pointer("/a/b/c::number")
    Pointer("/a/0/c::object"; shift_index = true)
"""
struct Pointer{T}
    tokens::Vector{Union{String, Int}}
end

function Pointer(token_string::AbstractString; shift_index::Bool = false)
    if startswith(token_string, "#")
        token_string = _unescape_jpath(token_string[2:end])
    end
    if isempty(token_string)
        return Pointer{Nothing}([""])
    end
    if !startswith(token_string, TOKEN_PREFIX)
        throw(ArgumentError("JSONPointer must starts with '$TOKEN_PREFIX' prefix"))
    end
    tokens = convert(
        Vector{Union{String, Int}},
        String.(split(token_string, TOKEN_PREFIX; keepempty = false)),
    )
    if length(tokens) == 0
        return Pointer{Any}([""])
    end
    T = _last_element_to_type!(tokens)
    for (i, token) in enumerate(tokens)
        if occursin(r"^\d+$", token) # index of a array
            tokens[i] = parse(Int, token)
            if shift_index
                tokens[i] += 1
            end
            if iszero(tokens[i])
                throw(ArgumentError("Julia uses 1-based indexing, use '1' instead of '0'"))
            end
        elseif occursin(r"^\\\d+$", token) # literal string for a number
            tokens[i] = String(chop(token; head = 1, tail = 0))
        elseif occursin("~", token)
            tokens[i] = replace(replace(token, "~0" => "~"), "~1" => "/")
        end
    end
    return Pointer{T}(tokens)
end

Base.length(x::Pointer) = length(x.tokens)

Base.eltype(::Pointer{T}) where {T} = T

function Base.show(io::IO, x::Pointer{T}) where {T}
    print(io, "JSONPointer{", T, "}(\"/", join(x.tokens, "/"), "\")")
end

function Base.show(io::IO, ::Pointer{Nothing})
    print(io, "JSONPointer{Nothing}(\"\")")
end

# This code block needs some explaining.
#
# Ideally, one would define methods like Base.haskey(::AbstractDict, ::Pointer).
# However, this causes an ambiguity with Base.haskey(::Dict, key), which has a
# more concrete first argument and a less concrete second argument. We could
# just define both methods to avoid the ambiguity with Dict, but this would
# probably break any package which defines an <:AbstractDict and fails to type
# the second argument to haskey, getindex, etc!
#
# To avoid the ambiguity issue, we have to manually encode each AbstractDict
# subtype that we support :(
for T in (Dict, OrderedCollections.OrderedDict)
    @eval begin
        # This method is used when creating new dictionaries from JSON pointers.
        function $T{K, V}(kv::Pair{<:Pointer, V}...) where {V, K<:Pointer}
            return $T{String, Any}()
        end

        _new_container(::$T) = $T{String, Any}()

        Base.haskey(dict::$T, p::Pointer) = _haskey(dict, p)
        Base.getindex(dict::$T, p::Pointer) = _getindex(dict, p)
        Base.setindex!(dict::$T, v, p::Pointer) = _setindex!(dict, v, p)
        Base.get(dict::$T, p::Pointer, default) = _get(dict, p, default)
    end
end

Base.getindex(A::AbstractArray, p::Pointer) = _getindex(A, p)
Base.haskey(A::AbstractArray, p::Pointer) = _haskey(A, p)

function Base.unique(arr::AbstractArray{<:Pointer, N}) where {N}
    out = deepcopy(arr)
    if isempty(arr)
        return out
    end
    pointers = getfield.(arr, :tokens)
    if allunique(pointers)
        return out
    end
    delete_target = Int[]
    for p in pointers
        indicies = findall(el -> el == p, pointers)
        if length(indicies) > 1
            append!(delete_target, indicies[1:end-1])
        end
    end
    deleteat!(out, unique(delete_target))
    return out
end

Base.:(==)(a::Pointer{U}, b::Pointer{U}) where {U} = a.tokens == b.tokens

# ==============================================================================

_checked_get(collection::AbstractArray, token::Int) = collection[token]

_checked_get(collection::AbstractDict, token::String) = collection[token]

function _checked_get(collection, token)
    error(
        "JSON pointer does not match the data-structure. I tried (and " *
        "failed) to index $(collection) with the key: $(token)"
    )
end

# ==============================================================================

_haskey(::Any, ::Pointer{Nothing}) = true

function _haskey(collection, p::Pointer)
    for token in p.tokens
        if !_haskey(collection, token)
            return false
        end
        collection = _checked_get(collection, token)
    end
    return true
end

_haskey(collection::AbstractDict, token::String) = haskey(collection, token)

function _haskey(collection::AbstractArray, token::Int)
    return 1 <= token <= length(collection)
end

_haskey(::Any, ::Any) = false

# ==============================================================================

_getindex(collection, ::Pointer{Nothing}) = collection

function _getindex(collection, p::Pointer)
    return _getindex(collection, p.tokens)
end

function _getindex(collection, tokens::Vector{Union{String, Int}})
    for token in tokens
        collection = _checked_get(collection, token)
    end
    return collection
end

# ==============================================================================

function _get(collection, p::Pointer, default)
    if _haskey(collection, p)
        return _getindex(collection, p)
    end
    return default
end

# ==============================================================================

_null_value(p::Pointer) = _null_value(eltype(p))
_null_value(::Type{String}) = ""
_null_value(::Type{<:Real}) = 0
_null_value(::Type{<:AbstractDict}) = OrderedCollections.OrderedDict{String, Any}()
_null_value(::Type{<:AbstractVector{T}}) where {T} = T[]
_null_value(::Type{Bool}) = false
_null_value(::Type{Nothing}) = nothing
_null_value(::Type{Missing}) = missing

_null_value(::Type{Any}) = missing

_convert_v(v::U, ::Pointer{U}) where {U} = v
function _convert_v(v::V, p::Pointer{U}) where {U, V}
    v = ismissing(v) ? _null_value(p) : v
    try
        return convert(eltype(p), v)
    catch
        error(
            "$(v)::$(typeof(v)) is not valid type for $(p). Remove type " *
            "assertion in the JSON pointer if you don't a need static type."
        )
    end
end

function _add_element_if_needed(prev::AbstractVector{T}, k::Int) where {T}
    x = k - length(prev)
    if x > 0
        append!(prev, [_null_value(T) for _ = 1:x])
    end
    return
end

function _add_element_if_needed(
    prev::AbstractDict{K, V}, k::String
) where {K, V}
    if !haskey(prev, k)
        prev[k] = _null_value(V)
    end
end

function _add_element_if_needed(collection, token)
    error(
        "JSON pointer does not match the data-structure. I tried (and " *
        "failed) to set $(collection) at the index: $(token)"
    )
end

_new_data(::Any, n::Int) = Vector{Any}(missing, n)
_new_data(::AbstractVector, ::String) = OrderedCollections.OrderedDict{String, Any}()
_new_data(x::AbstractDict, ::String) = _new_container(x)

function _setindex!(collection::AbstractDict, v, p::Pointer)
    prev = collection
    for (i, token) in enumerate(p.tokens)
        _add_element_if_needed(prev, token)
        if i != length(p)
            if ismissing(prev[token])
                prev[token] = _new_data(prev, p.tokens[i + 1])
            end
            prev = prev[token]
        end
    end
    prev[p.tokens[end]] = _convert_v(v, p)
    return v
end


text:
# A type to test unconventional indexing ranges

module ATs  # OffsetArrays

using Base: ReshapedArray, ReshapedArrayIterator
using Base.PermutedDimsArrays: PermutedDimsArray
import ArrayIteration: inds

struct OA{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
end

OA(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) where {T,N} = OA{T,N,typeof(A)}(A, offsets)

Base.parent(A::OA) = A.parent
Base.size(A::OA) = size(parent(A))
inds(A::OA, d) = (1:size(parent(A),d)) .+ A.offsets[d]
eachindexx(A::OA) = CartesianRange(inds(A))

Base.getindex(A::OA, inds::Int...) = parent(A)[offset(A.offsets, inds)...]
Base.setindex!(A::OA, val, inds::Int...) = parent(A)[offset(A.offsets, inds)...] = val

offset(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) where {N} = _offset((), offsets, inds)
_offset(out, ::Tuple{}, ::Tuple{}) = out
@inline _offset(out, offsets, inds) = _offset((out..., inds[1]-offsets[1]), Base.tail(offsets), Base.tail(inds))

# An iterator that deliberately makes PermutedDimsArrays more "dangerous"
# (sync to the rescue!)
struct PDAIterator
    iter::UnitRange{Int}
end
struct PDAIndex
    i::Int
end

Base.parent(A::PermutedDimsArray) = A.parent # move to Base
eachindexx(A::PermutedDimsArray{T,N,AA}) where {T,N,AA<:Array} = PDAIterator(eachindexx(A.parent))

function Base.iterate(iter::PDAIterator, s...)
    v = iterate(iter.iter)
    v === nothing && return nothing
    (i, s) = v
    PDAIndex(i), s
end

Base.getindex(A::PermutedDimsArray, i::PDAIndex) = parent(A)[i.i]
Base.setindex!(A::PermutedDimsArray, val, i::PDAIndex) = parent(A)[i.i] = val

# Turn on eachindex for ReshapedArrays
eachindexx(A::ReshapedArray) = ReshapedArrayIterator(A)

end


text:
struct LexicalException <: Exception
  msg::String
end

@enum TokenClass begin
  ident
  times
  plus
  minus
  divide
  log_not
  log_and
  equals
  less_than
  open_paren
  close_paren
  rng
  semicolon
  colon
  doublequote
  int_literal
  string_literal
  assign
  kw_var
  kw_for
  kw_end
  kw_in
  kw_do
  kw_read
  kw_print
  kw_int
  kw_string
  kw_bool
  kw_assert
  kw_true
  kw_false
  eoi
end

whitespace = [' ', '\t', '\n']
end_of_input_symbol = '$'
delimiters = union(whitespace, end_of_input_symbol)
ident_or_kw_initial = union('a':'z', 'A':'Z')
ident_or_kw_body = union(ident_or_kw_initial, '0':'9', '_')
keywords = Dict([
  "var" => kw_var,
  "for" => kw_for,
  "end" => kw_end,
  "in" => kw_in,
  "do" => kw_do,
  "read" => kw_read,
  "print" => kw_print,
  "int" => kw_int,
  "string" => kw_string,
  "bool" => kw_bool,
  "assert" => kw_assert,
  "true" => kw_true,
  "false" => kw_false
  ])
symbol_initials = ['*', '+', '-', '/', '(', ')', '.', ';', ':', '!', '=', '<', '&']
digits = '0':'9'
unary_ops = Dict(
  log_not => '!'
)
binary_ops = Dict(
  times => '*',
  plus => '+',
  minus => '-',
  divide => '/',
  equals => '=',
  less_than => '<',
  log_and => '&'
  )
syntactic_symbols = Dict(
  open_paren => '(',
  close_paren => ')',
  colon => ':',
  semicolon => ';',
  doublequote => '"'
)
symbol_to_character = union(unary_ops, binary_ops, syntactic_symbols)
character_to_symbol = Dict([(sym => op) for (op, sym) ∈ symbol_to_character])

mutable struct Token
  class::TokenClass
  lexeme::String
  line::Int
end

Token(class::TokenClass, lexeme::String) = Token(class, lexeme, 0)

"""
The main function of the scanner. The "next" index points at the
character to scan next. The outer while loop skips whitespace, keeps
track of the line number and the comment nesting level, calls the
getToken function, stores the resulting tokens and increments the index.
"""
function scanInput(input::AbstractString, next = 1)
  input *= end_of_input_symbol
  tokens = Array{Token,1}()
  lineNumber = 1
  while next <= length(input)
    commentNesting = 0
    while (input[next] in whitespace
      || commentNesting > 0
      || (input[next] == '/' && input[next+1] == '/')
      || (input[next] == '/' && input[next+1] == '*'))
      if input[next] == '\n'
        lineNumber += 1
        next += 1
      elseif input[next] == '/' && input[next+1] == '/' && commentNesting == 0
        while input[next] != '\n' && (next < length(input) - 1) next += 1 end
      elseif input[next] == '/' && input[next+1] == '*'
        commentNesting += 1
        next += 2
      elseif input[next] == '*' && input[next+1] == '/'
        commentNesting -= 1
        next += 2
      else
        next += 1
      end
    end

    if next <= length(input)
      token, next = getToken(input, next, lineNumber)
      token.line = lineNumber
      push!(tokens, token)
    end
  end
  return tokens
end

"""
Looks at the first character of a potential lexeme and
calls the corresponding function to produce a token.
"""
function getToken(input, next, lineNumber)
  c = input[next]
  if c ∈ ident_or_kw_initial
    return getIdentOrKw(input, next, lineNumber)
  end
  if c ∈ symbol_initials
    return getOperator(input, next, lineNumber)
  end
  if c ∈ digits
    return getInteger(input, next, lineNumber)
  end
  if c == '"'
    return getString(input, next, lineNumber)
  end
  if c == '$'
    return Token(eoi, string(end_of_input_symbol)), next+1
  end
  throw(LexicalException(
    "Characted $(c) on or near line $(lineNumber) is not part of a legal token."))
end

function getIdentOrKw(input, next, lineNumber)
  initial = next
  while input[next] ∉ union(delimiters, symbol_initials) next += 1 end
  str = input[initial:next-1]
  token = str ∈ keys(keywords) ? Token(keywords[str], str) : Token(ident, str)
  return token, next
end

function getOperator(input, next, lineNumber)
  if input[next] == '.'
    input[next+1] != '.' && throw(LexicalException(
      "A single dot (on line $(lineNumber)) is not a valid token. Did you mean '..'?"))
    return Token(rng, ".."), next+2
  end
  if input[next] == ':'
    input[next+1] == '=' && return Token(assign, ":="), next+2
    return Token(colon, ":"), next+1
  end
  tokenClass = character_to_symbol[input[next]]
  return Token(tokenClass, string(input[next])), next+1
end

function getInteger(input, next, lineNumber)
  initial = next
  while input[next] ∈ digits next += 1 end
  return Token(int_literal, input[initial:next-1]), next
end

function getString(input, next, lineNumber)
  initial = next
  next += 1
  str = Array{Char,1}()
  while (input[next] != '"' ||
    (input[next] == '"' && input[next-1] == '\\'))
    if input[next] == '\\'
      next += 1
      if input[next] == 'n'
        push!(str, '\n')
      end
      if input[next] == 't'
        push!(str, '\t')
      end
      if input[next] == '"'
        push!(str, '"')
      end
      if input[next] == '\\'
        push!(str, '\\')
      end
      next += 1
    else
      push!(str, input[next])
      next += 1
    end
    next >= length(input) && throw(LexicalException("Reached the end of the program while
      scanning a string literal. Did you forget the closing quote?"))
  end
  return Token(string_literal, join(str)), next+1
end

max_stars_repo_name: juhamyllari/mini-pl-julia

max_issues_repo_name: juhamyllari/mini-pl-julia

max_forks_repo_name: juhamyllari/mini-pl-julia


text:
create(c::Client, ::Type{Guild}; kwargs...) = create_guild(c; kwargs...)

retrieve(c::Client, ::Type{Guild}, guild::Integer) = get_guild(c, guild)
retrieve(c::Client, ::Type{Guild}; kwargs...) = get_current_user_guilds(c; kwargs...)

update(c::Client, g::AbstractGuild; kwargs...) = modify_guild(c, g.id; kwargs...)

delete(c::Client, g::AbstractGuild) = delete_guild(c, g.id)


text:
function Schnittpunkte(Dim,Winkel,X,Y,x0,y0,x1,y1,N,JJ,nx,ny,vx,vy,theta,Nullzeile)
SchnittVektorsort =  zeros(2*nx,JJ*N);
SchnittVektor = [];
HorzBody = zeros(nx+1,2);
VertBody = zeros(nx+1,2);
xS = x0:1:x1;
yS = y0:1:y1;
laengenmatrix=zeros(size(X,1),nx*ny);
Pixeln, Pixelx, Pixely = Pixelberechner(x0,x1,y0,y1,nx,ny);
 for i = 1:size(X,1) # für die Anzahl der parallelen Strahlen durchläuft Startvektoren
     counterhorz = 1;
     countervert = 1;
     VertBody = zeros(nx+1,2);
     HorzBody = zeros(nx+1,2);
     SchnittVektor = [];
     HorzSchnitt = zeros(nx+1,2);
     VertSchnitt = zeros(nx+1,2);
     Abstand = zeros(size(Pixelx,1),1);
     for k = 1:nx+1
        s = (X[i]-xS[k])/vx;
        t = (Y[i]-yS[k])/vy;

        HorzSchnitt[k,1] = X[i]-s*vx; #nächster Schnittpunkt mit Senkrechten, x-Koordinate
        HorzSchnitt[k,2] = Y[i]-s*vy; #nächster Schnittpunkt mit Senkrechten, y-Koordinate
        VertSchnitt[k,1] = X[i]-t*vx; #nächster Schnittpunkt mit Horizontalen, x-Koordinate
        VertSchnitt[k,2] = Y[i]-t*vy; #nächster Schnittpunkt mit Horizontalen, y-Koordinate

        #Überprüfung, dass Schnittpunkte im Körper liegen
        if (round(HorzSchnitt[k,1];digits=5) >= x0) && (round(HorzSchnitt[k,1];digits=5) <= x1) && (round(HorzSchnitt[k,2];digits=5) >= y0) && (round(HorzSchnitt[k,2];digits=5) <= y1)
            HorzBody[counterhorz,1] = HorzSchnitt[k,1];
            HorzBody[counterhorz,2] = HorzSchnitt[k,2];
            counterhorz = counterhorz + 1;
        end
        eps=1E-7;
        if (round(VertSchnitt[k,1];digits=5) >= x0) && (round(VertSchnitt[k,1];digits=5) <= x1) && (round(VertSchnitt[k,2];digits=5) >= y0) && (round(VertSchnitt[k,2];digits=5) <= y1)
            VertBody[countervert,1] = VertSchnitt[k,1];
            VertBody[countervert,2] = VertSchnitt[k,2];
            countervert = countervert+1;
        end
     end
     urHorzBody = HorzBody;
     urVertBody = VertBody;
     for r = 1:size(urVertBody,1);
         if HorzBody[size(urHorzBody,1)+1-r,1] == 0 && HorzBody[size(urHorzBody,1)+1-r,2] == 0
            if VertBody[size(urHorzBody,1)+1-r,2] == 0 && VertBody[size(urHorzBody,1)+1-r,1] == 0
                HorzBody1 = resize(HorzBody[:,1],size(urHorzBody,1)-r);
                HorzBody2 = resize(HorzBody[:,2],size(urHorzBody,1)-r);
                VertBody1 = resize(VertBody[:,1],size(urHorzBody,1)-r);
                VertBody2 = resize(VertBody[:,2],size(urHorzBody,1)-r);
                HorzBody = hcat(HorzBody1,HorzBody2);
                VertBody = hcat(VertBody1,VertBody2);
            else
                HorzBody1 = resize(HorzBody[:,1],size(urHorzBody,1)-r);
                HorzBody2 = resize(HorzBody[:,2],size(urHorzBody,1)-r);
                HorzBody = hcat(HorzBody1,HorzBody2);
            end
         else
            if VertBody[size(urHorzBody,1)+1-r,2] == 0 && VertBody[size(urHorzBody,1)+1-r,1] == 0
                VertBody1 = resize(VertBody[:,1],size(urHorzBody,1)-r);
                VertBody2 = resize(VertBody[:,2],size(urHorzBody,1)-r);
                VertBody = hcat(VertBody1,VertBody2);
            end
         end
     end
     if HorzBody == []
        SchnittVektor = VertBody;
     elseif VertBody == []
        SchnittVektor = HorzBody;
     else
         SchnittVektor = vcat(HorzBody, VertBody);
     end

     if theta <= pi/2
         SchnittVektorsort = sortieren(SchnittVektor,1,0);
     elseif theta == pi
         SchnittVektorsort = sortieren(SchnittVektor,1,1);
     else
         SchnittVektorsort = sortieren(SchnittVektor, 2,0);
     end
     # Berechnung der Mitten zwischen zwei Schnittpunkten und Bestimmung
     # des nächsten Mittelpunkts eines Pixels
     distance = zeros(size(SchnittVektorsort,1)-1,1);
     Mittex= zeros(size(SchnittVektorsort,1)-1,1);
     Mittey = zeros(size(SchnittVektorsort,1)-1,1);
     for l = 1:size(SchnittVektorsort,1)-1
         distance[l] = sqrt((SchnittVektorsort[l,1]-SchnittVektorsort[l+1,1])^2+(SchnittVektorsort[l,2]-SchnittVektorsort[l+1,2])^2);
         Mittex[l,1] = SchnittVektorsort[l,1]+distance[l]/2*vx;
         Mittey[l,1] = SchnittVektorsort[l,2]+distance[l]/2*vy;
         for m = 1:size(Pixelx,1)
             Abstand[m,1] = sqrt((Pixelx[m,1]-Mittex[l,1])^2+(Pixely[m,1]-Mittey[l,1])^2);
         end
         Sortiert= sortieren(Abstand,1,0); # suche Minimum
         Minimum = Sortiert[1,1];
         Ort = finde(Abstand, Minimum);
         x = Int(Ort[1,1]);
         laengenmatrix[i,x] = distance[l];
     end

     if size(SchnittVektorsort,1) == 0
        Nullzeile = Nullzeile+1;
     end
 end
return laengenmatrix,SchnittVektorsort,Nullzeile;
end


text:
using MLJNaiveBayesInterface
using MLJBase
using Test

import StableRNGs
rng = StableRNGs.StableRNG(1234)

## GAUSSIAN

gaussian_classifier = GaussianNBClassifier()

# gaussian classifier takes continuous features
X, y = @load_iris

train, test = partition(eachindex(y), 0.6)

fitresultG, cacheG, reportG = fit(gaussian_classifier, 1,
                    selectrows(X, train), y[train]);

params = fitted_params(gaussian_classifier, fitresultG)
@test params.n_obs == length(train)

gaussian_pred = predict(gaussian_classifier,
                        fitresultG,
                        selectrows(X, test));

yhat1 = gaussian_pred[1]
@test Set(classes(yhat1)) == Set(classes(y[1]))

# test with linear data:
x1 = randn(rng, 3000);
x2 = randn(rng, 3000);
x3 = randn(rng, 3000);
X = (x1=x1, x2=x2, x3=x3);
ycont = x1 - x2 -2x3;
ydiscrete = map(ycont) do η
    η > 0 ? "go" : "stop"
end
y = coerce(ydiscrete, Multiclass)
train, test = partition(eachindex(y), 0.8);

gaussian_classifier = GaussianNBClassifier()

fitresultG, cacheG, reportG = MLJBase.fit(gaussian_classifier, 1,
             selectrows(X, train), y[train])

params = fitted_params(gaussian_classifier, fitresultG)
@test params.n_obs == length(train)

gaussian_pred = MLJBase.predict_mode(gaussian_classifier,
                                     fitresultG, selectrows(X, test))

@test sum(gaussian_pred .!= y[test])/length(y) < 0.05

## MULTINOMIAL

# first contrive some test data

# some word counts in children's books about colours:
red = [2, 0, 1, 0, 1]
blue = [4, 1, 2, 3, 2]
green = [0, 2, 0, 6, 1]
X = (red=red, blue=blue, green=green)

# gender of author:
y = categorical(['m', 'f', 'm', 'f', 'm'])
male = y[1]
female = y[2]

# Note: The smoothing algorithm is to add to the training data, for
# each class observed, a row with every feature getting count of
# alpha. So smoothing also effects the class marginals (is this
# standard?) Only integer values of alpha allowed.

# computing conditional probabilities by hand with Lagrangian
# smoothing (alpha=1):
red_given_m = 5/16
blue_given_m = 9/16
green_given_m = 2/16
red_given_f = 1/15
blue_given_f = 5/15
green_given_f = 9/15

m_(red, blue, green) =
    4/7*(red_given_m^red)*(blue_given_m^blue)*(green_given_m^green)
f_(red, blue, green) =
    3/7*(red_given_f^red)*(blue_given_f^blue)*(green_given_f^green)
normalizer(red, blue, green) = m_(red, blue, green) + f_(red, blue, green)
m(a...) = m_(a...)/normalizer(a...)
f(a...) = f_(a...)/normalizer(a...)

Xnew = (red=[1, 1], blue=[1, 2], green=[1, 3])

# prediction by hand:

yhand =[MLJBase.UnivariateFinite([male, female], [m(1, 1, 1), f(1, 1, 1)]),
        MLJBase.UnivariateFinite([male, female], [m(1, 2, 3), f(1, 2, 3)])]

multinomial_classifier = MultinomialNBClassifier()

fitresultMLT, cacheMLT, reportMLT =
    MLJBase.fit(multinomial_classifier, 1, X, y)

params = fitted_params(multinomial_classifier, fitresultMLT)
@test params.n_obs == length(y) + 3 * 2 # 3 features, 2 classes

yhat = MLJBase.predict(multinomial_classifier, fitresultMLT, Xnew)

# see issue https://github.com/dfdx/NaiveBayes.jl/issues/42
@test pdf(yhand[1], 'm') ≈ pdf(yhat[1], 'm')
@test pdf(yhand[1], 'f') ≈ pdf(yhat[1], 'f')


text:
function absolute_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(abs.(y_true .- y_pred))
end

function percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return((y_true .- y_pred) ./ y_true)
end

function log_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(log(y_true - y_pred))
end

function squared_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return((y_true - y_pred) .^ 2)
end

function squared_log_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(log_error(y_true, y_pred) .^ 2)
end

function absolute_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(abs.(percent_error(y_true, y_pred)))
end

function mean_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(y_true - y_pred))
end

function mean_absolute_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(absolute_error(y_true, y_pred)))
end

function median_absolute_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(median(absolute_error(y_true, y_pred)))
end

function mean_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(percent_error(y_true, y_pred)))
end

function median_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(median(percent_error(y_true, y_pred)))
end

function mean_squared_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(squared_error(y_true, y_pred)))
end

function median_squared_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(median(squared_error(y_true, y_pred)))
end

function sum_squared_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(sum(squared_error(y_true, y_pred)))
end

function mean_squared_log_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(squared_log_error(y_true, y_pred)))
end

function mean_absolute_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(absolute_percent_error(y_true, y_pred)))
end

function median_absolute_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(median(absolute_percent_error(y_true, y_pred)))
end

function symmetric_mean_absolute_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(mean(abs.(y_pred .- y_true) ./ ((abs.(y_true) .+ abs.(y_pred)) ./ 3)))
end

function symmetric_median_absolute_percent_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(median(abs.(y_pred - y_true) ./ ((abs.(y_true) .+ abs.(y_pred)) ./ 2)))
end

function mean_absolute_scaled_error(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    n = max(length(y_true), length(y_pred))
    numerator = sum(abs.(y_true .- y_pred))
    denominator = (n / (n - 1)) * sum(abs.(y_true[2:n] .- y_pred[1:(n-1)]))
    return(numerator / denominator)
end

function total_variance_score(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(sum((y_true .- mean(y_true)) .^ 2))
end

function explained_variance_score(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(sum((y_pred .- mean(y_true)) .^ 2))
end

function unexplained_variance_score(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(sum((y_true .- y_pred) .^ 2))
end

function r2_score(y_true, y_pred)
    @_dimcheck length(y_true) == length(y_pred)
    return(explained_variance_score(y_true, y_pred) / total_variance_score(y_true, y_pred))
end


text:
using Cairo, Compose

function plot_network(graph::AbstractGraph, user_paths, locs_x, locs_y)
    used_edges = []
    used_by = Dict()

    colour_pal = [colorant"lightgrey", colorant"orange", colorant"lightslateblue", colorant"green"]

    # nodecolor = ["blue" for i in 1:nv(graph)]

    this_user = 0
    for paths in user_paths
        this_user += 1
        for path in paths
            for edge in path
                push!(used_edges, edge)
                used_by[(edge.src,edge.dst)] = this_user
            end
        end
    end

    colours = []
    widths = []
    for edge in edges(graph)
        if edge in used_edges
            push!(colours, used_by[(edge.src,edge.dst)] + 1)
            push!(widths, 5)
        else
            push!(colours, 1)
            push!(widths, 1)
        end
    end

    mygplot = gplot(graph, locs_x, locs_y, edgestrokec=colour_pal[colours],
    edgelinewidth=widths, arrowlengthfrac=0.04)#, layout=spring_layout)
    # Save to pdf
    draw(PDF("plots/network_drawing.pdf", 16cm, 16cm), mygplot)
end

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
    fix!(model::JuMP.Model, solution::AbstractDict)::Nothing

Fix the value of all binary variables to the ones specified by the given
solution. Useful for computing LMPs.
"""
function fix!(model::JuMP.Model, solution::AbstractDict)::Nothing
    instance, T = model[:instance], model[:instance].time
    is_on = model[:is_on]
    prod_above = model[:prod_above]
    reserve = model[:reserve]
    for g in instance.units
        for t in 1:T
            is_on_value = round(solution["Is on"][g.name][t])
            prod_value =
                round(solution["Production (MW)"][g.name][t], digits = 5)
            reserve_value =
                round(solution["Reserve (MW)"][g.name][t], digits = 5)
            JuMP.fix(is_on[g.name, t], is_on_value, force = true)
            JuMP.fix(
                prod_above[g.name, t],
                prod_value - is_on_value * g.min_power[t],
                force = true,
            )
            JuMP.fix(reserve[g.name, t], reserve_value, force = true)
        end
    end
    return
end


text:
module AcousticFeatures

using AlphaStableDistributions, DSP, FindPeaks1D, LinearAlgebra, Statistics, StatsBase, ProgressMeter

include("subsequences.jl")
include("utils.jl")

export Energy, Myriad, VMyriad, FrequencyContours, SoundPressureLevel, ImpulseStats, AlphaStableStats, Entropy, Score

export Subsequence

export spectrumflatten, myriadconstant, vmyriadconstant, pressure, envelope

abstract type AbstractAcousticFeature end

################################################################################
#
#   AcousticFeature types
#
################################################################################
struct Energy <: AbstractAcousticFeature end

struct Myriad{T<:Union{Nothing,Real}} <: AbstractAcousticFeature
    sqKscale::T
end
Myriad() = Myriad{Nothing}(nothing)

struct VMyriad{T<:Union{Nothing,Real},M<:Union{Nothing,AbstractMatrix}} <: AbstractAcousticFeature
    K²::T
    Rₘ::M
end
VMyriad() = VMyriad(nothing, nothing)

struct FrequencyContours{FT<:Real,T<:Real} <: AbstractAcousticFeature
    fs::FT
    n::Int
    tnorm::Union{Nothing,T} #time constant for normalization (sec)
    fd::T #frequency difference from one step to the next (Hz)
    minhprc::T
    minfdist::T
    mintlen::T
end

struct SoundPressureLevel{T<:Real} <: AbstractAcousticFeature
    ref::T
end
SoundPressureLevel() = SoundPressureLevel(1.0)

struct ImpulseStats{FT<:Real,T<:Real} <: AbstractAcousticFeature
    fs::FT
    k::Int
    tdist::T
    computeenvelope::Bool
end
ImpulseStats(fs) = ImpulseStats(fs, 10, 1e-3, true)
ImpulseStats(fs, k, tdist) = ImpulseStats(fs, k, tdist, true)

struct AlphaStableStats <: AbstractAcousticFeature end

struct Entropy{FT<:Real} <: AbstractAcousticFeature
    n::Int
    noverlap::Int
    fs::FT
end

mutable struct Score{VT1<:AbstractArray{<:Real},VT2<:AbstractRange{Int}}
    s::VT1
    indices::VT2
end
################################################################################
#
#   Implementations
#
################################################################################

outputndims(::Energy) = 1
outputeltype(::Energy) = Float64
"""
    Score of `x` based on mean energy.
"""
score(::Energy, x::AbstractVector{T}) where T<:Real = [mean(abs2, x)]

outputndims(::Myriad{S}) where S<:Real = 1
outputeltype(::Myriad{S}) where S<:Real = Float64
"""
    Score of `x` based on myriad.

    Reference:
    Mahmood et. al., "Optimal and Near-Optimal Detection in Bursty Impulsive Noise,"
    IEEE Journal of Oceanic Engineering, vol. 42, no. 3, pp. 639--653, 2016.
"""
function score(f::Myriad{S}, x::AbstractVector{T}) where {T<:Real, S<:Real}
    sqKscale = f.sqKscale
    [sum(x -> log(sqKscale + abs2(x)), x)]
end

score(f::Myriad{Nothing}, x) = score(Myriad(myriadconstant(x)), x)

outputndims(::VMyriad{S,M}) where {S<:Real,M<:AbstractMatrix} = 1
outputeltype(::VMyriad{S,M}) where {S<:Real,M<:AbstractMatrix} = Float64
"""
    Score of `x` based on vector myriad.

    Reference:
    Mahmood et. al., "Optimal and Near-Optimal Detection in Bursty Impulsive Noise,"
    IEEE Journal of Oceanic Engineering, vol. 42, no. 3, pp. 639--653, 2016.
"""
function score(f::VMyriad{S,M}, x::AbstractVector{T}) where {T<:Real,S<:Real,M<:AbstractMatrix}
    K², Rₘ = f.K², f.Rₘ
    m = size(Rₘ, 1)-1
    N = length(x)
    N < m && throw(ArgumentError("`m` has to be larger than length of `x`"))
    mplusonedivtwo = (m+1)/2
    mplustwodivtwo = (m+2)/2
    onetom = 1:m
    invRₘ = inv(cholesky(Rₘ).L)
    invRₘonetom = inv(cholesky(Rₘ[onetom, onetom]).L)
    s = mplusonedivtwo*log(K²+norm²(invRₘonetom*x[onetom]))
    for n in m+1:N
        s += @views mplustwodivtwo*log(K²+norm²(invRₘ*x[n-m:n])) - mplusonedivtwo*log(K²+norm²(invRₘonetom*x[n-m:n-1]))
    end
    s
end
score(f::VMyriad{Nothing,Nothing}, x) = score(VMyriad(vmyriadconstant(x)...), x)

outputndims(::FrequencyContours) = 1
outputeltype(::FrequencyContours) = Float64
"""
    Score of `x` based on frequency contours count.

    Reference:
    D. Mellinger, R. Morrissey, L. Thomas, J. Yosco, "A method for detecting whistles, moans, and other frequency
    contour sounds", 2011 J. Acoust. Soc. Am. 129 4055
"""
function score(f::FrequencyContours, x::AbstractVector{T}) where T<:Real
    spec = spectrogram(x, f.n, f.n÷2; fs=f.fs, window=DSP.hamming)
    p  = spec.power; frequency=spec.freq; t=spec.time
    δt = t[2]-t[1]
    δf = frequency[2]-frequency[1]
    f.tnorm === nothing ? Nnorm = size(p, 2) : Nnorm = f.tnorm÷(δt) |> Int
    p    = spectrumflatten(p, Nnorm) #noise-flattened spectrogram
    # crds,_ = @views peakprom(p[:, 1], Maxima(), trunc(Int, f.minfdist÷δf), eps(T)+percentile(p[:, 1], f.minhprc))
    crds, _ = findpeaks1d(p[:, 1]; height=eps(T)+percentile(p[:, 1], f.minhprc), distance=trunc(Int, f.minfdist/δf))
    ctrs = [[(crd, 1)] for crd in crds]
    for (i, col) in enumerate(eachcol(p[:, 2:end]))
        col = collect(col)
        # crds,_ = Peaks.peakprom(col, Maxima(), trunc(Int, f.minfdist/δf), eps(T)+percentile(col, f.minhprc))
        crds, _ = findpeaks1d(col; height=eps(T)+percentile(col, f.minhprc), distance=trunc(Int, f.minfdist/δf))
        for crd in crds
            if length(ctrs) == 0
                ctrs = [[(crd, 1)] for crd in crds]
            else
                idxselect = Int64[]
                costselect = Float64[]
                for (j, ctr) in enumerate(ctrs)
                    if (ctr[end][2] == i-1) && abs(frequency[ctr[end][1]]-frequency[crd]) <= f.fd
                        push!(idxselect, j)
                        push!(costselect, abs(frequency[ctr[end][1]]-frequency[crd]))
                    end
                end
                if isempty(idxselect)
                    push!(ctrs, [(crd, i)])
                else
                    idxopt = idxselect[argmin(costselect)]
                    push!(ctrs[idxopt], (crd, i))
                end
            end
        end
    end
    idxdelete = Int64[]
    for (i, ctr) in enumerate(ctrs)
        (length(ctr)-1)*(δt) < f.mintlen && push!(idxdelete, i)
    end
    deleteat!(ctrs, idxdelete)
    count = isempty(ctrs) ? 0 : sum(length, ctrs)
    [count/length(p)]
end

outputndims(::SoundPressureLevel) = 1
outputeltype(::SoundPressureLevel) = Float64
"""
Score of `x` based on Sound Pressure Level (SPL). `x` is in micropascal. In water, the common reference is 1 micropascal. In air, the common reference is 20 micropascal.
"""
function score(f::SoundPressureLevel, x::AbstractVector{T}) where T<:Real
    rmsx = sqrt(mean(abs2, x))
    [20*log10(rmsx/f.ref)]
end

outputndims(::ImpulseStats) = 3
outputeltype(::ImpulseStats) = Float64
"""
Score of `x` based on number of impulses, mean and variance of inter-impulse intervals. The minimum height of impulses is defined by `a+k*b` where `a` is median of the envelope of `x` and `b` is median absolute deviation (MAD) of the envelope of `x`.
"""
function score(f::ImpulseStats, x::AbstractVector{T}) where T<:Real
    if f.computeenvelope
        x = envelope(x)
    end
    center = Statistics.median(x)
    height = center+f.k*mad(x, center=center, normalize=false)
    distance = trunc(Int, f.tdist*f.fs)
    # crds, _ = Peaks.peakprom(x, Maxima(), distance, height)
    crds,_ = findpeaks1d(x; height=height, distance=distance)
    timeintervals = diff(crds)
    [length(crds) mean(timeintervals)/f.fs var(timeintervals)/f.fs]
end

outputndims(::AlphaStableStats) = 2
outputeltype(::AlphaStableStats) = Float64
"""
Score of `x` based on the parameters of Symmetric Alpha Stable Distributions. The parameter α measures the impulsiveness while the parameter scale measures the width of the distributions.
"""
function score(f::AlphaStableStats, x::AbstractVector{T}) where T<:Real
    d = fit(AlphaStable, x)
    [d.α d.scale]
end

outputndims(::Entropy) = 3
outputeltype(::Entropy) = Float64
"""
Score of `x` based on temporal entropy, spectral entropy and entropy index.

Reference:
J. Sueur, A. Farina, A. Gasc, N. Pieretti, S. Pavoine, Acoustic Indices for Biodiversity Assessment and Landscape Investigation, 2014.
"""
function score(f::Entropy, x::AbstractVector{T}) where T<:Real
    sp = spectrogram(x, f.n, f.noverlap; fs=f.fs).power

    ne = normalize_envelope(x)
    n = length(ne)
    Ht = -sum(ne .* log2.(ne)) ./ log2(n)
    ns = normalize_spectrum(sp)
    N = length(ns)
    Hf = -sum(ns .* log2.(ns)) ./ log2.(N)
    H = Ht*Hf
    [Ht Hf H]
end

function Score(f::AbstractAcousticFeature,
               x::AbstractVector{T};
               winlen::Int=length(x),
               noverlap::Int=0,
               subseqtype::DataType=Float64,
               preprocess::Function=x->x,
               map::Function=map,
               showprogress::Bool=true) where {T<:Real, N, L}
    xlen = length(x)
    if winlen < xlen
        (noverlap < 0) && throw(ArgumentError("`noverlap` must be larger or equal to zero."))
        subseqs = Subsequence(x, winlen, noverlap)
#        sc = Score(zeros(outputeltype(f), length(subseqs), outputndims(f)), 1:subseqs.step:xlen)
    elseif winlen == xlen
        stmp = score(f, preprocess(convert.(subseqtype, x)))
        if stmp isa Number
            return Score(reshape([stmp], (1, 1)), 1:1)
        else
            return Score(stmp, 1:1)
        end
    else
        throw(ArgumentError("`winlen` must be smaller or equal to the length of `x`."))
    end
    if showprogress
        s = @showprogress map(x -> score(f, preprocess(convert.(subseqtype, x))), subseqs)
    else
        s = map(x -> score(f, preprocess(convert.(subseqtype, x))), subseqs)
    end
    Score(reshape(vcat(s...), (length(s), length(s[1]))), 1:subseqs.step:xlen)
    # @inbounds for (i, subseq) in enumerate(subseqs)
    #     sc.s[i, :] = score(f, preprocess(convert.(subseqtype, subseq)))
    # end
    # sc
end

(f::AbstractAcousticFeature)(x) = Score(f, x)

end


text:
module P4est

using CBinding: @cstruct
using Reexport: @reexport

include("libp4est.jl")

end