text:
using Images, FileIO
using FixedPointNumbers

function initial_means(img_data, k_num)
    return rand(img_data, k_num)
end

function update_groups!(k_groups, img_data, means)
    for iter in eachindex(img_data)
        min_idx = 1
        min_val = 1000.0 # Marking this as float made a big difference as well
        for i = 1:16
            dist = norm(means[i] - img_data[iter])
            if dist < min_val
                min_val = dist
                min_idx = i
            end
        end
        k_groups[iter] = min_idx
        # This was the previous code that was not optimal and had a lot of gc
        # Got 2x speedup by refactoring to above
        # local distance = norm.(means .- img_data[iter])
        # _, idx = findmin(distance)
        # k_groups[iter] = idx
    end
end

function update_means!(means, img_data, k_groups, k_num::Int64)
    # means = Array(RGB{Float64}, k_num)
    for i = 1:k_num
        mask = k_groups .== i
        pixels = img_data[mask]
        means[i] = mean(pixels)
    end
end

function k_means(img_data, k_num)
    means = initial_means(img_data, k_num)
    k_groups = Array(Int64, (size(img_data, 1), size(img_data,2)))
    update_groups!(k_groups, img_data, means)
    loop_num = 0
    while loop_num < 5
        update_means!(means, img_data, k_groups, k_num)
        update_groups!(k_groups,img_data, means)
        loop_num += 1
        print("Iteration $loop_num \n")
    end
    return means
end

function reduce_image!(img2, means, k_groups, k_num)
    for i = 1:k_num
        mask = k_groups .== i
        img2[mask] = means[i]
    end
end
## Main script here....
srand(1234)  # random seed
const k= 16  # How many clusters in k-means

# Load Images, Small Mandrill and large Mandrill
img1_ = load("mandrill-small.tiff")
img2_ = load("mandrill-large.tiff")
# Convert the images to float values and assign with const
const img = float64.(img1_)  # convert to float 64
const img2 = float64.(img2_)

# Perform K Means
means = k_means(img, k)

# Use the Means to compress the larger madrill picture
groups = Array(Int64, (size(img2, 1), size(img2,2)))
update_groups!(groups, img2, means)
reduce_image!(img2, means, groups, k)
save("mandrill-large-reduced-julia.tiff", img2)


text:
module XGBoost

using Printf
using Random: randperm, seed!
using SparseArrays: SparseMatrixCSC, nnz
using Statistics: mean, std

export DMatrix, Booster
export xgboost, predict, save, nfold_cv, slice, get_info, set_info, dump_model, importance

include("xgboost_lib.jl")

end # module XGBoost

ext: jl


text:
#=
VecE2: a 2d euclidean vector
=#

struct VecE2{R<:Real} <: VecE{2,R}
    x::R
    y::R
end
VecE2() = VecE2(0.0,0.0)
VecE2(x::Integer, y::Integer) = VecE2(float(x), float(y))
VecE2(t::Tuple) = VecE2(promote(t...)...)

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

Base.show(io::IO, a::VecE2) = @printf(io, "VecE2(%.3f, %.3f)", a.x, a.y)

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

dist(a::VecE2, b::VecE2) = hypot(a.x-b.x, a.y-b.y)
function dist2(a::VecE2, b::VecE2)
    Δx = a.x-b.x
    Δy = a.y-b.y
    Δx*Δx + Δy*Δy
end

proj(a::VecE2, b::VecE2, ::Type{Float64}) = (a.x*b.x + a.y*b.y) / hypot(b.x, b.y) # dot(a,b) / |b|
function proj(a::VecE2, b::VecE2, ::Type{VecE2})
    # dot(a,b) / dot(b,b) ⋅ b
    s = (a.x*b.x + a.y*b.y) / (b.x*b.x + b.y*b.y)
    VecE2(s*b.x, s*b.y)
end

lerp(a::VecE2, b::VecE2, t::Real) = VecE2(a.x + (b.x-a.x)*t, a.y + (b.y-a.y)*t)

Base.rot180(a::VecE2) = VecE2(-a.x, -a.y)
Base.rotl90(a::VecE2) = VecE2(-a.y,  a.x)
Base.rotr90(a::VecE2) = VecE2( a.y, -a.x)
function rot(a::VecE2, θ::Float64)

    #=
    Rotates counter-clockwise about origin
    =#

    c = cos(θ)
    s = sin(θ)

    VecE2(a.x*c - a.y*s, a.x*s+a.y*c)
end


text:
function Plus(a, b)
    return a + b
end

function factorial(n)
    if n == 0 or n == 1
        return 1
    end
    return n * factorial(n - 1)
end

begin
    print(
        Plus(1, 1)
    )
end

max_stars_repo_name: groundf/course-julia

max_issues_repo_name: groundf/course-julia

max_forks_repo_name: groundf/course-julia


text:
using WordCloud
words = (1:200) .% 10 .|> string
weights = (1:200) .% 11 .+ 1
#md# ### SVG
#md# You already have or want to manually generate a SVG mask with outline, you should set a proper transparent region
svgmask = shape(squircle, 300, 200, outline=3, linecolor="navy", color="AliceBlue")
wc1 = wordcloud(
    words, weights,
    mask=svgmask,
    transparent=c -> c != WordCloud.torgba("AliceBlue"), # the outline should be regarded as transparent too
) |> generate!
#md# Or, if you set the `maskcolor` in `wordcloud`, the transparent will be automatically set correctly.
#md# ```julia
#md# wc1 = wordcloud(
#md#     words, weights,
#md#     maskshape = squircle, rt=0.5,
#md#     masksize = (300, 200),
#md#     maskcolor = "AliceBlue",
#md#     outline = 3, linecolor = "navy"
#md# ) |> generate!
#md# ```
paint(wc1, "outline.svg")
println("results are saved to outline.svg")
#md# ### Bitmap
#md# If you already have a bitmap mask without outline, you can outline it before painting
bitmapmask = WordCloud.svg2bitmap(shape(squircle, 300, 200, color="AliceBlue", backgroundsize=(312, 212)))
wc2 = wordcloud(
    words, weights,
    mask=bitmapmask,
) |> generate!
paint(wc2, "outline.png", background=outline(bitmapmask, color="navy", linewidth=3, smoothness=0.8))
println("results are saved to outline.png")
#md#
wc1, wc2
#eval# runexample(:outline)
#md# ![](outline.svg)
#md# ![](outline.png)

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
using Compat
import Base.Forward

using SortingAlgorithms, Base.Order, Compat, IndexedTables, DataFrames
import Base: ht_keyindex, rehash!, _setindex!
import SortingAlgorithms: uint_mapping, RADIX_SIZE, RADIX_MASK
import DataFrames: DataFrame, AbstractDataFrame
import IndexedTables: IndexedTable, column
import PooledArrays.PooledArray
import CategoricalArrays.CategoricalArray
import SplitApplyCombine.groupreduce

function sumby_radixsort_extra{T, S<:Number}(by::AbstractVector{T},  val::AbstractVector{S})::Dict{T,S}
  by_sim = similar(by)
  val1=similar(val)
  o = Forward
  lo = 1
  hi = length(by)

  if hi == 0;  return Dict{T,S}();
  elseif hi == 1;  return Dict{T,S}(by[1] => val[1]);  end

  # Make sure we're sorting a bits type
  #TT = Base.Order.ordtype(o, by)
  if !isbits(T)
      error("Radix sort only sorts bits types (got $T)")
  end

  # Init
  iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
  bin = zeros(UInt32, 2^RADIX_SIZE, iters)
  if lo > 1;  bin[1,:] = lo-1;  end

  # Histogram for each element, radix
  for i = lo:hi
    v = uint_mapping(o, by[i])
    for j = 1:iters
        idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
        @inbounds bin[idx,j] += 1
    end
  end

  # Sort!
  swaps = 0
  len = hi-lo+1
  for j = 1:iters
      tic()
      # Unroll first data iteration, check for degenerate case
      v = uint_mapping(o, by[hi])
      idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

      # are all values the same at this radix?
      if bin[idx,j] == len;  continue;  end

      cbin = cumsum(bin[:,j])
      ci = cbin[idx]
      by_sim[ci] = by[hi]
      val1[ci] = val[hi]

      cbin[idx] -= 1

      # Finish the loop...
      @inbounds for i in hi-1:-1:lo
          v = uint_mapping(o, by[i])
          idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
          ci = cbin[idx]
          by_sim[ci] = by[i]
          val1[ci] = val[i]
          cbin[idx] -= 1
      end
      by,by_sim = by_sim,by
      val,val1 = val1,val
      swaps += 1
      end

      @inbounds if isodd(swaps)
      by,by_sim = by_sim,by
      val,val1 = val1,val
      for i = lo:hi
          by[i] = by_sim[i]
          val[i] = val1[i]
      end
      toc()
  end

  sumby_contiguous(by, val)
end


text:
#####
##### `norm`
#####

function frule((_, Δx), ::typeof(norm), x)
    y = norm(x)
    return y, _norm2_forward(x, Δx, norm(x))
end
function frule((_, Δx), ::typeof(norm), x::Number, p::Real)
    y = norm(x, p)
    ∂y = if iszero(Δx) || iszero(p)
        zero(real(x)) * zero(real(Δx))
    else
        signx = x isa Real ? sign(x) : x * pinv(y)
        _realconjtimes(signx, Δx)
    end
    return y, ∂y
end

function rrule(
    ::typeof(norm),
    x::Union{StridedArray, LinearAlgebra.AbstractTriangular, Diagonal},
    p::Real,
)
    y = LinearAlgebra.norm(x, p)
    function norm_pullback(Δy)
        ∂x = Thunk() do
            return if isempty(x) || p == 0
                zero.(x) .* (zero(y) * zero(real(Δy)))
            elseif p == 2
                _norm2_back(x, y, Δy)
            elseif p == 1
                _norm1_back(x, y, Δy)
            elseif p == Inf
                _normInf_back(x, y, Δy)
            elseif p == -Inf
                _normInf_back(x, y, Δy)
            else
                _normp_back_x(x, p, y, Δy)
            end
        end
        ∂p = @thunk _normp_back_p(x, p, y, Δy)
        return (NO_FIELDS, ∂x, ∂p)
    end
    norm_pullback(::Zero) = (NO_FIELDS, Zero(), Zero())
    return y, norm_pullback
end
function rrule(
    ::typeof(norm),
    x::Union{StridedArray, LinearAlgebra.AbstractTriangular, Diagonal},
)
    y = LinearAlgebra.norm(x)
    function norm_pullback(Δy)
        ∂x = if isempty(x)
            zero.(x) .* (zero(y) * zero(real(Δy)))
        else
            _norm2_back(x, y, Δy)
        end
        return (NO_FIELDS, ∂x)
    end
    norm_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, norm_pullback
end
function rrule(
    ::typeof(norm),
    x::Union{LinearAlgebra.TransposeAbsVec, LinearAlgebra.AdjointAbsVec},
    p::Real,
)
    y, inner_pullback = rrule(norm, parent(x), p)
    function norm_pullback(Δy)
        (∂self, ∂x′, ∂p) = inner_pullback(Δy)
        fdual = x isa Transpose ? transpose : adjoint
        ∂x = @thunk fdual(unthunk(∂x′))
        return (∂self, ∂x, ∂p)
    end
    return y, norm_pullback
end
function rrule(::typeof(norm), x::Number, p::Real)
    y = norm(x, p)
    function norm_pullback(Δy)
        ∂x = if iszero(Δy) || iszero(p)
            zero(x) * zero(real(Δy))
        else
            signx = x isa Real ? sign(x) : x * pinv(y)
            signx * real(Δy)
        end
        return (NO_FIELDS, ∂x, Zero())
    end
    norm_pullback(::Zero) = (NO_FIELDS, Zero(), Zero())
    return y, norm_pullback
end

#####
##### `normp`
#####

function rrule(
    ::typeof(LinearAlgebra.normp),
    x::Union{StridedArray, LinearAlgebra.AbstractTriangular, Diagonal},
    p,
)
    y = LinearAlgebra.normp(x, p)
    function normp_pullback(Δy)
        ∂x = @thunk _normp_back_x(x, p, y, Δy)
        ∂p = @thunk _normp_back_p(x, p, y, Δy)
        return (NO_FIELDS, ∂x, ∂p)
    end
    normp_pullback(::Zero) = (NO_FIELDS, Zero(), Zero())
    return y, normp_pullback
end

function _normp_back_x(x, p, y, Δy)
    c = real(Δy) / y
    ∂x = similar(x)
    broadcast!(∂x, x) do xi
        a = norm(xi)
        ∂xi = xi * ((a / y)^(p - 2) * c)
        return ifelse(isfinite(∂xi), ∂xi, zero(∂xi))
    end
    return ∂x
end

function _normp_back_p(x, p, y, Δy)
    y > 0 && isfinite(y) && !iszero(p) || return zero(real(Δy)) * zero(y) / one(p)
    s = sum(x) do xi
        a = norm(xi)
        c = (a / y)^(p - 1) * a * log(a)
        return ifelse(isfinite(c), c, zero(c))
    end
    ∂p = real(Δy) * (s - y * log(y)) / p
    return ∂p
end

#####
##### `normMinusInf`/`normInf`
#####

function rrule(
    ::typeof(LinearAlgebra.normMinusInf),
    x::Union{StridedArray, LinearAlgebra.AbstractTriangular, Diagonal},
)
    y = LinearAlgebra.normMinusInf(x)
    normMinusInf_pullback(Δy) = (NO_FIELDS, _normInf_back(x, y, Δy))
    normMinusInf_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, normMinusInf_pullback
end

function rrule(
    ::typeof(LinearAlgebra.normInf),
    x::Union{StridedArray,LinearAlgebra.AbstractTriangular,Diagonal},
)
    y = LinearAlgebra.normInf(x)
    normInf_pullback(Δy) = (NO_FIELDS, _normInf_back(x, y, Δy))
    normInf_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, normInf_pullback
