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

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

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

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

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

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

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

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

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

# jo_full(jo)
jo_full(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = A*jo_eye(DDT,A.n)

# norm(jo)
norm(A::joAbstractLinearOperator,p::Real=2) = norm(elements(A),p)

# real(jo)
#real{DDT<:Real,RDT<:Real}(A::joAbstractLinearOperator{DDT,RDT}) = A
function real(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT}
    throw(joAbstractLinearOperatorException("real(jo) not implemented"))
end
joReal(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = real(A)

# imag(jo)
#imag{DDT<:Real,RDT<:Real}(A::joAbstractLinearOperator{DDT,RDT}) = joZeros(A.m,A.n,DDT,RDT)
function imag(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT}
    throw(joAbstractLinearOperatorException("imag(jo) not implemented"))
end
joImag(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} = imag(A)

# conj(jo)
joConj(A::joAbstractLinearOperator) = conj(A)
conj(A::joLinearOperator{DDT,RDT}) where {DDT,RDT} =
    joLinearOperator{DDT,RDT}("conj("*A.name*")",A.m,A.n,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop
        )
conj(A::joMatrix{DDT,RDT}) where {DDT,RDT} =
    joMatrix{DDT,RDT}("conj("*A.name*")",A.m,A.n,
        A.fop_C, A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop
        )
conj(A::joLinearFunction{DDT,RDT}) where {DDT,RDT} =
    joLinearFunction{DDT,RDT}("conj("*A.name*")",A.m,A.n,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop, A.fMVok,
        A.iop_C, A.iop_A, A.iop_T, A.iop, A.iMVok
        )

# transpose(jo)
transpose(A::joLinearOperator{DDT,RDT}) where {DDT,RDT} =
    joLinearOperator{RDT,DDT}("transpose("*A.name*")",A.n,A.m,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A
        )
transpose(A::joMatrix{DDT,RDT}) where {DDT,RDT} =
    joMatrix{RDT,DDT}("transpose("*A.name*")",A.n,A.m,
        A.fop_T, A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A
        )
transpose(A::joLinearFunction{DDT,RDT}) where {DDT,RDT} =
    joLinearFunction{RDT,DDT}("transpose("*A.name*")",A.n,A.m,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A, A.fMVok,
        A.iop_T, A.iop, A.iop_C, A.iop_A, A.iMVok
        )

# adjoint(jo)
adjoint(A::joLinearOperator{DDT,RDT}) where {DDT,RDT} =
    joLinearOperator{RDT,DDT}("adjoint("*A.name*")",A.n,A.m,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T
        )
adjoint(A::joMatrix{DDT,RDT}) where {DDT,RDT} =
    joMatrix{RDT,DDT}("adjoint("*A.name*")",A.n,A.m,
        A.fop_A, A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T
        )
adjoint(A::joLinearFunction{DDT,RDT}) where {DDT,RDT} =
    joLinearFunction{RDT,DDT}("adjoint("*A.name*")",A.n,A.m,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T, A.fMVok,
        A.iop_A, A.iop_C, A.iop, A.iop_T, A.iMVok
        )

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

# issymmetric(jo)
issymmetric(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} =
    (A.m == A.n && (norm(elements(A)-elements(transpose(A))) < joTol))

# ishermitian(jo)
ishermitian(A::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT} =
    (A.m == A.n && (norm(elements(A)-elements(adjoint(A))) < joTol))

# getindex(jo,...)
getindex(A::joAbstractLinearOperator{DDT,RDT},::Colon,i::Integer) where {DDT,RDT} =
    joHelpers_etc.elements_column_helper(A,i)
getindex(A::joAbstractLinearOperator{DDT,RDT},::Colon,r::UnitRange{URT}) where {DDT,RDT,URT<:Integer}  =
    joHelpers_etc.elements_column_helper(A,r)

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

# *(jo,jo)
function *(A::joAbstractLinearOperator{CDT,ARDT},B::joAbstractLinearOperator{BDDT,CDT}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joAbstractLinearOperatorException("shape 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
# in joAbstractSAparallelLinearOperator/base_functions.jl
#function *(A::joSAdistributingLinearOperator{CDT,ARDT,2},B::joAbstractLinearOperator{BDDT,CDT}) where {ARDT,BDDT,CDT}
#function *(A::joAbstractLinearOperator{CDT,ARDT},B::joSAgatheringLinearOperator{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
# in joAbstractDAparallelLinearOperator/base_functions.jl
#function *(A::joDAdistributingLinearOperator{CDT,ARDT,2},B::joAbstractLinearOperator{BDDT,CDT}) where {ARDT,BDDT,CDT}
#function *(A::joAbstractLinearOperator{CDT,ARDT},B::joDAgatheringLinearOperator{BDDT,CDT,2}) where {ARDT,BDDT,CDT}

# *(jo,mvec)
function *(A::joLinearOperator{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joLinearOperatorException("shape 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::joMatrix{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joMatrixException("shape 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::joLinearFunction{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joLinearFunctionException("shape mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    if A.fMVok
        MV=A.fop(mv)
        jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    else
        MV=Matrix{ARDT}(undef,A.m,size(mv,2))
        for i=1:size(mv,2)
            V=A.fop(mv[:,i])
            i==1 && jo_check_type_match(ARDT,eltype(V),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    end
    return MV
end

# *(mvec,jo)

# *(jo,vec)
function *(A::joLinearOperator{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.n == size(v,1) || throw(joLinearOperatorException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(jo,vec):",A.name,typeof(A),vDT]," / "))
    V=A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    return V
end
function *(A::joMatrix{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.n == size(v,1) || throw(joMatrixException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(jo,vec):",A.name,typeof(A),vDT]," / "))
    V=A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    return V
end
function *(A::joLinearFunction{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.n == size(v,1) || throw(joLinearFunctionException("shape mismatch"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(jo,vec):",A.name,typeof(A),vDT]," / "))
    V=A.fop(v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(vec,jo)

# *(num,jo)
function *(a::Number,A::joAbstractLinearOperator{ADDT,ARDT}) where {ADDT,ARDT}
    return joLinearOperator{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,
        v1->jo_convert(ARDT,a*(A*v1),false),
        v2->jo_convert(ADDT,a*(transpose(A)*v2),false),
        v3->jo_convert(ADDT,conj(a)*(adjoint(A)*v3),false),
        v4->jo_convert(ARDT,conj(a)*(conj(A)*v4),false),
        @joNF, @joNF, @joNF, @joNF
        )
end
function *(a::joNumber{ADDT,ARDT},A::joAbstractLinearOperator{ADDT,ARDT}) where {ADDT,ARDT}
    return joLinearOperator{ADDT,ARDT}("(N*"*A.name*")",A.m,A.n,
        v1->jo_convert(ARDT,a.rdt*(A*v1),false),
        v2->jo_convert(ADDT,a.ddt*(transpose(A)*v2),false),
        v3->jo_convert(ADDT,conj(a.ddt)*(adjoint(A)*v3),false),
        v4->jo_convert(ARDT,conj(a.rdt)*(conj(A)*v4),false),
        @joNF, @joNF, @joNF, @joNF
        )
end

# *(jo,num)
*(A::joAbstractLinearOperator{ADDT,ARDT},a::Number) where {ADDT,ARDT} = a*A
*(A::joAbstractLinearOperator{ADDT,ARDT},a::joNumber{ADDT,ARDT}) where {ADDT,ARDT} = a*A

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

# \(jo,jo)

# \(jo,mvec)
function \(A::joLinearOperator{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.m == size(mv,1) || throw(joLinearOperatorException("shape mismatch"))
    jo_check_type_match(ARDT,mvDT,join(["RDT for \\(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    if hasinverse(A)
        MV=get(A.iop)(mv)
        jo_check_type_match(ADDT,eltype(MV),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    elseif issquare(A)
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4square(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    elseif (istall(A) && !isnull(jo_iterative_solver4tall))
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4tall(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    elseif (iswide(A) && !isnull(jo_iterative_solver4wide))
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4wide(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    else
        throw(joLinearOperatorException("\\(jo,MultiVector) not supplied"))
    end
    return MV
end
function \(A::joMatrix{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.m == size(mv,1) || throw(joMatrixException("shape mismatch"))
    jo_check_type_match(ARDT,mvDT,join(["RDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    if hasinverse(A)
        MV=get(A.iop)(mv)
        jo_check_type_match(ADDT,eltype(MV),join(["DDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    else
        throw(joMatrixException("\\(jo,Vector) not supplied"))
    end
    return MV
end
function \(A::joLinearFunction{ADDT,ARDT},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.m == size(mv,1) || throw(joLinearFunctionException("shape mismatch"))
    jo_check_type_match(ARDT,mvDT,join(["RDT for \\(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    if hasinverse(A)
        if A.iMVok
            MV = get(A.iop)(mv)
        else
            MV=Matrix{ADDT}(undef,A.n,size(mv,2))
            for i=1:size(mv,2)
                V=get(A.iop)(mv[:,i])
                i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
                MV[:,i]=V
            end
        end
    elseif issquare(A)
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4square(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    elseif (istall(A) && !isnull(jo_iterative_solver4tall))
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4tall(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    elseif (iswide(A) && !isnull(jo_iterative_solver4wide))
        MV=Matrix{ADDT}(undef,A.n,size(mv,2))
        for i=1:size(mv,2)
            V=jo_convert(ADDT,jo_iterative_solver4wide(A,mv[:,i]))
            i==1 && jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,mvec):",A.name,typeof(A),eltype(V)]," / "))
            MV[:,i]=V
        end
    else
        throw(joLinearFunctionException("\\(jo,MultiVector) not supplied"))
    end
    return MV
end

# \(mvec,jo)

# \(jo,vec)
function \(A::joLinearOperator{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.m == size(v,1) || throw(joLinearOperatorException("shape mismatch"))
    jo_check_type_match(ARDT,vDT,join(["RDT for \\(jo,vec):",A.name,typeof(A),vDT]," / "))
    if hasinverse(A)
        V=get(A.iop)(v)
        jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    elseif issquare(A)
        V=jo_convert(ADDT,jo_iterative_solver4square(A,v))
    elseif (istall(A) && !isnull(jo_iterative_solver4tall))
        V=jo_convert(ADDT,jo_iterative_solver4tall(A,v))
    elseif (iswide(A) && !isnull(jo_iterative_solver4wide))
        V=jo_convert(ADDT,jo_iterative_solver4wide(A,v))
    else
        throw(joLinearOperatorException("\\(jo,Vector) not supplied"))
    end
    return V
end
function \(A::joMatrix{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.m == size(v,1) || throw(joMatrixException("shape mismatch"))
    jo_check_type_match(ARDT,vDT,join(["RDT for *(jo,vec):",A.name,typeof(A),vDT]," / "))
    if hasinverse(A)
        V=get(A.iop)(v)
        jo_check_type_match(ADDT,eltype(V),join(["DDT from *(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    else
        throw(joMatrixException("\\(jo,Vector) not supplied"))
    end
    return V
end
function \(A::joLinearFunction{ADDT,ARDT},v::LocalVector{vDT}) where {ADDT,ARDT,vDT<:Number}
    A.m == size(v,1) || throw(joLinearFunctionException("shape mismatch"))
    jo_check_type_match(ARDT,vDT,join(["RDT for \\(jo,vec):",A.name,typeof(A),vDT]," / "))
    if hasinverse(A)
        V=get(A.iop)(v)
        jo_check_type_match(ADDT,eltype(V),join(["DDT from \\(jo,vec):",A.name,typeof(A),eltype(V)]," / "))
    elseif issquare(A)
        V=jo_convert(ADDT,jo_iterative_solver4square(A,v))
    elseif (istall(A) && !isnull(jo_iterative_solver4tall))
        V=jo_convert(ADDT,jo_iterative_solver4tall(A,v))
    elseif (iswide(A) && !isnull(jo_iterative_solver4wide))
        V=jo_convert(ADDT,jo_iterative_solver4wide(A,v))
    else
        throw(joLinearFunctionException("\\(jo,Vector) not supplied"))
    end
    return V
end

# \(vec,jo)

# \(num,jo)

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

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

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

# +(jo,jo)
function +(A::joAbstractLinearOperator{DDT,RDT},B::joAbstractLinearOperator{DDT,RDT}) where {DDT,RDT}
    size(A) == size(B) || throw(joAbstractLinearOperatorException("shape mismatch"))
    return joLinearOperator{DDT,RDT}("("*A.name*"+"*B.name*")",size(A,1),size(B,2),
        v1->A*v1+B*v1,
        v2->transpose(A)*v2+transpose(B)*v2,
        v3->adjoint(A)*v3+adjoint(B)*v3,
        v4->conj(A)*v4+conj(B)*v4,
        @joNF, @joNF, @joNF, @joNF
        )
end

# +(jo,mvec)

# +(mvec,jo)

# +(jo,vec)

# +(vec,jo)

# +(jo,num)
function +(A::joAbstractLinearOperator{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return joLinearOperator{ADDT,ARDT}("("*A.name*"+N)",size(A,1),size(A,2),
        v1->A*v1+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
        v2->transpose(A)*v2+joConstants(A.n,A.m,b;DDT=ARDT,RDT=ADDT)*v2,
        v3->adjoint(A)*v3+joConstants(A.n,A.m,conj(b);DDT=ARDT,RDT=ADDT)*v3,
        v4->conj(A)*v4+joConstants(A.m,A.n,conj(b);DDT=ADDT,RDT=ARDT)*v4,
        @joNF, @joNF, @joNF, @joNF
        )
end
function +(A::joAbstractLinearOperator{ADDT,ARDT},b::joNumber{ADDT,ARDT}) where {ADDT,ARDT}
    return joLinearOperator{ADDT,ARDT}("("*A.name*"+N)",size(A,1),size(A,2),
        v1->A*v1+joConstants(A.m,A.n,b.rdt;DDT=ADDT,RDT=ARDT)*v1,
        v2->transpose(A)*v2+joConstants(A.n,A.m,b.ddt;DDT=ARDT,RDT=ADDT)*v2,
        v3->adjoint(A)*v3+joConstants(A.n,A.m,conj(b.ddt);DDT=ARDT,RDT=ADDT)*v3,
        v4->conj(A)*v4+joConstants(A.m,A.n,conj(b.rdt);DDT=ADDT,RDT=ARDT)*v4,
        @joNF, @joNF, @joNF, @joNF
        )
end

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

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

# -(jo)
-(A::joLinearOperator{DDT,RDT}) where {DDT,RDT} =
    joLinearOperator{DDT,RDT}("(-"*A.name*")",A.m,A.n,
        v1->-A.fop(v1),
        v2->-get(A.fop_T)(v2),
        v3->-get(A.fop_A)(v3),
        v4->-get(A.fop_C)(v4),
        v5->-get(A.iop)(v5),
        v6->-get(A.iop_T)(v6),
        v7->-get(A.iop_A)(v7),
        v8->-get(A.iop_C)(v8)
        )
-(A::joMatrix{DDT,RDT}) where {DDT,RDT} =
    joMatrix{DDT,RDT}("(-"*A.name*")",A.m,A.n,
        v1->-A.fop(v1),
        v2->-A.fop_T(v2),
        v3->-A.fop_A(v3),
        v4->-A.fop_C(v4),
        v5->-get(A.iop)(v5),
        v6->-get(A.iop_T)(v6),
        v7->-get(A.iop_A)(v7),
        v8->-get(A.iop_C)(v8)
        )
-(A::joLinearFunction{DDT,RDT}) where {DDT,RDT} =
    joLinearFunction{DDT,RDT}("(-"*A.name*")",A.m,A.n,
        v1->-A.fop(v1),
        v2->-get(A.fop_T)(v2),
        v3->-get(A.fop_A)(v3),
        v4->-get(A.fop_C)(v4),
        A.fMVok,
        v5->-get(A.iop)(v5),
        v6->-get(A.iop_T)(v6),
        v7->-get(A.iop_A)(v7),
        v8->-get(A.iop_C)(v8),
        A.iMVok
        )

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

# -(jo,mvec)

# -(mvec,jo)

# -(jo,vec)

# -(vec,jo)

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

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

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

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

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

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

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

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

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

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

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

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

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

# hcat(...jo...)
hcat(ops::joAbstractLinearOperator...) = joDict(ops...)

# vcat(...jo...)
vcat(ops::joAbstractLinearOperator...) = joStack(ops...)

# hvcat(...jo...)
hvcat(rows::Tuple{Vararg{Int}}, ops::joAbstractLinearOperator...) = joBlock(collect(rows),ops...)

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

# mul!(...,jo,...)
mul!(y::LocalVector{RDT},A::joAbstractLinearOperator{DDT,RDT},x::LocalVector{DDT}) where {DDT,RDT} = y[:] = A * x
mul!(y::LocalMatrix{RDT},A::joAbstractLinearOperator{DDT,RDT},x::LocalMatrix{DDT}) where {DDT,RDT} = y[:,:] = A * x

# ldiv!(...,jo,...)
ldiv!(y::LocalVector{DDT},A::joAbstractLinearOperator{DDT,RDT},x::LocalVector{RDT}) where {DDT,RDT} = y[:] = A \ x
ldiv!(y::LocalMatrix{DDT},A::joAbstractLinearOperator{DDT,RDT},x::LocalMatrix{RDT}) where {DDT,RDT} = y[:,:] = A \ x


text:
using MoodleQuestions
using Documenter

makedocs(;
    modules=[MoodleQuestions],
    authors="Daniel Molina <dmolina@decsai.ugr.es>",
    repo="https://github.com/dmolina/MoodleQuestions.jl/blob/{commit}{path}#L{line}",
    sitename="MoodleQuestions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dmolina.github.io/MoodleQuestions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Instructions" => "formato.md",
        "Server Mode" => "server.md"
    ],
)

deploydocs(;
    repo="github.com/dmolina/MoodleQuestions.jl",
)

ext: jl

max_stars_repo_name: dmolina/MoodleQuestions

max_issues_repo_name: dmolina/MoodleQuestions

max_forks_repo_name: dmolina/MoodleQuestions


text:
using Nystrom
using StaticArrays
using LinearAlgebra
using IterativeSolvers
using Plots

##
k = 3.3
pde = Maxwell(;dim=3,k)

# exact solution
G    = SingleLayerKernel(pde)
xs   = SVector(0.1,0.2,0.3)
x0   = SVector(5,5,5)
pol    = SVector(1+im,-2.,3.)
E    = (dof) -> G(dof,xs)*pol
exa  = E(x0)

# geometry and parameters
sph_radius = 2
geo = ParametricSurfaces.Sphere(;radius=sph_radius)
Ω = Geometry.Domain(geo)
Γ = boundary(Ω)

## CFIE
α = 1
β = im*k/sph_radius
order = 5         # quadrature order 1D
ndofs = Float64[]
errs = Float64[]
iterative = true;
#compress = Matrix
compress = Nystrom.IFGFCompressor(;p=(3,5,5),nmax=100);
##
for n in [4]
    M     = meshgen(Γ,(n,n))
    nmesh = NystromMesh(view(M,Γ);order)
    J     = Nystrom.diagonal_jacobian_matrix(nmesh)
    dualJ = Nystrom.diagonal_dualjacobian_matrix(nmesh)
    γ₀E   = γ₀(E,nmesh)
    rhs   = dualJ*ncross(γ₀E)

    @info "Computing operators..."
    T,K = Nystrom.maxwell_dim(pde,nmesh;compress)   # EFIE, MFIE
    @info "Assembling system matrix..."
    L   = Nystrom.assemble_maxwell_indirect_cfie(nmesh, α, β, K, T)
    @info "Solving..."
    if iterative
        ϕ_coeff = gmres(L,rhs; verbose=true,maxiter=600,restart=600,abstol=1e-6)
    else
        ϕ_coeff = L\rhs
    end
    ϕ     = J*ϕ_coeff
    Spot  = Nystrom.MaxwellEFIEPotential(pde,nmesh)
    Dpot  = Nystrom.MaxwellMFIEPotential(pde,nmesh)
    Eₐ    = (x) -> α*Dpot[ϕ](x) + β*Spot[ncross(ϕ)](x)
    er    = (Eₐ(x0) - exa)/norm(exa)

    ndof = length(Nystrom.dofs(nmesh))
    err = norm(er,Inf)

    push!(ndofs, ndof)
    push!(errs, err)
    @info "" ndof err
end

## Plot
sqrt_ndofs = sqrt.(ndofs)
fig = plot(sqrt_ndofs,errs,xscale=:log10,yscale=:log10,m=:o,label="error",lc=:black)
title = "k=$k, α=$α, β=$β, qorder=$order"
plot!(xlabel="√ndofs",ylabel="error",title=title)
for p in 1:5
    cc = errs[end]*sqrt_ndofs[end]^p
    plot!(fig,sqrt_ndofs,cc./sqrt_ndofs.^p,xscale=:log10,yscale=:log10,label="h^$p",ls=:dash)
end
display(fig)

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
"""
    infercontent(filename|io) -> (formats[], codings[])

Guess possible formats/codings from a filename or IO stream.

In the case of a filename, the extension is checked; for IO streams, the first
few bytes of the stream are examined.
"""
function infercontent(filename::AbstractString)
	format_guesses = []
	coding_guesses = []
	filename2, extension = splitext(filename)
	extension = lowercase(extension)
	for mime in get(registry.extensions, extension, [])
		if mime in registry.codings
			push!(coding_guesses, mime)
		elseif mime in registry.formats
			push!(format_guesses, mime)
		end
	end
	if length(coding_guesses) > 0 && length(format_guesses) > 0
		coding_guesses = []
	end
	if length(coding_guesses) > 0
		extension2 = lowercase(splitext(filename2)[2])
		for mime in get(registry.extensions, extension2, [])
			if mime in registry.formats
				push!(format_guesses, mime)
			end
		end
	end
	format_guesses, coding_guesses
end

function infercontent(io::IO)
	format_guesses = []
	coding_guesses = []
	mark(io)
	beginning = read(io, 512)
	reset(io)
	if length(beginning) > 0
		for (signature, mimes) in registry.signatures
			if length(beginning) >= length(signature)
				if view(beginning, axes(signature)...) == signature
					for mime in mimes
						if mime in registry.codings
							push!(coding_guesses, mime)
						elseif mime in registry.formats
							push!(format_guesses, mime)
						end
					end
				end
			end
		end
	end
	if length(format_guesses) > 0 && length(coding_guesses) > 0
		coding_guesses = []
	end
	format_guesses, coding_guesses
end


text:
using BenchmarkTools
using BasicBSpline
using GeometryBasics
using Random

const SUITE = BenchmarkGroup()

Random.seed!(42)

p = 3
k = Knots(rand(100)) + (p+1) * Knots(0,1)
P = FastBSplineSpace(p,k)
i = 1
t = 0.02

SUITE["fast-bsplinebasis"] = @benchmarkable bsplinebasis(i,P,t)

p1 = 2
k1 = Knots(rand(3)) + (p1 + 1) * Knots([0, 1])
P1 = FastBSplineSpace(p1, k1)
n1 = dim(P1)
a = [Point(i1, rand()) for i1 in 1:n1]
M = BSplineCurve([P1], a)

p1′ = p1 + 1
k1′ = k1 + unique(k1) + Knots(rand(2))
P1′ = FastBSplineSpace(p1′, k1′)
M′ = refinement(M, [P1′])

SUITE["refinement"] = @benchmarkable refinement(M, [P1′])
SUITE["fitting_I"] = @benchmarkable fittingcontrolpoints(M, [P1′])
SUITE["fitting_R"] = @benchmarkable fittingcontrolpoints(M, [P1′], domain=:R)


text:
function POMDPModelTools.gbmdp_handle_terminal(pomdp::Union{SubHuntPOMDP,DSubHuntPOMDP}, updater::Updater, b::ParticleCollection, s::SubState, a::Int, rng::AbstractRNG)
    @assert isterminal(pomdp, s)
    return ParticleCollection([s,s])
end


text:
using Base.Cartesian: @nexprs

macro multistep(fname,rule)
	@assert rule.head == :(=)
	y = rule.args[1].args[1]  # output variable
    vecs = Set{Symbol}()
    consts = Set{Symbol}()
    inds = Set{Symbol}()
    offsets = Set{Int}()
    extractvars!(rule,vecs,consts,inds,offsets)
    delete!(vecs,y)
    vecs = sort([vecs...])
    consts = sort([consts...])
    @assert length(inds) == 1 "only 1 index variable allowed"
    ind = pop!(inds)
    maxoffset = maximum(offsets)
    minoffset = minimum(offsets)
    nconsts = maxoffset - minoffset
    fname! = symbol(string(fname) * "!")

    yarg = :($y::AbstractVector{T})
    yconsts = [symbol(string(y) * "_$i") for i in 1:nconsts]
    yconstargs = [:($c::T) for c in yconsts]
    refargs = [Expr(:(::), v, :(AbstractVector{T})) for v in vecs]
    constargs = [Expr(:(::), c, :T) for c in consts]
    args = [yconstargs;refargs;constargs]
    fsig = esc(Expr(:call, Expr(:curly, fname, :T), args...))
    fsig! = esc(Expr(:call, Expr(:curly, fname!, :T), yarg, args...))
    fcall = esc(Expr(:call, Expr(:curly, fname!, :T), y, yconsts..., vecs..., consts...))
    ind = esc(ind)

    quote
        $fsig! = begin
            @nexprs $nconsts j -> $y[j] = $yconsts[j]
            @inbounds for $ind in $(1-minoffset):length($y)-$maxoffset
                $rule
            end
        end

        $fsig = begin
            $(esc(y)) = similar($(esc(vecs[1])))
            $(Expr(:escape, Expr(:call, fname!, y, yconsts..., vecs..., consts...)))
            $(esc(y))
        end
    end
end

function extractvars!(e, vecs, consts, inds, offsets)
    if isa(e, Symbol) && isalpha(string(e))
        push!(consts, e)
        return
    end

    if isa(e, Expr)
        if e.head == :ref
            push!(vecs, e.args[1])
            ind = e.args[2]

            if isa(ind, Symbol)
                push!(inds, ind)
            elseif isa(ind, Expr)
                op = ind.args[1]
                @assert ind.head == :call && (op == :+ || op == :-)
                offset = op == :+ ? ind.args[3] : -ind.args[3]
                push!(inds, ind.args[2])
                push!(offsets, offset)
            end

            return
        end

        for i in 1:length(e.args)
            extractvars!(e.args[i], vecs, consts, inds, offsets)
        end
    end
end

macro adamsmoulton3_multistep(name,kernel)
    vecs = Set()
    consts = Set()
    inds = Set()
    offsets = Set()
    extractvars!(kernel,vecs,consts,inds,offsets)
    delete!(vecs,:y)
    delete!(vecs,:f)
    push!(consts,:h)
    vecs = sort([vecs...])
    consts = sort([consts...])
    @assert length(inds) == 1 "only 1 index variable allowed"
    ind = pop!(inds)
    yarg = :(y::AbstractVector)
    vecargs = [Expr(:(::), v, :(AbstractVector)) for v in vecs]
    fsig = esc(Expr(:call, :($name), yarg, vecargs..., consts...))

	y = esc(:y)
	i = esc(ind)
	f = esc(:f)

    quote
        $fsig = begin
            $f = similar($y)

            $y[1] = 0.
            $f[1] = $(replace(kernel,ind,1))
            $y[2] = $y[1] + h*$f[1]
            $f[2] = $(replace(kernel,ind,2))
            $y[3] = $y[2] + h*($f[2]*3./2. - $f[1]/2.)
            $f[3] = $(replace(kernel,ind,3))
            for $i in 4:length($y)
                $y[$i] = $y[$i-1] + h*($f[$i-1]*23./12. - $f[$i-2]*4./3. + $f[$i-3]*5./12.)
                $f[$i] = $(replace(kernel,ind,i))
            end
        end
    end
end

function replace!(e,var,val)
    !isa(e,Expr) && return
    for i in 1:length(e.args)
        if isa(e.args[i],Expr)
            replace!(e.args[i],var,val)
        elseif e.args[i] == var
            e.args[i] = val
        end
    end
end

function replace(e,var,val)
    ec = copy(e)
    replace!(ec,var,val)
    ec
end

function isvariable(s)
    if isa(s, Symbol) && isalpha(string(s))
        return true
    end

    if isa(s, Expr)
        return reduce(|, map(isvariable, s.args))
    end

    return false
end

max_stars_repo_name: johncsnyder/Element


text:
"""
    search_server(data_config_path, io_port, search_server_ready; cache_path=nothing)

Search server for Garamond. It is a finite-state-machine that
when called, creates the searchers i.e. search objects using the
`data_config_path` and the proceeds to looping continuously
in order to asynchronously handle outside requests.

After the searchers are loaded, the search server sends a notification
using `search_server_ready` to any listening I/O servers.
"""
function search_server(data_config_path, io_port, search_server_ready; cache_path=nothing)

    # Build search environment
    env = build_search_env(data_config_path; cache_path=cache_path)

    # Start the search environment operator
    up_in_channel = Channel{String}(0)          # input: dictionary with keys the command and its argument
    up_out_channel = Channel{typeof(env)}(0)    # output: updated environment
    channels = (up_in_channel, up_out_channel)
    Threads.@spawn env_operator(env, channels)

    # Notify waiting I/O servers
    @info "• Notifying I/O servers..."
    notify(search_server_ready, true)

    # Start search server
    ipaddr = Sockets.localhost
    server = listen(ipaddr, io_port)
    @info "SEARCH server online @$ipaddr:$io_port..."

    # Main loop: accept connection, update (if updates available) and respond
    counter = [0]  # vector so the value is mutable
    while true
        sock = accept(server)
        if isready(up_out_channel)
            env = take!(up_out_channel)
        end
        @async respond(env, sock, counter, channels)
    end
end

"""
    respond(env, socket, counter, channels)

Responds to search server requests received on `socket` using
the search data from `searchers`. The requests are counted
through the variable `counter`.
"""
function respond(env, socket, counter, channels)
    # Channels for updating
    up_in_channel, up_out_channel = channels

    # Read and parse an outside request
    request = parse(InternalRequest, readline(socket))
    counter.+= 1
    @debug "* Received [#$(counter[1])]: $request."

    timer_start = time()
    if request.operation === :search
        ### Search ###
        search_results = search(env, request)
        ranked_search_results = rank(env, request, search_results)
        etime = time() - timer_start
        response = build_response(env.dbdata,
                                  request,
                                  ranked_search_results;
                                  id_key=env.id_key,
                                  elapsed_time=etime)
        write(socket, response * RESPONSE_TERMINATOR)
        @info "• Search [#$(counter[1])]: query='$(request.query)' completed in $etime(s)."

    elseif request.operation === :recommend
        ### Recommend ###
        recommendations = recommend(env, request)
        ranked_recommendations = rank(env, request, recommendations)
        etime = time() - timer_start
        response = build_response(env.dbdata,
                                  request,
                                  ranked_recommendations;
                                  id_key=env.id_key,
                                  elapsed_time=etime)
        write(socket, response * RESPONSE_TERMINATOR)
        @info "• Recommendation [#$(counter[1])]: completed in $etime(s)."

    elseif request.operation === :rank
        ranked_ids = rank(env, request, nothing)  #::Vector{SearchResult}
        etime = time() - timer_start
        response = build_response(env.dbdata, request, ranked_ids; id_key=env.id_key, elapsed_time=etime)
        @info "• Rank [#$(counter[1])]: completed in $etime(s)."
        write(socket, response * RESPONSE_TERMINATOR)

    elseif request.operation === :kill
        ### Kill the search server ###
        @info "• Kill [#$(counter[1])]: Exiting in 1(s)..."
        write(socket, RESPONSE_TERMINATOR)
        sleep(1)
        exit()

    elseif request.operation === :read_configs
        ### Read and return data configurations ###
        @info "• Get configuration [#$(counter[1])]."
        write(socket, read_configuration_to_json(env) * RESPONSE_TERMINATOR)

    elseif request.operation === :envop
        ### Read and return data configurations ###
        @info "• Environment operation [#$(counter[1])]."
        if !isready(up_in_channel)
            @async put!(up_in_channel, request.query)  # the take! is in the search server
        else
            @warn "Environment operation request ignored, another in progress."
        end
        write(socket, RESPONSE_TERMINATOR)  # send response asynchronously

    elseif request.operation === :error
        @info "• Errored request [#$(counter[1])]: Ignoring..."
        write(socket, RESPONSE_TERMINATOR)

    else
        @info "• Unknown request [#$(counter[1])]: Ignoring..."
        write(socket, RESPONSE_TERMINATOR)
    end

    close(socket)
    @debug "* Response [#$(counter[1])]: done after $(time()-timer_start)(s)."
end

"""
    build_response(dbdata, request, results, [; kwargs...])

Builds a response for an engine client using the data, request and results.
"""
function build_response(dbdata,
                        request,
                        results;
                        id_key=DEFAULT_DB_ID_KEY,
                        elapsed_time=-1.0)
    if !isempty(results)
        n_total_results = mapreduce(x->length(x.query_matches), +, results)
    else
        n_total_results = 0
    end

    response_results = Dict{String, Vector{Dict{Symbol, Any}}}()
    return_fields = Tuple(intersect!(unique!([request.return_fields..., id_key]),
                                     colnames(dbdata)))
    for result in results
        dict_vector = []
        offset = (request.response_page - 1) * request.response_size
        scores, indices = unzip(result.query_matches; n=request.response_size, offset=offset, ndims=2)
        dataresult = sort(rows(filter(in(indices), dbdata, select=id_key), return_fields),
                          by=row->getproperty(row, id_key))
        for (entry, score) in sort(collect(zip(dataresult, scores[sortperm(indices)])),
                                   by=zipped->zipped[2], rev=true)
            dict_entry = Dict{Symbol, Any}(pairs(entry))
            push!(dict_entry, :score => score)
            push!(dict_vector, dict_entry)
        end
        push!(response_results, result.id => dict_vector)
    end
    response = Dict("elapsed_time"=>elapsed_time,
                    "results" => response_results,
                    "n_total_results" => n_total_results,
                    "n_searchers" => length(results),
                    "n_searchers_w_results" => mapreduce(r->!isempty(r.query_matches), +, results),
                    "suggestions" => squash_suggestions(results, request.max_suggestions))
    return JSON.json(response)
 end


text:
using JuLIP, QMMM2, SHIPs, PrettyTables
using JuLIP.MLIPs
using NeighbourLists, LinearAlgebra

# fname = @__DIR__() * "/C_test_data.json"
fname = "/Users/ortner/Dropbox/Work/Projects/SiteEFitting/data/C_Dim3.json"
D = load_json(fname)

trans = PolyTransform(2, rnn(:C))
weights = Dict("Es" => 0.01, "dEs" => 10.0, "d2Esh" => 1.0)

err0 = Float64[]
err1 = Float64[]
err2 = Float64[]

@info("4B BASIS TEST")
for deg in [6, 8, 10, 12, 14, 16 ] #, 18, 20]
   basis = IPSuperBasis(OneBody(1.0),
                        SHIPBasis(3, deg, 1.0, trans, 2, 0.7*rnn(:C), 5.0) )
   @info("deg = $deg, len(B) = $(length(basis))")
   🚢, fitinfo = QMMM2.lsqfit(basis, D, weights; rtol=1e-5)
   @show fitinfo["maxe"]
   push!(err0, fitinfo["maxe"]["Es"])
   push!(err1, fitinfo["maxe"]["dEs"])
   push!(err2, fitinfo["maxe"]["d2Esh"])
end

@info("Summary of Errors:")
pretty_table([err0 err1 err2], ["Es", "dEs", "d2Es"],
             formatter=ft_printf("%.3e"))

# @info("5B BASIS TEST")
# for deg in [6, 8, 10, 12, 13, 14 ]
#    basis = IPSuperBasis(OneBody(1.0),
#                         SHIPBasis(4, deg, 1.5, trans, 2, 0.7*rnn(:C), 5.0) )
#    @info("deg = $deg, len(B) = $(length(basis))")
#    🚢, fitinfo = QMMM2.lsqfit(basis, D, weights; rtol=1e-5)
#    @show fitinfo["maxe"]
#    push!(err0, fitinfo["maxe"]["Es"])
#    push!(err1, fitinfo["maxe"]["dEs"])
#    push!(err2, fitinfo["maxe"]["d2Esh"])
# end
#
# @info("Summary of Errors:")
# pretty_table([err0 err1 err2], ["Es", "dEs", "d2Es"],
#              formatter=ft_printf("%.3e"))


text:
using DataDeps

register(DataDep(
  "seismic-bumps",
  "https://archive.ics.uci.edu/ml/datasets/seismic-bumps",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/00266/seismic-bumps.arff",
  "aabe512fab65b36d1dfb462650b75cfd8d99d8cc2723e8ecb4e6f5e1caccd5a7",
  post_fetch_method=(path -> begin
    # skipto causes an error here for some reason, just delete the lines
    lines = readlines(path)
    write(path, join(lines[155:end], "\n"))

    UCIData.process_dataset(path,
        target_index=19,
        feature_indices=1:18,
        categoric_indices=[1:3; 8],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using Test

using SIMD
using ChipSort

@testset "Sorting networks" begin

function test_sorting_networks(T, N)
    vv = rand(T, N)
    srtref = sort(vv)
    srthat = sort_net(vv...)

    @test all(srthat[2:end] .>= srthat[1:end-1])
    @test all(srthat .== srtref)
end

@testset for T in [Int8, Int16, Int32, Int64, Float32, Float64]
    for L in 2 .^ (2:5)
        test_sorting_networks(T, L)
    end
end

function test_sort_simd_vec(T, N)
    va = ntuple(n->Vec(tuple(rand(T, N)...)),N)
    srthat = sort_net(va...)

    xx = [va[j][k] for j in 1:N, k in 1:N]
    srtref = sort(xx; dims=1)
    srthat_arr = [srthat[j][k] for j in 1:N, k in 1:N]

    @test srthat_arr == srtref
end

@testset for T in [Int8, Int16, Int32, Int64, Float32, Float64]
    Niterations = 16
    for L in 2 .^ (2:5)
        for it in 1:Niterations
            test_sort_simd_vec(T, L)
        end
    end
end

@testset for N in [16, 32]
    for it in 1:2^16
        test_sort_simd_vec(UInt32, N)
    end
end

function test_sort_heaviside(N, Na)
    va = ntuple(n->if (Na&(1<<(n-1))) > 0 1 else 0 end, N)
    srthat = sort_net(va...)
    srtref = sort([va...])
    @test srthat == tuple(srtref...)
end

@testset for N in 2 .^ (2:4)
    for Na in 0:(2^N-1)
            test_sort_heaviside(N, Na)
    end
end

end


text:
# checks Fresnel's transmission coefficient numerically
n_incident = 1.
n_transmit = 1.45
ref_local_field = get_local_field(Float64[0.5], refractive_indexes=ones(3) * n_incident)[3]
interface_localfield = get_local_field(Float64[0.5],  refractive_indexes=Float64[n_transmit, n_incident, n_incident])[3]

transmission_numerical = abs2(interface_localfield)/abs2(ref_local_field)/n_incident*n_transmit
transmission_fresnel = ((2n_incident)/(n_incident+n_transmit))^2/n_incident*n_transmit

@test abs(transmission_fresnel - transmission_numerical)/abs(transmission_fresnel)< 1e-2


text:
"""
    AbstractCumulant
"""
abstract type AbstractCumulant end

function get(cumulant::AbstractCumulant, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("get(CumulantType, args...) not defined!"))
end

"""
    FeatureCumulant <: AbstractCumulant

Basic Cumulant which takes the value c_tp1 = s_tp1[idx] for 1<=idx<=length(s_tp1)
"""
struct FeatureCumulant <: AbstractCumulant
    idx::Int
end

Base.get(cumulant::FeatureCumulant, state_tp1, action_tp1, preds_tp1) =
    state_tp1[cumulant.idx]

"""
    PredictionCumulant <: AbstractCumulant

Basic cumulant which takes the value c_tp1 = preds_tp1[idx]
"""
struct PredictionCumulant <: AbstractCumulant
    idx::Int
end

Base.get(cumulant::PredictionCumulant, state_tp1, action_tp1, preds_tp1) =
    preds_tp1[cumulant.idx]

"""
    ScaledCumulant{F<:Number, T<:AbstractCumulant} <: AbstractCumulant

A cumulant which scales another AbstractCumulant
"""
struct ScaledCumulant{F<:Number, T<:AbstractCumulant} <: AbstractCumulant
    scale::F
    cumulant::T
end

Base.get(cumulant::ScaledCumulant, state_tp1, action_tp1, preds_tp1) =
    cumulant.scale*get(cumulant.cumulant, state_tp1, action_tp1, preds_tp1)

"""
    FunctionalCumulant{F<:Function} <: AbstractCumulant

A cumulant that has a user defined function c_tp1 = f(state_tp1, action_tp1, preds_tp1)
"""
struct FunctionalCumulant{F<:Function} <: AbstractCumulant
    f::F
end

Base.get(cumulant::FunctionalCumulant, state_tp1, action_tp1, preds_tp1) =
    cumulant.f(state_tp1, action_tp1, preds_tp1)

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
mutable struct CompoundXMLElement <: MyXMLElement
    els::Vector{MyXMLElement}
end

function make_xml(cp::CompoundXMLElement)
    for xml in cp.els
        make_xml(xml)
    end
end


text:
using Printf: @sprintf

import Base.copy
import LogicCircuits.IO: SDDElement,
    PSDDElement,
    save_lines,
    parse_psdd_file,
    PsddHeaderLine,
    LcHeaderLine,
    save_sdd_file,
    save_as_dot,
    get_nodes_level

# Saving psdd

#####################
# decompile for nodes
#####################

# decompile for psdd
decompile(n::ProbLiteral, node2id, vtree2id)::UnweightedLiteralLine =
    UnweightedLiteralLine(node2id[n], vtree2id[n.origin.vtree], literal(n), true)

make_element(n::Prob⋀, w::AbstractFloat, node2id) =
    PSDDElement(node2id[n.children[1]],  node2id[n.children[2]], w)

is_true_node(n)::Bool =
    GateType(n) isa ⋁Gate && num_children(n) == 2 && GateType(children(n)[1]) isa LiteralGate && GateType(children(n)[2]) isa LiteralGate &&
    positive(children(n)[1]) && negative(children(n)[2])

function decompile(n::Prob⋁, node2id, vtree2id)::Union{WeightedNamedConstantLine, DecisionLine{PSDDElement}}
    if is_true_node(n)
        WeightedNamedConstantLine(node2id[n], vtree2id[n.origin.vtree], lit2var(n.children[1].origin.literal), n.log_thetas[1]) # TODO
    else
        DecisionLine(node2id[n], vtree2id[n.origin.vtree], UInt32(num_children(n)), map(x -> make_element(x[1], x[2], node2id), zip(children(n), n.log_thetas)))
    end
end

#####################
# build maping
#####################

function get_node2id(ln::AbstractVector{X}, T::Type)where X #<: T#::Dict{T, ID}
    node2id = Dict{T, ID}()
    outnodes = filter(n -> !(GateType(n) isa ⋀Gate), ln)
    sizehint!(node2id, length(outnodes))
    index = ID(0) # node id start from 0
    for n in outnodes
        node2id[n] = index
        index += ID(1)
    end
    node2id
end

function get_vtree2id(ln::PlainVtree):: Dict{PlainVtreeNode, ID}
    vtree2id = Dict{PlainVtreeNode, ID}()
    sizehint!(vtree2id, length(ln))
    index = ID(0) # vtree id start from 0

    for n in ln
        vtree2id[n] = index
        index += ID(1)
    end
    vtree2id
end

#####################
# saver for circuits
#####################

function psdd_header()
    """
    c ids of psdd nodes start at 0
    c psdd nodes appear bottom-up, children before parents
    c
    c file syntax:
    c psdd count-of-sdd-nodes
    c L id-of-literal-sdd-node id-of-vtree literal
    c T id-of-trueNode-sdd-node id-of-vtree variable log(litProb)
    c D id-of-decomposition-sdd-node id-of-vtree number-of-elements {id-of-prime id-of-sub log(elementProb)}*
    c
    c File generated by Juice.jl
    c"""
end

function save_psdd_file(name::String, ln::ProbΔ, vtree::PlainVtree)
    # TODO add method isstructured
    @assert ln[end].origin isa StructLogicalΔNode "PSDD should decorate on StructLogicalΔ"
    @assert endswith(name, ".psdd")
    node2id = get_node2id(ln, ProbΔNode)
    vtree2id = get_vtree2id(vtree)
    formatlines = Vector{CircuitFormatLine}()
    append!(formatlines, parse_psdd_file(IOBuffer(psdd_header())))
    push!(formatlines, PsddHeaderLine(num_nodes(ln)))
    for n in filter(n -> !(GateType(n) isa ⋀Gate), ln)
        push!(formatlines, decompile(n, node2id, vtree2id))
    end
    save_lines(name, formatlines)
end

function lc_header()
    """
    c variables (from inputs) start from 1
    c ids of logistic circuit nodes start from 0
    c nodes appear bottom-up, children before parents
    c the last line of the file records the bias parameter
    c three types of nodes:
    c	T (terminal nodes that correspond to true literals)
    c	F (terminal nodes that correspond to false literals)
    c	D (OR gates)
    c
    c file syntax:
    c Logistic Circuit
    c T id-of-true-literal-node id-of-vtree variable parameters
    c F id-of-false-literal-node id-of-vtree variable parameters
    c D id-of-or-gate id-of-vtree number-of-elements (id-of-prime id-of-sub parameters)s
    c B bias-parameters
    c
    c File generated by Juice.jl
    c"""
end

function save_lc_file(name::String, ln::LogisticΔ, vtree)
    @assert ln[end].origin isa StructLogicalΔNode "LC should decorate on StructLogicalΔ"
    @assert endswith(name, ".circuit")
    node2id = get_node2id(ln, ProbΔNode)
    vtree2id = get_vtree2id(vtree)
    formatlines = Vector{CircuitFormatLine}()
    append!(formatlines, parse_lc_file(IOBuffer(lc_header())))
    push!(formatlines, LcHeaderLine())
    for n in filter(n -> !(GateType(n) isa ⋀Gate), ln)
        push!(formatlines, decompile(n, node2id, vtree2id))
    end
    save_lines(name, formatlines)
end

import LogicCircuits.save_circuit # make available for extension

function save_circuit(name::String, circuit, vtree=nothing)
    if endswith(name, ".circuit")
        save_lc_file(name, circuit, vtree)
    elseif endswith(name, ".psdd")
        save_psdd_file(name, circuit, vtree)
    elseif endswith(name, ".sdd")
        save_sdd_file(name, circuit, vtree)
    else
        error("Cannot save circuit to file with this extensions: $name")
    end
end

"Save prob circuit to .dot file"
function save_as_dot(root::ProbΔNode, file::String)
    return save_as_dot(node2dag(root), file)
end

"Save prob circuits to .dot file"
function save_as_dot(circuit::ProbΔ, file::String)
    # TODO (https://github.com/Juice-jl/LogicCircuits.jl/issues/7)
    node_cache = Dict{ProbΔNode, Int64}()
    for (i, n) in enumerate(circuit)
        node_cache[n] = i
    end

    levels = get_nodes_level(circuit)

    f = open(file, "w")
    write(f,"digraph Circuit {\nsplines=false\nedge[arrowhead=\"none\",fontsize=6]\n")

    for level in levels
        if length(level) > 1
            write(f,"{rank=\"same\";newrank=\"true\";rankdir=\"LR\";")
            rank = ""
            foreach(x->rank*="$(node_cache[x])->",level)
            rank = rank[1:end-2]
            write(f, rank)
            write(f,"[style=invis]}\n")
        end
    end

    for n in reverse(circuit)
        if n isa Prob⋀
            write(f, "$(node_cache[n]) [label=\"*$(node_cache[n])\"]\n")
        elseif n isa Prob⋁
            write(f, "$(node_cache[n]) [label=\"+$(node_cache[n])\"]\n")
        elseif n isa ProbLiteral && positive(n)
            write(f, "$(node_cache[n]) [label=\"+$(variable(n.origin))\"]\n")
        elseif n isa ProbLiteral && negative(n)
            write(f, "$(node_cache[n]) [label=\"-$(variable(n.origin))\"]\n")
        else
            throw("unknown ProbNode type")
        end
    end

    for n in reverse(circuit)
        if n isa Prob⋀
            for c in n.children
                write(f, "$(node_cache[n]) -> $(node_cache[c])\n")
            end
        elseif n isa Prob⋁
            for (c, p) in zip(n.children, exp.(n.log_thetas))
                prob = @sprintf "%0.1f" p
                write(f, "$(node_cache[n]) -> $(node_cache[c]) [label=\"$prob\"]\n")
            end
        else
        end
    end

    write(f, "}\n")
    flush(f)
    close(f)
end


text:
using JuDoc, Test, Markdown
const J = JuDoc
const D = joinpath(dirname(dirname(pathof(JuDoc))), "test", "_dummies")

# NOTE this first file MUST be included before running the rest of the tests
# otherwise you may get an error like "key 0x099191234..." was not found or
# saying that the key :in doesn't exist or something along those lines
include("jd_paths_vars.jl"); include("test_utils.jl")

include("misc.jl")

# MANAGER folder
include("manager/utils.jl")
println("🍺")

# PARSER folder
println("PARSER/MD+LX")
include("parser/markdown+latex.jl")
include("parser/markdown-extra.jl")
include("parser/footnotes.jl")
println("🍺")
println("PARSER/HTML")
include("parser/html.jl")
println("🍺")

# CONVERTER folder
println("CONVERTER/MD")
include("converter/markdown.jl")
include("converter/markdown2.jl")
include("converter/markdown3.jl")
include("converter/hyperref.jl")
println("🍺")
println("CONVERTER/HTML")
include("converter/html.jl")
println("🍺")

println("CONVERTER/LX")
include("converter/eval.jl")
include("converter/lx_input.jl")
include("converter/lx_simple.jl")
println("🍺")

println("INTEGRATION")
include("global/cases1.jl")
include("global/cases2.jl")
include("global/ordering.jl")

begin
    # create temp dir to do complete integration testing (has to be here in order
    # to locally play nice with node variables etc, otherwise it's a big headache)
    p = joinpath(D, "..", "__tmp");
    # after errors, this may not have been deleted properly
    isdir(p) && rm(p; recursive=true, force=true)
    # make dir, go in it, do the tests, then get completely out (otherwise windows
    # can't delete the folder)
    mkdir(p); cd(p);
    include("global/postprocess.jl");
    cd(joinpath(D, ".."))
    # clean up
    rm(p; recursive=true, force=true)
end
cd(dirname(dirname(pathof(JuDoc))))
println("😅 😅 😅 😅")


text:
struct PushRecoveryController{M<:MomentumBasedController, S, E<:SE3PDController, L}
    lowlevel::M
    robotmass::Float64
    gravitymag::Float64

    end_effector_tasks::Vector{SpatialAccelerationTask}
    linmomtask::LinearMomentumRateTask
    pelvistask::AngularAccelerationTask
    jointtasks::Dict{JointID, JointAccelerationTask{Revolute{Float64}}}

    statemachine::S
    end_effector_controllers::Vector{E}
    linear_momentum_controller::L
    pelvisgains::PDGains{Diagonal{Float64, SVector{3, Float64}}, Diagonal{Float64, SVector{3, Float64}}}
    jointgains::Dict{JointID, PDGains{Float64, Float64}}

    comref::Point3D{SVector{3, Float64}}
    jointrefs::Dict{JointID, Float64}

    active_contact_points::Dict{BodyID, Vector{SPoint3D{Float64}}}
end

function PushRecoveryController(
        lowlevel::MomentumBasedController,
        pelvis::RigidBody,
        nominalstate::MechanismState,
        statemachine,
        end_effector_controllers::Vector{<:SE3PDController},
        linear_momentum_controller,
        joint_regularization::Float64 = 0.05,
        linear_momentum_weight::Float64 = 1.0,
        pelvisgains::PDGains = default_pelvis_gains(),
        jointgains = Dict(JointID(j) => critically_damped_gains(100.0) for j in tree_joints(lowlevel.state.mechanism)),
        comref::Point3D = center_of_mass(nominalstate) - FreeVector3D(root_frame(lowlevel.state.mechanism), 0., 0., 0.05),
        jointrefs = Dict(JointID(j) => configuration(nominalstate, j)[1] for j in tree_joints(lowlevel.state.mechanism) if joint_type(j) isa Revolute))
    mechanism = lowlevel.state.mechanism
    world = root_body(mechanism)
    worldframe = root_frame(mechanism)
    m = mass(mechanism)
    g = norm(mechanism.gravitational_acceleration)

    regularize!.(Ref(lowlevel), tree_joints(mechanism), joint_regularization)

    end_effector_control_joints = Set{Joint}()
    end_effector_tasks = SpatialAccelerationTask[]
    for controller in end_effector_controllers
        base = findbody(mechanism, controller.base)
        body = findbody(mechanism, controller.body)
        task = SpatialAccelerationTask(mechanism, path(mechanism, base, body))
        push!(end_effector_tasks, task)
        for joint in task.path
            push!(end_effector_control_joints, joint)
        end
        addtask!(lowlevel, task, Parameter(lowlevel.qpmodel, val=controller.weight))
    end

    linmomtask = LinearMomentumRateTask(mechanism, centroidal_frame(lowlevel))
    addtask!(lowlevel, linmomtask, linear_momentum_weight)

    pelvistask = AngularAccelerationTask(mechanism, path(mechanism, world, pelvis))
    addtask!(lowlevel, pelvistask, 10.0)

    revolutejoints = filter(j -> joint_type(j) isa Revolute, tree_joints(mechanism))
    positioncontroljoints = setdiff(revolutejoints, end_effector_control_joints)
    jointtasks = Dict(JointID(j) => JointAccelerationTask(j) for j in positioncontroljoints)
    addtask!.(Ref(lowlevel), collect(values(jointtasks)), 10.0)

    active_contact_points = Dict{BodyID, Vector{SPoint3D{Float64}}}()

    PushRecoveryController(
        lowlevel, m, g,
        end_effector_tasks, linmomtask, pelvistask, jointtasks,
        statemachine, end_effector_controllers, linear_momentum_controller, pelvisgains, jointgains,
        comref, jointrefs,
        active_contact_points)
end

function (controller::PushRecoveryController)(τ::AbstractVector, t::Number, state::MechanismState)
    # Call state machine
    controller.statemachine(t, state)

    # Update active contact points
    update_active_contacts!(controller.active_contact_points, controller.lowlevel.contacts)

    # End effector control
    for i in eachindex(controller.end_effector_tasks)
        setdesired!(controller.end_effector_tasks[i], controller.end_effector_controllers[i](t, state))
    end

    # Linear momentum control
    c = center_of_mass(state)
    h = momentum(state)
    ċ = FreeVector3D(h.frame, linear(h) / controller.robotmass)
    l̇_des = controller.linear_momentum_controller(t, c, ċ, controller.active_contact_points, state)
    world_to_centroidal = Transform3D(c.frame, centroidal_frame(controller.lowlevel), -c.v)
    l̇_des = transform(l̇_des, world_to_centroidal)
    setdesired!(controller.linmomtask, l̇_des)

    # Pelvis orientation control
    pelvis = target(controller.pelvistask.path)
    Hpelvis = transform_to_root(state, pelvis)
    Tpelvis = transform(twist_wrt_world(state, pelvis), inv(Hpelvis))
    ωd_pelvis_des = FreeVector3D(Tpelvis.frame, pd(controller.pelvisgains, rotation(Hpelvis), Tpelvis.angular))
    setdesired!(controller.pelvistask, ωd_pelvis_des)

    # Joint position control
    for jointid in keys(controller.jointtasks)
        task = controller.jointtasks[jointid]
        gains = controller.jointgains[jointid]
        ref = controller.jointrefs[jointid]
        v̇_joint_des = pd(gains, configuration(state, jointid)[1], ref, velocity(state, jointid)[1], 0.0)
        setdesired!(task, v̇_joint_des)
    end

    controller.lowlevel(τ, t, state)
    τ
end

function update_active_contacts!(active_contact_points::AbstractDict{BodyID, <:AbstractVector{<:Point3D}}, contacts)
    for active_points in values(active_contact_points)
        empty!(active_points)
    end
    for (body, points) in contacts
        bodyid = BodyID(body)
        active_points_body = get!(valtype(active_contact_points), active_contact_points, bodyid)
        for contactpoint in points
            if QPControl.isenabled(contactpoint)
                push!(active_points_body, contactpoint.position)
            end
        end
    end
    nothing
end

function default_pelvis_gains()
    gains_xy = critically_damped_gains(100.0)
    gains_z = critically_damped_gains(50.0)
    PDGains(Diagonal(SVector(gains_xy.k, gains_xy.k, gains_z.k)), Diagonal(SVector(gains_xy.d, gains_xy.d, gains_z.d)))
end


text:
# ------------ pthread.h ------------
const libpthread = "libpthread"
# Justification: https://refspecs.linuxfoundation.org/LSB_3.1.0/LSB-Core-generic/LSB-Core-generic/libpthread-ddefs.html
const pthread_t = Culong

const __ULONG32_TYPE = Culong
const __t_uscalar_t = __ULONG32_TYPE
const __ULONGWORD_TYPE = __t_uscalar_t
const __CPU_MASK_TYPE = __ULONGWORD_TYPE
const __cpu_mask = __CPU_MASK_TYPE

const __CPU_SETSIZE = 1024
const __NCPUBITS = (8 * sizeof(__cpu_mask))

"""
Data structure to describe CPU mask.

Ref: [docs](https://github.com/lattera/glibc/blob/master/posix/bits/cpu-set.h)
"""
struct cpu_set_t
    __bits::NTuple{16,__cpu_mask}
end

"""
Returns the ID of the calling thread.
This is the same value that is returned in `*thread` in the
`pthread_create(3)`` call that created this thread.

Ref: [docs](https://man7.org/linux/man-pages/man3/pthread_self.3.html)
"""
pthread_self() = @ccall libpthread.pthread_self()::pthread_t

"""
Ref: [docs](https://man7.org/linux/man-pages/man3/pthread_setaffinity_np.3.html)
"""
function pthread_setaffinity_np(thread, cpussetsize, cpuset)
    @ccall libpthread.pthread_setaffinity_np(
        thread::pthread_t, cpussetsize::Csize_t, cpuset::Ptr{cpu_set_t}
    )::Cint
end

# function pinthread(processorId::Integer)
#     cpuset = cpu_set_t()
#     cpuset_ref = Ref(cpuset)
#     thread = zero(pthread_t)

#     thread = pthread_self()
#     # CPU_ZERO(&cpuset)
#     # CPU_SET(processorId, &cpuset);
#     pthread_setaffinity_np(thread, sizeof(cpu_set_t), cpuset_ref)
# end


text:
@doc raw"""
    βdifferential_geodesic_startpoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_x g(t;p,q)[X]$.
They are

```math
\beta(\kappa) = \begin{cases}
\frac{\sinh(d(1-t)\sqrt{-\kappa})}{\sinh(d\sqrt{-\kappa})}
&\text{ if }\kappa < 0,\\
1-t & \text{ if } \kappa = 0,\\
\frac{\sin((1-t)d\sqrt{\kappa})}{\sinh(d\sqrt{\kappa})}
&\text{ if }\kappa > 0.
\end{cases}
```

Due to a symmetry agrument, these are also used to compute $D_q g(t; p,q)[\eta]$

# See also

[`differential_geodesic_endpoint`](@ref), [`differential_geodesic_startpoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_geodesic_startpoint(κ,t,d)
    if (d==0) || (κ==0)
        return (1-t)
    else
        if κ < 0
            return sinh(sqrt(-κ)*(1-t)*d)/sinh(sqrt(-κ)*d)
        elseif κ > 0
            return sin(sqrt(κ)*(1-t)*d)/sin(sqrt(κ)*d)
        end
    end
end
@doc raw"""
    βdifferential_exp_basepoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_p \exp_p X [Y]$. They are

```math
\beta(\kappa) = \begin{cases}
\cosh(\sqrt{-\kappa})&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\cos(\sqrt{\kappa}) &\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_exp_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_exp_basepoint(κ,t,d)
    if κ < 0
        return cosh(sqrt(-κ)*d)
    elseif κ > 0
        return cos(sqrt(κ)*d)
    else
        return 1.0;
    end
end
@doc raw"""
    βdifferential_exp_argument(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_X \exp_p X[Y]$. They are

```math
$\beta(\kappa) = \begin{cases}
\frac{\sinh(d\sqrt{-\kappa})}{d\sqrt{-\kappa}}&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\frac{\sin(d\sqrt{\kappa})}{\sqrt{d\kappa}}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_exp_argument`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_exp_argument(κ,t,d)
    if (d==0) || (κ==0)
        return 1.0
    else
        if κ < 0
            return sinh(sqrt(-κ)*d)/( d*sqrt((-κ)) )
        elseif κ > 0
            return sin( sqrt(κ)*d )/( d*sqrt(κ) )
        end
    end
end
@doc raw"""
    βdifferential_log_basepoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_p \log_p q[X]$. They are

```math
\beta(\kappa) = \begin{cases}
-\sqrt{-\kappa}d\frac{\cosh(d\sqrt{-\kappa})}{\sinh(d\sqrt{-\kappa})}&\text{ if }\kappa < 0,\\
-1 & \text{ if } \kappa = 0,\\
-\sqrt{\kappa}d\frac{\cos(d\sqrt{\kappa})}{\sin(d\sqrt{\kappa})}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_log_argument`](@ref), [`differential_log_argument`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_log_basepoint(κ::Number,t::Number,d::Number)
    if (d==0) || (κ==0)
        return -1.0
    else
        if κ < 0
            return - sqrt(-κ)*d*cosh(sqrt(-κ)*d)/sinh(sqrt(-κ)*d)
        else #if κ > 0
            return - sqrt(κ)*d*cos(sqrt(κ)*d)/sin(sqrt(κ)*d)
        end
    end
end
@doc raw"""
    βdifferential_log_argument(κ,t,d)

weights for the JacobiField corresponding to the differential of the logarithmic
map with respect to its argument $D_q \log_p q[X]$. They are

```math
\beta(\kappa) = \begin{cases}
\frac{ d\sqrt{-\kappa} }{\sinh(d\sqrt{-\kappa})}&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\frac{ d\sqrt{\kappa} }{\sin(d\sqrt{\kappa})}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_log_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_log_argument(κ::Number,t::Number,d::Number)
    if (d==0) || (κ==0)
        return 1.0
    else
        if κ < 0
            return sqrt(-κ)*d/sinh(sqrt(-κ)*d)
        else #if κ > 0
            return sqrt(κ)*d/sin(sqrt(κ)*d)
        end
    end
end

@doc raw"""
    Y = adjoint_Jacobi_field(M, p, q, t, X, β)

Compute the AdjointJacobiField $J$ along the geodesic $γ_{p,q}$ on the manifold
$\mathcal M$ with initial conditions (depending on the application)
$X ∈ T_{γ_{p,q}(t)}\mathcal M$ and weights $β$. The result is a vector
$Y ∈ T_p\mathcal M$. The main difference to [`jacobi_field`](@ref) is the,
that the input `X` and the output `Y` switched tangent spaces.
For detais see [`jacobi_field`](@ref)
"""
function adjoint_Jacobi_field(M::Manifold, p, q, t, X, β::Function=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M,p,q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M),Ref(p), V, Ref(x), Ref(ParallelTransport()))
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (
            inner.(Ref(M),Ref(x),Ref(X),Θ)
        ) .* (
            β.(B.data.eigenvalues,Ref(t),Ref(distance(M,p,q)))
        ) .* V
    )
    return Y
end
function adjoint_Jacobi_field(M::Circle{ℝ}, p::Real, q::Real, t::Real, X::Real, β::Function=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M,p,q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M),Ref(p), V, Ref(x), Ref(ParallelTransport()))
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (
            inner.(Ref(M),Ref(x),Ref(X),Θ)
        ) .* (
            β.(B.data.eigenvalues,Ref(t),Ref(distance(M,p,q)))
        ) .* V
    )[1]
    return Y
end
function adjoint_Jacobi_field(M::AbstractPowerManifold, p, q, t, X, β::Function=βdifferential_geodesic_startpoint)
    rep_size = representation_size(M.manifold)
    Y = allocate_result(M, adjoint_Jacobi_field, p, X)
    for i in get_iterator(M)
        #lY = adjoint_Jacobi_field(M.manifold, p[i], q[i], t, X[i], β)
        lY = adjoint_Jacobi_field(
            M.manifold,
            get_component(M,p,i),
            get_component(M,q,i),
            t,
            get_component(M,X,i),
            β,
        )
        set_component!(M,Y,lY,i)
    end
    return Y
end

@doc raw"""
    Y = jacobi_field(M, p, q, t, X, β)
compute the Jacobi jield $J$ along the geodesic $γ_{p,q}$ on the manifold $\mathcal M$ with
initial conditions (depending on the application) $X ∈ T_p\mathcal M$ and weights $β$. The
result is a tangent vector `Y` from $T_{γ_{p,q}(t)}\mathcal M$.

# See also

[`adjoint_Jacobi_field`](@ref)
"""
function jacobi_field(M::Manifold, p, q, t, X, β::Function=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t);
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M,p,q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.( Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    Y = zero_tangent_vector(M,p)
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y .= sum(
        (
            inner.(Ref(M), Ref(p), Ref(X), V)
        ) .* (
            β.(B.data.eigenvalues,Ref(t),Ref(distance(M,p,q)))
        ) .* Θ
    )
    return Y
end
function jacobi_field(M::Circle{ℝ}, p::Real, q::Real, t::Real, X::Real, β::Function=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t);
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M,p,q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.( Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    Y = zero_tangent_vector(M,p)
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (
            inner.(Ref(M), Ref(p), Ref(X), V)
        ) .* (
            β.(B.data.eigenvalues,Ref(t),Ref(distance(M,p,q)))
        ) .* Θ
    )[1]
    return Y
end
function jacobi_field(M::AbstractPowerManifold, p, q, t, X, β::Function=βdifferential_geodesic_startpoint)
    rep_size = representation_size(M.manifold)
    Y = allocate_result(M, adjoint_Jacobi_field, p, X)
    for i in get_iterator(M)
        #lY = jacobi_field(M.manifold, p[i], q[i], t, X[i], β)
        set_component!(
            M,
            Y,
            jacobi_field(
                M.manifold,
                get_component(M,p,i),
                get_component(M,q,i),
                t,
                get_component(M,X,i),
                β,
            ),
            i,
        )
    end
    return Y
end


text:
# Hermitian lattices
@attributes mutable struct HermLat{S, T, U, V, W} <: AbsLat{S}
  space::HermSpace{S, T, U, W}
  pmat::V
  gram::U
  rational_span::HermSpace{S, T, U, W}
  base_algebra::S
  involution::W
  automorphism_group_generators::Vector{U}
  automorphism_group_order::fmpz
  generators
  minimal_generators
  norm
  scale

  function HermLat{S, T, U, V, W}() where {S, T, U, V, W}
    z = new{S, T, U, V, W}()
    return z
  end

  function HermLat(E::S, G::U, P::V) where {S, U, V}
    involutionL = involution(E)

    K = base_field(E)

    space = hermitian_space(E, G)

    z = new{S, typeof(K), U, V, typeof(involutionL)}(space, P)
    z.base_algebra = E
    z.involution = involutionL
    return z
  end

  function HermLat(E::S, G::U) where {S, U}
    n = nrows(G)
    M = pseudo_matrix(identity_matrix(E, n))
    return HermLat(E, G, M)
  end
end


text:
# Author: Mikael Koli
# All rights reserved

module BotMarket

    # Functions and structs for the module
    export AbstractInvestor, AbstractMarket,
        # Structs
        DoubleAuctionMarket,

        ZeroIntelligentInvestor,
        ExternalWorld,
        AskLimitOrder, BidLimitOrder,
        RandomWalkAsset, FixedAsset,

        # Functions
        place!,
        clear!,
        cancel_all!,
        release_all!,
        update!,
        pay_cashflows!,

        # Generics
        generic_currency,
        generic_stock,

        # For debugging
        get_order,
        get_order_quantity,
        get_order_price,
        get_side,
        get_unreserved,
        # Abstracts
        AbstractAsset,
        Trade

    # Including the module files where the
    # the above structs and functions come
    # from

    include("abstracts.jl")

    include("market/orders.jl")
    include("market/trade.jl")
    include("market/clear.jl")
    include("market/limitorderbook.jl")

    include("investor/api.jl")
    include("investor/generic.jl")
    include("investor/position.jl")

    include("asset/api.jl")
    include("external_world/external_world.jl")

    # Built in
    include("builtin/zero_intelligent.jl")
    include("builtin/asset.jl")
    include("builtin/double_auction.jl")
    include("builtin/currency.jl")

    # Defining some generics for easier use
    generic_currency = RandomWalkAsset(0, 0, 0, 0, "ccy")
    generic_stock = RandomWalkAsset(0, 0, 0, 0, "stock")

end


text:
using Pkg

Pkg.add("Publish")
Pkg.add("GraphPlot")

using Publish
using FluxTraining

deploy(FluxTraining; root = "/FluxTraining.jl", label="dev", force = true)


text:
"""
    mass_hu(ρ, I, V, S_BG, S_O)
    mass_hu(vol, ρ, S_BG, S_OBJ)

Calculates the mass of a known material given a CT volume segment.

# Arguments
- vol: segmented CT volume
- ρ: density of known material
- I: sum of total volume segment
- V: total volume of segment
    if voxels are 1 cm^3, then `V = length(vol)`
- S_BG: intensity of known background voxels
- S_O: intensity of known object voxels

# Citation
'Accurate quantification of vessel cross-sectional area using CT
angiography: a simulation study' [Molloi, Johnson, Ding, Lipinski]
(DOI: 10.1007/s10554-016-1007-9)

"""
function mass_hu(vol, ρ, S_BG, S_O)
	I = sum(vol)
	V = length(vol)
	mₒ = ρ * ((I - (V * S_BG)) / (S_O - S_BG))
	return mₒ
end

function mass_hu(ρ, I, V, S_BG, S_O)
	mₒ = ρ * ((I - (V * S_BG)) / (S_O - S_BG))
	return mₒ
end


text:
using HTTP
import HTTP:bytes
using JSON
using JSON2
"""
Routingを行う場合
"""
const localhost = "0.0.0.0"
const port = parse(Int, ENV["PORT"])
println("ポート番号: ", port)
const headers = ["Content-Type" => "application/json"]

struct TestStruct
    name::String
    number::Int
end

# apitest/nameのname部分とTestStructを紐づけ，uniqueなnameを持つ
test_structs = Dict{String, TestStruct}()

"""
辞書型をJSON文字列に変換する
"""
function dict2json_str(json_dict::AbstractDict)
    buf = IOBuffer()
    JSON.print(buf, json_dict, 4)
    return json_str = String(take!(buf))
end

"""
POSTされたJSONを保存する
"""
function postTest(req::HTTP.Request)::HTTP.Response
    println("createTest")
    request_body = String(HTTP.payload(req))
    try
        # String を parseして　TestStruct型変数作成
        test_struct = JSON2.read(request_body, TestStruct)
        # デフォルトはtest_structのnameをkeyとする
        push!(test_structs, test_struct.name=>test_struct)
        return HTTP.Response(200)
    catch
        return HTTP.Response(500)
    end
end

"""
GET，nameをpathから取得して，そのnameを持つデータをJSONで返す
"""
function getTest(req::HTTP.Request)::HTTP.Response
    println("getTest")
    try
        uri = HTTP.URI(req.target)
        params = HTTP.queryparams(uri)
        # testapi/name, get name
        name = HTTP.URIs.splitpath(uri.path)[2]
        test_struct::TestStruct = test_structs[name]
        json_str = JSON2.write(test_struct)
        json_dict = JSON2.read(json_str, Dict)
        json_str = dict2json_str(json_dict)
        return HTTP.Response(200, headers, body = bytes(json_str))
    catch
        body = """
        {
            "message": "No data"
        }
        """
        return HTTP.Response(404, headers, body = bytes(body))
    end
end

"""
PUT データを作成 or 更新する
nameを取得して，その他の値はbodyから取得する
あとはpostと同じ，request.bodyに"name":nameを追加してpostTestに投げる
"""
function putTest(req::HTTP.Request)::HTTP.Response
    println("putTest")
    try
        uri = HTTP.URI(req.target)
        params = HTTP.queryparams(uri)
        # testapi/name, get name
        name = HTTP.URIs.splitpath(uri.path)[2]
        json_dict = JSON2.read(String(req.body), Dict)
        # すでに"name"が存在する場合は，更新される
        push!(json_dict, "name"=>name)
        json_str = JSON2.write(json_dict)
        test_struct = JSON2.read(json_str, TestStruct)
        push!(test_structs, name=>test_struct)
        return HTTP.Response(200)
    catch
        body = """
        {
            "message": "No data"
        }
        """
        return HTTP.Response(404, headers, body = bytes(body))
    end
end

"""
DELETE データを削除する
pathからnameを取得して，そのnameを持つstructを削除する
"""
function deleteTest(req::HTTP.Request)::HTTP.Response
    println("deleteTest")
    try
        uri = HTTP.URI(req.target)
        params = HTTP.queryparams(uri)
        # testapi/name, get name
        name = HTTP.URIs.splitpath(uri.path)[2]
        # delete
        filter!(x->x.first != name, test_structs)
        println("deleted")
        return HTTP.Response(200)
    catch
        body = """
        {
            "message": "No data"
        }
        """
        return HTTP.Response(404, headers, body = bytes(body))
    end
end

const TEST_ROUTER = HTTP.Router()
HTTP.@register(TEST_ROUTER, "POST", "/testapi", postTest)
HTTP.@register(TEST_ROUTER, "GET", "/testapi/*", getTest)
HTTP.@register(TEST_ROUTER, "PUT", "/testapi/*", putTest)
HTTP.@register(TEST_ROUTER, "DELETE", "/testapi/*", deleteTest)
println("listen...")
HTTP.serve(TEST_ROUTER, localhost, port)

max_stars_repo_name: Soyukke/RESTAPIServer

max_issues_repo_name: Soyukke/RESTAPIServer

max_forks_repo_name: Soyukke/RESTAPIServer


text:
using EisensteinNumbers
using Base.Test: @test, @test_throws

# Eisenstein: same type arguments
@test begin
    l = Eisenstein(1, 2)
    r = Eisenstein{Int}(1, 2)
    l == r
end

# Eisenstein: mixed type arguments
@test begin
    l = Eisenstein(1, 1.2)
    r = Eisenstein{Float64}(1.0, 1.2)
    l == r
end

# Eisenstein: single argument
@test begin
    l = Eisenstein(3)
    r = Eisenstein(3, 0)
    l == r
end

# real
@test begin
    l = real(Eisenstein(1, 2))
    r = 0
    l == r
end

@test begin
    l = real(Eisenstein(1, 0))
    r = 1
    l == r
end

# imag
@test begin
    l = imag(Eisenstein(1, 2))
    r = sqrt(3)
    l == r
end

# complex
@test begin
    z = Eisenstein(1.2, 3.4)
    l = complex(z)
    r = complex(real(z), imag(z))
    l == r
end

# unreal: Eisenstein
@test begin
    l = unreal(Eisenstein(1.2, 3.4))
    r = 3.4
    l == r
end

# unreal: Real
@test begin
    l = unreal(Eisenstein(1, 0))
    r = 0
    l == r
end

# isreal: false case
@test begin
    l = isreal(Eisenstein(3, 4))
    r = false
    l == r
end

# isreal: true case
@test begin
    l = isreal(2 + 0ω)
    r = true
    l == r
end

# asarray
@test begin
    l = asarray(Eisenstein(5, 6))
    r = [5, 6]
    l == r
end

# iszero: false case
@test begin
    l = iszero(ω)
    r = false
    l == r
end

# iszero: true case
@test begin
    l = iszero(Eisenstein(0, 0))
    r = true
    l == r
end

# zero: type argument
@test begin
    l = zero(Eisenstein{Int})
    r = Eisenstein{Int}(0, 0)
    l == r
end

# zero: Eisenstein argument
@test begin
    l = zero(Eisenstein(1.2, 3.4))
    r = Eisenstein{Float64}(0.0, 0.0)
    l == r
end

# one: type argument
@test begin
    l = one(Eisenstein{Int})
    r = Eisenstein{Int}(1, 0)
    l == r
end

# one: Eisenstein argument
@test begin
    l = one(Eisenstein(1.2, 3.4))
    r = Eisenstein{Float64}(1.0, 0.0)
    l == r
end

# omega: type argument
@test begin
    l = omega(Eisenstein{Int})
    r = Eisenstein{Int}(0, 1)
    l == r
end

# omega: Eisenstein argument
@test begin
    l = omega(Eisenstein(1.2, 3.4))
    r = Eisenstein{Float64}(0.0, 1.0)
    l == r
end

# show
@test begin
    io = IOBuffer()
    show(io, Eisenstein(1, 2))
    l = String(take!(io))
    r = "[1: 1, ω: 2]"
    l == r
end

# conj
@test begin
    l = conj(Eisenstein(3, 4))
    r = Eisenstein(-1, -4)
    l == r
end

# Involutivity of conj
@test begin
    z = Eisenstein(2, 5)
    l = conj(conj(z))
    r = z
    l == r
end

# +: unary
@test begin
    l = +(Eisenstein(1, 2))
    r = Eisenstein(1, 2)
    l == r
end

# +: two Eisenstein numbers
@test begin
    x = Eisenstein(1, 2)
    y = Eisenstein(3, 4)
    l = x + y
    r = Eisenstein(4, 6)
    l == r
end

# +: add Eisenstein and Real
@test begin
    l = 2 + Eisenstein(5, 6)
    r = Eisenstein(7, 6)
    l == r
end

@test begin
    l = Eisenstein(7, 8) + 3
    r = Eisenstein(10, 8)
    l == r
end

# -: unary
@test begin
    l = -(Eisenstein(1.2, 3.4))
    r = Eisenstein(-1.2, -3.4)
    l == r
end

# -: two Eisenstein numbers
@test begin
    x = Eisenstein(1, 2)
    y = Eisenstein(3, 4)
    l = x - y
    r = Eisenstein(-2, -2)
    l == r
end

# -: add Eisenstein and Real
@test begin
    l = 2 - Eisenstein(5, 6)
    r = Eisenstein(-3, -6)
    l == r
end

@test begin
    l = Eisenstein(7, 8) - 3
    r = Eisenstein(4, 8)
    l == r
end

# ω is cubic root of unity
@test begin
    l = ω * ω * ω
    r = one(Eisenstein{Int})
    l == r
end

# *: multiply two Eisenstein numbers
@test begin
    l = Eisenstein(3, 1) * Eisenstein(2, -1)
    r = Eisenstein(7, 0)
    l == r
end

# *: commutativity
@test begin
    a = Eisenstein(1, 3)
    b = Eisenstein(2, 4)
    l = a * b
    r = b * a
    l == r
end

# *: multiply Eisenstein and Real
@test begin
    l = 2 * Eisenstein(3, 4)
    r = Eisenstein(6, 8)
    l == r
end

@test begin
    l = Eisenstein(3, 4) * 23
    r = Eisenstein(69, 92)
    l == r
end

# Positivity of abs2
@test begin
    l = abs2(Eisenstein(1, 2)) > 0
    r = true
    l == r
end

@test_throws ErrorException begin
    inv(zero(Eisenstein{Int}))
end

# /
@test begin
    x = Eisenstein(Int128(-50) // 9, -19 // 91)
    y = Eisenstein(Int128(-1) // 55, 91 // 8)
    z = x * y
    l = z / y
    r = x
    l == r
end

# \
@test begin
    x = Eisenstein(Int128(-50) // 9, -19 // 91)
    y = Eisenstein(Int128(-1) // 55, 91 // 8)
    z = x * y
    l = x \ z
    r = y
    l == r
end

# /: Eisenstein and Real
@test begin
    a = Int128(69) // 47
    z = Eisenstein(Int128(33) // 28, 13 // 35)
    l = a / z
    r = a * inv(z)
    l == r
end

@test begin
    a = Int128(6) // 19
    b = Int128(-25) // 71
    c = 34 // 1
    l = Eisenstein(b, c) / a
    r = Eisenstein(b / a, c / a)
    l == r
end

@test_throws ErrorException begin
    Eisenstein(3.1, 4.2) / 0.0
end

# \: Eisenstein and Real
@test begin
    a = Int128(69) // 47
    z = Eisenstein(Int128(33) // 28, 13 // 35)
    l = z \ a
    r = inv(z) * a
    l == r
end

@test begin
    a = Int128(6) // 19
    b = Int128(-25) // 71
    c = 34 // 1
    l = a \ Eisenstein(b, c)
    r = Eisenstein(a \ b, a \ c)
    l == r
end

@test_throws ErrorException begin
    0.0 \ Eisenstein(3.1, 4.2)
end

# associates
@test begin
    a, b, c, d, f, g = associates(Eisenstein(3, 7))
    abs2(a) == abs2(b) == abs2(c) == abs2(d) == abs2(f) == abs2(g)
end

# random
@test begin
    l = random(Eisenstein{Int64})
    r = random(Eisenstein{Int64})
    l != r
end


text:
using BenchmarkTools
using ArbFloats

round2(x) = round(x*100.0)/100.0

macro setprecisions(nbits)
  quote begin
    setprecision(BigFloat, $nbits)
    setprecision(ArbFloat, $nbits)
  end end
end

function bench(f,val)
    benchrunner = @benchmarkable ($f)($val)
    tune!(benchrunner)
    benchcatcher =  run(benchrunner).times
    return benchcatcher
end

function bench(f,val1,val2)
    benchrunner = @benchmarkable ($f)($val1,$val2)
    tune!(benchrunner)
    benchcatcher = run(benchrunner).times
    return benchcatcher
end

function benchbits_parsed_rel(f,val,nbits)
    setprecision(BigFloat, nbits)
    setprecision(ArbFloat, nbits)
    str = string(val)
    bigval = parse(BigFloat,str)
    arbval = parse(ArbFloat,str)
    bigbench = bench(f,bigval)
    arbbench = bench(f,arbval)
    big_slowdown = ratio( median(bigbench), median(arbbench) )
    return round2(big_slowdown)
  end

function benchbits_converted_rel(f,val,nbits)
    setprecision(BigFloat, nbits)
    setprecision(ArbFloat, nbits)
    bigval = convert(BigFloat,val)
    arbval = convert(ArbFloat,val)
    bigbench = bench(f,bigval)
    arbbench = bench(f,arbval)
    big_slowdown =  ratio( median(bigbench), median(arbbench) )
    return round2(big_slowdown)
end

function benchbits_converted_rel(f,val1,val2,nbits)
    setprecision(BigFloat, nbits)
    setprecision(ArbFloat, nbits)

    bigval1 = convert(BigFloat,val1)
    bigval2 = convert(BigFloat,val2)
    arbval1 = convert(ArbFloat,val1)
    arbval2 = convert(ArbFloat,val2)

    bigbench = bench(f,bigval1,bigval2)
    arbbench = bench(f,arbval1,arbval2)
    big_slowdown =  ratio( median(bigbench), median(arbbench) )
    return round2(big_slowdown)
end

function nbit_bigfloat_slowerby(fn, val, vecnbits)
    nprecisions = length(vecnbits)
    slowerby = zeros(Float64, nprecisions)
    for i in 1:nprecisions
         slowerby[i] = benchbits_converted_rel(fn, val, vecnbits[i])
    end
    return slowerby
end

function nbit_bigfloat_slowerby(fn, val1, val2, vecnbits)
    nprecisions = length(vecnbits)
    slowerby = zeros(Int, nprecisions)
    for i in 1:nprecisions
         slowerby[i] = benchbits_converted_rel(fn, val1, val2, vecnbits[i])
    end
    return slowerby
end

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
BenchmarkTools.DEFAULT_PARAMETERS.samples = 8
BenchmarkTools.DEFAULT_PARAMETERS.evals   = 1

#const fp_precise_bits = [128, 256, 512, 768, 1024, 1280, 1536, 2048, 2560, 3072, 3584];
# const fp_precise_bits = [128, 256, 512, 1024, 1280, 1536, 2048, 2560, 3072];
const fp_precise_bits = [128*i for i in 1:16];

println( "    bits: ",fp_precise_bits' );println()

g = Float64(golden); recipg = 1/g; ten=10.0; tenth = (1/ten);
#mul_slowerby = nbit_bigfloat_slowerby( (*), g, recipg, fp_precise_bits); println("    mul :",mul_slowerby' )
#div_slowerby = nbit_bigfloat_slowerby( (/), g, recipg, fp_precise_bits); println("    div :", div_slowerby' )
# sqrt_slowerby = nbit_bigfloat_slowerby( (sqrt), v1, v2, fp_precise_bits); println("    sqrt :", sqrt_slowerby' )
exp_slowerby = nbit_bigfloat_slowerby( (exp), g, fp_precise_bits);   println( "    exp : ", exp_slowerby' )
log_slowerby = nbit_bigfloat_slowerby( (log), g, fp_precise_bits);   println( "    log : ", log_slowerby' )
cos_slowerby = nbit_bigfloat_slowerby( (cos), g, fp_precise_bits);   println( "    cos : ", cos_slowerby' )
tan_slowerby = nbit_bigfloat_slowerby( (tan), g, fp_precise_bits);   println( "    tan : ", tan_slowerby' )
acos_slowerby = nbit_bigfloat_slowerby( (acos), recipg, fp_precise_bits); println( "    acos: ", acos_slowerby' )
atan_slowerby = nbit_bigfloat_slowerby( (atan), recipg, fp_precise_bits); println( "    atan: ", atan_slowerby' )
cosh_slowerby = nbit_bigfloat_slowerby( (cosh), g, fp_precise_bits);  println( "    cosh: ", cosh_slowerby' )
tanh_slowerby = nbit_bigfloat_slowerby( (tanh), g, fp_precise_bits);  println( "    tanh: ", tanh_slowerby' )
acosh_slowerby = nbit_bigfloat_slowerby( (acosh), g, fp_precise_bits);println( "    acosh: ", acosh_slowerby' )
atanh_slowerby = nbit_bigfloat_slowerby( (atanh), recipg, fp_precise_bits);println( "    atanh: ", atanh_slowerby' )

gamma_slowerby = nbit_bigfloat_slowerby( (gamma), g, fp_precise_bits);println( "    gamma:  ", gamma_slowerby' )
digamma_slowerby = nbit_bigfloat_slowerby( (digamma), g, fp_precise_bits);println( "    digamma:  ", digamma_slowerby' )
zeta_slowerby = nbit_bigfloat_slowerby( (zeta), g, fp_precise_bits);println( "     zeta:  ", zeta_slowerby' )

println();println( "    bits: ",fp_precise_bits' )

#println("     *      : ", mul_slowerby' );
#println("     /      : ", div_slowerby' );
#println("    cbrt    : ", cbrt_slowerby' );
println( "    exp    : ", exp_slowerby' );
println( "    log    : ", log_slowerby' );
println( "    cos    : ", cos_slowerby' );
println( "    tan    : ", tan_slowerby' );
println( "    acos   : ", acos_slowerby' );
println( "    atan   : ", atan_slowerby' );
println( "    cosh   : ", cosh_slowerby' );
println( "    tanh   : ", tanh_slowerby' );
println( "    acosh  : ", acosh_slowerby' );
println( "    atanh  : ", atanh_slowerby' );

println( "    gamma  : ", gamma_slowerby' );
println( "    digamma: ", digamma_slowerby' );
println( "    zeta   : ", zeta_slowerby' );

# fp_precise_bits and
bf_slower_by = hcat(  exp_slowerby ,log_slowerby,
           cos_slowerby ,tan_slowerby ,acos_slowerby, atan_slowerby,
           cosh_slowerby, tanh_slowerby, acosh_slowerby, atanh_slowerby,
           gamma_slowerby, digamma_slowerby ,zeta_slowerby)


text:
function OutlierFilter(arg0::jint, arg1::jdouble)
    return OutlierFilter((jint, jdouble), arg0, arg1)
end

function modify(obj::OutlierFilter, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
using DemeNet: Signer, Deme, ID, DemeID, Certificate, serialize, deserialize, datadir, Profile, Intent
using Base: UUID
using Sockets # do I need it?
#using ..Types: SystemConfig, PFID, BraidChain
#using ..Braiders: Mixer, Braider

#Recorder, register
#using ..Ledgers: serve # I could name it as ledger node or something
#using ..DataFormat
using PeaceVote: BraidChain, BraidChainServer, record
using Recruiters: Certifier
import Recruiters

import Recruiters: ticket, addtooken

### The configuration is stored in the ledger which can be transfered publically. One only needs to check that the configuration is signed by the server. So in the end one first downloads the ledger and then checks whether the configuration makes sense with serverid

### One should load system config with Deme. One thus would need to establish Ledger which thus would not require to have stuff. The constructor for Deme would be available here locally.

### Perhaps much better option would be to connect to the server over ssh since when one sets up the system one needs to have a key. We could use something like `scp` to copy a maintainerid and generate a server key executed by `julia "somefile.jl"` and returned that to a standart output.

# One would write it in the setup file as
# SERVER_ID = configure("pi@192.1.1.1",MAINTAINER_ID)
# And then latter on decide what duties the SERVER_ID would have
function configure(ssh,maintainerid)
    # (1) Installs julia (uses curl for downloads)
    # (2) Sets up packages PeaceFounder and etc...
    # (3) adds MAINTAINER_ID as an environment variable
    # (4) starts configurator server
    # (5) returns SERVER_ID
end

### The configurration file is constructed locally
function deploy(port,serverid,signedconfig)
    socket = connect(port)

    signature = signer.sign(config)
    serilaize(socket,(config,signature))

    serversignature = deserialize(socket)
    config, id = unwrap((config,signature))
    @assert id in serverid
end

function init(port,maintainerid,signer::Signer)
    server = listen(port)
    socket = accept(server)

    # I should use read wiht opcodes
    cmd = deserialize(socket)

    while true

        if cmd==:configure

            config = deserialize(socket)
            file, id = unwrap(config)
            @assert id in maintainerid
            serialize(socket,signer.sign(config))

        elseif cmd==:reset
            # verification
            tooken = rand(1:100)
            serialize(socket,tooken)
            signature = deserialize(socket)
            tooken, id = unwrap((tooken,signature))
            @assert id in maintainerid
            # removing all configuration except the serverkey
        elseif cmd==:pass
            return nothing
        end
    end
end

### We will need some improvement here on

### I might need to put this into Types.
#Port(port) = Port(getipaddr(),port)

struct PeaceFounderServer
    chain::BraidChain
    braidchain
    certifier
end

# I could perhaps start the server from DemeSpec file. In that way I could generate appropriate deme.

# This function starts the services which are specified fo this particular server.
# When onion sockets will become a thing the System should also be started on the participating members
# devices or some additionall dedicated mixers.

function storeprofile(cert::Certificate{Profile},deme::Deme)
    uuid = deme.spec.uuid
    intentprofile = Intent(cert,deme.notary)
    id = intentprofile.reference

    fname = datadir(uuid) * "/profiles/" * string(id,base=16)
    mkpath(dirname(fname))

    io = IOBuffer()
    serialize(io,cert)
    bytes = take!(io)

    write(fname,bytes)
end

### I could substitute chain with a deme and then create braidchain with in
function PeaceFounderServer(config::PeaceFounderConfig,chain::BraidChain,server::Signer)

    #config = deserialize(chain,PeaceFounderConfig)

    #config = SystemConfig(deme)

    #chain = BraidChain()
    braidchain = BraidChainServer(chain,server)

    certifier = Certifier(config.certifier,chain.deme,server)
    iddaemon = @async while true
        tooken,profile,cert = take!(certifier.tookencertifier.tickets)
        @show profile
        storeprofile(profile,chain.deme)
        record(config.braidchain.recorder,cert)
    end

    # if config.certifier==nothing
    #     certifier = nothing
    # else
    #     certifier = Certifier(config.certifier,chain.deme,server)
    #     iddaemon = @async while true
    #         tooken,profile,cert = take!(certifier.tookencertifier.tickets)
    #         @show profile
    #         storeprofile(profile,chain.deme)
    #         register(config.recorder,cert)
    #     end
    # end

    # mixer = Mixer(config.mixerport,chain.deme,server)
    # synchronizer = @async serve(config.syncport,chain.ledger)

    # braider = Braider(config.braider,chain.deme,server)

    ### So this is where the bug happens
    #recorder = Recorder(config.recorder,chain,braider,server)

    return PeaceFounderServer(chain,braidchain,certifier)
end

addtooken(config::PeaceFounderConfig,deme::Deme,tooken::Int,signer::Signer) = addtooken(config.certifier,deme,tooken,signer)

function ticket(config::PeaceFounderConfig,deme::Deme,tooken::Int)
    certifier = config.certifier
    cport = certifier.certifierport
    serverid = certifier.serverid
    ticket(deme.spec,cport,serverid,tooken) # could improve API
end

max_stars_repo_name: akels/DemocracyCredit

max_issues_repo_name: akels/DemocracyCredit

max_forks_repo_name: akels/DemocracyCredit


text:
# Define a dataset
a = "This is a string"

# Define lambda functions that perform actions

fsw = (x)->split(x," ")
fsl = (x)->split(x,"")

# Define the cells that will be nodes
w1 = j4pr.functioncell(x->x) 	# 1. copy-er
w2 = j4pr.functioncell(fsw)  	# 2. word splitter
w3 = j4pr.functioncell(fsl)  	# 3. letter counter
w41 = j4pr.functioncell(length)	# 4. length
w42 = j4pr.functioncell(length)  	# 5. length
w5 = j4pr.functioncell((x)->begin
		str = "the answer is: "*string(max(x...)/min(x...))
		@show str
		str
	end
) # 6. divide

# Define some form of connectivity
#	      / w2 -> w41 \
# Cin -> w1 ->		   -> w5 -> Cout
# 	      \ w3 -> w42 /

W=[w1,w2,w3,w41,w42,w5]
V= [ (0,1), (1,2), (1,3),(2,4), (3,5), (4,6), (5,6), (6,-1)] #vertices with indices from W

# Parse vertices and create async cells
V = unique(V)
C = [Channel(0) for i in 1:length(V)]
Cin = Channel(0);
Cout = Channel(0);

for (iw, w) in enumerate(W)
	ic = []; #input channels
	oc = []; #output channels
	for (iv, v) in enumerate(V)
		if iw == v[1] && v[2] != -1 oc=[oc...,C[iv]] end
		if iw == v[2] && v[1] != 0  ic=[ic...,C[iv]] end
		if iw == v[2] && v[1] == 0  ic=[ic...,Cin]  end
		if iw == v[1] && v[2] == -1 oc=[oc...,Cout] end
	end
	j4pr.@async_cell_exp w ic oc
end

# Put input, take output
print("Putting...")
tic()
put!(Cin, a)
println("OK")

println("The end with result $(take!(Cout))")
toc()

ps = w1+[w2+w41;w3+w42]+w5
@time r=a|> ps
@show r
[close(c) for c in C];


text:
############################################################################
# Find the ground state of the Ising Hamitlonian and get the magnetization #
############################################################################
using MatrixProductStates

let
    # Number of Spins
    N = 40
    # Bond Dimension
    D = 20
    # Physical dimension
    d = 2
    # Desired Accuracy
    acc = 1E-10
    # Maximum number of sweeps
    max_sweeps = 10

    # Coupling between nearest neigbours
    J = 1.0
    # Coupling to external field
    lambda = 1.0
    # Generate Ising Hamiltonian
    H = getIsingMPO(N, J, lambda)
    # Get the MPO for the total spin
    Sz = getTotalSpinMPO(N)

    # Run the variational ground state search
    E0, mps, num_of_sweeps = find_groundstate(H, D, d, acc, max_sweeps)

    # Print the results
    println("Groundstate energy:            ", E0)
    println("Groundstate energy density:    ", E0 / N)
    println("Number of sweeps:              ", num_of_sweeps)
    println("Total spin:                    ", expectation_value(mps, Sz))
    println("Entropy along center bond      ", compute_entropy(mps))

    nothing
end


text:
@doc """
Euler Equation -- test_viscous

The function tests for the correctness of the viscous terms computation.

This is a serial test and uses the input file called

`input_viscous_polynomial_ser.jl`

"""->

function test_viscous()

  @testset "--- Testing Viscous Terms ---" begin

    fname = "input_viscous_polynomial_ser.jl"
    # mesh, sbp, eqn, opts, pmesh = EulerEquationMod.createObjects(fname)
    mesh, sbp, eqn, opts = PDESolver.run_solver(fname)
    @assert mesh.isDG == true
    @assert opts["isViscous"] == true

    @testset "Checking Solution at Final Time Step" begin

      viscous_qvec_final_for_comparison = readdlm("viscous_files/solution_viscous_forcomparison_0.dat_keep")

      for dof_ix = 1:length(eqn.q_vec)
        diff = eqn.q_vec[dof_ix] - viscous_qvec_final_for_comparison[dof_ix]
        @test isapprox( diff, 0.0) atol=1e-13
      end

    end # End  testset("Checking Solution at Final Time Step")

    @testset "Checking Diffusion Tensor" begin

      Tsol = opts["Tsol"]
      Tdim = mesh.dim
      GtL = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
      GtR = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)

      # checking a couple interfaces
      interfaces = sview(mesh.interfaces, :)
      nfaces = length(mesh.interfaces)

      for interface = [3 4]
        fill!(GtL, 0.0)
        fill!(GtR, 0.0)

        q_faceL = Base.view(eqn.q_face, :, 1, :, interface)
        q_faceR = Base.view(eqn.q_face, :, 2, :, interface)

        NavierStokesMod.calcDiffusionTensor(eqn.params, q_faceL, GtL)
        NavierStokesMod.calcDiffusionTensor(eqn.params, q_faceR, GtR)

        GtL_vec = reshape(GtL, 128, )
        GtR_vec = reshape(GtR, 128, )

        GtL_file_to_check_against = string("viscous_files/viscous_test_face", interface, "_GtL_vec.dat_keep")
        GtR_file_to_check_against = string("viscous_files/viscous_test_face", interface, "_GtR_vec.dat_keep")

        GtL_vec_to_check_against = readdlm(GtL_file_to_check_against)
        GtR_vec_to_check_against = readdlm(GtR_file_to_check_against)

        for ii = 1:length(GtL_vec)

          diff = GtL_vec[ii] - GtL_vec_to_check_against[ii]
          @test isapprox( diff, 0.0) atol=1e-13

          diff = GtR_vec[ii] - GtR_vec_to_check_against[ii]
          @test isapprox( diff, 0.0) atol=1e-13

        end

      end   # end loop over interfaces

    end   # End  testset("Checking Diffusion Tensor")

  end   # End testset("--- Testing Viscous Terms ---")

end # End function test_viscous

add_func1!(NSTests, test_viscous, [TAG_SHORTTEST])


text:
using BenchmarkTools
using Random
using LinearAlgebra

BenchmarkTools.DEFAULT_PARAMETERS.evals = 20

# define arrays and matrices
const rows, cols = 500, 600
const reduceRows, reduceCols = Int(rows / 5), Int(cols / 6)

const float_arrayA = rand(rows * cols)
const float_arrayB = rand(rows * cols)
const small_int_matrixA = rand(1:10, reduceRows, reduceCols)
const small_int_matrixB = rand(1:10, reduceRows, reduceCols)
const float_matrixA = rand(rows, cols)
const float_matrixC = rand(cols, rows)
const small_float_matrixA = rand(reduceRows, reduceCols)
const small_float_matrixB = rand(reduceRows, reduceCols)

test1_name = "Element-wise sum of two $(reduceRows)x$(reduceCols) matrices (int), (1000 loops)"
function test1()
    for _ in 1:1000
        _ = small_int_matrixA + small_int_matrixB
    end
end

test2_name = "Element-wise multiplication of two $(reduceRows)x$(reduceCols) matrices (float64), (1000 loops)"
function test2()
    for _ in 1:1000
        _ = small_float_matrixA .* small_float_matrixB
    end
end

test3_name = "Dot (scalar) product of two $(rows * cols) arrays (float64), (1000 loops)"
function test3()
    for _ in 1:1000
        _ = dot(float_arrayA, float_arrayB)
    end
end

test4_name = "Matrix product of $(rows)x$(cols) and $(cols)x$(rows) matrices (float64)"
function test4()
    _ = float_matrixA * float_matrixC
end

test5_name = "L2 norm of $(rows)x$(cols) matrix (float64), (1000 loops)"
function test5()
    for _ in 1:1000
        _ = norm(float_matrixA)
    end
end

test6_name = "Sort of $(rows)x$(cols) matrix (float64)"
function test6()
    _ = sort(float_matrixA, dims = 2)
end

println(test1_name)
@btime test1()

println(test2_name)
@btime test2()

println(test3_name)
@btime test3()

println(test4_name)
@btime test4()

println(test5_name)
@btime test5()

println(test6_name)
@btime test6()