text:
const FN = typeof(identity) # `typeof(identity)` gives an object which is a DataType
const FC = typeof(conj)
const FA = typeof(adjoint)
const FT = typeof(transpose)
_conj(::FN) = conj
_conj(::FC) = identity
_conj(::FA) = transpose
_conj(::FT) = adjoint
_transpose(::FN) = transpose
_transpose(::FC) = adjoint
_transpose(::FA) = conj
_transpose(::FT) = identity
_adjoint(::FN) = adjoint
_adjoint(::FC) = transpose
_adjoint(::FA) = identity
_adjoint(::FT) = conj

"""
    abstract type AbstractStridedView{T,N,F<:Union{FN,FC,FA,FT}} <: AbstractArray{T,N} end

Abstract type for all kinds of StridedView.

# Type Parameters
- `T`: the type of the elements
- `N`: the dimension of the array
- `F<:Union{FN,FC,FA,FT}`: the function denotes the relation between the strided view and
                            its parent

The fields that needed in any concrete subtypes:
- `parent::A` where `A<:DenseArray`
- `size::NTuple{N,Int}`: size of the stridedview; default to be the size of `parent`
- `strides::NTuple{N,Int}`: strides of the stridedview; default to be strides of `parent`
- `offset::Int`: default to be 0; the position of the first element of the strided view is
                `offset + 1`.
- `op::F`: default to be `identity`
"""
abstract type AbstractStridedView{T,N,F<:Union{FN,FC,FA,FT}} <: AbstractArray{T,N} end

"""
    elsize(::Type{<:AbstractStridedView{T}})

Return the size of `T` which is the type of the elements.
"""
Base.elsize(::Type{<:AbstractStridedView{T}}) where {T} =
    Base.isbitstype(T) ? sizeof(T) :
        (Base.isbitsunion(T) ? Base.bitsunionsize(T) : sizeof(Ptr))

# Methods based on map!
Base.copy!(dst::AbstractStridedView{<:Any,N}, src::AbstractStridedView{<:Any,N}) where {N} =
    map!(identity, dst, src)
Base.conj!(a::AbstractStridedView{<:Real}) = a
Base.conj!(a::AbstractStridedView) = map!(conj, a, a)
LinearAlgebra.adjoint!(dst::AbstractStridedView{<:Any,N},
    src::AbstractStridedView{<:Any,N}) where {N} = copy!(dst, adjoint(src))
Base.permutedims!(dst::AbstractStridedView{<:Any,N},
    src::AbstractStridedView{<:Any,N}, p) where {N} = copy!(dst, permutedims(src, p))

"""
    _normalizestrides(size::Dims{N}, strides::Dims{N})

Return the nomalized strides with the following conventions:
If for `i`, `size[i] == 1`, set the corresponding `strides[i]` as `size[i-1]*strides[i-1]`.
If for `i`, `size[i] == 0`, set the `i`th and all the following strides elements after `i`
to be `size[i-1]*strides[i-1]`.

Used to construct general strided view type.
"""
_normalizestrides(size::Tuple{}, strides::Tuple{}) = strides
function _normalizestrides(size::Dims{N}, strides::Dims{N}) where {N}
    if size[1] == 1
        strides = Base.setindex(strides, 1, 1)
    elseif size[1] == 0
        strides = one.(strides)
    end
    for i = 2:N
        if size[i] == 1
            strides = Base.setindex(strides, size[i-1]*strides[i-1], i)
        elseif size[i] == 0
            for j = i:N
                strides = Base.setindex(strides, size[i-1]*strides[i-1], j)
            end
        end
    end
    return strides
end

"""
    _simplify(size::Dims{N}, strides::Dims{N}) where {N}

Simplify the `size` and `strides` based on the following rules:
(1) If `size[1] == 1`, move this dimension to the last one.
(2) If `size[1]*strides[1] == tailstrides[1]`, combine the second dimension to the first
and get the new first size as `size[1]*tailsize[1]` with `strides[1]`. We also add a one
dimensional space to the last, with strides `tailsize[1]*tailstrides[1]`.

Used in `sreshape` for the concrete type of strided view.

# Examples
```julia-repl
julia> _simplify((2, 1, 3, 4, 1), (1, 2, 2, 6, 24))
((24, 1, 1, 1, 1), (1, 24, 24, 2, 24))
```
"""
_simplify(size::Tuple{}, strides::Tuple{}) = size, strides
_simplify(size::Dims{1}, strides::Dims{1}) = size, strides
function _simplify(size::Dims{N}, strides::Dims{N}) where {N}
    tailsize, tailstrides = _simplify(tail(size), tail(strides))
    if size[1] == 1
        return (tailsize..., 1), (tailstrides..., strides[1])
    elseif size[1]*strides[1] == tailstrides[1]
        return (size[1]*tailsize[1], tail(tailsize)..., 1),
            (strides[1], tail(tailstrides)..., tailsize[1]*tailstrides[1])
    else
        return (size[1], tailsize...), (strides[1], tailstrides...)
    end
end

"""
    _computereshapestrides(newsize::Dims, oldsize::Dims{N}, oldstrides::Dims{N}) where {N}

Return the new strides after we reshape the strided view from the `oldsize` to the
`newsize`. The input `oldstrides` is the old strides before reshape.

Used in `sreshape` for the concrete type of strided view.
"""
_computereshapestrides(newsize::Tuple{}, oldsize::Tuple{}, oldstrides::Tuple{}) = oldstrides
function _computereshapestrides(newsize::Tuple{}, oldsize::Dims{N}, oldstrides::Dims{N}) where {N}
    all(isequal(1), oldsize) || throw(DimensionMismatch())
    return ()
end
function _computereshapestrides(newsize::Dims, oldsize::Tuple{}, oldstrides::Tuple{})
    all(isequal(1), newsize) || throw(DimensionMismatch())
    return newsize
end
function _computereshapestrides(newsize::Dims, oldsize::Dims{N}, oldstrides::Dims{N}) where {N}
    d,r = divrem(oldsize[1], newsize[1]) # The quotient and remainder (x÷y, x%y)
    # `oldsize[1]` should always be an integer multiple of `newsize[1]`.
    if r == 0
        s1 = oldstrides[1]
        if d == 1
            oldsize = (tail(oldsize)..., 1)
            oldstrides = (tail(oldstrides)..., newsize[1]*s1)
            return (s1, _computereshapestrides(tail(newsize), oldsize, oldstrides)...)
        else
            oldsize = (d, tail(oldsize)...)
            oldstrides = (newsize[1]*s1, tail(oldstrides)...)
            return (s1, _computereshapestrides(tail(newsize), oldsize, oldstrides)...)
        end
    else
        if prod(newsize) != prod(oldsize)
            throw(DimensionMismatch())
        else
            throw(ReshapeException())
        end
    end
end

"""
    sreshape(a::AbstractStridedView,newsize::Vararg{Int})

Return a strided view of that has the same parent with `a`, but with `newsize`.

This function is needed to be implemented for each concrete type of strided view.
"""
function sreshape end
sreshape(a::AbstractStridedView, newsize::Vararg{Int}) = sreshape(a, newsize)

struct ReshapeException <: Exception end
Base.show(io::IO, e::ReshapeException) = print(io, "Cannot produce a reshaped StridedView
    without allocating, try sreshape(copy(array), newsize) or fall back to reshape(array,
    newsize)")

# Used in getindex and setindex of the strided view
"""
    _computeind(indices::NTuple{N,Int}, strides::NTuple{N,Int}) where {N}

Return the index in `a.parent` that corresponds to the input index `indices` of the strided
view. `indices` is the index in the Cartesian form.
"""
@inline _computeind(indices::Tuple{}, strides::Tuple{}) = 1
@inline _computeind(indices::NTuple{N,Int}, strides::NTuple{N,Int}) where {N} =
    (indices[1]-1)*strides[1] + _computeind(tail(indices), tail(strides))

"""
    _computeviewsize(oldsize::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}

Return size after selecting the range `I` in the old strided view.

If `I[i]` is an Int, remove that dimension; if `I[i]` is a `Colon`, that dimension remain
invariant; if `I[i]` is a `RangeIndex`, the new size in that dimension is equal to
`length(I[i])`.
"""
@inline function _computeviewsize(oldsize::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Int)
        return _computeviewsize(tail(oldsize), tail(I))
    elseif isa(I[1], Colon)
        return (oldsize[1], _computeviewsize(tail(oldsize), tail(I))...)
    else
        return (length(I[1]), _computeviewsize(tail(oldsize), tail(I))...)
    end
end
_computeviewsize(::Tuple{}, ::Tuple{}) = ()

"""
    _computeviewstrides(oldstrides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}

Return the strides after selecting the range `I` in the old strided view.

If `I[i]` is an Int, remove that dimension; if `I[i]` is a `Colon`, leave that dimension
invariant; if `I[i]` is a `RangeIndex`, the new stride should be `oldstrides[1]*step(I[1])`.
"""
@inline function _computeviewstrides(oldstrides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Int)
        return _computeviewstrides(tail(oldstrides), tail(I))
    elseif isa(I[1], Colon)
        return (oldstrides[1], _computeviewstrides(tail(oldstrides), tail(I))...)
    else
        return (oldstrides[1]*step(I[1]), _computeviewstrides(tail(oldstrides), tail(I))...)
    end
end
_computeviewstrides(::Tuple{}, ::Tuple{}) = ()

"""
    _computeviewoffset(strides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}

Return the new offset after selecting range `I` from the old strided view.

If `I[i]` is a Colon, return `0` for that dimension; if `I[i]` is an Int or a RangeIndex,
return `first(I[1])-1)*strides[1]` for that dimension. Add all offsets for each dimension
to give the total new offset.
"""
@inline function _computeviewoffset(strides::NTuple{N,Int}, I::NTuple{N,SliceIndex}) where {N}
    if isa(I[1], Colon)
        return _computeviewoffset(tail(strides), tail(I))
    else
        return (first(I[1])-1)*strides[1]+_computeviewoffset(tail(strides), tail(I))
    end
end
_computeviewoffset(::Tuple{}, ::Tuple{}) = 0

"""
    struct ParentIndex
        i::Int
    end

The index directly into parent array. It is used to `getindex` from the `parent` of strided
view directly.
"""
struct ParentIndex
    i::Int
end

"""
    sview(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N}
    sview(a::AbstractStridedView, I::SliceIndex)
    sview(a::DenseArray{<:Any,N}, I::Vararg{SliceIndex,N}) where {N}
    sview(a::DenseArray, I::SliceIndex)

Return a view from index `I` of the strided view or a `DenseArray`as a new strided view.
"""
@inline sview(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(a, I...)
@inline sview(a::AbstractStridedView, I::SliceIndex) =
    getindex(sreshape(a, (length(a),)), I...)
@inline sview(a::DenseArray{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(StridedView(a), I...)
@inline sview(a::DenseArray, I::SliceIndex) =
    getindex(sreshape(StridedView(a), (length(a),)), I...)

"""
    Base.view(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N}

Return a new strided view from the index `I` of the old strided view `a`.
"""
@inline Base.view(a::AbstractStridedView{<:Any,N}, I::Vararg{SliceIndex,N}) where {N} =
    getindex(a, I...)

# linear algebra
"""
    LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N}, α::Number,
                        src::AbstractStridedView{<:Number,N}) where {N}

Return `dst = α * src`, where `α` is a number.
"""
function LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N}, α::Number,
                            src::AbstractStridedView{<:Number,N}) where {N}
    if α == 1
        copy!(dst, src)
    else
        dst .= α .* src
    end
    return dst
end

"""
    LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N},
                        src::AbstractStridedView{<:Number,N}, α::Number) where {N}

Return `dst = src * α`, where `α` is a number.
"""
function LinearAlgebra.mul!(dst::AbstractStridedView{<:Number,N},
                            src::AbstractStridedView{<:Number,N}, α::Number) where {N}
    if α == 1
        copy!(dst, src)
    else
        dst .= src .* α
    end
    return dst
end

"""
    LinearAlgebra.rmul!(dst::AbstractStridedView, α::Number)

Replace `dst` by `dst * α`, where `α` is a number.
"""
LinearAlgebra.rmul!(dst::AbstractStridedView, α::Number) = mul!(dst, dst, α)

"""
    LinearAlgebra.lmul!(α::Number, dst::AbstractStridedView)

Replace `dst` by `α * dst`, where `α` is a number.
"""
LinearAlgebra.lmul!(α::Number, dst::AbstractStridedView) = mul!(dst, α, dst)

"""
    isblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}

Check whether the strided view `A` is a blas matrix.

The definition of the blas matrix is that:

Case 1: if `A.op == identity`, it requires either `stride(A,1) == 1` or `stride(A,2) == 1`;

Case 2: if `A.op == conj`, it requires `stride(A, 2) == 1`.

Case 3: If `A.op` is neither `identity` nor `conj`, `A` is not a blas matrix.

The definition of blas matrix comes from the facts that the BLAS module requires contiguous
matrix columns.
"""
function isblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}
    if A.op == identity
        return stride(A,1) == 1 || stride(A,2) == 1
    elseif A.op == conj
        return stride(A, 2) == 1
    else # should never happen
        return false
    end
end

"""
    getblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}

Return the strided view with the compatible form of strides that can be used in BLAS
functions and the label which determine whether to transpose an argument in BLAS functions.

'N': The input matrix X is not transposed or conjugated.

'T': The input matrix X will be transposed.

'C': The input matrix X will be conjugated and transposed.
"""
function getblasmatrix(A::AbstractStridedView{T,2}) where {T<:LinearAlgebra.BlasFloat}
    if A.op == identity
        if stride(A,1) == 1
            return blasstrides(A), 'N'
        else
            return blasstrides(transpose(A)), 'T'
        end
    else
        return blasstrides(adjoint(A)), 'C'
    end
end

"""
    _threaded_blas_mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{T,2},
        B::AbstractStridedView{T,2}, α, β, nthreads) where {T<:LinearAlgebra.BlasFloat}

Return `C = A*B*α + C*β`. Use BLAS function `gemm!` and use multiple threads if
`nthreads > 1` and `size(C,1)*size(C,2) > 1024`.
"""
function _threaded_blas_mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{T,2},
        B::AbstractStridedView{T,2}, α, β, nthreads) where {T<:LinearAlgebra.BlasFloat}
    m, n = size(C)
    m == size(A, 1) && n == size(B, 2) || throw(DimensionMismatch())
    if nthreads == 1 || m*n < 1024
        A2, CA = getblasmatrix(A)
        B2, CB = getblasmatrix(B)
        C2 = blasstrides(C)
        LinearAlgebra.BLAS.gemm!(CA, CB, convert(T, α), A2, B2, convert(T, β), C2)
    else
        if m > n
            m2 = round(Int, m/16)*8
            nthreads2 = nthreads >> 1
            t = Threads.@spawn _threaded_blas_mul!(C[1:m2, :], A[1:m2, :], B, α, β, nthreads2)
            _threaded_blas_mul!(C[m2+1:m, :], A[m2+1:m, :], B, α, β, nthreads - nthreads2)
            wait(t)
            return C
        else
            n2 = round(Int, n/16)*8
            nthreads2 = nthreads >> 1
            t = Threads.@spawn _threaded_blas_mul!(C[:, 1:n2], A, B[:, 1:n2], α, β, nthreads2)
            _threaded_blas_mul!(C[:, n2+1:n], A, B[:, n2+1:n], α, β, nthreads - nthreads2)
            wait(t)
            return C
        end
    end
end

"""
    __mul!(C::AbstractStridedView{<:Any,2}, A::AbstractStridedView{<:Any,2},
            B::AbstractStridedView{<:Any,2}, α, β)

Return `C = A*B*α + C*β`.
This implementation is faster than LinearAlgebra.generic_matmatmul.
"""
function __mul!(C::AbstractStridedView{<:Any,2}, A::AbstractStridedView{<:Any,2},
                B::AbstractStridedView{<:Any,2}, α, β)
    (size(C,1) == size(A,1) && size(C,2) == size(B,2) && size(A,2) == size(B,1)) ||
        throw(DimensionMatch("A has size $(size(A)), B has size $(size(B)),
                                C has size $(size(C))"))
    m,n = size(C)
    k = size(A,2)
    A2 = sreshape(A, (m, 1, k))
    B2 = sreshape(permutedims(B,(2,1)), (1, n, k))
    C2 = sreshape(C, (m, n, 1))

    if α == 0 || k == 0
        rmul!(C, β)
    elseif α == 1
        if β == 0
            _mapreducedim!(*, +, zero, (m,n,k), (C2,A2,B2))
        elseif β == 1
            _mapreducedim!(*, +, nothing, (m,n,k), (C2,A2,B2))
        else
            _mapreducedim!(*, +, x->x*β, (m,n,k), (C2,A2,B2))
        end
    else
        f = (x,y)->(x*y*α)
        if β == 0
            _mapreducedim!(f, +, zero, (m,n,k), (C2,A2,B2))
        elseif β == 1
            _mapreducedim!(f, +, nothing, (m,n,k), (C2,A2,B2))
        else
            _mapreducedim!(f, +, x->x*β, (m,n,k), (C2,A2,B2))
        end
    end
    return C
end

# here we will have C.op == :identity && stride(C,1) < stride(C,2)
function _mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{T,2},
                B::AbstractStridedView{T,2}, α, β) where {T<:LinearAlgebra.BlasFloat}
    if stride(C,1) == 1 && isblasmatrix(A) && isblasmatrix(B)
        nthreads = use_threaded_mul() ? get_num_threads() : 1
        _threaded_blas_mul!(C, A, B, α, β, nthreads)
    else
        return __mul!(C, A, B, α, β)
    end
end

function LinearAlgebra.mul!(C::AbstractStridedView{T,2}, A::AbstractStridedView{<:Any,2},
    B::AbstractStridedView{<:Any,2}, α::Number = true, β::Number = false) where {T}
    if !(eltype(C) <: LinearAlgebra.BlasFloat && eltype(A) == eltype(B) == eltype(C))
        return __mul!(C, A, B, α, β)
    end
    # C.op is identity or conj
    if C.op == conj
        if stride(C,1) < stride(C,2)
            _mul!(conj(C), conj(A), conj(B), conj(α), conj(β))
        else
            _mul!(C', B', A', conj(α), conj(β))
        end
    elseif stride(C,1) > stride(C,2)
        _mul!(transpose(C), transpose(B), transpose(A), α, β)
    else
        _mul!(C, A, B, α, β)
    end
    return C
end

function LinearAlgebra.axpy!(a::Number, X::AbstractStridedView{<:Number,N},
                                Y::AbstractStridedView{<:Number,N}) where {N}
    if a == 1
        Y .= X .+ Y
    else
        Y .= a .* X .+ Y
    end
    return Y
end
function LinearAlgebra.axpby!(a::Number, X::AbstractStridedView{<:Number,N}, b::Number,
                                Y::AbstractStridedView{<:Number,N}) where {N}
    if b == 1
        axpy!(a, X, Y)
    elseif b == 0
        mul!(Y, a, X)
    else
        Y .= a .* X .+ b .* Y
    end
    return Y
end

# Converting back to other DenseArray type:
function Base.convert(T::Type{<:DenseArray}, a::AbstractStridedView)
    b = T(undef, size(a))
    copy!(StridedView(b), a)
    return b
end
# following method because of ambiguity warning
function Base.convert(::Type{T}, a::AbstractStridedView) where {T<:Array}
    b = T(undef, size(a))
    copy!(StridedView(b), a)
    return b
end
function Base.convert(::Type{Array}, a::AbstractStridedView{T}) where {T}
    b = Array{T}(undef, size(a))
    copy!(StridedView(b), a)
    return b
end


text:
using CompScienceMeshes
using Test

m1 = meshrectangle(1.0, 0.5, 0.1, 3)
m1 = CompScienceMeshes.translate(m1, [0.0, 0.5, 0.0])
m2 = meshrectangle(1.0, 0.5, 0.1, 3)

nc1 = numcells(m1)
nc2 = numcells(m2)

nv1 = numvertices(m1)
nv2 = numvertices(m2)

M1 = weld(m1,m2)
@test numcells(M1) == nc1 + nc2
@test numvertices(M1) == nv1 + nv2 - 11

M2 = weld(m1,-m2)
@test numcells(M2) == nc1 + nc2
@test numvertices(M2) == nv1 + nv2 - 11

ext: jl


text:
using ArgParse
include("../specfem_gll.jl/tasks/numerical_operation_on_bin_file.jl")
include("../specfem_gll.jl/src/utils/readfiles.jl")

function parse_command_line()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--target_basedir"
            help = "the target gll directory"
        "--reference_basedir"
            help = "the reference gll directory"
        "--mesh_basedir"
            help = "the mesh directory"
        "--tag"
            help = "the tag"
        "--iproc"
            help = "the bin file index"
    end
    return parse_args(s)
    end

function main()
    # parse args
    parsed_args = parse_command_line()
    target_basedir = parsed_args["target_basedir"]
    mesh_basedir = parsed_args["mesh_basedir"]
    reference_basedir = parsed_args["reference_basedir"]
    tag = parsed_args["tag"]
    iproc = parse(Int64, parsed_args["iproc"])
    # get nspec
    mesh_info = sem_mesh_read(mesh_basedir, 0)
    nspec = mesh_info.nspec
    # run generate_perturbation
    compare_model(target_basedir::String, reference_basedir::String, tag::String,  nspec::Int64, iproc::Int64)
end

main()

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
function sumdigit(n)
    numsum = 0
    for c in n
        numsum += parse(Int, c)
    end
    return numsum
end

while true
    n = readline()
    n == "0" && break
    println(sumdigit(n))
end


text:
# mean, median, std etc.. functions don't work on Distribution.Truncated instances,
# so we'll define them using a resampling approach

import Distributions.Truncated
import StatsBase:
    mean,
    median,
    middle,
    std,
    mode,
    quantile

mean(d::Truncated; n::Int = 10000) = mean(rand(d, n))
median(d::Truncated; n::Int = 10000) = median(rand(d, n))
middle(d::Truncated; n::Int = 10000) = median(rand(d, n))
std(d::Truncated; n::Int = 10000) = std(rand(d, n))
mode(d::Truncated; n::Int = 10000) = mode(rand(d, n))
quantile(d::Truncated, q; n::Int = 10000) = quantile(rand(d, n), q)


text:
using ImageFiltering, ImageFeatures, OffsetArrays, CoordinateTransformations, PaddedViews #Distances

const h_kernel = Kernel.DoG((5, 180), (5*sqrt(2), 180), (31, 901))
const α = pi/2 + linspace(-.02, .02, 10)
const wavelengths = (250, 600)

exiftool_base = joinpath(Pkg.dir("TrackRoots"), "deps", "src", "exiftool", "exiftool")
const exiftool = exiftool_base*(is_windows() ? ".exe" : "")

"""
File(path, time)
A type that holds the file path and time.
"""
struct DarkFile
    path::String
    time::Float64 # in hours
end

"""
Stage(timelapse)
A type that holds all the `FilePair`s, including how many there are, `n`.
"""
struct CalibStage
    timelapse::Vector{DarkFile}
    Δt::Float64 # in hours
    Δx::Float64 # in mm
    speed::Float64 # in mm/hour
    id::Int
    home::String
    base::String

    CalibStage(timelapse::Vector{DarkFile}, Δt::Float64, Δx::Float64, id::Int, home::String, base::String) = new(timelapse, Δt, Δx, Δx/Δt, id, home, base)
end

tohour{T}(t::T) = t/convert(T, Dates.Hour(1))

"""
parse2hours(txt)
Convert `txt` to milliseconds.
"""
function parse2hours(x::String)
    m = match(r"^(\d\d\d\d)(\d\d)(\d\d) (\d\d):(\d\d):(\d\d).?(\d?\d?\d?)", x)
    @assert m ≠ nothing "Failed to extract date and time from string $x"
    tohour(DateTime(parse.(Int, m.captures)...) - DateTime(0))
end

#="""
find_vertical_distances(file)
Find all the distances between vertical edges in an image. Helps detect the grid lines in the background of the light images.
"""
function find_vertical_distances(file::String)
    img = load(file)
    I = imfilter(img, h_kernel, "symmetric")
    threshold = quantile(vec(I), 0.99)
    img_edges = I .> threshold
    lines = hough_transform_standard(img_edges, 1, α, 40, 10)
    x = first.(lines)'
    return vec(pairwise(Cityblock(), x))
end=#

sort_sum(h::Vector{Float64}) = 1/sum(sort(h, rev=true)[1:4])
rotateim(img::Matrix{Float64}, θ::Float64) = warp(img, recenter(RotMatrix(θ), Images.center(img)), 0)
collapseim(img::OffsetArray{Float64,2,Matrix{Float64}}, dim::Int) = vec(parent(sum(img, dim)))
function rotate_sum(img::Matrix{Float64}, θ::Float64)
    imgw = rotateim(img, θ)
    h = collapseim(imgw, 1)
    # v = collapseim(imgw, 2)
    sort_sum(h)# + sort_sum(v)
end
# function shift_compare(v::Vector{Float64}, h::Vector{Float64}, step::Int)
    # v1 = circshift(v, step)
    # v1 + h
# end
# shift_sum(v::Vector{Float64}, h::Vector{Float64}, step::Int) = sort_sum(shift_compare(v, h, step))

function pixel_width(st::Stage)
    file1 = st.timelapse[1].light
    file2 = st.timelapse[end].light
    img1 = load(file1)
    img2 = load(file2)
    img = [(Float64(i1) + Float64(i2))/2 for (i1, i2) in zip(img1, img2)]
    _, i = findmin(rotate_sum(img, x) for x in α)
    imgw = rotateim(img, α[i])
    y = collapseim(imgw, 1)
    n = length(y)
    z = PaddedView(0, y', (1,n+1*wavelengths[2]))
    wls = wavelengths[1]:wavelengths[2]
    phs = 1:wavelengths[2]
    s = [sum(z[ph + i*wl]/n for i in 0:1) for wl in wls, ph in phs]
    _, i = findmax(s)
    iwl, iph = ind2sub(size(s), i)
    wl = wls[iwl]
    80/6/wl
end

# """
# pixel_width(stages)
# Return the first good-enough pixel width from all the images in this stack.
# """
# function pixel_width(stages::Vector{Stage})
    # for st in stages
        # x = pixel_width(st)
        # isempty(x) || return 80/6/x # in mm
        # file = st.timelapse[1].light
        # x = find_vertical_distances(file)
        # filter!(i -> 200 < i < 350, x)
        # isempty(x) || return 80/6/mean(x) # in mm
    # end
    # warn("failed to automatically calibrate the images. Assuming pixels are ~0.04 mm wide, this estimation would correspond to about 350 pixels between two grid lines (please check to make sure this estimation is not too far off)")
    # return 0.03819517804872148
# end

function stages2calib(stages::Vector{Stage}, Δx::Float64)
    nst = length(stages)
    ntl = length(first(stages).timelapse)
    files = Matrix{String}(ntl, nst)
    for j in 1:nst, i in 1:ntl
        files[i, j] = stages[j].timelapse[i].dark
    end
    ts = reshape([parse2hours(x) for x in eachline(`$exiftool -T -ModifyDate -n $files`)], ntl, nst)
    ts .-= RowVector(ts[1,:])
    Δt = mean(diff(ts, 1), 1)
    calibstages = Vector{CalibStage}(nst)
    for j in 1:nst
        dfs = Vector{DarkFile}(ntl)
        for i in 1:ntl
            dfs[i] = DarkFile(files[i,j], ts[i,j])
        end
        calibstages[j] = CalibStage(dfs, Δt[j], Δx, j, stages[j].home, stages[j].base)
    end
    return calibstages
end


text:
using ComponentArrays
using DifferentialEquations
using LabelledArrays
using Sundials
using Test
using Unitful

@testset "Performance" begin
    @testset "Issue 36" begin
        function f1(du,u,p,t)
            du.x .= -1 .* u.x .* u.y .* p[1]
            du.y .= -1 .* u.y .* p[2]
        end

        n = 100

        p = [0.1,0.1]

        lu_0 = @LArray fill(1000.0,2*n) (x=(1:n), y=(n+1:2*n))
        cu_0 = ComponentArray(x=fill(1000.0, n), y=fill(1000.0, n))

        lprob1 = ODEProblem(f1,lu_0,(0,100.0),p)
        cprob1 = ODEProblem(f1,cu_0,(0,100.0),p)

        solve(lprob1, Rodas5());
        solve(lprob1, Rodas5(autodiff=false));
        solve(cprob1, Rodas5());
        solve(cprob1, Rodas5(autodiff=false));

        ltime1 = @elapsed lsol1 = solve(lprob1, Rodas5());
        ltime2 = @elapsed lsol2 = solve(lprob1, Rodas5(autodiff=false));
        ctime1 = @elapsed csol1 = solve(cprob1, Rodas5());
        ctime2 = @elapsed csol2 = solve(cprob1, Rodas5(autodiff=false));

        @test (ctime1 - ltime1)/ltime1 < 0.3
        @test (ctime2 - ltime2)/ltime2 < 0.3
    end

    @testset "Slack Issue 2021-2-19" begin
        DifferentialEquations.DiffEqBase.diffeqbc(x::ComponentArray) = DifferentialEquations.DiffEqBase.DiffEqBC(x)

        function heat_conduction(du,u,p,t)
            nknots = 100
            h² = (1.0/(nknots+1))^2
            du .= zero(eltype(u))
            u₃ = @view u[3:end]
            u₂ = @view u[2:end-1]
            u₁ = @view u[1:end-2]
            @. du[2:end-1] = (u₃ - 2*u₂ + u₁)/h²
            nothing
        end

        t0, t1 = 0.0, 1.0
        u0 = randn(300)
        u0_ca = ComponentArray(a=u0[1:100],b=u0[101:200],c=u0[201:300])
        u0_la = @LArray u0 (a=1:100, b=101:200, c=201:300)

        cprob = ODEProblem(heat_conduction, u0_ca, (t0, t1))
        lprob = ODEProblem(heat_conduction, u0_la, (t0, t1))
        prob = ODEProblem(heat_conduction, u0, (t0, t1))

        solve(cprob, Tsit5(), saveat=0.2)
        solve(lprob, Tsit5(), saveat=0.2)
        solve(prob, Tsit5(), saveat=0.2)

        ctime = @elapsed solve(cprob, Tsit5(), saveat=0.2)
        ltime = @elapsed solve(lprob, Tsit5(), saveat=0.2)
        time = @elapsed solve(prob, Tsit5(), saveat=0.2)

        @test (ctime - time)/time < 0.1
        @test (ctime - ltime)/ltime < 0.05
    end
end

@testset "Issue 53" begin
    x0 = ComponentArray(x=ones(10))
    prob = ODEProblem((u,p,t)->u, x0, (0.,1.))
    sol = solve(prob, CVODE_BDF(linear_solver=:BCG), reltol=1e-15, abstol=1e-15)
    @test sol(1)[1] ≈ exp(1)
end

@testset "Issue 55" begin
    f!(D, x, p, t) = nothing
    x0 = ComponentArray(x=zeros(4))
    prob = ODEProblem(f!, x0, (0.0, 1.0), 0.0)
    sol = solve(prob, Rodas4())
    @test sol[1] == x0
end

@testset "Unitful" begin
    tspan = (0.0u"s", 10.0u"s")
    pos = 0.0u"m"
    vel = 0.0u"m/s"
    x0 = ComponentArray(pos=pos, vel=vel)
    F(t) = 1

    # double integrator in state-space form
    A = Union{typeof(0u"s^-1"), typeof(0u"s^-2"), Int}[0u"s^-1" 1; 0u"s^-2" 0u"s^-1"]
    B = Union{typeof(0u"m/s"), typeof(1u"m/s^2")}[0u"m/s"; 1u"m/s^2"]
    di(x,u,t) = A*x .+ B*u(t)

    prob = ODEProblem(di, x0, tspan, F)
    sol = solve(prob, Tsit5())
    @test unit(sol[end].pos) == u"m"
    @test unit(sol[end].vel) == u"m/s"
end


text:
using Metatheory
using Metatheory.Library

falseormissing(x) =
    x === missing || !x

r = @theory begin
    foo(~x, ~y) → 2 * ~x % ~y
    foo(~x, ~y) → sin(~x)
    sin(~x) → foo(~x, ~x)
end
@testset "Basic Equalities 1" begin
    @test (@areequal r foo(b, c) foo(d, d)) == false
end

r = @theory begin
    ~a * 1 → :foo
    ~a * 2 → :bar
    1 * ~a → :baz
    2 * ~a → :mag
end

@testset "Matching Literals" begin
    g = EGraph(:(a * 1))
    addexpr!(g, :foo)
    saturate!(g, r)
    display(g.classes); println()

    @test (@areequal r a * 1 foo) == true
    @test (@areequal r a * 2 foo) == false
    @test (@areequal r a * 1 bar) == false
    @test (@areequal r a * 2 bar) == true

    @test (@areequal r 1 * a baz) == true
    @test (@areequal r 2 * a baz) == false
    @test (@areequal r 1 * a mag) == false
    @test (@areequal r 2 * a mag) == true
    end

comm_monoid = @commutative_monoid (*) 1
@testset "Basic Equalities - Commutative Monoid" begin
    @test true == (@areequal comm_monoid a * (c * (1 * d)) c * (1 * (d * a)) )
    @test true == (@areequal comm_monoid x * y y * x )
    @test true == (@areequal comm_monoid (x * x) * (x * 1) x * (x * x) )
end

comm_group = @commutative_group (+) 0 inv
t = comm_monoid ∪ comm_group ∪ (@distrib (*) (+))

@testset "Basic Equalities - Comm. Monoid, Abelian Group, Distributivity" begin
    @test true == (@areequal t (a * b) + (a * c) a * (b + c) )
    @test true == (@areequal t a * (c * (1 * d)) c * (1 * (d * a)) )
    @test true == (@areequal t a + (b * (c * d)) ((d * c) * b) + a )
    @test true == (@areequal t (x + y) * (a + b) ((a * (x + y)) + b * (x + y)) ((x * (a + b)) + y * (a + b)) )
    @test true == (@areequal t (((x * a + x * b) + y * a) + y * b) (x + y) * (a + b) )
    @test true == (@areequal t a + (b * (c * d)) ((d * c) * b) + a )
    @test true == (@areequal t a + inv(a) 0 (x * y) + inv(x * y) 1 * 0 )
end

@testset "Basic Equalities - False statements" begin
    @test falseormissing(@areequal t (a * b) + (a * c) a * (b + a))
    @test falseormissing(@areequal t (a * c) + (a * c) a * (b + c))
    @test falseormissing(@areequal t a * (c * c) c * (1 * (d * a)))
    @test falseormissing(@areequal t c + (b * (c * d)) ((d * c) * b) + a)
    @test falseormissing(@areequal t (x + y) * (a + c) ((a * (x + y)) + b * (x + y)))
    @test falseormissing(@areequal t ((x * (a + b)) + y * (a + b)) (x + y) * (a + c))
    @test falseormissing(@areequal t (((x * a + x * b) + y * a) + y * b) (x + y) * (a + x))
    @test falseormissing(@areequal t a + (b * (c * a)) ((d * c) * b) + a)
    @test falseormissing(@areequal t a + inv(a) a)
    @test falseormissing(@areequal t (x * y) + inv(x * y) 1)
end

# Issue 21
simp_theory = @theory begin
    munit() => :foo
end
G = EGraph(:(munit()))
params = SaturationParams(timeout=1)
saturate!(G, simp_theory, params)

module Bar
foo = 42
export foo
    using Metatheory

    t = @theory begin
        :woo => foo
    end
    export t
end

module Foo
    foo = 12
    using Metatheory

    t = @theory begin
        :woo => foo
    end
    export t
end

g = EGraph(:woo);
saturate!(g, Bar.t);
saturate!(g, Foo.t);
foo = 12

@testset "Different modules" begin
    @test @areequalg g t 42 12
end

comm_monoid = @theory begin
    ~a * ~b --> ~b * ~a
    ~a * 1 --> ~a
    ~a * (~b * ~c) --> (~a * ~b) * ~c
    ~a::Number * ~b::Number => ~a * ~b
end

G = EGraph(:(3 * 4))
@testset "Basic Constant Folding Example - Commutative Monoid" begin
    @test (true == @areequalg G comm_monoid 3 * 4 12)
    @test (true == @areequalg G comm_monoid 3 * 4 12 4 * 3  6 * 2)
end

@testset "Basic Constant Folding Example 2 - Commutative Monoid" begin
    ex = :(a * 3 * b * 4)
    G = EGraph(ex)
    @test (true == @areequalg G comm_monoid (3 * a) * (4 * b) (12 * a) * b ((6 * 2) * b) * a)
end

@testset "Type Assertions in Ematcher" begin
    some_theory = @theory begin
        ~a * ~b --> ~b * ~a
        ~a::Number * ~b::Number --> matched(~a, ~b)
        ~a::Int64 * ~b::Int64 --> specific(~a, ~b)
        ~a * (~b * ~c) --> (~a * ~b) * ~c
    end

    g = EGraph(:(2 * 3))
    saturate!(g, some_theory)
    # display(g.classes)

    @test true == areequal(g, some_theory, :(2 * 3), :(matched(2, 3)))
    @test true == areequal(g, some_theory, :(matched(2, 3)), :(specific(3, 2)))
end

function Base.iszero(g::EGraph, ec::EClass)
    n = ENodeLiteral(0)
    return n ∈ ec
end

@testset "Predicates in Ematcher" begin
    some_theory = @theory begin
        ~a::iszero * ~b --> 0
        ~a * ~b --> ~b * ~a
    end

    g = EGraph(:(2 * 3))
    saturate!(g, some_theory)
    # display(g.classes)

    @test true == areequal(g, some_theory, :(a * b * 0), 0)
end


text:
quoted(x) = typof(x) === Quote || typof(x) === Quotenode
unquoted(x) = typof(x) === UnaryOpCall && typof(x.args[1]) === OPERATOR && kindof(x.args[1]) == CSTParser.Tokens.EX_OR

function get_ids(x, q = false, ids = [])
    if quoted(x)
        q = true
    end
    if q && unquoted(x)
        q = false
    end
    if isidentifier(x)
        !q && push!(ids, x)
    elseif x.args !== nothing
        for i in 1:length(x.args)
            get_ids(x.args[i], q, ids)
        end
    end
    ids
end

function collect_bindings_refs(x::EXPR, bindings = [], refs = [])
    if bindingof(x) !== nothing
        push!(bindings, x)
    end
    if StaticLint.hasref(x)
        push!(refs, x)
    end
    if x.args !== nothing
        for a in x.args
            collect_bindings_refs(a, bindings, refs)
        end
    end
    return bindings, refs
end

function remove_ref(x::EXPR)
    if hasref(x) && refof(x) isa Binding && refof(x).refs isa Vector
        for ia in enumerate(refof(x).refs)
            if ia[2] == x
                deleteat!(refof(x).refs, ia[1])
                setref!(x, nothing)
                return
            end
        end
        error()
    end
end

function clear_binding(x::EXPR)
    if bindingof(x) isa Binding
        for r in bindingof(x).refs
            if r isa EXPR
                setref!(r, nothing)
            elseif r isa Binding
                clear_binding(r)
            end
        end
        x.meta.binding = nothing
    end
end
function clear_scope(x::EXPR)
    if hasmeta(x) && scopeof(x) isa Scope
        setparent!(scopeof(x), nothing)
        empty!(scopeof(x).names)
        if typof(x) === FileH && scopeof(x).modules isa Dict && haskey(scopeof(x).modules, "Base") && haskey(scopeof(x).modules, "Core")
            m1 = scopeof(x).modules["Base"]
            m2 = scopeof(x).modules["Core"]
            empty!(scopeof(x).modules)
            scopeof(x).modules["Base"] = m1
            scopeof(x).modules["Core"] = m2
        else
            scopeof(x).modules = nothing
        end
    end
end

function clear_ref(x::EXPR)
    if refof(x) isa Binding
        if refof(x).refs isa Vector
            for i in 1:length(refof(x).refs)
                if refof(x).refs[i] == x
                    deleteat!(refof(x).refs, i)
                    break
                end
            end
        end
        setref!(x, nothing)
    elseif refof(x) !== nothing
        setref!(x, nothing)
    end
end
function clear_error(x::EXPR)
    if hasmeta(x) && x.meta.error !== nothing
        x.meta.error = nothing
    end
end
function clear_meta(x::EXPR)
    clear_binding(x)
    clear_ref(x)
    clear_scope(x)
    clear_error(x)
    if x.args !== nothing
        for a in x.args
            clear_meta(a)
        end
    end
end

function get_root_method(b, server)
    return b
end

function get_root_method(b::Binding, server, b1 = nothing, visited_bindings = Binding[])
    if b.prev === nothing || b == b.prev || !(b.prev isa Binding) || b in visited_bindings
        return b
    end
    push!(visited_bindings, b)
    if b.type == b.prev.type == CoreTypes.Function
        return get_root_method(b.prev, server, b, visited_bindings)
    elseif b.type == CoreTypes.Function && b.prev.type == CoreTypes.DataType
        return b.prev
    else
        return b
    end
end

function retrieve_delayed_scope(x)
    if (CSTParser.defines_function(x) || CSTParser.defines_macro(x)) && scopeof(x) !== nothing
        if parentof(scopeof(x)) !== nothing
            return parentof(scopeof(x))
        else
            return scopeof(x)
        end
    elseif typof(x) === Export
        return retrieve_scope(x)
    elseif parentof(x) !== nothing
        return retrieve_delayed_scope(parentof(x))
    end
    return nothing
end

function retrieve_scope(x)
    if scopeof(x) !== nothing
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_scope(parentof(x))
    end
    return
end

function find_return_statements(x::EXPR)
    rets = EXPR[]
    if CSTParser.defines_function(x)
        find_return_statements(x.args[3], true, rets)
    end
    return rets
end

function find_return_statements(x::EXPR, last_stmt, rets)
    if last_stmt && !(typof(x) === CSTParser.Block || typof(x) === CSTParser.If || typof(x) === CSTParser.KEYWORD)
        push!(rets, x)
        return rets, false
    end

    if typof(x) === CSTParser.Return
        push!(rets, x)
        return rets, true
    end

    if x.args isa Vector{EXPR}
        for i = 1:length(x.args)
            _, stop_iter = find_return_statements(x.args[i], last_stmt && (i == length(x.args) || (typof(x) === CSTParser.If && typof(x.args[i]) === CSTParser.Block)), rets)
            stop_iter && break
        end
    end
    return rets, false
end

function _expr_assert(x::EXPR, typ, nargs)
    typof(x) == typ && x.args isa Vector{EXPR} && length(x.args) == nargs
end

function _binary_assert(x, kind)
    typof(x) === CSTParser.BinaryOpCall && x.args isa Vector{EXPR} && length(x.args) == 3 && typof(x.args[2]) === CSTParser.OPERATOR && kindof(x.args[2]) === kind
end

# should only be called on Bindings to functions
function last_method(func::Binding)
    if func.next isa Binding && (func.next.type === CoreTypes.Function || (func.next.type === CoreTypes.DataType && func.type === CoreTypes.Function))
        return func.next
    else
        return func
    end
end

function prev_method(func::Binding)
    if func.prev isa Binding
        if func.prev.type === CoreTypes.Function
            return func.prev
        elseif func.prev.type === CoreTypes.DataType && func.prev.val isa EXPR && CSTParser.defines_struct(func.prev.val)
            return func.prev
        elseif (func.prev.val isa SymbolServer.FunctionStore || func.prev.val isa SymbolServer.DataTypeStore) && length(func.prev.val.methods) > 0
            return func.prev.val, 1
        end
    elseif (func.prev isa SymbolServer.FunctionStore || func.prev isa SymbolServer.DataTypeStore) && length(func.prev.methods) > 0
        return func.prev, 1
    end
    return nothing
end

function prev_method(func_iter::Tuple{T,Int}) where T <: Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore}
    func = func_iter[1]
    iter = func_iter[2]
    if iter < length(func.methods)
        return func, iter + 1
    else
        return nothing
    end
end

function find_exported_names(x::EXPR)
    exported_vars = EXPR[]
    for i in 1:length(x.args[3].args)
        expr = x.args[3].args[i]
        if typof(expr) == CSTParser.Export &&
            for j = 2:length(expr)
                if CSTParser.isidentifier(expr.args[j]) && StaticLint.hasref(expr.args[j])
                    push!(exported_vars, expr.args[j])
                end
            end
        end
    end
    return exported_vars
end

"""
    is_in_fexpr(x::EXPR, f)
Check whether `x` isa the child of an expression for which `f(parent) == true`.
"""
function is_in_fexpr(x::EXPR, f)
    if f(x)
        return true
    elseif parentof(x) isa EXPR
        return is_in_fexpr(parentof(x), f)
    else
        return false
    end
end

"""
    get_in_fexpr(x::EXPR, f)
Get the `parent` of `x` for which `f(parent) == true`. (is_in_fexpr should be called first.)
"""
function get_parent_fexpr(x::EXPR, f)
    if f(x)
        return x
    elseif parentof(x) isa EXPR
        return get_parent_fexpr(parentof(x), f)
    end
end

hasreadperm(p::String) = (uperm(p) & 0x04) == 0x04

# check whether a path is in (including subfolders) the julia base dir. Returns "" if not, and the path to the base dir if so.
function _is_in_basedir(path::String)
    i = findfirst(r".*base", path)
    i == nothing && return ""
    path1 = path[i]::String
    !hasreadperm(path1) && return ""
    !isdir(path1) && return ""
    files = readdir(path1)
    if all(f -> f in files, ["Base.jl", "coreio.jl", "essentials.jl", "exports.jl"])
        return path1
    end
    return ""
end


text:
# Starting point for dictionaries example for Learning Julia
# Dictionaries map unique keys to individual values

# Create a new dictionary
d1 = Dict("a" => 1, "b" => 2, "c" => 3)
println(d1, "\n")

d2 = Dict{String, String}()
d2["key1"] = "val1"
println(d2, "\n")

# look up an item
println(d1["b"], "\n")

# Add an item to a dictionary
d1["d"] = 4
println(d1, "\n")

# Test to see if an item is in a dictionary
println(haskey(d1, "e"), "\n")

# remove an item from a dictionary
delete!(d1, "c")
println(d1, "\n")

# get all the keys and values as arrays
println(keys(d1))
println(values(d1), "\n")

# Iterate over a dictionary's keys and values
for (key, val) in d1
    println(key, ", ", val)
end

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
function SingularValueDecomposition(arg0::RealMatrix)
    return SingularValueDecomposition((RealMatrix,), arg0)
end

function get_condition_number(obj::SingularValueDecomposition)
    return jcall(obj, "getConditionNumber", jdouble, ())
end

function get_covariance(obj::SingularValueDecomposition, arg0::jdouble)
    return jcall(obj, "getCovariance", RealMatrix, (jdouble,), arg0)
end

function get_inverse_condition_number(obj::SingularValueDecomposition)
    return jcall(obj, "getInverseConditionNumber", jdouble, ())
end

function get_norm(obj::SingularValueDecomposition)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_rank(obj::SingularValueDecomposition)
    return jcall(obj, "getRank", jint, ())
end

function get_s(obj::SingularValueDecomposition)
    return jcall(obj, "getS", RealMatrix, ())
end

function get_singular_values(obj::SingularValueDecomposition)
    return jcall(obj, "getSingularValues", Vector{jdouble}, ())
end

function get_solver(obj::SingularValueDecomposition)
    return jcall(obj, "getSolver", DecompositionSolver, ())
end

function get_u(obj::SingularValueDecomposition)
    return jcall(obj, "getU", RealMatrix, ())
end

function get_ut(obj::SingularValueDecomposition)
    return jcall(obj, "getUT", RealMatrix, ())
end

function get_v(obj::SingularValueDecomposition)
    return jcall(obj, "getV", RealMatrix, ())
end

function get_vt(obj::SingularValueDecomposition)
    return jcall(obj, "getVT", RealMatrix, ())
end


text:
# This file is a part of ElasticArrays.jl, licensed under the MIT License (MIT).

__precompile__(true)

module ElasticArrays

using Base: IteratorSize, HasLength, HasShape

include("util.jl")
include("elasticarray.jl")

end # module


text:
module TestSimpleEigensolver

using IterativeSolvers
using Test
using LinearMaps
using LinearAlgebra
using Random

@testset "Simple Eigensolvers" begin

Random.seed!(1234321)
n = 10

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)

    A = rand(T, n, n)
    A = A' * A
    λs = eigvals(A)

    tol = (eltype(T) <: Complex ? 2 : 1) * n^2 * cond(A) * eps(real(one(T)))

    ## Simple methods

    @testset "Power iteration" begin
        λ, x, history = powm(A; tol=tol, maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test λs[end] ≈ λ
        @test norm(A * x - λ * x) ≤ tol
    end

    @testset "Inverse iteration" begin
        # Set a target near the middle eigenvalue
        idx = div(n, 2)
        σ = T(0.75 * λs[idx] + 0.25 * λs[idx + 1])

        # Construct F = inv(A - σI) "matrix free"
        # Make sure we use complex arithmetic everywhere,
        # because of the follow bug in base: https://github.com/JuliaLang/julia/issues/22683
        F = lu(complex(A) - UniformScaling(σ))
        Fmap = LinearMap{complex(T)}((y, x) -> ldiv!(y, F, x), size(A, 1), ismutating = true)

        λ, x, history = invpowm(Fmap; shift=σ, tol=tol, maxiter=10n, log=true)

        @test isa(history, ConvergenceHistory)
        @test norm(A * x - λ * x) ≤ tol
        @test λ ≈ λs[idx]

    end
end
end

end # module


text:
#!/usr/bin/env julia
# File              : src/Airfoil.jl
# Author            : Denis N Dumoulin <denis.dumoulin@uclouvain.be>
# Date              : 18.03.2019
# Last Modified Date: 30.08.2021

"""
    gaw1(x)

Return the `y` coordinates of a GAW1 airfoil.
"""
function gaw1(x)
	ye = @. .278536sqrt(x) - .148567x + .006397x^2 - .220980x^3 + .081084x^4
	yi = @. -.190361sqrt(x)  + .161628x - .341176x^2 + .897631x^3 - .531252x^4
	return [ye, yi]
end

"""
    naca00(ZZ)(x)

Return the `y` coordinates of a NACA00ZZ airfoil.
"""
naca00(ZZ) = x -> begin
    ye = @. .05ZZ*(.2969sqrt(x) - .1260x - .3516x^2 + .2843x^3 - .1036x^4)
    yi = -ye
    return [ye, yi]
end

"""
    circle(R)(x)

Return the `y` coordinates of a circle of radius `R`.
"""
circle(R) = x -> begin
	ye = @. sqrt(R^2 - (x-R)^2)
	yi = -ye
	return [ye, yi]
end

"""
    ellipse(B)(x)

Return the `y` coordinates of an ellipse with half height `B`.
"""
ellipse(B) = x -> begin
	ye = @. .5b*sqrt(1 - 4(x-.5)^2)
	yi = -ye
	return [ye, yi]
end

"""
    uniform(vx, vy, α̇)(t)

Return the vector `[x, ẋ, ẍ]` corresponding to a uniform straight motion in
the `[-X, Y, -θ]` direction  with velocity `[vx, vy, α̇]`.
"""
uniform(vx, vy, α̇) = t -> ([-vx*t, vy*t, -α̇*t], [-vx, vy, -α̇], [0, 0, 0])

"""
    heavepitch(h0, αmax, ψ1, ψ2, strouhal)(t)

Return the vector `[x, ẋ, ẍ]` corresponding to a heaving & pitching motion where
`x = [-X, Y, α]`.
"""
heavepitch(h0, αmax, ψ1, ψ2, strouhal) = t -> begin
	ω = strouhal*pi/h0
	Y = h0*sin(ω*t + deg2rad(ψ1))
	Ẏ = ω*h0*cos(ω*t + deg2rad(ψ1))
	Ÿ = -(ω^2)*h0*sin(ω*t + deg2rad(ψ1))
	Y⃛ = -(ω^3)*h0*cos(ω*t + deg2rad(ψ1))

	αh = atan(-Ẏ)
	α̇h = -Ÿ / (1+Ẏ^2)
	α̈h = 2Ẏ*Ÿ^2/(1+Ẏ^2)^2 - Y⃛/(1 + Ẏ^2)

	α = deg2rad(αmax)*sin(ω*t + deg2rad(ψ2)) - αh
	α̇ = ω*deg2rad(αmax)*cos(ω*t + deg2rad(ψ2)) - α̇h
	α̈ = -ω^2*deg2rad(αmax)*sin(ω*t + deg2rad(ψ2)) - α̈h
	return [[-t, Y, α], [-1, Ẏ, α̇], [0, Ÿ, α̈]]
end

"""
    circularmotion(R, strouhal)(t)

Return the vector `[x, ẋ, ẍ]` corresponding to a circular motion of radius `R` where
`x = [X, Y, α]`.
"""
circularmotion(R, strouhal) = t -> begin
	ω = .5strouhal*pi/R
	ωt = ω*t

	X = R*cos(ωt)
	Ẋ = -ω*R*sin(ωt)
	Ẍ = -ω^2*R*cos(ωt)

	Y = R*sin(ωt)
	Ẏ = ω*R*cos(ωt)
	Ÿ = -ω^2*R*sin(ωt)

	return [[X, Y, -ωt+π/2], [Ẋ, Ẏ, -ω], [Ẍ, Ÿ, 0]]
end

#Constant values for the fish
const wh = .04
const sb = .04
const st = .95
const wt = .01

"""
    initiator(ξ, T)

Smooth increase from 0 to 1 in half a period `T` and its time derivative.
"""
function initiator(ξ, T)
	#= ξ >= T/2 && return [1, 0] =#
	#= return [.5 + .5sin(π*(2ξ/T - .5)), π/T*cos(π*(2ξ/T - .5))] =#
	ξ >= T && return [1, 0]
	return [.5 + .5sin(π*(ξ/T - .5)), .5π/T*cos(π*(ξ/T - .5))]
end

"""
    carlingfish(s, t, T)

Return the curvature of a curve of unit length associated with the motion of a
carling fish of period `T`. Its time derivative is computed too.
"""
function carlingfish(s, t, T)
	@assert 0 ≤ s ≤ 1 "Wrong coordinate 's'"
	a = .125
	b = .03125
	c = 1+b
	y = 2π*(s-t/T)
	κ1 = 4a/c*(π*cos(y) - π^2*(b+s)*sin(y))/(1+(a/c*(sin(y) + 2π*(b+s)*cos(y)))^2)^1.5
	κ2 = (75.3982*((1.0472a*c^2*(π*(b+s)*cos(y) + sin(y))) /((a^2*(2π*(b + s)*cos(y) + sin(y))^2)/c^2 + 1)^1.5
				   + (π*a^3*(2π*(b + s)*cos(y) + sin(y))*(cos(y) - 2π*(b + s)*sin(y))*(cos(y) - π*(b + s)*sin(y)))/((a^2*(2π*(b + s)*cos(y) + sin(y))^2)/c^2 + 1)^2.5))/(c^3*T)

	i1, i2 = initiator(t, T)
	return [κ1*i1; (κ1*i2 + κ2*i1)]
end

function anticarlingfish(s, t, T)
    -carlingfish(s, t, T)
end

function deadfish(s, t, T)
	zeros(2)
end

"""
    splinefish(s, t, T, X, K)

Return the curvature of a curve of unit length associated with the motion of a
fish following a cubic spline.
"""
function splinefish(s, t, T, X, K)
	@assert 0 ≤ s ≤ 1 "Wrong coordinate 's'"
	Ks = cubicspline(s, X, K)
	y = 2π*(t/T - 1.72s)
	κ1 = Ks*sin(y)
	κ2 = 2π/T*Ks*cos(y)

	i1, i2 = initiator(t, T)
	return [κ1*i1; κ1*i2 + κ2*i1]
end

function efficientfish(s, t, T)
	X = [0, 1/3, 2/3, 1]
	K = [3.34, 1.67, 2π, 2π]
	return splinefish(s, t, T, X, K)
end

function fastfish(s, t, T)
	X = [0, 1/3, 2/3, 1]
	K = [1.51, .48, 5.74, 2.73]
	return splinefish(s, t, T, X, K)
end

"""
    swimmerthickness(s)

Return the thickness of the swimming body along its coordinate `s`.
"""
function swimmerthickness(s)
	s < sb && return [sqrt(s*(2wh - s)), -sqrt(s*(2wh - s))]
	#= sb ≤ s < st && return [wh - (wh-wt)*((s-sb)/(st-sb))^2, -wh + (wh-wt)*((s-sb)/(st-sb))^2] =#
	sb ≤ s < st && return [wh - (wh-wt)*(s-sb)/(st-sb), -wh + (wh-wt)*(s-sb)/(st-sb)]
	return [wt*(1-s)/(1-st), -wt*(1-s)/(1-st)]
end


text:
using NamedVecs
using Test

@testset "NamedVecs.jl" begin
    v = NamedVec((a=[1], b=[2, 3], c=4))
    @test v == [1, 2, 3, 4]
    @test v.a == [1]
    @test v.b == [2, 3]
    @test v.c === 4

    v[2] = 5
    @test v == [1, 5, 3, 4]
    @test v.b == [5, 3]
    @test 2v == [2, 10, 6, 8]
    @test (2v).b == [10, 6]
    @test (2v).c === 8

    @test v .+ v == 2v
    @test 2 .* v == 2v
end


text:
using PatModules
@auxmodule one ("utils",) begin

    function say_hi()
        utils.say("hi")
    end
end


text:
e1 = Expr(:call, *, 3, 4) #> :((*)(3,4))
eval(e1) #> 12

# e1 = expr(:call, *, 3, a) #> ERROR: a not defined
e2 = Expr(:call, *, 3, :a) #> :((*)(3,a))
# eval(e2) #> ERROR: a not defined
a = 4 #> 4
eval(e2) #> 12

# b #> ERROR: b not defined
e3 = :(b = 1)
eval(e3) #> 1
b #> 1

e4 = :(a + b) #> :(a + b)
e5 = :($a + b) #> :(4 + b)
eval(e4) #> 5
eval(e5) #> 5

max_stars_repo_path: Module 1/Chapter07/eval.jl

max_issues_repo_path: Module 1/Chapter07/eval.jl

max_forks_repo_path: Module 1/Chapter07/eval.jl


text:
println("**********Testing Clustering**********")

facts("Run hclust on multiple FeatureVectors") do

  a = FeatureVector(Dict(["pickle"=>1, "wand"=>10]))
  b = FeatureVector(Dict(["word"=>131, "war"=>5]))
  c = FeatureVector(Dict(["saint"=>8, "waffle"=>44]))
  d = FeatureVector(Dict(["kid"=>97, "taco"=>18]))
  e = FeatureVector(Dict(["krap"=>49, "sock"=>22]))
  f = FeatureVector(Dict(["sack"=>201, "chartreuse"=>30]))
  g = FeatureVector(Dict(["knack"=>201, "balls"=>30]))
  i = FeatureVector(Dict(["snot"=>1, "falcon"=>2, "water"=>19]))
  j = FeatureVector(Dict(["fart"=>13, "snot"=>50, "calculator"=>11]))
  k = FeatureVector(Dict(["waffle"=>8, "feather"=>41, "water"=>19]))
  l = FeatureVector(Dict(["nard"=>197, "knock"=>45, "phanome"=>9]))
  m = FeatureVector(Dict(["knock"=>73, "sock"=>20, "butter"=>1]))
  n = FeatureVector(Dict(["knee"=>20, "falcon"=>300, "France"=>45]))
  o = FeatureVector(Dict(["snot"=>1, "falcon"=>2, "wrench"=>9, "taco"=>74]))
  p = FeatureVector(Dict(["fart"=>13, "snot"=>50, "water"=>74, "or"=>194]))
  q = FeatureVector(Dict(["waffle"=>8, "feather"=>41, "chartreuse"=>24, "pickle"=>6]))
  r = FeatureVector(Dict(["nard"=>197, "knock"=>45, "superhotdogmission"=>1, "androgyny"=>4]))
  s = FeatureVector(Dict(["knock"=>73, "sock"=>20, "taco"=>35, "hoop"=>16]))
  t = FeatureVector(Dict(["knee"=>20, "falcon"=>300, "fresh"=>88, "prince"=>88]))
  u = FeatureVector(Dict(["snot"=>1, "falcon"=>1, "wrench"=>1, "taco"=>1, "prince"=>1]))
  v = FeatureVector(Dict(["fart"=>1, "snot"=>1, "water"=>1, "or"=>1, "prince"=>1]))
  w = FeatureVector(Dict(["waffle"=>1, "feather"=>1, "chartreuse"=>1, "pickle"=>1, "prince"=>1]))
  x = FeatureVector(Dict(["nard"=>1, "knock"=>1, "superhotdogmission"=>1, "androgyny"=>1, "prince"=>1]))
  y = FeatureVector(Dict(["knock"=>1, "sock"=>1, "taco"=>35, "hoop"=>1, "prince"=>1]))
  z = FeatureVector(Dict(["knee"=>1, "falcon"=>1, "fresh"=>1, "prince"=>1, "princes"=>1]))

  data = [a,b,c,d,e,f,g,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z]
  h = hclust(data)
  @fact isa(h,BinaryTree) => true

end

facts("hclust deals returns nothing on incorrect input") do
  @fact hclust(3456787) => nothing
  @fact hclust("taco sauce") => nothing
  @fact hclust((4,5)) => nothing
end

facts("Random init provides k initial centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = random_init(features,2)

  @fact Base.length(centroids) => 2
  @fact centroids[1] != centroids[2] => true
end

facts("Max-min init provides k initial distinct centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = max_min_init(features,3,dist_cos)

  @fact Base.length(centroids) => 3
  @fact centroids[1] != centroids[2] => true
  @fact centroids[1] != centroids[3] => true
  @fact centroids[2] != centroids[3] => true
end

facts("Max-min init provides a random centroid and its furthest centroid") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["hello"=>8, "word"=>4]))
  d = FeatureVector(Dict(["hello"=>4, "word"=>8]))
  e = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  f = FeatureVector(Dict(["jello"=>4, "kim"=>0]))
  dict = Dict([1=>a,2=>b,3=>c,4=>d,5=>e,6=>f])
  features = collect(Base.values(dict))
  centroids = max_min_init(features,2,dist_cos)

  for fv in features
    @fact dist_cos(centroids[1], centroids[2]) >= dist_cos(centroids[1],fv) => true
  end
end

facts("kmeans returns correct amount of clusters") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  dict = Dict([1=>a,2=>b,3=>c])
  clusters = kmeans(dict,[],3)

  @fact Base.length(clusters) => 3
end

facts("kmeans returns no empty clusters") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>0]))
  b = FeatureVector(Dict(["hello"=>0, "word"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  dict = Dict([1=>a,2=>b,3=>c])
  clusters = kmeans(dict,[],3)

  for cluster in clusters
    @fact cluster.vectors != Cluster().vectors => true
  end
end

facts("kmeans returns the correct clustering of features") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>3]))
  b = FeatureVector(Dict(["part"=>3, "science"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  aa = FeatureVector(Dict(["hello"=>7, "word"=>5]))
  bb = FeatureVector(Dict(["part"=>8, "science"=>9]))
  cc = FeatureVector(Dict(["kim"=>2, "julia"=>3]))
  aaa = FeatureVector(Dict(["hello"=>1, "word"=>3]))
  bbb = FeatureVector(Dict(["part"=>3, "science"=>1]))
  ccc = FeatureVector(Dict(["kim"=>15, "julia"=>12]))
  dict = Dict([1=>a,2=>b,3=>c,4=>aa,5=>bb,6=>cc,7=>aaa,8=>bbb,9=>ccc])
  clusters = kmeans(dict,[],3)

  # what the clusters should be
  clust1_values = [a,aa,aaa]
  clust2_values = [b,bb,bbb]
  clust3_values = [c,cc,ccc]
  array = [clust1_values,clust2_values,clust3_values]

  # getting each cluster's FVs and making sure they are all there.
  for cluster in clusters
    fvs = collect(Base.values(cluster.vectors))
    for value in fvs
      @fact value in array => true
    end
  end
end

facts("kmeans lets you put in array of pre-defined centroids") do
  a = FeatureVector(Dict(["hello"=>4, "word"=>3]))
  b = FeatureVector(Dict(["part"=>3, "science"=>5]))
  c = FeatureVector(Dict(["kim"=>4, "julia"=>8]))
  aa = FeatureVector(Dict(["hello"=>7, "word"=>5]))
  bb = FeatureVector(Dict(["part"=>8, "science"=>9]))
  cc = FeatureVector(Dict(["kim"=>2, "julia"=>3]))
  aaa = FeatureVector(Dict(["hello"=>1, "word"=>3]))
  bbb = FeatureVector(Dict(["part"=>3, "science"=>1]))
  ccc = FeatureVector(Dict(["kim"=>15, "julia"=>12]))
  dict = Dict([1=>a,2=>b,3=>c,4=>aa,5=>bb,6=>cc,7=>aaa,8=>bbb,9=>ccc])
  array = [a,b,c]
  clusters = kmeans(dict,array)

  # what the clusters should be
  clust1_values = [a,aa,aaa]
  clust2_values = [b,bb,bbb]
  clust3_values = [c,cc,ccc]
  array = [clust1_values,clust2_values,clust3_values]

  # getting each cluster's FVs and making sure they are all there.
  for cluster in clusters
    fvs = collect(Base.values(cluster.vectors))
    for value in fvs
      @fact value in array => true
    end
  end
end


text:
# ccall wrapper for calling functions in libraries that might not be available

export @checked, with_workspace, @debug_ccall

"""
    @checked function foo(...)
        rv = ...
        return rv
    end

Macro for wrapping a function definition returning a status code. Two versions of the
function will be generated: `foo`, with the function body wrapped by an invocation of the
`@check` macro (to be implemented by the caller of this macro), and `unsafe_foo` where no
such invocation is present and the status code is returned to the caller.
"""
macro checked(ex)
    # parse the function definition
    @assert Meta.isexpr(ex, :function)
    sig = ex.args[1]
    @assert Meta.isexpr(sig, :call)
    body = ex.args[2]
    @assert Meta.isexpr(body, :block)

    # generate a "safe" version that performs a check
    safe_body = quote
        @check $body
    end
    safe_sig = Expr(:call, sig.args[1], sig.args[2:end]...)
    safe_def = Expr(:function, safe_sig, safe_body)

    # generate a "unsafe" version that returns the error code instead
    unsafe_sig = Expr(:call, Symbol("unsafe_", sig.args[1]), sig.args[2:end]...)
    unsafe_def = Expr(:function, unsafe_sig, body)

    return esc(:($safe_def, $unsafe_def))
end

"""
    with_workspace([eltyp=UInt8], size, [fallback::Int]; keep::Bool=false) do workspace
        ...
    end

Create a GPU workspace vector with element type `eltyp` and size in number of elements (in
the default case of an UInt8 element type this equals to the amount of bytes) specified by
`size`, and pass it to the do block. A fallback workspace size `fallback` can be specified
if the regular size would lead to OOM. Afterwards, the buffer is put back into the memory
pool for reuse (unless `keep` is set to `true`).

This helper protects against the rare but real issue of the workspace size getter returning
different results based on the GPU device memory pressure, which might change _after_
initial allocation of the workspace (which can cause a GC collection).
"""
@inline with_workspace(f, size::Union{Integer,Function}, fallback=nothing; kwargs...) =
    with_workspace(f, UInt8, isa(size, Integer) ? ()->size : size, fallback; kwargs...)

function with_workspace(f::Function, eltyp::Type{T}, size::Union{Integer,Function},
                        fallback::Union{Nothing,Integer}=nothing; keep::Bool=false) where {T}
    get_size() = Int(isa(size, Integer) ? size : size()::Integer)

    # allocate
    sz = get_size()
    workspace = nothing
    try
        while workspace === nothing || length(workspace) < sz
            workspace = CuVector{T}(undef, sz)
            sz = get_size()
        end
    catch ex
        fallback === nothing && rethrow()
        isa(ex, OutOfGPUMemoryError) || rethrow()
        workspace = CuVector{T}(undef, fallback)
    end
    workspace = workspace::CuVector{T}

    # use & free
    try
        f(workspace)
    finally
        keep || CUDA.unsafe_free!(workspace)
    end
end

macro debug_ccall(target, rettyp, argtyps, args...)
    @assert Meta.isexpr(target, :tuple)
    f, lib = target.args

    quote
        # get the call target, as e.g. libcuda() triggers initialization, even though we
        # can't use the result in the ccall expression below as it's supposed to be constant
        $(esc(target))

        # printing without task switches
        io = Core.stdout

        print(io, $f, '(')
        for (i, arg) in enumerate(($(map(esc, args)...),))
            i > 1 && print(io, ", ")
            render_arg(io, arg)
        end
        print(io, ')')
        rv = ccall($(esc(target)), $(esc(rettyp)), $(esc(argtyps)), $(map(esc, args)...))
        println(io, " = ", rv)
        for (i, arg) in enumerate(($(map(esc, args)...),))
            if arg isa Base.RefValue
                println(io, " $i: ", arg[])
            end
        end
        rv
    end
end

render_arg(io, arg) = print(io, arg)
render_arg(io, arg::Union{<:Base.RefValue, AbstractArray}) = summary(io, arg)


text:
function add_solver_xml!(::NEOSSolver{:CPLEX, :NL}, m::NLModel)
    # Add user options
    param_string = ""
    for key in keys(m.solver.params)
	param_string *= "set $(key) $(m.solver.params[key])\n"
    end
    m.xmlmodel = replace(m.xmlmodel, r"<options>.*</options>"s => "<options><![CDATA[$(param_string)]]></options>")
end


text:
#=
VecSE2: a 2d euclidean vector with an orientation
=#

immutable VecSE2 <: VecSE
    x :: Float64
    y :: Float64
    θ :: Float64

    VecSE2(x::Real=0.0, y::Real=0.0, θ::Real=0.0) = new(x, y, θ)
    VecSE2(a::VecE2, θ::Real=0.0) = new(a.x, a.y, θ)
end

polar(r::Real, ϕ::Real, θ::Real) = VecSE2(r*cos(ϕ), r*sin(ϕ), θ)

Base.length(::VecSE2) = 3
Base.copy(a::VecSE2) = VecSE2(a.x, a.y, a.θ)
Base.convert(::Type{Vector{Float64}}, a::VecSE2) = [a.x, a.y, a.θ]
Base.convert(::Type{VecE3}, a::VecSE2) = VecE3(a.x, a.y, a.θ)
Base.convert(::Type{VecE2}, a::VecSE2) = VecE2(a.x, a.y)
function Base.convert{R<:Real}(::Type{VecSE2}, a::AbstractArray{R})
    @assert(length(a) == 3)
    VecSE2(a[1], a[2], a[3])
end
Base.show(io::IO, a::VecSE2) = @printf(io, "VecSE2({%.3f, %.3f}, %.3f)", a.x, a.y, a.θ)

@compat Base.:+(b::Real, a::VecSE2) = VecSE2(a.x+b, a.y+b, a.θ)
@compat Base.:+(a::VecSE2, b::Real) = VecSE2(a.x+b, a.y+b, a.θ)
@compat Base.:+(a::VecSE2, b::VecE2) = VecSE2(a.x+b.x, a.y+b.y, a.θ)
@compat Base.:+(a::VecSE2, b::VecSE2) = VecSE2(a.x+b.x, a.y+b.y, a.θ+b.θ)

@compat Base.:-(b::Real, a::VecSE2) = VecSE2(b-a.x, b-a.y, a.θ)
@compat Base.:-(a::VecSE2, b::Real) = VecSE2(a.x-b, a.y-b, a.θ)
@compat Base.:-(a::VecSE2, b::VecE2) = VecSE2(a.x-b.x, a.y-b.y, a.θ)
@compat Base.:-(a::VecE2,  b::VecSE2) = VecE2(a.x-b.x, a.y-b.y)
@compat Base.:-(a::VecSE2, b::VecSE2) = VecSE2(a.x-b.x, a.y-b.y, a.θ-b.θ)

@compat Base.:*(b::Real, a::VecSE2) = VecSE2(b*a.x, b*a.y, a.θ)
@compat Base.:*(a::VecSE2, b::Real) = VecSE2(a.x*b, a.y*b, a.θ)

@compat Base.:/(a::VecSE2, b::Real) = VecSE2(a.x/b, a.y/b, a.θ)

@compat Base.:^(a::VecSE2, b::Integer) = VecSE2(a.x^b, a.y^b, a.θ)
@compat Base.:^(a::VecSE2, b::AbstractFloat) = VecSE2(a.x^b, a.y^b, a.θ)

@compat Base.:%(a::VecSE2, b::Real) = VecSE2(a.x%b, a.y%b, a.θ)

Base.clamp(a::VecSE2, lo::Real, hi::Real) = VecSE2(clamp(a.x, lo, hi), clamp(a.y, lo, hi), a.θ)

@compat Base.:(==)(a::VecSE2, b::VecSE2) = isequal(a.x, b.x) && isequal(a.y, b.y) && isequal(a.θ, b.θ)
Base.isequal(a::VecSE2, b::VecSE2) = isequal(a.x, b.x) && isequal(a.y, b.y) && isequal(a.θ, b.θ)
function Base.isapprox(a::VecSE2, b::VecSE2;
    _absa::Float64 = abs(a),
    _absb::Float64 = abs(b),
    _maxeps::Float64 = max(eps(_absa), eps(_absb)),
    rtol::Real=cbrt(_maxeps),
    atol::Real=sqrt(_maxeps)
    )

    isapprox(VecE2(a.x, a.y), VecE2(b.x, b.y), _absa=_absa, _absb=_absb, _maxeps=_maxeps, rtol=rtol, atol=atol) &&
    isapprox(a.θ, b.θ, rtol=rtol, atol=atol)
end

Base.isfinite(a::VecSE2) = isfinite(a.x) && isfinite(a.y) && isfinite(a.θ)
Base.isinf(a::VecSE2) = isinf(a.x) || isinf(a.y) || isinf(a.θ)
Base.isnan(a::VecSE2) = isnan(a.x) || isnan(a.y) || isnan(a.θ)

Base.abs(a::VecSE2) = hypot(a.x, a.y)
Base.hypot(a::VecSE2) = hypot(a.x, a.y)
Base.abs2(a::VecSE2) = a.x*a.x + a.y*a.y
function Base.norm(a::VecSE2)
    m = abs(a)
    VecSE2(a.x/m, a.y/m)
end

Base.atan2(a::VecSE2) = atan2(a.y, a.x)

function lerp(a::VecSE2, b::VecSE2, t::Real)
    x = a.x + (b.x-a.x)*t
    y = a.y + (b.y-a.y)*t
    θ = lerp_angle(a.θ, b.θ, t)
    VecSE2(x, y, θ)
end

Base.rot180(a::VecSE2) = VecSE2(a.x, a.y, a.θ+π)
Base.rotl90(a::VecSE2) = VecSE2(a.x, a.y, a.θ+0.5π)
Base.rotr90(a::VecSE2) = VecSE2(a.x, a.y, a.θ-0.5π)
rot(a::VecSE2, Δθ::Float64) = VecSE2(a.x, a.y, a.θ+Δθ)

Base.mod2pi(a::VecSE2) = VecSE2(a.x, a.y, mod2pi(a.θ))


text: abstract type AbstractRight <: AbstractThing end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
function resolve(op_sym, valstack)
	if op_sym == '+'
		operator = +
	elseif op_sym == '-'
		operator = -
	elseif op_sym == '*'
		operator = *
	end
	push!(valstack, operator(pop!(valstack), pop!(valstack)))
end

function eval_math(io)
	opstack = Char[]
	valstack = Int[]
	digits = IOBuffer()

	while !eof(io)
		c = read(io, Char)
		if c == '\n'
			#newline always results in returning the top of the valstack
			break
		elseif c == ')'
			resolve(pop!(opstack), valstack)
		elseif (c in ['+', '-', '*'] && peek(io, Char) == ' ')
			push!(opstack, c)
		elseif (isdigit(c) || c in ['-', '+'])
			write(digits, c)
			while isdigit(peek(io, Char))
				write(digits, read(io, Char))
			end
			push!(valstack, parse(Int, String(take!(digits))))
		end
		#ignoring any '(' and whitespaces except '\n'
	end

	return pop!(valstack)
end

#io = IOBuffer("(+ 2 (* (+ 23 -6) (* 3 7)))\n(+ 3 23)\n") # = 359
io = open("simple_program.scm", "r");

while !eof(io)
	result = eval_math(io)
	println(result)
end


text:
using Markdown # hide

# ## Automatic hyperparameter optimization
# By Eric S. Téllez
#

ENV["GKSwstype"] = "nul" # hide

using SimilaritySearch

const dim = 16

# This example optimizes different kinds of optimizations that allow different tradeoffs

X = rand(Float32, dim, 10^5)
Q = rand(Float32, dim, 10^3)
k = 10
nothing # hide

# We will use the squared L2, which preserves the order of L2 but is faster to compute.
dist = SqL2Distance()
db = MatrixDatabase(X)
queries = MatrixDatabase(Q)
nothing # hide
# Computing ground truth

goldI, goldD = searchbatch(ExhaustiveSearch(; db, dist), queries, k)  # `ExhaustiveSearch` solves with brute force
nothing # hide
# ## Different hyperparameter optimization strategies
# the way of specifying the hyperparameter optimization strategy and objective is with
# a `SearchGraphCallbacks` object, as follows:
G1 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks1 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=ParetoRecall()))  # ParetoRecall is the default
@elapsed index!(G1; callbacks=callbacks1)

# Using `ParetoRadius`: which should be faster since it doesn't needs a costly computation as the recall score
# but can be easily fool by distances distribution

G2 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks2 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=ParetoRadius()))  # it uses distances instead of recall, it will be faster but lower quality
@elapsed index!(G2; callbacks=callbacks2)

# Using `MinRecall`: It ensures a minimum quality in a small validation set
G3 = SearchGraph(; dist, db=MatrixDatabase(X))
callbacks3 = SearchGraphCallbacks(hyperparameters=OptimizeParameters(kind=MinRecall(0.95)))
@elapsed index!(G3; callbacks=callbacks3)

#
# `index!`, `append!`, and `push!` accept callbacks.
#
# ## Performances
# searching times
@time I1, D1 = searchbatch(G1, queries, k)
@time I2, D2 = searchbatch(G2, queries, k)
@time I3, D3 = searchbatch(G3, queries, k)
nothing #hide

# the recall (0-1, one is the best)
@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)
nothing #hide

#=
 here we can see smaller recalls than expected, but we also can appreciate the differences
 among them.

 ## Optimizing `SearchGraph` without inserting

 The hyperparameter optimization is performed in exponential stages while the `SearchGraph` is created,
 and therefore, the current hyperparameters could need an update. To optimize an already created `SearchGraph` we use `optimize` instead of `index`

=#
optimize!(G1, callbacks1.hyperparameters)
optimize!(G3, callbacks2.hyperparameters)
optimize!(G3, callbacks3.hyperparameters)

@time I1, D1 = searchbatch(G1, queries, k)
@time I2, D2 = searchbatch(G2, queries, k)
@time I3, D3 = searchbatch(G3, queries, k)
nothing #hide

# the recall ($0$ to  $1$, where one is the best)
@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)

#=
 The optimization is made using objects of the dataset as queries to compute the objective function.
 If it is possible to get external queries to optimize (not in the database and not in the query set,
 but also following the expected distrbiution), they can be provided as follows:
=#
equeries = MatrixDatabase(rand(dim, 64))
optimize!(G1, callbacks1.hyperparameters, queries=equeries)
optimize!(G2, callbacks2.hyperparameters, queries=equeries)
optimize!(G3, callbacks3.hyperparameters, queries=equeries)

@time I1, D1 = searchbatch(G1, queries, k)
@time I2, D2 = searchbatch(G2, queries, k)
@time I3, D3 = searchbatch(G3, queries, k)
nothing #hide

@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)

# Finally, we create our index with any hyperparameter optimization strategy, and optimize
# for quality, as follows:
minrecall = callbacks3.hyperparameters
optimize!(G1, minrecall, queries=equeries)
optimize!(G2, minrecall, queries=equeries)
optimize!(G3, minrecall, queries=equeries)

@time I1, D1 = searchbatch(G1, queries, k)
@time I2, D2 = searchbatch(G2, queries, k)
@time I3, D3 = searchbatch(G3, queries, k)
nothing #hide

@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)

#=
Please note that faster searches are expected for indexes created for higher qualities.

## `OptimizeParameters`
The structure is defined as follows

```
@with_kw mutable struct OptimizeParameters <: Callback
    kind::ErrorFunction = ParetoRecall()
    initialpopulation = 16
    params = SearchParams(maxpopulation=16, bsize=4, mutbsize=16, crossbsize=8, tol=-1.0, maxiters=16)
    ksearch::Int32 = 10
    numqueries::Int32 = 64
    space::BeamSearchSpace = BeamSearchSpace()
end
```

It allows to modify several parameters to adjusting parameters for a given objective,
these parameters can be tuned for faster constructions or more accurate ones.
=#

#- Using BeutifulMakie style
# ## Dependencies
using Pkg # hide
Pkg.status(["SimilaritySearch"]) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using MLBase

msg("BinaryClassEncoding: class hierachy")

@test BinaryClassEncoding <: ClassEncoding
@test ZeroOneClassEncoding <: BinaryClassEncoding
@test SignedClassEncoding <: BinaryClassEncoding
@test ClassEncodings.ZeroOneClassEncoding == ZeroOneClassEncoding
@test ClassEncodings.SignedClassEncoding == SignedClassEncoding

msg("BinaryClassEncoding: interface stability")

wrongDim1 = ["y","y","y"]
wrongDim2 = ["y","n","c"]