end

function _normInf_back(x, y, Δy)
    Δu = real(Δy)
    T = typeof(zero(float(eltype(x))) * zero(Δu))
    ∂x = fill!(similar(x, T), 0)
    # if multiple `xi`s have the exact same norm, then they must have been identically
    # produced, e.g. with `fill`. So we set only one to be non-zero.
    # we choose last index to match the `frule`.
    yind = findlast(xi -> norm(xi) == y, x)
    yind === nothing && throw(ArgumentError("y is not the correct norm of x"))
    @inbounds ∂x[yind] = sign(x[yind]) * Δu
    return ∂x
end

#####
##### `norm1`
#####

function rrule(
    ::typeof(LinearAlgebra.norm1),
    x::Union{StridedArray,LinearAlgebra.AbstractTriangular,Diagonal},
)
    y = LinearAlgebra.norm1(x)
    norm1_pullback(Δy) = (NO_FIELDS, _norm1_back(x, y, Δy))
    norm1_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, norm1_pullback
end

function _norm1_back(x, y, Δy)
    ∂x = similar(x)
    ∂x .= sign.(x) .* real(Δy)
    return ∂x
end

#####
##### `norm2`
#####

function frule((_, Δx), ::typeof(LinearAlgebra.norm2), x)
    y = LinearAlgebra.norm2(x)
    return y, _norm2_forward(x, Δx, y)
end

function rrule(
    ::typeof(LinearAlgebra.norm2),
    x::Union{StridedArray,LinearAlgebra.AbstractTriangular,Diagonal},
)
    y = LinearAlgebra.norm2(x)
    norm2_pullback(Δy) = (NO_FIELDS, _norm2_back(x, y, Δy))
    norm2_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, norm2_pullback
end

function _norm2_forward(x, Δx, y)
    ∂y = real(dot(x, Δx)) * pinv(y)
    return ∂y
end
function _norm2_back(x, y, Δy)
    ∂x = similar(x)
    ∂x .= x .* (real(Δy) * pinv(y))
    return ∂x
end

#####
##### `normalize`
#####

function rrule(::typeof(normalize), x::AbstractVector, p::Real)
    nrm, inner_pullback = rrule(norm, x, p)
    Ty = typeof(first(x) / nrm)
    y = copyto!(similar(x, Ty), x)
    LinearAlgebra.__normalize!(y, nrm)
    function normalize_pullback(Δy)
        invnrm = pinv(nrm)
        ∂nrm = -dot(y, Δy) * invnrm
        (_, ∂xnorm, ∂p) = inner_pullback(∂nrm)
        ∂x = @thunk unthunk(∂xnorm) .+ Δy .* invnrm
        return (NO_FIELDS, ∂x, ∂p)
    end
    normalize_pullback(::Zero) = (NO_FIELDS, Zero(), Zero())
    return y, normalize_pullback
end
function rrule(::typeof(normalize), x::AbstractVector)
    nrm = LinearAlgebra.norm2(x)
    Ty = typeof(first(x) / nrm)
    y = copyto!(similar(x, Ty), x)
    LinearAlgebra.__normalize!(y, nrm)
    function normalize_pullback(Δy)
        ∂x = (Δy .- real(dot(y, Δy)) .* y) .* pinv(nrm)
        return (NO_FIELDS, ∂x)
    end
    normalize_pullback(::Zero) = (NO_FIELDS, Zero())
    return y, normalize_pullback
end


text:
include("utils.jl")

@testset "fixed Rank" begin
    M = FixedRankMatrices(3, 2, 2)
    M2 = FixedRankMatrices(3, 2, 1)
    Mc = FixedRankMatrices(3, 2, 2, ℂ)
    x = SVDMPoint([1.0 0.0; 0.0 1.0; 0.0 0.0])
    x2 = SVDMPoint([1.0 0.0; 0.0 1.0; 0.0 0.0], 1)
    v = UMVTVector([0.0 0.0; 0.0 0.0; 1.0 1.0], [1.0 0.0; 0.0 1.0], zeros(2, 2))
    @test repr(M) == "FixedRankMatrices(3, 2, 2, ℝ)"
    @test repr(Mc) == "FixedRankMatrices(3, 2, 2, ℂ)"
    @test sprint(show, "text/plain", x) == """
    $(sprint(show, SVDMPoint{Matrix{Float64}, Vector{Float64}, Matrix{Float64}}))
    U factor:
     3×2 $(sprint(show, Matrix{Float64})):
      1.0  0.0
      0.0  1.0
      0.0  0.0
    singular values:
     2-element $(sprint(show, Vector{Float64})):
      1.0
      1.0
    Vt factor:
     2×2 $(sprint(show, Matrix{Float64})):
      1.0  0.0
      0.0  1.0"""
    @test sprint(show, "text/plain", v) == """
    $(sprint(show, UMVTVector{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}))
    U factor:
     3×2 $(sprint(show, Matrix{Float64})):
      0.0  0.0
      0.0  0.0
      1.0  1.0
    M factor:
     2×2 $(sprint(show, Matrix{Float64})):
      1.0  0.0
      0.0  1.0
    Vt factor:
     2×2 $(sprint(show, Matrix{Float64})):
      0.0  0.0
      0.0  0.0"""

    @test inner(M, x, v, v) == norm(M, x, v)^2
    @test x == SVDMPoint(x.U, x.S, x.Vt)
    @test v == UMVTVector(v.U, v.M, v.Vt)
    @testset "Fixed Rank Matrices – Basics" begin
        @test representation_size(M) == (3, 2)
        @test representation_size(Mc) == (3, 2)
        @test manifold_dimension(M) == 6
        @test manifold_dimension(Mc) == 12
        @test !is_manifold_point(M, SVDMPoint([1.0 0.0; 0.0 0.0], 2))
        @test_throws DomainError is_manifold_point(
            M,
            SVDMPoint([1.0 0.0; 0.0 0.0], 2),
            true,
        )
        @test is_manifold_point(M2, x2)

        @test !is_tangent_vector(
            M,
            SVDMPoint([1.0 0.0; 0.0 1.0; 0.0 0.0]),
            UMVTVector(zeros(2, 1), zeros(1, 2), zeros(2, 2)),
        )
        @test !is_tangent_vector(M, SVDMPoint([1.0 0.0; 0.0 0.0], 2), v)
        @test_throws DomainError is_tangent_vector(
            M,
            SVDMPoint([1.0 0.0; 0.0 0.0], 2),
            v,
            true,
        )
        @test !is_tangent_vector(M, x, UMVTVector(x.U, v.M, x.Vt, 2))
        @test_throws DomainError is_tangent_vector(
            M,
            x,
            UMVTVector(x.U, v.M, x.Vt, 2),
            true,
        )
        @test !is_tangent_vector(M, x, UMVTVector(v.U, v.M, x.Vt, 2))
        @test_throws DomainError is_tangent_vector(
            M,
            x,
            UMVTVector(v.U, v.M, x.Vt, 2),
            true,
        )

        @test is_manifold_point(M, x)
        @test is_tangent_vector(M, x, v)
    end
    types = [[Matrix{Float64}, Vector{Float64}, Matrix{Float64}]]
    TEST_FLOAT32 && push!(types, [Matrix{Float32}, Vector{Float32}, Matrix{Float32}])

    for T in types
        @testset "Type $T" begin
            y = retract(M, x, v, PolarRetraction())
            z = SVDMPoint([1/sqrt(2) 1/sqrt(2); 1/sqrt(2) -1/sqrt(2); 0.0 0.0])
            pts = []
            for p in [x, y, z]
                push!(pts, SVDMPoint(convert.(T, [p.U, p.S, p.Vt])...))
            end
            for p in pts
                @test is_manifold_point(M, p)
            end
            @testset "SVD MPoint Basics" begin
                s = svd(x.U * Diagonal(x.S) * x.Vt)
                x2 = SVDMPoint(s)
                x3 = SVDMPoint(s.U, s.S, s.Vt)
                @test SVDMPoint(x.U, x.S, x.Vt) == x
                @test x.S == x2.S
                @test x.U == x2.U
                @test x.Vt == x2.Vt
                @test x == x2
                @test x2.U == x3.U
                @test x2.S == x3.S
                @test x2.Vt == x3.Vt
                @test x2 == x3
                y = SVDMPoint([1.0 0.0; 0.0 0.0; 0.0 0.0], 1)
                s2 = svd([1.0 0.0; 0.0 0.0; 0.0 0.0])
                y2 = SVDMPoint(s2, 1)
                y3 = SVDMPoint(s2.U, s2.S, s2.Vt, 1)
                @test y.S == y2.S
                @test y.U == y2.U
                @test y.Vt == y2.Vt
                @test y == y2
                @test y2.U == y3.U
                @test y2.S == y3.S
                @test y2.Vt == y3.Vt
                @test y2 == y3

                @test is_manifold_point(M, x)
                xM = x.U * Diagonal(x.S) * x.Vt
                @test is_manifold_point(M, xM)
                @test !is_manifold_point(M, xM[1:2, :])
                @test_throws DomainError is_manifold_point(M, xM[1:2, :], true)
                @test_throws DomainError is_manifold_point(
                    FixedRankMatrices(3, 2, 1),
                    x,
                    true,
                )
                @test_throws DomainError is_manifold_point(
                    FixedRankMatrices(3, 2, 1),
                    xM,
                    true,
                )
                xF1 = SVDMPoint(2 * x.U, x.S, x.Vt)
                @test !is_manifold_point(M, xF1)
                @test_throws DomainError is_manifold_point(M, xF1, true)
                xF2 = SVDMPoint(x.U, x.S, 2 * x.Vt)
                @test !is_manifold_point(M, xF2)
                @test_throws DomainError is_manifold_point(M, xF2, true)
            end
            @testset "UMV TVector Basics" begin
                w = UMVTVector(v.U, 2 * v.M, v.Vt)
                @test v + w == UMVTVector(2 * v.U, 3 * v.M, 2 * v.Vt)
                @test v - w == UMVTVector(0 * v.U, -v.M, 0 * v.Vt)
                @test 2 * v == UMVTVector(2 * v.U, 2 * v.M, 2 * v.Vt)
                @test v * 2 == UMVTVector(v.U * 2, v.M * 2, v.Vt * 2)
                @test 2 \ v == UMVTVector(2 \ v.U, 2 \ v.M, 2 \ v.Vt)
                @test v / 2 == UMVTVector(v.U / 2, v.M / 2, v.Vt / 2)
                @test +v == v
                @test -v == UMVTVector(-v.U, -v.M, -v.Vt)
                w = UMVTVector(v.U, v.M, v.Vt)
                @test v == w
                w = allocate(v, number_eltype(v))
                zero_tangent_vector!(M, w, x)
                oneP = SVDMPoint(one(zeros(3, 3)), ones(2), one(zeros(2, 2)), 2)
                @test oneP == one(x)
                oneV = UMVTVector(one(zeros(3, 3)), one(zeros(2, 2)), one(zeros(2, 2)), 2)
                @test oneV == one(v)
            end
            test_manifold(
                M,
                pts,
                test_exp_log=false,
                default_inverse_retraction_method=nothing,
                test_injectivity_radius=false,
                default_retraction_method=PolarRetraction(),
                test_is_tangent=false,
                test_project_tangent=true,
                test_default_vector_transport=false,
                test_forward_diff=false,
                test_reverse_diff=false,
                test_vector_spaces=false,
                test_vee_hat=false,
                test_tangent_vector_broadcasting=false, #broadcast not so easy for 3 matrix type
                projection_atol_multiplier=15,
                retraction_methods=[PolarRetraction()],
                mid_point12=nothing,
            )
        end
    end
end


text:
export convGEMMKernel,Amv,ATmv,transposeTest,getConvGEMMKernel
using DistributedArrays
mutable struct convGEMMKernel{T} <: AbstractConvKernel{T}
    nImg    		:: Array{Int,1}
    sK      		:: Array{Int,1}
	shiftX  		:: Array{Int,1}
	shiftT  		:: Array{Int,1}
	aux_sk3 		:: Array{T, 3}
	aux_sk4 		:: Array{T, 3}
	KK 				:: Array{Array{T,2}}

	prev_nex_AMV	:: Int
	AY				:: Array{T,3}
	prev_nex_ATmv 	:: Int
	ATZ 			:: Array{T,3}
	dtheta 			:: Array{T,4}

end
function getConvGEMMKernel(TYPE::Type,nImg,sK)

	if sK[1] == 1 && sK[2] == 1
		shiftX = [0;0];
		shiftT = [0;0];
	elseif sK[1] == 3 && sK[2] == 3
		shiftX = [0;-1;0;0;1;0];
		shiftT = [1;0;0;0;0;-1];
	else
		error("Code only supports 1X1 and 3X3 convolutions");
	end

	aux_sk3 = zeros(TYPE,nImg[1],nImg[2],sK[3]);
	aux_sk4 = zeros(TYPE,nImg[1],nImg[2],sK[4]);

	KK = Array{Array{TYPE,2}}(sK[1],sK[2]);
	prev_nex_AMV = -1
	AY = zeros(TYPE,0,0,0)
	prev_nex_ATmv = -1
	ATZ = zeros(TYPE,0,0,0)
	dtheta = zeros(TYPE, sK[1], sK[2], sK[3], sK[4])

	return convGEMMKernel{TYPE}(copy(nImg),copy(sK),shiftX,shiftT,aux_sk3,aux_sk4,KK,prev_nex_AMV,AY,prev_nex_ATmv,ATZ,dtheta);
end

function Amv(this::convGEMMKernel{T},theta::Array{T},Y::Array{T}) where {T<:Number}
    ## We assume that the data Y is held in the order XYCN.

	sK = this.sK;
	nImg = this.nImg;
	nex   = div(numel(Y),prod(nImgIn(this)))
	KK = this.KK

	# compute convolution
	Y     = reshape(Y,nImg[1],nImg[2],this.sK[3],nex);

	if nex != this.prev_nex_AMV
		this.AY    = Array{T, 3}(nImg[1]*nImg[2],this.sK[4],nex);
	end

	aux   = this.aux_sk3;
    AYk   = reshape(this.aux_sk4,nImg[1]*nImg[2],sK[4]);
	### reshape the kernels for gemm!:
	K = reshape(theta, sK[1], sK[2], sK[3], sK[4])

	for k1 = 1:sK[1]
		for k2 = 1:sK[2]
			@inbounds KK[k1,k2] = K[k1,k2,:,:]';
		end
	end

	for k = 1:nex
		AYk[:] = zero(T)
		AYk = multConv2Dblock(Y,KK, AYk,aux,this.shiftX,this.shiftT,k);
		@inbounds this.AY[:,:,k] = AYk;
	end
	AY_out = reshape(this.AY,:,nex)
	this.prev_nex_AMV = nex

    return AY_out
end

function ATmv(this::convGEMMKernel{T},theta::Array{T},Zin::Array{T}) where {T<:Number}

	nImg  = this.nImg;
	sK    = this.sK;
    nex   =  div(numel(Zin),prod(nImgOut(this)));
    K     = reshape(theta, sK[1], sK[2], sK[3], sK[4]);
	Z     = reshape(Zin,nImg[1],nImg[2],sK[4],nex);
	aux   = this.aux_sk4;
	# ATZ = this.ATZ
	ATZk  = reshape(this.aux_sk3,nImg[1]*nImg[2],sK[3]);

	if nex != this.prev_nex_ATmv
		this.ATZ   = zeros(T,nImg[1]*nImg[2],sK[3],nex);
	end

	### reshape the kernels for gemm!:
	KK = this.KK

	for k1 = 1:sK[1]
		for k2 = 1:sK[2]
			@inbounds KK[k1,k2] = K[k1,k2,:,:];
		end
	end
	## flipping:
	KK = flipdim(flipdim(KK,2),1);
	for k = 1:nex
		ATZk[:] = zero(T)
		ATZk = multConv2Dblock(Z,KK, ATZk,aux,this.shiftX,this.shiftT,k);
		@inbounds this.ATZ[:,:,k] = ATZk;
	end
	ATZ_out = reshape(this.ATZ,:,nex);
	this.prev_nex_ATmv = nex

    return ATZ_out
end

function Jthetamv(this::convGEMMKernel{T},dtheta::Array{T},dummy::Array{T},Y::Array{T},temp=nothing) where {T<:Number}
    nex    =  div(numel(Y),nFeatIn(this));
    Z      = Amv(this,dtheta,Y);
    return Z
end

function JthetaTmv(this::convGEMMKernel{T}, Zin::Array{T}, dummy::Array{T}, Yin::Array{T}) where {T<:Number}
     # derivative of Z*(A(theta)*Y) w.r.t. theta
	sK = this.sK
	nImg = this.nImg
	nex   = div(numel(Yin),prod(nImgIn(this)))
    # compute convolution
	Y     = reshape(Yin, nImg[1], nImg[2], this.sK[3], nex)
	Z	  = reshape(Zin, nImg[1]*nImg[2], this.sK[4], nex)
	Zk    = reshape(this.aux_sk4, nImg[1]*nImg[2], this.sK[4]);
	aux   = this.aux_sk3;
	dtheta = zeros(T, sK[1], sK[2], sK[3], sK[4])
	### reshape the kernels for gemm!:
	KK = this.KK
	for k1 = 1:sK[1]
		for k2 = 1:sK[2]
			@inbounds KK[k1, k2] = zeros(T, sK[3], sK[4])
		end
	end
    for k = 1:nex
		getColumn!(Z, Zk, k)
		multConv2Dblock(Y, KK,  Zk, aux, this.shiftX, this.shiftT, k, doDerivative = 1)
	end
	### Assemble the kernels from gemm!:
	for k1 = 1:sK[1]
		for k2 = 1:sK[2]
			@inbounds dtheta[k1, k2, :, :] = KK[k1, k2]
		end
	end
    dtheta_out = reshape(dtheta, sK[1], sK[2], sK[3], sK[4])
    return dtheta_out
end

function getColumn!(Z::Array{T},Zk::Array{T},k::Int64) where {T<:Number}
for c=1:size(Z,2)
	for j=1:size(Z,1)
		@inbounds	Zk[j,c] = Z[j,c,k];
	end
end
end

function multConv2Dblock(x::Array{T},K::Array{Array{T,2},2}, y::Array{T}, tin::Array{T},shiftX,shiftT,imIdx;doDerivative = 0) where {T<:Number}
	## y = K*x
	## K - 3X3 array of Arrays
	## x - a vector of length |nImgag+2|*cin (zero padded)
	## y - a vector of length |nImgag|*cout
	nImg1 = size(x,1);
	nImg2 = size(x,2);
	cin = size(x,3);
	cout = size(y,2);
	OneType = one(T);
	t = reshape(tin,nImg1,nImg2,cin);
	kernelWidth = size(K,1);
	# y = reshape(y,nImg1*nImg2,cout); # it is supposed to be of this shape...
	k=1;
	jt=0;it=0;jt=0;jx=0;

	for p = 1:2:2*kernelWidth
		for q = 1:2:2*kernelWidth
			lower = nImg2+shiftT[p+1]  # Move outside of the forloop for increased speed
			upper = nImg1+shiftT[q+1]  # Move outside of the forloop for increased speed
			for cc = 1:cin
				jx = 1+shiftX[p];  # Moving these outside didn't seem to help
				jt = 1+shiftT[p];
				if jt > 1
					###################### Dirichlet #######################
					@inbounds t[:,1:(jt-1),cc] = zero(T);
					###################### Periodic #######################
					# ix = 1+shiftX[q];
					# if shiftT[q] > 0
						#@inbounds t[1,1,cc] = x[end,end,cc,imIdx];
					# end
					# for it = (1+shiftT[q]):upper
						#@inbounds t[it,1,cc] = x[ix,end,cc,imIdx];
						# ix +=1;
					# end
					# if shiftT[q+1] < 0
						#@inbounds t[end,1,cc] = x[1,end,cc,imIdx];
					# end
					###################### End Periodic #######################
				end
				while jt <= lower
					it = 1+shiftT[q];
					ix = 1+shiftX[q];
					if it > 1
						for ii = 1:(it-1)
							###################### Dirichlet #######################
							@inbounds t[ii,jt,cc] = zero(T)   #@inbounds t[1:(it-1),jt,cc] = 0.0 - faster unvectorized
							###################### Periodic #######################
							#@inbounds t[ii,jt,cc] = x[end,jx,cc,imIdx];
						end
					end
					while it <= upper
						@inbounds t[it,jt,cc] = x[ix,jx,cc,imIdx];
						it+=1;ix+=1;
					end
					if it <= nImg1
						for ii = it:nImg1
							###################### Dirichlet #######################
							@inbounds t[ii,jt,cc] = zero(T)	#@inbounds t[it:nImg1,jt,cc] = 0.0 - faster unvectorized
							###################### Periodic #######################
							# @inbounds t[ii,jt,cc] = x[1,jx,cc,imIdx];
						end
					end
					jt+=1;jx+=1;

				end
				if jt <= nImg2
					###################### Dirichlet #######################
					@inbounds t[:,jt:nImg2,cc] = zero(T);
					###################### Periodic #######################
					# if shiftT[q] > 0
						# @inbounds t[1,end,cc] = x[end,1,cc,imIdx];
					# end
					# ix = ix = 1+shiftX[q];
					# for it = (1+shiftT[q]):upper
						# @inbounds t[it,end,cc] = x[ix,1,cc,imIdx];
						# ix +=1;
					# end
					# if shiftT[q+1] < 0
						# @inbounds t[end,end,cc] = x[1,1,cc,imIdx];
					# end
					###################### End Periodic #######################
				end
			end
			if doDerivative == 0

				BLAS.gemm!('N','T',OneType,reshape(t,nImg1*nImg2,cin),K[k],OneType,y);
			else
				BLAS.gemm!('T','N',OneType,reshape(t,nImg1*nImg2,cin),y,OneType,K[k]);
			end
			k+=1;
		end
	end
	return y;
end


text:
include("dependencies_for_runtests.jl")

using Statistics
using Oceananigans.Architectures: arch_array
using Oceananigans.Fields: ReducedField, CenterField, ZFaceField, compute_at!, @compute
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Grids: halo_size

@testset "Fields computed by Reduction" begin
    @info "Testing Fields computed by reductions..."

    for arch in archs
        arch_str = string(typeof(arch))

        grid = RectilinearGrid(arch, size = (2, 2, 2),
                               x = (0, 2), y = (0, 2), z = (0, 2),
                               topology = (Periodic, Periodic, Bounded))

        Nx, Ny, Nz = size(grid)

        w = ZFaceField(grid)
        T = CenterField(grid)

        trilinear(x, y, z) = x + y + z

        set!(T, trilinear)
        set!(w, trilinear)

        @compute Txyz = Field(Average(T, dims=(1, 2, 3)))

        # Note: halo regions must be *filled* prior to computing an average
        # if the average within halo regions is to be correct.
        fill_halo_regions!(T, arch)
        @compute Txy = Field(Average(T, dims=(1, 2)))

        fill_halo_regions!(T, arch)
        @compute Tx = Field(Average(T, dims=1))

        @compute wxyz = Field(Average(w, dims=(1, 2, 3)))
        @compute wxy = Field(Average(w, dims=(1, 2)))
        @compute wx = Field(Average(w, dims=1))

        @testset "Averaged fields [$arch_str]" begin
            @info "  Testing averaged Fields [$arch_str]"

            @test Txyz[1, 1, 1] ≈ 3
            @test Array(interior(Txy))[1, 1, :] ≈ [2.5, 3.5]
            @test Array(interior(Tx))[1, :, :] ≈ [[2, 3] [3, 4]]
            @test wxyz[1, 1, 1] ≈ 3
            @test Array(interior(wxy))[1, 1, :] ≈ [2, 3, 4]
            @test Array(interior(wx))[1, :, :] ≈ [[1.5, 2.5] [2.5, 3.5] [3.5, 4.5]]

            # Test whether a race condition gets hit for averages over large fields
            big_grid = RectilinearGrid(arch,
                                       topology = (Periodic, Periodic, Bounded),
                                       size = (256, 256, 128),
                                       x = (0, 2),
                                       y = (0, 2),
                                       z = (0, 2))

            c = CenterField(big_grid)
            c .= 1

            C = Field(Average(c, dims=(1, 2)))

            # Test that the mean consistently returns 1 at every z for many evaluations
            # Warm up
            for i = 1:10
                sum!(C, C.operand.operand)
            end

            results = []
            for i = 1:10
                sum!(C, C.operand.operand)
                push!(results, all(interior(C) .== 1))
            end

            @test mean(results) == 1.0
        end

        @testset "Allocating reductions [$arch_str]" begin
            @info "  Testing allocating reductions"

            # Mean
            @test Txyz[1, 1, 1] == mean(T)
            @test interior(Txy) == interior(mean(T, dims=(1, 2)))
            @test interior(Tx) == interior(mean(T, dims=1))

            @test wxyz[1, 1, 1] == mean(w)
            @test interior(wxy) == interior(mean(w, dims=(1, 2)))
            @test interior(wx) == interior(mean(w, dims=1))

            # Maximum and minimum
            @test maximum(T) == maximum(interior(T))
            @test minimum(T) == minimum(interior(T))

            for dims in (1, (1, 2))
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
                @test interior(minimum(T; dims)) == minimum(interior(T); dims)
            end
        end

        @testset "Conditional computation of averaged Fields [$(typeof(arch))]" begin
            @info "  Testing conditional computation of averaged Fields [$(typeof(arch))]"
            for FT in float_types
                grid = RectilinearGrid(arch, FT, size=(2, 2, 2), extent=(1, 1, 1))
                c = CenterField(grid)

                for dims in (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))
                    C = Field(Average(c, dims=dims))

                    @test !isnothing(C.status)

                    # Test conditional computation
                    set!(c, 1)
                    compute_at!(C, FT(1)) # will compute
                    @test all(interior(C) .== 1)
                    @test C.status.time == FT(1)

                    set!(c, 2)
                    compute_at!(C, FT(1)) # will not compute because status == 1
                    @test C.status.time == FT(1)
                    @test all(interior(C) .== 1)

                    compute_at!(C, FT(2))
                    @test C.status.time == FT(2)
                    @test all(interior(C) .== 2)
                end
            end
        end
    end
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: DoublePipeIncr.mso								$
#*------------------------------------------------------------------
type Properties_Wall
	Properties_Wall()=begin
		N=outers.N
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			fill(viscosity(Dict{Symbol,Any}(
				:Brief=>"Stream Incremental Viscosity",
				:Default=>1,
				:Lower=>1e-5,
				:Upper=>1e5,
				:Symbol=>"\\mu"
			)),(N)),
			fill(temperature(Dict{Symbol,Any}(
				:Brief=>"Incremental Wall Temperature",
				:Lower=>50
			)),(N)),
			[:N,],
			[:Mu,:Twall,]
		)
	end
	N::DanaInteger
	Mu::Array{viscosity}
	Twall::Array{temperature}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Properties_Wall
