text:
############################################################
## joAbstractDAparallelLinearOperator(s) - overloaded Base functions

# eltype(jo)
eltype(A::joAbstractDAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = promote_type(DDT,RDT)

# deltype(jo)
deltype(A::joAbstractDAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = DDT

# reltype(jo)
reltype(A::joAbstractDAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = RDT

# show(jo)
function show(A::joAbstractDAparallelLinearOperator)
    println("Type: $(typeof(A).name)")
    println("Name: $(A.name)")
    println("Size: $(size(A))")
    println(" NVC: $(A.nvc)")
    println(" DDT: $(deltype(A))")
    println(" RDT: $(reltype(A))")
    return nothing
end

# display(jo)
display(A::joAbstractDAparallelLinearOperator) = println((typeof(A),((A.m,A.n),A.nvc),A.name))

# size(jo)
size(A::joAbstractDAparallelLinearOperator) = A.m,A.n

# size(jo,1/2)
function size(A::joAbstractDAparallelLinearOperator,ind::Integer)
    if ind==1
		return A.m
	elseif ind==2
		return A.n
	elseif ind==3
		return A.nvc
	else
		throw(joAbstractDAparallelLinearOperatorException("invalid index"))
	end
end

# length(jo)
length(A::joAbstractDAparallelLinearOperator) = A.m*A.n

# jo_full(jo)

# norm(jo)

# real(jo)

# imag(jo)

# conj(jo)
conj(A::joDALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDALinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.PAs_out, A.fclean, A.rclean
        )
conj(A::joDAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributedLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.PAs_out, A.fclean, A.rclean
        )
conj(A::joDAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributingLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_out, A.gclean
        )
conj(A::joDAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAgatheringLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.gclean
        )

# transpose(jo)
transpose(A::joDALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDALinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
transpose(A::joDAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributedLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
transpose(A::joDAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAgatheringLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.gclean
        )
transpose(A::joDAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributingLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_in, A.gclean
        )

# adjoint(jo)
adjoint(A::joDALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDALinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
adjoint(A::joDAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributedLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
adjoint(A::joDAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAgatheringLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.gclean
        )
adjoint(A::joDAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joDAdistributingLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_in, A.gclean
        )

# isreal(jo)
isreal(A :: joAbstractDAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = (DDT<:Real && RDT<:Real)

# issymmetric(jo)

# ishermitian(jo)

# getindex(jo,...)

############################################################
## overloaded Base *(...jo...)

# *(jo,jo)
function *(A::joDALinOpsUnion{CDT,ARDT,2},B::joDALinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joDALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joDALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.fclean,A.rclean)
end
function *(A::joDALinOpsUnion{CDT,ARDT,2},B::joDAdistributeLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joDALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joDAdistributingLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               A.PAs_out,B.gclean)
end
function *(A::joDAgatherLinOpsUnion{CDT,ARDT,2},B::joDALinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joDAgatheringLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joDAgatheringLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joDAgatheringLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joDAgatheringLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,B.fclean)
end
function *(A::joDAdistributeLinOpsUnion{CDT,ARDT,2},B::joDAgatherLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joDALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    return joDALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.gclean,A.gclean)
end
function *(A::joDAgatherLinOpsUnion{CDT,ARDT,2},B::joDAdistributeLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joLinearOperator{BDDT,ARDT}("($(A.name)*$(B.name))",size(A,1),size(B,2),
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF)
end
# extras with warnings
function *(A::joDAdistribute{CDT,ARDT,2},B::joDAgather{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joDALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.PAs_out.dims[2]==B.PAs_in.dims[2] || throw(joDALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    @warn "*($(typeof(A)),$(typeof(B))) is a senseless operation. Get rid of it!"
    return joDALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.PAs_out.dims[2],
               v1->v1, v2->v2, v3->v3, v4->v4,
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.gclean,A.gclean)
end
function *(A::joDAgather{CDT,ARDT,2},B::joDAdistribute{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.PAs_in.dims[2]==B.PAs_out.dims[2] || throw(joLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    @warn "*($(typeof(A)),$(typeof(B))) is a senseless operation. Get rid of it!"
    return joLinearOperator{BDDT,ARDT}("($(A.name)*$(B.name))",size(A,1),size(B,2),
               v1->v1, v2->v2, v3->v3, v4->v4,
               @joNF, @joNF, @joNF, @joNF)
end
# joDA with joAbstractLinearOperator
function *(A::joDAdistributeLinOpsUnion{CDT,ARDT,2},B::joAbstractLinearOperator{BDDT,CDT}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joAbstractLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    return joDAdistributingLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               A.PAs_out,A.gclean)
end
function *(A::joAbstractLinearOperator{CDT,ARDT},B::joDAgatherLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joAbstractLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    return joDAgatheringLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),B.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,B.gclean)
end

# *(jo,mvec)
function *(A::joDALinearOperator{ADDT,ARDT,2},mv::DArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joDALinearOperatorException("shape mismatch"))
    A.nvc == size(mv,2) || throw(joDALinearOperatorException("shape mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joDALinearOperatorException("*($(A.name),mv::DArray): input distributor mismatch"))
    MV=A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joDAdistributedLinearOperator{ADDT,ARDT,2},mv::DArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joDAdistributedLinearOperatorException("shape mismatch"))
    A.nvc == size(mv,2) || throw(joDAdistributedLinearOperatorException("shape mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joDAdistributedLinearOperatorException("*($(A.name),mv::DArray): input distributor mismatch"))
    MV=dalloc(A.PAs_out)
    spmd(joDAutils.jo_x_mv!,A.fop,mv,MV,pids=A.PAs_out.procs)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joDAdistributingLinearOperator{ADDT,ARDT,2},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joDAdistributingLinearOperatorException("shape mismatch in A$(size(A))*v$(size(mv))"))
    A.nvc == size(mv,2) || throw(joDAdistributingLinearOperatorException("nvc size mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    MV = A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joDAgatheringLinearOperator{ADDT,ARDT,2},mv::DArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joDAgatheringLinearOperatorException("shape mismatch in A$(size(A))*v$(size(mv))"))
    A.nvc == size(mv,2) || throw(joDAgatheringLinearOperatorException("nvc size mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joDAgatheringLinearOperatorException("*($(A.name),mv::DArray): input distributor mismatch"))
    MV = A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end

# *(mvec,jo)

# *(jo,vec)

# *(vec,jo)

# *(num,jo)

# *(jo,num)

############################################################
## overloaded Base \(...jo...)

# \(jo,jo)

# \(jo,mvec)

# \(mvec,jo)

# \(jo,vec)

# \(vec,jo)

# \(num,jo)

# \(jo,num)
#\{ADDT,ARDT}(A::joAbstractDAparallelLinearOperator{ADDT,ARDT},a::Number) = inv(a)*A
#\{ADDT,ARDT}(A::joAbstractDAparallelLinearOperator{ADDT,ARDT},a::joNumber{ADDT,ARDT}) = inv(a)*A

############################################################
## overloaded Base +(...jo...)

# +(jo)
+(A::joAbstractDAparallelLinearOperator) = A

# +(jo,jo)

# +(jo,mvec)

# +(mvec,jo)

# +(jo,vec)

# +(vec,jo)

# +(jo,num)

# +(num,jo)
+(b::Number,A::joAbstractDAparallelLinearOperator{ADDT,ARDT}) where {ADDT,ARDT} = A+b
+(b::joNumber{ADDT,ARDT},A::joAbstractDAparallelLinearOperator{ADDT,ARDT}) where {ADDT,ARDT} = A+b

############################################################
## overloaded Base -(...jo...)

# -(jo)

# -(jo,jo)
-(A::joAbstractDAparallelLinearOperator,B::joAbstractDAparallelLinearOperator) = A+(-B)

# -(jo,mvec)

# -(mvec,jo)

# -(jo,vec)

# -(vec,jo)

# -(jo,num)
-(A::joAbstractDAparallelLinearOperator,b::Number) = A+(-b)
-(A::joAbstractDAparallelLinearOperator,b::joNumber) = A+(-b)

# -(num,jo)
-(b::Number,A::joAbstractDAparallelLinearOperator) = -A+b
-(b::joNumber,A::joAbstractDAparallelLinearOperator) = -A+b

############################################################
## overloaded Base .*(...jo...)

# .*(num,jo)
Base.Broadcast.broadcasted(::typeof(*),a::Number,A::joAbstractDAparallelLinearOperator) = a*A
Base.Broadcast.broadcasted(::typeof(*),a::joNumber,A::joAbstractDAparallelLinearOperator) = a*A

# .*(jo,num)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractDAparallelLinearOperator,a::Number) = a*A
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractDAparallelLinearOperator,a::joNumber) = a*A

############################################################
## overloaded Base .\(...jo...)

############################################################
## overloaded Base .+(...jo...)

# .+(jo,num)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractDAparallelLinearOperator,b::Number) = A+b
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractDAparallelLinearOperator,b::joNumber) = A+b

# .+(num,jo)
Base.Broadcast.broadcasted(::typeof(+),b::Number,A::joAbstractDAparallelLinearOperator) = A+b
Base.Broadcast.broadcasted(::typeof(+),b::joNumber,A::joAbstractDAparallelLinearOperator) = A+b

############################################################
## overloaded Base .-(...jo...)

# .-(jo,num)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractDAparallelLinearOperator,b::Number) = A+(-b)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractDAparallelLinearOperator,b::joNumber) = A+(-b)

# .-(num,jo)
Base.Broadcast.broadcasted(::typeof(-),b::Number,A::joAbstractDAparallelLinearOperator) = -A+b
Base.Broadcast.broadcasted(::typeof(-),b::joNumber,A::joAbstractDAparallelLinearOperator) = -A+b

############################################################
## overloaded Base block methods

# hcat(...jo...)

# vcat(...jo...)

# hvcat(...jo...)

############################################################
## overloaded LinearAlgebra functions

# mul!(...,jo,...)

# ldiv!(...,jo,...)


text:
insts = OrderedDict()
insts["minimal"] = [
    ((4, true, false, false),),
    ((4, false, true, false),),
    ((4, false, false, true),),
    ]
insts["fast"] = [
    ((10, true, false, false),),
    ((10, false, true, false),),
    ((10, false, false, true),),
    ((100, true, false, false),),
    ((100, false, true, false), (default_tol_relax = 100,)),
    ((100, false, false, true),),
    ((1000, true, false, false),),
    ]
insts["various"] = [
    ((750, true, false, false),),
    ((750, false, false, true),),
    ((1500, true, false, false),),
    ((1500, false, false, true),),
    ]
return (MaxVolumeNative, insts)

ext: jl


text:
##################################################
#   Override default handling (default = inline) #
##################################################

get_latex_function(r::DataFrames.DataFrame) = mdtable
get_md_function(args::DataFrames.DataFrame) = mdtable

###############################################
#         Overload environment functions      #
###############################################

function mdtable(d::DataFrames.DataFrame; kwargs...)
    body = Matrix(d)
    head = propertynames(d)
    mdtable(body; head=head, kwargs...)
end

function _latextabular(d::DataFrames.DataFrame; kwargs...)
    body = Matrix(d)
    head = propertynames(d)
    _latextabular(body; head=head, kwargs...)
end

function _latexarray(d::DataFrames.DataFrame; kwargs...)
    body = Matrix(d)
    head = permutedims(propertynames(d))
    result = vcat(head, body)
    _latexarray(result; kwargs...)
end


text:
#

julia> using Checkers

#tests the conditional property of Julia
 julia> @test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#100 test will be run
#Passes if 100 tests meet the crierea of holding true. Will fail becuase generating 100 pairs of x,y, some will bound to fail condition

#another test the conditional property of Julia
false.
#100,000 tests will performed until 100 passing test
julia> @test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#Same as top but will pass because 100000test will be performed.

#another test the conditional property of Julia

julia> @test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#Same as top but will pass because 100000test will be performed.


text:
#################################
"""
    fgmax = loadfgmax(loaddir::AbstractString, fg::VisClaw.FixedGrid; nval_save::Integer=fg.nval)

Function: fgmaxXXXX.txt reader
"""
function loadfgmax(loaddir::AbstractString, fg::VisClaw.FixedGrid; nval_save::Integer=fg.nval)
    ## check nval_save
    nval_save > fg.nval && (nval_save=fg.nval)

    ## fgmaxXXXX.txt
    filename = "fgmax"*@sprintf("%04d",fg.id)*".txt"
    ## load
    dat = readdlm(joinpath(loaddir, filename), Float64)
    dat[dat.<-1e10] .= NaN

    ncol = 4 + 2(fg.nval) + 1
    ## assign
    if fg.style == 0 || fg.style == 1
        topo = dat[:,4]
        D = dat[:,5]
        tD = dat[:,5+fg.nval]
        tarrival = dat[:,end]
        if fg.nval >= 2
            v = dat[:,6]
            tv = dat[:,6+fg.nval]
        end
        if fg.nval >= 5
            M = dat[:,7]
            tM = dat[:,7+fg.nval]
            Mflux = dat[:,8]
            tMflux = dat[:,8+fg.nval]
            hmin = dat[:,9]
            thmin = dat[:,9+fg.nval]
        end
    elseif fg.style == 2 || fg.style == 3
        valall = permutedims(reshape(dat, (fg.nx, fg.ny, ncol)), [2 1 3])
        topo = valall[:,:,4]
        D = valall[:,:,5]
        tD = valall[:,:,5+fg.nval]
        tarrival = valall[:,:,end]
        if fg.nval >= 2
            v = valall[:,:,6]
            tv = valall[:,:,6+fg.nval]
        end
        if fg.nval >= 5
            M = valall[:,:,7]
            tM = valall[:,:,7+fg.nval]
            Mflux = valall[:,:,8]
            tMflux = valall[:,:,8+fg.nval]
            hmin = valall[:,:,9]
            thmin = valall[:,:,9+fg.nval]
        end
    elseif fg.style == 4
        indc = [[fg.flag[i][1] fg.flag[i][2]] for i=1:fg.npts]
        indc = vcat(indc...)
        ind = sortperm(indc[:,1])

        topo = dat[ind,4]
        D = dat[ind,5]
        tD = dat[ind,5+fg.nval]
        tarrival = dat[ind,end]
        if fg.nval >= 2
            v = dat[ind,6]
            tv = dat[ind,6+fg.nval]
        end
        if fg.nval >= 5
            M = dat[ind,7]
            tM = dat[ind,7+fg.nval]
            Mflux = dat[ind,8]
            tMflux = dat[ind,8+fg.nval]
            hmin = dat[ind,9]
            thmin = dat[ind,9+fg.nval]
        end
    end

    if nval_save == 1;     fgmax = VisClaw.FGmax(topo,D,tD,tarrival)
    elseif nval_save == 2; fgmax = VisClaw.FGmax(topo,D,v,tD,tv,tarrival)
    elseif nval_save == 5; fgmax = VisClaw.FGmax(topo,D,v,M,Mflux,hmin,tD,tv,tM,tMflux,thmin,tarrival)
    end

    # return
    return fgmax

end
#################################


text:
# Source/Sink with NamedTuple, both row and column oriented
"A default row-oriented \"table\" that supports both the `Data.Source` and `Data.Sink` interfaces. Can be used like `Data.stream!(source, Data.RowTable)`. It is represented as a Vector of NamedTuples."
const RowTable{T} = Vector{T} where {T <: NamedTuple}
Data.isdone(source::RowTable, row, col, rows, cols) = row > rows || col > cols
function Data.isdone(source::RowTable, row, col)
    rows = length(source)
    return Data.isdone(source, row, col, rows, rows > 0 ? length(rows[1]) : 0)
end
Data.streamtype(::Type{Array}, ::Type{Data.Field}) = true
Data.streamfrom(source::RowTable, ::Type{Data.Field}, ::Type{T}, row, col) where {T} = source[row][col]
Data.streamtypes(::Type{Array}) = [Data.Row]
Data.accesspattern(::RowTable) = Data.RandomAccess()

Data.streamto!(sink::RowTable, ::Type{Data.Row}, val, row, col) =
    row > length(sink) ? push!(sink, val) : setindex!(sink, val, row)
Data.streamto!(sink::RowTable, ::Type{Data.Row}, val, row, col::Int, ::Type{Val{false}}) =
    push!(sink, val)
Data.streamto!(sink::RowTable, ::Type{Data.Row}, val, row, col::Int, ::Type{Val{true}}) =
    setindex!(sink, val, row)

# Basically, a NamedTuple with any # of AbstractVector elements, accessed by column name
"A default column-oriented \"table\" that supports both the `Data.Source` and `Data.Sink` interfaces. Can be used like `Data.stream!(source, Data.Table). It is represented as a NamedTuple of AbstractVectors.`"
const Table = NamedTuple{names, T} where {names, T <: NTuple{N, AbstractVector{S} where S}} where {N}

function Data.schema(rt::RowTable{NamedTuple{names, T}}) where {names, T}
    return Data.Schema(Type[A for A in T.parameters],
                        collect(map(string, names)), length(rt))
end
# NamedTuple Data.Source implementation
# compute Data.Schema on the fly
function Data.schema(df::NamedTuple{names, T}) where {names, T}
    return Data.Schema(Type[eltype(A) for A in T.parameters],
                        collect(map(string, names)), length(df) == 0 ? 0 : length(getfield(df, 1)))
end

Data.isdone(source::Table, row, col, rows, cols) = row > rows || col > cols
function Data.isdone(source::Table, row, col)
    cols = length(source)
    return Data.isdone(source, row, col, cols == 0 ? 0 : length(getfield(source, 1)), cols)
end

# We support both kinds of streaming
Data.streamtype(::Type{<:NamedTuple}, ::Type{Data.Column}) = true
Data.streamtype(::Type{<:NamedTuple}, ::Type{Data.Field}) = true
Data.accesspattern(::Table) = Data.RandomAccess()

# Data.streamfrom is pretty simple, just return the cell or column
@inline Data.streamfrom(source::Table, ::Type{Data.Column}, T, row::Integer, col::Integer) = source[col]
@inline Data.streamfrom(source::Table, ::Type{Data.Field}, T, row::Integer, col::Integer) = source[col][row]

# NamedTuple Data.Sink implementation
# we support both kinds of streaming to our type
Data.streamtypes(::Type{<:NamedTuple}) = [Data.Column, Data.Field]
# we support streaming WeakRefStrings
Data.weakrefstrings(::Type{<:NamedTuple}) = true

# convenience methods for "allocating" a single column for streaming
allocate(::Type{T}, rows, ref) where {T} = Vector{T}(undef, rows)
# allocate(::Type{T}, rows, ref) where {T <: Union{CategoricalValue, Missing}} =
#     CategoricalArray{CategoricalArrays.unwrap_catvalue_type(T)}(rows)
# special case for WeakRefStrings
allocate(::Type{WeakRefString{T}}, rows, ref) where {T} = WeakRefStringArray(ref, WeakRefString{T}, rows)
allocate(::Type{Union{WeakRefString{T}, Missing}}, rows, ref) where {T} = WeakRefStringArray(ref, Union{WeakRefString{T}, Missing}, rows)
allocate(::Type{Missing}, rows, ref) = fill(missing, rows)

# NamedTuple doesn't allow duplicate names, so make sure there are no duplicates in our column names
function makeunique(names::Vector{String})
    nms = [Symbol(nm) for nm in names]
    seen = Set{Symbol}()
    for (i, x) in enumerate(nms)
        x in seen ? setindex!(nms, Symbol("$(x)_$i"), i) : push!(seen, x)
    end
    return (nms...,)
end

function Array(sch::Data.Schema{R}, ::Type{Data.Row}, append::Bool=false, args...) where {R}
    types = Data.types(sch)
    # check if we're dealing with an existing NamedTuple sink or not
    if !isempty(args) && args[1] isa RowTable && types == Data.types(Data.schema(args[1]))
        sink = args[1]
        sinkrows = size(Data.schema(sink), 1)
        if append && !R
            sch.rows = sinkrows
        else
            newsize = ifelse(!R, 0, ifelse(append, sinkrows + sch.rows, sch.rows))
            resize!(sink, newsize)
            sch.rows = newsize
        end
    else
        rows = ifelse(!R, 0, sch.rows)
        names = makeunique(Data.header(sch))
        # @show rows, names, types
        sink = Vector{NamedTuple{names, Tuple{types...}}}(undef, rows)
        sch.rows = rows
    end
    return sink
end

function Array(sink::RowTable, sch::Data.Schema, ::Type{Data.Row}, append::Bool)
    return Array(sch, Data.Row, append, sink)
end

# Construct or modify a NamedTuple to be ready to stream data from a source with a schema of `sch`
# We support WeakRefString streaming, so we include the `reference` keyword
function NamedTuple(sch::Data.Schema{R}, ::Type{S}=Data.Field,
                    append::Bool=false, args...; reference::Vector{UInt8}=UInt8[]) where {R, S <: StreamType}
    types = Data.types(sch)
    # check if we're dealing with an existing NamedTuple sink or not
    if !isempty(args) && args[1] isa Table && types == Data.types(Data.schema(args[1]))
        # passing in an existing NamedTuple Sink w/ same types as source (as indicated by `sch`)
        sink = args[1]
        sinkrows = size(Data.schema(sink), 1)
        if append && (S == Data.Column || !R) # are we appending and either column-streaming or there are an unknown # of rows
            sch.rows = sinkrows
            # dont' need to do anything because:
              # for Data.Column, we just append columns anyway (see Data.streamto! below)
              # for Data.Field, the # of rows in the source are unknown (ismissing(rows)), so we'll just push! in streamto!
        else
            # need to adjust the existing sink
            # similar to above, for Data.Column or unknown # of rows for Data.Field, we'll append!/push!, so we empty! the columns
            # if appending, we want to grow our columns to be able to include every row in source (sinkrows + sch.rows)
            # if not appending, we're just "re-using" a sink, so we just resize it to the # of rows in the source
            newsize = ifelse(S == Data.Column || !R, 0, ifelse(append, sinkrows + sch.rows, sch.rows))
            foreach(col->resize!(col, newsize), sink)
            sch.rows = newsize
        end
        # take care of a possible reference from source by letting WeakRefStringArrays hold on to them
        if !isempty(reference)
            foreach(col-> col isa WeakRefStringArray && push!(col.data, reference), sink)
        end
    else
        # allocating a fresh NamedTuple Sink; append is irrelevant
        # for Data.Column or unknown # of rows in Data.Field, we only ever append!, so just allocate empty columns
        rows = ifelse(S == Data.Column, 0, ifelse(!R, 0, sch.rows))
        names = makeunique(Data.header(sch))

        sink = NamedTuple{names}(Tuple(allocate(types[i], rows, reference) for i = 1:length(types)))
        sch.rows = rows
    end
    return sink
end

# Constructor that takes an existing NamedTuple sink, just pass it to our mega-constructor above
function NamedTuple(sink::Table, sch::Data.Schema, ::Type{S}, append::Bool; reference::Vector{UInt8}=UInt8[]) where {S}
    return NamedTuple(sch, S, append, sink; reference=reference)
end

# Data.streamto! is easy-peasy, if there are known # of rows from source, we pre-allocated
# so we can just set the value; otherwise (didn't pre-allocate), we push!/append! the values
Data.streamto!(sink::Table, ::Type{Data.Field}, val, row, col::Int) =
    (A = getfield(sink, col); row > length(A) ? push!(A, val) : setindex!(A, val, row))
Data.streamto!(sink::Table, ::Type{Data.Field}, val, row, col::Int, ::Type{Val{false}}) =
    push!(getfield(sink, col), val)
Data.streamto!(sink::Table, ::Type{Data.Field}, val, row, col::Int, ::Type{Val{true}}) =
    getfield(sink, col)[row] = val
Data.streamto!(sink::Table, ::Type{Data.Column}, column, row, col::Int, knownrows) =
    append!(getfield(sink, col), column)

# Row iteration for Data.Sources
struct Rows{S, NT}
    source::S
end

Base.eltype(rows::Rows{S, NT}) where {S, NT} = NT
function Base.length(rows::Rows)
    sch = Data.schema(rows.source)
    return size(sch, 1)
end

"Returns a NamedTuple-iterator of any `Data.Source`"
function rows(source::S) where {S}
    sch = Data.schema(source)
    names = makeunique(Data.header(sch))
    types = Data.types(sch)
    return Rows{S, NamedTuple{names, Tuple{types...}}}(source)
end

function Base.iterate(rows::Rows{S, NamedTuple{names, types}}, row::Int=1) where {S, names, types}
    if @generated
        cols = length(names)
        vals = Tuple(:(Data.streamfrom(rows.source, Data.Field, $typ, row, $col)) for (col, typ) in zip(1:cols, types.parameters) )
        r = :(($(NamedTuple{names, types})(($(vals...),)), row + 1))
        # println(r)
        return quote
            Data.isdone(rows.source, row, $cols) && return nothing
            $r
        end
    else
        cols = length(names)
        return NamedTuple{names, types}(Tuple(Data.streamfrom(rows.source, Data.Field, T, row, i) for (i, T) in enumerate(types.parameters)))
    end
end


text:
function generate_html(integrators, figures, figures_detail, prefix, suffix, fig_dir, output_dir, output_file;
                        title="", params=nothing, max_height="200px")
    header = """
    <!DOCTYPE html>

            <title>GeometricExamples.jl $title</title>
            <meta charset="utf-8" />
            <meta name="authors" content="Michael Kraus" />

    """

    footer = """

    """

    # 1. Header
    # 2. Simulation information (runid, Δt, nt, ...)
    # 3. Figures with caption: integrator, tableau
    # 4. Libray settings (solver, atol, rtol, stol, ...)
    # 5. Footer

    # generate empty output string
    html = ""

    # add header
    html = html * header

    # add title
    if title != ""
        html = html * "
" * title * "
\n"
    end

    # write simulation parameters
    if params != nothing
        html = html * "
\n"
        html = html * "
\n"
        html = html * "
\n"
        html = html * "
\n"
        html = html * "Simulation Parameters	\n"
        html = html * "	\n"
        html = html * "
\n"
        html = html * "
\n"
        html = html * "
\n"

        for (key, value) in params
            html = html * "
\n"
            html = html * "" * string(key) * "	\n"
            html = html * "" * string(value) * "	\n"
            html = html * "
\n"
        end

        html = html * "
\n"
        html = html * "
\n"
        html = html * "

\n"
    end

    # start table
    html = html * "
\n"

    # table header
    html = html * "
\n"
    html = html * "
\n"
    html = html * "Run ID	\n"
    html = html * "Integrator	\n"
    html = html * "Tableau	\n"

    for fig in figures
        html = html * "	\n"
    end

    html = html * "
\n"
    html = html * "
\n"

    # plots for each integrator/tableau combination
    html = html * "
\n"
    for run in integrators

        int  = string(run[1])[(lastindex(findlast("Integrator", string(run[1]))) + 1):end]
        tab  = string(run[2].name)
        rid  = prefix * run[3]

        html = html * "
\n"
        html = html * "" * "" * rid * "" * "	\n"
        html = html * "" * int * "	\n"
        html = html * "" * tab * "	\n"

        for fig in figures
            fig_file = fig_dir * "/" * rid * fig * suffix

            html = html * ""
            if isfile(output_dir * "/" * fig_file)
                html = html * ""
                html = html * ""
                html = html * ""
            end
            html = html * "	\n"
        end

        html = html * "
\n"

        # generate empty output string for run file
        html_run = ""

        # add header
        html_run = html_run * header

        # add title
        if title != ""
            html_run = html_run * "
" * title * "
\n"
        end

        for fig in figures_detail
            fig_file = rid * fig * suffix

            if isfile(output_dir * "/" * fig_dir * "/" * fig_file)
                html_run = html_run * "\n"
                html_run = html_run * "
\n"
                html_run = html_run * fig_file * "

\n\n"
            end
        end

        # add footer
        html_run = html_run * footer

        # write file
        open(output_dir * "/" * rid * ".html", "w") do f
            write(f, html_run)
        end

    end
    html = html * "
\n"

    # finish table
    html = html * "
\n"

    # write library settings
    settings = get_config_dictionary()

    html = html * "
\n"
    html = html * "
\n"
    html = html * "
\n"
    html = html * "
\n"
    html = html * "Library Settings	\n"
    html = html * "	\n"
    html = html * "
\n"
    html = html * "
\n"
    html = html * "
\n"

    for key in sort(collect(keys(settings)))
        html = html * "
\n"
        html = html * "" * string(key) * "	\n"
        html = html * "" * string(settings[key]) * "	\n"
        html = html * "
\n"
    end

    html = html * "
\n"
    html = html * "
\n"
    html = html * "

\n"

    # add footer
    html = html * footer

    # write file
    open(output_dir * "/" * output_file, "w") do f
        write(f, html)
    end

end


text:
"""
Package: Forecast

    co2(full = false)

Return dataset with atmospheric Carbon Dioxide Dry Air Mole Fractions from quasi-continuous measurements at Mauna Loa, Hawaii.

K.W. Thoning, A.M. Crotwell, and J.W. Mund (2020), Atmospheric Carbon Dioxide Dry Air Mole Fractions from continuous measurements at Mauna Loa, Hawaii, Barrow, Alaska, American Samoa and South Pole. 1973-2019, Version 2020-08 National Oceanic and Atmospheric Administration (NOAA), Global Monitoring Laboratory (GML), Boulder, Colorado, USA https://doi.org/10.15138/yaf1-bk21 FTP path: ftp://aftp.cmdl.noaa.gov/data/greenhouse_gases/co2/in-situ/surface/

# Arguments
- `full`: if `true` Returns full original dataset from 1973 to 2020 in a DataFrame, otherwise returns the subset used in "STL: A Seasonal-Trend Decomposition Procedure Based on Loess" from Cleveland et. al. Its default value is `false`.

# Returns
 Dataframe or TimeArray containing the descrived dataset.

# Examples
```julia-repl
julia> co2()
[ Info: Dataset used in Cleveland et al. paper
4609×1 TimeArray{Union{Missing, Float64},1,Date,Array{Union{Missing, Float64},1}}
1974-05-17 to 1986-12-31
│            │ co2     │
├────────────┼─────────┤
│ 1974-05-17 │ 333.38  │
│ 1974-05-18 │ 333.11  │
│ 1974-05-19 │ 333.46  │
   [...]
```
"""
function co2(full = false)

    data = "data/co2.csv.gz"
    path = joinpath(splitdir(@__DIR__)[1], data)

    co2_df = GZip.open(path, "r") do io
        CSV.read(io,DataFrame)
    end

    if full
        @info "Full dataset from 1973 to 2020"
        return co2_df
    else
        @info "Dataset used in Cleveland et al. paper"
    end

    # recode missing values
    co2_df.value = replace(co2_df.value, -999.99 => missing)

    dates_co2 = Dates.Date(1973,1,1):Dates.Day(1):Dates.Date(2019,12,31)
    co2_ts = TimeSeries.TimeArray(dates_co2, co2_df.value)

    # paper dates
    # using 17 may instead 17 april (as the paper claims)
    # since there seems to be no data for april in 1974
    dates_co2_stl = Dates.Date(1974,5,17):Dates.Day(1):Dates.Date(1986,12,31)

    # revome leap year day
    dates_co2_stl = filter(dates_co2_stl) do x
        (Dates.isleapyear(x) & (Dates.month(x) == 2)) ? Dates.day(x) != 29 : true
    end

    co2_stl_ts = co2_ts[dates_co2_stl]
    TimeSeries.rename!(co2_stl_ts,:co2)
    co2_stl_ts

end


text:
using DashBootstrapComponents

tabs = html_div([
    dbc_tabs([
        dbc_tab(label = "Tab 1", tab_style = Dict("marginLeft" => "auto")),
        dbc_tab(label = "Tab 2", label_style = Dict("color" => "#00AEF9")),
    ]),
    html_br(),
    dbc_tabs([
        dbc_tab(label = "Tab 1", tabClassName = "ms-auto"),
        dbc_tab(label = "Tab 2", labelClassName = "text-success"),
    ]),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module EmpiricalMeasures

using Distributions
import Distributions:DiscreteMultivariateDistribution
import Distributions:AliasTable
using Random
using ArraysOfArrays
using StatsBase

export MvDiscreteNonParametric
export MvDiscreteNonParametricSampler
export empiricalmeasure
export support, probs
export mean, var, cov, pdf, entropy

struct MvDiscreteNonParametric{T <: Real,
    P <: Real,Ts <:  ArrayOfSimilarArrays{T},Ps <: AbstractVector{P}} <: DiscreteMultivariateDistribution

    support::Ts
    p::Ps

    function MvDiscreteNonParametric{T,P,Ts,Ps}(support::Ts,
        p::Ps) where {T <: Real,P <: Real,Ts <: AbstractVector{<:AbstractVector{T}},Ps <: AbstractVector{P}}

        length(support) == length(p) || error("length of `support` and `p` must be equal")
        isprobvec(p) || error("`p` must be a probability vector")
        allunique(support) || error("`support` must contain only unique value")
        new{T,P,Ts,Ps}(support, p)
    end
end

"""
    MvDiscreteNonParametric(
        support::AbstractVector,
        p::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
    )

Construct a multivariate discrete nonparametric probability distribution with `support` and corresponding
probabilities `p`. If the probability vector argument is not passed, then
equal probability is assigned to each entry in the support.

# Examples
```julia
using ArraysOfArrays
# rows correspond to samples
μ = MvDiscreteNonParametric(nestedview(rand(7,3)'))

# columns correspond to samples
ν = MvDiscreteNonParametric(nestedview(rand(7,3)))
```
"""
function MvDiscreteNonParametric(
    support::AbstractVector{<:AbstractVector{<:Real}},
    p::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
)
    return MvDiscreteNonParametric{eltype(eltype(support)),eltype(p),typeof(ArrayOfSimilarArrays(support)),typeof(p)}(
        ArrayOfSimilarArrays(support), p)
end

"""
    MvDiscreteNonParametric(
        support::Matrix{<:Real},
        p::AbstractVector{<:Real}=fill(inv(length(support)), length(support)
    )

Construct a multivariate discrete nonparametric probability distribution
from a matrix as `support` where each row is a sample, and corresponding
probabilities `p`. If the probability vector argument is not passed, then
equal probability is assigned to each entry in the support.

# Examples
```julia
# the rows correspond to the samples
using LinearAlgebra
μ = MvDiscreteNonParametric(rand(10,3), normalize!(rand(10),1))
```
"""
function MvDiscreteNonParametric(
    support::Matrix{<:Real},
    p::AbstractVector{<:Real}=fill(inv(size(support)[1]), size(support)[1])
)
    return MvDiscreteNonParametric(nestedview(support'), p)
end

Base.eltype(::Type{<:MvDiscreteNonParametric{T}}) where T = T

"""
    support(d::MvDiscreteNonParametric)
Get a sorted AbstractVector defining the support of `d`.
"""
Distributions.support(d::MvDiscreteNonParametric) = d.support

"""
    probs(d::MvDiscreteNonParametric)
Get the vector of probabilities associated with the support of `d`.
"""
Distributions.probs(d::MvDiscreteNonParametric) = d.p

# It would be more intuitive if length was the
#
"""
    length(d::MvDiscreteNonParametric)
Retunrs the dimension of the mass points (samples).
It corresponds to `innersize(d.support)[1]`,
where `innersize` is a function from from ArraysOfArrays.jl.
"""
Base.length(d::MvDiscreteNonParametric) = innersize(d.support)[1]

"""
    size(d::MvDiscreteNonParametric)
Returns the size of the support as a tuple where
the first value is the number of points (samples)
and the second is the dimension of the samples (e.g ℝⁿ).
It corresponds to `size(flatview(d.support)')`
where `flatview` is a function from from ArraysOfArrays.jl
that turns an Array of Arrays into a matrix.
"""
Base.size(d::MvDiscreteNonParametric) = size(flatview(d.support)')

"""
    empiricalmeasure(
        support::AbstractVector,
        probs::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
    )

Construct a finite discrete probability measure with `support` and corresponding
`probabilities`. If the probability vector argument is not passed, then
equal probability is assigned to each entry in the support.

# Examples
```julia
using ArraysOfArrays
# rows correspond to samples
μ = empiricalmeasure(nestedview(rand(7,3)'), normalize!(rand(10),1))

# columns correspond to samples, each with equal probability
ν = empiricalmeasure(nestedview(rand(3,12)))
```

!!! note
    If `support` is a 1D vector, the constructed measure will be sorted,
    e.g. for `mu = empiricalmeasure([3, 1, 2],[0.5, 0.2, 0.3])`, then
    `mu.support` will be `[1, 2, 3]` and `mu.p` will be `[0.2, 0.3, 0.5]`.
    Also, avoid passing 1D distributions as `RowVecs(rand(3))` or `[[1],[3],[4]]`,
    since this will be dispatched to the multivariate case instead
    of the univariate case for which the algorithm is more efficient.

"""
function empiricalmeasure(
    support::AbstractVector{<:Real},
    p::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
    )
    return DiscreteNonParametric(support, p)
end
function empiricalmeasure(
    support::AbstractVector{<:AbstractVector{<:Real}},
    p::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
    )
    return MvDiscreteNonParametric(support, p)
    # return MvDiscreteNonParametric{eltype(eltype(support)),eltype(p),typeof(support),typeof(p)}(support, p)
end

"""
    empiricalmeasure(
        support::Matrix{<:Real},
        probs::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
    )
Construct a multivariate empirical measure from a matrix
as `support` and sample by rows.
# Examples
```julia
using ArraysOfArrays
# the rows correspond to the samples
μ = empiricalmeasure(rand(7,3), normalize!(rand(10),1))
```
"""
function empiricalmeasure(
    support::Matrix{<:Real},
    p::AbstractVector{<:Real}=fill(inv(size(support)[1]), size(support)[1])
)
    return MvDiscreteNonParametric(nestedview(support'), p)
end

"""
    MvDiscreteNonParametricSampler(support, p)
Data structure for efficiently sampling from an arbitrary probability mass
function defined by `support` and probabilities `p`.
"""
struct MvDiscreteNonParametricSampler{T <: Real,S <: AbstractVector{<:AbstractVector{T}},A <: AliasTable} <: Sampleable{Multivariate,Discrete}
    support::S
    aliastable::A

    function MvDiscreteNonParametricSampler{T,S}(support::S, probs::AbstractVector{<:Real}
                                                ) where {T <: Real,S <: AbstractVector{<:AbstractVector{T}}}
        aliastable = AliasTable(probs)
        new{T,S,typeof(aliastable)}(support, aliastable)
    end
end

MvDiscreteNonParametricSampler(support::S, p::AbstractVector{<:Real}
                              ) where {T <: Real,S <: AbstractVector{<:AbstractVector{T}}} =
    MvDiscreteNonParametricSampler{T,S}(support, p)

# Sampling

sampler(d::MvDiscreteNonParametric) =
    MvDiscreteNonParametricSampler(support(d), probs(d))

function Distributions._rand!(rng::AbstractRNG, d::MvDiscreteNonParametric, x::AbstractVector{T}) where T <: Real

    length(x) == length(d) || throw(DimensionMismatch("Invalid argument dimension."))
    s = d.support
    p = d.p

    n = length(p)
    draw = Base.rand(rng, float(eltype(p)))
    cp = p[1]
    i = 1
    while cp <= draw && i < n
        @inbounds cp += p[i += 1]
    end
    for (j, v) in enumerate(s[i])
        x[j] = v
    end
    return x
end

function Distributions._logpdf(d::MvDiscreteNonParametric, x::AbstractVector{T}) where T <: Real
    s = support(d)
    p = probs(d)
    for i in 1:length(p)
        if s[i] == x
            return log(p[i])
        end
    end
    return log(zero(eltype(p)))
end

function Distributions.mean(d::MvDiscreteNonParametric)
    return StatsBase.mean(d.support, weights(d.p))
end

function Distributions.var(d::MvDiscreteNonParametric)
    x = support(d)
    p = probs(d)
    return StatsBase.var(x, Weights(p, one(eltype(p))), corrected=false)
end

function Distributions.cov(d::MvDiscreteNonParametric)
    x = support(d)
    p = probs(d)
    return cov(x, Weights(p, one(eltype(p))), corrected=false)
end

Distributions.entropy(d::MvDiscreteNonParametric) = entropy(probs(d))
Distributions.entropy(d::MvDiscreteNonParametric, b::Real) = entropy(probs(d), b)

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.jl

A pure Julia implementation of the [apache arrow](https://arrow.apache.org/) memory format specification.

This implementation supports the 1.0 version of the specification, including support for:
  * All primitive data types
  * All nested data types
  * Dictionary encodings, nested dictionary encodings, and messages
  * Extension types
  * Streaming, file, record batch, and replacement and isdelta dictionary messages
  * Buffer compression/decompression via the standard LZ4 frame and Zstd formats

It currently doesn't include support for:
  * Tensors or sparse tensors
  * Flight RPC
  * C data interface

Third-party data formats:
  * csv and parquet support via the existing [CSV.jl](https://github.com/JuliaData/CSV.jl) and [Parquet.jl](https://github.com/JuliaIO/Parquet.jl) packages
  * Other [Tables.jl](https://github.com/JuliaData/Tables.jl)-compatible packages automatically supported ([DataFrames.jl](https://github.com/JuliaData/DataFrames.jl), [JSONTables.jl](https://github.com/JuliaData/JSONTables.jl), [JuliaDB.jl](https://github.com/JuliaData/JuliaDB.jl), [SQLite.jl](https://github.com/JuliaDatabases/SQLite.jl), [MySQL.jl](https://github.com/JuliaDatabases/MySQL.jl), [JDBC.jl](https://github.com/JuliaDatabases/JDBC.jl), [ODBC.jl](https://github.com/JuliaDatabases/ODBC.jl), [XLSX.jl](https://github.com/felipenoris/XLSX.jl), etc.)
  * No current Julia packages support ORC or Avro data formats

See docs for official Arrow.jl API with the [User Manual](@ref) and reference docs for [`Arrow.Table`](@ref), [`Arrow.write`](@ref), and [`Arrow.Stream`](@ref).
"""
module Arrow

using Mmap
import Dates
using DataAPI, Tables, SentinelArrays, PooledArrays, CodecLz4, CodecZstd, TimeZones, BitIntegers

export ArrowTypes

using Base: @propagate_inbounds
import Base: ==

const DEBUG_LEVEL = Ref(0)

function setdebug!(level::Int)
    DEBUG_LEVEL[] = level
    return
end

function withdebug(f, level)
    lvl = DEBUG_LEVEL[]
    try
        setdebug!(level)
        f()
    finally
        setdebug!(lvl)
    end
end

macro debug(level, msg)
    esc(quote
        if DEBUG_LEVEL[] >= $level
            println(string("DEBUG: ", $(QuoteNode(__source__.file)), ":", $(QuoteNode(__source__.line)), " ", $msg))
        end
    end)
end

const FILE_FORMAT_MAGIC_BYTES = b"ARROW1"
const CONTINUATION_INDICATOR_BYTES = 0xffffffff

# vendored flatbuffers code for now
include("FlatBuffers/FlatBuffers.jl")
using .FlatBuffers

include("metadata/Flatbuf.jl")
using .Flatbuf; const Meta = Flatbuf

include("ArrowTypes/src/ArrowTypes.jl")
using .ArrowTypes
include("utils.jl")
include("arraytypes/arraytypes.jl")
include("eltypes.jl")
include("table.jl")
include("write.jl")
include("append.jl")
include("show.jl")

const LZ4_FRAME_COMPRESSOR = LZ4FrameCompressor[]
const ZSTD_COMPRESSOR = ZstdCompressor[]

function __init__()
    for _ = 1:Threads.nthreads()
        zstd = ZstdCompressor(; level=3)
        CodecZstd.TranscodingStreams.initialize(zstd)
        push!(ZSTD_COMPRESSOR, zstd)
        lz4 = LZ4FrameCompressor(; compressionlevel=4)
        CodecLz4.TranscodingStreams.initialize(lz4)
        push!(LZ4_FRAME_COMPRESSOR, lz4)
    end
    OBJ_METADATA_LOCK[] = ReentrantLock()
    return
end

end  # module Arrow


text:
# log2(x) =  vgetexp(x) + log2(vgetmant8(x))
@inline function vlog_v1(x1::VectorizationBase.AbstractSIMD{W,Float64}) where {W} # Testing if an assorted mix of operations
    x14 = vgetmant(x1)
    x8 = vgetexp(1.3333333333333333*x1)
    # x2 = reinterpret(UInt64, x1)
    # x3 = x2 >>> 0x0000000000000020
    # greater_than_zero = x1 > zero(x1)
    # alternative = VectorizationBase.ifelse(x1 == zero(x1), -Inf, NaN)
    # isinf = x1 == Inf
    # x5 = x3 + 0x0000000000095f62
    # x6 = x5 >>> 0x0000000000000014
    # x7 = x6 - 0x00000000000003ff
    # x8 = convert(Float64, x7 % Int)
    # @show x3 x5 x6 x7 x8
    # x9 = x5 << 0x0000000000000020
    # x10 = x9 & 0x000fffff00000000
    # x11 = x10 + 0x3fe6a09e00000000
    # x12 = x2 & 0x00000000ffffffff
    # x13 = x11 | x12
    # x14 = reinterpret(Float64, x13)
    x15 = x14 - 1.0
    x18 = x14 + 1.0
    x16 = x15 * x15
    x17 = 0.5 * x16
    x19 = x15 / x18
    x20 = x19 * x19
    x21 = x20 * x20
    x22 = vfmadd(x21, 0.15313837699209373, 0.22222198432149784)
    x23 = vfmadd(x21, x22, 0.3999999999940942)
    x24 = x23 * x21
    x25 = vfmadd(x21, 0.14798198605116586, 0.1818357216161805)
    x29 = x24 + x17
    x26 = vfmadd(x21, x25, 0.2857142874366239)
    x27 = vfmadd(x21, x26, 0.6666666666666735)
    x28 = x27 * x20
    x30 = x29 + x28
    x31 = x8 * 1.9082149292705877e-10
    x32 = vfmadd(x19, x30, x31)
    x33 = x15 - x17
    x34 = x33 + x32
    x35 = vfmadd(x8, 0.6931471803691238, x34)
    return x35
    # @show x35
    # x36 = VectorizationBase.ifelse(greater_than_zero, x35, alternative)
    # VectorizationBase.ifelse(isinf, Inf, x36)
end

const LOGTABLE = [0.0,0.007782140442054949,0.015504186535965254,0.02316705928153438,
    0.030771658666753687,0.0383188643021366,0.0458095360312942,0.053244514518812285,
    0.06062462181643484,0.06795066190850775,0.07522342123758753,0.08244366921107459,
    0.08961215868968714,0.09672962645855111,0.10379679368164356,0.11081436634029011,
    0.11778303565638346,0.12470347850095724,0.13157635778871926,0.13840232285911913,
    0.1451820098444979,0.15191604202584197,0.15860503017663857,0.16524957289530717,
    0.17185025692665923,0.1784076574728183,0.184922338494012,0.19139485299962947,
    0.19782574332991987,0.2042155414286909,0.21056476910734964,0.21687393830061436,
    0.22314355131420976,0.22937410106484582,0.2355660713127669,0.24171993688714516,
    0.24783616390458127,0.25391520998096345,0.25995752443692605,0.26596354849713794,
    0.27193371548364176,0.2778684510034563,0.2837681731306446,0.28963329258304266,
    0.2954642128938359,0.3012613305781618,0.3070250352949119,0.3127557100038969,
    0.3184537311185346,0.324119468654212,0.329753286372468,0.3353555419211378,
    0.3409265869705932,0.34646676734620857,0.3519764231571782,0.3574558889218038,
    0.3629054936893685,0.3683255611587076,0.37371640979358406,0.37907835293496944,
    0.38441169891033206,0.3897167511400252,0.394993808240869,0.4002431641270127,
    0.4054651081081644,0.4106599249852684,0.415827895143711,0.42096929464412963,
    0.4260843953109001,0.4311734648183713,0.43623676677491807,0.4412745608048752,
    0.44628710262841953,0.45127464413945856,0.4562374334815876,0.46117571512217015,
    0.46608972992459924,0.470979715218791,0.4758459048699639,0.4806885293457519,
    0.4855078157817008,0.4903039880451938,0.4950772667978515,0.4998278695564493,
    0.5045560107523953,0.5092619017898079,0.5139457511022343,0.5186077642080457,
    0.5232481437645479,0.5278670896208424,0.5324647988694718,0.5370414658968836,
    0.5415972824327444,0.5461324375981357,0.5506471179526623,0.5551415075405016,
    0.5596157879354227,0.564070138284803,0.5685047353526688,0.5729197535617855,
    0.5773153650348236,0.5816917396346225,0.5860490450035782,0.5903874466021763,
    0.5947071077466928,0.5990081896460834,0.6032908514380843,0.6075552502245418,
    0.6118015411059929,0.616029877215514,0.6202404097518576,0.6244332880118935,
    0.6286086594223741,0.6327666695710378,0.6369074622370692,0.6410311794209312,
    0.6451379613735847,0.6492279466251099,0.6533012720127457,0.65735807270836,
    0.661398482245365,0.6654226325450905,0.6694306539426292,0.6734226752121667,
    0.6773988235918061,0.6813592248079031,0.6853040030989194,0.689233281238809,
    0.6931471805599453];

logb(::Type{Float32},::Val{2})  = 1.4426950408889634
logb(::Type{Float32},::Val{:ℯ}) = One()
logb(::Type{Float32},::Val{10}) = 0.4342944819032518
logbU(::Type{Float64},::Val{2})  = 1.4426950408889634
logbL(::Type{Float64},::Val{2})  = 2.0355273740931033e-17
logbU(::Type{Float64},::Val{:ℯ}) = One()
logbL(::Type{Float64},::Val{:ℯ}) = Zero()
logbU(::Type{Float64},::Val{10}) = 0.4342944819032518
logbL(::Type{Float64},::Val{10}) = 1.098319650216765e-17

@inline function vlog_base(v::AbstractSIMD{W,Float64}) where {W}
    y = vgetmant12(v)
    mf = vgetexp(v)

    y128 = 128.0 * y
    f128 = vsreduce(y128, Val(0))
    F128 = y128 - f128# - 128.0
    jp = convert(UInt,F128) - 0x0000000000000080 # - 128
    # jp = convert(UInt,F128 - 128.0) # - 128
    # jp = convert(UInt,F128) - 0x000000000000007f # - 127

    hi = vload(zstridedpointer(LOGTABLE), (jp,))
    # l_hi = muladd(0.6931471805601177, mf, hi)
    l_hi = muladd(0.6931471805599453, mf, hi)

    u = (2.0*f128)/(y128+F128)
    # u = (2.0*f128)*vinv_fast(y128+F128)
    v = u*u
    q = u*v*muladd(0.012500053168098584, v, 0.08333333333303913)

    ## Step 4
    m_hi = logbU(Float64, Val(:ℯ))
    m_lo = logbL(Float64, Val(:ℯ))
    return fma(m_hi, l_hi, fma(m_hi, (u + q), m_lo*l_hi))
end

const LOG_TABLE_1 = Vec((
    Core.VecElement(-0.6931471805599453), Core.VecElement(-0.6325225587435105), Core.VecElement(-0.5753641449035618), Core.VecElement(-0.5212969236332861),
    Core.VecElement(-0.4700036292457356), Core.VecElement(-0.42121346507630353), Core.VecElement(-0.3746934494414107), Core.VecElement(-0.33024168687057687)
))
const LOG_TABLE_2 = Vec((
    Core.VecElement(-0.2876820724517809), Core.VecElement(-0.24686007793152578), Core.VecElement(-0.2076393647782445), Core.VecElement(-0.16989903679539747),
    Core.VecElement(-0.13353139262452263), Core.VecElement(-0.09844007281325252), Core.VecElement(-0.06453852113757118), Core.VecElement(-0.0317486983145803)
))
const LOG2_TABLE_1 = Vec((
    Core.VecElement(-1.0), Core.VecElement(-0.9125371587496606), Core.VecElement(-0.8300749985576876), Core.VecElement(-0.7520724865564145),
    Core.VecElement(-0.6780719051126377),Core.VecElement(-0.6076825772212398), Core.VecElement(-0.5405683813627028),Core.VecElement(-0.4764380439429871)
))
const LOG2_TABLE_2 = Vec((
    Core.VecElement(-0.4150374992788438),Core.VecElement(-0.3561438102252753),Core.VecElement(-0.2995602818589078), Core.VecElement(-0.24511249783653147),
    Core.VecElement(-0.19264507794239588),Core.VecElement(-0.14201900487242788),Core.VecElement(-0.09310940439148147),Core.VecElement(-0.04580368961312479)
))

@inline function logkern_5(x)
    # c5 = 1.442695040910480151153610022148007414933621426294604120324193438507895065110221
    # c4 = -0.72134752044717115185862113438227555323895953631804827391698847566536447158971
    # c3 = 0.4808981720904872241142144122992626169360496403735493070672998980156544684558863
    # c2 = -0.360673694648970133332498335009073996964643657971752173300330580172269152150781
    # c1 = 0.288894281140228357047827707854224408808934973997456175085011693761970856737952
    # c0 = -0.2406712024835816804281715317202731627491220621210104128295077305793044245650939
    # c5 = 1.000000000004163234554825918991486005405540498226964649982446185635197444806964
    # c4 = -0.5000000000218648309257667907050168790946684146132007063399540958947352513227305
    # c3 = 0.3333332566080562266977844691232384812586343066695585749657855454332169568535108
    # c2 = -0.2499998582263749067900056202432438929369945607183640079498376650719435347943278
    # c1 = 0.2002094490921391961924703669397489236496809044839667749701044680362284866088264
    # c0 = -0.1669110816794161712236614512036245351950242871767130610770379611369082994542161
    c5 = 1.442695040894969685206229605569847427762291596137745860941223179497930808392396
    c4 = -0.7213475204760259868264417322489103259481854442059435485772582433076339753304261
    c3 = 0.4808982362718110098786870255779326384270826960362097549483653590161925169305216
    c2 = -0.3606735556861350010148197635190486357722791582949450655689225850005157925747294
    c1 = 0.2888411793443405953870653225945083997249994212173033378047806528684724370915621
    c0 = -0.2408017898083064242160316054616972180231905999993934707245482866675760802361285
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
end
@inline function logkern_8(x)
    c8 = 1.442695040888963407816086035203962841602383176205260263758551036257573600888366
    c7 = -0.7213475204444770815507673261267005370242876005369959065788590988093119140411675
    c6 = 0.4808983469629686528945376971132664305807537133089662498292066034260696561965772
    c5 = -0.360673760285340441741375667112629490399269780822657391947176916709523570522447
    c4 = 0.2885390083116887609918806358733183651489107200432558827730599551390737454696275
    c3 = -0.2404489409194524614303773871505793252012695382714584724417731408820873961880021
    c2 = 0.20609895474146851962537542011143008396458423224465328706920087552924788983302
    c1 = -0.180654273148821532237028827440264938703627661154871403363931834846318112573294
    c0 = 0.1606521926313221723948643034905571707676419356400609690156705532267686492918868
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6),x,c7),x,c8)
end
@inline function logkern_9(x)
    c9 = 1.442695040888963407587958876357516873223860332904434711906078224014586198339343
    c8 = -0.7213475204444817057711203479236040464424888160067823280588699959470059217833771
    c7 = 0.480898346962976128170017519528240603877475137860740135494448418612334521964808
    c6 = -0.3606737602222044528926485178925468152698544632914287402177742224438034341355266
    c5 = 0.2885390082734153359755295484533201676037233143490110170185338945645707433116462
    c4 = -0.2404491736638805517094434611266213341107664043116792359649609367170386448753448
    c3 = 0.2060990174486639706727491347594731948792605354335332681880496880277396108033051
    c2 = -0.1803364995021978052782243890647317433194251060644842283875753942662516598663569
    c1 = 0.1606200865414345528500280443621751155379142697854250695344975262087060612426921
    c0 = -0.1446222903636375678552838798467050326611974502493221959482776774515654211467748
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6),x,c7),x,c8),x,c9)
end
@inline function logkern_7(x)
    c7 = 1.442695040888962266311320415441181965256745101864924017704730474290825179713386
    c6 = -0.7213475204444734849561916107798427186725025515630934214552341603017790847666743
    c5 = 0.4808983470003734483568608991659360631657097156445816019788134536213757250523677
    c4 = -0.3606737603148037624174654830835163329225507561910479104865837838596482628691627
    c3 = 0.2885388167992364548202154388340846682584745440281798969580530860342799719518947
    c2 = -0.2404488805785497044335968015014331102915764909299868970339917618519049206072763
    c1 = 0.2064127286335890836287774329056506792116396261746810041599847620260037837288407
    c0 = -0.1806895812056507095939337124225687685035687284931465242170306355644855623118632
    # c7 = 1.442695040888962266311320415440673226691509089399247319200433080789551133947874
    # c6 = -0.7213475204444734849561916106463536879222835723990934589483454213207143172699565
    # c5 = 0.4808983470003734483568609116686822486155905301580259810441289744019151890859293
    # c4 = -0.3606737603148037624174668433501054879718352818773442452949795967039102009618056
    # c3 = 0.2885388167992364548201727628199776711909314676794911979761485405529591339257216
    # c2 = -0.2404488805785497044302647238681858349839732924674030562755638231437855429968676
    # c1 = 0.2064127286335890836637376093489490181845837682333985654448331246682934490703086
    # c0 = -0.1806895812056507118633103913062383463992924441916921974716424429045368328209497
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6),x,c7)
end
@inline function logkern_6(x)
    c6 = 1.44269504088896112481574913926580271258562668309654298701097420877374680031478
    c5 = -0.7213475204628791729494574154147857065200754915843724028524643433219030089552769
    c4 = 0.4808983470214135184345854322580473442742138878734777669129377480116915710051349
    c3 = -0.3606736095355376179878548883602333628571188158521650924319649055538192403728833
    c2 = 0.2885387593457860540619732799086814986482101737149964973666904711488634555687257
    c1 = -0.2407576763559816924463087687841796294903125418141104155734727426854999960321288
    c0 = 0.2064519502023243513328580164119739425968179159770070826936661162470222233739379
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5),x,c6)
end

# Probably the best
@inline function vlog2_fast(a::AbstractSIMD{W,Float64}) where {W}
    # log2(vgetmant(a, Val(8))) + vgetexp(a) == log2(a)
    m = vgetmant(a, Val(8)) # m ∈ [1,2)
    e = vgetexp(a)
    # m ∈ [1,2); We reduce further...
    y = inv_approx(m) # y ∈ (0.5,1]
    # return y
    y₀ = vroundscale(y, Val(80)) # 80 = 16*(1+4)
    # y₀ is a multiple of 32; 0.5:(1/32):1.0
    # r + 1 = m*y;
    # log(r+1) = log(m) + log(y₀)
    # log(m) = log(1+r) - log(y₀)
    r = vfmsub(m, y₀, 1.0)
    # log1pr = logkern_5(r)
    # log1pr = logkern_6(r)
    log1pr = logkern_7(r)
    # log1pr = logkern_9(r)
    # @show r log1pr m
    # log1pr = logkern_8(r)
    y₀notone = y₀ ≠ 1.0
    inds = (reinterpret(UInt, y₀) >>> 48) & 0x0f
    # @show y₀ inds r
    # logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
    logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
    # emlogy₀ = e - logy₀
    emlogy₀ = ifelse(y₀notone, e - logy₀, e)
    log1pr + emlogy₀
    # logm = ifelse(y₀notone, log1pr - logy₀, log1pr)
    # @show r y₀ e
    # return r, m, y₀, logy₀, e
    # logkern_5(r) - logy₀ + e
    # logm = ifelse(y₀isone, log1pr, log1pr - logy₀)
    # log1pr += e

    # logm = ifelse(y₀notone, log1pr - logy₀, log1pr)
    # @show logm log1pr logy₀ e
    # logm + e
end
@inline vlog_fast(x::T) where {T} = vlog2_fast(x) * convert(T, 0.6931471805599453)
@inline vlog10_fast(x::T) where {T} = vlog2_fast(x) * convert(T, 0.3010299956639812)
# @inline Base.FastMath.log_fast(v::AbstractSIMD) = vlog_fast(float(v))
# @inline Base.FastMath.log2_fast(v::AbstractSIMD) = vlog2_fast(float(v))
# @inline Base.FastMath.log10_fast(v::AbstractSIMD) = vlog10_fast(float(v))

@inline function log2_kern_5_256(x)
    c5 = 1.442695040888963430242018860033667730529246266516984575206349078995811595149434
    c4 = -0.7213475204444818238116304451613321925305725926194819917643836452178866375548279
    c3 = 0.480898346935995048804761640936874918715608085587068841980211938177173076051027
    c2 = -0.3606737601723789903264415189400867418476075783130968675245051405508304609911771
    c1 = 0.2885437254815682617764304037319501313717679206671517967342968352103864725323379
    c0 = -0.2404546770228689226321452773936263107841762894627528631584290286000941540854746
    x * vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0,x,c1),x,c2),x,c3),x,c4),x,c5)
end

const LOG2_TABLE_128 = Float64[log2(x) for x ∈  range(big"0.5", step = 1/256, length = 128)]

@inline function vlog2_v2(a::AbstractSIMD{W,Float64}) where {W}
    # log2(vgetmant(a, Val(8))) + vgetexp(a) == log2(a)
    m = vgetmant(a, Val(8)) # m ∈ [1,2)
    e = vgetexp(a)
    # m ∈ [1,2); We reduce further...
    y = inv_approx(m) # y ∈ (0.5,1]
    # return y
    y₀ = vroundscale(y, Val(128)) # 80 = 16*(1+7)

    # y₀ is a multiple of 32; 0.5:
    # r + 1 = m*y;
    # log(r+1) = log(m) + log(y₀)
    # log(m) = log(1+r) - log(y₀)
    r = vfmsub(m, y₀, 1.0)
    # log1pr = logkern_5(r)
    log1pr = log2_kern_5_256(r)
    y₀notone = y₀ ≠ 1.0
    inds = (reinterpret(Int, y₀) >> 45) & 0xff
    # inds = (reinterpret(UInt, y₀) >>> 48) & 0x0f
    # @show y₀ inds r
    # logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
    logy₀ = vload(zstridedpointer(LOG2_TABLE_128), (inds,))

    emlogy₀ = ifelse(y₀notone, e - logy₀, e)
    log1pr + emlogy₀
end

# @inline function Base.log(x1::AbstractSIMD{W,Float64}) where {W}
# @inline function vlog(x1::Float64)
# @inline Base.log(v::AbstractSIMD) = log(float(v))

# @inline function vlog_fast(x1::AbstractSIMD{W,Float32}) where {W}
#     notzero = x1 != zero(x1)
#     greater_than_zero = x1 > zero(x1)
#     x3 = x1 < 1.1754944f-38#3.4332275f-5
#     # x6 = true if x3 entirely false
#     x7 = x1 * 8.388608f6#14.0f0
#     x8 = ifelse(x3, x7, x1)
#     x10 = reinterpret(UInt32, x8)
#     x11 = x10 + 0x004afb0d
#     x12 = x11 >>> 0x00000017
#     x13 = ifelse(x3, 0xffffff6a, 0xffffff81)
#     x15 = x12 + x13
#     x16 = x11 & 0x007fffff
#     x17 = x16 + 0x3f3504f3
#     x18 = reinterpret(Float32, x17)
#     x19 = x18 - 1f0
#     x20 = x18 + 1f0
#     x21 = x19 / x20
#     x22 = x21 * x21
#     x23 = x22 * x22
#     x24 = vfmadd(x23, 0.24279079f0, 0.40000972f0)
#     x25 = x24 * x23
#     x26 = vfmadd(x23, 0.6666666f0, 0.6666666f0)
#     x27 = x26 * x22
#     x28 = x19 * x19
#     x29 = x28 * 5f-1
#     x30 = convert(Float32, x15 % Int32)
#     x31 = x27 + x29
#     x32 = x31 + x25
#     x33 = x30 * -0.00021219444f0
#     x34 = vfmadd(x21, x32, x33)
#     x35 = x19 - x29
#     x36 = x35 + x34
#     x37 = vfmadd(x30, 0.6933594f0, x36)
#     # x37 = vfmadd(x30, 0.6931472f0, x36)
#     x39 = ifelse(x1 == Inf32, Inf32, x37)
#     x40 = ifelse(notzero, x39, -Inf32)
#     ifelse(x1 < zero(x1), NaN32, x40)
# end
# @inline Base.FastMath.log_fast(v::AbstractSIMD) = vlog_fast(float(v))


text:
using Distributions

# ρ(x) in integrand
ρ(x) = 1/sqrt(π)*exp(-x^2)

#=
Implementation of random walk sampler Metropolis-Hastings.
Inputs:
- n: The number of states to return
- b: The standard deviation of the normal distribution used in the proposal
     distribution
Outputs:
- Returns an array of length n denoting states S_i in Metropolis-Hastings alg.
=#
function metropolis_hastings(n, b)

    # Samples
    samples = zeros(n)
    samples[1] = rand()

    # Conditional Q distribution
    d = Normal(0,b)

    acceptances = 0

    for i=2:n
        # Propose Y
        ξ = rand(d)
        Y = samples[i-1] + ξ

        # Acceptance probability
        α = min(1, ρ(Y)/ρ(samples[i-1]))

        # Next state
        if rand() < α
            samples[i] = Y
            acceptances += 1
        else
            samples[i] = samples[i-1]
        end
    end

    println(acceptances)

    return samples
end

# Obtain 1000 samples with b = 2
S = metropolis_hastings(10000,1)

# f(x) in integrand
f(x) = -1/2 < x < 2 ? abs(cos(x)) : 0

# Estimate integral
int_estimate = sqrt(π)*mean([f(x) for x in S])

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
## cholesky.jl : Gaussian random field generator based on the Cholesky decomposition of the covariance matrix

"""
    Cholesky()

Retuns a [`GaussianRandomFieldGenerator`](@ref) based on a Cholesky factorization of the covariance matrix.

# Examples
```jldoctest
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, Cholesky(), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a Cholesky decomposition

julia> heatmap(grf)
[...]

```

!!! warning

    The Cholesky factorization requires the covariance matrix to be symmetric and positive definite. If the covariance matrix is not positive definite, an error will be thrown. Try using the `Spectral` method in that case.

See also: [`Spectral`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
struct Cholesky <: GaussianRandomFieldGenerator end

# compose a GaussianRandomField using Cholesky factorization
function _GaussianRandomField(mean, cov::CovarianceFunction, method::Cholesky, pts...)
    C = apply(cov, pts...)

    # compute Cholesky factorization
    isposdef!(Symmetric(C)) || throw(ArgumentError("to use a Cholesky factorization, the covariance matrix must be positive definite"))
    L = LowerTriangular(C')

    GaussianRandomField{Cholesky,typeof(cov),typeof(pts),typeof(mean),typeof(L)}(mean,cov,pts,L)
end

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{Cholesky}) = size(grf.data, 1)

# sample from the GaussianRandomField using Cholesky factorization
function _sample(grf::GaussianRandomField{Cholesky}, xi)
	grf.mean + std(grf.cov) * reshape(grf.data * xi, size(grf.mean))
end

Base.show(io::IO, ::Cholesky) = print(io, "Cholesky decomposition")


text:
"""
bijective function {1, ..., k} -> {1, ..., k}
"""
struct Perm{T <: Integer} <: BasicDiscreteFunction
    m :: Vector{T}
    dom :: ShiftedSegment{Int}
    rng :: ShiftedSegment{Int}
end

"""
bijective function {0, .. k-1} -> {0, ..., k-1}
"""
struct OffsetPerm{T <: Integer} <: BasicDiscreteFunction
    m :: Vector{T}
    dom :: Segment{Int}
    rng :: Segment{Int}
end

function Perm(v)
    k = maximum(v)
    dom = Segment(1, k)
    return Perm(v, dom, dom)
end

function OffsetPerm(v)
    k = maximum(v)
    dom = Segment(k)
    return OffsetPerm(v, dom, dom)
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions