text:
module VectorizationBase

import ArrayInterface, LinearAlgebra, Libdl, Hwloc, IfElse, LayoutPointers
using ArrayInterface:
    StaticInt, Zero, One, StaticBool, True, False,
    contiguous_axis, contiguous_axis_indicator, contiguous_batch_size, stride_rank,
    device, CPUPointer, CPUIndex, eq, ne, lt, le, gt, ge,
    known_length, known_first, known_last, strides, offsets,
    static_first, static_last, static_length
import IfElse: ifelse

using CPUSummary:
  cache_type, num_cache, num_threads, num_cache_levels, num_cores, num_l1cache,
  num_l2cache, cache_associativity, num_l3cache, sys_threads, cache_inclusive,
  num_l4cache, cache_linesize, num_machines, cache_size, num_sockets
using HostCPUFeatures:
  register_size, static_sizeof, fast_int64_to_double, pick_vector_width,
  pick_vector_width_shift, prevpow2, simd_integer_register_size, fma_fast,
  smax, smin, has_feature, has_opmask_registers, register_count, static_sizeof,
  cpu_name, register_size, unwrap, intlog2, nextpow2, fast_half

using SIMDTypes:
  Bit, FloatingTypes, SignedHW, UnsignedHW, IntegerTypesHW, NativeTypesExceptBitandFloat16,
  NativeTypesExceptBit, NativeTypesExceptFloat16, NativeTypes, _Vec
using LayoutPointers:
  AbstractStridedPointer, StridedPointer, StridedBitPointer, memory_reference,
  stridedpointer, zstridedpointer, similar_no_offset, similar_with_offset,
  grouped_strided_pointer, stridedpointers, bytestrides, DensePointerWrapper,
  zero_offsets

asbool(::Type{True}) = true
asbool(::Type{False}) = false
# TODO: see if `@inline` is good enough.
# @inline asvalbool(r) = Val(map(Bool, r))
# @inline asvalint(r) = Val(map(Int, r))
@generated function asvalint(r::T) where {T<:Tuple{Vararg{StaticInt}}}
    t = Expr(:tuple)
    for s ∈ T.parameters
        push!(t.args, s.parameters[1])
    end
    Expr(:call, Expr(:curly, :Val, t))
end
@generated function asvalbool(r::T) where {T<:Tuple{Vararg{StaticBool}}}
    t = Expr(:tuple)
    for b ∈ T.parameters
        push!(t.args, b === True)
    end
    Expr(:call, Expr(:curly, :Val, t))
end
@inline val_stride_rank(A) = asvalint(stride_rank(A))
@inline val_dense_dims(A) = asvalbool(ArrayInterface.dense_dims(A))

# doesn't export `Zero` and `One` by default, as these names could conflict with an AD library
export Vec, Mask, EVLMask, MM, stridedpointer, vload, vstore!, StaticInt, True, False, Bit,
    vbroadcast, mask, vfmadd, vfmsub, vfnmadd, vfnmsub,
    VecUnroll, Unroll, pick_vector_width

using Base: llvmcall, VecElement, HWReal, tail
const LLVMCALL = GlobalRef(Base, :llvmcall)

const Boolean = Union{Bit,Bool}

abstract type AbstractSIMD{W,T <: Union{<:StaticInt,NativeTypes}} <: Real end
abstract type AbstractSIMDVector{W,T} <: AbstractSIMD{W,T} end
"""
    VecUnroll{N,W,T,V<:Union{NativeTypes,AbstractSIMD{W,T}}} <: AbstractSIMD{W,T}

`VecUnroll` supports optimizations when interleaving instructions across different memory storage schemes.
`VecUnroll{N,W,T} is typically a tuple of `N+1` `AbstractSIMDVector{W,T}`s. For example, a `VecUnroll{3,8,Float32}`
is a collection of 4× `Vec{8,Float32}`.

# Examples

```jldoctest; setup=:(using VectorizationBase)
julia> rgbs = [(R = Float32(i)/255, G = Float32(i+100)/255, B = Float32(i+200)/255) for i in 0:7:49]
8-element Vector{NamedTuple{(:R, :G, :B), Tuple{Float32, Float32, Float32}}}:
 (R = 0.0, G = 0.39215687, B = 0.78431374)
 (R = 0.02745098, G = 0.41960785, B = 0.8117647)
 (R = 0.05490196, G = 0.44705883, B = 0.8392157)
 (R = 0.08235294, G = 0.4745098, B = 0.8666667)
 (R = 0.10980392, G = 0.5019608, B = 0.89411765)
 (R = 0.13725491, G = 0.5294118, B = 0.92156863)
 (R = 0.16470589, G = 0.5568628, B = 0.9490196)
 (R = 0.19215687, G = 0.58431375, B = 0.9764706)

julia> ret = vload(stridedpointer(reinterpret(reshape, Float32, rgbs)), Unroll{1,1,3,2,8,zero(UInt),1}((1,1)))
3 x Vec{8, Float32}
Vec{8, Float32}<0.0f0, 0.02745098f0, 0.05490196f0, 0.08235294f0, 0.10980392f0, 0.13725491f0, 0.16470589f0, 0.19215687f0>
Vec{8, Float32}<0.39215687f0, 0.41960785f0, 0.44705883f0, 0.4745098f0, 0.5019608f0, 0.5294118f0, 0.5568628f0, 0.58431375f0>
Vec{8, Float32}<0.78431374f0, 0.8117647f0, 0.8392157f0, 0.8666667f0, 0.89411765f0, 0.92156863f0, 0.9490196f0, 0.9764706f0>

julia> typeof(ret)
VecUnroll{2, 8, Float32, Vec{8, Float32}}
```
While the `R`, `G`, and `B` are interleaved in `rgb`s, they have effectively been split out in `ret`
(the first contains all 8 `R` values, with `G` and `B` in the second and third, respectively).

To optimize for the user's CPU, in real code it would typically be better to use `Int(pick_vector_width(Float32))`
in place of `8` (`W`) in the `Unroll` construction.
"""
struct VecUnroll{N,W,T,V<:Union{NativeTypes,AbstractSIMD{W,T}}} <: AbstractSIMD{W,T}
    data::Tuple{V,Vararg{V,N}}
    @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMD{W,T}} = new{N,W,T,V}(data)
    @inline (VecUnroll(data::Tuple{T,Vararg{T,N}})::VecUnroll{N,T,T}) where {N,T<:NativeTypes} = new{N,1,T,T}(data)
    # # following two definitions are for checking that you aren't accidentally creating `VecUnroll{0}`s.
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMD{W,T}} = (@assert(N > 0); new{N,W,T,V}(data))
    # @inline (VecUnroll(data::Tuple{T,Vararg{T,N}})::VecUnroll{N,T,T}) where {N,T<:NativeTypes} = (@assert(N > 0); new{N,1,T,T}(data))

    # @inline VecUnroll{N,W,T,V}(data::Tuple{V,Vararg{V,N}}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,Vec{W,T}}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
end
# @inline VecUnroll(data::Tuple) = VecUnroll(promote(data))
# const AbstractSIMD{W,T} = Union{AbstractSIMDVector{W,T},VecUnroll{<:Any,W,T}}
const IntegerTypes = Union{StaticInt,IntegerTypesHW}
const VecOrScalar = Union{AbstractSIMDVector,NativeTypes}
const NativeTypesV = Union{AbstractSIMD,NativeTypes,StaticInt}
# const NativeTypesV = Union{AbstractSIMD,NativeTypes,StaticInt}
const IntegerTypesV = Union{AbstractSIMD{<:Any,<:IntegerTypes},IntegerTypesHW}
struct Vec{W,T} <: AbstractSIMDVector{W,T}
    data::NTuple{W,Core.VecElement{T}}
    @inline Vec{W,T}(x::NTuple{W,Core.VecElement{T}}) where {W,T<:NativeTypes} = new{W,T}(x)
    @generated function Vec(x::Tuple{Core.VecElement{T},Vararg{Core.VecElement{T},_W}}) where {_W,T<:NativeTypes}
        W = _W + 1
        # @assert W === pick_vector_width(W, T)# || W === 8
        vtyp = Expr(:curly, :Vec, W, T)
        Expr(:block, Expr(:meta,:inline), Expr(:(::), Expr(:call, vtyp, :x), vtyp))
    end
    # @inline function Vec(x::NTuple{W,<:Core.VecElement}) where {W}
    #     T = eltype(x)
    #     @assert W === pick_vector_width(W, T)
    #     # @assert ispow2(W) && (W ≤ max(pick_vector_width(W, T), 8))
    #     new{W,T}(x)
    # end
end

@inline Base.copy(v::AbstractSIMDVector) = v
@inline asvec(x::_Vec) = Vec(x)
@inline asvec(x) = x
@inline data(vu::VecUnroll) = getfield(vu, :data)

@inline unrolleddata(x) = x
@inline unrolleddata(x::VecUnroll) = getfield(x, :data)

@inline _demoteint(::Type{T}) where {T} = T
@inline _demoteint(::Type{Int64}) = Int32
@inline _demoteint(::Type{UInt64}) = UInt32

# abstract type AbstractMask{W,U<:Union{UnsignedHW,UInt128,UInt256,UInt512,UInt1024}} <: AbstractSIMDVector{W,Bit} end
abstract type AbstractMask{W,U<:Union{UnsignedHW,UInt128}} <: AbstractSIMDVector{W,Bit} end
struct Mask{W,U} <: AbstractMask{W,U}
    u::U
    @inline function Mask{W,U}(u::Unsigned) where {W,U} # ignores U...
        U2 = mask_type(StaticInt{W}())
        new{W,U2}(u % U2)
    end
end
struct EVLMask{W,U} <: AbstractMask{W,U}
    u::U
    evl::UInt32
    @inline function EVLMask{W,U}(u::Unsigned, evl) where {W,U} # ignores U...
        U2 = mask_type(StaticInt{W}())
        new{W,U2}(u % U2, evl % UInt32)
    end
end
const AnyMask{W} = Union{AbstractMask{W},VecUnroll{<:Any,W,Bit,<:AbstractMask{W}}}
@inline Mask{W}(u::U) where {W,U<:Unsigned} = Mask{W,U}(u)
@inline EVLMask{W}(u::U, i) where {W,U<:Unsigned} = EVLMask{W,U}(u, i)
@inline Mask{1}(b::Bool) = b
@inline EVLMask{1}(b::Bool, i) = b
@inline Mask(m::EVLMask{W,U}) where {W,U} = Mask{W,U}(getfield(m,:u))
# Const prop is good enough; added an @inferred test to make sure.
# Removed because confusion can cause more harm than good.

@inline Base.broadcastable(v::AbstractSIMDVector) = Ref(v)

Vec{W,T}(x::Vararg{NativeTypes,W}) where {W,T<:NativeTypes} = Vec(ntuple(w -> Core.VecElement{T}(x[w]), Val{W}()))
Vec{1,T}(x::Union{Float32,Float64}) where {T<:NativeTypes} = T(x)
Vec{1,T}(x::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Bool}) where {T<:NativeTypes} = T(x)

@inline Base.length(::AbstractSIMDVector{W}) where W = W
@inline Base.size(::AbstractSIMDVector{W}) where W = (W,)
@inline Base.eltype(::AbstractSIMD{W,T}) where {W,T} = T
@inline Base.conj(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.adjoint(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.transpose(v::AbstractSIMDVector) = v # so that things like dot products work.

# Not using getindex/setindex as names to emphasize that these are generally treated as single objects, not collections.
@generated function extractelement(v::Vec{W,T}, i::I) where {W,I <: IntegerTypesHW,T}
    typ = LLVM_TYPES[T]
    instrs = """
        %res = extractelement <$W x $typ> %0, i$(8sizeof(I)) %1
        ret $typ %res
    """
    call = :($LLVMCALL($instrs, $T, Tuple{_Vec{$W,$T},$I}, data(v), i))
    Expr(:block, Expr(:meta, :inline), call)
end
@generated function insertelement(v::Vec{W,T}, x::T, i::I) where {W,I <: IntegerTypesHW,T}
    typ = LLVM_TYPES[T]
    instrs = """
        %res = insertelement <$W x $typ> %0, $typ %1, i$(8sizeof(I)) %2
        ret <$W x $typ> %res
    """
    call = :(Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T},$T,$I}, data(v), x, i)))
    Expr(:block, Expr(:meta, :inline), call)
end
@inline (v::AbstractSIMDVector)(i::IntegerTypesHW) = extractelement(v, i - one(i))
@inline (v::AbstractSIMDVector)(i::Integer) = extractelement(v, Int(i) - 1)
Base.@propagate_inbounds (vu::VecUnroll)(i::Integer, j::Integer) = getfield(vu, :data)[j](i)

@inline Base.Tuple(v::Vec{W}) where {W} = ntuple(v, Val{W}())

# Use with care in function signatures; try to avoid the `T` to stay clean on Test.detect_unbound_args

@inline data(v) = v
@inline data(v::Vec) = getfield(v, :data)

function Base.show(io::IO, v::AbstractSIMDVector{W,T}) where {W,T}
    name = typeof(v)
    print(io, "$(name)<")
    for w ∈ 1:W
        print(io, repr(extractelement(v, w-1)))
        w < W && print(io, ", ")
    end
    print(io, ">")
end
Base.bitstring(m::AbstractMask{W}) where {W} = bitstring(data(m))[end-W+1:end]
function Base.show(io::IO, m::AbstractMask{W}) where {W}
  bits = data(m)
  if m isa EVLMask
    print(io, "EVLMask{$W,Bit}<")
  else
    print(io, "Mask{$W,Bit}<")
  end
  for w ∈ 0:W-1
    print(io, (bits & 0x01) % Int)
    bits >>= 0x01
    w < W-1 && print(io, ", ")
  end
  print(io, ">")
end
function Base.show(io::IO, vu::VecUnroll{N,W,T,V}) where {N,W,T,V}
    println(io, "$(N+1) x $V")
    d = data(vu)
    for n in 1:N+1
        show(io, d[n]);
        n > N || println(io)
    end
end

"""
The name `MM` type refers to _MM registers such as `XMM`, `YMM`, and `ZMM`.
`MMX` from the original MMX SIMD instruction set is a [meaningless initialism](https://en.wikipedia.org/wiki/MMX_(instruction_set)#Naming).

The `MM{W,X}` type is used to represent SIMD indexes of width `W` with stride `X`.
"""
struct MM{W,X,I<:Union{HWReal,StaticInt}} <: AbstractSIMDVector{W,I}
    i::I
    @inline MM{W,X}(i::T) where {W,X,T<:Union{HWReal,StaticInt}} = new{W,X::Int,T}(i)
end
@inline MM(i::MM{W,X}) where {W,X} = MM{W,X}(getfield(i, :i))
@inline MM{W}(i::Union{HWReal,StaticInt}) where {W} = MM{W,1}(i)
@inline MM{W}(i::Union{HWReal,StaticInt}, ::StaticInt{X}) where {W,X} = MM{W,X}(i)
@inline data(i::MM) = getfield(i, :i)

@inline extractelement(i::MM{W,X,I}, j) where {W,X,I} = getfield(i, :i) + (X % I) * (j % I)

function Base.getproperty(::AbstractSIMD, ::Symbol)
    throw(ErrorException("""
`Base.getproperty` not defined on AbstractSIMD.
If you wish to work with the data as a tuple, it is recommended to use `Tuple(v)`. Once you have an ordinary tuple, you can access
individual elements normally. Alternatively, you can index using parenthesis, e.g. `v(1)` indexes the first element.
Parenthesis are used instead of `getindex`/square brackets because `AbstractSIMD` objects represent a single number, and
for `x::Number`, `x[1] === x`.

If you wish to perform a reduction on the collection, the naming convention is prepending the base function with a `v`. These functions
are not overloaded, because for `x::Number`, `sum(x) === x`. Functions include `vsum`, `vprod`, `vmaximum`, `vminimum`, `vany`, and `vall`.

If you wish to define a new operation applied to the entire vector, do not define it in terms of operations on the individual eleemnts.
This will often lead to bad code generation -- bad in terms of both performance, and often silently producing incorrect results!
Instead, implement them in terms of existing functions defined on `::AbstractSIMD`. Please feel free to file an issue if you would like
clarification, and especially if you think the function may be useful for others and should be included in `VectorizationBase.jl`.
"""))
end

"""
  pause()

For use in spin-and-wait loops, like spinlocks.
"""
@inline pause() = ccall(:jl_cpu_pause, Cvoid, ())

# notinthreadedregion() = iszero(ccall(:jl_in_threaded_region, Cint, ()))
# function assert_init_has_finished()
#     _init_has_finished[] || throw(ErrorException("bad stuff happened"))
#     return nothing
# end

include("static.jl")
include("cartesianvindex.jl")
include("early_definitions.jl")
include("promotion.jl")
include("llvm_types.jl")
include("lazymul.jl")
include("strided_pointers/stridedpointers.jl")
# include("strided_pointers/bitpointers.jl")
include("strided_pointers/cartesian_indexing.jl")
include("strided_pointers/cse_stridemultiples.jl")
include("llvm_intrin/binary_ops.jl")
include("vector_width.jl")
include("ranges.jl")
include("llvm_intrin/conversion.jl")
include("llvm_intrin/masks.jl")
include("llvm_intrin/intrin_funcs.jl")
include("llvm_intrin/memory_addr.jl")
include("llvm_intrin/unary_ops.jl")
include("llvm_intrin/vbroadcast.jl")
include("llvm_intrin/vector_ops.jl")
include("llvm_intrin/nonbroadcastingops.jl")
include("llvm_intrin/integer_fma.jl")
include("llvm_intrin/vfmaddsub.jl")
include("vecunroll/memory.jl")
include("vecunroll/fmap.jl")
include("base_defs.jl")
include("alignment.jl")
include("special/misc.jl")
include("special/double.jl")
include("special/exp.jl")
include("special/verf.jl")
# include("special/log.jl")

demoteint(::Type{T}, W) where {T} = False()
demoteint(::Type{UInt64}, W::StaticInt) = gt(W, pick_vector_width(UInt64))
demoteint(::Type{Int64}, W::StaticInt) = gt(W, pick_vector_width(Int64))

@generated function simd_vec(::DemoteInt, y::_T, x::Vararg{_T,_W}) where {DemoteInt<:StaticBool,_T,_W}
    W = 1 + _W
    T = DemoteInt === True ? _demoteint(_T) : _T
    trunc = T !== _T
    Wfull = nextpow2(W)
    ty = LLVM_TYPES[T]
    init = W == Wfull ? "undef" : "zeroinitializer"
    instrs = ["%v0 = insertelement <$Wfull x $ty> $init, $ty %0, i32 0"]
    Tup = Expr(:curly, :Tuple, T)
    for w ∈ 1:_W
        push!(instrs, "%v$w = insertelement <$Wfull x $ty> %v$(w-1), $ty %$w, i32 $w")
        push!(Tup.args, T)
    end
    push!(instrs, "ret <$Wfull x $ty> %v$_W")
    llvmc = :($LLVMCALL($(join(instrs,"\n")), _Vec{$Wfull,$T}, $Tup))
    trunc ? push!(llvmc.args, :(y % $T)) : push!(llvmc.args, :y)
    for w ∈ 1:_W
        ref = Expr(:ref, :x, w)
        trunc && (ref = Expr(:call, :%, ref, T))
        push!(llvmc.args, ref)
    end
    quote
        $(Expr(:meta,:inline))
        Vec($llvmc)
    end
end
function vec_quote(demote, W, Wpow2, offset::Int = 0)
  call = Expr(:call, :simd_vec, Expr(:call, demote ? :True : :False)); Wpow2 += offset
  iszero(offset) && push!(call.args, :y)
  foreach(w -> push!(call.args, Expr(:call, getfield, :x, w, false)), max(1,offset):min(W,Wpow2)-1)
  foreach(w -> push!(call.args, Expr(:call, :zero, :T)), W+1:Wpow2)
  call
end
@generated function _vec(::StaticInt{_Wpow2}, ::DemoteInt, y::T, x::Vararg{T,_W}) where {DemoteInt<:StaticBool, _Wpow2, _W, T <: NativeTypes}
    1+3
  W = _W + 1
  demote = DemoteInt === True
  Wpow2 = demote ? 2_Wpow2 : _Wpow2
  if W ≤ Wpow2
    vec_quote(demote, W, Wpow2)
  else
    tup = Expr(:tuple)
    offset = 0
    while offset < W
      push!(tup.args, vec_quote(demote, W, Wpow2, offset)); offset += Wpow2
    end
    Expr(:call, :VecUnroll, tup)
  end
end
@inline function Vec(y::T, x::Vararg{T,_W}) where {_W, T <: NativeTypes}
  W = StaticInt{_W}() + One()
  _vec(pick_vector_width(W, T), demoteint(T, W), y, x...)
end

@inline reduce_to_onevec(f::F, vu::VecUnroll) where {F} = ArrayInterface.reduce_tup(f, data(vu))

include("precompile.jl")
_precompile_()

end # module


text:
# K-means algorithm

#### Interface

type KmeansResult{T<:AbstractFloat} <: ClusteringResult
    centers::Matrix{T}         # cluster centers (d x k)
    assignments::Vector{Int}   # assignments (n)
    costs::Vector{T}           # costs of the resultant assignments (n)
    counts::Vector{Int}        # number of samples assigned to each cluster (k)
    cweights::Vector{Float64}  # cluster weights (k)
    totalcost::Float64         # total cost (i.e. objective) (k)
    iterations::Int            # number of elapsed iterations
    converged::Bool            # whether the procedure converged
end

const _kmeans_default_init = :kmpp
const _kmeans_default_maxiter = 100
const _kmeans_default_tol = 1.0e-6
const _kmeans_default_display = :none

function kmeans!{T<:AbstractFloat}(X::Matrix{T}, centers::Matrix{T};
                                   weights=nothing,
                                   maxiter::Integer=_kmeans_default_maxiter,
                                   tol::Real=_kmeans_default_tol,
                                   display::Symbol=_kmeans_default_display,
                                   distance::SemiMetric=SqEuclidean())

    m, n = size(X)
    m2, k = size(centers)
    m == m2 || throw(DimensionMismatch("Inconsistent array dimensions."))
    (2 <= k < n) || error("k must have 2 <= k < n.")

    assignments = zeros(Int, n)
    costs = zeros(T, n)
    counts = Vector{Int}(k)
    cweights = Vector{Float64}(k)

    _kmeans!(X, conv_weights(T, n, weights), centers,
             assignments, costs, counts, cweights,
             round(Int, maxiter), tol, display_level(display), distance)
end

function kmeans{T<:AbstractFloat}(X::Matrix{T}, k::Int;
                weights=nothing,
                init=_kmeans_default_init,
                maxiter::Integer=_kmeans_default_maxiter,
                tol::Real=_kmeans_default_tol,
                display::Symbol=_kmeans_default_display,
                distance::SemiMetric=SqEuclidean())

    m, n = size(X)
    (2 <= k < n) || error("k must have 2 <= k < n.")
    iseeds = initseeds(init, X, k)
    centers = copyseeds(X, iseeds)
    kmeans!(X, centers;
            weights=weights,
            maxiter=maxiter,
            tol=tol,
            display=display,
            distance=distance)
end

#### Core implementation

# core k-means skeleton
function _kmeans!{T<:AbstractFloat}(
    x::Matrix{T},                   # in: sample matrix (d x n)
    w::Union{Void, Vector{T}},      # in: sample weights (n)
    centers::Matrix{T},             # in/out: matrix of centers (d x k)
    assignments::Vector{Int},       # out: vector of assignments (n)
    costs::Vector{T},               # out: costs of the resultant assignments (n)
    counts::Vector{Int},            # out: the number of samples assigned to each cluster (k)
    cweights::Vector{Float64},      # out: the weights of each cluster
    maxiter::Int,                   # in: maximum number of iterations
    tol::Real,                      # in: tolerance of change at convergence
    displevel::Int,                 # in: the level of display
    distance::SemiMetric)             # in: function to calculate the distance with

    # initialize

    k = size(centers, 2)
    to_update = Vector{Bool}(k) # indicators of whether a center needs to be updated
    unused = Int[]
    num_affected::Int = k # number of centers, to which the distances need to be recomputed

    dmat = pairwise(distance, centers, x)
    dmat = convert(Array{T}, dmat) #Can be removed if one day Distance.result_type(SqEuclidean(), T, T) == T
    update_assignments!(dmat, true, assignments, costs, counts, to_update, unused)
    objv = w == nothing ? sum(costs) : dot(w, costs)

    # main loop
    t = 0
    converged = false
    if displevel >= 2
        @printf "%7s %18s %18s | %8s \n" "Iters" "objv" "objv-change" "affected"
        println("-------------------------------------------------------------")
        @printf("%7d %18.6e\n", t, objv)
    end

    while !converged && t < maxiter
        t = t + 1

        # update (affected) centers

        update_centers!(x, w, assignments, to_update, centers, cweights)

        if !isempty(unused)
            repick_unused_centers(x, costs, centers, unused)
        end

        # update pairwise distance matrix

        if !isempty(unused)
            to_update[unused] = true
        end

        if t == 1 || num_affected > 0.75 * k
            pairwise!(dmat, distance, centers, x)
        else
            # if only a small subset is affected, only compute for that subset
            affected_inds = find(to_update)
            dmat_p = pairwise(distance, centers[:, affected_inds], x)
            dmat[affected_inds, :] = dmat_p
        end

        # update assignments

        update_assignments!(dmat, false, assignments, costs, counts, to_update, unused)
        num_affected = sum(to_update) + length(unused)

        # compute change of objective and determine convergence

        prev_objv = objv
        objv = w == nothing ? sum(costs) : dot(w, costs)
        objv_change = objv - prev_objv

        if objv_change > tol
            warn("The objective value changes towards an opposite direction")
        end

        if abs(objv_change) < tol
            converged = true
        end

        # display iteration information (if asked)

        if displevel >= 2
            @printf("%7d %18.6e %18.6e | %8d\n", t, objv, objv_change, num_affected)
        end
    end

    if displevel >= 1
        if converged
            println("K-means converged with $t iterations (objv = $objv)")
        else
            println("K-means terminated without convergence after $t iterations (objv = $objv)")
        end
    end

    return KmeansResult(centers, assignments, costs, counts, cweights,
                        Float64(objv), t, converged)
end

#
#  Updates assignments, costs, and counts based on
#  an updated (squared) distance matrix
#
function update_assignments!{T<:AbstractFloat}(
    dmat::Matrix{T},            # in:  distance matrix (k x n)
    is_init::Bool,              # in:  whether it is the initial run
    assignments::Vector{Int},   # out: assignment vector (n)
    costs::Vector{T},           # out: costs of the resultant assignment (n)
    counts::Vector{Int},        # out: number of samples assigned to each cluster (k)
    to_update::Vector{Bool},    # out: whether a center needs update (k)
    unused::Vector{Int})        # out: the list of centers get no samples assigned to it

    k::Int, n::Int = size(dmat)

    # re-initialize the counting vector
    fill!(counts, 0)

    if is_init
        fill!(to_update, true)
    else
        fill!(to_update, false)
        if !isempty(unused)
            empty!(unused)
        end
    end

    # process each sample
    @inbounds for j = 1 : n

        # find the closest cluster to the i-th sample
        a::Int = 1
        c::T = dmat[1, j]
        for i = 2 : k
            ci = dmat[i, j]
            if ci < c
                a = i
                c = ci
            end
        end

        # set/update the assignment
        if is_init
            assignments[j] = a
        else  # update
            pa = assignments[j]
            if pa != a
                # if assignment changes,
                # both old and new centers need to be updated
                assignments[j] = a
                to_update[a] = true
                to_update[pa] = true
            end
        end

        # set costs and counts accordingly
        costs[j] = c
        counts[a] += 1
    end

    # look for centers that have no associated samples

    for i = 1 : k
        if counts[i] == 0
            push!(unused, i)
            to_update[i] = false # this is handled using different mechanism
        end
    end
end

#
#  Update centers based on updated assignments
#
#  (specific to the case where samples are not weighted)
#
function update_centers!{T<:AbstractFloat}(
    x::Matrix{T},                   # in: sample matrix (d x n)
    w::Void,                        # in: sample weights
    assignments::Vector{Int},       # in: assignments (n)
    to_update::Vector{Bool},        # in: whether a center needs update (k)
    centers::Matrix{T},             # out: updated centers (d x k)
    cweights::Vector)               # out: updated cluster weights (k)

    d::Int = size(x, 1)
    n::Int = size(x, 2)
    k::Int = size(centers, 2)

    # initialize center weights
    for i = 1 : k
        if to_update[i]
            cweights[i] = 0.
        end
    end

    # accumulate columns
    @inbounds for j = 1 : n
        cj = assignments[j]
        1 <= cj <= k || error("assignment out of boundary.")
        if to_update[cj]
            if cweights[cj] > 0
                for i = 1:d
                    centers[i, cj] += x[i, j]
                end
            else
                for i = 1:d
                    centers[i, cj] = x[i, j]
                end
            end
            cweights[cj] += 1
        end
    end

    # sum ==> mean
    for j = 1:k
        if to_update[j]
            @inbounds cj::T = 1 / cweights[j]
            vj = view(centers,:,j)
            for i = 1:d
                @inbounds vj[i] *= cj
            end
        end
    end
end

#
#  Update centers based on updated assignments
#
#  (specific to the case where samples are weighted)
#
function update_centers!{T<:AbstractFloat}(
    x::Matrix{T},                   # in: sample matrix (d x n)
    weights::Vector{T},             # in: sample weights (n)
    assignments::Vector{Int},       # in: assignments (n)
    to_update::Vector{Bool},        # in: whether a center needs update (k)
    centers::Matrix{T},             # out: updated centers (d x k)
    cweights::Vector)               # out: updated cluster weights (k)

    d::Int = size(x, 1)
    n::Int = size(x, 2)
    k::Int = size(centers, 2)

    # initialize center weights
    for i = 1 : k
        if to_update[i]
            cweights[i] = 0.
        end
    end

    # accumulate columns
    # accumulate_cols_u!(centers, cweights, x, assignments, weights, to_update)
    for j = 1 : n
        @inbounds wj = weights[j]

        if wj > 0
            @inbounds cj = assignments[j]
            1 <= cj <= k || error("assignment out of boundary.")

            if to_update[cj]
                rj = view(centers, :, cj)
                xj = view(x, :, j)
                if cweights[cj] > 0
                    for i = 1:d
                        @inbounds rj[i] += xj[i] * wj
                    end
                else
                    for i = 1:d
                        @inbounds rj[i] = xj[i] * wj
                    end
                end
                cweights[cj] += wj
            end
        end
    end

    # sum ==> mean
    for j = 1:k
        if to_update[j]
            @inbounds cj::T = 1 / cweights[j]
            vj = view(centers,:,j)
            for i = 1:d
                @inbounds vj[i] *= cj
            end
        end
    end
end

#
#  Re-picks centers that get no samples assigned to them.
#
function repick_unused_centers{T<:AbstractFloat}(
    x::Matrix{T},           # in: the sample set (d x n)
    costs::Vector{T},       # in: the current assignment costs (n)
    centers::Matrix{T},     # to be updated: the centers (d x k)
    unused::Vector{Int})    # in: the set of indices of centers to be updated

    # pick new centers using a scheme like kmeans++
    ds = similar(costs)
    tcosts = copy(costs)
    n = size(x, 2)

    for i in unused
        j = wsample(1:n, tcosts)
        tcosts[j] = 0
        v = x[:,j]
        centers[:,i] = v

        colwise!(ds, SqEuclidean(), v, x)
        tcosts = min(tcosts, ds)
    end
end

ext: jl


text:
"""
A type for storing a frequency histogram of MerCounts, also sometimes referred
to as a Kmer spectra.
"""
struct MerCountHist
    data::Vector{UInt64}
    min::UInt8
end

function Base.summary(io::IO, hist::MerCountHist)
    print(io, "Count histogram of motifs appearing more than ", hist.min, "times")
end

"""
    hist!(freqs::Vector{MerCount{M}}, min_count::Integer = 0) where {M<:AbstractMer}

Build a histogram of kmer frequencies, excluding any kmer counts that don't meet
`min_count`.
"""
function hist(freqs::Vector{MerCount{M}}, min_count::Integer = 0) where {M<:AbstractMer}
    hist = zeros(UInt64, 256)
    for x in freqs
        f = freq(x)
        if f ≥ min_count
            old = hist[f]
            hist[f] = old + 1
        end
    end
    return MerFreqHist(hist, convert(UInt8, min_count))
end

"""
    hist!(freqs::Vector{MerCount{M}}, min_count::Integer = 0) where {M<:AbstractMer}

A version of `hist` that builds a MerCountHist, AND mutates the input `freqs` by
filtering it to remove MerCounts that did not meet `min_count` and so were not
included in the constructed MerCountHist.
"""
function hist!(freqs::Vector{MerCount{M}}, min_count::Integer = 0) where {M<:AbstractMer}
    hist = zeros(UInt64, 256)
    wi = firstindex(freqs)
    used = 0
    for x in freqs
        f = freq(x)
        if f ≥ min_count
            old = hist[f]
            hist[f] = old + 1
            freqs[wi] = x
            wi = wi + 1
            used = used + 1
        end
        resize!(freqs, used)
    end
    return MerCountHist(hist, convert(UInt8, min_count))
end


text:
module TestRandomInitialParameters
    verbose = false

    using Base.Test

    import IntertemporalChoiceHeuristics:
        load,
        fit,
        models,
        initial_parameters,
        training_losses

    import Distributions: Uniform

    import Optim

    inputs, weights = load()

    λ = 0.0

    @testset "Identical optima found from different starting points" begin
        for model in models
            Θ = initial_parameters(model)
            K = length(Θ)
            for loss in training_losses
                res = fit(model, inputs, weights, loss, λ, Θ)

                fₓ = Optim.minimum(res)

                n_points = 10
                costs = Array(Float64, n_points)
                for i in 1:n_points
                    ϵ = rand(Uniform(-0.1, +0.1), K)
                    Θ′ = Θ + ϵ
                    try
                        res = fit(model, inputs, weights, loss, λ, Θ′)
                        fₓ′ = Optim.minimum(res)
                        costs[i] = fₓ′
                        δ = (fₓ′ - fₓ) / fₓ
                        if verbose
                            @printf(
                                "%s\t%s\t%d\t%.4f\t%.4f\t%.4f\n",
                                model,
                                loss,
                                i,
                                fₓ,
                                fₓ′,
                                δ,
                            )
                        end
                    catch
                        costs[i] = fₓ
                    end
                end
                # TODO: Get this to pass reliably.
                # @test all(x -> x < 1e-10, (costs ./ minimum(costs)) - 1.0)
            end
        end
    end
end


text:
@doc raw"""
`DAE`: Differential Algebraic Equation

Defines a differential algebraic initial value problem
```math
\begin{aligned}
\dot{q} (t) &= v(t, q(t)) + u(t, q(t), \lambda(t)) , & q(t_{0}) &= q_{0} , \\
0 &= \phi (t, q(t), \lambda(t)) , & \lambda(t_{0}) &= \lambda_{0} ,
\end{aligned}
```
with vector field ``v``, projection ``u``, algebraic constraint ``\phi=0``,
initial conditions ``q_{0}`` and ``\lambda_{0}``, the dynamical variable ``q``
taking values in ``\mathbb{R}^{d}`` and the algebraic variable ``\lambda``
taking values in ``\mathbb{R}^{m}``.

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `vType <: Function`: type of `v`
* `uType <: Function`: type of `u`
* `ūType <: OptionalFunction`: type of `ū`
* `ϕType <: Function`: type of `ϕ`
* `ψType <: OptionalFunction`: type of `ψ`
* `v̄Type <: Function`: type of `v̄`
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `d`: dimension of dynamical variable ``q`` and the vector field ``v``
* `m`: dimension of algebraic variable ``\lambda`` and the constraint ``\phi``
* `v`: function computing the vector field
* `u`: function computing the projection
* `ū`: function computing the secondary projection field ``\bar{u}`` (optional)
* `ϕ`: algebraic constraint
* `ψ`: secondary constraints (optional)
* `v̄`: function computing an initial guess for the velocity field ``v`` (defaults to `v`)
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `λ₀`: initial condition for algebraic variable ``\lambda``
* `μ₀`: initial condition for algebraic variable ``μ`` (optional)
* `invariants`: either a `NamedTuple` containing the equation's invariants or `nothing`
* `parameters`: either a `NamedTuple` containing the equation's parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

The function `v`, providing the vector field, takes three arguments,
`v(t, q, v)`, the functions `u` and `ϕ`, providing the projection and the
algebraic constraint take four arguments, `u(t, q, λ, u)` and `ϕ(t, q, λ, ϕ)`,
where `t` is the current time, `q` and `λ` are the current solution vectors,
and `v`, `u` and `ϕ` are the vectors which hold the result of evaluating the
vector field ``v``, the projection ``u`` and the algebraic constraint ``\phi``
on `t`, `q` and `λ`.

### Constructors

```julia
DAE(v, u, ū, ϕ, ψ, v̄, t₀, q₀, λ₀, invariants, parameters, periodicity)

DAE(v, u, ϕ, t₀, q₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
DAE(v, u, ϕ, q₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
DAE(v, u, ϕ, t₀, q₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
DAE(v, u, ϕ, q₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)

DAE(v, u, ū, ϕ, ψ, t₀, q₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
DAE(v, u, ū, ϕ, ψ, q₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
DAE(v, u, ū, ϕ, ψ, t₀, q₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
DAE(v, u, ū, ϕ, ψ, q₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
```

### Example

```julia
    function v(t, q, v)
        v[1] = q[1]
        v[2] = q[2]
    end

    function u(t, q, λ, u)
        u[1] = +λ[1]
        u[2] = -λ[1]
    end

    function ϕ(t, q, λ, ϕ)
        ϕ[1] = q[2] - q[1]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ₀ = [0.]

    dae = DAE(v, u, ϕ, t₀, q₀, λ₀)

```

"""
struct DAE{vType <: Callable,
           uType <: Callable, ϕType <: Callable,
           ūType <: OptionalCallable, ψType <: OptionalCallable,
           v̄Type <: Callable,
           invType <: OptionalInvariants,
           parType <: OptionalParameters,
           perType <: OptionalPeriodicity} <: AbstractEquationDAE{invType,parType,perType,ψType}

    v::vType
    u::uType
    ϕ::ϕType
    ū::ūType
    ψ::ψType
    v̄::v̄Type

    invariants::invType
    parameters::parType
    periodicity::perType

    function DAE(v, u, ϕ, ū, ψ, v̄, invariants, parameters, periodicity)
        @assert !isempty(methods(v))
        @assert !isempty(methods(u))
        @assert !isempty(methods(ϕ))
        @assert !isempty(methods(ū)) || ū === nothing
        @assert !isempty(methods(ψ)) || ψ === nothing
        @assert !isempty(methods(v̄))

        new{typeof(v), typeof(u), typeof(ϕ), typeof(ū), typeof(ψ), typeof(v̄),
            typeof(invariants), typeof(parameters), typeof(periodicity)}(
                v, u, ϕ, ū, ψ, v̄, invariants, parameters, periodicity)
    end
end

DAE(v, u, ϕ, ū, ψ; v̄=v, invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = DAE(v, u, ϕ, ū, ψ, v̄, invariants, parameters, periodicity)
DAE(v, u, ϕ; kwargs...) = DAE(v, u, ϕ, nothing, nothing; kwargs...)

GeometricBase.invariants(equation::DAE) = equation.invariants
GeometricBase.parameters(equation::DAE) = equation.parameters
GeometricBase.periodicity(equation::DAE) = equation.periodicity

hasvectorfield(::DAE) = true

function check_initial_conditions(equ::DAE, ics::NamedTuple)
    haskey(ics, :q) || return false
    haskey(ics, :λ) || return false
    if hassecondary(equ)
        haskey(ics, :μ) || return false
        eltype(ics.λ) == eltype(ics.μ) || return false
        typeof(ics.λ) == typeof(ics.μ) || return false
        axes(ics.λ) == axes(ics.μ) || return false
    end
    return true
end

function check_methods(equ::DAE, tspan, ics::NamedTuple, params)
    applicable(equ.v, tspan[begin], ics.q, zero(ics.q), params) || return false
    applicable(equ.u, tspan[begin], ics.q, ics.λ, zero(ics.q), params) || return false
    applicable(equ.ϕ, tspan[begin], ics.q, zero(ics.λ), params) || return false
    return true
end

function datatype(equ::DAE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return eltype(ics.q)
end

function arrtype(equ::DAE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return typeof(ics.q)
end

_get_v(equ::DAE, params) = (t,q,v)   -> equ.v(t, q, v, params)
_get_u(equ::DAE, params) = (t,q,λ,u) -> equ.u(t, q, λ, u, params)
_get_ϕ(equ::DAE, params) = (t,q,ϕ)   -> equ.ϕ(t, q, ϕ, params)
_get_ū(equ::DAE, params) = (t,q,λ,u) -> equ.ū(t, q, λ, u, params)
_get_ψ(equ::DAE, params) = (t,q,v,ψ) -> equ.ψ(t, q, v, ϕ, params)
_get_v̄(equ::DAE, params) = (t,q,v)   -> equ.v̄(t, q, v, params)

function _functions(equ::DAE)
    if hassecondary(equ)
        (v = equ.v, u = equ.u, ϕ = equ.ϕ, ū = equ.ū, ψ = equ.ψ, v̄ = equ.v̄)
    else
        (v = equ.v, u = equ.u, ϕ = equ.ϕ, v̄ = equ.v̄)
    end
end

function _functions(equ::DAE, params::OptionalParameters)
    if hassecondary(equ)
        (v = _get_v(equ, params), u = _get_u(equ, params), ϕ = _get_ϕ(equ, params), ū = _get_ū(equ, params), ψ = _get_ψ(equ, params), v̄ = _get_v̄(equ, params))
    else
        (v = _get_v(equ, params), u = _get_u(equ, params), ϕ = _get_ϕ(equ, params), v̄ = _get_v̄(equ, params))
    end
end


text:
using Documenter
using StridedTensorXD

makedocs(;
    modules=[StridedTensorXD],
    sitename="StridedTensorXD.jl",
    authors = "PhysicsCodesLab",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax()
    ),
    pages = [
        "Home" => "index.md",
         "Index" => ["index/index.md"]
    ],
)

deploydocs(;
   repo = "github.com/PhysicsCodesLab/StridedTensorXD.jl.git",
   devbranch="main"
)


text:
@test DiscretizeQuantile(2) == DiscretizeQuantile(2, false)

# NOTE: Test cases are generated by sklearn.preprocessing.KBinsDiscretizer
@test array_matches(binedges(DiscretizeQuantile(1), [1.0,2.0,3.0,4.0,5.0,6.0]), [1.0, 6.0], 1e-8)
@test array_matches(binedges(DiscretizeQuantile(2), [1.0,2.0,3.0,4.0,5.0,6.0]), [1.0, 3.5, 6.0], 1e-8)
@test array_matches(binedges(DiscretizeQuantile(2), [6.0,2.0,30.0,4.0,1.0,1.0]), [1.0, 3.0, 30.0], 1e-8)
@test array_matches(binedges(DiscretizeQuantile(2), [1.0,2.0,3.0,4.0,5.0]), [1.0, 3.0, 5.0], 1e-8)
@test array_matches(binedges(DiscretizeQuantile(3), [1.0,2.0,3.0,4.0,5.0]), [1.0, 2+1/3, 3+2/3, 5.0], 1e-8)
@test_throws ErrorException binedges(DiscretizeQuantile(3, false), [1.0, 1.0])
@test_throws ErrorException binedges(DiscretizeQuantile(3), [1.0,1.0,1.0,1.0,1.0,1.0])

@test array_matches(binedges(DiscretizeQuantile(2, true), [1,2]), [1.0, 1.5, 2.0])
@test array_matches(binedges(DiscretizeQuantile(3, true), [1,2,2,2,3]), [1.0, 2.0, 3.0])
@test array_matches(binedges(DiscretizeQuantile(3, true), [1,2,2,2,50]), [1.0, 2.0, 50.0])
@test array_matches(binedges(DiscretizeQuantile(3, true), [1,2,2,2,51]), [1.0, 2.0, 51.0])
@test array_matches(binedges(DiscretizeQuantile(3, true), [1,1,1,2,3,4,50,51,52]), [1.0, 1+2/3, 19+1/3, 52.0], 1e-8)
@test_throws ErrorException binedges(DiscretizeQuantile(3, false), [1,2,2,2,3])
@test_throws ErrorException binedges(DiscretizeQuantile(3), [2,2,2,2,2])


text:
# Generate some test data
t = collect(0:0.1:5.0)
f_(t) = [sin(t); cos(t)]
df_(t) = [cos(t); -sin(t)]
X = hcat(map(i -> f_(i), t)...)
Y = randn(2,2)*X
DX = hcat(map(i -> df_(i), t)...)
p = [2.7]
u_(u, p, t) = [p[1] * u[1] - 3.0 * t]
U = hcat(map(i -> u_(X[:, i], p, t[i]), 1:length(t))...)

@testset "DiscreteProblem" begin
    p1 = DiscreteDataDrivenProblem(X)
    p2 = DiscreteDataDrivenProblem(X, t)
    p3 = DiscreteDataDrivenProblem(X, t, U)
    p4 = DiscreteDataDrivenProblem(X, t, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_discrete(p_)
        @test_nowarn @is_applicable p_
    end
    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
end

@testset "Collocation" begin
    t_ = 0.0:0.001:10.0
    x_ = permutedims(sin.(t_))
    y_ = permutedims(cos.(t_))
    for m in [EpanechnikovKernel, UniformKernel,TriangularKernel,QuarticKernel,
        TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel,
        LogisticKernel, SigmoidKernel, SilvermanKernel]
        p_ = ContinuousDataDrivenProblem(x_, t_)
        @test norm(y_ .- p_.DX) < 1e-1
    end
end

@testset "ContinuousProblem" begin
    p2 = ContinuousDataDrivenProblem(X, t, GaussianKernel())
    p3 = ContinuousDataDrivenProblem(X, DX)
    p4 = ContinuousDataDrivenProblem(X, t, U)
    p5 = ContinuousDataDrivenProblem(X, t, DX, U)
    p6 = ContinuousDataDrivenProblem(X, t, DX, u_, p = p)

    for p_ in [p2;p3;p4;p5;p6]
        @test is_valid(p_)
        @test is_continuous(p_)
    end

    @test isapprox(p3.X, p2.X, atol = 0.3)
    @test isapprox(p3.DX, p2.DX, atol = 0.3)
    @test isequal(p5.U, p6.U)
    @test is_parametrized(p6)
    @test isequal(p6.p, p)
    @test has_timepoints(p4)
end

@testset "DirectProblem" begin
    p1 = DirectDataDrivenProblem(X, Y)
    p2 = DirectDataDrivenProblem(X, t, Y)
    p3 = DirectDataDrivenProblem(X, t, Y, U)
    p4 = DirectDataDrivenProblem(X, t, Y, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_direct(p_)
    end

    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
    @test !is_autonomous(p3)
end

@testset "Problem Basis Interaction" begin
    @variables x y z t α β u
    b1 = Basis([α*x; β*y; z*t^2+u], [x; y; z], iv = t, parameters = [α; β], controls = [u])
    b2 = Basis([α*x; β*y; z*t; α], [x; y; z], iv = t, parameters = [α; β])
    sample_size = 100
    X1 = randn(3, sample_size)
    DX1 = randn(3, sample_size)
    DX2 = randn(3, sample_size-1)
    X2 = randn(4, sample_size)
    Y1 = randn(3, sample_size)
    Y2 = randn(4, sample_size)
    t = randn(sample_size)
    U = randn(1, sample_size)
    ps = randn(3)

    p1 = DirectDataDrivenProblem(X1, Y1, t = t, p = ps, U = U)
    p2 = DirectDataDrivenProblem(X1, Y2, t = t, p = ps, U = U)
    p3 = ContinuousDataDrivenProblem(X1, t, DX1,  Y = Y1, p = ps, U = U)
    p4 = DiscreteDataDrivenProblem(X1, Y = Y2, t = t, p = ps, U = U)

    @testset "Check validity" begin
        @test_throws AssertionError @is_applicable p2 b2
        @test_throws AssertionError @is_applicable p1 b2
        @test_throws AssertionError @is_applicable p3 b2
        @test_throws AssertionError @is_applicable p4 b2
        @test_throws AssertionError @is_applicable p4 b1 DX1

        @test_nowarn @is_applicable p1 b1
        @test_nowarn @is_applicable p3 b1
        @test_nowarn @is_applicable p4 b1
        @test_nowarn @is_applicable p1 b1 DX1
        @test_nowarn @is_applicable p3 b1 DX1
    end

    @testset "Apply problem" begin
        DX1 .= 0.0
        @test_nowarn b1(DX1, p1)
        @test iszero(norm(DX1 - b1(p1)))

        DX1 .= 0.0
        @test_nowarn b1(DX1, p3)
        @test iszero(norm(DX1 - b1(p3)))

        DX2 .= 0.0
        @test_nowarn b1(DX2, p4)
        @test iszero(norm(DX2 - b1(p4)))

        DX1 .= 0.0
        @test_nowarn @views b1(DX1[:,1:end-1], p4)
        @test iszero(norm(DX1[:, 1:end-1] - b1(p4)))
    end
end


text:
using LShapedSolvers
using Base.Test
using JuMP
using StochasticPrograms
using Gurobi

solver = GurobiSolver(OutputFlag=0)
lsolvers = [(LShapedSolver(:ls,solver,log=false),"L-Shaped"),
            (LShapedSolver(:rd,solver,crash=Crash.EVP(),autotune=true,log=false),"RD L-Shaped"),
            (LShapedSolver(:tr,solver,crash=Crash.EVP(),autotune=true,log=false),"TR L-Shaped"),
            (LShapedSolver(:lv,solver,crash=Crash.EVP(),λ=0.6,log=false),"Leveled L-Shaped")]

problems = Vector{Tuple{JuMP.Model,String}}()
info("Loading test problems...")
info("Loading simple...")
include("simple.jl")
info("Loading farmer...")
include("farmer.jl")

info("Test problems loaded. Starting test sequence.")
@testset "$lsname Solver: $name" for (lsolver,lsname) in lsolvers, (sp,name) in problems
    solve(sp,solver=solver)
    x̄ = copy(sp.colVal)
    Q̄ = copy(sp.objVal)
    solve(sp,solver=lsolver)
    @test norm(optimal_decision(sp) - x̄) <= 1e-2
    @test abs(optimal_value(sp) - Q̄) <= 1e-2
end

info("Starting distributed tests...")

include("/opt/julia-0.6/share/julia/test/testenv.jl")
push!(test_exeflags.exec,"--color=yes")
cmd = `$test_exename $test_exeflags run_dtests.jl`

if !success(pipeline(cmd; stdout=STDOUT, stderr=STDERR)) && ccall(:jl_running_on_valgrind,Cint,()) == 0
    error("Distributed test failed, cmd : $cmd")
end


text:
#=
This file benchmarks the parallel and serial computation of a recurrence matrix
=#
using DynamicalSystemsBase, RecurrenceAnalysis
using Statistics, BenchmarkTools, Base.Threads

function pretty_time(b::BenchmarkTools.Trial)
    med = median(b)

    return string(
        BenchmarkTools.prettytime(BenchmarkTools.time(med)),
        " (",
        BenchmarkTools.prettypercent(BenchmarkTools.gcratio(med)),
        " GC)"
    )

end

Ns = round.(Int, 10.0 .^ (2:0.5:4.5))
ro = Systems.roessler()

println("I am using $(nthreads()) threads.")
for N in Ns
    # set up datasets
    tr = trajectory(ro, N*0.1; dt = 0.1, Tr = 10.0)
    printstyled("For N = $(length(tr))..."; color = :blue, bold = true)
    println()
    x = tr[:, 1]

    # calculate only the lower triangle
    b_serial   = @benchmark RecurrenceMatrix($(tr), 5.0; metric = Euclidean(), parallel = false)
    b_parallel = @benchmark RecurrenceMatrix($(tr), 5.0; metric = Euclidean(), parallel = true)
    b_s_1 = @benchmark RecurrenceMatrix($(x), 2.0; parallel = false)
    b_p_1 = @benchmark RecurrenceMatrix($(x), 2.0; parallel = true)

    # calculate triangular metrics
    threeD = median(b_serial.times)/median(b_parallel.times)
    oneD = median(b_s_1.times)/median(b_p_1.times)
    threeD, oneD = round.((threeD, oneD); sigdigits = 3)
    printstyled("Lower triangle only"; color = :green)
    println("\nSpeedups: \n    3D=$threeD\n    1D=$oneD")
    println(
        """
        Raw timings:
            3D serial:   $(pretty_time(b_serial))
            3D parallel: $(pretty_time(b_parallel))
            1D serial:   $(pretty_time(b_s_1))
            1D parallel: $(pretty_time(b_p_1))
        """
    )

    # calculate the full matrix
    b_serial_full   = @benchmark CrossRecurrenceMatrix($(tr), $(tr), 5.0; metric = Euclidean(), parallel = false)
    b_parallel_full = @benchmark CrossRecurrenceMatrix($(tr), $(tr), 5.0; metric = Euclidean(), parallel = true)
    b_s_1_full = @benchmark CrossRecurrenceMatrix($(x), $(x), 2.0; parallel = false)
    b_p_1_full = @benchmark CrossRecurrenceMatrix($(x), $(x), 2.0; parallel = true)

    # calculate full metrics
    threeD = median(b_serial_full.times)/median(b_parallel_full.times)
    oneD = median(b_s_1_full.times)/median(b_p_1_full.times)
    threeD, oneD = round.((threeD, oneD); sigdigits = 3)
    printstyled("Full matrix"; color = :green)
    println("\nSpeedups: \n    3D=$threeD\n    1D=$oneD")
    println(
        """
        Raw timings:
            3D serial:   $(pretty_time(b_serial_full))
            3D parallel: $(pretty_time(b_parallel_full))
            1D serial:   $(pretty_time(b_s_1_full))
            1D parallel: $(pretty_time(b_p_1_full))
        """
    )
    println()
end


text:
######################################################################
# Pia UTs
# -----
# Licensed under Apache License 2.0
module TestPia
using Test
using ExtraFun
using Pia
import Base.Threads: Atomic, atomic_add!
Pia.set_default_timeout(1)

function app(port::Integer)
    inst = application("TestApp", v"0.1.0", transport=PiaLink(ip"::1", port))
    Pia.bootstrap(inst, ip"::1")
    return inst
end

@testset "Pia" begin
    @testset "Address" begin
        @test Address{40}(tuple((0x1:0x28)...)) isa Address{40}
        @test Address("test") === Address((0x74, 0x65, 0x73, 0x74))
        @test_throws TypeError Address{40}("test")
    end

    @testset "Link" begin
        @testset "connect" begin
            @with app1 = app(Pia.DEFAULT_PORT) app2 = app(Pia.DEFAULT_PORT+1) begin
                @assert app1.protocols.transport isa PiaLink && app2.protocols.transport isa PiaLink

                @test_throws ConnectionError Pia.connect(app1)
                @test Pia.connect(app2) !== nothing
                sleep(0.01) # Mandatory sleep to allow app1 thread to catch up
                let peers1 = peers(app1), peers2 = peers(app2)
                    @test haskey(peers1, endpoint(app2)) && haskey(peers2, endpoint(app1))
                end
            end
        end

        @testset "connectionless" begin
            @with app1 = app(Pia.DEFAULT_PORT) app2 = app(randport(40000)) begin
                @assert app1.protocols.transport isa PiaLink && app2.protocols.transport isa PiaLink
                Pia.connect(app2)
                sleep(0.01)
                @assert Pia.isconnected(app1) && Pia.isconnected(app2)

                @test begin
                    let success = false
                        task1 = @async success = unsafe_string(pointer(read(app2, endpoint(app1))::Vector{UInt8})) == "foobar"
                        task2 = @async write(app1, endpoint(app2), "foobar")
                        wait(task2); wait(task1)
                        success
                    end
                end
            end
        end
    end
end

end


text:
using Distributions
using Bijectors
using Transits: Kipping13Transform
using HypothesisTests

@testset "Kipping13" begin
    N = 10000
    @test length(Kipping13()) == 2
    samps = rand(rng, Kipping13(), N)
    u1, u2 = eachrow(samps)

    @test all(@.(u1 + u2 < 1)) # flux is always positive
    @test all(@.(u1 + 2 * u2 > 0)) # flux is strictly monotonic increasing
    @test all(@.(u1 > 0)) # u1 always positive (flux is always positive)
    @test loglikelihood(Kipping13(), samps) == 0.0

    # make sure qs are Uniform
    q1 = @. (u1 + u2)^2
    q2 = @. 0.5 * u1 / (u1 + u2)
    for q in (q1, q2)
        t = ApproximateOneSampleKSTest(q, Uniform(0, 1))
        @test t.δ < 0.05
    end

end

@testset "Kipping13Transform" begin
    N = 1000
    @test bijector(Kipping13()) isa Kipping13Transform

    𝔅 = Kipping13Transform()
    xs = rand(Kipping13(), N)
    for x in eachcol(xs)
        y = 𝔅(x)
        xi = inv(𝔅)(y)
        @test xi ≈ x
    end

    ys = 2 .* randn(2, N)
    @test all(eachcol(ys)) do y
        u1, u2 = inv(𝔅)(y)

        u1 + u2 ≤ 1 &&
        u1 + 2 * u2 ≥ 0 &&
        u1 ≥ 0
    end
end


text:
module Keymods

export Keymod

@enum Keymod::UInt16 begin
  NONE = 0
  LSHIFT = 1
  RSHIFT = 2
  LCTRL = 64
  RCTRL = 128
  LALT = 256
  RALT = 512
  LMETA = 1024
  RMETA = 2048
  NUM = 4096
  CAPS = 8192
  MODE = 16384
  CTRL = 192
  SHIFT = 3
  ALT = 768
  META = 3072
end

end

using .Keymods


text:
@testset "462.minimum-moves-to-equal-array-elements-ii.jl" begin
    @test min_moves2([-9, -17, -18, 12, -22, 3, -23, -15, 27, 23, 14, -16, 7, -7, -24, 22, -6, -10, -6, 2, -19, 18, 7, 16, 12, 9, -28, -25, 2, 4, 27, -19, -3, 21, 10, 14, -8, -11, 14, -17, -7, -22, -3, -25, 18, -29, -24, -13, 3, -4]) == 695
    @test min_moves2([1, 2, 3]) == 2
end


text:
# Gurobi model

#################################################
#
#  model type & constructors
#
#################################################

type Model
    env::Env
    ptr_model::Ptr{Void}
    callback::Any
    finalize_env::Bool

    function Model(env::Env, p::Ptr{Void}; finalize_env::Bool=false)
        model = new(env, p, nothing, finalize_env)
        finalizer(model, m -> (free_model(m); if m.finalize_env; free_env(m.env); end))
        model
    end
end

# If environment is tied to this model, the, we can safely finalize
# both at the same time, working around the Julia GC.
function Model(env::Env, name::String; finalize_env::Bool=false)

    @assert is_valid(env)
    @assert isascii(name)

    a = Array(Ptr{Void}, 1)
    ret = @grb_ccall(newmodel, Cint, (
        Ptr{Void},  # env
        Ptr{Ptr{Void}},  # pointer to model pointer
        Ptr{UInt8},  # name
        Cint,  # numvars
        Ptr{Float64}, # obj coeffs
        Ptr{Float64}, # lbounds
        Ptr{Float64}, # ubounds
        Ptr{UInt8},   # var types,
        Ptr{Ptr{UInt8}} # varnames
        ),
        env, a, name, 0,
        C_NULL, C_NULL, C_NULL, C_NULL, C_NULL
    )

    if ret != 0
        throw(GurobiError(env, ret))
    end

    Model(env, a[1]; finalize_env=finalize_env)
end

function Model(env::Env, name::String, sense::Symbol)
    @assert isascii(name)
    model = Model(env, name)
    if sense != :minimize
        set_sense!(model, sense)
    end
    model
end

#################################################
#
#  model manipulation
#
#################################################

Base.unsafe_convert(ty::Type{Ptr{Void}}, model::Model) = model.ptr_model::Ptr{Void}

function free_model(model::Model)
    if model.ptr_model != C_NULL
        @grb_ccall(freemodel, Void, (Ptr{Void},), model.ptr_model)
        model.ptr_model = C_NULL
    end
end

function copy(model::Model)
    pm::Ptr{Void} = C_NULL
    if model.ptr_model != C_NULL
        pm = @grb_ccall(copymodel, Ptr{Void}, (Ptr{Void},), model.ptr_model)
        if pm == C_NULL
            error("Failed to copy a Gurobi model.")
        end
    end
    Model(model.env, pm)
end

function update_model!(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(updatemodel, Cint, (Ptr{Void},), model.ptr_model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function reset_model!(model::Model)
    @assert model.ptr_model != C_NULL
    ret = @grb_ccall(resetmodel, Cint, (Ptr{Void},), model.ptr_model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# read / write file

function read_model(model::Model, filename::String)
    @assert isascii(filename) # TODO: support non-ascii file names
    @assert is_valid(model.env)
    a = Array(Ptr{Void}, 1)
    ret = @grb_ccall(readmodel, Cint,
        (Ptr{Void}, Ptr{UInt8}, Ptr{Ptr{Void}}),
        model.env, filename, a)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    model.ptr_model = a[1]
    nothing
end

function write_model(model::Model, filename::String)
    @assert isascii(filename) # TODO: support non-ascii file names
    ret = @grb_ccall(write, Cint, (Ptr{Void}, Ptr{UInt8}),
        model.ptr_model, filename)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function tune_model(model::Model)
    ret = @grb_ccall(tunemodel, Cint, (Ptr{Void},), model.ptr_model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function get_tune_result!(model::Model,i::Int)
    ret = @grb_ccall(gettuneresult, Cint, (Ptr{Void}, Cint), model.ptr_model, i)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

terminate(model::Model) = @grb_ccall(terminate, Void, (Ptr{Void},), model.ptr_model)

# Presolve the model but don't solve. For some reason this is not
# documented for the C interface, but is for all the other interfaces.
# Source: https://twitter.com/iaindunning/status/519620465992556544
function presolve_model(model::Model)
    ret = @grb_ccall(presolvemodel, Ptr{Void}, (Ptr{Void},), model.ptr_model)
    if ret == C_NULL
        # Presumably failed to return a model
        error("presolve_model did not return a model")
    end
    return Model(model.env, ret)
end

function fixed_model(model::Model)
    @assert model.ptr_model != C_NULL
    fixed::Ptr{Void} = C_NULL
    fixed = @grb_ccall(fixedmodel, Ptr{Void}, (Ptr{Void},), model.ptr_model)
    if fixed == C_NULL
        error("Unable to create fixed model")
    end
    Model(model.env, fixed)
end


text:
@testset "6.4 Hyperbolic cotangent" begin
    include("6.4.1 (c+d x)^m (a+b coth)^n.jl")
    include("6.4.2 Hyperbolic cotangent functions.jl")
    include("6.4.7 (d hyper)^m (a+b (c coth)^n)^p.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/.jl


text:
using SafeTestsets

@safetestset "Rotations" begin include("rotations.jl") end
@safetestset "Ellipsoid" begin include("ellipsoid.jl") end
@safetestset "Atmosphere" begin include("atmosphere.jl") end
@safetestset "Kinematics" begin include("kinematics.jl") end
@safetestset "Mechanics" begin include("mechanics.jl") end
@safetestset "Anemometry" begin include("anemometry.jl") end


text:
#push!(Base.LOAD_PATH,joinpath(dirname(@__FILE__),"..",".."))

using Xpress, Base.Test, MathOptInterface, MathOptInterface.Test, MathOptInterfaceXpress

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIXPR = MathOptInterfaceXpress

@testset "MathOptInterfaceXpress" begin
    @testset "Unit Tests" begin
        config = MOIT.TestConfig()
        solver = XpressOptimizer()

        MOIT.basic_constraint_tests(solver, config;
            exclude = [
            ]
        )

        MOIT.unittest(solver, config, [
            "solve_affine_interval"
        ])
    end
    @testset "Linear tests" begin
        linconfig = MOIT.TestConfig()
        solver = XpressOptimizer()
        MOIT.contlineartest(solver , linconfig, ["linear10","linear12","linear8a","linear8b","linear8c"])

        solver_nopresolve = XpressOptimizer(PRESOLVE = 0)
        MOIT.contlineartest(solver_nopresolve , linconfig, ["linear10","linear12"])

        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        MOIT.linear12test(solver, linconfig_nocertificate)

        # Intervals
        linconfig_noquery = MOIT.TestConfig(query=false)
        MOIT.linear10test(solver, linconfig_noquery)
    end

    @testset "Quadratic tests" begin
        quadconfig = MOIT.TestConfig(atol=1e-5, rtol=1e-5, duals=false, query=false)
        solver = XpressOptimizer()
        MOIT.contquadratictest(solver, quadconfig)
    end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig()
        solver = XpressOptimizer()
        MOIT.lintest(solver, linconfig, ["lin3","lin4"])

        solver_nopresolve = XpressOptimizer(PRESOLVE=0)
        MOIT.lintest(solver_nopresolve, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig()
        solver = XpressOptimizer()
        MOIT.intlineartest(solver, intconfig)
    end

    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig()
        solver = XpressOptimizer()
        MOIT.nametest(solver)
        @testset "validtest" begin
            MOIT.validtest(solver)
        end
        @testset "emptytest" begin
            MOIT.emptytest(solver)
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(solver)
        end
        @testset "canaddconstrainttest" begin
            MOIT.canaddconstrainttest(solver, Float64, Complex{Float64})
        end
        @testset "copytest" begin
            solver2 = XpressOptimizer()
            MOIT.copytest(solver,solver2)
        end
    end
end
;


text:
function parareal_serial!(cache::PararealCache, solution::PararealSolution, problem::AbstractInitialValueProblem, parareal::Parareal)
    @↓ U, T, F, G = cache
    @↓ errors, saveiterates = solution
    @↓ u0, (t0, tN) ← tspan = problem
    @↓ finesolver, coarsolver, P, K = parareal
    @↓ weights, ψ, ϵ = parareal.control
    T[1] = t0
    for n = 1:P
        # stable sum
        T[n+1] = (1-n) * t0 / P + n * tN / P
    end
    F[1] = G[1] = U[1] = u0
    # coarse guess
    for n = 1:P
        chunkproblem = subproblemof(problem, U[n], T[n], T[n+1])
        chunksolution = coarsolver(chunkproblem)
        U[n+1] = chunksolution.u[end]
    end
    G .= U
    for k in 1:K
        if saveiterates
            for n = 1:(k - 1)
                solution[k][n] = solution[k-1][n]
            end
        end
        # fine run (parallelisable)
        for n = k:P
            chunkproblem = subproblemof(problem, U[n], T[n], T[n+1])
            chunksolution = finesolver(chunkproblem)
            F[n + 1] = chunksolution.u[end]
            if saveiterates
                solution[k][n] = chunksolution
            else
                solution[n] = chunksolution
            end
        end
        # check convergence
        update!(weights, U, F)
        errors[k] = ψ(cache, solution, k, weights)
        if errors[k] ≤ ϵ
            resize!(errors, k)
            resize!(solution, k)
            break
        end
        # serial update
        for n = k:P
            chunkproblem = subproblemof(problem, U[n], T[n], T[n+1])
            chunksolution = coarsolver(chunkproblem)
            U[n+1] = chunksolution.u[end] + F[n+1] - G[n+1]
            G[n+1] = chunksolution.u[end]
        end
    end
    return solution
end


text:
module test_solidProperties

using Modia3D
using Modia3D.ModiaMath
using Modia3D.Unitful
using Modia3D.Test
using Modia3D.LinearAlgebra

filename = joinpath(Modia3D.path, "objects", "primitiveFiles", "cube.obj")

# Properties
vmat1 = Modia3D.Material(color="LightBlue", transparency=0.5, wireframe=false)
vmat2 = deepcopy(vmat1)
vmat2.transparency = 0.7
smat1 = Modia3D.MassProperties(m=1.0)
font  = Modia3D.Font(fontFamily="TimesNewRoman", color="Black", bold=false, charSize=0.2)
rtext = [0.0, 0.0, 1.0]
ltext = [0.0, 0.0, -0.5]

# Dimensions
LxBox = 0.9; LyBox = 0.5; LzBox = 0.3
DxSphere = 0.7
DxEll = 0.7; DyEll = 0.3; DzEll = 0.4
DxCyl = 0.5; DyCyl = 0.6; LzCyl = 0.8
DxCone = 0.3; DyCone = 0.4; LzCone = 0.7
DxCap = 0.4; DyCap = 0.3; LzCap = 0.45
DxPipe = 0.5; DyPipe = 0.8; LzPipe = 0.4
LxBeam = 0.4; DyBeam = 0.5; LzBeam = 0.3

# Solid objects
@assembly Solids begin
   world     = Modia3D.Object3D(visualizeFrame=false)
   box       = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(LxBox,LyBox,LzBox)      , "Aluminium", vmat1); r=[ 4.5,0.0,0.0])
   sphere    = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(DxSphere)           , "Aluminium", vmat1); r=[ 3.0,0.0,0.0])
   ellipsoid = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidEllipsoid(DxEll,DyEll,DzEll), "Aluminium", vmat1); r=[ 1.5,0.0,0.0])
   cylinder  = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCylinder(DxCyl,LzCyl; Dy=DyCyl)     , smat1      , vmat1); r=[ 0.0,0.0,0.0])
   cone      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCone(DxCone,LzCone; Dy=DyCone)         , smat1      , vmat1); r=[-1.5,0.0,0.0])
   capsule   = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidCapsule(DxCap,LzCap; Dy=DyCap)     , smat1      , vmat1); r=[-3.0,0.0,0.0])
   pipe      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidPipe(DxPipe,LzPipe; Dy=DyPipe)         , smat1      , vmat1); r=[ 4.5,0.0,-2.5])
   beam      = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBeam(LxBeam,DyBeam,LzBeam)     , smat1      , vmat1); r=[ 3.0,0.0,-2.5])

   fileMesh = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidFileMesh(filename,1.0), smat1      , vmat1);  R=ModiaMath.rot2(180u"°")*ModiaMath.rot1(90u"°"), r=[ 1.5,0.0,-2.5])

   # coord Systems
   coordBox     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 4.5,0.0,0.0])
   coordSphere  = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 3.0,0.0,0.0])
   coordCyl     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 1.5,0.0,0.0])
   coordCone    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 0.0,0.0,0.0])
   coordCap     = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ -1.5,0.0,0.0])
   coord2       = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ -3.0,0.0,0.0])
   coordpipe    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 4.5,0.0,-2.5])
   coordBeam    = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 3.0,0.0,-2.5])
   coordMesh1   = Modia3D.Object3D(world,  Modia3D.CoordinateSystem(0.7) ; r=[ 1.5,0.0,-2.5])
   # Place text above the shapes
   boxText       = Modia3D.Object3D(box      , Modia3D.TextShape("SolidBox"     ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   sphereText    = Modia3D.Object3D(sphere   , Modia3D.TextShape("SolidSphere"  ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   ellipsoidText = Modia3D.Object3D(ellipsoid, Modia3D.TextShape("SolidEllipsoid";font=font); r=rtext, visualizeFrame=Modia3D.False)
   cylinderText  = Modia3D.Object3D(cylinder , Modia3D.TextShape("SolidCylinder"; font=font); r=rtext, visualizeFrame=Modia3D.False)
   coneText      = Modia3D.Object3D(cone     , Modia3D.TextShape("SolidCone"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   capsuleText   = Modia3D.Object3D(capsule  , Modia3D.TextShape("SolidCapsule" ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   pipeText      = Modia3D.Object3D(pipe     , Modia3D.TextShape("SolidPipe"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   beamText      = Modia3D.Object3D(beam     , Modia3D.TextShape("SolidBeam"    ; font=font); r=rtext, visualizeFrame=Modia3D.False)
   meshText1     = Modia3D.Object3D(coordMesh1, Modia3D.TextShape("SolidFileMesh"; font=font); r=rtext, visualizeFrame=Modia3D.False)
   meshText4     = Modia3D.Object3D(coordMesh1, Modia3D.TextShape("(.obj)"        ; font=font); r=ltext, visualizeFrame=Modia3D.False)
end
solid = Solids(sceneOptions=Modia3D.SceneOptions(visualizeFrames=false, defaultFrameLength=0.7) )
Modia3D.visualizeAssembly!( solid )

#=
println("Modia3D.volume(box) = ",       Modia3D.volume(solid.box.data.geo))
println("Modia3D.volume(sphere) = ",    Modia3D.volume(solid.sphere.data.geo))
println("Modia3D.volume(ellipsoid) = ", Modia3D.volume(solid.ellipsoid.data.geo))
println("Modia3D.volume(cylinder) = ",  Modia3D.volume(solid.cylinder.data.geo))
println("Modia3D.volume(cone) = ",      Modia3D.volume(solid.cone.data.geo))
println("Modia3D.volume(capsule) = ",   Modia3D.volume(solid.capsule.data.geo))
println("Modia3D.volume(pipe) = ",      Modia3D.volume(solid.pipe.data.geo))
println("Modia3D.volume(beam) = ",      Modia3D.volume(solid.beam.data.geo))
println("Modia3D.volume(fileMesh) = ", Modia3D.volume(solid.fileMesh.data.geo))
=#

#=
# println("Modia3D.bottomArea(sphere) = ",    Modia3D.bottomArea(solid.sphere.data.geo))
# println("Modia3D.bottomArea(ellipsoid) = ", Modia3D.bottomArea(solid.ellipsoid.data.geo))
# println("Modia3D.bottomArea(capsule) = ",   Modia3D.bottomArea(solid.capsule.data.geo))
# println("Modia3D.bottomArea(fileMesh) = ", Modia3D.bottomArea(solid.fileMesh.data.geo))
println("A sphere, an ellipsoid, a capsule, and a fileMesh do not have a bottom area.")
println("Modia3D.bottomArea(box) = ",       Modia3D.bottomArea(solid.box.data.geo))
println("Modia3D.bottomArea(cylinder) = ",  Modia3D.bottomArea(solid.cylinder.data.geo))
println("Modia3D.bottomArea(cone) = ",      Modia3D.bottomArea(solid.cone.data.geo))
println("Modia3D.bottomArea(pipe) = ",      Modia3D.bottomArea(solid.pipe.data.geo))
println("Modia3D.bottomArea(beam) = ",      Modia3D.bottomArea(solid.beam.data.geo))
=#

#=
#println("Modia3D.topArea(sphere) = ",    Modia3D.topArea(solid.sphere.data.geo))
#println("Modia3D.topArea(ellipsoid) = ", Modia3D.topArea(solid.ellipsoid.data.geo))
#println("Modia3D.topArea(capsule) = ",   Modia3D.topArea(solid.capsule.data.geo))
println("A sphere, an ellipsoid, and a capsule do not have a top area.")
println("Modia3D.topArea(box) = ",       Modia3D.topArea(solid.box.data.geo))
println("Modia3D.topArea(cylinder) = ",  Modia3D.topArea(solid.cylinder.data.geo))
println("Modia3D.topArea(cone) = ",      Modia3D.topArea(solid.cone.data.geo))
println("Modia3D.topArea(pipe) = ",      Modia3D.topArea(solid.pipe.data.geo))
println("Modia3D.topArea(beam) = ",      Modia3D.topArea(solid.beam.data.geo))

println("Modia3D.longestEdge(box) = ",       Modia3D.longestEdge(solid.box.data.geo))
println("Modia3D.longestEdge(sphere) = ",    Modia3D.longestEdge(solid.sphere.data.geo))
println("Modia3D.longestEdge(ellipsoid) = ", Modia3D.longestEdge(solid.ellipsoid.data.geo))
println("Modia3D.longestEdge(cylinder) = ",  Modia3D.longestEdge(solid.cylinder.data.geo))
println("Modia3D.longestEdge(cone) = ",      Modia3D.longestEdge(solid.cone.data.geo))
println("Modia3D.longestEdge(capsule) = ",   Modia3D.longestEdge(solid.capsule.data.geo))
println("Modia3D.longestEdge(pipe) = ",      Modia3D.longestEdge(solid.pipe.data.geo))
println("Modia3D.longestEdge(beam) = ",      Modia3D.longestEdge(solid.beam.data.geo))
println("Modia3D.longestEdge(fileMesh) = ", Modia3D.longestEdge(solid.fileMesh.data.geo))

println("Modia3D.lengthGeo(box) = ",       Modia3D.lengthGeo(solid.box.data.geo))
println("Modia3D.lengthGeo(sphere) = ",    Modia3D.lengthGeo(solid.sphere.data.geo))
println("Modia3D.lengthGeo(ellipsoid) = ", Modia3D.lengthGeo(solid.ellipsoid.data.geo))
println("Modia3D.lengthGeo(cylinder) = ",  Modia3D.lengthGeo(solid.cylinder.data.geo))
println("Modia3D.lengthGeo(cone) = ",      Modia3D.lengthGeo(solid.cone.data.geo))
println("Modia3D.lengthGeo(capsule) = ",   Modia3D.lengthGeo(solid.capsule.data.geo))
println("Modia3D.lengthGeo(pipe) = ",      Modia3D.lengthGeo(solid.pipe.data.geo))
println("Modia3D.lengthGeo(beam) = ",      Modia3D.lengthGeo(solid.beam.data.geo))
#println("Modia3D.lengthGeo(fileMesh) = ", Modia3D.lengthGeo(solid.fileMesh.data.geo))

println("Modia3D.centroid(box) = ",       Modia3D.centroid(solid.box.data.geo))
println("Modia3D.centroid(sphere) = ",    Modia3D.centroid(solid.sphere.data.geo))
println("Modia3D.centroid(ellipsoid) = ", Modia3D.centroid(solid.ellipsoid.data.geo))
println("Modia3D.centroid(cylinder) = ",  Modia3D.centroid(solid.cylinder.data.geo))
println("Modia3D.centroid(cone) = ",      Modia3D.centroid(solid.cone.data.geo))
println("Modia3D.centroid(capsule) = ",   Modia3D.centroid(solid.capsule.data.geo))
println("Modia3D.centroid(pipe) = ",      Modia3D.centroid(solid.pipe.data.geo))
println("Modia3D.centroid(beam) = ",      Modia3D.centroid(solid.beam.data.geo))
println("Modia3D.centroid(fileMesh) = ", Modia3D.centroid(solid.fileMesh.data.geo))
=#
#=
println("Modia3D.inertiaMatrix(box) = ",       Modia3D.inertiaMatrix(solid.box.data.geo,1))
println("Modia3D.inertiaMatrix(sphere) = ",    Modia3D.inertiaMatrix(solid.sphere.data.geo,1))
println("Modia3D.inertiaMatrix(ellipsoid) = ", Modia3D.inertiaMatrix(solid.ellipsoid.data.geo,1))
println("Modia3D.inertiaMatrix(cylinder) = ",  Modia3D.inertiaMatrix(solid.cylinder.data.geo,1))
#println("Modia3D.inertiaMatrix(cone) = ",      Modia3D.inertiaMatrix(solid.cone.data.geo,1))
#println("Modia3D.inertiaMatrix(capsule) = ",   Modia3D.inertiaMatrix(solid.capsule.data.geo,1))
println("Modia3D.inertiaMatrix(pipe) = ",      Modia3D.inertiaMatrix(solid.pipe.data.geo,1))
#println("Modia3D.inertiaMatrix(beam) = ",      Modia3D.inertiaMatrix(solid.beam.data.geo,1))
println("Modia3D.inertiaMatrix(fileMesh) = ", Modia3D.inertiaMatrix(solid.fileMesh.data.geo,1))
=#

#=
box
sphere
ellipsoid
cylinder
cone
capsule
pipe
beam
=#

function volumeCone(geo)
    h = LzCone
     # bottom area
    A = DxCone/2
    B = DyCone/2
    # top area
    a = A * geo.extras[1]
    b = B * geo.extras[1]
    # for cylinder (also top area)
    a_inner = A * geo.extras[1] * geo.extras[2]
    b_inner = B * geo.extras[1] * geo.extras[2]
    if geo.extras[1] == 0.0   # total cone
        return (h*pi/3)*A*B
    else    # frustum of a cone (with elliptic ground area): (h*pi/3) * (A*B^2 - a*b^2)/(B - b)
        # solid frustum of a cone - cylinder
        return (h*pi/3) * (A*B^2 - a*b^2)/(B - b) - pi*h*a_inner*b_inner
    end
end

@testset "Solids: volume test" begin
   @test isapprox(Modia3D.volume(solid.box.data.geo),       LxBox*LyBox*LzBox)
   @test isapprox(Modia3D.volume(solid.sphere.data.geo),    4/3*pi*(DxSphere/2)^3)
   @test isapprox(Modia3D.volume(solid.ellipsoid.data.geo), 4/3*pi*DxEll/2 * DyEll/2 * DzEll/2)
   @test isapprox(Modia3D.volume(solid.cylinder.data.geo),  pi * DxCyl/2 * DyCyl/2 * LzCyl)
   @test isapprox(Modia3D.volume(solid.cone.data.geo),      volumeCone(solid.cone.data.geo))
   @test isapprox(Modia3D.volume(solid.capsule.data.geo),   pi*LzCap*DxCap/2*DyCap/2 + 4/3*pi*(DyCap/2)^2*DxCap/2)
   @test isapprox(Modia3D.volume(solid.pipe.data.geo),      pi*LzPipe*(DxPipe*DyPipe - DxPipe*DyPipe*solid.pipe.data.geo.extras[2]^2)/4)
   @test isapprox(Modia3D.volume(solid.beam.data.geo),      LzBeam*(LxBeam*DyBeam + pi*(DyBeam/2)^2))
   @test isapprox(Modia3D.volume(solid.beam.data.geo),      LzBeam*(LxBeam*DyBeam + pi*(DyBeam/2)^2))
   @test isapprox(Modia3D.volume(solid.fileMesh.data.geo),  1.0)

   @test isapprox(Modia3D.bottomArea(solid.box.data.geo),   solid.box.data.geo.Lx*solid.box.data.geo.Ly)
   @test isapprox(Modia3D.bottomArea(solid.cylinder.data.geo), pi*DxCyl/2*DyCyl/2)
   @test isapprox(Modia3D.bottomArea(solid.cone.data.geo), pi*DxCone/2*DyCone/2)
end

println("... success of test_solidProperties.jl!")
end


text:
using PracticePackage
using Documenter

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

makedocs(;
    modules=[PracticePackage],
    authors="Dale-Black <djblack@uci.edu> and contributors",
    repo="https://github.com/Dale-Black/PracticePackage.jl/blob/{commit}{path}#{line}",
    sitename="PracticePackage.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Dale-Black.github.io/PracticePackage.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Dale-Black/PracticePackage.jl",
)


text:
cube, angles = BetaPictoris[:cube, :pa]

@testset "interface" begin
    X = flatten(cube)
    med = median(X, dims=1)

    d = ADI.fit(Classic(), cube)
    @test ADI.design(d) ≈ med

    t = reconstruct(d)
    @test size(t) == size(X)
    @test X .- t ≈ X .- med

end

@testset "evaluation" begin
    res = Classic()(cube, angles)
    res2 = Classic()(cube, angles, ref=cube)
    @test res == res2
    @test median(res) ≈ 0 atol = 1 / sqrt(length(res)) # within 1σ

    # using different statistic
    res = Classic(mean)(cube, angles)
    res2 = Classic(mean)(cube, angles, ref=cube)
    @test res == res2
    @test median(res) ≈ 0 atol = 1 / sqrt(length(res)) # within 1σ
end


text:
using CompScienceMeshes
using BEAST

sphere = CompScienceMeshes.tetmeshsphere(1.0, 0.35)
hemi = submesh(tet -> cartesian(center(chart(sphere,tet)))[3] < 0, sphere)

X = BEAST.nedelecd3d(sphere)
Y = BEAST.restrict(X, hemi)


text:
include("../src/Operation.jl")
using .Operation
using Test

@testset "Operation.jl" begin
    @test Operation.str2num("1 1 1 1", 2, 2) == [1 1; 1 1] #Квадратная
    @test Operation.str2num("1 1 1 1 1 1", 3, 2) == [1 1; 1 1; 1 1] #Для не квадратной
    @test Operation.transH("1 2 3 4", 2, 2) == [1 3; 2 4]#Транспонирование квадратной
    @test Operation.transH("1 2 3 4 5 6", 3, 2) == [1 3 5; 2 4 6]#Транспанирование не квадратной
    @test Operation.dethelp([1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0], 1, 1) == [5 6; 8 9]#Убрать строку и столбец
    @test Operation.det("1 1 1 1", 2, 2) == 0#Определитель единичной
    @test Operation.det("1 1 1 1 1 1", 3, 2) == "Матрица должна быть квадратной"#Ошибка ввода
    @test Operation.det("2 2 3 1 5 6 2 4 3", 3, 3) == -18#Произвольный определитель
    @test Operation.tr("1 1 1 1", 2, 2) == 2#След единичной
    @test Operation.tr("1 1 1 1 1 1", 3, 2) == "Матрица должна быть квадратной"#Ошибка ввода
    @test Operation.tr("1 2 3 4", 2, 2) == 5#След произвольной
    @test Operation.rang("2 2 3 1 5 6 2 4 3", 3, 3) == 3#Ранг произвольной
    @test Operation.rang("1 1 1 1", 2, 2) == 1#Ранг единичной
    @test Operation.obr("1 1 1 1", 2, 2) == "Определитель равен нулю"# Ошибка ввода
    @test Operation.obr("1 1 1 2", 2, 2) == [2 -1; -1 1]#Обратная произвольной
    @test Operation.obr("1 1 1 1 1 1", 3, 2) == "Матрица должна быть квадратной"#Ошибка ввода
    @test Operation.mul("1 1 1 1", 2, 2, "1 1 1 1", 2, 2) == [2 2; 2 2]#Умножение единичных
    @test Operation.mul("1 1 1 1 1 1", 3, 2, "1 1 1 1", 2, 2) == [2 2; 2 2; 2 2]#Умножение столбцы первой равны строкам второй
    @test Operation.mul("1 1 1 1 1 1", 2, 3, "1 1 1 1", 2, 2) == "Количество столбцов первой матрицы должно быть равно количеству строк второй"#Ошибка ввода
    @test Operation.mul("1 2 3 4", 2, 2, "2 11 5 5", 2, 2) == [12 21; 26 53]#Умножение произвольных
    @test Operation.plus("1 1 1 1", 2, 2, "1 1 1 1", 2, 2) == [2 2; 2 2]#Сложение единичных
    @test Operation.plus("1 1 1 1 1 1", 3, 3, "1 1 1 1", 2, 2) == "Размеры матриц должны совпадать"#Ошибка ввода
    @test Operation.plus("1 2 3 4", 2, 2, "1 2 3 4", 2, 2) == [2 4; 6 8]#Сложение произвольных
    @test Operation.minus("1 1 1 1", 2, 2, "1 1 1 1", 2, 2) == [0 0; 0 0]#Разность единичных
    @test Operation.minus("1 1 1 1 1 1", 3, 3, "1 1 1 1", 2, 2) == "Размеры матриц должны совпадать"#Ошибка ввода
    @test Operation.minus("2 2 13 5", 2, 2, "2 1 5 45", 2, 2) == [0 1; 8 -40]#Разность произвольных
    @test Operation.float_to_frac(5) == 5#Преобразование обычного числа
    @test Operation.float_to_frac(1.5) == "3/2"#Преобразование дроби > 1
    @test Operation.float_to_frac(0.000000001) == 0#Преобразование с выходом за погрешность
    @test Operation.float_to_frac(0.005) == "1/200"#Преобразование дроби < 1
    @test Operation.float_to_frac(0.6666666666) == "2/3"#преобразование периодической лроби
end

max_stars_repo_name: NikitaKutelev/kurs

max_issues_repo_name: NikitaKutelev/kurs

max_forks_repo_name: NikitaKutelev/kurs


text:
#Test code
#-------------------------------------------------------------------------------

using MiscTools

#No real test code yet... just demonstrate use:

#==Input data
===============================================================================#
filename = "./outfile.bin"

#==Tests
===============================================================================#
open(filename, "w") do io
	writeword(io, UInt32(5))
	writeword(io, Int32(-5))
	writeword(io, UInt8(5))
	writeword(io, Int8(-5))
	writeword(io, convert(Vector{UInt8}, "Output String"))
	writeword(io, UInt8(5))
	writeword(io, UInt8(5))
end

:Test_Complete


text:
###########################################################################
#
#  NUTS variant Hamiltonian Monte Carlo
#
#  Reference: Hoffman M.D, Gelman A.The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.
#  arXiv, 2011
#
#  Parameters :
#    - maxdoublings : max number of doublings in inner steps
#    - tuner: for tuning the NUTS parameters
#
###########################################################################

export NUTS

println("Loading NUTS(maxdoublings, tuner) sampler")

###########################################################################
#                  NUTS specific 'tuners'
###########################################################################
abstract NUTSTuner <: MCMCTuner

# TODO : describe here the leapSteps adaptation

###########################################################################
#                  NUTS type
###########################################################################
immutable NUTS <: MCMCSampler
  maxdoublings::Int
  tuner::Union(Nothing, NUTSTuner)

  function NUTS(i::Int, t::Union(Nothing, NUTSTuner))
    @assert i>0 "max doublings should be > 0"
    @assert i<20 "max doublings reasonably be < 20"
    new(i,t)
  end
end
NUTS(i::Int=5) = NUTS(i , nothing)
NUTS(t::NUTSTuner) = NUTS(5, t)
# keyword args version
NUTS(;maxdoublings::Int=5, tuner::Union(Nothing, NUTSTuner)=nothing) = NUTS(maxdoublings, tuner)

###########################################################################
#                  HMC task
###########################################################################

####  Helper functions and types for HMC sampling task
uturn(s1::HMCSample, s2::HMCSample) = dot(s2.pars-s1.pars, s1.m) < 0. || dot(s2.pars-s1.pars, s2.m) < 0.

####  HMC task
function SamplerTask(model::MCMCModel, sampler::NUTS, runner::MCMCRunner)
    local epsilon, u_slice
    local state0  # starting state of each loop
    local scale

 	@assert hasgradient(model) "NUTS sampler requires model with gradient function"

	# hook inside Task to allow remote resetting
	task_local_storage(:reset,
	             (resetPars::Vector{Float64}) -> (state0 = HMCSample(copy(resetPars));
	                                              calc!(state0, model.evalallg)) )

	# initialization
	state0 = HMCSample(copy(model.init))
	calc!(state0, model.evalallg)
  @assert isfinite(state0.logTarget) "Initial values out of model support, try other values"
	scale = model.scale

	# find initial value for epsilon
	epsilon = 1.
	state0.m = randn(model.size) .* scale
	state1 = leapfrog(state0, epsilon, model.evalallg)

	ratio = exp(state0.H-state1.H)
	a = 2*(ratio>0.5)-1.
	while ratio^a > 2^-a
		epsilon *= 2^a
		state1 = leapfrog(state0, epsilon, model.evalallg)
		ratio = exp(state0.H-state1.H)
	end

	# buidtree function
	function buildTree(state, dir, j, ll)  # TODO : pass espilon, u_slice, state0 ?
		local state1, n1, s1, alpha1, nalpha1
		local state2, n2, s2, alpha2, nalpha2
		local state_plus, state_minus
		local dummy
		const deltamax = 100

		if j == 0
			state1 = leapfrog(state, dir*epsilon, ll)
			n1 = ( u_slice <= -state1.H ) + 0
			s1 = u_slice < ( deltamax-state1.H )

			return state1, state1, state1, n1, s1, min(1., exp(state0.H-state1.H)), 1
		else
			state_minus, state_plus, state1, n1, s1, alpha1, nalpha1 = buildTree(state, dir, j-1, ll)
			if s1
				if dir == -1
					state_minus, dummy, state2, n2, s2, alpha2, nalpha2 = buildTree(state_minus, dir, j-1, ll)
	 			else
	 				dummy, state_plus, state2, n2, s2, alpha2, nalpha2 = buildTree(state_plus, dir, j-1, ll)
	 			end
	 			if rand() <= n2/(n2+n1)
	 				state1 = deepcopy(state2)
	 			end

	 			alpha1 += alpha2
	 			nalpha1 += nalpha2
	 			s1 = s2 && !uturn(state_minus, state_plus)
	 			n1 += n2
	 		end

	 		return state_minus, state_plus, state1, n1, s1, alpha1, nalpha1
		end
	end

	### adaptation parameters  # TODO : should be part of the tuner
	const delta = 0.7  # target acceptance
	const nadapt = 1000  # nb of steps to adapt epsilon
	const gam = 0.05
	const kappa = 0.75
	const t0 = 10
	### adaptation inital values
	hbar = 0.
	mu = log(10*epsilon)
	lebar = 0.0

	### main loop
	# using a 'for' loop instead of 'while true' for this sampler
	#  because leapStep adaptation requires a sampling counter
  i = 1
  while true
 		local alpha, nalpha, n, s, j, n1, s1
 		local dummy, state_minus, state_plus, state, state1

 		state0.m = randn(model.size) .* scale
 		update!(state0)

 		u_slice  = log(rand())-state0.H # use log ( != paper) to avoid underflow

 		state = deepcopy(state0)
 		state_minus = deepcopy(state0)
 		state_plus = deepcopy(state0)

 		# inner loop
 		j, n = 0, 1
 		s = true
 		while s && j < sampler.maxdoublings
 			dir = randbool() ? 1 : -1
 			if dir == -1
 				state_minus, dummy, state1, n1, s1, alpha, nalpha = buildTree(state_minus, dir, j, model.evalallg)
 			else
 				dummy, state_plus, state1, n1, s1, alpha, nalpha = buildTree(state_plus, dir, j, model.evalallg)
 			end
 			if s1 && rand() < n1/n  # accept
 				state = deepcopy(state1)
 			end
 			n += n1
 			j += 1
 			s = s1 && !uturn(state_minus, state_plus)
 		end

 		# epsilon adjustment
 		if i <= nadapt  # warming up period
 			hbar = hbar * (1-1/(i+t0)) + (delta-alpha/nalpha)/(i+t0)
			le = mu-sqrt(i)/gam*hbar
			lebar = i^(-kappa) * le + (1-i^(-kappa)) * lebar
			epsilon = exp(le)
		else # post warm up, keep dual epsilon
			epsilon = exp(lebar)
		end

		ms = MCMCSample(state.pars, state.logTarget, state.grad,
		                state0.pars, state0.logTarget, state0.grad,
		                {"epsilon" => epsilon, "ndoublings" => j} )
		produce(ms)

		state0 = deepcopy(state)

    i += 1
	end

end


text:
#######################
# recursive fibonacci #
#######################

perf_micro_fib(n) = n < 2 ? n : perf_micro_fib(n-1) + perf_micro_fib(n-2)

############
# parseint #
############

function perf_micro_parseint(t)
    local n, m
    for i=1:t
        n = rand(UInt32)
        s = hex(n)
        m = UInt32(parse(Int64,s,16))
    end
    @assert m == n
    return n
end

##############
# mandelbrot #
##############

function mandel(z)
    c = z
    maxiter = 80
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

perf_micro_mandel() = [mandel(complex(r,i)) for i=-1.:.1:1., r=-2.0:.1:0.5]

#############
# quicksort #
#############

function quicksort!(a, lo, hi)
    i, j = lo, hi
    while i < hi
        pivot = a[(lo+hi)>>>1]
        while i <= j
            while a[i] < pivot; i += 1; end
            while a[j] > pivot; j -= 1; end
            if i <= j
                a[i], a[j] = a[j], a[i]
                i, j = i+1, j-1
            end
        end
        if lo < j; quicksort!(a,lo,j); end
        lo, j = i, hi
    end
    return a
end

perf_micro_quicksort(n) = quicksort!(rand(n), 1, n)

########
# πsum #
########

function perf_micro_pisum()
    sum = 0.0
    for j = 1:500
        sum = 0.0
        for k = 1:10000
            sum += 1.0/(k*k)
        end
    end
    return sum
end

###############
# randmatstat #
###############

function perf_micro_randmatstat(t)
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i=1:t
        a = randn(n,n)
        b = randn(n,n)
        c = randn(n,n)
        d = randn(n,n)
        P = [a b c d]
        Q = [a b; c d]
        v[i] = trace((P.'*P)^4)
        w[i] = trace((Q.'*Q)^4)
    end
    return (std(v)/mean(v), std(w)/mean(w))
end

##############
# randmatmul #
##############

perf_micro_randmatmul(t) = rand(t,t)*rand(t,t)


text:
# Run `julia make.jl` in this folder to generate .html pages in a build/
# directory, then open index.md for docs

using Documenter, DPMeansClustering
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/MurrellGroup/DPMeansClustering.jl.git",
    julia  = "nightly",
    osname = "osx"
makedocs()


text:
module rollingBall1

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmatGraphics = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of Graphics
vmatSolids = Modia3D.Material(color="Red" , transparency=0.0)         # material of solids
vmatTable = Modia3D.Material(color="Green", transparency=0.1)         # material of table

cmatTable = "BilliardTable"
cmatBall = "BilliardBall"

LxGround = 3.0
LyBox = 0.5
LzBox = 0.02
diameter = 0.06
@assembly Table(world) begin
  withBox = Modia3D.Solid(Modia3D.SolidBox(LxGround, LyBox, LzBox) , "BilliardTable", vmatTable; contactMaterial = cmatTable)
  box1 = Modia3D.Object3D(world, withBox, r=[1.5, 0.0, -LzBox/2], fixed=true)
end

@assembly OneRollingBall() begin
  world = Modia3D.Object3D(visualizeFrame=false)
  table = Table(world)
  ball1 = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(diameter), "BilliardBall", vmatSolids ;
                           contactMaterial = cmatBall), fixed = false, r=[0.2, 0.0, diameter/2], v_start=[3.0, 0.0, 0.0], visualizeFrame=true )
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
bill = OneRollingBall(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=false,
                       defaultFrameLength=0.1,nz_max = 100, enableContactDetection=true, visualizeContactPoints=false, visualizeSupportPoints=false))
model = Modia3D.SimulationModel( bill )

using PyPlot
using PyCall

pyplot_rc = PyCall.PyDict(PyPlot.matplotlib["rcParams"])
pyplot_rc["font.family"]      = "sans-serif"
pyplot_rc["font.sans-serif"]  = ["Calibri", "Arial", "Verdana", "Lucida Grande"]
pyplot_rc["font.size"]        = 12.0
pyplot_rc["lines.linewidth"]  = 1.5
pyplot_rc["grid.linewidth"]   = 0.5
pyplot_rc["axes.grid"]        = true
pyplot_rc["axes.titlesize"]   = "medium"
pyplot_rc["figure.titlesize"] = "medium"

result = ModiaMath.simulate!(model; stopTime=0.5, tolerance=1e-8, log=false)

ModiaMath.plot(result, [("ball1.r[1]"),
                        ("ball1.r[3]"),
                        ("ball1.v[1]"),
                        ("ball1.w[2]")])

ball_r = result.series["ball1.r"][:,1]
t      = result.series["time"]
l_t    = length(t)

dist = ( ball_r[2] - ball_r[1] ) / (t[2] - t[1] )
distEnd = ( ball_r[l_t] - ball_r[l_t - 1] ) / (t[l_t] - t[l_t-1] )

println("... success of rollingBall1.jl!")
end


text:
#test-infrastructure.jl
#tests unums infrastructure (everything except mathematical operations).

################################################################################
## UNIT TESTING

#integer operations and helpers for the type constructor
include("./unit/unit-int64-op.jl")
include("./unit/unit-int64-cmp.jl")

#functions which directly participate in construction and validation of unums
#include("./unit/unit-helpers.jl")

#type constructor testing
include("./unit/unit-unum.jl")
include("./unit/unit-ubound.jl")

#testing things that get done to unums
include("./unit/unit-constants.jl")
include("./unit/unit-convert.jl")
include("./unit/unit-comparison.jl")
include("./unit/unit-operations.jl")

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
# Install the HurdleDMR package
#Pkg.clone("https://github.com/AsafManela/HurdleDMR.jl")

# Add parallel workers and make package available to workers
using Distributed
import HurdleDMR
addprocs(Sys.CPU_THREADS-2)
@everywhere using HurdleDMR

# Setup your data into an n-by-p covars matrix, and a (sparse) n-by-d counts matrix
# Here we generate some random data
using CSV, GLM, DataFrames, Distributions, Random, SparseArrays, CategoricalArrays, InteractiveUtils #, StatsModels
n = 200
p = 3
d = 4

Random.seed!(13)
m = 1 .+ rand(Poisson(5),n)
vs = rand(n,p)
covarsdf = DataFrame(vs,[:y, :x, :z])
covarsdf[:cat] = CategoricalArray(rand(["1","2","3"], n))
ηfn(vi,g) = exp.([0 + i*sum(vi) - parse(Float64, g) for i=1:d])
q = [ηfn(vs[i,:], covarsdf[i,:cat]) for i=1:n]
for i=1:n
  q[i] ./= sum(q[i])
end
counts = convert(SparseMatrixCSC{Float64,Int}, hcat(broadcast((qi,mi)->rand(Multinomial(mi, qi)),q,m)...)')

# construct equivalent covars matrix so we can show how that api works too
covars = ModelMatrix(ModelFrame(@formula(y ~ cat + x + y + z), covarsdf)).m[:,2:end]

## To fit a hurdle distribtued multiple regression (hdmr):
m = hdmr(covars, counts)

prem = vec(sum(counts, dims=2))
m1 = hdmr(covars, counts[:,1:2]; m=prem)
m2 = hdmr(covars, counts[:,3:4]; m=prem)

hcat(m1.coefspos, m2.coefspos) == m.coefspos
hcat(m1.coefszero, m2.coefszero) == m.coefszero

# or with a dataframe and formula
mf = @model(h ~ cat + x + y + z, c ~ cat + x + y)
m2 = fit(HDMR, mf, covarsdf, counts)
coef(m2.model)

using GLM
lm1 = fit(LinearModel, @formula(y ~ 0 + cat + x + z), covarsdf)
yhat = predict(lm1, covarsdf[1:10,:])
mm = lm1
mdf = covarsdf[1:10,:]
mm.mf.terms
mm.mf.terms.intercept = false
newTerms = StatsModels.dropresponse!(mm.mf.terms)
# create new model frame/matrix
mf = ModelFrame(newTerms, mdf; contrasts = mm.mf.contrasts)
newX = ModelMatrix(mf).m

# m2.mm.m
# trmszero = HurdleDMR.getrhsterms(mf, :h)
# trmspos = HurdleDMR.getrhsterms(mf, :c)
# trms, inzero, inpos = HurdleDMR.mergerhsterms(trmszero,trmspos)
# mf, mm, counts = HurdleDMR.createmodelmatrix(trms, covarsdf, counts, Dict())
# inzero
#
# inzero, inpos = HurdleDMR.mapins(inzero, inpos, mm)
# projdir = HurdleDMR.ixprojdir(trms, :y, mm)
# where the h ~ equation is the model for zeros (hurdle crossing) and c ~ is the model for positive counts

# in either case we can get the coefficients matrix for each variable + intercept as usual with
coefspos, coefszero = coef(m)

# By default we only return the AICc maximizing coefficients.
# To also get back the entire regulatrization paths, run
paths = fit(HDMRPaths, mf, covarsdf, counts)

coef(paths; select=AllSeg())

# To get a sufficient reduction projection in direction of vy
z = srproj(m,counts,1,1)

# Counts inverse regression (cir) allows us to predict a covariate with the counts and other covariates
# Here we use hdmr for the backward regression and another model for the forward regression
# This can be accomplished with a single command, by fitting a CIR{HDMR,FM} where the forward model is FM <: RegressionModel
cir = fit(CIR{HDMR,LinearModel},mf,covarsdf,counts,:y; nocounts=true)
# where the argument nocounts=true means we also fit a benchmark model without counts

# we can get the forward and backward model coefficients with
coefbwd(cir)
coeffwd(cir)

# the fitted model can be used to predict vy with new data
yhat = predict(cir, covarsdf[1:10,:], counts[1:10,:])

# and we can also predict only with the other covariates, which in this case
# is just a linear regression
yhat_nocounts = predict(cir, covarsdf[1:10,:], counts[1:10,:]; nocounts=true)

# To fit a distribtued multinomial regression (dmr):
m = dmr(covars, counts)

# or with a dataframe and formula
mf = @model(c ~ vy + v1 + v2)
m = fit(DMR, mf, covarsdf, counts)

# in either case we can get the coefficients matrix for each variable + intercept as usual with
coef(m)

# By default we only return the AICc maximizing coefficients.
# To also get back the entire regulatrization paths, run
paths = fit(DMRPaths, mf, covarsdf, counts)

# we can now select, for example the coefficients that minimize CV mse (takes a while)
coef(paths; select=MinCVmse(Kfold(size(covars,1), 5)))

# To get a sufficient reduction projection in direction of vy
z = srproj(m,counts,1)

# A multinomial inverse regression (mnir) uses dmr for the backward regression and another model for the forward regression
# This can be accomplished with a single command, by fitting a CIR{DMR,FM} where FM is the forward RegressionModel
mnir = fit(CIR{DMR,LinearModel},mf,covarsdf,counts,:y)

# we can get the forward and backward model coefficients with
coefbwd(mnir)
coeffwd(mnir)

# the fitted model can be used to predict vy with new data
yhat = predict(mnir, covarsdf[1:10,:], counts[1:10,:])

# Suppose instead we want to predict a discrete variable, then perhaps use a Poisson GLM as follows
mnir = fit(CIR{DMR,GeneralizedLinearModel},mf,covarsdf,counts,:y,Gamma())

using Coverage
# defaults to src/; alternatively, supply the folder name as argument
coverage = process_folder(srcfolder)
# Get total coverage for all Julia files
covered_lines, total_lines = get_summary(coverage)
# Or process a single file
@show get_summary(process_file("src/HurdleDMR.jl"))

function focus(X; j=:)
  X[:,j]
end
function focus2(X; j=axes(X,2))
  X[:,j]
end

A = sprand(5000,4000,0.2)
focus(A)
focus2(A)
focus(A; j=3)
focus2(A; j=3)

j = axes(A,2)
j == 3:3
A
in(4000,axes(A,2))
in(4000,3)
A[axes(A,1),axes(A,2)]
j = 4
focusj = to_axes(A,(:,4))
A[axes(A,1),axes(A,2)]
A[:,4]
A[focusj[1],focusj[2]]

typeof(j)
typeof(:)
colon(1,1,3)
j == 3:3

# clean_folder("/home/amanela/.julia/v0.4")