function atributes(in::Properties_Wall,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Incremental Physical properties of the streams at wall temperature."
	fields[:Info]="to be documented."
	drive!(fields,_)
	return fields
end
Properties_Wall(_::Dict{Symbol,Any})=begin
	newModel=Properties_Wall()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Catlab.Theories: adom, attr, attrtype, attr, adom, acodom
using Catlab.Present
using Catlab.CategoricalAlgebra.CSetDataStructures: struct_acset

const VPSI = Vector{Pair{Symbol, Int}}

"""
To compute automorphisms of Attributed CSets, we create a pseudo CSet which has
additional components for each data type.

This is inefficient for attributes which have a total order on them
(e.g. integers/strings) since we solve for a canonical permutation of the
attributes. Future work could address this by initializing the coloring with
the 'correct' canonical order.
"""
function pseudo_cset(g::StructACSet{S}
                    )::Tuple{StructACSet, Dict{Symbol,Vector{Any}}} where {S}
  tabs, arrs = collect(ob(S)), collect(hom(S))
  src, tgt = dom(S), codom(S)
  dtabs, darrs = collect(attrtype(S)), collect(attr(S))
  dsrc, dtgt = adom(S), acodom(S)

  # Create copy of schema (+ an extra component for each datatype)
  pres = Presentation(FreeSchema)
  xobs = [Ob(FreeSchema,t) for t in vcat(tabs,dtabs)]
  xobsdic = Dict([t=>Ob(FreeSchema,t) for t in vcat(tabs,dtabs)])
  n = length(tabs)
  for x in xobs add_generator!(pres, x) end

  for (arr, s, t) in zip(map(collect, [arrs, src, tgt])...)
    add_generator!(pres, Hom(arr, xobsdic[s], xobsdic[t]))
  end

  # Add an arrow for each attr and store all possible Data values
  attrvals = Dict([t=>Set() for t in dtabs])
  for (arr, s, t) in zip(darrs, dsrc, dtgt)
    add_generator!(pres, Hom(arr, xobsdic[s], xobsdic[t]))
    union!(attrvals[t], Set(g[arr]))
  end

  # Use Julia ordering to give each value an index
  attrvals = Dict([k=>sort(collect(x)) for (k,x) in collect(attrvals)])

  # Create and populate pseudo-cset
  name = Symbol("Pseudo_$(typeof(g).name.name)")
  expr = struct_acset(name, StructACSet, pres, index=vcat(arrs,darrs))
  eval(expr)
  csettype = eval(name)
  res = Base.invokelatest(eval(name))

  # Copy the original non-attr data
  for t in tabs
    add_parts!(res, t, nparts(g,t))
  end
  for a in arrs
    set_subpart!(res, a, g[a])
  end

  # initialize pseudo data components
  for t in dtabs
    add_parts!(res, t, length(attrvals[t]))
  end

  # Replace data value with an index for each attribute
  for (a,t) in zip(darrs, dtgt)
    fks = [findfirst(==(v), attrvals[t]) for v in g[a]]
    set_subpart!(res, a, fks)
  end

  return res, attrvals
end

pluspair(x::Pair{Int, Int}, y::Pair{Int,Int})::Pair{Int,Int} = x[1]+y[1] => x[2]+y[2]

"""
Heuristic for total ordering of objects of a CSet. We want things that will
be easily distinguished to be first. Things are very easily distinguished by
in-arrows, and somewhat distinguished by out arrows. So each object gets a
score based on this criteria. We then iterate until convergence.

We use this total ordering to order arrows by their target's score.

As an example, consider the triangle:
     f
  A --> B
  g ↘  ↙ h
     C
The stable order of the objects is:
  [A (low, high), B (mid mid), C (high, low)]
This induces an order on the arrows:
  [f (low+mid, high+mid), g (low+high, high+low), h (mid+high, mid+low)]

TO DO test this heuristic vs random heuristics to see that it's actually
effective.
"""
function order_syms(::StructACSet{S})::Vector{Symbol} where {S}
  os, arrs, srcs, tgts = ob(S), hom(S), dom(S), codom(S)
  function score_obj(obj::Symbol, scores::Dict{Symbol, Pair{Int,Int}})::Pair{Int,Int}
    arr_in = sum([scores[s][1] for (s, t) in zip(srcs, tgts) if t == obj])
    arr_ot = sum([scores[t][2] for (s, t) in zip(srcs, tgts) if s == obj])
    return arr_in => arr_ot
  end
  getorder(scores)::Vector{Symbol} = map(last, sort([(b,a) for (a,b)
                                                     in collect(scores)]))
  scores = Dict(zip(os, [(1=>1) for _ in os]))
  oldorder = []
  neworder = getorder(scores)
  while oldorder != neworder
    oldorder = neworder
    scores = Dict([o=> pluspair(scores[o],score_obj(o, scores)) for o in os])
    neworder = getorder(scores)
  end
  ordarrs = [(h, pluspair(scores[s], scores[t]))
             for (h, s, t) in zip(arrs, srcs, tgts)]
  return reverse(getorder(ordarrs))
end

"""
Inverse of pseudo_cset. Requires mapping (generated by `pseudo_cset`) of indices
for each Data to the actual data values.
"""
function pseudo_cset_inv(g::StructACSet,
                         orig::StructACSet{S},
                         attrvals::Dict{Symbol,Vector{Any}}
                        )::StructACSet{S} where {S}
  orig = deepcopy(orig)
  arrs = hom(S)
  darrs, dtabs = attr(S), acodom(S)
  for arr in arrs
    set_subpart!(orig, arr, g[arr])
  end
  for (darr,tgt) in zip(darrs, dtabs)
    set_subpart!(orig, darr, attrvals[tgt][g[darr]])
  end
  return orig
end

"""Lexicographic minimum of all automorphisms"""
function canonical_iso(g::StructCSet{S})::StructCSet{S} where {S}
  os = order_syms(g)
  opt = sort(collect(autos(g)[1]), by=(γ->order_perms(g, γ, os)))
  ord(x) = [x[s] for s in os]
  isos = sort([apply_automorphism(g, Dict(a)) for a in autos(g)[1]], by=ord)
  return isempty(isos) ? g : isos[1]
end

function isos(g::StructACSet{S})::Vector{StructACSet{S}} where {S}
  p, avals = pseudo_cset(g)
  return [pseudo_cset_inv(apply_automorphism(p, Dict(a)), g, avals)
               for a in autos(p)[1]]
end

"""
Compute automorphisms for the pseudo-cset, but then substitute in
the actual attribute values before evaluating the lexicographic order
"""
function canonical_iso(g::StructACSet{S}; pres::Union{Nothing,Presentation})::StructACSet{S} where {S}
  if !isnothing(pres) return canonical_iso_nauty(g, pres) end
  os = order_syms(g)
  ord(x) = vcat([x[a] for a in attr(S)],[x[s] for s in os])

  is = isos(g)
  !isempty(is) || error("Empty isos?")
  return sort(is, by=ord)[1]
end

"""Hash of canonical isomorphism."""
function canonical_hash(g::StructACSet;
                        pres::Union{Nothing,Presentation}=nothing)::UInt64
  hash(string(canonical_iso(g; pres=pres)))
end

"""Find index at which two vectors diverge (used in `search_tree`)"""
function common(v1::Vector{T}, v2::Vector{T})::Int where {T}
  for (i, (x, y)) in enumerate(zip(v1, v2))
    if x != y
      return i-1
    end
  end
  return min(length(v1), length(v2))
end

"""Search tree explored by Nauty"""
mutable struct Tree
  coloring::CDict
  saturated::CDict
  indicator::UInt64
  children::Dict{Pair{Symbol, Int}, Tree}
  function Tree()
    return new(CDict(),CDict(),UInt64(0),
               Dict{Pair{Symbol, Int},Tree}())
  end
end

"""
Keep track of what happens during the automorphism search. Possible actions:
- start_iter: val is a search tree + a location in the search tree
- add_leaf
- auto_prune
- order_prune
- orbit_prune
- flag_skip
- new_child
"""
mutable struct History
  action::String
  val::Any
end

"""
Container for the optimal automorphism seen so far
Store values of indicator function for each node in search tree
en route to the leaf node
"""
mutable struct Indicator
  val::Union{Nothing,Vector{UInt64}}
  function Indicator()
    return new(nothing)
  end
end

hpush!(hist::Vector{History}, act::String, val::Any) = (isempty(hist) ?
  nothing : push!(hist, History(act, val)))

function Base.size(t::Tree)::Int
  return 1 + (isempty(t.children) ? 0 : sum(map(Base.size, values(t.children))))
end

"""Get a node via a sequence of edges from the root"""
function Base.getindex(t::Tree, pth::VPSI)::Tree
  ptr = t
  for p in pth
    ptr = ptr.children[p]
  end
  return ptr
end

function get_indicators(t::Tree, pth::VPSI)::Vector{UInt64}
  ptr, inds = t, UInt64[]
  for p in pth
    ptr = ptr.children[p]
    push!(inds, ptr.indicator)
  end
  return inds
end

"""Automorphism based pruning when we've found a new leaf node (τ @ t)"""
function compute_auto_prune(tree::Tree, t::VPSI, τ::CDict, leafnodes::Set{VPSI})::Set{VPSI}
  skip = Set{VPSI}()
  for p in filter(!=(t), leafnodes)
    π = tree[p].saturated
    i = common(p, t)
    a, b, c = abc = p[1:i], p[1:i+1], t[1:i+1] # == t[1:i]
    a_, b_, c_ = [tree[x].saturated for x in abc]
    γ = compose_perms(π, invert_perms(τ))
    if (compose_perms(γ, a_) == a_ &&
        compose_perms(γ, b_) == c_)
      # skip everything from c to a
      for i in length(c):length(t)
        push!(skip, t[1:i])
      end
      break
    end
  end
  return filter(!=(t), skip) # has something gone wrong if we need to do this?
end

"""To reduce branching factor, split on the SMALLEST nontrivial partition"""
function split_data(coloring::CDict)::Tuple{Symbol, Int, Vector{Int}}
  colors_by_size = sort(get_colors_by_size(coloring), rev=true)
  if isempty(colors_by_size)
    return :_nothing, 0, []
  end
  split_tab, split_color = colors_by_size[1][2]
  colors = coloring[split_tab]
  split_inds = findall(==(split_color), colors)
  return split_tab, split_color, split_inds
end

"""
DFS tree of colorings, with edges being choices in how to break symmetry
Goal is to acquire all leaf nodes.

Algorithm from "McKay’s Canonical Graph Labeling Algorithm" by Hartke and
Radcliffe (2009).

McKay's "Practical Graph Isomorphism" (Section 2.29: "storage of identity
nodes") warns that it's not a good idea to check for every possible automorphism
pruning (for memory and time concerns). To do: look into doing this in a more
balanced way. Profiling code will probably reveal that checking for automorphism
pruning is a bottleneck.

Inputs:
 - g: our structure that we are computing automorphisms for
 - res: all automorphisms found so far
 - split_seq: sequence of edges (our current location in the tree)
 - tree: all information known so far - this gets modified
 - leafnodes: coordinates of all automorphisms found so far
 - perm_worsts: data used for order pruning
 - skip: flagged coordinates which have been pruned
 - history: log all the transformations made to a data structure for analysis
 - auto_prune:  whether or not to use automorphisms to prune
 - orbit_prune: whether or not to use orbit information to prune
 - order_prune: whether or not to use "indicator functions" to prune
"""
function search_tree!(g::StructACSet{S},
                      init_coloring::CDict,
                      split_seq::VPSI,
                      tree::Tree,
                      leafnodes::Set{VPSI},
                      indicator::Indicator,
                      skip::Set{VPSI};
                      history::Vector{History}=History[],
                      auto_prune::Bool=true,
                      orbit_prune::Bool=true,
                      order_prune::Bool=true,
                     )::Nothing where {S}
  # Perform color saturation
  color_seq, new_ind = color_saturate(g; init_color=init_coloring, history=!isempty(history))
  curr_tree = tree[split_seq]
  curr_tree.coloring = init_coloring
  curr_tree.indicator = new_ind
  curr_tree.saturated = coloring = color_seq[end]
  curr_inds = get_indicators(tree, split_seq)
  ci_len = length(curr_inds)

  if order_prune && !isnothing(indicator.val) && (
      curr_inds < let v=indicator.val; v[1:min(length(v), ci_len)] end)
    hpush!(history, "order_prune", (curr_inds, indicator.val))
    return nothing
  end

  split_tab, split_color, split_inds = sd = split_data(coloring)
  hpush!(history, "start_iter", (split_seq, color_seq, new_ind, sd))

  # Check if we are now at a leaf node
  if isempty(split_inds)
    # Add result to list of results
    push!(leafnodes, split_seq)
    invert_perms(coloring) # fail if not a perm
    #hpush!(history, "add_leaf", split_seq)

    if isnothing(indicator.val) || indicator.val < curr_inds
      indicator.val = curr_inds
    end

    # Prune with automorphisms
    if auto_prune
      pruned = compute_auto_prune(tree, split_seq, coloring, leafnodes)
      if !isempty(pruned)
        hpush!(history, "auto_prune", pruned)
        union!(skip, pruned)
      end
    end

    return nothing
  end
  # Branch on this leaf
  for split_ind in split_inds
    if split_ind != split_inds[1]
      hpush!(history, "return", split_seq)
    end

    if split_seq ∈ skip
      hpush!(history, "flag_skip", split_ind)
    elseif  orbit_prune && orbit_check(curr_tree, split_tab, split_color)
      hpush!(history, "orbit_prune", nothing)
    else
      # Construct arguments for recursive call to child
      new_coloring = deepcopy(coloring)
      new_seq = vcat(split_seq, [split_tab => split_ind])
      new_coloring[split_tab][split_ind] = maximum(coloring[split_tab]) + 1
      curr_tree.children[split_tab => split_ind] = Tree()
      search_tree!(g, new_coloring, new_seq, tree, leafnodes, indicator, skip;
                   auto_prune=auto_prune, orbit_prune=orbit_prune,
                   order_prune=order_prune, history=history)
    end
  end

  return nothing
end

"""
Get coordinates of all nodes in a tree that have no children
"""
function get_leaves(t::Tree)::Vector{VPSI}
  if isempty(t.children)
    return [Pair{Symbol,Int}[]]
  else
    res = []
    for (kv, c) in t.children
      for pth in get_leaves(c)
        push!(res, vcat([kv],pth))
      end
    end
    return res
  end
end

"""
Look at all automorphisms between leaf nodes. If they generate a group with
just one orbit in the specific part where the current tree got differentiated,
we need not explore further (returns true)
"""
function orbit_check(t::Tree, pₛ::Symbol, pᵢ::Int; verbose::Bool=false)::Bool
  leaves, auts = filter(l->is_perms(t[l].saturated),get_leaves(t)), Perm[]
  for i in 1:length(leaves)
    icolor = t[leaves[i]].saturated[pₛ]
    for j in i+1:length(leaves)
      jcolor = t[leaves[j]].saturated[pₛ]
      push!(auts,Perm(compose_comp(icolor, Base.invperm(jcolor))))
    end
  end
  all_i = findall(==(pᵢ), t.saturated[pₛ])
  length(all_i) > 1 || error("branched on nontrivial color")
  o = Orbit(auts, all_i[1])
  if verbose
    println("pₛpᵢ=$pₛ$pᵢ all_i=$all_i \nauts=$auts\norb=$(collect(o))") end
  return isempty(setdiff(all_i, collect(o)))
end

"""Compute the automorphisms of a CSet"""
function autos(g::StructACSet;
               history::Bool=false,
               auto_prune::Bool=true,
               orbit_prune::Bool=true,
               order_prune::Bool=true)::Tuple{Set{CDict},Tree, Vector{History}}
  tree, leafnodes = Tree(), Set{VPSI}()
  hist = history ? [History("start", g)] : History[]
  search_tree!(g, nocolor(g), VPSI(), tree, leafnodes,Indicator(),Set{VPSI}();
               history=hist, auto_prune=auto_prune, orbit_prune=orbit_prune,
               order_prune=order_prune)
  return Set([tree[ln].saturated for ln in leafnodes]), tree, hist
end


text:
p = plot(df, x=:comp1, y=:comp2, color=:group, Geom.point(tag=:lda));
figure()
hfig = display(p)
hit((hfig,:lda), @guarded (figtag, i, event, dist) -> if dist < 2 showimg(imgs, i) end)
lasso_initialize(hfig, (figno,selections) -> showimgs(imgs, first(values(selections))))
nothing


text:
@doc """
	trapezquad(x::Vector, y::Vector) -> Number

Approximate the integral of `y` over the the interval `[min(x), max(x)]`, using the divisions in `x`.
"""->
function trapezquad(x::AbstractVector, y::AbstractVector)
	2 <= (Nx = length(x)) == length(y) || throw(DimensionMismatch())
	!issorted(x) && throw(AssertionError("x must be sorted"))

	integral = 0.0
	for n in 2:Nx
		@inbounds integral += (x[n]-x[n-1])*0.5*(y[n]+y[n-1])
	end

	return integral
end

@doc """
	inner(x::Vector, y::Vector, z::Vector) -> Number

Approximate the inner product between `y` and `z` over the the interval `[min(x), max(x)]`, using the divisions in `x`.
The inner product is conjugate linear in `z`.
"""->
function inner(x::AbstractVector, y::AbstractVector, z::AbstractVector)
	trapezquad( x, y.*conj(z) )
end

@doc """
	l2norm(x::Vector, y::Vector) -> Number

Approximate the L2 norm squared of `y` over the the interval `[min(x), max(x)]`, using the divisions in `x`.
"""->
function l2norm(x::AbstractVector, y::AbstractVector)
	trapezquad( x, abs2(y) )
end


text:
# Noop Codec
# ==========

"""
    Noop()

Create a noop codec.

Noop (no operation) is a codec that does nothing. The data read from or written
to the stream are kept as-is without any modification. This is often useful as a
buffered stream or an identity element of a composition of streams.

The implementations are specialized for this codec. For example, a `Noop` stream
uses only one buffer rather than a pair of buffers, which avoids copying data
between two buffers and the throughput will be larger than a naive
implementation.
"""
struct Noop <: Codec end

const NoopStream{S} = TranscodingStream{Noop,S} where S<:IO

"""
    NoopStream(stream::IO)

Create a noop stream.
"""
function NoopStream(stream::IO; kwargs...)
    return TranscodingStream(Noop(), stream; kwargs...)
end

function TranscodingStream(codec::Noop, stream::IO;
                           bufsize::Integer=DEFAULT_BUFFER_SIZE,
                           sharedbuf::Bool=(stream isa TranscodingStream))
    checkbufsize(bufsize)
    checksharedbuf(sharedbuf, stream)
    if sharedbuf
        buffer = stream.state.buffer1
    else
        buffer = Buffer(bufsize)
    end
    return TranscodingStream(codec, stream, State(buffer, buffer))
end

"""
    position(stream::NoopStream)

Get the current poition of `stream`.

Note that this method may return a wrong position when
- some data have been inserted by `TranscodingStreams.unread`, or
- the position of the wrapped stream has been changed outside of this package.
"""
function Base.position(stream::NoopStream)
    return position(stream.stream) - buffersize(stream.state.buffer1)
end

function Base.seek(stream::NoopStream, pos::Integer)
    seek(stream.stream, pos)
    initbuffer!(stream.state.buffer1)
    return
end

function Base.seekstart(stream::NoopStream)
    seekstart(stream.stream)
    initbuffer!(stream.state.buffer1)
    return
end

function Base.seekend(stream::NoopStream)
    seekend(stream.stream)
    initbuffer!(stream.state.buffer1)
    return
end

function Base.unsafe_read(stream::NoopStream, output::Ptr{UInt8}, nbytes::UInt)
    changemode!(stream, :read)
    buffer = stream.state.buffer1
    p = output
    p_end = output + nbytes
    while p < p_end && !eof(stream)
        if buffersize(buffer) > 0
            m = min(buffersize(buffer), p_end - p)
            copydata!(p, buffer, m)
        else
            # directly read data from the underlying stream
            m = p_end - p
            Base.unsafe_read(stream.stream, p, m)
        end
        p += m
    end
    if p < p_end && eof(stream)
        throw(EOFError())
    end
    return
end

function Base.unsafe_write(stream::NoopStream, input::Ptr{UInt8}, nbytes::UInt)
    changemode!(stream, :write)
    buffer = stream.state.buffer1
    if marginsize(buffer) ≥ nbytes
        copydata!(buffer, input, nbytes)
        return Int(nbytes)
    else
        flushbuffer(stream)
        # directly write data to the underlying stream
        return unsafe_write(stream.stream, input, nbytes)
    end
end

function Base.transcode(::Type{Noop}, data::ByteData)
    return transcode(Noop(), data)
end

function Base.transcode(::Noop, data::ByteData)
    # Copy data because the caller may expect the return object is not the same
    # as from the input.
    return Vector{UInt8}(data)
end

# Stats
# -----

function stats(stream::NoopStream)
    state = stream.state
    mode = state.mode
    @checkmode (:idle, :read, :write)
    buffer = state.buffer1
    @assert buffer == stream.state.buffer2
    if mode == :idle
        consumed = supplied = 0
    elseif mode == :read
        supplied = buffer.transcoded
        consumed = supplied - buffersize(buffer)
    elseif mode == :write
        supplied = buffer.transcoded + buffersize(buffer)
        consumed = buffer.transcoded
    else
        @assert false "unreachable"
    end
    return Stats(consumed, supplied, supplied, supplied)
end

# Buffering
# ---------
#
# These methods are overloaded for the `Noop` codec because it has only one
# buffer for efficiency.

function fillbuffer(stream::NoopStream; eager::Bool = false)
    changemode!(stream, :read)
    buffer = stream.state.buffer1
    @assert buffer === stream.state.buffer2
    if stream.stream isa TranscodingStream && buffer === stream.stream.state.buffer1
        # Delegate the operation when buffers are shared.
        return fillbuffer(stream.stream, eager = eager)
    end
    nfilled::Int = 0
    while ((!eager && buffersize(buffer) == 0) || (eager && makemargin!(buffer, 0, eager = true) > 0)) && !eof(stream.stream)
        makemargin!(buffer, 1)
        nfilled += readdata!(stream.stream, buffer)
    end
    buffer.transcoded += nfilled
    return nfilled
end

function flushbuffer(stream::NoopStream, all::Bool=false)
    changemode!(stream, :write)
    buffer = stream.state.buffer1
    @assert buffer === stream.state.buffer2
    nflushed::Int = 0
    if all
        while buffersize(buffer) > 0
            nflushed += writedata!(stream.stream, buffer)
        end
    else
        nflushed += writedata!(stream.stream, buffer)
        makemargin!(buffer, 0)
    end
    buffer.transcoded += nflushed
    return nflushed
end

function flushuntilend(stream::NoopStream)
    stream.state.buffer1.transcoded += writedata!(stream.stream, stream.state.buffer1)
    return
end


text:
# first idea of implementation
# Very basic example, with signals, only with updatePosition!
# two frames are externaly driven

module Test_BasicExBox2
using Modia3D

vgreen  = Modia3D.Material(color="Green", transparency=0.5)
vyellow = Modia3D.Material(color="Yellow", transparency=0.5)

lx = 1.0
ly = 0.5
lz = 0.3

@enum WaveType sinus cosinus sinush cosinush

@assembly Signals begin
  sinWave = Wave(;wave=sinus)
  cosWave = Wave(;wave=cosinus)
end

mutable struct Wave
  waveType::WaveType
  amplitude::Float64
  frequency::Float64
  phaseShift::Float64
  output
  function Wave(; wave::WaveType=sinus, amplitude::Float64=1.0, frequency::Float64=1/(2*pi), phaseShift::Float64=0.0)
    new(wave,amplitude,frequency,phaseShift)
  end
end

function setSignal!(signals::Signals, time)
  signals.sinWave.output = setOutput(signals.sinWave,time)
  signals.cosWave.output = setOutput(signals.cosWave,time)
end

function setOutput(wave::Wave, time)
  if wave.waveType == sinus
    return wave.amplitude*sin(2*pi*wave.frequency*time + wave.phaseShift)
  elseif wave.waveType == cosinus
    return wave.amplitude*cos(2*pi*wave.frequency*time + wave.phaseShift)
  end
end

@assembly Pendulum(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D()
   green = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vgreen),[ [-lx/2,0,0],[ lx/2,0,0] ])
   rev1 = Modia3D.Revolute(world, green.frames[1])
   gelb = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vyellow),[ [-lx/2,0,0],[ lx/2,0,0] ])
   rev2 = Modia3D.Revolute(world, gelb.frames[1])
end
pendulum = Pendulum()
waves = Signals()

Modia3D.initAnalysis!(pendulum, Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3, enableContactDetection=false))

for time = range(0.0, stop=2*pi, length=101)
  setSignal!(waves,time)
  Modia3D.setAngle!(pendulum.rev1, waves.sinWave.output)
  Modia3D.setAngle!(pendulum.rev2, waves.cosWave.output)
  Modia3D.updatePosition!(pendulum)

  Modia3D.visualize!(pendulum,time)
end
Modia3D.closeAnalysis!(pendulum)

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


text:
module SampleConventionalApp

function run_app()
	@info("Running application...")
	println("Hello world!")
end

end # module

max_stars_repo_name: ma-laforge/SampleConventionalApp

max_issues_repo_name: ma-laforge/SampleConventionalApp

max_forks_repo_name: ma-laforge/SampleConventionalApp


text:
"""
This package is intended to provide an abstract interface for working
with tree structures.
Though the package itself is not particularly sophisticated, it defines
the interface that can be used by other packages to talk about trees.
"""
module AbstractTrees

export children, ischild, intree, isdescendant, treesize, treebreadth, treeheight
export print_tree, TreeCharSet
export TreeIterator, Leaves, PostOrderDFS, PreOrderDFS, StatelessBFS
export Tree, ShadowTree, AnnotationNode, treemap, treemap!

import Base: getindex, setindex!, iterate, nextind, print, show,
    eltype, IteratorSize, IteratorEltype, length, push!, pop!
using Base: SizeUnknown, EltypeUnknown

abstract type AbstractShadowTree end

include("traits.jl")
include("base.jl")
include("implicitstacks.jl")
include("printing.jl")
include("indexing.jl")
include("iteration.jl")
include("builtins.jl")

end # module


text:
# FIXME: Add link to https://github.com/paranumal/libparanumal here and in
# advection (also update the license)

# {{{ Fill sendbuf on device with buf (for all dimensions)
function knl_fillsendbuf!(::Val{Np}, ::Val{nvar}, sendbuf, buf,
                          sendelems) where {Np, nvar}
  n = threadIdx().x
  e = blockIdx().x

  @inbounds if n <= Np && e <= length(sendelems)
    re = sendelems[e]
    for s = 1:nvar
      sendbuf[n, s, e] = buf[n, s, re]
    end
  end
  nothing
end
# }}}

# {{{ Fill buf on device with recvbuf (for all dimensions)
function knl_transferrecvbuf!(::Val{Np}, ::Val{nvar}, buf, recvbuf, nelem,
                              nrealelem) where {Np, nvar}
  n = threadIdx().x
  e = blockIdx().x

  @inbounds if n <= Np && e <= nelem
    for s = 1:nvar
      buf[n, s, nrealelem + e] = recvbuf[n, s, e]
    end
  end
  nothing
end
# }}}

# {{{ MPI Buffer handling
function fillsendbuf!(sendbuf, d_sendbuf::CuArray, d_buf::CuArray, d_sendelems)
  nsendelem = length(d_sendelems)
  Np = size(d_buf, 1)
  nvar = size(d_buf, 2)
  if nsendelem > 0
    @cuda(threads=Np, blocks=nsendelem,
          knl_fillsendbuf!(Val(Np), Val(nvar), d_sendbuf, d_buf, d_sendelems))
    sendbuf .= d_sendbuf
  end
end

function transferrecvbuf!(d_recvbuf::CuArray, recvbuf, d_buf::CuArray, nrealelem)
  nrecvelem = size(recvbuf)[end]
  Np = size(d_buf, 1)
  nvar = size(d_buf, 2)
  if nrecvelem > 0
    d_recvbuf .= recvbuf
    @cuda(threads=Np, blocks=nrecvelem,
          knl_transferrecvbuf!(Val(Np), Val(nvar), d_buf, d_recvbuf,
                               nrecvelem, nrealelem))
  end
