text:
"""
Provide interpolation and integration.
"""
module Interp

# using StaticArrays, FastGaussQuadrature, CompositeGrids
using StaticArrays, FastGaussQuadrature
using ..SimpleG
using ..CompositeG

#include("chebyshev.jl")

# include("simple.jl")
# using .SimpleG
# include("composite.jl")
# using .CompositeG

abstract type InterpStyle end
struct LinearInterp <: InterpStyle end
struct ChebInterp <: InterpStyle end
struct CompositeInterp <: InterpStyle end
const LINEARINTERP = LinearInterp()
const CHEBINTERP = ChebInterp()
const COMPOSITEINTERP = CompositeInterp()

InterpStyle(::Type) = LinearInterp()
InterpStyle(::Type{<:SimpleG.BaryCheb}) = ChebInterp()
InterpStyle(::Type{<:CompositeG.Composite}) = CompositeInterp()

# return types of findneighbor, contains complete information for interp
abstract type InterpNeighbor end

# for linear interp, two nearest neighbor are needed
struct LinearNeighbor{T<:AbstractFloat} <: InterpNeighbor
    index::UnitRange{Int}
    grid::SVector{2,T}
    x::T
end

# for cheb interp, whole cheb grid is needed, together with interp weight
# index is recorded as first and last index involved, useful for composite grids
struct ChebNeighbor{T<:AbstractFloat} <: InterpNeighbor
    index::UnitRange{Int}
    grid::Vector{T}
    weight::Vector{T}
    x::T
end

"""
    function findneighbor(xgrid::T, x; method=:default) where {T}

Find neighbor grid points and related information for extrapolating
the value of x on xgrid.

#Members:
- xgrid: grid to be interpolated
- x: value to be interpolated
- method: :default use optimized method, :linear use linear interp.
"""
function findneighbor(xgrid::T, x; method=InterpStyle(T)) where {T}
    return findneighbor(method, xgrid, x)
end

function findneighbor(::LinearInterp, xgrid, x)
    T = eltype(xgrid.grid)
    xi0, xi1 = 0, 0
    if (x <= xgrid[firstindex(xgrid)])
        xi0 = 1
        xi1 = 2
    elseif (x >= xgrid[lastindex(xgrid)])
        xi0 = lastindex(xgrid) - 1
        xi1 = xi0 + 1
    else
        xi0 = floor(xgrid, x)
        xi1 = xi0 + 1
    end

    x0, x1 = xgrid[xi0], xgrid[xi1]

    return LinearNeighbor{T}(xi0:xi1, [x0, x1], x)
end

function findneighbor(::ChebInterp, xgrid, x)
    T = eltype(xgrid.grid)
    return ChebNeighbor{T}(1:xgrid.size, xgrid.grid, xgrid.weight, x)
end

function findneighbor(::CompositeInterp, xgrid, x)
    if CompositeG.getbottomtype(xgrid) <: SimpleG.BaryCheb
        T = eltype(xgrid.grid)
        curr = xgrid
        xi0 = 1
        while !(typeof(curr) <: SimpleG.BaryCheb)
            i = floor(curr.panel, x)
            xi0 += curr.inits[i] - 1
            curr = curr.subgrids[i]
        end
        return ChebNeighbor{T}(xi0:(curr.size-1+xi0), curr.grid, curr.weight, x)
    else
        return findneighbor(LinearInterp(), xgrid, x)
    end
end

"""
    function dataslice(data, axes, indices)

Return a view of data sliced on given axes with given indices.
Works like view(data, (:, ..., :, i_1:f_1,  :, ..., i_n:f_n, :, ..., :)).
Type unstable unless slice dims are constant.

#Members:
- data: data to be sliced.
- axes: axes to be sliced. accept Int or NTuple{DIM, Int} for single or multiple axes. when omitted, assume all axes.
- indices: indices of slicing. accept UnitRange{Int} or Vector of UnitRange{Int} like 2:8 or [2:8, 3:7]
"""
function dataslice(data, axes::Int, indices)
    return selectdim(data, axes, indices[1])
end

function dataslice(data, axes::Int, indices::UnitRange{Int})
    return selectdim(data, axes, indices)
end

function dataslice(data, axes::NTuple{DIM,Int}, indices) where {DIM}
    # @assert DIM == length(indices)
    slice = data
    for i in 1:DIM
        slice = selectdim(slice, axes[i], indices[i])
    end
    return slice
end

@inline function dataslice(data, indices::UnitRange{Int})
    return view(data, indices)
end

@inline function dataslice(data, indices::Vector{UnitRange{Int}})
    return view(data, indices...)
end

"""
    function interpsliced(neighbor, data; axis=1)
Interpolate with given neighbor and sliced data. Assume data already sliced on given axis.

#Members:
- neighbor: neighbor from findneighbor()
- data: sliced data
- axis: axis sliced and to be interpolated
"""
function interpsliced(neighbor, data; axis=1)
    return dropdims(mapslices(u -> _interpsliced(neighbor, u), data, dims=axis), dims=axis)
end

function interpsliced(neighbor, data::AbstractMatrix; axis=1)
    # trying to make it type stable for matrix
    if axis == 1
        return map(u -> _interpsliced(neighbor, u), eachcol(data))
    elseif axis == 2
        return map(u -> _interpsliced(neighbor, u), eachrow(data))
    else
        throw(DomainError(axis, "axis should be 1 or 2 for Matrix"))
    end
end

function interpsliced(neighbor, data::AbstractVector; axis=1)
    return _interpsliced(neighbor, data)
end

# function interpsliced(neighbor, data; axis=1)
#     if ndims(data) == 1
#         return _interpsliced(neighbor, data)
#     else
#         return dropdims(mapslices(u->_interpsliced(neighbor, u), data, dims=axis), dims=axis)
#     end
# end

function _interpsliced(neighbor::LinearNeighbor, data)
    # data should be sliced priorly
    dx0, dx1 = neighbor.x - neighbor.grid[1], neighbor.grid[2] - neighbor.x

    d0, d1 = data[1], data[2]

    g = d0 * dx1 + d1 * dx0

    gx = g / (dx0 + dx1)
    return gx
end

function _interpsliced(neighbor::ChebNeighbor, data)
    # data should be sliced priorly
    return SimpleG.barycheb(length(neighbor.grid), neighbor.x, data, neighbor.weight, neighbor.grid)
end

function interpND(data, xgrids, xs; method=LINEARINTERP)
    #WARNING: This function works but is not type stable
    dim = length(xs)
    neighbors = [findneighbor(xgrids[i], xs[i]; method) for i in 1:dim]
    indices = [nei.index for nei in neighbors]

    # data_slice = dataslice(data, indices)
    data_slice = view(data, indices...)
    curr_data_slice = copy(data_slice)

    for i in 1:dim-1
        curr_data_slice = interpsliced(neighbors[i], curr_data_slice)
    end

    return interpsliced(neighbors[dim], curr_data_slice)
end

"""
    function linearND(data, xgrids, xs)

linear interpolation of data(xs)

#Arguments:
- xgrids: n-dimensional grids, xgrids[i] is a 1D grid
- data: n-dimensional array of data
- xs: list of x, x[i] corresponds to xgrids[i]
"""
function linearND(data, xgrids, xs)
    @inline function enumX(xs)
        # dim = length(xs)
        # result = ones(Float64, 2^dim)
        # for i in 1:2^dim
        #     for j in 1:dim
        #         result[i] *= xs[j]^digits(i-1,base=2,pad=dim)[j]
        #     end
        # end
        # return result
        return [prod(xs .^ digits(i - 1, base=2, pad=length(xs))) for i in 1:2^length(xs)]
    end
    @inline function f(as, xs)
        # dim = length(xs)
        # result = 0.0
        # ex = enumX(xs)
        # for i in 1:2^dim
        #     result += as[i]*ex[i]
        # end
        # return result
        return sum(as .* enumX(xs))
    end

    dim = length(xs)

    # find grid points below and above xs
    xis = zeros(Int, (dim, 2))
    for i in 1:dim
        xi0, xi1 = 0, 0
        x = xs[i]
        if (x <= xgrids[i].grid[firstindex(xgrids[i])])
            xi0 = 1
            xi1 = 2
        elseif (x >= xgrids[i].grid[lastindex(xgrids[i])])
            xi0 = lastindex(xgrid) - 1
            xi1 = xi0 + 1
        else
            xi0 = floor(xgrids[i], xs[i])
            xi1 = xi0 + 1
        end
        xis[i, 1], xis[i, 2] = xi0, xi1
    end

    # data value at nearby grid points
    datas = [data[[xis[j, 1+digits(i - 1, base=2, pad=dim)[j]] for j in 1:dim]...] for i in 1:2^dim]
    # datas = zeros(Float64, 2^dim)
    # for i in 1:2^dim
    #     datas[i] = data[[xis[j , 1+digits(i,base=2,pad=dim)[j]]  for j in 1:dim]...]
    # end

    # create x matrix
    xmat = ones(Float64, (2^dim, 2^dim))
    for i in 1:2^dim
        ii = digits(i - 1, base=2, pad=dim)
        xxs = [xgrids[j][xis[j, 1+ii[j]]] for j in 1:dim]
        xmat[:, i] = enumX(xxs)
    end

    xtran = copy(transpose(xmat))
    as = xtran \ datas
    #println(xtran)
    #println(as)
    #println(xtran*as)
    #println(datas)
    #@assert xtran*as == datas "$(xtran*as), $(datas), $(xtran)"

    return f(as, xs)
end

@inline function linearND(data::Matrix, xgrids, xs)
    # for 2d, use linear2D
    @assert length(xgrids) == length(xs) == 2
    return linear2D(data, xgrids[1], xgrids[2], xs[1], xs[2])
end

@inline function linearND(data::Vector, xgrids, xs)
    # for 2d, use linear2D
    @assert length(xgrids) == length(xs) == 1
    return linear2D(data, xgrids[1], xs[1])
end

"""
   linear2D(data, xgrid, ygrid, x, y)

linear interpolation of data(x, y)

#Arguments:
- xgrid: one-dimensional grid of x
- ygrid: one-dimensional grid of y
- data: two-dimensional array of data
- x: x
- y: y
"""
@inline function linear2D(data, xgrid, ygrid, x, y)

    xarray, yarray = xgrid.grid, ygrid.grid

    # if (
    #     x <= xarray[firstindex(xgrid)] ||
    #     x >= xarray[lastindex(xgrid)] ||
    #     y <= yarray[firstindex(ygrid)] ||
    #     y >= yarray[lastindex(ygrid)]
    # )
    #     return 0.0
    # end
    # xi0, yi0 = floor(xgrid, x), floor(ygrid, y)
    # xi1, yi1 = xi0 + 1, yi0 + 1

    xi0, xi1, yi0, yi1 = 0, 0, 0, 0
    if (x <= xarray[firstindex(xgrid)])
        xi0 = 1
        xi1 = 2
    elseif (x >= xarray[lastindex(xgrid)])
        xi0 = lastindex(xgrid) - 1
        xi1 = xi0 + 1
    else
        xi0 = floor(xgrid, x)
        xi1 = xi0 + 1
    end

    if (y <= yarray[firstindex(ygrid)])
        yi0 = 1
        yi1 = 2
    elseif (y >= yarray[lastindex(ygrid)])
        yi0 = lastindex(ygrid) - 1
        yi1 = yi0 + 1
    else
        yi0 = floor(ygrid, y)
        yi1 = yi0 + 1
    end

    dx0, dx1 = x - xarray[xi0], xarray[xi1] - x
    dy0, dy1 = y - yarray[yi0], yarray[yi1] - y

    d00, d01 = data[xi0, yi0], data[xi0, yi1]
    d10, d11 = data[xi1, yi0], data[xi1, yi1]

    g0 = data[xi0, yi0] * dx1 + data[xi1, yi0] * dx0
    g1 = data[xi0, yi1] * dx1 + data[xi1, yi1] * dx0

    gx = (g0 * dy1 + g1 * dy0) / (dx0 + dx1) / (dy0 + dy1)
    return gx
end

"""
    function linear1D(data, xgrid, x)

linear interpolation of data(x)

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
"""
@inline function linear1D(data, xgrid, x)

    xarray = xgrid.grid

    xi0, xi1 = 0, 0
    if (x <= xarray[firstindex(xgrid)])
        xi0 = 1
        xi1 = 2
    elseif (x >= xarray[lastindex(xgrid)])
        xi0 = lastindex(xgrid) - 1
        xi1 = xi0 + 1
    else
        xi0 = floor(xgrid, x)
        xi1 = xi0 + 1
    end

    dx0, dx1 = x - xarray[xi0], xarray[xi1] - x

    d0, d1 = data[xi0], data[xi1]

    g = d0 * dx1 + d1 * dx0

    gx = g / (dx0 + dx1)
    return gx
end

"""
    function interp1D(data, xgrid, x; axis=1, method=InterpStyle(T))

linear interpolation of data(x) with single or multiple dimension.
For 1D data, return a number; for multiple dimension, reduce the given axis.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
- axis: axis to be interpolated in data
- method: by default use optimized method; use linear interp if Interp.LinearInterp()
"""
function interp1D(data, xgrid::T, x; axis=1, method=InterpStyle(T)) where {T}
    return dropdims(mapslices(u -> interp1D(method, u, xgrid, x), data, dims=axis), dims=axis)
end

function interp1D(data::AbstractMatrix, xgrid::T, x; axis=1, method=InterpStyle(T)) where {T}
    # trying to make it type stable for matrix
    if axis == 1
        return map(u -> interp1D(method, u, xgrid, x), eachcol(data))
    elseif axis == 2
        return map(u -> interp1D(method, u, xgrid, x), eachrow(data))
    else
        throw(DomainError(axis, "axis should be 1 or 2 for Matrix"))
    end
end

function interp1D(data::AbstractVector, xgrid::T, x; axis=1, method=InterpStyle(T)) where {DT,T}
    return interp1D(method, data, xgrid, x)
end

"""
    function interp1D(::LinearInterp,data, xgrid, x)

linear interpolation of data(x), use floor and linear1D

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
"""
function interp1D(::LinearInterp, data, xgrid, x)
    return linear1D(data, xgrid, x)
end

"""
    function interp1D(::ChebInterp, data, xgrid, x)

linear interpolation of data(x), barycheb for BaryCheb grid

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
"""
function interp1D(::ChebInterp, data, xgrid, x)
    return SimpleG.barycheb(xgrid.size, x, data, xgrid.weight, xgrid.grid)
end

"""
    function interp1D(::CompositeInterp,data, xgrid, x)

linear interpolation of data(x),
first floor on panel to find subgrid, then call interp1D on subgrid

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: x
"""
function interp1D(::CompositeInterp, data, xgrid, x)
    i = floor(xgrid.panel, x)
    head, tail = xgrid.inits[i], xgrid.inits[i] + xgrid.subgrids[i].size - 1
    return interp1D(view(data, head:tail), xgrid.subgrids[i], x)
end

"""
    function interp1DGrid(data, xgrid, grid; axis=1, method=InterpStyle(T))
For 1D data, do interpolation of data(grid[1:end]), return a Vector.
For ND data, do interpolation of data(grid[1:end]) at given axis, return data of same dimension.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- grid: points to be interpolated on
- axis: axis to be interpolated in data
- method: by default use optimized method; use linear interp if :linear
"""
function interp1DGrid(data, xgrid::T, grid; axis=1, method=InterpStyle(T)) where {T}
    return mapslices(u -> interp1DGrid(method, u, xgrid, grid), data, dims=axis)
end

function interp1DGrid(data::AbstractVector, xgrid::T, grid; axis=1, method=InterpStyle(T)) where {T}
    return interp1DGrid(method, data, xgrid, grid)
end

"""
    function interp1DGrid(::Union{LinearInterp,ChebInterp}, data, xgrid, grid)

linear interpolation of data(grid[1:end]), return a Vector
simply call interp1D on each points

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- grid: points to be interpolated on
"""
function interp1DGrid(::Union{LinearInterp,ChebInterp}, data, xgrid, grid)
    ff = zeros(eltype(data), length(grid))
    for (xi, x) in enumerate(grid)
        ff[xi] = interp1D(data, xgrid, x)
        # if x == 1.1386851268496132
        #     println(xgrid.bound)
        #     println(xgrid.grid)
        # end
    end
    return ff
end

"""
    function interp1DGrid(::CompositeInterp, data, xgrid, grid)

linear interpolation of data(grid[1:end]), return a Vector
grid should be sorted.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- grid: points to be interpolated on
"""
function interp1DGrid(::CompositeInterp, data, xgrid, grid)
    ff = zeros(eltype(data), length(grid))

    init, curr = 1, 1
    for pi in 1:xgrid.panel.size-1
        if grid[init] < xgrid.panel[pi+1]
            head, tail = xgrid.inits[pi], xgrid.inits[pi] + xgrid.subgrids[pi].size - 1
            while grid[curr] < xgrid.panel[pi+1] && curr < length(grid)
                curr += 1
            end
            if grid[curr] < xgrid.panel[pi+1] && curr == length(grid)
                @assert xgrid.subgrids[pi].bound[1] <= grid[init] <= grid[curr] <= xgrid.subgrids[pi].bound[2]
                ff[init:curr] = interp1DGrid(view(data, head:tail), xgrid.subgrids[pi], view(grid, init:curr))
                return ff
            else
                @assert xgrid.subgrids[pi].bound[1] <= grid[init] <= grid[curr-1] <= xgrid.subgrids[pi].bound[2]
                ff[init:curr-1] = interp1DGrid(view(data, head:tail), xgrid.subgrids[pi], view(grid, init:curr-1))
            end
            # println(view(data, head:tail))
            # println(xgrid.subgrids[pi].grid)
            # println(grid[init:curr-1])
            # println(ff[init:curr-1])
            init = curr
        end
    end
    return ff
end

abstract type IntegrateStyle end
struct WeightIntegrate <: IntegrateStyle end
struct NoIntegrate <: IntegrateStyle end
struct CompositeIntegrate <: IntegrateStyle end
struct ChebIntegrate <: IntegrateStyle end
const WEIGHTINTEGRATE = WeightIntegrate()
const NOINTEGRATE = NoIntegrate()
const COMPOSITEINTEGRATE = CompositeIntegrate()
const CHEBINTEGRATE = ChebIntegrate()

IntegrateStyle(::Type) = NoIntegrate()
IntegrateStyle(::Type{<:SimpleG.BaryCheb}) = ChebIntegrate()
IntegrateStyle(::Type{<:SimpleG.GaussLegendre}) = WeightIntegrate()
IntegrateStyle(::Type{<:SimpleG.Uniform}) = WeightIntegrate()
IntegrateStyle(::Type{<:SimpleG.Arbitrary}) = WeightIntegrate()
IntegrateStyle(::Type{<:SimpleG.Log}) = WeightIntegrate()
IntegrateStyle(::Type{<:CompositeG.Composite}) = CompositeIntegrate()

"""
    function integrate1D(data, xgrid; axis=1)

calculate integration of data[i] on xgrid.
For 1D data, return a number; for multiple dimension, reduce the given axis.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- axis: axis to be integrated in data
"""
function integrate1D(data, xgrid::T; axis=1) where {T}
    return dropdims(mapslices(u -> integrate1D(IntegrateStyle(T), u, xgrid), data, dims=axis), dims=axis)
end

function integrate1D(data::AbstractMatrix, xgrid::T; axis=1) where {T}
    if axis == 1
        return map(u -> integrate1D(IntegrateStyle(T), u, xgrid), eachcol(data))
    elseif axis == 2
        return map(u -> integrate1D(IntegrateStyle(T), u, xgrid), eachrow(data))
    else
        throw(DomainError(axis, "axis should be 1 or 2 for Matrix"))
    end
end

function integrate1D(data::AbstractVector, xgrid::T; axis=1) where {T}
    return integrate1D(IntegrateStyle(T), data, xgrid)
end

"""
    function integrate1D(::NoIntegrate, data, xgrid)

calculate integration of data[i] on xgrid
works for grids that do not have integration weight stored

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
"""
function integrate1D(::NoIntegrate, data, xgrid)
    result = eltype(data)(0.0)

    grid = xgrid.grid
    for i in 1:xgrid.size
        if i == 1
            weight = 0.5 * (grid[2] - xgrid.bound[1])
        elseif i == xgrid.size
            weight = 0.5 * (xgrid.bound[2] - grid[end-1])
        else
            weight = 0.5 * (grid[i+1] - grid[i-1])
        end
        result += data[i] * weight
    end
    return result
end

"""
    function integrate1D(::WeightIntegrate, data, xgrid)

calculate integration of data[i] on xgrid
works for grids that have integration weight stored

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
"""
function integrate1D(::WeightIntegrate, data, xgrid)
    result = eltype(data)(0.0)

    for i in 1:xgrid.size
        result += data[i] * xgrid.weight[i]
    end
    return result
end

"""
    function integrate1D(::ChebIntegrate, data, xgrid)

calculate integration of data[i] on xgrid
works for grids that have integration weight stored

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
"""
function integrate1D(::ChebIntegrate, data, xgrid)
    a, b = xgrid.bound[1], xgrid.bound[2]
    return SimpleG.chebint(xgrid.size, -1.0, 1.0, data, xgrid.invVandermonde) * (b - a) / 2.0
end

"""
    function integrate1D(::CompositeIntegrate, data, xgrid)

calculate integration of data[i] on xgrid
call integrate1D for each subgrid and return the sum.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
"""
function integrate1D(::CompositeIntegrate, data, xgrid)
    result = eltype(data)(0.0)

    for pi in 1:xgrid.panel.size-1
        head, tail = xgrid.inits[pi], xgrid.inits[pi] + xgrid.subgrids[pi].size - 1
        result += integrate1D(view(data, head:tail), xgrid.subgrids[pi])
        currgrid = xgrid.subgrids[pi]
    end
    return result

end

"""
    function integrate1D(data, xgrid, range; axis=1)

calculate integration of data[i] on xgrid.
For 1D data, return a number; for multiple dimension, reduce the given axis.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- range: range of integration, [init, fin] within bound of xgrid.
- axis: axis to be integrated in data
"""
function integrate1D(data, xgrid::T, range; axis=1) where {T}
    return dropdims(mapslices(u -> integrate1D(IntegrateStyle(T), u, xgrid, range), data, dims=axis), dims=axis)
end

function integrate1D(data::AbstractMatrix, xgrid::T, range; axis=1) where {T}
    if axis == 1
        return map(u -> integrate1D(IntegrateStyle(T), u, xgrid, range), eachcol(data))
    elseif axis == 2
        return map(u -> integrate1D(IntegrateStyle(T), u, xgrid, range), eachrow(data))
    else
        throw(DomainError(axis, "axis should be 1 or 2 for Matrix"))
    end
end

function integrate1D(data::AbstractVector, xgrid::T, range; axis=1) where {T}
    return integrate1D(IntegrateStyle(T), data, xgrid, range)
end

@inline function trapezoidInt(data, grid, range)
    return (0.5 * (range[1] + range[2]) * (data[2] - data[1]) - grid[1] * data[2] + grid[2] * data[1]) / (grid[2] - grid[1]) * (range[2] - range[1])
end

function integrate1D(::NoIntegrate, data, xgrid, range)
    sign, x1, x2 = (range[1] < range[2]) ? (1.0, range[1], range[2]) : (-1.0, range[2], range[1])
    @assert xgrid.bound[1] <= x1 <= x2 <= xgrid.bound[2]
    xi1, xi2 = floor(xgrid, x1), floor(xgrid, x2)

    result = eltype(data)(0.0)

    grid = xgrid.grid
    # g[xi1], x1, g[xi1+1], g[xi1+2], ... , g[xi2], x2
    if xi1 == xi2
        result += trapezoidInt(view(data, xi1:xi1+1), view(grid, xi1:xi1+1), [x1, x2])
    else
        for i in xi1:xi2
            if i == xi1
                result += trapezoidInt(view(data, i:i+1), view(grid, i:i+1), [x1, grid[i+1]])
            elseif i == xi2
                result += trapezoidInt(view(data, i:i+1), view(grid, i:i+1), [grid[i], x2])
            else
                result += 0.5 * (data[i] + data[i+1]) * (grid[i+1] - grid[i])
            end
        end
    end
    return result * sign
end

function integrate1D(::ChebIntegrate, data, xgrid, range)
    a, b = xgrid.bound[1], xgrid.bound[2]
    x1, x2 = range[1], range[2]
    c1, c2 = (2x1 - a - b) / (b - a), (2x2 - a - b) / (b - a)
    return SimpleG.chebint(xgrid.size, c1, c2, data, xgrid.invVandermonde) * (b - a) / 2.0
end

function integrate1D(::WeightIntegrate, data, xgrid, range)
    return integrate1D(NoIntegrate(), data, xgrid, range)
    # may give bad result for GaussLegendre grid
    # sign, x1, x2 = (range[1]<range[2]) ? (1.0, range[1], range[2]) : (-1.0, range[2], range[1])
    # @assert xgrid.bound[1] <= x1 <= x2 <= xgrid.bound[2]
    # xi1, xi2 = floor(xgrid, x1), floor(xgrid, x2)

    # result = eltype(data)(0.0)
    # grid=xgrid.grid
    # for i in xi1+1:xi2
    #     if i==xi1+1
    #         weight = 0.5*(grid[i+1]-x1)
    #     elseif i==xi2
    #         weight = 0.5*(x2-grid[xi2])
    #     else
    #         weight = xgrid.weight[i]
    #     end
    #     result += data[i]*weight
    # end
    # return result*sign
end

function integrate1D(::CompositeIntegrate, data, xgrid, range)
    sign, x1, x2 = (range[1] < range[2]) ? (1.0, range[1], range[2]) : (-1.0, range[2], range[1])
    @assert xgrid.bound[1] <= x1 <= x2 <= xgrid.bound[2]
    pi1, pi2 = floor(xgrid.panel, x1), floor(xgrid.panel, x2)

    result = eltype(data)(0.0)
    if pi1 == pi2
        pi = pi1
        head, tail = xgrid.inits[pi], xgrid.inits[pi] + xgrid.subgrids[pi].size - 1
        result += integrate1D(view(data, head:tail), xgrid.subgrids[pi], [x1, x2])
    else
        for pi in pi1:pi2
            head, tail = xgrid.inits[pi], xgrid.inits[pi] + xgrid.subgrids[pi].size - 1
            if pi == pi1
                result += integrate1D(view(data, head:tail), xgrid.subgrids[pi], [x1, xgrid.subgrids[pi].bound[2]])
            elseif pi == pi2
                result += integrate1D(view(data, head:tail), xgrid.subgrids[pi], [xgrid.subgrids[pi].bound[1], x2])
            else
                result += integrate1D(view(data, head:tail), xgrid.subgrids[pi])
            end
        end
    end
    return result * sign
end

abstract type DifferentiateStyle end
struct NoDifferentiate <: DifferentiateStyle end
struct ChebDifferentiate <: DifferentiateStyle end
struct CompositeDifferentiate <: DifferentiateStyle end
const NODIFFERENTIATE = NoDifferentiate()
const COMPOSITEDIFFERENTIATE = CompositeDifferentiate()
const CHEBDIFFERENTIATE = ChebDifferentiate()

DifferentiateStyle(::Type) = NoDifferentiate()
DifferentiateStyle(::Type{<:SimpleG.BaryCheb}) = ChebDifferentiate()
DifferentiateStyle(::Type{<:CompositeG.Composite}) = CompositeDifferentiate()

"""
    function differentiate1D(data, xgrid, x; axis=1)

calculate integration of data[i] on xgrid.
For 1D data, return a number; for multiple dimension, reduce the given axis.

#Arguments:
- xgrid: one-dimensional grid of x
- data: one-dimensional array of data
- x: point to differentiate
- axis: axis to be differentiated in data
"""
function differentiate1D(data, xgrid::T, x; axis=1) where {T}
    return dropdims(mapslices(u -> differentiate1D(DifferentiateStyle(T), u, xgrid, x), data, dims=axis), dims=axis)
end

function differentiate1D(data::AbstractMatrix, xgrid::T, x; axis=1) where {T}
    if axis == 1
        return map(u -> differentiate1D(DifferentiateStyle(T), u, xgrid, x), eachcol(data))
    elseif axis == 2
        return map(u -> differentiate1D(DifferentiateStyle(T), u, xgrid, x), eachrow(data))
    else
        throw(DomainError(axis, "axis should be 1 or 2 for Matrix"))
    end
end

function differentiate1D(data::AbstractVector, xgrid::T, x; axis=1) where {T}
    return differentiate1D(DifferentiateStyle(T), data, xgrid, x)
end

function differentiate1D(::NoDifferentiate, data, xgrid, x)
    #simple numerical differentiate
    xi = floor(xgrid, x)
    return (data[xi+1] - data[xi]) / (xgrid[xi+1] - xgrid[xi])
end

function differentiate1D(::ChebDifferentiate, data, xgrid, x)
    a, b = xgrid.bound[1], xgrid.bound[2]
    c = (2x - a - b) / (b - a)
    return SimpleG.chebdiff(xgrid.size, c, data, xgrid.invVandermonde) / (b - a) * 2.0
end

function differentiate1D(::CompositeDifferentiate, data, xgrid, x)
    i = floor(xgrid.panel, x)
    head, tail = xgrid.inits[i], xgrid.inits[i] + xgrid.subgrids[i].size - 1
    return differentiate1D(view(data, head:tail), xgrid.subgrids[i], x)
end

end


text:
using HDF5

"""
`savedaqdata(h5, X::MeasData; kw...)`

Saves data read by a DAQ device in a HDF5 group

The group where the data will be stored is specified by the parameter `h5`.

The data is stored under the device name of the daq device (`devname(dev)`).

Other parameters related to the data acquisition can be provided using the key word
arguments `kw`. They will be stored under the keyword name.

The user should ensure that the data being stored (including attributes) have types
compatible with the HDF5.jl package.

"""
function savedaqdata(h5, X::MeasData; kw...)
    dname = devname(X)
    h5[dname] = measdata(X)
    d = h5[dname]
    attributes(d)["devname"] = dname
    attributes(d)["devtype"] = devtype(X)
    attributes(d)["time"] = AbstractDAQs.time2ms(meastime(X))
    attributes(d)["rate"] = samplingrate(X)
    attributes(d)["chans"] = collect(keys(X.chans))

    for (k,v) in kw
        attributes(d)[string(k)] = v
    end
    return
end

"""
`savedaqconfig(h5, dev::AbstractDAQ; kw...)`

Save device configuration. Stores information contained in the `conf::DAQConfig`
field of the daq device. Other information can be stored using the keyword parameters
where the keyword names are used as names in the HDF5 files.

The configuration is stored under the device name (`devname(dev)`). User should make
sure that different devices do not have the same name.

The user should ensure that the data being stored ( have types
compatible with the HDF5.jl package.

"""
function savedaqconfig(h5, dev::AbstractDAQ; kw...)

    device = string(typeof(dev))
    dname = devname(dev)
    g = create_group(h5, dname)
    g["device"] = device
    g["devname"] = dname
    g["ip"] = daqdevip(dev)
    g["model"] = daqdevmodel(dev)
    g["sn"] = daqdevserialnum(dev)
    g["tag"] = daqdevtag(dev)

    g["iparameters_names"] = collect(keys(iparameters(dev)))
    g["iparameters"] = collect(values(iparameters(dev)))

    g["sparameters_names"] = collect(keys(sparameters(dev)))
    g["sparameters"] = collect(values(sparameters(dev)))

    g["fparameters_names"] = collect(keys(fparameters(dev)))
    g["fparameters"] = collect(values(fparameters(dev)))

    g["daqchannels"] = daqchannels(dev)

    for (k,v) in kw
        g[string(k)] = v
    end

    return
end

function readdaqdata_dev(h, dtype, attr)

    if "devname" ∈ attr
        dname = read(attributes(h)["devname"])
    else
        dname = "unknown"
    end

    if "time" ∈ attr
        t = ms2time(read(attributes(h)["time"]))
    else
        t = now()
    end

    if "rate" ∈ attr
        rate = read(attributes(h)["rate"])
    elseif "fs" ∈ attr
        rate = read(attributes(h)["fs"])
    else
        rate = 1.0
    end

    data = read(h)

    nchans = size(data,1)

    if "chans" ∈ attr
        chans = read(attributes(h)["chans"])
    else
        chans = "C" .* string.(1:nchans)
    end

    ch = OrderedDict{String,Int}()
    for (i,c) in enumerate(chans)
        ch[c] = i
    end

    return MeasData(dname, dtype, t, rate, data, ch)
end

function readdaqdata_devset(h, dtype, attr)

    if !("devices" ∈ attr)
        error("There should be a 'devices' attribute in HDF5 file for it to be a valid MeasDataSet object")
    else
        devices = read(attributes(h)["devices"])
    end

    if "devname" ∈ attr
        dname = read(attributes(h)["devname"])
    else
        dname = "unknown"
    end

    if "time" ∈ attr
        t = ms2time(read(attributes(h)["time"]))
    else
        t = now()
    end

    data = OrderedDict{String,MeasData}()

    for dev in devices
        data[dev] = readdaqdata(h[dev])
    end

    return MeasDataSet(dname, dtype, t, data)
end

"""
`readdaqdata(h)`

Reads data acquisition data from a HDF5 file located at `h`. It will recognize
whether the data is a `MeasData` or `MeasDataSet` and act accordingly.

"""
function readdaqdata(h)

    attr = collect(keys(attributes(h)))

    if "devtype" ∈ attr
        dtype = read(attributes(h)["devtype"])
    elseif "type" ∈ attr
        dtype = read(attributes(h)["type"])
    else
        error("HDF5 section NOT and MeasData! No 'devtype' attribute!")
    end

    if dtype == "DeviceSet"
        return readdaqdata_devset(h, dtype, attr)
    else
        return readdaqdata_dev(h, dtype, attr)
    end

end

ext: jl


text:
module AngleBetweenVectors

import Base: angle

import LinearAlgebra: norm

Floats = Union{AbstractFloat, Complex{T} where T<:AbstractFloat}

@inline unitize(p) = p ./ norm(p)

"""
    angle(point1::T, point2::T) where {T}

Accurately ascertains the undirected angle (0 <= radians < pi)
between two points given in N-dimensional Cartesian coordinates.

Prefer this to `acos` alternatives
- more reliably accurate
- more consistently stable

Suggested when any |coordinate| of either point may be outside 2^±20 or [1/1_000_000, 1_000_000].
Strongly recommended when any |coordinate| is outside 2^±24 or [1/16_000_000, 16_000_000].

If either of the points is at the origin, the result is `NaN` because angle is undefined in that case.

You *must* define a tuple constructor `Tuple(x::YourPointType) = ...` if one does not already exist.
"""
angle(tuple1::NTuple{N,T}, tuple2::NTuple{N,T}) where {N, T<:Floats} = angle(real(T), tuple1, tuple2)

# because of the broadcasts .- and .+ below, it is essential to check size compatibility
# for arrays, whereas for tuples the consistency ensured by the "N" in the type
function angle(a1::AbstractArray{T}, a2::AbstractArray{T}) where {T<:Floats}
	size(a1) == size(a2) || throw(DimensionMismatch())
    return angle(real(T), a1, a2)
end

function angle(::Type{T}, point1, point2) where {T<:AbstractFloat}
    unitpoint1 = unitize(point1)
    unitpoint2 = unitize(point2)

    y = unitpoint1 .- unitpoint2
    x = unitpoint1 .+ unitpoint2

    a = 2 * atan(norm(y) / norm(x))

    !(signbit(a) || signbit(T(pi) - a)) ? a : (signbit(a) ? zero(T) : T(pi))
end

# this method allows the arrays to have different types, even non-float types
function angle(a1::AbstractArray{T1}, a2::AbstractArray{T2}) where {T1 <: Number, T2 <: Number}
    T = float(promote_type(T1, T2)) # the "T" for T(pi) must be a float type
    return angle(convert(AbstractArray{T}, a1), convert(AbstractArray{T}, a2))
end

@inline angle(point1::T, point2::T) where {T} = angle(Tuple(point1), Tuple(point2))

end # AngleBetweenVectors


text:
struct InputPair
    v :: Int32
    function InputPair(s::AbstractString)
        v = CoolProp.get_input_pair_index(s)
        new(v)
    end
end
Base.convert(::Type{Int32}, x::InputPair) = x.v

struct Property
    v :: Int32
    function Property(s::AbstractString)
        v = CoolProp.get_param_index(s)
        new(v)
    end
end
Base.convert(::Type{Int32}, x::Property) = x.v

struct WorkingFluid
    handle :: Int32
    function WorkingFluid(propertylibrary::AbstractString, fluid::AbstractString)
        handle = CoolProp.AbstractState_factory(propertylibrary, fluid)
        new(handle)
    end
end
Base.convert(::Type{Int32}, x::WorkingFluid) = x.handle


text:
module FluidFields

using InplaceRealFFT, LinearAlgebra, FFTW
using FluidTensors
import FluidTensors

export ScalarField, VectorField, SymTenField, SymTrTenField, isrealspace, setreal!, setfourier!, real!, fourier!, AbstractField, AbstractTensorField
export AntiSymTenField

export dealias!

const Float3264 = InplaceRealFFT.Float3264

include("kvec.jl")
include("scalar.jl")
include("vector.jl")
include("tensor.jl")

const AbstractField{T} = Union{<:ScalarField{T},<:VectorField{T},<:AbstractTensorField{T}}

include("functions.jl")

function __init__()
    path = @__DIR__
    fftwplan = joinpath(path,"fftw_wisdom")
    FFTW.set_num_threads(Threads.nthreads())
    if isfile(fftwplan)
        try
            FFTW.import_wisdom(fftwplan)
        catch
        end
    end
    atexit(()->(FFTW.export_wisdom(fftwplan)))
end

end # module


text:
include("../src/numder_module.jl")

module dtest
using PyPlot
PyPlot.pygui(true)
# using ..dev
using ..numder

function test_estimate_noise()
    # define function + number of evals
    fcn(x) = cos(x) + sin(x) + 1e-6 * rand() # sin(x) + 0.00001 * randn()
    n_eval = 7

    # initial point + step length
    x0 = 0.0 # 1 # pi/4
    h  = 0.01 # 1e3 * sqrt(eps(x0))
    # h = 1e3 .* sqrt.(eps.(x0))
    noise_est = numder.estimate_noise(
        fcn,
        x0;
        n_eval = n_eval)
end

function test_estimate_noise_nd()
    fcn(x) = x[1].^2 + x[2].^2 + 0.5 * randn()
    n_eval = 7

    x0 = [1.0, 1.0]

    noise_est = numder.estimate_noise(
            fcn,
            x0)
end

fcn(x) = sin(x) + 0.001 * randn()
fcn_theo(x) = sin(x)

x0 = 0.0
noise_est = numder.estimate_noise(fcn, x0)
h = noise_est

xs = collect( -pi/2:0.1:2*pi )
ys = fcn.(xs)
ys_theo = fcn_theo.(xs)
x_mids = 0.5 .* (xs[1:end-1] + xs[2:end])

dfdx = numder.numdiff_fw.(fcn, xs, sqrt(h) ) #, h
dfdx_theo = cos.(xs)

dfdx_from_signal = (ys[2:end] - ys[1:end-1]) ./ diff(xs)

fig = PyPlot.figure()

ax1 = PyPlot.subplot(2,1,1)
PyPlot.grid()
PyPlot.plot( xs, ys_theo )
PyPlot.plot( xs, ys )

ax1 = PyPlot.subplot(2,1,2, sharex = ax1)
PyPlot.grid()
PyPlot.plot( xs, dfdx_theo )
PyPlot.plot( xs, dfdx )
PyPlot.plot( x_mids, dfdx_from_signal )

end # dtest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module miraccess1
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using Test
function test()
    a = VecAttrib([1, 2, 4]);
    # @test show(a) != ""
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    # @show ir[1]
    # @show ir[1][3]
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test summary(ir.left) == "q4s = 6 x Q4"
    @test summary(ir.right) == "vrts = 12 x P1"
    @test summary(ir) == "(q4s, vrts): q4s = 6 x Q4, vrts = 12 x P1"
    # @show Q4
    # @show q4s
    true
end
end
using .miraccess1
miraccess1.test()

module mtattr10
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using Test
function test()
    a = VecAttrib([1, 2, 4]);
    # @test show(a) != ""
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test summary(ir.left) == "q4s = 6 x Q4"
    @test summary(ir.right) == "vrts = 12 x P1"
    @test summary(ir) == "(q4s, vrts): q4s = 6 x Q4, vrts = 12 x P1"
    # @show Q4
    # @show q4s
    true
end
end
using .mtattr10
mtattr10.test()

module mod1ame
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test summary(ir.left) == "q4s = 6 x Q4"
    @test summary(ir.right) == "vrts = 12 x P1"
    @test summary(ir) == "(q4s, vrts): q4s = 6 x Q4, vrts = 12 x P1"

    true
end
end
using .mod1ame
mod1ame.test()

module mtest1
using MeshCore: NoSuchShape, P1, L2, L3, SHAPE_DESC
using MeshCore: NoSuchShape, ShapeColl, manifdim, nfacets
using MeshCore: shapedesc, nshapes, manifdim, nvertices, facetdesc, nfacets, ridgedesc, nridges
using Test
function test()
    @test manifdim(P1) == 0
    @test nfacets(L2) == 2
    vrts = ShapeColl(P1, 12)
    lines = ShapeColl(L2, 32)
    @test shapedesc(vrts) == P1
    @test nshapes(vrts) == 12
    @test manifdim(vrts) == manifdim(P1)
    @test nvertices(lines) == nvertices(L2)
    @test facetdesc(lines) == P1
    @test nfacets(lines) == nfacets(L2)
    @test ridgedesc(lines) == NoSuchShape
    @test nridges(lines) == 0

    lines3 = ShapeColl(L3, 13)
    @test nvertices(lines3) == nvertices(L3)
    @test nvertices(lines3) == 3
    @test facetdesc(lines3) == P1
    @test nfacets(lines3) == nfacets(L3)
    @test ridgedesc(lines3) == NoSuchShape
    @test nridges(lines3) == 0
    @test manifdim(lines3) == 1

    @test SHAPE_DESC["P1"].name == "P1"
    @test SHAPE_DESC["L2"].name == "L2"
    @test SHAPE_DESC["T3"].name == "T3"
    @test SHAPE_DESC["T4"].name == "T4"
    @test SHAPE_DESC["Q4"].name == "Q4"
    @test SHAPE_DESC["H8"].name == "H8"
    @test SHAPE_DESC["L3"].name == "L3"
    true
end
end
using .mtest1
mtest1.test()

module mtattr11
using MeshCore: VecAttrib
using Test
function test()
    a = VecAttrib([1, 2, 4]);
    @test getindex(a, 3)     == 4
    setindex!(a, 6, 3)
    @test  getindex(a, 3) == 6
    @test a[3] == 6
    a[3] = -3
    @test a[3] == -3
    # @show a
    true
end
end
using .mtattr11
mtattr11.test()

module mtest1a1
using MeshCore: _sense
using Test
function test()
    oc = [4, 1, 5, 13]
    fc = [4, 1, 5, 13]
    @test _sense(fc, oc, 4) == +1
    fc = [13, 4, 1, 5]
    @test _sense(fc, oc, 4) == +1
    fc = [1, 5, 13, 4]
    @test _sense(fc, oc, 4) == +1
    fc = [1, 4, 13, 5]
    @test _sense(fc, oc, 4) == -1
    fc = [5, 1, 4, 13]
    @test _sense(fc, oc, 4) == -1
    fc = [13, 5, 1, 4]
    @test _sense(fc, oc, 4) == -1
    fc = [4, 13, 5, 1]
    @test _sense(fc, oc, 4) == -1
    true
end
end
using .mtest1a1
mtest1a1.test()

module mtest2
using StaticArrays
using MeshCore: VecAttrib
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    a =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])

    @test length(a) == 12
    x = a[SVector{2}([2, 4])]
    @test x[1] == SVector{2}([633.3333333333334 0.0])
    @test length(a[1]) == 2
    @test eltype(a[1]) == Float64

    true
end
end
using .mtest2
mtest2.test()

module mtest3
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test ir.name == "(q4s, vrts)"
    @test manifdim(q4s) == 2
    @test nvertices(q4s) == 4
    @test facetdesc(q4s) == L2
    @test nfacets(q4s) == 4
    @test nshapes(q4s) == 6
    @test manifdim(shapedesc(q4s)) == 2
    @test nvertices(shapedesc(q4s)) == 4
    @test nfacets(shapedesc(q4s)) == 4
    @test nridges(shapedesc(q4s)) == 4
    @test n1storderv(shapedesc(q4s)) == 4
    @test nshifts(shapedesc(q4s)) == 4
    @test ir[3, 3] == 7
    @test ir[1, 4] == 5
    @test ir[6, 1] == 7

    true
end
end
using .mtest3
mtest3.test()

module mtest5
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nrelations, facetconnectivity
using MeshCore: IncRel
using MeshCore: ir_skeleton
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)

    facemesh = ir_skeleton(mesh)
    @test nrelations(facemesh) == 17
    true
end
end
using .mtest5
mtest5.test()

module mtest6
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nrelations, facetconnectivity
using MeshCore: VecAttrib, IncRel
using MeshCore: ir_skeleton
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)

    facemesh = ir_skeleton(mesh)
    @test nrelations(facemesh) == 17
    # @show facemeshx
    for i in 1:nrelations(facemesh)
        x = locs[facemesh[i]]
    end #
    @test locs[facemesh[17]] == StaticArrays.SArray{Tuple{2},Float64,1,2}[[1266.6666666666667, 800.0], [1900.0, 800.0], ]
    @test locs[facemesh[17][1]] == [1266.6666666666667, 800.0]
    true
end
end
using .mtest6
mtest6.test()

module mtest7
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nrelations, facetconnectivity
using MeshCore: VecAttrib, IncRel
using MeshCore: ir_skeleton, ir_boundary, nshapes
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)

    bdrymesh = ir_boundary(mesh)
    @test nrelations(bdrymesh) == 10
    # @show bdrymesh.left, bdrymesh.right
    @test nshapes(bdrymesh.left) == 10
    @test nshapes(bdrymesh.right) == 12
    true
end
end
using .mtest7
mtest7.test()

module mtopoop1
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nshapes
using MeshCore: IncRel, ir_transpose, nshapes, VecAttrib
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)
    # Test the transpose incidence relation
    tmesh = ir_transpose(mesh)
    shouldget = Array{Int64,1}[[1], [1, 3], [3, 5], [5], [1, 2], [1, 2, 3, 4], [3, 4, 5, 6], [5, 6], [2], [2, 4], [4, 6], [6]]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (tmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    true
end
end
using .mtopoop1
mtopoop1.test()

module mtopoop2
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nshapes
using MeshCore: ir_bbyfacets, ir_skeleton, ir_transpose
using MeshCore: IncRel, VecAttrib
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)
    fmesh = ir_skeleton(mesh)
    # Test the transpose incidence relation
    tmesh = ir_bbyfacets(mesh, fmesh)
    shouldget = StaticArrays.SArray{Tuple{4},Int64,1,4}[
    [2, 4, 1, 8],
    [9, 11, 8, 15],
    [4, 6, 3, 10],
    [11, 13, 10, 16],
    [6, 7, 5, 12],
    [13, 14, 12, 17]   ]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (tmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    true
end
end
using .mtopoop2
mtopoop2.test()

module mtopoop3
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nshapes
using MeshCore: ir_bbyfacets, ir_skeleton, ir_transpose
using MeshCore: IncRel, VecAttrib
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
        locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)
    fmesh = ir_skeleton(mesh)
    # Test the bounded-by incidence relation
    bbmesh = ir_bbyfacets(mesh, fmesh)
    shouldget = StaticArrays.SArray{Tuple{4},Int64,1,4}[
    [2, 4, 1, 8],
    [9, 11, 8, 15],
    [4, 6, 3, 10],
    [11, 13, 10, 16],
    [6, 7, 5, 12],
    [13, 14, 12, 17]    ]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (bbmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    # Now test the transposed incidence relation
    tbbmesh = ir_transpose(bbmesh)
    shouldget = Array{Int64,1}[[1], [1], [3], [1, 3], [5], [3, 5], [5], [1, 2], [2], [3, 4], [2, 4], [5, 6],
    [4, 6], [6], [2], [4], [6]]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (tbbmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    true
end
end
using .mtopoop3
mtopoop3.test()

module mtopoop4
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nvertices, nfacets, facetdesc, nshapes
using MeshCore: ir_bbyridges, ir_skeleton, ir_transpose
using MeshCore: IncRel, VecAttrib
using Test
function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
        locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    mesh = IncRel(q4s, vrts, cc)
    # fmesh = skeleton(mesh)
    vmesh = IncRel(ShapeColl(P1, 12), vrts, [[i] for i in 1:nshapes(vrts)] )
    # Test the bounded-by incidence relation
    bbmesh = ir_bbyridges(mesh, vmesh)
    shouldget = Array{Int64,1}[[1, 2, 6, 5], [5, 6, 10, 9], [2, 3, 7, 6], [6, 7, 11, 10], [3, 4, 8, 7], [7, 8, 12, 11]]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (bbmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    # Now test the transposed incidence relation
    tbbmesh = ir_transpose(bbmesh)
    shouldget = Array{Int64,1}[[1], [1, 3], [3, 5], [5], [1, 2], [1, 2, 3, 4], [3, 4, 5, 6], [5, 6], [2], [2, 4], [4, 6], [6]]
    allmatch = true
    for j in 1:length(shouldget)
        for k in 1:length(shouldget[j])
            allmatch = allmatch && (tbbmesh[j, k] == shouldget[j][k])
        end
    end
    @test allmatch
    true
end
end
using .mtopoop4
mtopoop4.test()

include("samplet4.jl")

module mt4topo1
using StaticArrays
using MeshCore: P1, T4, ShapeColl,  manifdim, nvertices, nridges, nshapes
using MeshCore: ir_bbyridges, ir_skeleton, ir_bbyfacets, nshifts, _sense
using MeshCore: IncRel, ir_transpose, nrelations, nentities, VecAttrib
using ..samplet4: samplet4mesh
using Test
function test()
    xyz, cc = samplet4mesh()
    # Construct the initial incidence relation
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
    vrts = ShapeColl(P1, length(locs))
    tets = ShapeColl(T4, size(cc, 1))
    ir30 = IncRel(tets, vrts, cc)
    # Test the incidence relations 3 -> 0 & 0 -> 3
    ir03 = ir_transpose(ir30)
    allmatch = true
    for j in 1:nrelations(ir03)
        for k in 1:nentities(ir03, j)
            f = ir03[j, k]
            allmatch = allmatch && (j in ir30[f])
        end # k
    end # j
    @test allmatch

    # Test the incidence relations 2 -> 0 & 0 -> 2
    ir20 = ir_skeleton(ir30)
    ir02 = ir_transpose(ir20)
    allmatch = true
    for j in 1:nrelations(ir02)
        for k in 1:nentities(ir02, j)
            f = ir02[j, k]
            allmatch = allmatch && (j in ir20[f])
        end # k
    end # j
    @test allmatch

    # Test the incidence relations 1 -> 0 & 0 -> 1
    ir10 = ir_skeleton(ir20)
    ir01 = ir_transpose(ir10)
    allmatch = true
    for j in 1:nrelations(ir01)
        for k in 1:nentities(ir01, j)
            f = ir01[j, k]
            allmatch = allmatch && (j in ir10[f])
        end # k
    end # j
    @test allmatch

    # Test the incidence relations 3 -> 2 & 2 -> 3
    ir20 = ir_skeleton(ir30)
    ir32 = ir_bbyfacets(ir30, ir20)
    # Check that the tetrahedra connectivity refers to all three vertices of the faces
    allmatch = true
    @test nrelations(ir30) == nrelations(ir32)
    for j in 1:nrelations(ir32)
        for k in 1:nentities(ir32, j)
            f = ir32[j, k]
            for m in 1:3
                allmatch = allmatch && (ir20[abs(f), m] in ir30[j])
            end # m
        end # k
    end # j
    @test allmatch
    # Check that the orientations of the faces is correct
    allmatch = true
    @test nrelations(ir30) == nrelations(ir32)
    for j in 1:nrelations(ir32)
        for k in 1:nentities(ir32, j)
            for f in ir32[j]
                orientation = sign(f)
                f = abs(f)
                fc = ir20[f]
                matchone = false
                for m in 1:4
                    oc = ir30[j][ir30.left.shapedesc.facets[m, :]]
                    if length(intersect(fc, oc)) == 3
                        s = _sense(fc, oc, nshifts(ir20.left.shapedesc))
                        matchone = (s == orientation)
                    end
                    if matchone
                        break
                    end
                end # m
                allmatch = allmatch && matchone
            end
        end # k
    end # j
    @test allmatch
    # Check the transpose incidence relation
    ir23 = ir_transpose(ir32)
    allmatch = true
    for j in 1:nrelations(ir23)
        for k in 1:nentities(ir23, j)
            e = ir23[j, k]
            allmatch = allmatch && (j in abs.(ir32[e]))
        end # k
    end # j
    @test allmatch
    true
end
end
using .mt4topo1
mt4topo1.test()

module mtestattr1
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute
using Test
# using BenchmarkTools

function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])

    @test locs[10] == [633.3333333333334, 800.0]

    # @show typeof(Dict(a.name=>a))
    # @show typeof(Dict("geom"=>locs))
    vertices = ShapeColl(P1, size(xyz, 1))
    vertices.attributes["geom"] = locs
    locs = vertices.attributes["geom"]
    @test locs[10] == [633.3333333333334, 800.0]
    @test length(locs) == 12

    locs = attribute(vertices, "geom")
    @test locs[10] == [633.3333333333334, 800.0]

    true
end
end
using .mtestattr1
mtestattr1.test()

module m113651
using Test
using StaticArrays
using MeshCore: VecAttrib, ShapeColl, P1, attribute
function test()
	xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
	N, T = size(xyz, 2), eltype(xyz)
	locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])
	vertices = ShapeColl(P1, size(xyz, 1))
    vertices.attributes["geom"] = locs
	a = attribute(vertices, "geom")
	@test a[2] == [633.3333333333334, 0.0]

    dofn =  VecAttrib([i for i in 1:size(xyz, 1)])
    vertices.attributes["dofn"] = dofn
    # @show keys(vertices.attributes)
    @test "geom" in keys(vertices.attributes)
    @test "dofn" in keys(vertices.attributes)
end
end
using .m113651
m113651.test()

# module m1136193
# using Test
# using StaticArrays
# using MeshCore: VecAttrib, ShapeColl, P1, attribute
# function test()
# 	xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
# 	a = Attrib(i -> 1)
# 	vertices = ShapeColl(P1, size(xyz, 1), Dict("label"=>a))
# 	a = attribute(vertices, "label")
# 	@test a.co(3) == 1
# end
# end
# using .m1136193
# m1136193.test()

module mtest3a1
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test ir.name == "(q4s, vrts)"
    @test manifdim(q4s) == 2
    @test nvertices(q4s) == 4
    @test facetdesc(q4s) == L2
    @test nfacets(q4s) == 4
    @test nshapes(q4s) == 6
    @test manifdim(shapedesc(q4s)) == 2
    @test nvertices(shapedesc(q4s)) == 4
    @test nfacets(shapedesc(q4s)) == 4
    @test nridges(shapedesc(q4s)) == 4
    @test n1storderv(shapedesc(q4s)) == 4
    @test nshifts(shapedesc(q4s)) == 4
    @test ir[3, 3] == 7
    @test ir[1, 4] == 5
    @test ir[6, 1] == 7

    ir = IncRel(q4s, vrts, cc, "vertices")
    @test ir.name == "vertices"

    true
end
end
using .mtest3a1
mtest3a1.test()

module mtest3a2
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, ir_subset, nrelations
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    sir = ir_subset(ir, [1, 2])
    @test nrelations(sir) == 2
    @test isapprox(sir[1], [1, 2, 6, 5])

    true
end
end
using .mtest3a2
mtest3a2.test()

module mfeat1
using StaticArrays
using MeshCore: P1, L2, Q4, H8, T3, T4
using MeshCore: nfeatofdim
# using BenchmarkTools
using Test
function test()
    @test nfeatofdim(P1, 0) == 1
    @test nfeatofdim(P1, 1) == 0
    @test nfeatofdim(P1, 2) == 0
    @test nfeatofdim(P1, 3) == 0

    @test nfeatofdim(L2, 0) == 2
    @test nfeatofdim(L2, 1) == 1
    @test nfeatofdim(L2, 2) == 0
    @test nfeatofdim(L2, 3) == 0

    @test nfeatofdim(T3, 0) == 3
    @test nfeatofdim(T3, 1) == 3
    @test nfeatofdim(T3, 2) == 1
    @test nfeatofdim(T3, 3) == 0

    @test nfeatofdim(T4, 0) == 4
    @test nfeatofdim(T4, 1) == 6
    @test nfeatofdim(T4, 2) == 4
    @test nfeatofdim(T4, 3) == 1

    @test nfeatofdim(Q4, 0) == 4
    @test nfeatofdim(Q4, 1) == 4
    @test nfeatofdim(Q4, 2) == 1
    @test nfeatofdim(Q4, 3) == 0

    @test nfeatofdim(H8, 0) == 8
    @test nfeatofdim(H8, 1) == 12
    @test nfeatofdim(H8, 2) == 6
    @test nfeatofdim(H8, 3) == 1

    # @btime nfeatofdim($H8, 0)
    # @btime nfeatofdim($H8, 1)
    # @btime nfeatofdim($H8, 2)
    # @btime nfeatofdim($H8, 3)
    true
end
end
using .mfeat1
mfeat1.test()

module mtest4a5
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, ir_subset, nrelations, ir_identity
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    sir = ir_subset(ir, [1, 2])
    @test nrelations(sir) == 2
    @test isapprox(sir[1], [1, 2, 6, 5])
    ir1 = ir_identity(ir, :left)
    @test nrelations(ir1) == 6
    @test nshapes(ir1.left) == 6
    @test nshapes(ir1.right) == 6
    true
end
end
using .mtest4a5
mtest4a5.test()

module mtest1q
using MeshCore: NoSuchShape, P1, L3, T6, SHAPE_DESC
using MeshCore: NoSuchShape, ShapeColl, manifdim, nfacets
using MeshCore: shapedesc, nshapes, manifdim, nvertices, facetdesc, nfacets, ridgedesc, nridges, n1storderv
using Test
function test()
    @test manifdim(P1) == 0
    @test nfacets(L3) == 2
    @test nvertices(L3) == 3
    @test n1storderv(L3) == 2
    @test nfacets(T6) == 3

    vrts = ShapeColl(P1, 12)
    lines = ShapeColl(L3, 32)
    @test shapedesc(vrts) == P1
    @test nshapes(vrts) == 12
    @test manifdim(vrts) == manifdim(P1)
    @test nvertices(lines) == nvertices(L3)
    @test facetdesc(lines) == P1
    @test nfacets(lines) == nfacets(L3)
    @test ridgedesc(lines) == NoSuchShape
    @test nridges(lines) == 0

    tris6 = ShapeColl(T6, 13)
    @test nvertices(tris6) == nvertices(T6)
    @test nvertices(tris6) == 6
    @test n1storderv(T6) == 3
    @test facetdesc(tris6) == L3
    @test nfacets(tris6) == nfacets(T6)
    @test ridgedesc(tris6) == P1
    @test nridges(tris6) == 3
    @test manifdim(tris6) == 2

    @test SHAPE_DESC["P1"].name == "P1"
    @test SHAPE_DESC["L2"].name == "L2"
    @test SHAPE_DESC["T3"].name == "T3"
    @test SHAPE_DESC["T6"].name == "T6"
    @test SHAPE_DESC["T4"].name == "T4"
    @test SHAPE_DESC["Q4"].name == "Q4"
    @test SHAPE_DESC["H8"].name == "H8"
    @test SHAPE_DESC["L3"].name == "L3"
    true
end
end
using .mtest1q
mtest1q.test()

module mtest2q
using MeshCore: NoSuchShape, P1, L3, Q8, SHAPE_DESC
using MeshCore: NoSuchShape, ShapeColl, manifdim, nfacets
using MeshCore: shapedesc, nshapes, manifdim, nvertices, facetdesc, nfacets, ridgedesc, nridges, n1storderv
using Test
function test()
    @test manifdim(P1) == 0
    @test nfacets(L3) == 2
    @test nvertices(L3) == 3
    @test n1storderv(L3) == 2
    @test nfacets(Q8) == 4
    @test nridges(Q8) == 4

    vrts = ShapeColl(P1, 12)
    lines = ShapeColl(L3, 32)
    @test shapedesc(vrts) == P1
    @test nshapes(vrts) == 12
    @test manifdim(vrts) == manifdim(P1)
    @test nvertices(lines) == nvertices(L3)
    @test facetdesc(lines) == P1
    @test nfacets(lines) == nfacets(L3)
    @test ridgedesc(lines) == NoSuchShape
    @test nridges(lines) == 0

    q8s = ShapeColl(Q8, 13)
    @test nvertices(q8s) == nvertices(Q8)
    @test nvertices(q8s) == 8
    @test n1storderv(Q8) == 4
    @test facetdesc(q8s) == L3
    @test nfacets(q8s) == nfacets(Q8)
    @test ridgedesc(q8s) == P1
    @test nridges(q8s) == 4
    @test manifdim(q8s) == 2

    @test SHAPE_DESC["P1"].name == "P1"
    @test SHAPE_DESC["L2"].name == "L2"
    @test SHAPE_DESC["T3"].name == "T3"
    @test SHAPE_DESC["T6"].name == "T6"
    @test SHAPE_DESC["T4"].name == "T4"
    @test SHAPE_DESC["Q4"].name == "Q4"
    @test SHAPE_DESC["Q8"].name == "Q8"
    @test SHAPE_DESC["H8"].name == "H8"
    @test SHAPE_DESC["L3"].name == "L3"
    true
end
end
using .mtest2q
mtest2q.test()

module mtestbb1
using MeshCore: @_check
using Test
function test()
    @_check 1 == 1 "Really?"
    @_check 1 == 1 "Why " * "not?"
    a = 13
    b = 13
    @_check a == b "Why is $a " * "not equal to $b?"
    a = 61
    @_check a > b "Why is $a " * "not greater than $b?"
    Test.@test_throws AssertionError @_check(a < b, "Why is $a " * "not greater than $b?")
    true
end
end
using .mtestbb1
mtestbb1.test()

module mtestattr6
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute
using Test
# using BenchmarkTools

function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])

    @test locs[10] == [633.3333333333334, 800.0]
    locs[10] = -locs[10]
    @test locs[10] == [-633.3333333333334, -800.0]
    @test isapprox(locs[8:10], StaticArrays.SArray{Tuple{2},Float64,1,2}[[1900.0, 400.0], [0.0, 800.0], [-633.3333333333334, -800.0]])
    locs[8:10] .= StaticArrays.SArray{Tuple{2},Float64,1,2}[[1.0, 1.0], [1.0, 1.0], [1.0, 1.0]]
    @test isapprox(locs[8:10], StaticArrays.SArray{Tuple{2},Float64,1,2}[[1.0, 1.0], [1.0, 1.0], [1.0, 1.0]])

    true
end
end
using .mtestattr6
mtestattr6.test()

module mtesat4a5
using StaticArrays
using MeshCore: P1, L2, Q4, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel, ir_subset, nrelations, ir_identity, indextype, ir_code
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    @test indextype(ir) == Int64
    @test ir_code(ir) == (2, 0)
    sir = ir_subset(ir, [1, 2])
    @test nrelations(sir) == 2
    @test isapprox(sir[1], [1, 2, 6, 5])
    ir1 = ir_identity(ir, :left)
    @test nrelations(ir1) == 6
    @test nshapes(ir1.left) == 6
    @test nshapes(ir1.right) == 6
    true
end
end
using .mtesat4a5
mtesat4a5.test()

module mfeatx1
using StaticArrays
using MeshCore: P1, L2, Q4, H8, T3, T4, T6, L3
using MeshCore: nfeatofdim
# using BenchmarkTools
using Test
function test()

    @test nfeatofdim(L3, 0) == 3
    @test nfeatofdim(L3, 1) == 1
    @test nfeatofdim(L3, 2) == 0
    @test nfeatofdim(L3, 3) == 0

    @test nfeatofdim(T6, 0) == 6
    @test nfeatofdim(T6, 1) == 3
    @test nfeatofdim(T6, 2) == 1
    @test nfeatofdim(T6, 3) == 0

    true
end
end
using .mfeatx1
mfeatx1.test()

module mtest3y1
using StaticArrays
using MeshCore.Exports
using Test
function test()
    c = [(1, 2, 6, 5), (5, 6, 10, 9), (2, 3, 7, 6), (6, 7, 11, 10), (3, 4, 8, 7), (7, 8, 12, 11)]
    cc = [SVector{nvertices(Q4)}(c[idx]) for idx in 1:length(c)]
    q4s = ShapeColl(Q4, 6)
    vrts = ShapeColl(P1, 12)
    ir = IncRel(q4s, vrts, cc)
    q4s = ShapeColl(Q4, 6, "q4s")
    vrts = ShapeColl(P1, 12, "vrts")
    ir = IncRel(q4s, vrts, cc)
    @test ir.name == "(q4s, vrts)"
    @test manifdim(q4s) == 2
    @test nvertices(q4s) == 4
    @test facetdesc(q4s) == L2
    @test nfacets(q4s) == 4
    @test nshapes(q4s) == 6
    @test manifdim(shapedesc(q4s)) == 2
    @test nvertices(shapedesc(q4s)) == 4
    @test nfacets(shapedesc(q4s)) == 4
    @test nridges(shapedesc(q4s)) == 4
    @test n1storderv(shapedesc(q4s)) == 4
    @test nshifts(shapedesc(q4s)) == 4
    @test ir[3, 3] == 7
    @test ir[1, 4] == 5
    @test ir[6, 1] == 7

    true
end
end
using .mtest3y1
mtest3y1.test()

module mtest1y2
using MeshCore.Exports
using Test
function test()
    @test manifdim(P1) == 0
    @test nfacets(L2) == 2
    vrts = ShapeColl(P1, 12)
    lines = ShapeColl(L2, 32)
    @test shapedesc(vrts) == P1
    @test nshapes(vrts) == 12
    @test manifdim(vrts) == manifdim(P1)
    @test nvertices(lines) == nvertices(L2)
    @test facetdesc(lines) == P1
    @test nfacets(lines) == nfacets(L2)
    @test ridgedesc(lines) == NoSuchShape
    @test nridges(lines) == 0

    lines3 = ShapeColl(L3, 13)
    @test nvertices(lines3) == nvertices(L3)
    @test nvertices(lines3) == 3
    @test facetdesc(lines3) == P1
    @test nfacets(lines3) == nfacets(L3)
    @test ridgedesc(lines3) == NoSuchShape
    @test nridges(lines3) == 0
    @test manifdim(lines3) == 1

    @test SHAPE_DESC["P1"].name == "P1"
    @test SHAPE_DESC["L2"].name == "L2"
    @test SHAPE_DESC["T3"].name == "T3"
    @test SHAPE_DESC["T4"].name == "T4"
    @test SHAPE_DESC["Q4"].name == "Q4"
    @test SHAPE_DESC["H8"].name == "H8"
    @test SHAPE_DESC["L3"].name == "L3"
    true
end
end
using .mtest1y2
mtest1y2.test()

module mtestattrbc1
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute, FunAttrib, datavaluetype
using Test
# using BenchmarkTools

function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0; 1900.0 0.0; 0.0 400.0; 633.3333333333334 400.0; 1266.6666666666667 400.0; 1900.0 400.0; 0.0 800.0; 633.3333333333334 800.0; 1266.6666666666667 800.0; 1900.0 800.0]
    N, T = size(xyz, 2), eltype(xyz)
    locs =  VecAttrib([SVector{N, T}(xyz[i, :]) for i in 1:size(xyz, 1)])

    @test locs[10] == [633.3333333333334, 800.0]

    # @show typeof(Dict(a.name=>a))
    # @show typeof(Dict("geom"=>locs))
    vertices = ShapeColl(P1, size(xyz, 1))
    vertices.attributes["geom"] = locs
    locs = vertices.attributes["geom"]
    @test locs[10] == [633.3333333333334, 800.0]
    @test length(locs) == 12

    locs = attribute(vertices, "geom")
    @test locs[10] == [633.3333333333334, 800.0]
    @test locs[[3, 5]] == StaticArrays.SArray{Tuple{2},Float64,1,2}[[1266.6666666666667, 0.0], [0.0, 400.0]]

    @test datavaluetype(locs) == typeof(SVector{N, T}([0.0 0.0]))
    @test size(locs) == (size(xyz, 1),)

    flocs = FunAttrib(SVector{N, T}([0.0 0.0]), size(xyz, 1), i -> SVector{N, T}(xyz[i, :]))
    vertices.attributes["geom"] = flocs
    flocs = vertices.attributes["geom"]
    @test flocs[10] == [633.3333333333334, 800.0]

    @test datavaluetype(flocs) == typeof(SVector{N, T}([0.0 0.0]))
    @test size(flocs) == (size(xyz, 1),)
    @test flocs[[3, 5]] == StaticArrays.SArray{Tuple{2},Float64,1,2}[[1266.6666666666667, 0.0], [0.0, 400.0]]

    true
end
end
using .mtestattrbc1
mtestattrbc1.test()

module mtestattrbc2
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute, FunAttrib, datavaluetype
using Test
# using BenchmarkTools

function test()
    list = [3, 7, 1, 2, 5, 4, 6]
    N, T = size(list, 2), eltype(list)
    la =  VecAttrib(list)

    @test la[4] == 2

    la = FunAttrib(0, size(list, 1), i -> list[i])
    @test la[4] == 2

    list = collect(1:7)
    N, T = size(list, 2), eltype(list)
    la =  VecAttrib(list)

    @test la[4] == 4

    la = FunAttrib(0, size(list, 1), i -> i)
    @test la[4] == 4

    true
end
end
using .mtestattrbc2
mtestattrbc2.test()

module mtestattrbc3
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute, FunAttrib, datavaluetype
using Test
# using BenchmarkTools

function test()
    list = [3, 7, 1, 2, 5, 4, 6]
    N, T = size(list, 2), eltype(list)
    la =  VecAttrib(list)

    @test la[4] == 2

    la = FunAttrib(0, size(list, 1), i -> list[i])
    @test la[4] == 2
    @test_throws ErrorException la[4] = -4
    @test la[4] == 2

    list = collect(1:7)
    N, T = size(list, 2), eltype(list)
    la =  VecAttrib(list)

    @test la[4] == 4

    la = FunAttrib(0, size(list, 1), i -> i)
    @test la[4] == 4
    @test_throws ErrorException la[4] = -4
    @test la[4] == 4

    true
end
end
using .mtestattrbc3
mtestattrbc3.test()

module mtestattrbc4
using StaticArrays
using MeshCore: VecAttrib, P1, ShapeColl, attribute, FunAttrib, datavaluetype
using Test
# using BenchmarkTools

function test()
    xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
    a = FunAttrib(0, size(xyz, 1), i -> 1)
    vertices = ShapeColl(P1, size(xyz, 1))
    vertices.attributes["label"] = a
    @test a[3] == 1

    true
end
end
using .mtestattrbc4
mtestattrbc4.test()

module mfunattraccess1
using StaticArrays
using MeshCore: VecAttrib
using MeshCore: P1, L2, H8, ShapeColl, manifdim, nfacets, facetdesc, nshapes
using MeshCore: Q4ShapeDesc, shapedesc, n1storderv, nridges, nshifts, nvertices
using MeshCore: IncRel
using MeshCore: FunAttrib
using Test
function test()
    L, W, H = 1.5*pi, 3.0, 1.5
    R = 4.0
    M, N, K = 13, 4, 2
    linix = LinearIndices((M+1, N+1, K+1))
    sn(m, n, k) = linix[m, n, k]

    C = SVector{8, Int}[]
    for k in 1:K, n in 1:N, m in 1:M
        c = (
        sn(m, n, k), sn(m+1, n, k), sn(m+1, n+1, k), sn(m, n+1, k),
        sn(m, n, k+1), sn(m+1, n, k+1), sn(m+1, n+1, k+1), sn(m, n+1, k+1)
        )
        push!(C, SVector{8, Int}(c))
    end
    loc = SVector{3, Float64}[]
    for k in 1:K+1, n in 1:N+1, m in 1:M+1
        a = (m-1)/M*L
        x, y, z = (R+(n-1)/N*W) * cos(a), (R+(n-1)/N*W) * sin(a), (k-1)/K*H * (1 + a)
        push!(loc, SVector{3, Float64}(x, y, z))
    end

    elements = ShapeColl(H8, length(C), "elements")
    vertices = ShapeColl(P1, length(loc), "vertices")

    vertices.attributes["geom"] = VecAttrib(loc)

    connectivities = IncRel(elements, vertices, C)

    cart = CartesianIndices((M+1, N+1, K+1))
    access_location(i) = let
        ci  = cart[i]
        m, n, k = ci[1], ci[2], ci[3]
        a = (m-1)/M*L
        x, y, z = (R+(n-1)/N*W) * cos(a), (R+(n-1)/N*W) * sin(a), (k-1)/K*H * (1 + a)
        SVector{3, Float64}(x, y, z)
    end

    vertices2 = ShapeColl(P1, length(loc), "vertices2")
    vertices2.attributes["geom"] = FunAttrib(0.0, prod((M+1, N+1, K+1)), access_location);

    @test size(vertices.attributes["geom"]) == size(vertices2.attributes["geom"])

    n = 0.0
    for i in 1:size(vertices2.attributes["geom"], 1)
        n += sum((vertices.attributes["geom"][i] - vertices2.attributes["geom"][i]).^2)
    end
    @test n ≈ 0

    true
end
end
using .mfunattraccess1
mfunattraccess1.test()


text:
using MAT
using LinearAlgebra
using SparseArrays
using MatrixNetworks
include("../src/Graph_Louvain.jl")
mat = matread("../data/emailEUcore.mat")
A = mat["A"]            # Largest connected component of email EU
ASim = mat["Asim"]      # Undirected version of same graph: ASim = spones(A+A')
m = sum(A)
label = vec(round.(Int64,mat["labels"]))

## Performs weighted triangle clique expansion on an already undirected graph
include("../src/Motif_Adj_Mat.jl")
Ata = ExtractTriMat(ASim)

## The clique expanded graph isn't connected, so extract the largest connected component
Acc, p = largest_component(Ata)
core_inds = findall(x->x==true,p)
matwrite("../data/EmailTriCore.mat", Dict("core_inds"=>core_inds))
Ata = Ata[p,p]
ASim = ASim[p,p]    # Symmetric version of the graph
label = label[p]

matwrite("../data/EmailTriangle.mat",Dict("Ata"=>Ata))
n = size(ASim,1)

# Get total weight of graphs
mta = sum(Ata)
mSim = sum(A)

# Normalize so that total weight is the same for both graphs
ASim = mta/mSim*ASim

## Get weighted degree vectors
dta = vec(sum(Ata,dims = 2))
dSim = vec(sum(ASim,dims = 2))

numpoints = 20
Lams = range(0.05/m,3/m,length=numpoints)

# Store clusterings
SimC  = zeros(n,numpoints)
taC = zeros(n,numpoints)    # triangle motifs

##
maxits = 1000           # Louvain maximum iterations
numtimes = 1            # Number of times to run Lambda-Lovain for each
totaltimes = zeros(length(Lams),2)
for i = 1:length(Lams)
    lam = Lams[i]

    println("Lambda = $lam")
    start = time()
    cSim, BestObj = Many_Louvain(ASim,dSim,lam,numtimes,maxits)
    totaltimes[i,1] = time()-start
    SimC[:,i] = cSim

    start = time()
    cta, BestObj = Many_Louvain(Ata,dta,lam,numtimes,maxits)
    totaltimes[i,2] = time()-start
    taC[:,i] = cta
end

matwrite("Output/Louvain_Email_clusterings.mat", Dict("taC" => taC,
    "SimC" => SimC,"Lams" => collect(Lams),"totaltimes"=>totaltimes))

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
function base_to_conservative(Qbase, cellxmax, cellymax, specific_heat_ratio, nval)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    Qcon = zeros(cellxmax, cellymax, nval)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qcon[i,j,1] = Qbase[i,j,1]
            Qcon[i,j,2] = Qbase[i,j,1]*Qbase[i,j,2]
            Qcon[i,j,3] = Qbase[i,j,1]*Qbase[i,j,3]
            Qcon[i,j,4] = Qbase[i,j,4]/(specific_heat_ratio-1)+Qbase[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2
        end
    end
    return Qcon
end

function conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval)
    """
    Qbase=[rho,u,v,p]
    Qcon=[rho,rhou,rhov,e]
    """
    Qbase = zeros(cellxmax, cellymax, nval)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            Qbase[i,j,1] = Qcon[i,j,1]
            Qbase[i,j,2] = Qcon[i,j,2]/Qcon[i,j,1]
            Qbase[i,j,3] = Qcon[i,j,3]/Qcon[i,j,1]
            Qbase[i,j,4] = (Qcon[i,j,4]-Qcon[i,j,1]*(Qbase[i,j,2]^2+Qbase[i,j,3]^2)/2)*(specific_heat_ratio-1)
        end
    end
    return Qbase
end

function setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval)
    Qcon_hat = zeros(cellxmax, cellymax, nval)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            for l in 1:nval
                Qcon_hat[i,j,l] = Qcon[i,j,l] * volume[i,j]
            end
        end
    end
    return Qcon_hat
end

function Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval)
    Qcon = zeros(cellxmax, cellymax, nval)
    Threads.@threads for i in 1:cellxmax
        for j in 1:cellymax
            for l in 1:nval
                Qcon[i,j,l] = Qcon_hat[i,j,l] / volume[i,j]
            end
        end
    end
    return Qcon
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using DelimitedFiles

function storycloze_init()
    register(DataDep(
        "StoryClozeTest",
        """
        Story Cloze Test and ROCStories Corpora

        http://cs.rochester.edu/nlp/rocstories/

        Please fill in the form in the url,
        the dataset owner will send the download link to you with email.
        please copy-paste the link after you get the link.

        """,
        "http://cs.rochester.edu/nlp/rocstories/",
        "27eb558d129b922ddb7d5427ee140a21d9d636fbd46c90343fd800ab533f9cb6";
        fetch_method= (dummyremote, localdir) -> begin
            println(dummyremote)
            println("Please paste link you get after recieve the email")
            print("ROCStories winter 2017 set: ")
            rocw_url = readline()
            print("ROCStories spring 2016 set: ")
            rocs_url = readline()
            println("Story Cloze Test Spring 2016 set:")
            print("* val set: ")
            valset_url = readline()
            print("* test set: ")
            testset_url = readline()
            gdownload(rocw_url, localdir)
            gdownload(rocs_url, localdir)
            gdownload(valset_url, localdir)
            gdownload(testset_url, localdir)
        end
    ))
end

struct StoryCloze <: Dataset end

function testfile(::StoryCloze)
    sets, headers = readdlm(datadep"StoryClozeTest/cloze_test_test__spring2016-cloze_test_ALL_test.csv", ',', String; header=true)

    [selectdim(sets, 2, i) for i = 2:8]
end

function trainfile(::StoryCloze)
    sets, headers = readdlm(datadep"StoryClozeTest/cloze_test_val__spring2016-cloze_test_ALL_val.csv", ',', String; header=true)

    [selectdim(sets, 2, i) for i = 2:8]
end

get_labels(::StoryCloze) = ("1", "2")


text:
using TestPackage2

TestPackage2.greet()


text:
export DisplayArray, render

struct DisplayArray
    size::Tuple{Int, Int}
    background::Char
    content::Vector{DisplayRow}
end

function DisplayArray(h, w; background=' ')
    (textwidth(background) != 1) && (throw(
        "Invalid background character. " *
        "Text width of background character must be 1."
    ))

    content = Array{DisplayRow}(undef, h)
    foreach(i->content[i]=DisplayRow(w, background=background), 1:h)

    return DisplayArray(
        (h, w),
        background,
        content
    )
end

function Base.show(io::IO, row::DisplayArray)
    Base.print(io, "DisplayArray(" *
        "size=$(size(row)), " *
        "background char=Char($(convert(UInt16, row.background)))" *
        ")")
end

Base.size(array::DisplayArray) = array.size

Base.lastindex(array::DisplayArray, d::Int) = size(array)[d]

function Base.setindex!(
    array::DisplayArray,
    c::Char,
    display_row::Int, display_col::Int
)
    array.content[display_row][display_col] = c
end

function Base.setindex!(
    array::DisplayArray,
    str::AbstractString,
    display_row::Int, display_col_range::UnitRange{Int64}
)
    array.content[display_row][display_col_range] = str
end

function Base.setindex!(
    array::DisplayArray,
    str::AbstractString,
    display_row_range::UnitRange{Int64}, display_col_range::UnitRange{Int64}
)
    str = padding_vertical(str, length(display_row_range))
    split_str = split(str, '\n')

    for (i, s) in zip(display_row_range, split_str)
        array[i, display_col_range] = s
    end
end

function Base.setindex!(
    array::DisplayArray,
    str::AbstractString,
    display_row_range::UnitRange{Int64}, display_col::Int
)
    str = padding_vertical(str, length(display_row_range))
    split_str = split(str, '\n')

    for (i, s) in zip(display_row_range, split_str)
        (length(s) > 1) && throw(BoundsError)
        array[i, display_col:(display_col+textwidth(s))] = s
    end
end

function render(stream::IO, array::DisplayArray; pos=(-1, -1))
    T.buffered(stream, array, pos) do buffer, array, pos
        (pos != (-1, -1)) && T.cmove(buffer, pos[1], pos[2])
        foreach(
            row->(
                T.csave(buffer);
                render(buffer, row);
                T.crestore(buffer);
                T.cmove_down(buffer)
            ),
            array.content
        )
    end
end

render(array::DisplayArray; pos=(-1, -1)) = render(T.out_stream, array, pos=pos)


text:
module RandTest

using QuasiCopula, Random, Statistics, Test, LinearAlgebra, StatsFuns

# Normal
@testset "Normal(0,1) * (1 + 0.5 x^2)" begin ## where x = mean(dist)
dist = Normal(0, 1)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_normal = pdf_constants(γ, dist)

# test if the CDF is a proper density (integrates to 1)
inv(d_normal.c0 + d_normal.c1 * d_normal.μ + d_normal.c2 * (d_normal.σ2 + abs2(d_normal.μ)))
@test d_normal.c * (d_normal.c0 + d_normal.c1  * d_normal.μ + d_normal.c2 * (d_normal.σ2 + abs2(d_normal.μ))) == 1.0
@test minimum(d_normal) == -Inf
@test maximum(d_normal) == Inf
@test cdf(d_normal, -Inf) == 0
@test cdf(d_normal, Inf) == 1

Random.seed!(1234)
nsample = 10_000
@info "sample $nsample points for the $dist distribution"
s = Vector{Float64}(undef, nsample)
rand!(d_normal, s) # compile
@time rand!(d_normal, s)
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_normal))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_normal))")
end

# GAMMA
# when α = 2.5, 1.0 = θ
@testset "Gamma(α = 2.5, θ = 1.0) * (c0 + c1 * x + c2 * x^2)" begin ## where x = mean(dist)
Random.seed!(1234)
α, θ = (2.5, 1.0)
# res = Vector{Float64}(undef, 5)
# Γ = rand(5, 5)
γ = 1.0
dist = Gamma(α, θ)

# test if the proper c0, c1, c2 constants are stored.
d_gamma = pdf_constants(γ, dist)
@test d_gamma.c0 == 1 + (mean(d_gamma.d)^2 * inv(var(d_gamma.d)) * 0.5 * γ)
@test d_gamma.c1 == 0.5 * γ * (-2 * mean(d_gamma.d) * inv(var(d_gamma.d)))
@test d_gamma.c2 == 0.5 * γ * (inv(var(d_gamma.d)))

# test if the CDF is a proper density (integrates to 1)
α, θ = Distributions.params(dist)
normalizing_c1 = (StatsFuns.gamma(α + 1)/ StatsFuns.gamma(α))
normalizing_c2 = (StatsFuns.gamma(α + 2)/ StatsFuns.gamma(α))
@test d_gamma.c * (d_gamma.c0 + normalizing_c1 * d_gamma.c1 + normalizing_c2 * d_gamma.c2) == 1.0

# @test d_gamma.c ≈ inv(1+ 0.5 * tr(Γ))
@test minimum(d_gamma) == 0
@test maximum(d_gamma) == Inf
@test cdf(d_gamma, 0) == 0
@test cdf(d_gamma, Inf) == 1

###
Random.seed!(1234)
nsample = 10_000 #
@info "sample $nsample points for the $dist distribution"
s = Vector{Float64}(undef, nsample)
rand!(d_gamma, s) # compile
@time rand!(d_gamma, s)
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_gamma))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_gamma))")
end

# n = 10,0000 ::: time = 0.028 seconds
# sample mean = 3.1510221305244226; theoretical mean = 3.1666666666666665
# sample var = 5.67119904573935; theoretical var = 5.722222222222223

# ## when there is a lot of noise there is some deviation from truth
# @testset "Gamma(α = 3.0, θ = 2.0) * (c0 + c1 * x + c2 * x^2); Lots of nosie ==> deviation" begin ## where x = mean(dist)
# α, θ = (3.0, 2.0)
# Random.seed!(1234)
# dist = Gamma(α, θ)

# # test if the proper c0, c1, c2 constants are stored.
# d_gamma = pdf_constants(γ, dist)
# @test d_gamma.c0 == 1 + (mean(d_gamma.d)^2 * inv(var(d_gamma.d)) * 0.5 * γ)
# @test d_gamma.c1 == 0.5 * γ * (-2 * mean(d_gamma.d) * inv(var(d_gamma.d)))
# @test d_gamma.c2 == 0.5 * γ * (inv(var(d_gamma.d)))

# # test if the CDF is a proper density (integrates to 1)
# α, θ = params(dist)
# normalizing_c1 = (StatsFuns.gamma(α + 1)/ StatsFuns.gamma(α))
# normalizing_c2 = (StatsFuns.gamma(α + 2)/ StatsFuns.gamma(α))
# @test d_gamma.c * (d_gamma.c0 + d_gamma.c1  * d_gamma.μ + d_gamma.c2 * (d_gamma.σ2 + abs2(d_gamma.μ))) == 1.0

# @test d_gamma.c ≈ inv(1+ 0.5 * γ)
# @test minimum(d_gamma) == 0
# @test maximum(d_gamma) == Inf
# @test cdf(d_gamma, 0) == 0
# @test cdf(d_gamma, Inf) == 1

# ###
# Random.seed!(1234)
# nsample = 10_000 #
# @info "sample $nsample points for the $dist distribution"
# s = Vector{Float64}(undef, nsample)
# rand!(d_gamma, s) # compile
# @time rand!(d_gamma, s)
# println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_gamma))")
# println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_gamma))")
# end
# # n = 10,0000 ::: time = 0.049717 seconds
# # sample mean = 5.663051756167577; theoretical mean = 6.632654124304445
# # sample var = 11.845482796775215; theoretical var = 19.191598250653946

### EXPONENTIAL test
@testset "Exponential(θ = 3) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(1234)
dist = Exponential(3)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_exp = pdf_constants(γ, dist)
@test d_exp.c0 == 1  + (mean(d_exp.d)^2 * inv(var(d_exp.d)) * 0.5 * γ)
@test d_exp.c1 == 0.5 * γ * (-2 * mean(d_exp.d) * inv(var(d_exp.d)))
@test d_exp.c2 == 0.5 * γ * (inv(var(d_exp.d)))

# test if the CDF is a proper density (integrates to 1)
θ = Distributions.params(d_exp.d)[1]
normalizing_c1 = θ * (StatsFuns.gamma(2)/ StatsFuns.gamma(1))
normalizing_c2 = θ^2 * (StatsFuns.gamma(3)/ StatsFuns.gamma(1))
@test d_exp.c * (d_exp.c0 + normalizing_c1 * d_exp.c1 + normalizing_c2 * d_exp.c2) == 1.0

@test d_exp.c ≈ inv(1+ 0.5 * γ)
@test minimum(d_exp) == 0
@test maximum(d_exp) == Inf
@test cdf(d_exp, 0) == 0
@test cdf(d_exp, Inf) == 1

###
Random.seed!(1234)
nsample = 10_000 #_000
@info "sample $nsample points for the $dist distribution"
s = Vector{Float64}(undef, nsample)
rand!(d_exp, s) # compile
@time rand!(d_exp, s)
println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_exp))")
println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_exp))")
end
# n = 10,0000 ::: time = 0.011543 seconds
# sample mean = 4.9647655815252465; theoretical mean = 5.0
# sample var = 28.68846969375196; theoretical var = 29.0

# ### EXPONENTIAL test
# @testset "Exponential(θ = 20) * (c0 + c1 * x + c2 * x^2);" begin
# Random.seed!(1234)
# Γ = rand(5, 5)
# res = Vector{Float64}(undef, 5)
# dist = Exponential(20)

# # test if the proper c0, c1, c2 constants are stored.
# d_exp = pdf_constants(Γ, res, 1, dist)
# @test d_exp.c0 == 1 + 0.5tr(Γ[2:end, 2:end]) + (mean(d_exp.d)^2 * inv(var(d_exp.d)) * 0.5 * Γ[1,1])
# @test d_exp.c1 == 0.5 * Γ[1, 1] * (-2 * mean(d_exp.d) * inv(var(d_exp.d)))
# @test d_exp.c2 == 0.5 * Γ[1, 1] * (inv(var(d_exp.d)))

# # test if the CDF is a proper density (integrates to 1)
# θ = params(d_exp.d)[1]
# normalizing_c1 = θ * (StatsFuns.gamma(2)/ StatsFuns.gamma(1))
# normalizing_c2 = θ^2 * (StatsFuns.gamma(3)/ StatsFuns.gamma(1))
# @test d_exp.c * (d_exp.c0 + normalizing_c1 * d_exp.c1 + normalizing_c2 * d_exp.c2) == 1.0

# @test d_exp.c ≈ inv(1+ 0.5 * tr(Γ))
# @test minimum(d_exp) == 0
# @test maximum(d_exp) == Inf
# @test cdf(d_exp, 0) == 0
# @test cdf(d_exp, Inf) == 1

# ###
# Random.seed!(1234)
# nsample = 10_000 #_000
# @info "sample $nsample points for the $dist distribution"
# s = Vector{Float64}(undef, nsample)
# rand!(d_exp, s) # compile
# @time rand!(d_exp, s)
# println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_exp))")
# println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_exp))")
# end
# # n = 10,0000 ::: time = 0.012240 seconds
# # sample mean = 26.110450150884134; theoretical mean = 26.326541243044463
# # sample var = 856.9798254204945; theoretical var = 866.0981753436141

####
#########  Beta test using bisection ############
@testset "Beta(α = 2.0, β = 3.0) * (c0 + c1 * x + c2 * x^2);" begin
Random.seed!(123)
α, β = (2.0, 3.0)
dist = Beta(α, β)

# test if the proper c0, c1, c2 constants are stored.
γ = 1.0
d_beta = pdf_constants(γ, dist)
@test d_beta.c ≈ inv(1 + 0.5 * γ)
@test d_beta.c0 == 1 + (mean(d_beta.d)^2 * inv(var(d_beta.d)) * 0.5 * γ)
@test d_beta.c1 == 0.5 * γ * (-2 * mean(d_beta.d) * inv(var(d_beta.d)))
@test d_beta.c2 == 0.5 * γ * (inv(var(d_beta.d)))

# test if the CDF is a proper density (integrates to 1)
α, β = Distributions.params(d_beta.d)
normalizing_c1 = inv(StatsFuns.gamma(α) * StatsFuns.gamma(α + β + 1)) * (StatsFuns.gamma(α + β) * StatsFuns.gamma(α + 1))
normalizing_c2 = inv(StatsFuns.gamma(α) * StatsFuns.gamma(α + β + 2)) * (StatsFuns.gamma(α + β) * StatsFuns.gamma(α + 2))
@test d_beta.c * (d_beta.c0 + normalizing_c1 * d_beta.c1 + normalizing_c2 * d_beta.c2) ≈  1.0
@test minimum(d_beta) == 0
@test maximum(d_beta) ≈ 1
@test cdf(d_beta, 0) == 0
@test cdf(d_beta, 1) ≈ 1

###
# Random.seed!(1234)
# nsample = 10_000 #
# @info "sample $nsample points for the $dist distribution using the Bisection method."
# s = Vector{Float64}(undef, nsample)
# rand!(d_beta, s) # compile
# @time rand!(d_beta, s)
# println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_beta))")
# println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_beta))")
end

#########  Beta test ############
# @testset "Beta(α = 0.5, β = 3.0) * (c0 + c1 * x + c2 * x^2);" begin
# Random.seed!(123)
# α, β = (0.5, 3.0)
# dist = Beta(α, β)
#
# # test if the proper c0, c1, c2 constants are stored.
# γ = 1.0
# d_beta = pdf_constants(γ, dist)
# @test d_beta.c ≈ inv(1 + 0.5 * γ)
# @test d_beta.c0 == 1 + (mean(d_beta.d)^2 * inv(var(d_beta.d)) * 0.5 * γ)
# @test d_beta.c1 == 0.5 * γ * (-2 * mean(d_beta.d) * inv(var(d_beta.d)))
# @test d_beta.c2 == 0.5 * γ * (inv(var(d_beta.d)))
#
# # test if the CDF is a proper density (integrates to 1)
# α, β = Distributions.params(d_beta.d)
# normalizing_c1 = inv(StatsFuns.gamma(α) * StatsFuns.gamma(α + β + 1)) * (StatsFuns.gamma(α + β) * StatsFuns.gamma(α + 1))
# normalizing_c2 = inv(StatsFuns.gamma(α) * StatsFuns.gamma(α + β + 2)) * (StatsFuns.gamma(α + β) * StatsFuns.gamma(α + 2))
# @test d_beta.c * (d_beta.c0 + normalizing_c1 * d_beta.c1 + normalizing_c2 * d_beta.c2) == 1.0
# @test minimum(d_beta) == 0
# @test maximum(d_beta) ≈ 1
# @test cdf(d_beta, 0) == 0
# @test cdf(d_beta, 1) ≈ 1
#
# ###
# Random.seed!(1234)
# nsample = 10_000 #
# @info "sample $nsample points for the $dist distribution using the Bisection method."
# s = Vector{Float64}(undef, nsample)
# rand!(d_beta, s) # compile
# @time rand!(d_beta, s)
# println("sample mean = $(Statistics.mean(s)); theoretical mean = $(QuasiCopula.mean(d_beta))")
# println("sample var = $(Statistics.var(s)); theoretical var = $(QuasiCopula.var(d_beta))")
end


text:
"""
"""
function filename_extension(filename::AbstractString)
    result = lowercase(strip(splitext(strip(filename))[2]))
    return result
end