text:
using BasisFunctionExpansions
using Test, DSP, LinearAlgebra, Statistics, Random
Random.seed!(0)

# write your own tests here
@test BasisFunctionExpansions.get_centers_automatic(1:10,5,false)[1] |> length == 5
@test BasisFunctionExpansions.get_centers_automatic(1:10,5,true)[1] |> length == 10

@test BasisFunctionExpansions.get_centers_automatic([1:10 2:11],[5,6])[1] |> size == (2,5*6)

rbf = UniformRBFE(1:3, 5)
@test size(rbf(randn(10))) == (10,5)
@test size(rbf(randn())) == (5,)

x = range(0, stop=4, length=50)
a = rbf(x)
@test a[1] == a[end] # Should be completely symmetric in this case
@test isa(rbf, BasisFunctionExpansion)

v   = [1:3 2:4]
Nv  = [2,3]
rbf = MultiUniformRBFE(v,Nv)
@test rbf(randn(2)) |> size == (prod(Nv),)
@test rbf(randn(10,2)) |> size == (10,prod(Nv))

# Single dim
N    = 1000
v    = range(0, stop=10, length=N)
y    = randn(N)
y    = filt(ones(500)/500,[1],y)
Nv   = 10
rbf  = UniformRBFE(v,Nv, normalize=true)
bfa  = BasisFunctionApproximation(y,v,rbf,1)
yhat = bfa(v)
e = y-yhat
@test √(mean(e.^2)) < 0.012

# Multidim
N    = 1000
x    = range(0, stop=2pi-0.2, length=N)
v    = [cos.(x) sin.(x)].*x
y    = randn(N)
y    = filt(ones(500)/500,[1],y)
Nv   = [10,10]
rbf  = MultiUniformRBFE(v,Nv, normalize=true)
bfa  = BasisFunctionApproximation(y,v,rbf,0.0001)
yhat = bfa(v)
e = y-yhat
@test isapprox.(sum(rbf(randn(10,2)), dims=2), 1, atol=1e-7) |> all
@test √(mean(e.^2)) < 0.08
# plot(rbf)

# Multidim Diagonal
N    = 1000
x    = range(0, stop=2pi-0.2, length=N)
v    = [cos.(x) sin.(x)].*x
y    = randn(N)
y    = filt(ones(500)/500,[1],y)
Nc   = 8
rbf  = BasisFunctionExpansions.MultiDiagonalRBFE(v,Nc, normalize=true)
bfa  = BasisFunctionApproximation(y,v,rbf,0.0001)
yhat = bfa(v)
e = y-yhat
√(mean(e.^2))

# @test isapprox.(sum(rbf(randn(10,2)), dims=2), 1, atol=1e-7) |> all
@test √(mean(e.^2)) < 0.02

# Multidim
# using BasisFunctionExpansions, Plots
N    = 1000
x    = range(0, stop=2pi-0.2, length=N)
v    = [cos.(x) sin.(x)].*x
y    = randn(N)
y    = filt(ones(500)/500,[1],y)
Nc   = 8
rbf  = MultiRBFE(v,Nc, normalize=true)
bfa  = BasisFunctionApproximation(y,v,rbf,0.0001)
yhat = bfa(v)
e = y-yhat
√(mean(e.^2))

# @test isapprox.(sum(rbf(randn(10,2)), dims=2), 1, atol=1e-7) |> all
@test √(mean(e.^2)) < 0.02
# scatter3d(v[:,1],v[:,2],y, lab="Signal")
# scatter3d!(v[:,1],v[:,2],yhat, lab="Reconstruction")

# Dynamics
# using BasisFunctionExpansions
# using Plots; plotly()
T = 100
A = [1,2*0.7*1,1]
B = [10,5]
u = 1randn(T)
y = filt(B,A,sqrt.(abs.(u)))

yr,A = getARXregressor(y,u,3,2)

rbf = MultiUniformRBFE(A,[2,2,4,4,4], normalize=true)
bfa = BasisFunctionApproximation(yr,A,rbf, 1e-4)
e = √(mean((yr - bfa(A)).^2))
x = A\yr
e2 = √(mean((yr - A*x).^2))

# plot([yr bfa(A) A*x]); gui()

@test e < e2
@test e < 0.04

include("test_LPVSS.jl")


text:
function V=Rosenbrock(V,dt,Fcn,Jac,CG,Param)
Vn=V;
ROS=Param.ROS;
nV1=size(V,1);
nV2=size(V,2);
nV3=size(V,3);
nJ=nV1*nV2*nV3;
nStage=ROS.nStage;
k=zeros([size(V) nStage]);
J=Jac(V,CG,Param);
[L,U,p,q]=lu((1/(dt*ROS.d))*speye(nJ)-J,'vector');
for iStage=1:nStage
  V=Vn;
  for jStage=1:iStage-1
    V=V+ROS.a(iStage,jStage)*k(:,:,:,jStage);
  end
  fV(:,:,:)=Fcn(V,CG,Param);
  for jStage=1:iStage-1
    fV=fV+(ROS.c(iStage,jStage)/(ROS.d*dt))*k(:,:,:,jStage);
  end
  v=reshape(permute(fV,[2 1 3]),nJ,1);
  v(q)=U\(L\v(p));
  k(:,:,:,iStage)=permute(reshape(v,nV2,nV1,nV3),[2 1 3]);
end
V=Vn;
for iStage=1:nStage
  V=V+ROS.m(iStage)*k(:,:,:,iStage);
end
end

ext: jl


text:
firstPayment = 3000
i = 0.12
m = 12

# (A/P, i, N)
function capitalRecoveryFactor(i, N)
	(i * (1 + i)^N) / ((i+1)^N - 1)
end

# Find: How much the guy paid the student
P = firstPayment / (1 + i/m)^(m * 1/12)
A = P * capitalRecoveryFactor(i/m, 12)
C = A * 1 / capitalRecoveryFactor(i/m, 5) + A + 1000