end
# }}}

# {{{ Kernel wrappers
function volumegrad!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                     d_grad::CuArray, d_Q, d_vgeo, gravity, d_D,
                     elems) where {dim, N, nmoist, ntrace}
  ngrad = _nstategrad + 3nmoist
  Qshape    = (ntuple(j->N+1, dim)..., size(d_Q, 2), size(d_Q, 3))
  gradshape = (ntuple(j->N+1, dim)..., ngrad, size(d_Q, 3))
  vgeoshape = (ntuple(j->N+1, dim)..., size(d_vgeo, 2), size(d_Q, 3))

  d_gradC = reshape(d_grad, gradshape)
  d_QC = reshape(d_Q, Qshape)
  d_vgeoC = reshape(d_vgeo, vgeoshape)

  nelem = length(elems)
  @cuda(threads=ntuple(j->N+1, dim), blocks=nelem,
        knl_volumegrad!(Val(dim), Val(N), Val(nmoist), Val(ntrace), d_gradC,
                        d_QC, d_vgeoC, gravity, d_D, nelem))
end

function facegrad!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                   d_grad::CuArray, d_Q, d_vgeo, d_sgeo, gravity, elems,
                   d_vmapM, d_vmapP, d_elemtobndy) where {dim, N, nmoist,
                                                          ntrace}
  nelem = length(elems)
  @cuda(threads=(ntuple(j->N+1, dim-1)..., 1), blocks=nelem,
        knl_facegrad!(Val(dim), Val(N), Val(nmoist), Val(ntrace), d_grad, d_Q,
                      d_vgeo, d_sgeo, gravity, nelem, d_vmapM, d_vmapP,
                      d_elemtobndy))
end

function volumerhs!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                    d_rhs::CuArray, d_Q, d_grad, d_vgeo, gravity, viscosity,
                    d_D, elems) where {dim, N, nmoist, ntrace}
  ngrad = _nstategrad + 3nmoist

  Qshape    = (ntuple(j->N+1, dim)..., size(d_Q, 2), size(d_Q, 3))
  gradshape = (ntuple(j->N+1, dim)..., ngrad, size(d_Q, 3))
  vgeoshape = (ntuple(j->N+1, dim)..., size(d_vgeo,2), size(d_Q, 3))

  d_rhsC = reshape(d_rhs, Qshape...)
  d_QC = reshape(d_Q, Qshape)
  d_gradC = reshape(d_grad, gradshape)
  d_vgeoC = reshape(d_vgeo, vgeoshape)

  nelem = length(elems)
  @cuda(threads=ntuple(j->N+1, dim), blocks=nelem,
        knl_volumerhs!(Val(dim), Val(N), Val(nmoist), Val(ntrace), d_rhsC, d_QC,
                       d_gradC, d_vgeoC, gravity, viscosity, d_D, nelem))
end

function facerhs!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                  d_rhs::CuArray, d_Q, d_grad, d_vgeo, d_sgeo, gravity,
                  viscosity, elems, d_vmapM, d_vmapP,
                  d_elemtobndy) where {dim, N, nmoist, ntrace}
  nelem = length(elems)
  @cuda(threads=(ntuple(j->N+1, dim-1)..., 1), blocks=nelem,
        knl_facerhs!(Val(dim), Val(N), Val(nmoist), Val(ntrace), d_rhs, d_Q,
                     d_grad, d_vgeo, d_sgeo, gravity, viscosity, nelem,
                     d_vmapM, d_vmapP, d_elemtobndy))
end

# }}}

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
# These functions relate to both an ansiotropic and an isotropic grid

function smolyak_grid(node_type::Function,d::S,mu::Union{S,Array{S,1}}) where {S<:Integer}

  T = typeof(1.0)

  multi_index        = generate_multi_index(d,mu)
  unique_multi_index = sort(unique(multi_index))
  unique_node_number = m_i(unique_multi_index)

  # Create base nodes to be used in the sparse grid

  base_nodes   = Array{Array{T,1},1}(undef,length(unique_node_number))
  base_weights = Array{Array{T,1},1}(undef,length(unique_node_number))
  for i = 1:length(unique_node_number)
    base_nodes[i], base_weights[i] = node_type(unique_node_number[i])
  end

  # Determine the unique nodes introduced at each higher level

  unique_base_nodes = Array{Array{T,1},1}(undef,length(unique_node_number))
  unique_base_nodes[1] = base_nodes[1]
  for i = 2:length(unique_base_nodes)
    unique_base_nodes[i] = setdiff(base_nodes[i],base_nodes[i-1])
  end

  # Construct the sparse grid from the unique nodes

  nodes = Array{T,2}(undef,determine_grid_size(multi_index))
  l = 1
  @inbounds for j = 1:size(multi_index,1)
    new_nodes = unique_base_nodes[multi_index[j,1]]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,unique_base_nodes[multi_index[j,i]])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  # Eventually this function should also return the weights at each node on the grid
  # so that it can be used for numerical integration.

  if d == 1
    nodes = nodes[:]
  end

  return nodes, multi_index

end

function smolyak_grid(node_type::Function,d::S,mu::Union{S,Array{S,1}},domain::Union{Array{T,1},Array{T,2}}) where {S<:Integer, T<:AbstractFloat}

  if size(domain,2) != d
    error("domain is inconsistent with the number of dimensions")
  end

  (nodes, multi_index) = smolyak_grid(node_type,d,mu)

  # Now scale the nodes to the desired domain

  nodes = scale_nodes(nodes,domain)

  return nodes, multi_index

  # Eventually this function should also return the weights at each node on the grid
  # so that it can be used for numerical integration.

end

function determine_grid_size(mi)

  temp = similar(mi)

  for i = 1:size(mi,1)
    for j = 1:size(mi,2)
      if mi[i,j] == 1
        temp[i,j] = 1
      elseif mi[i,j] == 2
        temp[i,j] = 2^(mi[i,j]-1)
      else
        temp[i,j] = 2^(mi[i,j]-1)+1 - (2^(mi[i,j]-2)+1)
      end
    end
  end

  s = 0
  for i = 1:size(mi,1)
    t = 1
    for j = 1:size(mi,2)
      t *= temp[i,j]
    end
    s += t
  end

  return (s, size(mi,2))

end


text:
module Manifests

import ..ManifestFormatEntry
import ..hpf_manifest_format
import ..base_path, ..base_path_default, ..manifest_subdir_default, ..data_subdir_default

export read
export get_data_filenames
export get_base_path, set_base_path!, get_manifest_path, set_manifest_path!, get_data_path, set_data_path!
#export get_format, set_format!

using DataFrames
using Printf

# Specify paths to manifest and data files
manifest_path_default = joinpath(base_path_default,manifest_subdir_default)
data_path_default = joinpath(base_path_default,data_subdir_default)

get_manifest_path() = manifest_path_default
set_manifest_path!(path::String) = manifest_path_default = path

get_data_path() = data_path_default
set_data_path!(path::String) = data_path_default = path

get_base_path() = base_path
function set_base_path!(path::String)
   base_path = path
   set_manifest_path!( joinpath(base_path,manifest_subdir_default) )
   set_data_path!( joinpath(base_path,data_subdir_default) )
end

# Initialize base_path
set_base_path!(base_path_default)

#manifest_format = hpf_manifest_format
#get_format() = manifest_format
#set_format!(format::Vector{ManifestFormatEntry}) = manifest_formation = format

# Read in manifest files
function read(filename::String; filter_object::Regex = r"", filter_obstype::Regex = r"", verbose=false, manifest_format = hpf_manifest_format )
    df_types = map(t->Union{t,Missing},getproperty.(manifest_format,:type) )
    manifest = DataFrame(df_types,getproperty.(manifest_format,:name),0)

    manifest_lines = readlines(filename)
    for line in manifest_lines[2:end]
        line_parsed = Dict()
        for entry in manifest_format
            name = entry.name
            type = entry.type
            if entry.type == String
                line_parsed[entry.name] = strip(line[entry.cols])
            elseif entry.type <: Number
                val = tryparse(entry.type,line[entry.cols])
                if val == nothing
                   val = missing
                end
                line_parsed[entry.name] = val
            else
                @error("What am I supposed to do with an entry of type" * string(entry.type) * "?")
            end
        end # for entry
        if occursin(filter_object,line_parsed[:Object]) && occursin(filter_obstype,line_parsed[:ObsType])
            push!(manifest,line_parsed)
        end
    end # for line
   return manifest
end

function read(; path::String = get_manifest_path(), filter_files::Regex = r"", filter_object::Regex = r"", filter_obstype::Regex = r"", manifest_format = hpf_manifest_format, verbose=false )
    df_types = map(t->Union{t,Missing},getproperty.(manifest_format,:type) )
    manifest = DataFrame(df_types,getproperty.(manifest_format,:name),0)
    manifest_files = readdir(path)
    for manifest_filename in manifest_files
        if ! occursin(filter_files,manifest_filename )
            continue
        end
        manifest_this = read(joinpath(path,manifest_filename), filter_object=filter_object, filter_obstype=filter_obstype, verbose=verbose)
        append!(manifest,manifest_this)
    end
    return manifest
end

"""Get data filenames from a manifest DataFrame.  Optioanlly specify path to data."""
function get_data_filenames(df::DataFrame; path::String = get_data_path(), verbose=false )
    num_files = size(df,1)
    files = Vector{String}()
    for i in 1:num_files
       obsnum = Printf.@sprintf("%04d",df[i,:ObsNum])
       subdir = joinpath(df[i,:UTDate],obsnum)
       dirname = joinpath(path,subdir)
       file_pattern = Regex("(" * df[i,:Frame][1:end-5] * ".*fits)")
       local files_here
       try
          files_here = readdir(dirname)
       catch
          if verbose
             println("# Skipping " * string(file_pattern) * " not found in " * subdir )
          end
          continue
       end
       files_matched = map(f->match(file_pattern,f),files_here)
       idx_keep = findall(x->x!=nothing,files_matched)
       append!(files,map(m->joinpath(dirname,m.match),files_matched[idx_keep]))

    end
    return files
end

"""Get data filenames from searching all manifest entries.  Optioanlly specify path to data, filter for manifest filenames, and/or filter for object names."""
function get_data_filenames(; manifest_path::String = get_manifest_path(), data_path::String = get_data_path(), filter_files::Regex = r"", filter_object::Regex = r"", filter_obstype::Regex = r"", manifest_format = hpf_manifest_format, verbose = false  )
     df = read(path=manifest_path, filter_files=filter_files, filter_object=filter_object, filter_obstype=filter_obstype, manifest_format=manifest_format, verbose=verbose)
     get_data_filenames(df, path=data_path, verbose=verbose)
end

end


text:
abstract type CompiledState <: State end

function generate_field_type(domain::Domain, sig::Signature{N}) where {N}
    dtype = get(get_datatypes(domain), sig.type, GLOBAL_DATATYPES[sig.type])
    return N == 0 ? dtype : (dtype == Bool ? BitArray{N} : Array{dtype,N})
end

function generate_field_type(domain::AbstractedDomain, sig::Signature{N}) where {N}
    dtype = domain.interpreter.abstractions[sig.type]
    return N == 0 ? dtype : Array{dtype,N}
end

function generate_pred_init(domain::Domain, state::State,
                            sig::Signature{N}) where {N}
    if N == 0 return false end
    dims = generate_fluent_dims(domain, state, sig)
    return :(falses($(dims...)))
end

function generate_func_init(domain::Domain, state::State,
                            sig::Signature{N}) where {N}
    default = QuoteNode(defaultval(Val(sig.type)))
    if N == 0 return default end
    dims = generate_fluent_dims(domain, state, sig)
    return :(fill($default, $(dims...)))
end

function generate_func_init(domain::AbstractedDomain, state::State,
                            sig::Signature{N}) where {N}
    abstype = domain.interpreter.abstractions[sig.type]
    default = QuoteNode(abstype(defaultval(Val(sig.type))))
    if N == 0 return default end
    dims = generate_fluent_dims(domain, state, sig)
    return :(fill($default, $(dims...)))
end

function generate_state_type(domain::Domain, state::State, domain_type::Symbol)
    # Generate type definition
    state_fields = Expr[]
    for (_, pred) in sortedpairs(get_predicates(domain))
        type = generate_field_type(domain, pred)
        field = Expr(:(::), pred.name, QuoteNode(type))
        push!(state_fields, field)
    end
    for (_, fn) in sortedpairs(get_functions(domain))
        type = generate_field_type(domain, fn)
        field = Expr(:(::), fn.name, QuoteNode(type))
        push!(state_fields, field)
    end
    if domain isa AbstractedDomain
        name = "CompiledAbstracted" * pddl_to_type_name(get_name(domain)) * "State"
    else
        name = "Compiled" * pddl_to_type_name(get_name(domain)) * "State"
    end
    state_type = gensym(name)
    state_typesig = Expr(:(<:), state_type, QuoteNode(CompiledState))
    state_typedef = :(@auto_hash_equals $(Expr(:struct, true, state_typesig,
                                               Expr(:block, state_fields...))))
    state_constructor_defs =
        generate_state_constructors(domain, state, domain_type, state_type)
    state_method_defs =
        generate_state_methods(domain, state, domain_type, state_type)
    state_defs = Expr(:block, state_constructor_defs, state_method_defs)
    return (state_type, state_typedef, state_defs)
end

function generate_state_constructors(domain::Domain, state::State,
                                     domain_type::Symbol, state_type::Symbol)
    # Generate constructor with no arguments
    state_inits = []
    state_copies = Expr[]
    for (name, pred) in sortedpairs(get_predicates(domain))
        push!(state_inits, generate_pred_init(domain, state, pred))
        push!(state_copies, :(copy(state.$name)))
    end
    for (name, fn) in sortedpairs(get_functions(domain))
        push!(state_inits, generate_func_init(domain, state, fn))
        push!(state_copies, :(copy(state.$name)))
    end
    state_constructor_defs = quote
        $state_type() = $state_type($(state_inits...))
        $state_type(state::$state_type) = $state_type($(state_copies...))
        function $state_type(state::State)
            new = $state_type()
            for (term, val) in get_fluents(state)
                if val === false continue end
                set_fluent!(new, val, term)
            end
            return new
        end
    end
    return state_constructor_defs
end

function generate_state_methods(domain::Domain, state::State,
                                domain_type::Symbol, state_type::Symbol)
    # Construct domaintype and statetype methods
    types_def = quote
        domaintype(::Type{$state_type}) = $domain_type
        domaintype(::$state_type) = $domain_type
        statetype(::Type{$domain_type}) = $state_type
        statetype(::$domain_type) = $state_type
        get_domain(::$state_type) = $domain_type()
    end
    # Fluent name listing
    get_fluent_names_def = quote
        function get_fluent_names(state::$state_type)
            domain = get_domain(state)
            f = name -> begin
                grounded_args = groundargs(domain, state, name)
                return (!isempty(args) ? Compound(name, collect(args)) :
                        Const(name) for args in grounded_args)
            end
            fluents = Base.Generator(f, fieldnames($state_type))
            return Iterators.flatten(fluents)
        end
    end
    # Fluent accessors
    fluent_conds, get_fluent_brs, set_fluent_brs, groundargs_defs =
        Expr[], Expr[], Expr[], Expr[]
    for (name, sig) in pairs(get_fluents(domain))
        # Ground args
        def = generate_groundargs(domain, state, domain_type, state_type, name)
        push!(groundargs_defs, def)
        # Generate accessor branch conditions
        push!(fluent_conds, :(term.name == $(QuoteNode(name))))
        # Generate accessor branch defs for 0-arity fluents
        if length(sig.args) == 0
            push!(get_fluent_brs, :(state.$name))
            push!(set_fluent_brs, :(state.$name = val))
            continue
        end
        # Generate accessor branch defs multiple-arity fluents
        term = convert(Term, sig)
        varmap = Dict(a => :(term.args[$i].name)
                      for (i, a) in enumerate(sig.args))
        idxs = generate_fluent_ids(domain, state, term, sig, varmap)
        push!(get_fluent_brs, :(@inbounds state.$name[$(idxs...)]))
        push!(set_fluent_brs, :(@inbounds state.$name[$(idxs...)] = val))
    end
    err_br = :(error("Unrecognized fluent: $(term.name)"))
    get_fluent_def = quote
        function get_fluent(state::$state_type, term::Term)
            return $(generate_switch_stmt(fluent_conds, get_fluent_brs, err_br))
        end
    end
    set_fluent_def = quote
        function set_fluent!(state::$state_type, val, term::Term)
            return $(generate_switch_stmt(fluent_conds, set_fluent_brs, err_br))
        end
    end
    state_method_defs = Expr(:block, types_def,
        get_fluent_names_def, groundargs_defs...,
        get_fluent_def, set_fluent_def)
    return state_method_defs
end

function generate_groundargs(domain::Domain, state::State,
                             domain_type::Symbol, state_type::Symbol,
                             fluent::Symbol)
    sig = get_fluent(domain, fluent)
    argtypes = sig.argtypes
    if get_requirements(domain)[:typing]
        objs_exprs = [:(get_objects(state, $(QuoteNode(ty)))) for ty in argtypes]
    else
        objs_exprs = [:(get_objects(state)) for ty in argtypes]
    end
    iter_expr = :(Iterators.product($(objs_exprs...)))
    valtype = QuoteNode(Val{fluent})
    groundargs_def = quote
        function groundargs(domain::$domain_type, state::$state_type,
                            fluent::$valtype)
            return $iter_expr
        end
    end
    return groundargs_def
end

function (::Type{S})(state::State) where {S <: CompiledState}
    new = S()
    for (term, val) in get_fluents(state)
        if val === false continue end
        set_fluent!(new, val, term)
    end
    return new
end

Base.copy(state::S) where {S <: CompiledState} =
    S(state)

groundargs(domain::CompiledDomain, state::State, fluent::Symbol) =
    groundargs(domain, state, Val(fluent))

get_fluents(state::CompiledState) =
    (term => get_fluent(state, term) for term in get_fluent_names(state))

get_facts(state::CompiledState) =
    (term for term in get_fluent_names(state)
     if get_fluent(state, term) in (true, both))


text:
using Test

function Λ1(x, T)
    return (T-x)/T
end

#add basis for the final value
function Λ2(x,T)
    return x/T
end

function Λ(x, T)
    return max(0.0, sqrt(T)/2 - abs(x/sqrt(T) - sqrt(T)/2))
end
#T = 3
#plot([Λ(x, T) for x=0:0.001:T])

function Λ(x, l, k, T)
    x = mod(x, T)   #ok
    2^(-l/2)*Λ((x)*(2)^(l) - k*T, T)
end

function c(k::Fs, ξ, t, u, v, T)
    c = 0.0
    for n in k.nhb
        i,j = Faber(n)
        c += Λ(t, i, j, T)*ξ[n]
    end
    return c + Λ1(t, T)*u + Λ2(t, T)*v
end

L = 1
T = 10.0
ϕ = Fs(1, L, T)
ξ = fill(1.0, 2<<L - 1)
u = -1.0
v = -1.0

S = System(L, T, ξ)

fs_expansion(S, t, u, v)
c(ϕ, ξ, t, u, v, T)
for k in S.ϕ
    @test MCintegration(k) < k.ub
    @test MCintegration(k) > k.lb
    @test k.lb + range = k.ub
end

fs_expansion()
fs_expansion(ϕ, ξ, u, v, L, T)

for i
function MCintegration(k::Fs)
end

using Test
for i in 1:2^20
    @test i >= 1 << binlog(i)
    @test i < 2 << binlog(i)
    @test binlog(i) == Faber(i)[1]
end

function Λ(x, l, k, T)
   x = mod(x, T)   #ok
   2^(-l/2)*Λ((x)*1<<(l) -T*(k),T)
end

using Plots

T=1
dt= 0:0.0001:T
plot(dt, [Λ(t, 2, 0, T) for t in dt])

function num_int(i1,j1,i2,j2, T, dt)
    t = 0:dt:T
    r = 0
    for i in t[2:end]
        r += Λ(i + dt/2, i1, j1, T)*Λ(i + dt/2, i2, j2, T)*dt
    end
    return r
end

function num_int(i1,j1, T, dt)
    t = 0:dt:T
    r = 0
    for i in t[2:end]
        r += Λ(i + dt/2, i1, j1, T)*dt
    end
    return r
end

function num_int_bound1(n, T, dt)
    i1,j1 = Faber(n)
    t = 0:dt:T
    r = 0.0
    for i in t[2:end]
        r += Λ(i + dt/2, i1, j1, T)*(i + dt/2)/T*dt
    end
    return r
end

function num_int_bound2(n, T, dt)
    i1,j1 = Faber(n)
    t = 0:dt:T
    r = 0.0
    for i in t[2:end]
        r += Λ(i + dt/2, i1, j1, T)*(1 -(i + dt/2)/T)*dt
    end
    return r
end

num_int(0, 0, 2, 2 , 100, 0.00001)
num_int(0, 0 , 10, 0.00001)

generate_matrix(2, 10)
generate_vector(2, 10)

Λ(2.5, 1, 0, 1)
sqrt(10)*2^-(1/2 + 1)

bound2(3, 100)
num_int_bound2(3, 100.0, 0.001)

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
#= DirectToXYZ.jl
this script converts a matrix of direct, fractional coordinates
such as you might find in a POSCAR (VASP) file
 into standard cartesian coordinates

 INPUT: files: lattice-vectors, direct-list, atomlist
               & it will ask you for the value of Factor
 Files: lattice-vectors contains lattice vectors, a 3x3 list of their xyz components
        direct-list is a nx3 list of direct coordinates
 Factor = the scaling factor for the fractional coordinates
 OUTPUT: file coords.txt containing standard cartesian coordinates

 type info: All numbers are Float64.
=#
# get the scaling factor for the fractional/direct coordinates.
println("What is the scaling factor? Enter a Float64")
Factor=parse(Float64, readline())
println("Factor=", Factor)

# get the atomlist file into an array
AtomList=readdlm("atomlist")

# set the Lattice Vectors a, b, and c. Each is a 1x3 vector.
# having components e.g. ax, ay, az.

LatticeVectors=readdlm("lattice-vectors")
a=LatticeVectors[1:1,:]
b=LatticeVectors[2:2,:]
c=LatticeVectors[3:3,:]
println("...the lattice vectors are:")
println("a=", a)
println("b=", b)
println("c=", c)

DirectList=readdlm("direct-list")

n=size(DirectList,1)
XList=DirectList[1:n,1]
YList=DirectList[1:n,2]
ZList=DirectList[1:n,3]

# this is the old version of the functions. they were all identical except
# for interchanging the name of the coordinate.
#function UnscaleX(Factor::Float64, XList::Array{Float64}, a::Array{Float64})
# CarteX=zeros(XList)
# n=size(XList,1)
# for i in 1:n
#  CarteX[i]= Factor*(XList[i]*a[1] + XList[i]*a[2] + XList[i]*a[3])
# end
# CarteX
#end

function UnscaleX(Factor::Float64, XList::Array{Float64}, YList::Array{Float64}, ZList::Array{Float64}, a::Array{Float64}, b::Array{Float64}, c::Array{Float64})
 CarteX=zeros(XList)
 n=size(XList, 1)
 for i in 1:n
   CarteX[i]=Factor*(XList[i]*a[1] + YList[i]*b[1] + ZList[i]*c[1])
 end
 CarteX
end

function UnscaleY(Factor::Float64, XList::Array{Float64}, YList::Array{Float64}, ZList::Array{Float64}, a::Array{Float64}, b::Array{Float64}, c::Array{Float64})
 CarteY=zeros(YList)
 n=size(YList, 1)
 for i in 1:n
  CarteY[i]=Factor*(XList[i]*a[2] + YList[i]*b[2] + ZList[i]*c[2])
 end
 CarteY
end

function UnscaleZ(Factor::Float64, XList::Array{Float64}, YList::Array{Float64}, ZList::Array{Float64}, a::Array{Float64}, b::Array{Float64}, c::Array{Float64})
 CarteZ=zeros(ZList)
 n=size(ZList,1)
 for i in 1:n
  CarteZ[i]=Factor*(XList[i]*a[3] + YList[i]*b[3] + ZList[i]*c[3])
 end
 CarteZ
end

FinalX=zeros(XList)
FinalY=zeros(YList)
FinalZ=zeros(ZList)

FinalX=UnscaleX(Factor, XList, YList, ZList, a, b, c)
FinalY=UnscaleY(Factor, XList, YList, ZList, a, b, c)
FinalZ=UnscaleZ(Factor, XList, YList, ZList, a, b, c)

# Format results and save to an XYZ file
n=size(YList, 1)
println("Enter the text (< 80 chars) for the comment line of the xyz file:")
println("This is commonly the empirical formula of the system, or something descriptive")
CommentLine=readline()
write("commentline.txt", "$CommentLine")
writedlm("natom.txt", n)
writedlm("coords.txt", [AtomList FinalX FinalY FinalZ], '\t')
println("coordinates have been converted to XYZ format and saved in coords2.xyz.")
run(pipeline(`cat natom.txt commentline.txt coords.txt`, stdout="coords2.xyz"))

# clean up temp and .txt files before exiting
  run(`\rm natom.txt commentline.txt coords.txt`)
# run(`\rm lattice-vectors direct-list atomlist`)
#


text:
module MueLu
  using CxxWrap, MPI
  import ..libjltrilinos
  import ..Teuchos
  import ..Tpetra

  @wrapmodule(libjltrilinos, :register_muelu)

  """
  Return the "easy" parameters from the MueLu tutorial
  """
  function parameters_easy()
    pl = Teuchos.ParameterList("MueLu")
    parameters_easy!(pl)
    return pl
  end
  function parameters_easy!(pl)
    pl["verbosity"] = "low"
    pl["max levels"] = 3
    pl["coarse: max size"] = 10
    pl["multigrid algorithm"] = "sa"

    pl["smoother: type"] = "RELAXATION"
    smoother_params = Teuchos.sublist(pl, "smoother: params")
    smoother_params["relaxation: type"] = "Jacobi"
    smoother_params["relaxation: sweeps"] = 1
    smoother_params["relaxation: damping factor"] = 0.9

    pl["aggregation: type"] = "uncoupled"
    pl["aggregation: min agg size"] = 3
    pl["aggregation: max agg size"] = 9

    return
  end

end


text:
using WordTokenizers
using Test

files = ["simple",
         "sedbased",
         "sentence_splitting",
         "set_method_api",
         "tweet_tokenize",
         "reversible_tok",
         "toktok",
         "sp_unigram"
        ]

@testset "$file" for file in files
    include(file * ".jl")
end


text:
@testset "star attack" begin
    occ = [0,0,0,0,0,0,0,0,
           0,0,1,0,0,0,1,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,1,0,
           0,1,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]

    att = [0,0,0,0,1,0,0,0,
           0,0,0,0,1,0,0,0,
           0,0,0,1,1,1,0,0,
           1,1,1,1,0,1,1,1,
           0,0,0,1,1,1,0,0,
           0,0,1,0,1,0,0,0,
           0,0,0,0,1,0,0,0,
           0,0,0,0,1,0,0,0]

    edg = [0,0,0,0,0,0,0,0,
           0,0,1,0,0,0,1,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,1,0,
           0,1,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]
    occ_uint = Bobby.cvt_to_uint(BitArray(occ))
    ui = Bobby.INT2UINT[29]
    m, e = Bobby.star_attack(occ_uint, ui)
    @test all(Int.(Bobby.cvt_to_bitarray(m)) .== att)
    @test all(Int.(Bobby.cvt_to_bitarray(e)) .== edg)