@test_throws ArgumentError ZeroOneClassEncoding(wrongDim1)
@test_throws ArgumentError ZeroOneClassEncoding(wrongDim2)
@test_throws ArgumentError SignedClassEncoding(wrongDim1)
@test_throws ArgumentError SignedClassEncoding(wrongDim2)

t = ["y", "n", "n", "y", "n"]
wrongLabel1 = ["a","b","a"]
wrongLabel2 = [1,2,2]

ce = ZeroOneClassEncoding(t)
@test_throws KeyError labelencode(ce, wrongLabel1)
@test_throws MethodError labelencode(ce, wrongLabel2)

ce = SignedClassEncoding(t)
@test_throws KeyError labelencode(ce, wrongLabel1)
@test_throws MethodError labelencode(ce, wrongLabel2)

#-----------------------------------------------------------

msg("BinaryClassEncoding: encoding and decoding")

t = ["y", "n", "n", "y", "n"]
y = ["n","y", "y", "n"]

ce = ZeroOneClassEncoding(t)
pred = labelencode(ce, y)
idx = groupindices(ce, t)
@test idx[1] == [1,4]
@test idx[2] == [2,3,5]
@test classDistribution(ce, t) == (["y","n"], [2, 3])
@test classDistribution(ce.labelmap, t) == (["y","n"], [2, 3])
@test pred == [1., 0, 0, 1]
@test labeldecode(ce, pred) == y