println("Guy paid: ", C)


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03xz(datfile, io=stdout)
    ZER = 0.0e0
    ZERO = 0.0+0.0im
    ONE = 1.0+0.0im
    NIN = 5
    NOUT = 6
    NMAX = 100
    LDA   = 2*NMAX
    LDAE  = 2*NMAX
    LDAS  = NMAX
    LDGE  = 2*NMAX
    LDQG  = 2*NMAX
    LDQGE = 2*NMAX
    LDQGS  = NMAX
    LDRES = 2*NMAX
    LDU1  = 2*NMAX
    LDU2   = 2*NMAX
    LDWORK = 20*NMAX*NMAX + 12*NMAX + 2
    LZWORK = 12*NMAX - 2
    A = Array{ComplexF64,2}(undef, LDA,2*NMAX)
    QG = Array{ComplexF64,2}(undef, LDQG,2*NMAX+1)
    U1 = Array{ComplexF64,2}(undef, LDU1,2*NMAX)
    U2 = Array{ComplexF64,2}(undef, LDU2,2*NMAX)
    WR = Array{Float64,1}(undef, 2*NMAX)
    WI = Array{Float64,1}(undef, 2*NMAX)
    SCALE = Array{Float64,1}(undef, NMAX)
    BWORK = Array{BlasBool,1}(undef, 2*NMAX)
    AS = Array{ComplexF64,2}(undef, LDAS,NMAX)
    QGS = Array{ComplexF64,2}(undef, LDQGS,NMAX+1)
    DWORK = Array{Float64,1}(undef, LDWORK)
    QGE = Array{ComplexF64,2}(undef, LDQGE,2*NMAX+1)
    GE = Array{ComplexF64,2}(undef, LDGE,2*NMAX)
    AE = Array{ComplexF64,2}(undef, LDAE,2*NMAX)
    RES = Array{ComplexF64,2}(undef, LDRES,2*NMAX)
    ZWORK = Array{ComplexF64,1}(undef, LZWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    BALANC =  vs[2][1]
    JOB =  vs[3][1]
    JOBU =  vs[4][1]
    if ( N<=0 || N>NMAX )
        @error "Illegal N=$N"
    end

    M = 2*N
    A[1:M,1:M] .= ZERO
    QG[1:M,1:2M] .= ZERO

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
    end
    AS[1:N,1:N] .= A[1:N,1:N]

    vs = String[]
    _isz,_jsz = (N,N+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       QG[i,1:_jsz] .= parsex.(ComplexF64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 1

    #FOREIGN.zlacpy!( 'All', N, N+1, QG, LDQG, QGS, LDQGS )
    QGS[1:N,1:N+1] .= QG[1:N,1:N+1]

# interp call 2

    ILO, INFO = SLICOT.mb03xz!(BALANC, JOB, JOBU, N, A, QG, U1, U2, WR, WI, SCALE, BWORK)
    println(io, "ILO = $ILO")
    @test INFO == 0
    INFO == 0 || return

    println(io, "eigvals:")
    show(io, "text/plain", WR[1:M] + im*WI[1:M])
    println(io)
    if ( LSAME( JOB, 'S' )||LSAME( JOB, 'G' ) )
# interp output 1
    println(io, "A:")
    _nc = M
    _nr = M
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

    end # if
    if ( LSAME( JOB, 'G' ) )
# interp output 2
    println(io, "QG:")
    _nc = M
    _nr = M
    show(io, "text/plain", QG[1:_nr,1:_nc])
    println(io)

    end # if
    if ( LSAME( JOB, 'G' )&&LSAME( JOBU, 'U' ) )
# interp call 3

        ILO, INFO = SLICOT.mb04dz!(BALANC, N, AS, QGS, DWORK)
        println(io, "ILO = $ILO")
        @test INFO == 0
    INFO == 0 || return

# interp call 4

        #FOREIGN.zlaset!( 'Lower', M-1, M-1, ZERO, ZERO, A(2,1), LDA )
        triu!(view(A,1:M,1:M-1))

# interp call 5

        SLICOT.ma02ez!('U', 'C', 'N', M, QG)

        # compute Ae,Ge,Qe
        AE[1:N,1:N] .= im * imag.(AS[1:N,1:N])
        AE[N+1:2N,1:N] .= -im * real.(AS[1:N,1:N])
        for j in 1:N
            AE[1:N,j+N] .= -AE[N+1:2N,j]
            AE[N+1:2N,j+N] .= AE[1:N,j]
        end
        QGE[1:N,1:N+1] .= im * imag.(QGS[1:N,1:N+1])
        for j in 1:N+1
            QGE[N+j:2N,j] .= -im * real.(QGS[j:N,j])
        end
        for j in 1:N
            QGE[1:j,j+N+1] .= -QGE[1:j,j+1]
            QGE[N+1:2N,j+N+1] .= QGE[1:N,j+1]
        end
        QGE[N+1:2N,N+1] .= QGE[1:N,1]

# interp call 7

        SLICOT.ma02ez!('L', 'T', 'N', N, view(QGE,N+1:2N,1:N+1))

# interp call 8

        SLICOT.ma02ez!('U', 'T', 'N', N, view(QGE,1:N,N+2:2N+2))

# interp call 9

        # FOREIGN.zlacpy!( 'Upper', M, M, QGE(1,2), LDQGE, GE, LDGE )
        GE[1:M,1:M] .= triu(QGE[1:M,2:M+1])

# interp call 10

        SLICOT.ma02ez!('U', 'T', 'S', M, GE)

# interp call 11

        SLICOT.ma02ez!('L', 'T', 'S', M, QGE)

# interp call 12

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, AE, LDAE, U1, LDU1, ZERO, RES, LDRES )
        BLAS.gemm!('N','N',true,AE[1:M,1:M],U1[1:M,1:M],false,view(RES,1:M,1:M))

# interp call 13

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, -ONE, GE, LDGE, U2, LDU2, ONE, RES, LDRES )
        BLAS.gemm!('N','N',-ONE,GE[1:M,1:M],U2[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 14

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, -ONE, U1, LDU1, A, LDA, ONE, RES, LDRES )
        BLAS.gemm!('N','N',-ONE,U1[1:M,1:M],A[1:M,1:M],true,view(RES,1:M,1:M))

        # TEMP = ZLANGE( 'Frobenius', M, M, RES, LDRES, DWORK )
        TEMP = norm(RES[1:M,1:M])
# interp call 15

        #FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, AE, LDAE, U2, LDU2, ZERO, RES, LDRES )
        BLAS.gemm!('N','N',true,AE[1:M,1:M],U2[1:M,1:M],false,view(RES,1:M,1:M))

# interp call 16

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, GE, LDGE, U1, LDU1, ONE, RES, LDRES )
        BLAS.gemm!('N','N',true,GE[1:M,1:M],U1[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 17

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, -ONE, U1, LDU1, QG, LDQG, ONE, RES, LDRES )
        BLAS.gemm!('N','N',-ONE,U1[1:M,1:M],QG[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 18

        # FOREIGN.zgemm!( 'No Transpose', 'Conj Transpose', M, M, M, ONE, U2, LDU2, A, LDA, ONE, RES, LDRES )
        BLAS.gemm!('N','C',true,U2[1:M,1:M],A[1:M,1:M],true,view(RES,1:M,1:M))

        # TEMP = DLAPY2( TEMP, ZLANGE( 'Frobenius', M, M, RES, LDRES, DWORK ) )
        TEMP = hypot(TEMP,norm(RES[1:M,1:M]))
# interp call 19

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, true, QGE, LDQGE, U1, LDU1, ZERO, RES, LDRES )
        BLAS.gemm!('N','N',true,QGE[1:M,1:M],U1[1:M,1:M],false,view(RES,1:M,1:M))

# interp call 20

        # FOREIGN.zgemm!( 'Transpose', 'No Transpose', M, M, M, -ONE, AE, LDAE, U2, LDU2, ONE, RES, LDRES )
        BLAS.gemm!('T','N',-ONE,AE[1:M,1:M],U2[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 21

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, U2, LDU2, A, LDA, ONE, RES, LDRES )
        BLAS.gemm!('N','N',true,U2[1:M,1:M],A[1:M,1:M],true,view(RES,1:M,1:M))

        # TEMP = DLAPY2( TEMP, ZLANGE( 'Frobenius', M, M, RES, LDRES, DWORK ) )
        TEMP = hypot(TEMP,norm(RES[1:M,1:M]))
# interp call 22

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, QGE, LDQGE, U2, LDU2, ZERO, RES, LDRES )
        BLAS.gemm!('N','N',true,QGE[1:M,1:M],U2[1:M,1:M],false,view(RES,1:M,1:M))

# interp call 23

        # FOREIGN.zgemm!( 'Transpose', 'No Transpose', M, M, M, ONE, AE, LDAE, U1, LDU1, ONE, RES, LDRES )
        BLAS.gemm!('T','N',true,AE[1:M,1:M],U1[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 24

        # FOREIGN.zgemm!( 'No Transpose', 'No Transpose', M, M, M, ONE, U2, LDU2, QG, LDQG, ONE, RES, LDRES )
        BLAS.gemm!('N','N',true,U2[1:M,1:M],QG[1:M,1:M],true,view(RES,1:M,1:M))

# interp call 25

        # FOREIGN.zgemm!( 'No Transpose', 'Conj Transpose', M, M, M, ONE, U1, LDU1, A, LDA, ONE, RES, LDRES )
        BLAS.gemm!('N','C',true,U1[1:M,1:M],A[1:M,1:M],true,view(RES,1:M,1:M))

        #TEMP = DLAPY2( TEMP, ZLANGE( 'Frobenius', M, M, RES, LDRES, DWORK ) )
        TEMP = hypot(TEMP,norm(RES[1:M,1:M]))
        println(io, "residual norm = $TEMP")
    end # if
    if ( !LSAME( JOB, 'E' )&&LSAME( JOBU, 'U' ) )
        # interp output 3
        println(io, "U1:")
        _nc = M
        _nr = M
        show(io, "text/plain", U1[1:_nr,1:_nc])
        println(io)
        println(io, "U2:")
        _nc = M
        _nr = M
        show(io, "text/plain", U1[2:_nr,1:_nc])
        println(io)

        # unable to translate write statement:
        #  write MA02JZ( .FALSE., .FALSE., M, U1, LDU1, U2, LDU2, RES, LDRES )
        resid = SLICOT.ma02jz!(false,false,M,U1,U2)
        println(io, "U residual = $resid")
    end # if
    if ( LSAME( BALANC, 'S' )||LSAME( BALANC, 'B' ) )
        println(io, "SCALE:")
        show(io, "text/plain", SCALE[1:N])
        println(io)
    end # if
end # run_mb03xz()


text:
###
### TIndex
###

struct TIndex{N,ipos,T<:Tuple}
    I::T
    TIndex{N,ipos}(t::T) where {N,ipos,T<:Tuple} = new{N,ipos,T}(t)
end

TIndex{N,ipos}(t...) where {N,ipos} = TIndex{N,ipos}(t)

function Base.show(io::IO, ti::TIndex{N,ipos}) where {N,ipos}
    print(io, "TIndex{", N, ",", ipos, "}", ti.I)
end

###
### (::Tensor)[::TIndex]( = x)
###

@inline function Base.to_indices(
    t::Tensor{<:Any,N,ipos},
    inds,
    I::Tuple{TIndex{N,ipos}}
) where {N,ipos}
    #to_indices(t, inds, first(I).I)
    to_indices(t, inds, expand_tbeginends(t, first(I)).I)
end

@noinline function Base.to_indices(
    ::Tensor{<:Any,<:Any,ipos},
    _,
    ::Tuple{TIndex{<:Any,_ipos}},
) where {ipos,_ipos}
    throw(ArgumentError("Index positions $_ipos don't match indices of Tensor $ipos"))
end

# needed for views
Base.to_indices(t::Tensor, I::Tuple{TIndex}) = to_indices(t, axes(t), I)

###
### TI[1; 2 3; 4]
###

struct _TI end
const TI = _TI()

### `@generated` typed_hvcat

function static_index_pos_from_alternations(alt, ipos=true)
    i = 0
    x = UInt(0)
    while !isempty(alt)
        if first(alt) > 0
            x |= UInt(ipos) << i
            i += 1
            alt = (first(alt) - 1, Base.tail(alt)...)
        else
            ipos = !ipos
            alt = Base.tail(alt)
        end
    end
    :(IndexPos{$i}($x))
end

expanded_dims(::Type{<:Union{Integer,AbstractUnitRange{<:Integer},Colon}}) = 1
expanded_dims(::Type{<:Union{CartesianIndex{N},CartesianIndices{N}}}) where {N} = N
expanded_dims(::Type{BitArray{N}}) where {N} = N
expanded_dims(::Type{AbstractArray{T,N}}) where {T,N} = expanded_dims(T)
expanded_dims(::Type) = -1

function static_expand_indices(alt::Tuple, inds::Tuple{Vararg{Type}})
    N = 0
    new_alt = (0,)
    alt_ptr = 1
    for T_i in inds
        n = expanded_dims(T_i)::Int
        n == -1 && return nothing
        N += n
        new_alt = (Base.front(new_alt)..., last(new_alt) + n)
        if alt_ptr == alt[length(new_alt)]
            new_alt = (new_alt..., 0)
            alt_ptr = 1
        else
            alt_ptr += 1
        end
    end
    new_alt, N
end

function Base.typed_hvcat(
    ::_TI,
    ::Val{alt},
    inds...;
    start_ipos::Val{_start_ipos} = Val(true)
) where {alt,_start_ipos}
    if @generated
        e = static_expand_indices(alt, inds)
        e === nothing && return :(
            _typed_hvcat_fallback(TI, $alt, $(inds...), start_ipos=_start_ipos)
        )
        alt, N = e
        ipos = static_index_pos_from_alternations(alt, _start_ipos)
        :(Base.@_inline_meta; TIndex{$N,$ipos}(inds))
    else
        _typed_hvcat_fallback(TI, alt, inds..., start_ipos=_start_ipos)
    end
end
@inline function Base.typed_hvcat(::_TI, ::Val{alt}, ::typeof(\), inds...) where {alt}
    Base.typed_hvcat(TI, Val{Base.tail(alt)}(), inds...; start_ipos=Val(false))
end
@inline function Base.typed_hvcat(::_TI, alt::T, inds...) where {T<:Tuple}
    Base.typed_hvcat(TI, Val{alt}(), inds...)
end

### non-`@generated` typed_hvcat

function index_pos_from_alternations_fallback(alt::Tuple, ipos=true)
    i = 0
    x = UInt(0)
    while !isempty(alt)
        if first(alt) > 0
            x |= UInt(ipos) << i
            i += 1
            alt = (first(alt) - 1, Base.tail(alt)...)
        else
            ipos = !ipos
            alt = Base.tail(alt)
        end
    end
    IndexPos{i}(x)
end

function expand_indices_fallback(alt, inds)
    new_inds = ()
    new_alt = (0,)
    alt_ptr = 1
    for i in inds
        t = to_indices(AbstractArray, (), (i,))
        new_inds = (new_inds..., t...)
        new_alt = (Base.front(new_alt)..., last(new_alt) + length(t))
        if alt_ptr == alt[length(new_alt)]
            new_alt = (new_alt..., 0)
            alt_ptr = 1
        else
            alt_ptr += 1
        end
    end
    return new_alt, new_inds
end

function _typed_hvcat_fallback(::_TI, alt::Tuple, inds...; start_ipos=true) where {N}
    alt, _inds = expand_indices_fallback(alt, inds)
    ipos = index_pos_from_alternations_fallback(alt, start_ipos)
    TIndex{length(_inds),ipos}(inds)
end

for T in [:_TI, :Tensor]
    @eval begin
        @inline function Base.typed_vcat(x::$T, i::Vararg{Any,N}) where {N}
            Base.typed_hvcat(x, Val{ntuple(one, N)}(), i...)
        end
        @inline function Base.typed_hcat(x::$T, i::Vararg{Any,N}) where {N}
            Base.typed_hvcat(x, Val{(N,)}(), i...)
        end
    end
end
# for TI[] and TI[1]
@inline Base.getindex(::_TI) = Base.typed_hvcat(TI, Val{()}())
@inline Base.getindex(::_TI, i) = Base.typed_hvcat(TI, Val{(1,)}(), i)

###
### (::Tensor)[1; 2 3; 4]
###

# a[1; 2 3; 4] gets transformed into a[TI[1; 2 3; 4]]
@inline function Base.typed_hvcat(t::Tensor, alt::T, i...) where {T<:Tuple}
    Base.typed_hvcat(t, Val{alt}(), i...)
end
@inline function Base.typed_hvcat(
    t::Tensor{T,N,ipos},
    alt::Val{_alt},
    i...,
) where {T,N,ipos,_alt}
    t[Base.typed_hvcat(TI, alt, i...)]
end
@inline function Base.typed_hvcat(
    t::Tensor{T,N,ipos},
    alt::Val{_alt},
    ::typeof(\),
    i...,
) where {T,N,ipos,_alt}
    t[Base.typed_hvcat(TI, alt, \, i...)]
end


text:
mutable struct SimpleZone{T<:Integer} <: AbstractZone{T}
    id::T
    issource::Bool
    issink::Bool
    thruallowed::Bool
    props::Union{Nothing, Dict}
end

function SimpleZone{T}(d::Dict; idkey, issrckey, issnkkey, thrukey) where T <: Integer
    (id = d[idkey]; delete!(d, idkey))
    (issource = d[issrckey]; delete!(d, issrckey))
    (issink = d[issnkkey]; delete!(d, issnkkey))
    (thruallowed = d[thrukey]; delete!(d, thrukey))
    SimpleZone{T}(id, issource, issink, thruallowed, d)
end

SimpleZone{T}(dfrow::DataFrameRow, args...; kwargs...) where T <: Integer = SimpleZone{T}(Dict(pairs(dfrow)); idkey=kwargs[:idkey], issrckey=kwargs[:issrckey], issnkkey=kwargs[:issnkkey], thrukey=kwargs[:thrukey])

id(z::SimpleZone{T}) where T = z.id
issource(z::SimpleZone{T}) where T = z.issource
issink(z::SimpleZone{T}) where T = z.issink
throughflowallowed(z::SimpleZone{T}) where T = z.thruallowed
props(z::SimpleZone{T}) where T = z.props


text:
module OpenEXR

export load_exr, save_exr

using FileIO

module C
using OpenEXR_jll
using Colors
const ImfHalf = Float16
const ImfRgba = RGBA{ImfHalf}
include("OpenEXR_common.jl")
include("OpenEXR_api.jl")
const IMF_WRITE_RGBA = IMF_WRITE_RGB + IMF_WRITE_A
end  # module C

function check(ret)
    ret == 0 && error(unsafe_string(C.ImfErrorMessage()))
end

"""
    load_exr(filename)

If `filename` is a file in OpenEXR format, return a tuple containing the corresponding
RGBA{Float16} image along with in integer encoding the populated channels.
"""
function load_exr(filename)
    infile = C.ImfOpenInputFile(filename)  # open the file
    check(infile)
    try
        # get the header
        hdr = C.ImfInputHeader(infile)

        # read its data window
        xmin = Ref{Cint}()
        ymin = Ref{Cint}()
        xmax = Ref{Cint}()
        ymax = Ref{Cint}()
        C.ImfHeaderDataWindow(hdr, xmin, ymin, xmax, ymax)

        # compute the window size
        width = xmax[]-xmin[]+1
        height = ymax[]-ymin[]+1

        # allocate space for the result and get its strides
        data = Array{C.ImfRgba, 2}(undef, height, width)
        (xstride, ystride) = strides(data)

        # get the pointer to the data, shifting it according to the expected window
        dataptr = Base.unsafe_convert(Ptr{C.ImfRgba}, data) - xmin[] * xstride - ymin[] * ystride

        # copy the data
        check(C.ImfInputSetFrameBuffer(infile, dataptr, ystride, xstride))
        check(C.ImfInputReadPixels(infile, ymin[], ymax[]))

        # return the loaded raster along with the channels
        return (data, C.ImfInputChannels(infile))
    finally
        check(C.ImfCloseInputFile(infile))
    end
end

"""
    save_exr(filename, image[, channels])

Save the channels of `image` indicated by `channels` (by default, all are saved) into a
file in OpenEXR format named `filename`.
"""
function save_exr(filename, image::AbstractArray{C.ImfRgba, 2}, channels = C.IMF_WRITE_RGBA)
    # get the size of the data
    (height, width) = size(image)

    # create a new header
    hdr = C.ImfNewHeader()
    check(hdr)
    try
        # set the correct window sizes
        C.ImfHeaderSetDataWindow(hdr, 0, 0, width-1, height-1)
        C.ImfHeaderSetDisplayWindow(hdr, 0, 0, width-1, height-1)

        # open the output file
        outfile = C.ImfOpenOutputFile(filename, hdr, channels)
        check(outfile)
        try
            # get the strides and a pointer to the raster
            (xstride, ystride) = strides(image)
            dataptr = Base.unsafe_convert(Ptr{C.ImfRgba}, image)

            # copy the data
            check(C.ImfOutputSetFrameBuffer(outfile, dataptr, ystride, xstride))
            check(C.ImfOutputWritePixels(outfile, height))
        finally
            check(C.ImfCloseOutputFile(outfile))
        end
    finally
        C.ImfDeleteHeader(hdr)
    end
    nothing
end

function save_exr(f, image::AbstractArray{T, 2}, channels = C.IMF_WRITE_RGBA) where T
    save_exr(f, (c->convert(C.ImfRgba, c)).(image), channels)
end

# FileIO interface
load(f::File{DataFormat{:EXR}}) = load_exr(f.filename)[1]
save(f::File{DataFormat{:EXR}}, args...) = save_exr(f.filename, args...)

end  # module OpenEXR


text:
@inline get_xlsxfile(wb::Workbook) :: XLSXFile = wb.package
@inline get_xlsxfile(ws::Worksheet) :: XLSXFile = ws.package
@inline get_workbook(ws::Worksheet) :: Workbook = get_xlsxfile(ws).workbook
@inline get_workbook(xl::XLSXFile) :: Workbook = xl.workbook

function check_for_xlsx_file_format(filepath::AbstractString)
    @assert isfile(filepath) "File $filepath not found."
    io = open(filepath, "r")
    header = Base.read(io, 4)
    close(io)

    if header == [ 0x50, 0x4b, 0x03, 0x04 ] # valid Zip file header
        return
    elseif header == [ 0xd0, 0xcf, 0x11, 0xe0 ] # old XLS file
        error("$filepath looks like an old XLS file (not XLSX). This package does not support XLS file format.")
    else
        error("$filepath is not a valid XLSX file.")
    end
end

"""
    readxlsx(filepath) :: XLSXFile

Main function for reading an Excel file.
This function will read the whole Excel file into memory
and return a closed XLSXFile.

Consider using `openxlsx` for lazy loading of Excel file contents.
"""
@inline readxlsx(filepath::AbstractString) :: XLSXFile = open_or_read_xlsx(filepath, true, true, false)

"""
    openxlsx(filepath; [enable_cache]) :: XLSXFile

Open a XLSX file for reading. The user must close this file after using it with `close(xf)`.
XML data will be fetched from disk as needed.

If `enable_cache=true`, all read worksheet cells will be cached.
If you read a worksheet cell twice it will use the cached value instead of reading from disk
in the second time.

If `enable_cache=false`, worksheet cells will always be read from disk.
This is useful when you want to read a spreadsheet that doesn't fit into memory.

The default value is `enable_cache=true`.

The following example shows how you would read worksheet cells, one row at a time,
where `filename.xlsx` is a spreadsheet that doesn't fit into memory.

```julia
julia> f = XLSX.openxlsx("filename.xlsx", enable_cache=false)

julia for r in XLSX.eachrow(f["sheetname"])
        # do something with r
      end
```

See also `readxlsx` method.
"""
@inline openxlsx(filepath::AbstractString; enable_cache::Bool=true) :: XLSXFile = open_or_read_xlsx(filepath, false, enable_cache, false)

function open_or_read_xlsx(filepath::AbstractString, read_files::Bool, enable_cache::Bool, read_as_template::Bool) :: XLSXFile
    # sanity check
    if read_as_template
        @assert read_files && enable_cache
    end

    xf = XLSXFile(filepath, enable_cache, read_as_template)

    try
        for f in xf.io.files

            # ignore xl/calcChain.xml in any case (#31)
            if f.name == "xl/calcChain.xml"
                continue
            end

            if endswith(f.name, ".xml") || endswith(f.name, ".rels")
                # XML file
                internal_xml_file_add!(xf, f.name)
                if read_files

                    # ignore worksheet files because they'll be read thru streaming
                    # If reading as template, it will be loaded in two places: here and WorksheetCache.
                    if !read_as_template && startswith(f.name, "xl/worksheets") && endswith(f.name, ".xml")
                        continue
                    end

                    internal_xml_file_read(xf, f.name)
                end
            elseif read_as_template

                # Binary file
                # we only read binary files to save the Excel file later
                xf.binary_data[f.name] = ZipFile.read(f, UInt8, f.uncompressedsize)
            end
        end

        check_minimum_requirements(xf)
        parse_relationships!(xf)
        parse_workbook!(xf)

        # read data from Worksheet streams
        if read_files
            for sheet_name in sheetnames(xf)
                sheet = getsheet(xf, sheet_name)

                # to read sheet content, we just need to iterate a SheetRowIterator and the data will be stored in cache
                for r in eachrow(sheet)
                    nothing
                end
            end
        end

        if read_as_template
            wb = get_workbook(xf)
            if has_sst(wb)
                sst_load!(wb)
            end
        end

    finally
        if read_files
            close(xf)
        end
    end

    return xf
end

function get_default_namespace(r::EzXML.Node) :: String
    for (prefix, ns) in EzXML.namespaces(r)
        if prefix == ""
            return ns
        end
    end

    error("No default namespace found.")
end

# See section 12.2 - Package Structure
function check_minimum_requirements(xf::XLSXFile)
    mandatory_files = ["_rels/.rels",
                       "xl/workbook.xml",
                       "[Content_Types].xml",
                       "xl/_rels/workbook.xml.rels"
                       ]

    for f in mandatory_files
        @assert in(f, filenames(xf)) "Malformed XLSX File. Couldn't find file $f in the package."
    end

    nothing
end

"""
Parses package level relationships defined in `_rels/.rels`.
Prases workbook level relationships defined in `xl/_rels/workbook.xml.rels`.
"""
function parse_relationships!(xf::XLSXFile)

    # package level relationships
    xroot = get_package_relationship_root(xf)
    for el in EzXML.eachelement(xroot)
        push!(xf.relationships, Relationship(el))
    end
    @assert !isempty(xf.relationships) "Relationships not found in _rels/.rels!"

    # workbook level relationships
    wb = get_workbook(xf)
    xroot = get_workbook_relationship_root(xf)
    for el in EzXML.eachelement(xroot)
        push!(wb.relationships, Relationship(el))
    end
    @assert !isempty(wb.relationships) "Relationships not found in xl/_rels/workbook.xml.rels"

    nothing
end

"""
  parse_workbook!(xf::XLSXFile)

Updates xf.workbook from xf.data[\"xl/workbook.xml\"]
"""
function parse_workbook!(xf::XLSXFile)
    xroot = xmlroot(xf, "xl/workbook.xml")
    @assert EzXML.nodename(xroot) == "workbook" "Malformed xl/workbook.xml. Root node name should be 'workbook'. Got '$(EzXML.nodename(xroot))'."

    # workbook to be parsed
    workbook = get_workbook(xf)

    # workbookPr
    local foundworkbookPr::Bool = false
    for node in EzXML.eachelement(xroot)

        if EzXML.nodename(node) == "workbookPr"
            foundworkbookPr = true

            # read date1904 attribute
            if haskey(node, "date1904")
                attribute_value_date1904 = node["date1904"]

                if attribute_value_date1904 == "1" || attribute_value_date1904 == "true"
                    workbook.date1904 = true
                elseif attribute_value_date1904 == "0" || attribute_value_date1904 == "false"
                    workbook.date1904 = false
                else
                    error("Could not parse xl/workbook -> workbookPr -> date1904 = $(attribute_value_date1904).")
                end
            else
                # does not have attribute => is not date1904
                workbook.date1904 = false
            end

            break
        end
    end
    @assert foundworkbookPr "Malformed: couldn't find workbookPr node element in 'xl/workbook.xml'."

    # sheets
    sheets = Vector{Worksheet}()
    for node in EzXML.eachelement(xroot)
        if EzXML.nodename(node) == "sheets"

            for sheet_node in EzXML.eachelement(node)
                @assert EzXML.nodename(sheet_node) == "sheet" "Unsupported node $(EzXML.nodename(sheet_node)) in 'xl/workbook.xml'."
                worksheet = Worksheet(xf, sheet_node)
                push!(sheets, worksheet)
            end

            break
        end
    end
    workbook.sheets = sheets

    # named ranges
    for node in EzXML.eachelement(xroot)
        if EzXML.nodename(node) == "definedNames"
            for defined_name_node in EzXML.eachelement(node)
                @assert EzXML.nodename(defined_name_node) == "definedName"
                defined_value_string = EzXML.nodecontent(defined_name_node)
                name = defined_name_node["name"]

                local defined_value::DefinedNameValueTypes

                if is_valid_fixed_sheet_cellname(defined_value_string) || is_valid_sheet_cellname(defined_value_string)
                    defined_value = SheetCellRef(defined_value_string)
                elseif is_valid_fixed_sheet_cellrange(defined_value_string) || is_valid_sheet_cellrange(defined_value_string)
                    defined_value = SheetCellRange(defined_value_string)
                elseif ismatch(r"^\".*\"$", defined_value_string) # is enclosed by quotes
                    defined_value = defined_value_string[2:end-1] # remove enclosing quotes
                    if isempty(defined_value)
                        defined_value = Missings.missing
                    end
                elseif tryparse(Int, defined_value_string)
                    defined_value = parse(Int, defined_value_string)
                elseif tryparse(Float64, defined_value_string)
                    defined_value = parse(Float64, defined_value_string)
                elseif isempty(defined_value_string)
                    defined_value = Missings.missing
                else

                    # Couldn't parse definedName. Will silently ignore it, since this is not a critical feature.
                    continue

                    # debug
                    #error("Could not parse value $(defined_value_string) for definedName $name.")
                end

                if haskey(defined_name_node, "localSheetId")
                    # is a Worksheet level name

                    # localSheetId is the 0-based index of the Worksheet in the order
                    # that it is displayed on screen.
                    # Which is the order of the elements under <sheets> element in workbook.xml .
                    localSheetId = parse(Int, defined_name_node["localSheetId"]) + 1
                    sheetId = workbook.sheets[localSheetId].sheetId
                    workbook.worksheet_names[(sheetId, name)] = defined_value
                else
                    # is a Workbook level name
                    workbook.workbook_names[name] = defined_value
                end
            end

            break
        end
    end

    nothing
end

function tryparse(t::Type, s::String)
    try
        parse(t, s)
        return true
    catch
        return false
    end
end

# Lazy loading of XML files

"""
Lists internal files from the XLSX package.
"""
@inline filenames(xl::XLSXFile) = keys(xl.files)

"""
Returns true if the file data was read into xl.data.
"""
@inline internal_xml_file_isread(xl::XLSXFile, filename::String) :: Bool = xl.files[filename]
@inline internal_xml_file_exists(xl::XLSXFile, filename::String) :: Bool = haskey(xl.files, filename)

function internal_xml_file_add!(xl::XLSXFile, filename::String)
    @assert endswith(filename, ".xml") || endswith(filename, ".rels")
    xl.files[filename] = false
    nothing
end

function internal_xml_file_read(xf::XLSXFile, filename::String) :: EzXML.Document
    @assert internal_xml_file_exists(xf, filename) "Couldn't find $filename in $(xf.filepath)."

    if !internal_xml_file_isread(xf, filename)
        @assert isopen(xf) "Can't read from a closed XLSXFile."
        file_not_found = true
        for f in xf.io.files
            if f.name == filename
                xf.files[filename] = true # set file as read
                xf.data[filename] = EzXML.readxml(f)
                file_not_found = false
                break
            end
        end

        if file_not_found
            # shouldn't happen
            error("$filename not found in XLSX package.")
        end
    end

    return xf.data[filename]
end

function Base.close(xl::XLSXFile)
    xl.io_is_open = false
    close(xl.io)

    # close all internal file streams from worksheet caches
    for sheet in xl.workbook.sheets
        if !isnull(sheet.cache)
            cache = get(sheet.cache)
            close(cache.stream_state)
        end
    end
end

Base.isopen(xl::XLSXFile) = xl.io_is_open

"""
    xmldocument(xl::XLSXFile, filename::String) :: EzXML.Document

Utility method to find the XMLDocument associated with a given package filename.
Returns xl.data[filename] if it exists. Throws an error if it doesn't.
"""
@inline xmldocument(xl::XLSXFile, filename::String) :: EzXML.Document = internal_xml_file_read(xl, filename)

"""
    xmlroot(xl::XLSXFile, filename::String) :: EzXML.Node

Utility method to return the root element of a given XMLDocument from the package.
Returns EzXML.root(xl.data[filename]) if it exists.
"""
@inline xmlroot(xl::XLSXFile, filename::String) :: EzXML.Node = EzXML.root(xmldocument(xl, filename))


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reverse_diff)
#nb Turing.turnprogress(false);

d = DataFrame(CSV.read(joinpath(@__DIR__, "..", "..", "data", "reedfrogs.csv"),
  delim=';'));

size(d) |> display # Should be 48x5

# Set number of tanks

d[:tank] = 1:size(d,1);

# Thanks to Kai Xu!

@model m12_2(density, tank, surv) = begin

    # Separate priors on α and σ for each tank
    σ ~ Truncated(Cauchy(0, 1), 0, Inf)
    α ~ Normal(0, 1)

    # Number of unique tanks in the data set
    N_tank = length(tank)

    # vector for the priors for each tank
    a_tank = Vector{Real}(undef, N_tank)

    # For each tank [1,..,48] set a prior. Note the [] around Normal().
    a_tank ~ [Normal(α, σ)]

    # Observation
    logitp = [a_tank[tank[i]] for i = 1:N_tank]
    surv ~ VecBinomialLogit(density, logitp)

end

# Sample

chns = sample(m12_2(Vector{Int64}(d[:density]), Vector{Int64}(d[:tank]),
    Vector{Int64}(d[:surv])), Turing.NUTS(0.8), 1000);

# CmdStan result

m122rethinking = "
                mean   sd  5.5% 94.5% n_eff Rhat
a               1.30 0.25  0.90  1.70 11662    1
sigma       1.62 0.22  1.30  1.99  6556    1
a_tank[1]   2.12 0.88  0.84  3.60 16091    1
a_tank[2]   3.05 1.10  1.52  4.92 10962    1
a_tank[3]   1.00 0.66 -0.02  2.10 18175    1
a_tank[4]   3.05 1.11  1.47  4.96 10181    1
a_tank[5]   2.13 0.87  0.85  3.62 13720    1
a_tank[6]   2.12 0.86  0.86  3.59 11628    1
a_tank[7]   3.07 1.13  1.47  5.03 10315    1
a_tank[8]   2.13 0.87  0.86  3.60 13754    1
a_tank[9]  -0.18 0.60 -1.14  0.76 18218    1
a_tank[10]  2.11 0.86  0.83  3.58 15121    1
a_tank[11]  1.00 0.67 -0.04  2.09 17390    1
a_tank[12]  0.58 0.62 -0.41  1.60 17209    1
a_tank[13]  0.99 0.66 -0.04  2.09 15225    1
a_tank[14]  0.19 0.62 -0.79  1.20 18293    1
a_tank[15]  2.13 0.89  0.83  3.63 12445    1
a_tank[16]  2.11 0.87  0.87  3.61 12385    1
a_tank[17]  2.89 0.80  1.76  4.29 12583    1
a_tank[18]  2.38 0.66  1.43  3.49 14437    1
a_tank[19]  2.00 0.58  1.12  2.99 13959    1
a_tank[20]  3.67 1.03  2.20  5.44 10629    1
a_tank[21]  2.38 0.65  1.42  3.47 15309    1
a_tank[22]  2.39 0.66  1.42  3.49 13614    1
a_tank[23]  2.40 0.67  1.41  3.53 11868    1
a_tank[24]  1.69 0.52  0.90  2.55 18468    1
a_tank[25] -1.00 0.45 -1.74 -0.30 18153    1
a_tank[26]  0.16 0.40 -0.47  0.81 21895    1
a_tank[27] -1.44 0.50 -2.28 -0.69 16718    1
a_tank[28] -0.47 0.41 -1.15  0.17 20160    1
a_tank[29]  0.15 0.40 -0.48  0.80 19401    1
a_tank[30]  1.44 0.49  0.70  2.24 15407    1
a_tank[31] -0.64 0.42 -1.33 -0.01 15356    1
a_tank[32] -0.31 0.40 -0.95  0.32 19130    1
a_tank[33]  3.18 0.78  2.06  4.55 10894    1
a_tank[34]  2.70 0.66  1.75  3.84 13573    1
a_tank[35]  2.69 0.64  1.74  3.78 13876    1
a_tank[36]  2.06 0.53  1.26  2.92 16329    1
a_tank[37]  2.06 0.51  1.29  2.91 14672    1
a_tank[38]  3.88 0.97  2.52  5.57  9349    1
a_tank[39]  2.70 0.64  1.77  3.78 13444    1
a_tank[40]  2.34 0.56  1.49  3.31 14966    1
a_tank[41] -1.82 0.48 -2.61 -1.10 14214    1
a_tank[42] -0.58 0.36 -1.16 -0.02 17203    1
a_tank[43] -0.46 0.35 -1.02  0.08 17762    1
a_tank[44] -0.34 0.34 -0.90  0.20 16740    1
a_tank[45]  0.58 0.35  0.02  1.14 18946    1
a_tank[46] -0.57 0.34 -1.13 -0.03 19761    1
a_tank[47]  2.05 0.51  1.30  2.90 15122    1
a_tank[48]  0.00 0.33 -0.53  0.53 18236    1
";

# Draw summary

describe(chns)

# End of `12/m12.2t.jl`


text:
# Test if tf listner works correctly
using RobotOS
@rosimport std_msgs.msg: Bool
rostypegen(@__MODULE__)
import .std_msgs.msg: BoolMsg

using ROSTransforms.TF
using Test

# Initialize node
init_node("listner", anonymous=true) # delete

# Setup the transform listener
tl = TransformListener()
# Wait and collect the transform
waitForTransform(tl, "parent_link", "child_link", Time(), Duration(1.0))
transform = lookupTransform(tl, "parent_link", "child_link", Time())

@testset "TransformListener" begin
    tf_correct = Transform([1, 2, 3.], [0, 0, 0, 1.])
    @test transform == tf_correct
end

# Publsih success
pub = Publisher{BoolMsg}("/success", queue_size=10)
for i ∈ 1:100
    publish(pub, BoolMsg(true))
    rossleep(.1)
end


text:
################################################################
################################################################
# Entering_Commands
################################################################
################################################################
# # Copy of Julia Functions/Types

################################################################
# ansM can't be defined. use ans instead.
################################################################
"""
    clcM()

Clear Command Window. Works for any IDE or console.
"""
function clcM()

    if !isdefined(Main, :Juno)
        println("\33[2J")
    end

end


text:
module MaybesTest

using Funs, Maybes
using Base.Test

function test_basic()
    m0 = Maybe{Int}()
    m1 = Maybe{Int}(42)
    m2 = just(42)

    @test !isjust(m0)
    @test isnothing(m0)
    @test eltype(m0) == Int
    @test m0 == Maybe{Int}()
    @test m0 != Maybe{Int}(43)
    @test m0 != m1
    @test length(m0) == 0
    ms = [m for m in m0]
    @test isempty(ms)

    @test isjust(m1)
    @test !isnothing(m1)
    @test eltype(m1) == Int
    @test m1 == Maybe{Int}(42)
    @test m1 != Maybe{Int}(43)
    @test m1 != m0
    @test m1 == m2
    @test length(m1) == 1
    ms = [m for m in m1]
    @test !isempty(ms)
    @test ms == [42]

    @test typeof(m2) == Maybe{Int}

    @test maybe(0, x->x+1, m0) == 0
    @test maybe(0, x->x+1, m1) == 43
    @test fromjust(m1) == 42
    @test frommaybe(43, m0) == 43
    @test frommaybe(43, m1) == 42

    @test arrayToMaybe(Int[]) == Maybe{Int}()
    @test arrayToMaybe([42]) == just(42)
    @test arrayToMaybe([42, 43]) == just(42)
    @test maybeToArray(m0) == Int[]
    @test maybeToArray(m1) == Int[42]

    @test catMaybes([m0, m1, m2]) == [42, 42]
    @test mapMaybe(R=Int, x->x+1, [m0, m1, m2]) == [43, 43]

    @test string(m0) == "(?)"
    @test string(m1) == "(?42)"
end

function test_foldable()
    m0 = Maybe{Int}()
    m1 = Maybe{Int}(2)
    m2 = Maybe{Int}(3)
    @test freduce(+, 0, m0) == 0
    @test freduce(+, 0, m1) == 2
    @test freduce(+, 0, m1, m2) == 5
end

function test_functor()
    m0 = Maybe{Int}()
    m1 = Maybe{Int}(2)
    m2 = Maybe{Int}(3)
    @test fmap(R=Int, x->2*x, m0) == Maybe{Int}()
    @test fmap(R=Int, x->2*x, m1) == Maybe{Int}(4)
    @test fmap(R=Int, (x,y)->2*x+y, m1, m2) == Maybe{Int}(7)
end

function test_monad()
    u = munit(Maybe{Int}, 42)
    @test u == just(42)

    uu = munit(Maybe{Maybe{Int}}, u)
    @test uu == just(just(42))

    j = mjoin(uu)
    @test j == just(42)

    b = mbind(x::Int->just(x+1.0), just(1), R=Maybe{Float64})
    @test b == just(2.0)

    b = mbind(Fun{Maybe{Float64}}(x->just(x+1.0)), just(1))
    @test b == just(2.0)

    z = mzero(Maybe{Int})
    @test isempty(z)
    @test z == Maybe{Int}()

    p1 = mplus(Set([1,2]))
    p2 = mplus(Set([1,2]), Set([3,4]))
    p3 = mplus(Set([1,2]), Set([3,4]), Set([5,6]))
    @test p1 == Set([1,2])
    @test p2 == Set([1,2,3,4])
    @test p3 == Set([1,2,3,4,5,6])
end

function main()
    test_basic()
    test_foldable()
    test_functor()
    test_monad()
end

end


text:
function uniquearrayrefs(ls::LoopSet)
    uniquerefs = ArrayReferenceMeta[]
    # for arrayref ∈ ls.refs_aliasing_syms
    for op ∈ operations(ls)
        arrayref = op.ref
        arrayref === NOTAREFERENCE && continue
        any(ref -> sameref(arrayref, ref), uniquerefs) || push!(uniquerefs, arrayref)
    end
    uniquerefs
end

otherindexunrolled(loopsym::Symbol, ind::Symbol, loopdeps::Vector{Symbol}) = (loopsym !== ind) && (loopsym ∈ loopdeps)
function otherindexunrolled(ls::LoopSet, ind::Symbol, ref::ArrayReferenceMeta)
    us = ls.unrollspecification[]
    u₁sym = names(ls)[us.u₁loopnum]
    u₂sym = us.u₂loopnum > 0 ? names(ls)[us.u₂loopnum] : Symbol("##undefined##")
    otherindexunrolled(u₁sym, ind, loopdependencies(ref)) || otherindexunrolled(u₂sym, ind, loopdependencies(ref))
end
multiple_with_name(n::Symbol, v::Vector{ArrayReferenceMeta}) = sum(ref -> n === vptr(ref), v) > 1
# TODO: DRY between indices_calculated_by_pointer_offsets and use_loop_induct_var
function indices_calculated_by_pointer_offsets(ls::LoopSet, ar::ArrayReferenceMeta)
    indices = getindices(ar)
    ls.isbroadcast[] && return fill(false, length(indices))
    looporder = names(ls)
    offset = isdiscontiguous(ar)
    gespinds = Expr(:tuple)
    out = Vector{Bool}(undef, length(indices))
    li = ar.loopedindex
    # @show ls.vector_width[]
    for i ∈ eachindex(li)
        ii = i + offset
        ind = indices[ii]
        # j = findfirst(isequal(ind), view(indices, offset+1:ii-1))
        # if !isnothing(j)
        #     out[i] = out[j - offset]
        #     continue
        # end
        if (!li[i]) || multiple_with_name(vptr(ar), ls.lssm[].uniquearrayrefs) ||
            (iszero(ls.vector_width[]) && isstaticloop(getloop(ls, ind)))# ||
            # ((ls.align_loops[] > 0) && (first(names(ls)) == ind))
            out[i] = false
        elseif (isone(ii) && (first(looporder) === ind))
            out[i] = otherindexunrolled(ls, ind, ar)
        else
            out[i] = true
        end
    end
    out
end

"""
Returns a vector of length equal to the number of indices.
A value > 0 indicates which loop number that index corresponds to when incrementing the pointer.
A value < 0 indicates that abs(value) is the corresponding loop, and a `loopvalue` will be used.
"""
function use_loop_induct_var!(ls::LoopSet, q::Expr, ar::ArrayReferenceMeta, allarrayrefs::Vector{ArrayReferenceMeta})
    us = ls.unrollspecification[]
    li = ar.loopedindex
    looporder = reversenames(ls)
    uliv = Vector{Int}(undef, length(li))
    indices = getindices(ar)
    offset = first(indices) === Symbol("##DISCONTIGUOUSSUBARRAY##")
    if length(indices) != offset + length(li)
        println(ar)
        throw("Length of indices and length of offset do not match!")
    end
    isbroadcast = ls.isbroadcast[]
    gespinds = Expr(:tuple)
    offsetprecalc_descript = Expr(:tuple)
    use_offsetprecalc = false
    for i ∈ eachindex(li)
        ii = i + offset
        ind = indices[ii]
        # j = findfirst(isequal(ind), view(indices, offset+1:ii-1))
        # if !isnothing(j)
        #     j -= offset
        #     push!(gespinds.args, gespinds.args[j])
        #     uliv[i] = uliv[j]
        # else
        if (!li[i])
            uliv[i] = 0
            push!(gespinds.args, Expr(:call, lv(:Zero)))
            push!(offsetprecalc_descript.args, 0)
        elseif isbroadcast ||
            ((isone(ii) && (last(looporder) === ind)) && !(otherindexunrolled(ls, ind, ar)) ||
             multiple_with_name(vptr(ar), allarrayrefs)) ||
             (iszero(ls.vector_width[]) && isstaticloop(getloop(ls, ind)))# ||
             # ((ls.align_loops[] > 0) && (first(names(ls)) == ind))

            # Not doing normal offset indexing
            uliv[i] = -findfirst(isequal(ind), looporder)::Int
            push!(gespinds.args, Expr(:call, lv(:Zero)))
            push!(offsetprecalc_descript.args, 0) # not doing offset indexing, so push 0
        else
            uliv[i] = findfirst(isequal(ind), looporder)::Int
            loop = getloop(ls, ind)
            if loop.startexact
                push!(gespinds.args, Expr(:call, Expr(:curly, lv(:Static), loop.starthint - 1)))
            else
                push!(gespinds.args, Expr(:call, lv(:staticm1), loop.startsym))
            end
            if ind === names(ls)[us.vectorizedloopnum]
                push!(offsetprecalc_descript.args, 0)
            elseif (ind === names(ls)[us.u₁loopnum]) & (us.u₁ > 3)
                use_offsetprecalc = true
                push!(offsetprecalc_descript.args, us.u₁)
            elseif (ind === names(ls)[us.u₂loopnum]) & (us.u₂ > 3)
                use_offsetprecalc = true
                push!(offsetprecalc_descript.args, us.u₂)
            else
                push!(offsetprecalc_descript.args, 0)
            end
        end
    end
    if use_offsetprecalc
        push!(q.args, Expr(:(=), vptr(ar), Expr(:call, lv(:offsetprecalc), Expr(:call, lv(:gesp), vptr(ar), gespinds), Expr(:call, Expr(:curly, :Val, offsetprecalc_descript)))))
    else
        push!(q.args, Expr(:(=), vptr(ar), Expr(:call, lv(:gesp), vptr(ar), gespinds)))
    end
    uliv
end

# Plan here is that we increment every unique array
function add_loop_start_stop_manager!(ls::LoopSet)
    q = Expr(:block)
    us = ls.unrollspecification[]
    # Presence of an explicit use of a loopinducation var means we should use that, so we look for one
    # TODO: replace first with only once you add Compat as a dep or drop support for older Julia versions
    loopinductvars = map(op -> first(loopdependencies(op)), filter(isloopvalue, operations(ls)))
    # Filtered ArrayReferenceMetas, we must increment each
    arrayrefs = uniquearrayrefs(ls)
    use_livs = map(ar -> use_loop_induct_var!(ls, q, ar, arrayrefs), arrayrefs)
    # loops, sorted from outer-most to inner-most
    looporder = reversenames(ls)
    # For each loop, we need to choose an induction variable
    nloops = length(looporder)
    # loopstarts = Vector{Vector{ArrayReferenceMeta}}(undef, nloops)
    loopstarts = fill(ArrayReferenceMeta[], nloops)
    terminators = Vector{Int}(undef, nloops) # zero is standard loop induct var
    # loopincrements = Vector{Vector{ArrayReferenceMeta}}(undef, nloops) # Not needed; copy of loopstarts
    # The purpose of the reminds thing here is to pick which of these to use for the terminator
    # We want to pick the one with the fewest outer loops with respect to this one, to minimize
    # the number of redefinitions of max-pointer used for the termination comparison.
    reached_indices = zeros(Int, length(arrayrefs))
    for (i,loopsym) ∈ enumerate(looporder) # iterates from outer to inner
        loopstartᵢ = ArrayReferenceMeta[]; arⱼ = 0; minrem = typemax(Int);
        ric = Tuple{Int,Int}[]
        for j ∈ eachindex(use_livs) # j is array ref number
            for (l,k) ∈ enumerate(use_livs[j])# l is index number, k is loop number
                if k == i
                    push!(loopstartᵢ, arrayrefs[j])
                    push!(ric, ((reached_indices[j] += 1), length(loopstartᵢ)))
                end
            end
        end
        loopstarts[i] = loopstartᵢ
        terminators[i] = if (loopsym ∈ loopinductvars) || (any(r -> any(isequal(-i), r), use_livs)) || iszero(length(loopstartᵢ))
            0
        else
            # @show i, loopsym loopdependencies.(operations(ls)) operations(ls)
            # @assert !iszero(length(loopstartᵢ))
            last(ric[argmin(first.(ric))]) # index corresponds to array ref's position in loopstart
        end
    end
    ls.lssm[] = LoopStartStopManager(
        reverse!(terminators), reverse!(loopstarts), arrayrefs
    )
    q
end
maxsym(ptr, sub) = Symbol(ptr, "##MAX##", sub, "##")
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool)::Expr
    pointermax(ls, ar, n, sub, isvectorized, getloop(ls, names(ls)[n]))
end
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, loop::Loop)::Expr
    pointermax(ls, ar, n, sub, isvectorized, looplengthexpr(loop, n))::Expr
end
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, stophint::Int)::Expr
    # @unpack u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂ = us
    loopsym = names(ls)[n]
    index = Expr(:tuple)
    for i ∈ getindicesonly(ar)
        if i === loopsym
            if iszero(sub)
                push!(index.args, stophint)
            elseif isvectorized
                if isone(sub)
                    push!(index.args, Expr(:call, lv(:valsub), stophint, VECTORWIDTHSYMBOL))
                else
                    push!(index.args, Expr(:call, lv(:vsub), stophint, Expr(:call, lv(:valmul), VECTORWIDTHSYMBOL, sub)))
                end
            else
                push!(index.args, stophint - sub)
            end
            ptr = vptr(ar)
            return Expr(:call, lv(:pointerforcomparison), ptr, index)
        else
            push!(index.args, Expr(:call, lv(:Zero)))
        end
    end
    @show ar, loopsym
end
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, stopsym)::Expr
    # @unpack u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂ = us
    loopsym = names(ls)[n]
    index = Expr(:tuple)
    for i ∈ getindicesonly(ar)
        if i === loopsym
            if iszero(sub)
                push!(index.args, stopsym)
            elseif isvectorized
                if isone(sub)
                    push!(index.args, Expr(:call, lv(:valsub), stopsym, VECTORWIDTHSYMBOL))
                else
                    push!(index.args, Expr(:call, lv(:vsub), stopsym, Expr(:call, lv(:valmul), VECTORWIDTHSYMBOL, sub)))
                end
            else
                push!(index.args, Expr(:call, lv(:vsub), stopsym, sub))
            end
            return Expr(:call, lv(:pointerforcomparison), vptr(ar), index)
        else
            push!(index.args, Expr(:call, lv(:Zero)))
        end
    end
    @show ar, loopsym
end

function defpointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool)::Expr
    Expr(:(=), maxsym(vptr(ar), sub), pointermax(ls, ar, n, sub, isvectorized))
end

function maxunroll(us::UnrollSpecification, n)
    @unpack u₁loopnum, u₂loopnum, u₁, u₂ = us
    if n == u₁loopnum# && u₁ > 1
        u₁
    elseif n == u₂loopnum# && u₂ > 1
        u₂
    else
        1
    end
end

function startloop(ls::LoopSet, us::UnrollSpecification, n::Int, submax = maxunroll(us, n))
    @unpack u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂ = us
    lssm = ls.lssm[]
    termind = lssm.terminators[n]
    ptrdefs = lssm.incrementedptrs[n]
    loopstart = Expr(:block)
    firstloop = n == num_loops(ls)
    for ar ∈ ptrdefs
        ptr = vptr(ar)
        push!(loopstart.args, Expr(:(=), ptr, ptr))
    end
    if iszero(termind)
        loopsym = names(ls)[n]
        push!(loopstart.args, startloop(getloop(ls, loopsym), loopsym))
    else
        isvectorized = n == vectorizedloopnum
        for sub ∈ 0:submax
            push!(loopstart.args, defpointermax(ls, ptrdefs[termind], n, sub, isvectorized))
        end
    end
    loopstart