end


text:
module Triangulate

type triangulateio64
  pointlist::Array{Float64,1}
  numberofpoints::Int32
  trianglelist::Array{Int32,1}
  numberoftriangles::Int32
  edgelist::Array{Int32,1}
  numberofedges::Int32
end

const TrianglelibPath = abspath(joinpath(splitdir(Base.source_path())[1],"..","lib","libtriangle.so"))

function triangulate(size::Float64,numberofpoints::Int64,pointlist::Array{Float64, 1})
  if numberofpoints * 2 != length(pointlist)
    error("Triangulate: Inconsistent Mesh Input.")
  end
  ret = ccall((:triangleio, TrianglelibPath), Any, (Ptr{Void}, Ptr{UInt8}, Ptr{UInt8}, Int64, Ptr{Float64}),(pointer_from_objref(Triangulate)),"QpcAevn",string("Qprq34.0","a",size,"BPe"), Int32(numberofpoints), pointlist)
  return ret
end

function UnitSquareMesh(tol::Float64)
  ret = triangulate(tol, 4, [0., 0., 1., 0., 1., 1., 0., 1.])
  return ret
end

function UnitCircle(tol::Float64, size::Int64)
  Θ = 2 * π/size
  pointlist = zeros(2 * size)
  for idx in 1:size
    pointlist[2 * idx - 1] = cos(idx * Θ)
    pointlist[2 * idx    ] = sin(idx * Θ)
  end
  ret = triangulate(tol, size, pointlist)
end

export triangulate, UnitSquareMesh, UnitCircle

end


text:
struct TensorInformation
    indices::Vector{Tuple{Symbol,Int}}
    isconj::Bool
end

function parse_tensors_rhs(s::AbstractString)
    rhs_matched = collect(eachmatch(
        r"\[\s*@??\s*(?<label>[0-9]+)\s*(?<isconj>[\*⋆]??)\s*\|" *
        r"(?<index>[,_\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\s]+?)\]",
        s;
        overlap = false,
    ))
    isempty(rhs_matched) && error("no argument tensors")
    tensors_rhs = Vector{TensorInformation}(undef, length(rhs_matched))
    for x in rhs_matched
        label = parse(Int, x[:label])
        isconj = !isempty(x[:isconj])
        indices = Vector{Tuple{Symbol,Int}}()
        for n in split(x[:index], ",")
            sym = Meta.parse(n)
            sym isa Symbol || error("invalid index syntax @$(label)")
            push!(indices, (sym, label))
        end
        tensors_rhs[label] = TensorInformation(indices, isconj)
    end
    isassigned(tensors_rhs) || error("labels of tensors should be in 1:(#tensor)")
    return tensors_rhs
end

function parse_indices_lhs(s::AbstractString)
    lhs_matched = collect(eachmatch(
        r"\[\s*@??\s*>\s*\|(?<index>[0-9@,_\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\s]+?)\]",
        s;
        overlap = false,
    ))
    isempty(lhs_matched) && error("no definition of return tensor")
    length(lhs_matched) > 1 && error("duplicate definition of return tensor")
    lhs_matched = first(lhs_matched)
    indices_lhs = Vector{Tuple{Symbol,Int}}()
    for n in split(lhs_matched[:index], ",")
        index, label = split(n, "@")
        sym = Meta.parse(index)
        sym isa Symbol || error("invalud index syntax")
        label = parse(Int, label)
        push!(indices_lhs, (sym, label))
    end
    return indices_lhs
end

function parse_pairs_index(s::AbstractString)
    pairs_matched = collect(eachmatch(
        r"\(\s*(?<kind>[_\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}]+)\s*@\s*(?<klab>[0-9]+)\s*" *
        r"\s*=\s*(?<vind>[_\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\s]+)\s*@\s*(?<vlab>[0-9]+)\s*\)",
        s;
        overlap = false,
    ))
    pairs_index = Dict{Tuple{Symbol,Int},Tuple{Symbol,Int}}()
    for bond in pairs_matched
        kind = Meta.parse(bond[:kind])
        klab = parse(Int, bond[:klab])
        kind isa Symbol || error("invalid index syntax")
        vind = Meta.parse(bond[:vind])
        vlab = parse(Int, bond[:vlab])
        vind isa Symbol || error("invalid index syntax")
        pairs_index[(kind, klab)] = (vind, vlab)
    end
    length(union(keys(pairs_index), values(pairs_index))) == 2 * length(pairs_matched) ||
        error("duplication in index pairs")
    return pairs_index
end

function parse_tree_contract(s::AbstractString)
    t = collect(eachmatch(r"\[(?<tree>[0-9\s,\(\)]+?)\]", s; overlap = false))
    isempty(t) && return nothing
    length(t) > 1 && error("duplicate definition of contraction tree")
    return Meta.eval(Meta.parse(t[1][:tree]))
end


text:
# cheers_1 = NVModel(1, 5, truncated(Normal(90, 30), 0, 180) )
@test leftover(cheers_1, 86) ≈  9.97386724034834
@test sales(cheers_1, 86) ≈  76.02613275965166
@test lost_sales(cheers_1, 86) ≈  13.97386724034834

@test isapprox(leftover(cheers_1, 200), 110.0)
@test sales(cheers_1, 200) ≈  90.0

@test isapprox(lost_sales(cheers_1, 200),  0.0, atol=1/1_000_000)
@test profit(cheers_1, 86) ≈  294.1306637982583
@test q_opt(cheers_1) ≈ 115
@test q_opt(cheers_1, rounded = false) ≈ 115.16195141374729
@test q_scarf(cheers_1) ≈ 112.19801383255745
@test profit_shift(cheers_1) == 0
@test q_min(cheers_1) == 0
@test q_max(cheers_1) == Inf

# cheers_2 = cheers_2 = NVModel(4.0, 5.5, Uniform(0, 60))
@test q_opt(cheers_2) ≈ 16
@test q_opt(cheers_2, rounded = false) ≈ 16.363636363636363
@test q_scarf(cheers_2) ≈ 21.161165235168156
@test profit_shift(cheers_1) == 0.0

# cheers_4 = NVModel(10, 42, DiscreteNonParametric([0,1,2],[.3,.5,.2]))
@test isapprox(lost_sales(cheers_4, 2),  0.0, atol=1/1_000_000)
@test isapprox(lost_sales(cheers_4, 3),  0.0, atol=1/1_000_000)
@test sales(cheers_4, 2) ≈  0.9
@test leftover(cheers_4, 2) ≈  1.1
@test profit(cheers_4, 2) ≈  17.8
@test q_opt(cheers_4) ≈ 1
@test q_opt(cheers_4, rounded = false) ≈ 1
@test q_scarf(cheers_4) ≈ 1.3304430856687095

# bagel = bagel = NVModel(1, 4, Normal(150, 50))
@test q_opt(bagel) ≈ 184

# xs = vec([5_000	10_000	15_000	20_000	25_000	30_000	35_000	40_000	45_000	50_000	55_000	60_000	65_000	70_000	75_000])
# ps = vec([0.0181	0.0733	0.1467	0.1954	0.1954	0.1563	0.1042	0.0595	0.0298	0.0132	0.0053	0.0019	0.0006	0.0002	0.0001])
# elvis  = NVModel(cost = 6, price = 12, salvage = 2.5, demand = DiscreteNonParametric(xs,ps))
@test lost_sales(elvis, 50_000) ≈  61.0
@test sales(elvis, 50_000) ≈   24939.0
@test leftover(elvis, 50_000) ≈  25_061.0
@test profit(elvis, 75_000) ≈ -25000.0
@test q_opt(elvis) ≈ 30000
@test q_opt(elvis, rounded = false) ≈ 30_000
@test q_scarf(elvis) ≈ 27726.15473567755

# ksweet = NVModel(d = Binomial(2000, 0.02), c = 28, p = 45, salvage = 27)
@test q_opt(ksweet) == 50
@test q_opt(ksweet, rounded = false)  == 50
@test q_scarf(ksweet) ≈ 52.14810563784474
@test profit(ksweet, 50)  ≈  666.8822263692309

# poisson distribution
m = NVModel(Poisson(90), 0, 3, q_min=90, q_max = 120)
@test q_opt(m) == 120.0
@test leftover(m, 120.0) == 30.003682468511524
@test profit(m, 120.0) == 269.98895259446545

# lognormal distribution
x̄ = 1000
var = 600^2
σ = log(var/ x̄^2 + 1)
μ = log(x̄) - σ^2 / 2
flextrola_demand = LogNormal(μ, σ)
flextrola = NVModel(cost = 72,
                    price = 121,
                    salvage = 50,
                    demand = flextrola_demand
                    )
@test q_opt(flextrola) == 1111

solve(flextrola)

## backorder penalty
beer = NVModel(Uniform(0, 300), cost = 0, price = 0, backorder = 1, salvage = -0.5)
@test q_opt(beer) == 200

# holding cost and substitute margin
thesaurus_nvm = NVModel(cost = 1.15, price = 2.75, demand = Normal(18, 6),
                        substitute = 2.75-1.15, backorder = 0.5, salvage = 1.15,
                        holding = 1.15 * 0.2 / 12)
solve(thesaurus_nvm)
@test q_opt(thesaurus_nvm) == 29

## Incpomplete struct

struct IncompleteScenario <: AbstractNewsvendorProblem end
concreteIncompleteScenario = IncompleteScenario()
@test_throws ErrorException underage_cost(concreteIncompleteScenario)
@test_throws ErrorException overage_cost(concreteIncompleteScenario)
@test_throws ErrorException distr(concreteIncompleteScenario)
@test_throws ErrorException profit_shift(concreteIncompleteScenario)
@test_throws ErrorException q_min(concreteIncompleteScenario)
@test q_max(concreteIncompleteScenario) == Inf
@test_throws ErrorException profit(concreteIncompleteScenario)

## Boundary cases

# If there is no penalty in ordering too much
m = NVModel(Normal(90, 30), 0, 3)
@test q_opt(m) == Inf
@test profit(m, q_opt(m)) == 270.00

# If there is no penalty in ordering too much but a limit
m = NVModel(Normal(90, 30), 0, 3, q_min=90, q_max = 120)
@test q_opt(m) == 120.00

# If there is no point in ordering
m = NVModel(Normal(90, 30), 4, 3, fixcost = 100)
@test q_opt(m) == 0
@test profit(m, q_opt(m)) == -100.0343938885343

# If there is no point in ordering but we are forced to
m = NVModel(Normal(90, 30), 4, 3, q_min = 20, fixcost = 100)
@test q_opt(m) == 20
@test profit(m, q_opt(m)) ≈ -120.29875100641834
@test profit(m) ≈ -120.29875100641834

# If rounding exceeds q_min
m = NVModel(Normal(90, 30), 1, 5, q_min = 115.1)
@test q_opt(m, rounded = false) == 115.24863700718743
@test q_opt(m) == q_min(m)

# If rounding exceeds q_max
m = NVModel(Normal(90, 28), 1, 5, q_min = 20, q_max = 113.9)
@test q_opt(m, rounded = false) == 113.56539454004161
@test q_opt(m) == q_max(m)


text:
# Structure.jl

# [[file:~/Documents/Julia/scrap.org::*Structure.jl][Structure.jl:1]]
Base.length(A::UpTuple) = length(A.data)
Base.eachindex(A::UpTuple) = eachindex(A.data)

square(arr::UpTuple) = [arr.data...]' * [arr.data...]
square(a::SymExpr) = a^2

# addition and substraction
for op in [:+, :-]
    @eval Base.$op(x::UpTuple, y::UpTuple) = UpTuple($op(x.data, y.data))
    @eval Base.$op(A::UpTuple, b) = UpTuple([$op(a, b) for a in A.data])
    @eval Base.$op(a, B::UpTuple) = UpTuple([$op(a, b) for b in B.data])
end

Base.:(==)(x::UpTuple, y::UpTuple) = x.data == y.data

Base.iterate(x::UpTuple) = iterate(x.data)
Base.iterate(x::UpTuple, n::Int64) = iterate(x.data, n)

# Structure.jl:1 ends here


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    vcat(𝒟₁, 𝒟₂)

Concatenate spatial domains `𝒟₁` and `𝒟₂` vertically, i.e.
concatenate the coordinates into a new point set containing
all points.
"""
Base.vcat(𝒟₁::AbstractDomain, 𝒟₂::AbstractDomain) =
  PointSet(hcat(coordinates(𝒟₁), coordinates(𝒟₂)))

"""
    vcat(𝒮₁, 𝒮₂)

Concatenate spatial data `𝒮₁` and `𝒮₂` vertically, i.e.
concatenate the underlying domains and additionally the
table of variables.
"""
function Base.vcat(𝒮₁::AbstractData, 𝒮₂::AbstractData)
  𝒯 = vcat(values(𝒮₁), values(𝒮₂), cols=:union)
  𝒟 = vcat(domain(𝒮₁), domain(𝒮₂))
  georef(𝒯, 𝒟)
end

"""
    hcat(𝒮₁, 𝒮₂)

Concatenate spatial data `𝒮₁` and `𝒮₂` horizontally, i.e.
concatenate the columns of the underlying table assuming
`nelms(𝒮₁) == nelms(𝒮₂)`.
"""
function Base.hcat(𝒮₁, 𝒮₂)
  @assert nelms(𝒮₁) == nelms(𝒮₂) "cannot join different number of points"
  𝒯 = hcat(values(𝒮₁), values(𝒮₂), makeunique=true)
  𝒟 = domain(𝒮₁)
  georef(𝒯, 𝒟)
end


text:
module Entropia

include("funciones.jl")

export bandt_and_pompe
export bandt_and_pompe_normal
export weight_entropy
export renyi_entropy
export complexity_entropy

end # module