ce = ZeroOneClassEncoding(labelmap(["F","T"]))
@test labeldecode(ce, [1., 0, 1]) == ["T", "F", "T"]
@test labeldecode(ce, [0., 1, 0]) == ["F", "T", "F"]

ce = SignedClassEncoding(t)
pred = labelencode(ce, y)
idx = groupindices(ce, t)
@test idx[1] == [1,4]
@test idx[2] == [2,3,5]
@test classDistribution(ce, t) == (["y","n"], [2, 3])
@test classDistribution(ce.labelmap, t) == (["y","n"], [2, 3])
@test pred == [1., -1, -1, 1]
@test labeldecode(ce, pred) == y

ce = SignedClassEncoding(labelmap(["F","T"]))
@test labeldecode(ce, [1., -1, 1]) == ["T", "F", "T"]
@test labeldecode(ce, [-1., 1, -1]) == ["F", "T", "F"]

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

msg("MultinomialClassEncoding: class hierachy")

@test MultinomialClassEncoding <: ClassEncoding
@test MultivalueClassEncoding <: MultinomialClassEncoding
@test OneOfKClassEncoding <: MultinomialClassEncoding
@test OneHotClassEncoding <: MultinomialClassEncoding
@test OneHotClassEncoding == OneOfKClassEncoding
@test ClassEncodings.MultivalueClassEncoding == MultivalueClassEncoding
@test ClassEncodings.OneOfKClassEncoding == OneOfKClassEncoding

msg("MultinomialClassEncoding: interface stability")

wrongDim1 = ["y","y","y"]

@test_throws ArgumentError MultivalueClassEncoding(wrongDim1)
@test_throws ArgumentError OneOfKClassEncoding(wrongDim1)

t = ["y", "n", "k", "y", "n"]
wrongLabel1 = ["a","b","a"]
wrongLabel2 = [1,2,2]

ce = MultivalueClassEncoding(t)
@test_throws KeyError labelencode(ce, wrongLabel1)
@test_throws MethodError labelencode(ce, wrongLabel2)

ce = OneOfKClassEncoding(t)
@test_throws KeyError labelencode(ce, wrongLabel1)
@test_throws MethodError labelencode(ce, wrongLabel2)

#-----------------------------------------------------------

msg("MultinomialClassEncoding: encoding and decoding")

t = ["y", "n", "k", "y", "n"]
y = ["n", "y", "y", "k", "n"]

ce = MultivalueClassEncoding(t)
pred = labelencode(ce, y)
idx = groupindices(ce, t)
@test idx[1] == [1,4]
@test idx[2] == [2,5]
@test idx[3] == [3]
@test classDistribution(ce, t) == (["y","n","k"], [2, 2, 1])
@test classDistribution(ce.labelmap, t) == (["y","n","k"], [2, 2, 1])
@test pred == [2., 1, 1, 3, 2]
@test labeldecode(ce, pred) == y

ce = MultivalueClassEncoding(t, zero_based=true)
pred = labelencode(ce, y)
idx = groupindices(ce, t)
@test idx[1] == [1,4]
@test idx[2] == [2,5]
@test idx[3] == [3]
@test classDistribution(ce, t) == (["y","n","k"], [2, 2, 1])
@test classDistribution(ce.labelmap, t) == (["y","n","k"], [2, 2, 1])
@test pred == [1., 0, 0, 2, 1]
@test labeldecode(ce, pred) == y

ce = MultivalueClassEncoding(labelmap(["F","T", "W"]))
@test labeldecode(ce, [2., 1, 3]) == ["T", "F", "W"]
@test labeldecode(ce, [3., 2, 1]) == ["W", "T", "F"]

ce = OneOfKClassEncoding(t)
pred = labelencode(ce, y)
@test pred ==
  [0.  1  0;
   1  0  0;
   1  0  0;
   0  0  1;
   0  1  0]'
@test labeldecode(ce, pred) == y


text:
struct pwmEVA{T, U<:Distribution} <: fittedEVA{T}
    "Extreme value model definition"
    model::T
    "Maximum likelihood estimate"
    θ̂::Vector{Float64}
end

"""
    getdistribution(fittedmodel::pwmEVA)::Vector{<:Distribution}

Return the fitted distribution for the model fitted with the probability weigthed moments.

"""
function getdistribution(fittedmodel::pwmEVA)::Vector{<:Distribution}

    model = fittedmodel.model
    θ̂ = fittedmodel.θ̂

    res = getdistribution(model, θ̂)

    return res

end

"""
    quantile(fm::pwmEVA, p::Real)::Vector{<:Real}

Compute the quantile of level `p` from the fitted model. For the probability weighted moment method, the model has to be stationary.

"""
function quantile(fm::pwmEVA, p::Real)::Vector{<:Real}

    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    q = quantile(fm.model, fm.θ̂, p)

    return q

end