end
function offset_ptr(ar::ArrayReferenceMeta, us::UnrollSpecification, loopsym::Symbol, n::Int, UF::Int, offsetinds::Vector{Bool})
    indices = getindices(ar)
    offset = first(indices) === DISCONTIGUOUS
    gespinds = Expr(:tuple)
    li = ar.loopedindex
    for i ∈ eachindex(li)
        ii = i + offset
        ind = indices[ii]
        if !offsetinds[i]
            push!(gespinds.args, Expr(:call, lv(:Zero)))
        elseif ind == loopsym
            incrementloopcounter!(gespinds, us, n, UF)
        else
            push!(gespinds.args, Expr(:call, lv(:Zero)))
        end
        ind == loopsym && break
    end
    Expr(:(=), vptr(ar), Expr(:call, lv(:gesp), vptr(ar), gespinds))
end
function incrementloopcounter(ls::LoopSet, us::UnrollSpecification, n::Int, UF::Int)
    @unpack u₁loopnum, u₂loopnum, vectorizedloopnum, u₁, u₂ = us
    lssm = ls.lssm[]
    ptrdefs = lssm.incrementedptrs[n]
    looporder = names(ls)
    loopsym = looporder[n]
    q = Expr(:block)
    termind = lssm.terminators[n]
    if iszero(termind) # increment liv
        push!(q.args, incrementloopcounter(us, n, loopsym, UF))
    end
    for (j,ar) ∈ enumerate(ptrdefs)
        offsetinds = indices_calculated_by_pointer_offsets(ls, ar)
        push!(q.args, offset_ptr(ar, us, loopsym, n, UF, offsetinds))
    end
    q
end
function terminatecondition(ls::LoopSet, us::UnrollSpecification, n::Int, inclmask::Bool, UF::Int)
    lssm = ls.lssm[]
    termind = lssm.terminators[n]
    if iszero(termind)
        loop = getloop(ls, names(ls)[n])
        return terminatecondition(loop, us, n, loop.itersymbol, inclmask, UF)
    end

    termar = lssm.incrementedptrs[n][termind]
    ptr = vptr(termar)
    # @show UF, isvectorized(us, n)
    if inclmask && isvectorized(us, n)
        Expr(:call, :<, callpointerforcomparison(ptr), maxsym(ptr, 0))
    else
        Expr(:call, :≤, callpointerforcomparison(ptr), maxsym(ptr, UF))
    end
end


text:
using Documenter, SciMLBase

makedocs(
    sitename="SciMLBase.jl",
    authors="Chris Rackauckas",
    modules=[SciMLBase],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://scimlbase.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Fundamentals" => Any[
            "fundamentals/Problems.md",
            "fundamentals/PDE.md",
            "fundamentals/SciMLFunctions.md",
            "fundamentals/Differentiation.md",
            "fundamentals/FAQ.md"
        ]
    ]
)

deploydocs(
   repo = "github.com/SciML/SciMLBase.jl.git";
   push_preview = true
)


text:
using Term
import Term: Renderable

pprint(x) = begin
    println(x)
    println(x.measure)
end

p1 = Panel()
p2 = Panel(; width = 24, height = 3)
p3 = Panel("test[red]aajjaja[/red]"^5; width = 12)

# pprint(p1 * p2)
# pprint(p1 / p2)

a = p1 / p2


text:
##