"""
    parametervar(fm::pwmEVA, nboot::Int=1000)

Estimate the parameter estimates covariance matrix by bootstrap.
"""
function parametervar(fm::pwmEVA, nboot::Int=1000)::Array{Float64, 2}

    @assert nboot>0 "the number of bootstrap samples should be positive."

    y = fm.model.data.value
    n = length(y)

    θ̂ = Array{Float64}(undef, nboot, length(fm.θ̂))

    fitfun = fitpwmfunction(fm)

    for i=1:nboot
        ind = rand(1:n, n)            # Generate a bootstrap sample
        θ̂[i,:] = fitfun(y[ind]).θ̂   # Compute the parameter estimates
    end

    V = cov(θ̂)                    # Compute the approximate covariance matrix

    return V

end

function cint(fm::pwmEVA, clevel::Real=.95, nboot::Int=5000)::Array{Array{Float64,1},1}

    @assert 0<clevel<1 "the confidence level should be between 0 and 1."
    @assert nboot>0 "the number of bootstrap samples should be positive."

    α = 1-clevel

    y = fm.model.data.value
    n = length(y)

    θ̂ = Array{Float64}(undef, nboot, length(fm.θ̂))

    fitfun = fitpwmfunction(fm)

    for i=1:nboot
        ind = rand(1:n, n)            # Generate a bootstrap sample
        θ̂[i,:] = fitfun(y[ind]).θ̂   # Compute the parameter estimates
    end

    confint = Vector{Vector{Float64}}()

    for c in eachcol(θ̂)
        push!(confint, quantile(c,[α/2; 1-α/2]))
    end

    return confint

end

"""
    fitpwmfunction(fm::pwmEVA{BlockMaxima, GeneralizedExtremeValue})::Function

Returns the corresponding fitpwm function.

"""
function fitpwmfunction(fm::pwmEVA{BlockMaxima, GeneralizedExtremeValue})::Function
    return gevfitpwm
end

"""
    fitpwmfunction(fm::pwmEVA{ThresholdExceedance, GeneralizedPareto})::Function

Returns the corresponding fitpwm function.

"""
function fitpwmfunction(fm::pwmEVA{ThresholdExceedance, GeneralizedPareto})::Function
    return gpfitpwm
end

"""
    fitpwmfunction(fm::pwmEVA{BlockMaxima, Gumbel})::Function

Returns the corresponding fitpwm function.

"""
function fitpwmfunction(fm::pwmEVA{BlockMaxima, Gumbel})::Function
    return gumbelfitpwm
end

"""
    quantilevar(fm::pwmEVA, level::Real, nboot::Int=1000)::Vector{<:Real}

Compute the  approximate variance of the quantile of level `level` from the fitted model `fm` by bootstrap.

"""
function quantilevar(fm::pwmEVA, level::Real, nboot::Int=1000)::Vector{<:Real}

    # Compute the approximate covariance matrice of the parameter estimates by bootstrap
    V = parametervar(fm, nboot)

    # Compute the approximate quantile variance by the delta method
    f(θ::DenseVector) = quantile(fm.model,θ,level)[]  # With the pwm method, the model is stationary
    Δf(θ::DenseVector) = ForwardDiff.gradient(f, θ)
    G = Δf(fm.θ̂)

    qv = G'*V*G

    return [qv]

end

"""
    returnlevel(fm::pwmEVA{BlockMaxima, T} where T<:Distribution, returnPeriod::Real)::ReturnLevel

Compute the confidence intervel for the return level corresponding to the return period
`returnPeriod` from the fitted model `fm` with confidence level `confidencelevel`.

"""
function returnlevel(fm::pwmEVA{BlockMaxima, T} where T<:Distribution, returnPeriod::Real)::ReturnLevel

      @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

      # quantile level
      p = 1-1/returnPeriod

      return ReturnLevel(BlockMaximaModel(fm), returnPeriod, quantile(fm, p))

end

function cint(rl::ReturnLevel{pwmEVA{BlockMaxima, T}} where T<:Distribution, confidencelevel::Real=.95, nboot::Int=1000)::Vector{Vector{Real}}

      @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
      @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

      # quantile level
      p = 1-1/rl.returnperiod

      # Compute the credible interval
      α = (1 - confidencelevel)

      y = rl.model.fm.model.data.value
      n = length(y)

      qboot = Array{Float64}(undef, nboot)

      fitfun = Extremes.fitpwmfunction(rl.model.fm)

      for i=1:nboot
          ind = rand(1:n, n)            # Generate a bootstrap sample
          θ̂ = fitfun(y[ind]).θ̂          # Compute the parameter estimates
          qboot[i] = quantile(rl.model.fm.model, θ̂, p)[]
      end

      return [quantile(qboot,[α/2, 1-α/2])]

end

"""
    returnlevel(fm::pwmEVA{ThresholdExceedance, T} where T<:Distribution, threshold::Real, nobservation::Int,
        nobsperblock::Int, returnPeriod::Real)::ReturnLevel

Compute the return level corresponding to the return period `returnPeriod` from the fitted model `fm`.

The threshold should be a scalar. A varying threshold is not yet implemented.

"""
function returnlevel(fm::pwmEVA{ThresholdExceedance, T} where T<:Distribution, threshold::Real, nobservation::Int,
    nobsperblock::Int, returnPeriod::Real)::ReturnLevel

    @assert returnPeriod > zero(returnPeriod) "the return period should be positive."

    # Exceedance probability
    ζ = length(fm.model.data.value)/nobservation

    # Appropriate quantile level given the probability exceedance and the number of obs per year
    p = 1-1/(returnPeriod * nobsperblock * ζ)

    return ReturnLevel(PeakOverThreshold(fm, threshold, nobservation, nobsperblock),
        returnPeriod, threshold .+ quantile(fm, p))

end

function cint(rl::ReturnLevel{pwmEVA{ThresholdExceedance, T}} where T<:Distribution, confidencelevel::Real=.95)::Vector{Vector{Real}}

    @assert rl.returnperiod > zero(rl.returnperiod) "the return period should be positive."
    @assert zero(confidencelevel)<confidencelevel<one(confidencelevel) "the confidence level should be in (0,1)."

    # Exceedance probability
    ζ = length(rl.model.fm.model.data.value)/rl.model.nobservation

    # Appropriate quantile level given the probability exceedance and the number of obs per year
    p = 1-1/(rl.returnperiod * rl.model.nobsperblock * ζ)

    # Compute the credible interval

    nboot = 5000
    α = (1 - confidencelevel)

    y = rl.model.fm.model.data.value
    n = length(y)

    qboot = Array{Float64}(undef, nboot)

    fitfun = Extremes.fitpwmfunction(rl.model.fm)

    for i=1:nboot
        ind = rand(1:n, n)            # Generate a bootstrap sample
        θ̂ = fitfun(y[ind]).θ̂          # Compute the parameter estimates
        qboot[i] = quantile(rl.model.fm.model, θ̂, p)[]
    end

    return [rl.model.threshold .+ quantile(qboot,[α/2, 1-α/2])]

end

"""
    showfittedEVA(io::IO, obj::pwmEVA; prefix::String = "")

Displays a pwmEVA with the prefix `prefix` before every line.

"""
function showfittedEVA(io::IO, obj::pwmEVA; prefix::String = "")

    println(io, prefix, "pwmEVA")
    println(io, prefix, "model :")
    showEVA(io, obj.model, prefix = prefix*"\t")
    println(io)
    println(io, prefix, "θ̂  :\t", obj.θ̂)

end


text:
function MidpointIntegrator(arg0::jdouble)
    return MidpointIntegrator((jdouble,), arg0)
end

function get_a(obj::MidpointIntegrator)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::MidpointIntegrator)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::MidpointIntegrator)
    return jcall(obj, "getC", Vector{jdouble}, ())
end


text:
module DeltaSigma

export
    calculateSNR,
    calculateTF,
    dbp,
    dbv,
    designHBF,
    ds_hann,
    dsdemo1,
    dsdemo2,
    dsdemo3,
    dsdemo5,
    dsdemo6,
    dsdemo7,
    exampleHBF,
    find2dPIS,
    mapABCD,
    partitionABCD,
    peakSNR,
    predictSNR,
    realizeNTF,
    rmsGain,
    scaleABCD,
    simulateDSM,
    simulateHBF,
    simulateMS,
    simulateSNR,
    stuffABCD,
    synthesizeNTF,
    undbp,
    undbv

# Key functions
include("calculateTF.jl")
include("mapABCD.jl")
include("realizeNTF.jl")
include("scaleABCD.jl")
include("simulateDSM.jl")
include("simulateMS.jl")
include("simulateSNR.jl")
include("stuffABCD.jl")
include("synthesizeNTF.jl")

# Speciality Functions
include("designHBF.jl")
include("exampleHBF.jl")
include("predictSNR.jl")
include("simulateHBF.jl")
include("PosInvSet/find2dPIS.jl")

# Utility functions
include("calculateSNR.jl")
include("dbp.jl")
include("dbv.jl")
include("ds_hann.jl")
include("partitionABCD.jl")
include("peakSNR.jl")
include("rmsGain.jl")
include("undbp.jl")
include("undbv.jl")

# Demos
include("demos/dsdemo1.jl")
include("demos/dsdemo2.jl")
include("demos/dsdemo3.jl")
include("demos/dsdemo5.jl")
include("demos/dsdemo6.jl")
include("demos/dsdemo7.jl")

# Internal helper functions
include("bquantize.jl")
include("bunquantize.jl")
include("delay.jl")
include("ds_optzeros.jl")
include("ds_therm.jl")
include("evalF0.jl")
include("evalF1.jl")
include("evalRPoly.jl")
include("evalTF.jl")
include("frespHBF.jl")
include("selectElement.jl")
include("PosInvSet/dscut.jl")
include("PosInvSet/dsmap.jl")
include("PosInvSet/dssplit2d.jl")
include("PosInvSet/hull2d.jl")
include("PosInvSet/leftof.jl")
include("PosInvSet/outconvex2d.jl")

# Internal graphing functions
include("axisLabels.jl")
include("figureMagic.jl")
include("logsmooth.jl")
include("plotPZ.jl")
include("plotSpectrum.jl")
include("plotUsage.jl")
include("PosInvSet/dotplot.jl")
include("PosInvSet/polyplot.jl")

end # module

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
"""
`init_subspec!(m::SmetsWouters)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::SmetsWouters)

    if subspec(m) == "ss0"
        return
    else
        error("This subspec is not defined.")
    end

end


text:
@testset "HTTP.Form for multipart/form-data" begin
    headers = Dict("User-Agent" => "HTTP.jl")
    body = HTTP.Form(Dict())
    uri = "https://httpbin.org/post"
    uri_put = "https://httpbin.org/put"
    @testset "Setting of Content-Type" begin
        for r in (HTTP.request("POST", uri, headers, body), HTTP.post(uri, headers, body),
                  HTTP.request("PUT", uri_put, headers, body), HTTP.put(uri_put, headers, body))
            @test r.status == 200
            json = JSON.parse(IOBuffer(HTTP.payload(r)))
            @test startswith(json["headers"]["Content-Type"], "multipart/form-data; boundary=")
        end
    end
    @testset "HTTP.Multipart ensure show() works correctly" begin
        # testing that there is no error in printing when nothing is set for filename
        str = sprint(show, (HTTP.Multipart(nothing, IOBuffer("some data"), "plain/text", "", "testname")))
        @test findfirst("contenttype=\"plain/text\"", str) !== nothing
    end
    @testset "HTTP.Multipart test constructor" begin
        @test_nowarn HTTP.Multipart(nothing, IOBuffer("some data"), "plain/text", "", "testname")
        @test_throws MethodError HTTP.Multipart(nothing, "some data", "plain/text", "", "testname")
    end

    @testset "Boundary" begin
        @test HTTP.Form(Dict()) isa HTTP.Form
        @test HTTP.Form(Dict(); boundary="a") isa HTTP.Form
        @test HTTP.Form(Dict(); boundary=" Aa1'()+,-.:=?") isa HTTP.Form
        @test HTTP.Form(Dict(); boundary='a'^70) isa HTTP.Form
        @test_throws ArgumentError HTTP.Form(Dict(); boundary="")
        @test_throws ArgumentError HTTP.Form(Dict(); boundary='a'^71)
        @test_throws ArgumentError HTTP.Form(Dict(); boundary="a ")
    end
end


text:
push!(LOAD_PATH,"../src/")
using NeatText
using Documenter

DocMeta.setdocmeta!(NeatText, :DocTestSetup, :(using NeatText); recursive=true)

makedocs(;
    modules=[NeatText],
    authors="JCharis <jcharistech.com> and contributors",
    repo="https://github.com/jcharistech/NeatText.jl/blob/{commit}{path}#{line}",
    sitename="NeatText.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jcharistech.github.io/NeatText.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "User Guide" => "userguide.md",
        "API Reference" => "apireference.md",

    ],
)

deploydocs(;
    repo="github.com/jcharistech/NeatText.jl",
)