##
"""
    function skeleton()

This 'Scowleton' brought to you by one of
{appel,kube,rowe}@csua.berkeley.edu

# Example
```jldoctest
julia> cowsay("All T-Bone, no steak", cow=Cowsay.skeleton)
 ______________________
< All T-Bone, no steak >
 ----------------------
          \\      (__)
           \\     /oo|
            \\   (_"_)*+++++++++*
                   //I#\\\\\\\\\\\\\\\\I\\
                   I[I|I|||||I I `
                   I`I'///'' I I
                   I I       I I
                   ~ ~       ~ ~
                     Scowleton
```
"""
function skeleton(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
          $thoughts      (__)
           $thoughts     /$eyes|
            $thoughts   (_"_)*+++++++++*
                   //I#\\\\\\\\\\\\\\\\I\\
                   I[I|I|||||I I `
                   I`I'///'' I I
                   I I       I I
                   ~ ~       ~ ~
                     Scowleton
"""
return the_cow
end


text:
"""
    rf = downscale_spaceonly(r, f, weight=1.; fglob=false, fsmooth=false)

Downscale precipitation field `r` using metagaussian spectral field `f`. An optional weights array can be specified. Precipitation can be conserved globally (`fglob`) or using convolution (`fsmooth`).
"""
function downscale_spaceonly(r, f, weight=1.; fglob=false, fsmooth=false)

    (nas, nas) = size(r)
    (ns, ns) = size(f)

    # No negative rain
    ii = findall(r[:].<=0)
    r[ii] .= 0
    rs = r

    # Gaussianize
    # Alternative: rg=log(rs)
    rg = gaussianize(rs)

    # Prepare small scale random fields
    g = metagauss(f)

    # fix the amplitude of f
    (pstg, pstgt) = fft3d(g)
    pstg = pstg*((ns*ns)^2)
    (pstrg, pstrgt) = fft3d(rg)
    pstrg = pstrg*((nas*nas)^2)

    varg = pstg[div(nas, 2)]
    varrg = pstrg[div(nas, 2)]

    c = varrg/varg
    g = g*sqrt(c)

    gm = mergespec_spaceonly(rg, g, div(nas, 2))

    st = Statistics.std(gm)
    if st==0.
        st=1.
    end
    gm = gm/st
    fm = exp.(gm)
    fm = fm.*weight

    # We want the aggregated field to be the same as pa
    ii = findall(isnan.(r)); fm[ii] .= NaN

    if fglob
        ii = findall(.~isnan.(fm))
        ri = interpola(r, ns, 1)
        fm = fm*Statistics.mean(ri[ii])/Statistics.mean(fm[ii])
    elseif fsmooth
        ii = findall(isnan.(fm))
        fmi = interpola(agg(fm, nas, 1), ns, 1)
        fmi[ii] .= NaN
        fma = smoothconv(fmi, nas)
        ri = interpola(r, ns, 1)
        ri[ii] .= NaN
        raa = smoothconv(ri, nas)
        fm = raa./fma.*fm  #pa=p aggregated at L0 and T0
    else
        #fm=mergespec_spaceonly(fm,raa,div(nas,2))
        raa = agg(r, nas, 1)
        fma = agg(fm, nas, 1)
        ca=raa./fma #pa=p aggregated at L0 and T0
        cai = interpola(ca,ns,1)
        fm = cai.*fm
    end
    return fm
end


text:
module FixedSizeStrings

export FixedSizeString

import Base: iterate, lastindex, getindex, sizeof, length, ncodeunits, codeunit, isvalid, read, write

struct FixedSizeString{N} <: AbstractString
    data::NTuple{N,UInt8}
    FixedSizeString{N}(itr) where {N} = new(NTuple{N,UInt8}(itr))
end

FixedSizeString(s::AbstractString) = FixedSizeString{length(s)}(s)

function iterate(s::FixedSizeString{N}, i::Int = 1) where N
    i > N && return nothing
    return (Char(s.data[i]), i+1)
end

lastindex(s::FixedSizeString{N}) where {N} = N

getindex(s::FixedSizeString, i::Int) = Char(s.data[i])

sizeof(s::FixedSizeString) = sizeof(s.data)

length(s::FixedSizeString) = length(s.data)

ncodeunits(s::FixedSizeString) = length(s.data)

codeunit(::FixedSizeString) = UInt8
codeunit(s::FixedSizeString, i::Integer) = s.data[i]

isvalid(s::FixedSizeString, i::Int) = checkbounds(Bool, s, i)

function read(io::IO, T::Type{FixedSizeString{N}}) where N
    return read!(io, Ref{T}())[]::T
end

function write(io::IO, s::FixedSizeString{N}) where N
    return write(io, Ref(s))
end

end


text:
function iterative_probabilistic_improvement(tuning_run::Run,
                                             reference::RemoteChannel;
                                             threshold::AbstractFloat = 2.)
    cost_calls         = tuning_run.cost_evaluations
    iteration          = 1
    name               = "Iterative Probabilistic Improvement"
    stopping_criterion = @task tuning_run.stopping_criterion(tuning_run.duration)
    stop               = consume(stopping_criterion)

    while !stop
        iteration                 += 1
        result                     = probabilistic_improvement(tuning_run,
                                                               threshold = threshold)
        cost_calls                += result.cost_calls
        result.cost_calls          = cost_calls
        result.start               = tuning_run.starting_point
        result.technique           = name
        result.iterations          = iteration
        result.current_iteration   = iteration
        tuning_run.starting_point  = result.minimum
        tuning_run.starting_cost   = result.cost_minimum
        stop                       = consume(stopping_criterion)
        put!(reference, result)
    end
end


text:
const MAX_OR_MIN = Union{typeof(max),typeof(min)}

"""
    extreme_filter(f, A; r=1, [dims]) -> out
    extreme_filter(f, A, Ω) -> out

Filter the array `A` using select function `f(x, y)` for each Ω-neighborhood. The name
"extreme" comes from the fact that typical select function `f` choice is `min` and `max`.

For each pixel `p` in `A`, the select function `f` is applied to its Ω-neighborhood
iteratively in a `f(...(f(f(A[p], A[p+Ω[1]]), A[p+Ω[2]]), ...)` manner. For instance, in the
1-dimensional case, `out[p] = f(f(A[p], A[p-1]), A[p+1])` for each `p` is the default
behavior.

The Ω-neighborhood is defined by the `dims` or `Ω` argument. The `r` and `dims` keywords
specifies the box shape neighborhood `Ω` using [`strel_box`](@ref). The `Ω` is also known as
structuring element (SE), it can be either displacement offsets or bool array mask, please
refer to [`strel`](@ref) for more details.

# Examples

```jldoctest extreme_filter; setup=:(using ImageMorphology)
julia> M = [4 6 5 3 4; 8 6 9 4 8; 7 8 4 9 6; 6 2 2 1 7; 1 6 5 2 6]
5×5 Matrix{$Int}:
 4  6  5  3  4
 8  6  9  4  8
 7  8  4  9  6
 6  2  2  1  7
 1  6  5  2  6

julia> extreme_filter(max, M) # max-filter using 4 direct neighbors along both dimensions
5×5 Matrix{$Int}:
 8  9  9  9  8
 8  9  9  9  9
 8  9  9  9  9
 8  8  9  9  9
 6  6  6  7  7

julia> extreme_filter(max, M; dims=1) # max-filter along the first dimension (column)
5×5 Matrix{$Int}:
 8  6  9  4  8
 8  8  9  9  8
 8  8  9  9  8
 7  8  5  9  7
 6  6  5  2  7
```

`Ω` can be either an `AbstractArray{Bool}` mask array with `true` element indicating
connectivity, or a `AbstractArray{<:CartesianIndex}` array with each element indicating the
displacement offset to its center element.

```jldoctest extreme_filter
julia> Ω_mask = centered(Bool[1 1 0; 1 1 0; 1 0 0]) # custom neighborhood in mask format
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 1  1  0
 1  1  0
 1  0  0

julia> out = extreme_filter(max, M, Ω_mask)
5×5 Matrix{$Int}:
 4  8  6  9  4
 8  8  9  9  9
 8  8  9  9  9
 7  8  8  9  9
 6  6  6  5  7

julia> Ω_offsets = strel(CartesianIndex, Ω_mask) # custom neighborhood as displacement offset
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(-1, -1)
 CartesianIndex(0, -1)
 CartesianIndex(1, -1)
 CartesianIndex(-1, 0)

julia> out == extreme_filter(max, M, Ω_offsets) # both versions work equivalently
true
```

See also the in-place version [`extreme_filter!`](@ref). Another function in ImageFiltering
package `ImageFiltering.mapwindow` provides similar functionality.
"""
extreme_filter(f, A; r=nothing, dims=coords_spatial(A)) = extreme_filter(f, A, strel_box(A, dims; r))
extreme_filter(f, A, Ω::AbstractArray) = extreme_filter!(f, similar(A), A, Ω)

"""
    extreme_filter!(f, out, A; [r], [dims])
    extreme_filter!(f, out, A, Ω)

The in-place version of [`extreme_filter`](@ref) where `out` is the output array that gets
modified.
"""
extreme_filter!(f, out, A; r=nothing, dims=coords_spatial(A)) = extreme_filter!(f, out, A, strel_box(A, dims; r))
function extreme_filter!(f, out, A, Ω)
    axes(out) == axes(A) || throw(DimensionMismatch("axes(out) must match axes(A)"))
    require_select_function(f, eltype(A))
    return _extreme_filter!(strel_type(Ω), f, out, A, Ω)
end

_extreme_filter!(::MorphologySE, f, out, A, Ω) = _extreme_filter_generic!(f, out, A, Ω)
_extreme_filter!(::SEDiamond, f, out, A, Ω) = _extreme_filter_diamond!(f, out, A, Ω)
function _extreme_filter!(::MorphologySE, f::MAX_OR_MIN, out, A::AbstractArray{T}, Ω) where {T<:Union{Gray{Bool},Bool}}
    # NOTE(johnnychen94): empirical choice based on benchmark results (intel i9-12900k)
    true_ratio = gray(sum(A) / length(A)) # this usually takes <2% of the time but gives a pretty good hint to do the decision
    true_ratio == 1 && (out .= true; return out)
    true_ratio == 0 && (out .= false; return out)
    use_bool = prod(strel_size(Ω)) > 9 || (f === max && true_ratio > 0.8) || (f === min && true_ratio < 0.2)
    if use_bool
        return _extreme_filter_bool!(f, out, A, Ω)
    else
        return _extreme_filter_generic!(f, out, A, Ω)
    end
end
function _extreme_filter!(::SEDiamond, f::MAX_OR_MIN, out, A::AbstractArray{T}, Ω) where {T<:Union{Gray{Bool},Bool}}
    rΩ = strel_size(Ω) .÷ 2
    if ndims(A) == 2 && all(rΩ[1] .== rΩ)
        # super fast implementation and wins in all cases
        return _extreme_filter_diamond!(f, out, A, Ω)
    end

    # NOTE(johnnychen94): empirical choice based on benchmark results (intel i9-12900k)
    true_ratio = gray(sum(A) / length(A)) # this usually takes <2% of the time but gives a pretty good hint to do the decision
    true_ratio == 1 && (out .= true; return out)
    true_ratio == 0 && (out .= false; return out)
    use_bool = (f === max && true_ratio > 0.4) || (f === min && true_ratio < 0.6)
    if use_bool
        return _extreme_filter_bool!(f, out, A, Ω)
    else
        return _extreme_filter_diamond!(f, out, A, Ω)
    end
end

###
# Implementation details
###

function _extreme_filter_generic!(f, out, A, Ω)
    @debug "call the generic `extreme_filter` implementation" fname = _extreme_filter_generic!
    Ω = strel(CartesianIndex, Ω)
    δ = CartesianIndex(strel_size(Ω) .÷ 2)

    R = CartesianIndices(A)
    R_inner = (first(R) + δ):(last(R) - δ)

    # NOTE(johnnychen94): delibrately duplicate the kernel codes here for inner loop and
    # boundary loop, because keeping the big function body allows Julia to do further
    # optimization and thus significantly improves the overall performance.
    @inbounds for p in R_inner
        # for interior points, boundary check is unnecessary
        s = A[p]
        for o in Ω
            v = A[p + o]
            s = f(s, v)
        end
        out[p] = s
    end
    @inbounds for p in EdgeIterator(R, R_inner)
        # for edge points, skip if the offset exceeds the boundary
        s = A[p]
        for o in Ω
            q = p + o
            checkbounds(Bool, A, q) || continue
            s = f(s, A[q])
        end
        out[p] = s
    end
    return out
end

# optimized implementation for SEDiamond -- a typical case of separable filter
function _extreme_filter_diamond!(f, out, A, Ω::SEDiamondArray{N}) where {N}
    rΩ = strel_size(Ω) .÷ 2
    if N == 2 && f isa MAX_OR_MIN && all(rΩ[1] .== rΩ)
        iter = rΩ[1]
        return _extreme_filter_diamond_2D!(f, out, A, iter)
    else
        return _extreme_filter_diamond_generic!(f, out, A, Ω)
    end
end

function _extreme_filter_diamond_generic!(f, out, A, Ω::SEDiamondArray)
    @debug "call the optimized `extreme_filter` implementation for SEDiamond SE" fname = _extreme_filter_diamond!
    rΩ = strel_size(Ω) .÷ 2
    r = maximum(rΩ)

    # To avoid the result affected by loop order, we need two arrays
    src = (out === A) || (r > 1) ? copy(A) : A
    out .= src

    # applying radius=r filter is equivalent to applying radius=1 filter r times
    for i in 1:r
        Ω = strel_diamond(A, Ω.dims)
        rΩ = strel_size(Ω) .÷ 2
        inds = axes(A)
        for d in 1:ndims(A)
            # separately apply to each dimension
            if size(out, d) == 1 || rΩ[d] == 0
                continue
            end
            Rpre = CartesianIndices(inds[1:(d - 1)])
            Rpost = CartesianIndices(inds[(d + 1):end])
            _extreme_filter_C2!(f, out, src, Rpre, inds[d], Rpost)
        end

        if r > 1 && i < r
            src .= out
        end
    end
    return out
end

# fast version for C2 connectivity along each dimension
@noinline function _extreme_filter_C2!(f, dst, src, Rpre, inds, Rpost)
    # manually unroll the loop for r=1 case
    # for r>1 case, it is equivalently to apply r times
    @inbounds for Ipost in Rpost, Ipre in Rpre
        # loop inner region
        for i in (first(inds) + 1):(last(inds) - 1)
            a1 = src[Ipre, i - 1, Ipost]
            a2 = dst[Ipre, i, Ipost]
            a3 = src[Ipre, i + 1, Ipost]
            dst[Ipre, i, Ipost] = f(f(a1, a2), a3)
        end
        # process two edge points
        i = first(inds)
        a2 = dst[Ipre, i, Ipost]
        a3 = src[Ipre, i + 1, Ipost]
        dst[Ipre, i, Ipost] = f(a2, a3)

        i = last(inds)
        a1 = src[Ipre, i - 1, Ipost]
        a2 = dst[Ipre, i, Ipost]
        dst[Ipre, i, Ipost] = f(a1, a2)
    end
    return dst
end

# Shift vector of size N up or down by 1 accorded to dir, pad with v
const SafeArrayTypes{T,N,NA} = Union{Array{T,N},Base.SubArray{T,N,Array{T,NA}}}
function _unsafe_padded_copyto!(dest::SafeArrayTypes{T,1}, src::SafeArrayTypes{T,1}, dir, N, v) where {T}
    # ptr level optimized implementation for Real types
    # short-circuit all check so unsafe
    if dir
        unsafe_store!(pointer(dest), v)
        unsafe_copyto!(pointer(dest, 2), pointer(src, 1), N - 1)
    else
        unsafe_copyto!(pointer(dest, 1), pointer(src, 2), N - 1)
        unsafe_store!(pointer(dest), v, N)
    end
    return dest
end
function _unsafe_padded_copyto!(dest::AbstractVector, src::AbstractVector, dir, N, v)
    if dir
        dest[begin] = v
        copyto!(dest, 2, src, 1, N - 1)
    else
        dest[end] = v
        copyto!(dest, 1, src, 2, N - 1)
    end
    return dest
end

# ptr level optimized implementation for Real types
# short-circuit all check
# call LoopVectorization directly
function _unsafe_shift_arith!(
    f::MAX_OR_MIN,
    out::AbstractVector,
    tmp::AbstractVector,
    tmp2::AbstractVector,
    A::AbstractVector,
) #tmp external to reuse external allocation
    if f === min
        padd = typemax(eltype(out))
    else
        padd = typemin(eltype(out))
    end
    N = length(out)
    #in  src  = [1,2,3,4], padd, N=4
    #out tmp  = [padd,1,2,3]
    #out tmp2 = [1,2,3,padd]
    _unsafe_padded_copyto!(tmp, A, true, N, padd)
    _unsafe_padded_copyto!(tmp2, A, false, N, padd)
    @turbo warn_check_args = false for i in eachindex(out, A, tmp, tmp2)
        out[i] = f(f(A[i], tmp[i]), tmp2[i])
    end
    return out
end

# optimized `extreme_filter` implementation for SEDiamond SE and 2D images
# Similar approach could be found in
# Žlaus, Danijel & Mongus, Domen. (2018). In-place SIMD Accelerated Mathematical Morphology. 76-79. 10.1145/3206157.3206176.
# see https://www.researchgate.net/publication/325480366_In-place_SIMD_Accelerated_Mathematical_Morphology
function _extreme_filter_diamond_2D!(f::MAX_OR_MIN, out, A, iter)
    @debug "call the AVX-enabled `extreme_filter` implementation for 2D SEDiamond" fname = _extreme_filter_diamond_2D!
    out_actual = out

    out = OffsetArrays.no_offset_view(out)
    A = OffsetArrays.no_offset_view(A)

    # To avoid the result affected by loop order, we need two arrays
    src = (out === A) || (iter > 1) ? copy(A) : A
    # NOTE(johnnychen94): we don't need to do `out .= src` here because it is write-only;
    # all values are generated from the read-only `src`.

    # LoopVectorization currently doesn't understand Gray and N0f8 types, thus we
    # reinterpret to its raw data type UInt8
    src = rawview(channelview(src))
    out = rawview(channelview(out))

    # creating temporaries
    ySize, xSize = size(A)
    tmp = similar(out, eltype(out), ySize)
    tmp2 = similar(tmp)
    tmp3 = similar(tmp)

    # applying radius=r filter is equivalent to applying radius=1 filter r times
    for i in 1:iter
        # NOTE(johnnychen94): this implementation is essentially equivalent to the
        # following loop. Here we buffer the getindex to further improve the performance
        # by explicitly introducing SIMD buffers.

        # for y in 2:(size(src, 2) - 1), x in 2:(size(src, 1) - 1)
        #     tmp = f(f(src[x, y], src[x - 1, y]), src[x + 1, y])
        #     out[x, y] = f(f(tmp, src[x, y - 1]), src[x, y + 1])
        # end

        #compute first edge column
        #translate to clipped connection, x neighborhood of interest, ? neighborhood we don't care, . center
        # x ?
        # . x
        # x ?
        viewprevious = view(src, :, 1)
        # dilate/erode col 1
        _unsafe_shift_arith!(f, tmp2, tmp, tmp3, viewprevious)
        viewnext = view(src, :, 2)
        viewout = view(out, :, 1)
        # inf/sup between dilate/erode col 1 0 and col 2
        LoopVectorization.vmap!(f, viewout, viewnext, tmp2)
        #next->current
        viewcurrent = view(src, :, 2)
        for c in 3:xSize
            # viewprevious c-2
            # viewcurrent  c-1
            # viewnext     c
            # ? x ?
            # x . x
            # ? x ?
            viewout = view(out, :, c - 1)
            viewnext = view(src, :, c)
            # dilate(x-1)/erode(x-1)
            _unsafe_shift_arith!(f, tmp2, tmp, tmp3, viewcurrent)
            @turbo warn_check_args = false for i in eachindex(viewout, viewprevious, tmp2)
                #sup(x-2,dilate(x-1)),inf(x-2,erode(x-1))
                #sup(sup(x-2,dilate(x-1),x) || inf(inf(x-2,erode(x-1),x)
                viewout[i] = f(f(viewprevious[i], tmp2[i]), viewnext[i])
            end
            #current->previous
            viewprevious = view(src, :, c - 1)
            #next->current
            viewcurrent = view(src, :, c)
        end
        #end last column
        #translate to clipped connection
        # ? x
        # x .
        # ? x
        # dilate/erode col x
        viewout = view(out, :, xSize)
        _unsafe_shift_arith!(f, tmp2, tmp, tmp3, viewcurrent)
        LoopVectorization.vmap!(f, viewout, tmp2, viewprevious)
        if iter > 1 && i < iter
            src .= out
        end
    end

    return out_actual
end

# optimized implementation for Bool inputs with max/min select function
# 1) use &&, || instead of max, min
# 2) short-circuit the result to avoid unnecessary indexing and computation
function _extreme_filter_bool!(f, out, A::AbstractArray{Bool}, Ω)
    @debug "call the optimized max/min `extreme_filter` implementation for boolean array" fname = _extreme_filter_bool!
    Ω = strel(CartesianIndex, Ω)
    δ = CartesianIndex(strel_size(Ω) .÷ 2)

    R = CartesianIndices(A)
    R_inner = (first(R) + δ):(last(R) - δ)

    select = _fast_select(f)
    @inbounds for p in R_inner
        out[p] = select(A, p, Ω)
    end
    for p in EdgeIterator(R, R_inner)
        out[p] = select(A, p, Ω)
    end
    return out
end
function _extreme_filter_bool!(f, out, A::AbstractArray{Gray{Bool}}, Ω)
    return _extreme_filter_bool!(f, out, reinterpret(Bool, A), Ω)
end
_fast_select(::typeof(max)) = _maximum_fast
_fast_select(::typeof(min)) = _minimum_fast
Base.@propagate_inbounds function _maximum_fast(A::AbstractArray{Bool}, p, Ω)
    rst = A[p]
    for o in Ω
        # the complicated control flow ruins the SIMD and thus for small Ω, the performance
        # will be worse
        q = p + o
        @boundscheck checkbounds(Bool, A, q) || continue
        x = A[q]
        x && return true
        rst = rst || x
    end
    return rst
end
Base.@propagate_inbounds function _minimum_fast(A::AbstractArray{Bool}, p, Ω)
    rst = A[p]
    for o in Ω
        # the complicated control flow ruins the SIMD and thus for small Ω, the performance
        # will be worse
        q = p + o
        @boundscheck checkbounds(Bool, A, q) || continue
        x = A[q]
        x || return false
        rst = rst && x
    end
    return rst
end


text:
"""This function returns the successor of any given positive integer."""
function succ(n :: Int64) :: Int64
    return n + 1;
end;

"""This function implements the fast exponentiation algorithm"""
function fast_expo(n :: Int64, exp :: Int64) :: Int64
    if expo == 0
        return n
    elseif expo % 2 == 0
        let res = fast_expo(n, div(exp, 2))
        return (res * res)
    else
        let res = fast_expo(n, div(exp, 2))
        return (res * res * n)
    end
end


text:
using OptionPrice
using SDEModels
using Test

# write your own tests here

analytic = Analytic()
lewis = Lewis()
carr_madan = CarrMadan()

@test price(BlackScholes(0.01, 0.1), analytic, 100.0, 90.0, 0.01, 2.0) ≈ 13.147550110559813
@test price(BlackScholes(0.01, 0.1), carr_madan, 100.0, 90.0, 0.01, 2.0) ≈ 13.147550110559813
@test price(BlackScholes(0.01, 0.1), lewis, 100.0, 90.0, 0.01, 2.0) ≈ 13.147550110559813

@test price(Heston(0.01, 10, 0.16, 0.1, -0.8), carr_madan, [90.0, 0.16], 80.0, 0.01, 1.0) ≈ 19.489631391639065
@test price(Heston(0.01, 10, 0.16, 0.1, -0.8), lewis, [90.0, 0.16], 80.0, 0.01, 1.0) ≈ 19.489631391639065


text:
module TestFileFormats

using MathOptInterface
using Test

const MOI = MathOptInterface

function test_CBF()
    return include(joinpath(@__DIR__, "CBF", "CBF.jl"))
end

function test_LP()
    return include(joinpath(@__DIR__, "LP", "LP.jl"))
end

function test_MOF()
    return include(joinpath(@__DIR__, "MOF", "MOF.jl"))
end

function test_MPS()
    return include(joinpath(@__DIR__, "MPS", "MPS.jl"))
end

function test_NL()
    return include(joinpath(@__DIR__, "NL", "NL.jl"))
end

function test_SDPA()
    return include(joinpath(@__DIR__, "SDPA", "SDPA.jl"))
end

function test_copying()
    models = [
        MOI.FileFormats.CBF.Model,
        MOI.FileFormats.LP.Model,
        MOI.FileFormats.MOF.Model,
        MOI.FileFormats.MPS.Model,
        MOI.FileFormats.SDPA.Model,
    ]
    for src in models
        model_src = src()
        for dest in models
            model_dest = dest()
            MOI.copy_to(model_dest, model_src)
            @test !isempty(sprint(write, model_dest))
        end
        model_dest = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
        MOI.copy_to(model_dest, model_src)
    end
end

function test_instantiate()
    models = [
        MOI.FileFormats.CBF.Model,
        MOI.FileFormats.LP.Model,
        MOI.FileFormats.MOF.Model,
        MOI.FileFormats.MPS.Model,
        MOI.FileFormats.SDPA.Model,
    ]
    for src in models
        model = MOI.instantiate(src)
        @test model isa src
        bridged = MOI.instantiate(src; with_bridge_type = Float64)
        @test bridged isa MOI.Bridges.LazyBridgeOptimizer
    end
    return
end

function test_compressed_open()
    for cs in [MOI.FileFormats.Bzip2(), MOI.FileFormats.Gzip()]
        for open_type in ["a", "r+", "w+", "a+"]
            @test_throws ArgumentError MOI.FileFormats.compressed_open(
                (x) -> nothing,
                "dummy.gz",
                open_type,
                cs,
            )
        end
    end
end

function test_compression()
    filename = joinpath(@__DIR__, "MPS", "free_integer.mps")
    model = MOI.FileFormats.Model(filename = filename)
    MOI.read_from_file(model, filename)
    filename = joinpath(@__DIR__, "free_integer.mps")
    MOI.write_to_file(model, filename * ".garbage")
    for ext in ["", ".bz2", ".gz"]
        MOI.write_to_file(model, filename * ext)
        model2 = MOI.FileFormats.Model(filename = filename * ext)
        MOI.read_from_file(model2, filename)
    end

    sleep(1.0)  # Allow time for unlink to happen.
    for ext in ["", ".garbage", ".bz2", ".gz"]
        rm(filename * ext, force = true)
    end
    return
end

function test_Model()
    for (format, model) in [
        (MOI.FileFormats.FORMAT_CBF, MOI.FileFormats.CBF.Model()),
        (MOI.FileFormats.FORMAT_LP, MOI.FileFormats.LP.Model()),
        (MOI.FileFormats.FORMAT_MOF, MOI.FileFormats.MOF.Model()),
        (MOI.FileFormats.FORMAT_MPS, MOI.FileFormats.MPS.Model()),
        (MOI.FileFormats.FORMAT_NL, MOI.FileFormats.NL.Model()),
        (MOI.FileFormats.FORMAT_SDPA, MOI.FileFormats.SDPA.Model()),
    ]
        @test typeof(
            MOI.FileFormats.Model(format = format, filename = "foo.bar"),
        ) == typeof(model)
    end
    @test_throws(
        ErrorException(
            "When `format==FORMAT_AUTOMATIC` you must pass a `filename`.",
        ),
        MOI.FileFormats.Model(format = MOI.FileFormats.FORMAT_AUTOMATIC)
    )
    for (ext, model) in [
        (".cbf", MOI.FileFormats.CBF.Model()),
        (".lp", MOI.FileFormats.LP.Model()),
        (".mof.json", MOI.FileFormats.MOF.Model()),
        (".mps", MOI.FileFormats.MPS.Model()),
        (".nl", MOI.FileFormats.NL.Model()),
        (".sdpa", MOI.FileFormats.SDPA.Model()),
    ]
        @test typeof(MOI.FileFormats.Model(filename = "a$(ext)")) ==
              typeof(model)
        @test typeof(MOI.FileFormats.Model(filename = "a$(ext).gz")) ==
              typeof(model)
    end
    @test_throws(
        ErrorException("Unable to automatically detect format of a.b."),
        MOI.FileFormats.Model(filename = "a.b")
    )
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

end

TestFileFormats.runtests()


text:
function ChewMandestam(s,m1sq,m2sq)
    m1, m2 = sqrt(m1sq), sqrt(m2sq)
    #
    sth,spth = (m1+m2)^2, (m1-m2)^2;
    λh = sqrt(spth-s)*sqrt(sth-s)
    #
    val = 1/(π) * (
        λh/s*log((m1sq+m2sq-s+λh)/(2*m1*m2))-
            (m1sq-m2sq) * (1.0/s - 1/sth) * log(m1/m2)
    )
    return val
end

function Rho(s,m1sq,m2sq)
    m1, m2 = sqrt(m1sq), sqrt(m2sq)
    #
    sth,spth = (m1+m2)^2, (m1-m2)^2;
    λh = sqrt(s-sth)*sqrt(s-spth)
    return λh/s
end
iRho(s,m1sq,m2sq) = 1im*Rho(s,m1sq,m2sq)

#  _|  _|                                _|
#  _|      _|_|_|      _|_|      _|_|_|  _|_|_|      _|_|_|  _|_|_|      _|_|
#  _|  _|  _|    _|  _|_|_|_|  _|_|      _|    _|  _|    _|  _|    _|  _|_|_|_|
#  _|  _|  _|    _|  _|            _|_|  _|    _|  _|    _|  _|    _|  _|
#  _|  _|  _|    _|    _|_|_|  _|_|_|    _|    _|    _|_|_|  _|_|_|      _|_|_|
#                                                            _|
#                                                            _|

pole(σ,mcsq) = 1.0/(mcsq - σ)
BW(σ,m,Γ) = pole(σ,m^2-1im*m*Γ)

function BWdw(σ,m,Γ,m1,m2)
    √σ <= (m1+m2) && return 0.0im
    ρσ = sqrt(Kallen(σ,m1^2,m2^2))/σ
    ρ0 = sqrt(Kallen(m^2,m1^2,m2^2))/m^2
    pole(σ,m^2-1im*m*Γ*ρσ/ρ0)
end

struct Lineshape{T<:Real}
    itype::Int
    type::String
    pars::Vector{T}
end

# Breit-Wigner with constant width
Lineshape(m,Γ) = Lineshape(0,"BW",[m,Γ])
BreitWigner(m,Γ) = Lineshape(0,"BW",[m,Γ])

function amp(s,lsh::Lineshape)
    (lsh.itype == 0) && return BW(s,lsh.pars[1],lsh.pars[2])
    (lsh.itype == 9) && return ScattLenT11(s,lsh.pars[1:3],lsh.pars[4:5],lsh.pars[6:7])
    error("No itype #$(lsh.itype) found")
    return 0.0im
end

# Scattering length approximation
ScattLenApproximation(;
    Ms::Tuple{Real,Real,Real} = error("give Ms = (c11, c12, c22)"),
    ms1::Tuple{Real,Real} = error("give ms1 = (m1, m2)"),
    ms2::Tuple{Real,Real}) = Lineshape(9,"ScattLen",[Ms...,ms1...,ms2...])

function ScattLenMatrix(s,Ms,ms1,ms2)
    M  = SMatrix{2,2}(Ms[1],Ms[2],Ms[2],Ms[3])
    iρ = SMatrix{2,2}(-sqrt((ms1[1]+ms1[2])^2-s)*sqrt(s-(ms1[1]-ms1[2])^2)/s, 0.0im, 0.0im,
                      -sqrt((ms2[1]+ms2[2])^2-s)*sqrt(s-(ms2[1]-ms2[2])^2)/s)
    return inv(M - iρ)
end
ScattLenT11(s,Ms,ms1,ms2) = ScattLenMatrix(s,Ms,ms1,ms2)[1,1]

#            _|
#  _|_|_|    _|_|_|          _|_|_|  _|_|_|
#  _|    _|  _|    _|      _|_|      _|    _|
#  _|    _|  _|    _|          _|_|  _|    _|
#  _|_|_|    _|    _|  _|  _|_|_|    _|_|_|    _|
#  _|                                _|
#  _|                                _|

function RhoQTB(s,Xlineshape,msq; channel::Int=1)
    m1, m2, m3 = sqrt.(msq)
    #
    k = channel
    i,j = ij_from_k(k)
    #
    (√s < (m1+m2+m3) || √s ≈ (m1+m2+m3)) && return 0.0
    val = quadgk(σ->abs2(Xlineshape(σ))*sqrt(Kallen(s,σ,msq[k])*Kallen(σ,msq[i],msq[j]))/σ, (√msq[i]+√msq[j])^2, (√s-√msq[k])^2)[1]
    return val / s # /(2π*(8π)^2)
end

RhoQTB(s,m1,m2,Γ1,m1th) = RhoQTB(s,σ->BW(σ,m1,Γ1),[m2,m1th/2,m1th/2].^2)
iRhoQTB(s,m1,m2,Γ1,m1th) = 1im*RhoQTB(s,m1,m2,Γ1,m1th)

function three_body_phase_space_integral(function_σs, ms)
    m1sq,m2sq,m3sq,s = ms^2
    #
    σ1min,σ1max = lims1(ms)
    function integrand(x,f)
        σ1 = σ1min + x[1]*(σ1max-σ1min)
        z = 2*x[2]-1
        σs = Invariants(ms;σ1=σ1,σ3=σ3of1(z,σ1,ms^2))
        #
        r = function_σs(σs) *
            sqrt(Kallen(s,σ1, m1sq)*Kallen(σ1,m2sq,m3sq))/σ1 *
            (σ1max-σ1min) # jacobians
        f[1],f[2] = reim(r)
    end
    return complex(cuhre(integrand,2,2)[1]...) / s # /(2π*(8π)^2)
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module TestDecisionTree

using Test
import CategoricalArrays
using MLJBase

# load code to be tested:
import MLJModels
import DecisionTree
using MLJModels.DecisionTree_

# get some test data:
using RDatasets
iris = dataset("datasets", "iris")
const X = iris[:, 1:4]
const y = iris[:, 5]

baretree = DecisionTreeClassifier()

baretree.max_depth = 1
fitresult, cache, report = MLJBase.fit(baretree, 2, X, y);
baretree.max_depth = -1 # no max depth
fitresult, cache, report =
    MLJBase.update(baretree, 1, fitresult, cache, X, y);

# in this case decision tree is a perfect predictor:
yhat = MLJBase.predict_mode(baretree, fitresult, X);
@test yhat == y

# but pruning upsets this:
baretree.post_prune = true
baretree.merge_purity_threshold=0.1
fitresult, cache, report =
    MLJBase.update(baretree, 2, fitresult, cache, X, y)
yhat = MLJBase.predict_mode(baretree, fitresult, X);
@test yhat != y
yhat = MLJBase.predict(baretree, fitresult, X);

# check preservation of levels:
yyhat = predict_mode(baretree, fitresult, MLJBase.selectrows(X, 1:3))
@test MLJBase.classes(yyhat[1]) == MLJBase.classes(y[1])

info_dict(baretree)

# # testing machine interface:
# tree = machine(baretree, X, y)
# fit!(tree)
# yyhat = predict_mode(tree, MLJBase.selectrows(X, 1:3))
using Random: seed!
seed!(0)

n,m = 10^3, 5 ;
raw_features = rand(n,m);
weights = rand(-1:1,m);
labels = raw_features * weights;
features = MLJBase.table(raw_features);

R1Tree = DecisionTreeRegressor(min_samples_leaf=5, pruning_purity_threshold=0.1)
R2Tree = DecisionTreeRegressor(min_samples_split=5)
model1, = MLJBase.fit(R1Tree,1, features, labels)

vals1 = MLJBase.predict(R1Tree,model1,features)
R1Tree.post_prune = true
model1_prune, = MLJBase.fit(R1Tree,1, features, labels)
vals1_prune = MLJBase.predict(R1Tree,model1_prune,features)
@test vals1 !=vals1_prune

@test DecisionTree.R2(labels, vals1) > 0.8

model2, = MLJBase.fit(R2Tree, 1, features, labels)
vals2 = MLJBase.predict(R2Tree, model2, features)
@test DecisionTree.R2(labels, vals2) > 0.8

info_dict(R1Tree)

end
true


text:
# RolloutSimulator
# maintained by @zsunberg

"""
    RolloutSimulator(rng, max_steps)
    RolloutSimulator(; <keyword arguments>)

A fast simulator that just returns the reward

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

# Keyword arguments:
- rng: A random number generator to use.
- eps: A small number; if γᵗ where γ is the discount factor and t is the time step becomes smaller than this, the simulation will be terminated.
- max_steps: The maximum number of steps to simulate.

# Usage (optional arguments in brackets):

    ro = RolloutSimulator()
    history = simulate(ro, pomdp, policy, [updater [, init_belief [, init_state]]])

See also: [`HistoryRecorder`](@ref), [`run_parallel`](@ref)
"""
struct RolloutSimulator{RNG<:AbstractRNG} <: Simulator
    rng::RNG

    # optional: if these are null, they will be ignored
    max_steps::Union{Nothing,Int}
    eps::Union{Nothing,Float64}
end

RolloutSimulator(rng::AbstractRNG, d::Int=typemax(Int)) = RolloutSimulator(rng, d, nothing)
function RolloutSimulator(;rng=Random.GLOBAL_RNG,
                           eps=nothing,
                           max_steps=nothing)
    return RolloutSimulator{typeof(rng)}(rng, max_steps, eps)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy) begin
    @req updater(::typeof(policy))
    bu = updater(policy)
    @subreq simulate(sim, pomdp, policy, bu)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, bu::Updater) begin
    @req initialstate_distribution(::typeof(pomdp))
    dist = initialstate_distribution(pomdp)
    @subreq simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initialstate_distribution(pomdp)
    return @inferred simulate(sim, pomdp, policy, bu, dist)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief) begin
    @req rand(::typeof(sim.rng), ::typeof(initial_belief))
    @subreq simulate(sim, pomdp, policy, updater, initial_belief, s)
end

function simulate(sim::RolloutSimulator, pomdp::POMDP{S}, policy::Policy, updater::Updater, initial_belief) where {S}
    s = rand(sim.rng, initial_belief)::S
    return @inferred simulate(sim, pomdp, policy, updater, initial_belief, s)
end

@POMDP_require simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s) begin
    P = typeof(pomdp)
    S = statetype(P)
    A = actiontype(P)
    O = obstype(P)
    @req initialize_belief(::typeof(updater), ::typeof(initial_belief))
    @req isterminal(::P, ::S)
    @req discount(::P)
    @req generate_sor(::P, ::S, ::A, ::typeof(sim.rng))
    b = initialize_belief(updater, initial_belief)
    @req action(::typeof(policy), ::typeof(b))
    @req update(::typeof(updater), ::typeof(b), ::A, ::O)
end

function simulate(sim::RolloutSimulator, pomdp::POMDP, policy::Policy, updater::Updater, initial_belief, s)
    if isterminal(pomdp, s)
        return 0.0
    end

    if sim.eps == nothing
        eps = 0.0
    else
        eps = sim.eps
    end

    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end

    disc = 1.0
    r_total = 0.0

    b = initialize_belief(updater, initial_belief)

    step = 1

    while disc > eps && step <= max_steps
        a = action(policy, b)

        sp, o, r = generate_sor(pomdp, s, a, sim.rng)

        r_total += disc*r

        s = sp

        bp = update(updater, b, a, o)
        b = bp

        disc *= discount(pomdp)
        step += 1

        if isterminal(pomdp, s, a, o)
            break
        end
    end

    return r_total
end

@POMDP_require simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy) begin
    istate = initialstate(mdp, sim.rng)
    @subreq simulate(sim, mdp, policy, istate)
end

@POMDP_require simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy, initialstate) begin
    P = typeof(mdp)
    S = typeof(initialstate)
    A = actiontype(mdp)
    @req isterminal(::P, ::S)
    @req action(::typeof(policy), ::S)
    @req generate_sr(::P, ::S, ::A, ::typeof(sim.rng))
    @req discount(::P)
end

function simulate(sim::RolloutSimulator, mdp::MDP, policy::Policy)
    istate = initialstate(mdp, sim.rng)
    @inferred simulate(sim, mdp, policy, istate)
end

# This function is also defined in Entropy for an InformationRewardBeliefMDP
function simulate(sim::RolloutSimulator, mdp::Union{MDP{S}, POMDP{S}}, policy::Policy, initialstate::S) where {S}
    if isterminal(mdp, initialstate)
        return 0.0
    end

    if sim.eps == nothing
        eps = 0.0
    else
        eps = sim.eps
    end

    if sim.max_steps == nothing
        max_steps = typemax(Int)
    else
        max_steps = sim.max_steps
    end

    s = initialstate

    disc = 1.0
    r_total = 0.0
    step = 1

    while disc > eps && step <= max_steps
        a = action(policy, s)

        sp, r = generate_sr(mdp, s, a, sim.rng)

        r_total += disc*r

        s = sp

        disc *= discount(mdp)
        step += 1

        if isterminal(mdp, s, a)
            break
        end
    end

    return r_total
end