text:
@testset "Comparisons" begin
  grid2D = RegularGrid{Float64}(100,100)
  problem2D = EstimationProblem(data2D, grid2D, :value)
  solver₁ = Kriging(:value => (variogram=GaussianVariogram(range=35.),))
  solver₂ = Kriging(:value => (variogram=SphericalVariogram(range=35.),))

  @testset "Visual comparison" begin
    if ismaintainer || istravis
      function plot_comparison(fname)
        compare([solver₁, solver₂], problem2D, VisualComparison())
        png(fname)
      end
      refimg = joinpath(datadir,"VisualComparison.png")
      @test test_images(VisualTest(plot_comparison, refimg), popup=!istravis, tol=0.1) |> success
    end
  end

  @testset "Cross-validation" begin
    result = compare([solver₁, solver₂], problem2D, CrossValidation(3))
    errors4var = result.errors4var

    # data set with 3 points + 3-fold CV => 3 error values per solver
    @test length(errors4var[:value][1]) == 3
    @test length(errors4var[:value][2]) == 3

    # number of folds must be smaller than number of points, throw an error
    @test_throws AssertionError compare([solver₁, solver₂], problem2D, CrossValidation())
  end
end


text:
module KMeansController
using DataFrames, Query
using ..KMeans
using JSON

export call

"""
call(message)

Computes k-means from records.

# Arguments
- `message::Dict`: the arguments dict
  - `"records"`: the array of array of rows
  - `"cluster_num"``: cluster number

# Return
the json string of result object

# Examples
```jldoctest
julia> KMeansController.call(Dict("cluster_num" => 2, "records" => [[1.0, 2.0, 0.0], [2.0, 1.0, 2.0], [3.0, 0.0, 9.0]]))
"{\"centers\":[[1.5,1.5,1.0],[3.0,0.0,9.0]],\"clusters\":[[[1.0,2.0,0.0],[2.0,1.0,2.0]],[[3.0,0.0,9.0]]]}"
```
"""
function call(message)
  records = message["records"] |> collect
  cluster_num = message["cluster_num"]
  kmeans_result = KMeans.call(records, cluster_num)

  Dict(
    :clusters => kmeans_result[:clusters],
    :centers => kmeans_result[:centers]
  ) |> JSON.json
end
end

ext: jl


text:
using GAP, Hecke, Printf, ArgParse

include(joinpath(Hecke.pkgdir,"examples/FieldEnumeration/FieldEnumeration.jl"))

###############################################################################
#
#  Read-Write
#
###############################################################################

function ArgParse.parse_item(::Type{fmpz}, x::AbstractString)
  if in('^', x)
    l = split(x, '^')
    if length(l) != 2
      throw(error("Could not parse $x as fmpz"))
    end
    l = string.(l)
    return (parse(fmpz, l[1]))^parse(Int, l[2])
  else
    return parse(fmpz, string(x))
  end
end

function ArgParse.parse_item(::Type{Vector{Int}}, x::AbstractString)
  return map(y -> parse(Int, y), split(x, ','))
end

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--input"
      help = "Input file"
      arg_type = String
      default = ""
    "--only-real"
      help = "Only totally real fields"
      action = :store_true
    "--disc-bound"
      help = "Discriminant bound"
      arg_type = fmpz
    "--rootdisc-bound"
      help = "Root discriminant bound"
      arg_type = Float64
    "--only-complex"
      help = "Only totally complex fields"
      action = :store_true
    "--simplify"
      help = "Simplify the field"
      action = :store_true
    "--output"
      help = "Output file"
      arg_type = String
      default = ""
    "--type"
      help = "Type of the abelian extension"
      arg_type = Vector{Int}
      required = true
    "--silent"
      help = "No verbose out"
      action = :store_true
  end

  return parse_args(s)
end

function main()
  parsed_args = parse_commandline()

  local dbound::Union{fmpz, Nothing}
  local rdbound::Union{Nothing, Float64}
  local only_real::Bool
  local only_complex::Bool
  local input::String
  local output::String
  local simplify::Bool
  local type::Vector{Int}
  local silent::Bool

  for (arg, val) in parsed_args
    if arg == "disc-bound"
      dbound = val
    elseif arg == "rootdisc-bound"
      rdbound = val
    elseif arg == "only-real"
      only_real = val
    elseif arg == "only-complex"
      only_complex = val
    elseif arg == "simplify"
      simplify = val
    elseif arg == "input"
      input = val
    elseif arg == "output"
      output = val
    elseif arg == "type"
      type = val
    elseif arg == "silent"
      silent = val
    end
  end

  if !silent
    if dbound !== nothing
      @info "Discriminant bound $dbound"
    end

    if rdbound !== nothing
      @info "Root discriminant bound $rdbound"
    end

    @info "Extensions of type: $type"

    @info "Simplifying fields: $simplify"

    old_output = output

    if output == ""
      output = input * ".out"
    end

    @info "Only real fields: $only_real"
    @info "Only complex fields: $only_complex"
    @info "Input from: $input"

    if old_output == ""
      @info "Output to: $output (default)"
    else
      @info "Output to: $output"
    end

    set_verbose_level(:AbExt, 1)
  end

  if input == ""
    error("Input file must be specified by --input")
  end

  if dbound == nothing && rdbound == nothing
    error("One of --disc-bound or --rootdisc-bound must be specified")
  elseif dbound != nothing && rdbound != nothing
    error("Only one of --disc-bound or --rootdisc-bound can be set")
  end

  if only_complex && only_real
    error("Only one of --only-real or --only-complex can be set")
  end

  sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

  if output == ""
    output = input * ".out"
  end

  if isfile(output)
    error("Output file $output already exists. You may specify the output file by setting --output")
  end

  ending = String(split(basename(input), ".")[2])

  relative_deg = prod(type)

  if ending == "nfdb"
    DBin = Base.read(input, Hecke.NFDB)

    res = eltype(DBin)[]

    for _K in DBin
      K = Hecke.field(_K, cached = false)
      if dbound isa fmpz
        _dbound = dbound
      else
        @assert rdbound isa Float64
        _dbound = fmpz(BigInt(ceil(BigFloat(rdbound)^(degree(K) * relative_deg))))
      end

      if only_real
        if !is_totally_real(K)
          continue
        end

        l = Hecke.abelian_extensions(K, type, _dbound, ramified_at_inf_plc = (true, InfPlc[]))
      elseif only_complex
        l = Hecke.abelian_extensions(K, type, _dbound, ramified_at_inf_plc = (true, real_places(K)))
      else
        l = Hecke.abelian_extensions(K, type, _dbound)
      end

      for C in l
        if dbound isa fmpz
          @assert absolute_discriminant(C) <= dbound
        elseif rdbound isa Float64
          @assert Float64(absolute_discriminant(C)) <= rdbound^(degree(K) * relative_deg)
        end
        if only_real
          @assert signature(C)[2] == 0
        elseif only_complex
          @assert signature(C)[1] == 0
        end

        L = absolute_simple_field(number_field(C), cached = false, simplify = simplify)[1]
        r = Hecke._create_record(L)
        r[:discriminant] = absolute_discriminant(C)
        push!(res, r)
      end
    end

    DBout = Hecke.NFDB(res)
    t = "Abelian extensions with type $type"
    if dbound isa fmpz
      t = t * ", dbound $dbound"
    else
      t = t * ", rdbound $rdbound"
    end
    t = t * ", of $input"
    Hecke.add_meta!(DBout, :title => t)
    open(output, "w") do f
      Base.write(f, DBout)
    end
  end
end

main()


text:
# averaging functions, to operate on resistances
# cond_* functions convert to conductance, calulate mean, then convert to resistance
# res_* functions directly operate on the resistances
cond_cardinal_avg(x, y) = 1 / ((1/x + 1/y) / 2)
cond_diagonal_avg(x, y) = 1 / ((1/x + 1/y) / (2 * √2))
res_cardinal_avg(x, y) = (x + y) / 2
res_diagonal_avg(x, y) = ((x + y) * √2) / 2

function get_lat_lon_dims(A::GeoData.GeoArray)
    dim_names = String.(name(dims(A)))
    lat_first = findall(dim_names .== "Latitude")[1] < findall(dim_names .== "Longitude")[1]
    first_dim_type = lat_first ? Lat : Lon
    second_dim_type = lat_first ? Lon : Lat

    (dims(A, first_dim_type), dims(A, second_dim_type), Band(1:1, Categorical(order = Ordered())))
end

max_stars_repo_name: Circuitscape/SpatialGraphs-OLD

max_issues_repo_name: Circuitscape/SpatialGraphs-OLD

max_forks_repo_name: Circuitscape/SpatialGraphs-OLD


text:
struct BetterFeaturization <: SeaPearl.AbstractFeaturization end

function SeaPearl.featurize(sr::SeaPearl.DefaultStateRepresentation{BetterFeaturization})
    g = sr.cplayergraph
    features = zeros(Float32, nqueens_generator.board_size+3, nv(g))
    for i in 1:nv(g)
        cp_vertex = SeaPearl.cpVertexFromIndex(g, i)
        if isa(cp_vertex, SeaPearl.VariableVertex)
            if isa(cp_vertex.variable, SeaPearl.IntVarViewOffset)
                features[1, i] =0
            else
                id = cp_vertex.variable.id
                features[1, i] = string_to_queen(id)/nqueens_generator.board_size
            end
        end
        if isa(cp_vertex, SeaPearl.ConstraintVertex)
            features[2, i] = 1
        end
        if isa(cp_vertex, SeaPearl.ValueVertex)
            features[3, i] = 1.
            value = cp_vertex.value
            features[3+value, i] = 1.
        end
    end
    features
end

function string_to_queen(id::String)::Int
    parse(Int, split(id,"_")[end])
end

#initializing phase
#update_with_cp_model
#

function SeaPearl.feature_length(::Type{SeaPearl.DefaultStateRepresentation{BetterFeaturization, TS}}) where TS
    return 3 + gen.board_size
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using CairoMakie
using CSV

f = CSV.File("HopfSampling.csv")

random_sampling_trace = f.random_sampling_trace .+ eps()
optimal_sampling_trace = f.optimal_sampling_trace .+ eps()
samples = f.samples

noto_sans = "../resources/NotoSans-Regular.ttf"

tickfontsize    = 28
labelfontsize   = 32

f = Figure(font=noto_sans, resolution=(800,600))

ax = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Sample", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
            yscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10))

l1 = lines!(ax, samples, optimal_sampling_trace, color=:black)
l2 = lines!(ax, samples, random_sampling_trace, color=:black, linestyle=:dash)

Legend(f[1,2], [l1,l2], ["Our Method", "Random"])

ylims!(ax,(1e-3,1))
xlims!(ax,(0,samples[end]))

save("HopfSampling.pdf",f)


text:
using .PDMats: PDMat

export kernelpdmat

"""
    kernelpdmat(k::Kernel, X::AbstractVector)

Compute a positive-definite matrix in the form of a `PDMat` matrix (see [PDMats.jl](https://github.com/JuliaStats/PDMats.jl)),
with the Cholesky decomposition precomputed.
The algorithm adds a diagonal "nugget" term to the kernel matrix which is increased until positive
definiteness is achieved. The algorithm gives up with an error if the nugget becomes larger than 1% of the largest value in the kernel matrix.
"""
function kernelpdmat(κ::Kernel, X::AbstractVector)
    K = kernelmatrix(κ, X)
    Kmax = maximum(K)
    α = eps(eltype(K))
    while !isposdef(K + α * I) && α < 0.01 * Kmax
        α *= 2.0
    end
    if α >= 0.01 * Kmax
        error(
            "Adding noise on the diagonal was not sufficient to build a positive-definite" *
            " matrix:\n\t- Check that your kernel parameters are not extreme\n\t- Check" *
            " that your data is sufficiently sparse\n\t- Maybe use a different kernel",
        )
    end
    return PDMat(K + α * I)
end

"""
    kernelpdmat(k::Kernel, X::AbstractMatrix; obsdim)

If `obsdim=1`, equivalent to `kernelpdmat(k, RowVecs(X))`.
If `obsdim=2`, equivalent to `kernelpdmat(k, ColVecs(X))`.

See also: [`ColVecs`](@ref), [`RowVecs`](@ref)
"""
function kernelpdmat(κ::Kernel, X::AbstractMatrix; obsdim::Union{Int,Nothing}=defaultobs)
    return kernelpdmat(κ, vec_of_vecs(X; obsdim=obsdim))
end


text:
__precompile__(true)

module VectorField

import Plots: quiver, quiver!
import LinearAlgebra: adjoint
export meshgrid, vectorfield, vectorfield!

function meshgrid(x, y)
    x, y = float.(x), float.(y)
    [repeat(x, inner=length(y))'; repeat(y, outer=length(x))']
end

function vectorfield(points, field::Function; arrowscale=0.1)
    vectors = similar(points)
    d = size(points, 2)
    for i in 1:d
        vectors[:, i] .= collect(field(points[:, i]...))
    end
    vectors .*= arrowscale
    quiver(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end
# field is a Function that takes d inputs. d = dimension = 2 (3d not supported by Plots - use Makie.jl)
# points is an Array of size 2xd, each column represents a point in space
# vectorfield2d(field::Function, points <: AbstractArray{<:Number,2}, arrowscale::Real=0.1)

function vectorfield!(points, field::Function; arrowscale=0.1)
    vectors = similar(points)
    d = size(points, 2)
    for i in 1:d
        vectors[:, i] .= collect(field(points[:, i]...))
    end
    vectors .*= arrowscale
    quiver!(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end
# identical, but plots over existing items.
# to keep in line with Plots formatting

function vectorfield(xs, ys, field::Function; arrowscale=0.1)
    points = meshgrid(xs, ys)
    vectors = similar(points)
    d = size(points, 2)
    for i in 1:d
        vectors[:, i] .= collect(field(points[:, i]...))
    end
    vectors .*= arrowscale
    quiver(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end
# input x and y ranges, grid is interpolated

function vectorfield!(xs, ys, field::Function; arrowscale=0.1)
    points = meshgrid(xs, ys)
    vectors = similar(points)
    d = size(points, 2)
    for i in 1:d
        vectors[:, i] .= collect(field(points[:, i]...))
    end
    vectors .*= arrowscale
    quiver!(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end

function vectorfield(xs, ys, u::T, v::T; arrowscale=0.1) where {T <: AbstractArray{<:Number,N} where N}
    points = meshgrid(xs, ys)
    d = size(points, 2)
    d == length(u) == length(v) || error("Dimensions do not match. Ensure length(u)==length(v)==length(xs)*length(ys)")
    vectors = [reshape(u', 1, :); reshape(v', 1, :)]
    vectors .*= arrowscale
    quiver(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end
# u, v are 2 matrices, the value of u[i, j] is the x component of the arrow,
# the value of v[i, j] is the y component of the arrow
# Grid will be interpolated as:
#     ys - >           ys - >
# xs                xs
# |    u            |    v
# v           ,     v

function vectorfield!(xs, ys, u::T, v::T; arrowscale=0.1) where {T <: AbstractArray{<:Number,N} where N}
    points = meshgrid(xs, ys)
    d = size(points, 2)
    d == length(u) == length(v) || error("Dimensions do not match. Ensure length(u)==length(v)==length(xs)*length(ys)")
    vectors = [reshape(u', 1, :); reshape(v', 1, :)]
    vectors .*= arrowscale
    quiver!(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end

function vectorfield(xs, ys, uv::T; arrowscale=0.1) where {T <: AbstractArray{S,N} where S where N}
    points = meshgrid(xs, ys)
    d = size(points, 2)
    u = [uv[i][1] for i in CartesianIndices(uv)]
    v = [uv[i][2] for i in CartesianIndices(uv)]
    d == length(u) == length(v) || error("Dimensions do not match. Ensure length(u)==length(v)==length(xs)*length(ys)")
    vectors = [reshape(u', 1, :); reshape(v', 1, :)]
    vectors .*= arrowscale
    quiver(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end
# uv is a matrix of 2-element arrays [x, y]. Arrows will be interpreted componentwise.
# Grid will be interpolated as:
#     ys - >
# xs
# |    uv
# v

function vectorfield!(xs, ys, uv::T; arrowscale=0.1) where {T <: AbstractArray{S,N} where S where N}
    points = meshgrid(xs, ys)
    d = size(points, 2)
    u = [uv[i][1] for i in CartesianIndices(uv)]
    v = [uv[i][2] for i in CartesianIndices(uv)]
    d == length(u) == length(v) || error("Dimensions do not match. Ensure length(u)==length(v)==length(xs)*length(ys)")
    vectors = [reshape(u', 1, :); reshape(v', 1, :)]
    vectors .*= arrowscale
    quiver!(points[1, :], points[2, :], quiver=(vectors[1, :], vectors[2, :]))
end

end


text:
@testset "linear algebra" begin
    @test issquare(reshape(rand(Double32,3*3), 3, 3))

    # issue #77
    M = Double64.([1 0; 1 1])
    @test_throws ErrorException exp(M)
    M = Complex{Double64}.([1 0; 1 1])
    @test_throws ErrorException exp(M)
end


text:
module NFFT

using Base.Cartesian
using FFTW
using Distributed
using SparseArrays
using LinearAlgebra
using CUDA
using Graphics: @mustimplement

export AbstractNFFTPlan, plan_nfft, nfft, nfft_adjoint, ndft, ndft_adjoint
import Base.size

@enum PrecomputeFlags begin
  LUT = 1
  FULL = 2
end

@enum Device begin
  CPU = 1
  CUDAGPU = 2
end

#########################
# abstract NFFT interface
#########################
include("Abstract.jl")

#################################
# currently implemented NFFTPlans
#################################
include("windowFunctions.jl")
include("precomputation.jl")

include("CpuNFFT.jl")
if CUDA.functional()
    using CUDA.CUSPARSE
    include("CuNFFT.jl")
end

#################
# factory methods
#################
"""
        plan_nfft(x::Matrix{T}, N::NTuple{D,Int}, rest...; device::Device=CPU, kargs...)

compute a plan for the NFFT of a size-`N` array at the nodes contained in `x`.

The computing device (CPU or GPU) can be set using the keyworkd argument `device`
to NFFT.CPU or NFFT.CUDAGPU
"""
function plan_nfft(x::Matrix{T}, N::NTuple{D,Int}, rest...; device::Device=CPU, kargs...) where {T,D}
    if device==CPU
        p = NFFTPlan(x, N, rest...; kargs...)
    elseif device==CUDAGPU
        @assert CUDA.functional()==true "device $(dev) requires a functional CUDA setup"
        p = CuNFFTPlan(x, N, rest...; kargs...)
    else
        error("device $(dev) is not yet supported by NFFT.jl")
    end
    return p
end

"""
        plan_nfft(x::Vector{T}, dim::Integer, N::NTuple{D,Int}, rest...; device::Device=CPU, kargs...)

compute a plan for the NFFT of a size-`N` array at the nodes contained in `x`
and along the direction `dim`.
The computing device (CPU or GPU) can be set using the keyworkd argument `device`.
Currently only NFFT.CPU is supported.
"""
function plan_nfft(x::Vector{T}, dim::Integer, N::NTuple{D,Int}, rest...; device::Device=CPU, kargs...) where {T,D}
    if device==CPU
        p = NFFTPlan(x, dim, N, rest...; kargs...)
    elseif device==CUDAGPU
        error("device $(dev) doest not yet support directional NFFTs.")
    else
        error("device $(dev) is not yet supported by NFFT.jl")
    end
    return p
end

function plan_nfft(x::Matrix{T}, dim::Integer, N::NTuple{D,Int}, rest...; device::Device=CPU, kargs...) where {T,D}
    if size(x,1) != 1 && size(x,2) != 1
        throw(DimensionMismatch())
    end
    if device==CPU
        p = NFFTPlan(vec(x), dim, N, rest...; kargs...)
    elseif device==CUDAGPU
        error("device $(dev) doest not yet support directional NFFTs.")
    else
        error("device $(dev) is not yet supported by NFFT.jl")
    end
    return p
end

plan_nfft(x::AbstractMatrix{T}, N::NTuple{D,Int}, rest...; kwargs...) where {D,T} =
    plan_nfft(collect(x), N, rest...; kwargs...)

plan_nfft(x::AbstractVector, N::Integer, rest...; kwargs...) =
    plan_nfft(reshape(x,1,length(x)), (N,), rest...; kwargs...)

#########################
# high-level NFFT methods
#########################
"""
        nfft(x, f::AbstractArray{T,D}, rest...; device::Device=CPU, kwargs...)

calculates the NFFT of the array `f` for the nodes contained in the matrix `x`
The output is a vector of length M=`size(nodes,2)`
"""
function nfft(x, f::AbstractArray{T,D}, rest...; device::Device=CPU, kwargs...) where {T,D}
    p = plan_nfft(x, size(f), rest...; device=device, kwargs... )
    return nfft(p, f)
end

"""
        nfft_adjoint(x, N, fHat::AbstractArray{T,D}, rest...; device::Device=CPU, kwargs...)

calculates the adjoint NFFT of the vector `fHat` for the nodes contained in the matrix `x`.
The output is an array of size `N`
"""
function nfft_adjoint(x, N, fHat::AbstractVector{T}, rest...; device::Device=CPU, kwargs...) where T
    p = plan_nfft(x, N, rest...; device=device, kwargs...)
    return nfft_adjoint(p, fHat)
end

include("directional.jl")
include("multidimensional.jl")
include("samplingDensity.jl")
include("NDFT.jl")

end


text:
include("generate_side_nav.jl")
@info "Generated side nav"
include("generate_toc_pages.jl")
@info "Generated TOC pages"
include("generate_papers.jl")
@info "Generated 'Papers' page"


text:
"""Julia program to implement Comb Sort
Comb sort is an improvised Bubble sort. Bubble sort removes its inversions one by one whereas In Comb sort,
more than one inversion is removed by one swap and hence it is more efficient than Bubble sort.
"""

function comb_sort(arr, n)
    flag = 1
    space = n
    while (space > 1 || flag == 1)
        # Here empirically found shrink factor 1.3 is used.
        space = (space * 10) ÷ 13

        # Assigning space its minimum value of 1
        if (space < 1)
            space = 1
        end
        flag = 0

        # Compare elements with the space value
        for i in 1:(n - space)
            if (arr[i] > arr[i + space])
                arr[i], arr[i + space] = arr[i + space], arr[i]
                flag = 1
            end
        end
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int64, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
comb_sort(arr, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers do you want to sort? 6
Enter the numbers: 123 13 42 324 153 43
The numbers in sorted order is [13, 42, 43, 123, 153, 324]
"""


text:
using Documenter
using QuantumCumulants

ENV["GKSwstype"] = "100" # enable headless mode for GR to suppress warnings when plotting

pages = [
        "index.md",
        "theory.md",
        "tutorial.md",
        "correlation.md",
        "implementation.md",
        "api.md",
        "Examples" => [
            "examples/single-atom-laser-spectrum.md"
            "examples/mollow.md"
            "examples/many-atom-laser.md"
            "examples/optomechanical-cooling.md"
            "examples/superradiant-laser.md"
            "examples/excitation-transport-chain.md"
            ]
    ]

makedocs(
    sitename = "QuantumCumulants.jl",
    modules = [QuantumCumulants],
    pages = pages,
    checkdocs=:exports,
    format = Documenter.HTML(
                            mathengine=MathJax(),
                            footer="[**Back to GitHub**](https://github.com/qojulia/QuantumCumulants.jl)"
                            )
    )

deploydocs(
    repo = "github.com/qojulia/QuantumCumulants.jl.git",
    push_preview = false,
    )


text:
module Record
using MLStyle
using MLStyle.MatchCore
using MLStyle.MatchImpl
using MLStyle.AbstractPatterns
using MLStyle.AbstractPatterns.BasicPatterns
using MLStyle.Qualification

export @as_record, record_def

@nospecialize
function P_partial_struct_decons(t, partial_fields, ps, prepr::AbstractString="$t")
    function tcons(_...)
        t
    end

    comp = PComp(
        prepr, tcons;
    )
    function extract(sub::Any, i::Int, ::Any, ::Any)
        :($sub.$(partial_fields[i]))
    end
    decons(comp, extract, ps)
end

function record_def(@nospecialize(Struct), line::LineNumberNode, ::Module)
    quote
        $line
        function $MLStyle.pattern_uncall(t::Type{$Struct}, self::Function, type_params, type_args, args)
            $line
            isempty(type_params) || return begin
                call = Expr(:call, t, args...)
                ann = Expr(:curly, t, type_args...)
                self(Where(call, ann, type_params))
            end
            all_field_names = fieldnames(t)
            partial_field_names = Symbol[]
            patterns = Function[]
            $MLStyle.@switch args begin
            @case [Expr(:parameters, kwargs...), args...]
            @case let kwargs = [] end
            end
            n_args = length(args)
            if all(Meta.isexpr(arg, :kw) for arg in args)
                for arg in args
                    field_name = arg.args[1]
                    field_name in all_field_names || error("$t has no field $field_name.")
                    push!(partial_field_names, field_name)
                    push!(patterns, self(arg.args[2]))
                end
            elseif length(all_field_names) === n_args
                append!(patterns, map(self, args))
                append!(partial_field_names, all_field_names)
            elseif n_args === 1 && args[1] === :_
            elseif n_args !== 0
                error("count of positional fields should be 0 or the same as the fields($all_field_names)")
            end
            for e in kwargs
                $MLStyle.@switch e begin
                @case :: Symbol
                    e in all_field_names || error("unknown field name $e for $t when field punnning.")
                    push!(partial_field_names, e)
                    push!(patterns, $P_capture(e))
                    continue
                @case Expr(:kw, key::Symbol, value)
                    key in all_field_names || error("unknown field name $key for $t when field punnning.")
                    push!(partial_field_names, key)
                    push!(patterns, $and([$P_capture(key), self(value)]))
                    continue
                @case _

                    error("unknown sub-pattern $e in " * string(t) * ".")
                end
            end

            ret = $P_partial_struct_decons(t, partial_field_names, patterns)
            isempty(type_args) && return ret
            $and([self(Expr(:(::), Expr(:curly, t, type_args...))) , ret])
        end
    end
end

function as_record(@nospecialize(n), line::LineNumberNode, __module__::Module)
    @switch n begin
    @case ::Symbol
        return record_def(n, line, __module__)
    @case :(struct $hd{$(_...)}
                $(_...)
            end) ||
          :(struct $hd{$(_...)} <: $_
            $(_...)
            end) ||
          :(struct $hd <: $_
                $(_...)
            end) ||
          :(struct $hd
                $(_...)
            end)
        return Expr(
            :block,
            n,
            record_def(hd, line, __module__)
        )
    @case _
        error("malformed structure $n")
    end
end

macro as_record(qualifier, n)
    deprecate_qualifier_macro(qualifier, __source__)
    esc(as_record(n, __source__, __module__))
end

macro as_record(n)
    esc(as_record(n, __source__, __module__))
end
@specialize
end


text:
tests =
  ["test_empmctuner",
  "test_ARS1",
  "test_ARS2",
  "test_slice_sampler"]

println("Running tests:")

for t in tests
  tfile = t*".jl"
  println("  * $(tfile) *")
  include(tfile)
end


text:
include("/home/bendavid/rRNA_16s_Analysis/OligoCompressor.jl")

using FASTX
#using Plots

function edit_distance(seq1::LongSequence{DNAAlphabet{4}},seq2::LongSequence{DNAAlphabet{4}})
    m=length(seq1);
    n=length(seq2);
    d=zeros(m+1,n+1)
    for i = 1:m
        d[i+1,1]=i
    end
    for j = 1:n
        d[1,j+1]=j
    end
    for j = 1:n
        for i=1:m
            substitution_cost=1
            if iscompatible(seq1[i],seq2[j])
                substitution_cost=0;
            end
            d[i+1,j+1]= min(d[i,j+1]+1,d[i+1,j]+1,d[i,j]+substitution_cost)
        end
    end
    return d[m+1,n+1]
end

function hamming_distance(seq1::LongSequence{DNAAlphabet{4}},seq2::LongSequence{DNAAlphabet{4}})
    dist_counter=0
    if  length(seq1) != length(seq2)
        error("Sequences must be the same length.")
    end
    n=length(seq1)
    for i=1:n
        if !iscompatible(seq1[i],seq2[i])
            dist_counter+=1
        end
    end
    return dist_counter
end

function hamming_distance_slider(longseq::LongSequence{DNAAlphabet{4}},shortseq::LongSequence{DNAAlphabet{4}})
    sliderlength=length(longseq)-length(shortseq)+1
    dists=zeros(Int64,sliderlength)
    for i=1:sliderlength
        dists[i]=hamming_distance(longseq[i:i+length(shortseq)-1],shortseq)
    end
    return dists
end

seqs=String[]
reader= FASTA.Reader(open("/home/bendavid/rRNA_16s_Analysis/rRNA_16s.fasta","r"))

for record in reader
    push!(seqs,sequence(record))
end
close(reader)

seqs=LongDNASeq.(seqs)
logics=trues(length(seqs))
for i=1:length(seqs)
    if sum(isambiguous.(seqs[i]))>0
        logics[i]=false
    end
end
seqs=seqs[logics]

primer=dna"AGAGTTTGATCCTGGCTCAG"

k=collect(0:length(primer))
n=length(seqs)
counts=zeros(length(k))
idxs=zeros(Int64,n,length(k))
for i=1:length(k)
    for j=1:n
        idx=approxsearchindex(seqs[j],primer,i)
        if idx!=0
            counts[i]+=1
            idxs[j,i]=idx
        end
    end
end

pcts=counts./n
#plot(k,pcts)
#histogram(idxs[:,5])
#histogram(idxs[:,6])
#histogram(idxs[:,7])

mismatches=4
valids=idxs[:,mismatches+1].>0 .& idxs[:,mismatches+1].<2
valid_seqs=seqs[valids]

seqs_100=[dna"-" for i=1:n]

for i=1:n
    seqs_100[i]=seqs[i][1:100]
end

valid_primers=LongDNASeq[]
for i=1:n
    dists=hamming_distance_slider(seqs_100[i],primer)
    val,idx=findmin(dists)
    if val <=4
        push!(valid_primers,seqs_100[i][idx:idx+length(primer)-1])
    end
end

valid_primers
unique(valid_primers)

compressed_primers=oligo_pool_compressor(unique(valid_primers))

DF1=DataFrame(sequence=unique(valid_primers))
DF2=DataFrame(sequence=compressed_primers)

CSV.write("./rRNA16s_uncompressed.csv",DF1)
CSV.write("./rRNA16s_compressed.csv",DF2)


text:
@testset "FVP" begin
    # test for correct output and no side effects
    d1, d2 = 5, 10
    @test let glls = rand(d1, d2), x = rand(d1)
        glls2 = copy(glls)
        x2 = copy(x)
        y1 = LyceumAI.FVP(glls2) * x2
        y2 = glls * transpose(glls) * x
        isapprox(y1, y2) && eltype(y1) === eltype(glls) && x == x2 && glls == glls2
    end
    @test let glls = rand(d1, d2), x = rand(Float32, d1)
        glls2 = copy(glls)
        x2 = copy(x)
        y1 = LyceumAI.FVP{Float32}(glls2) * x2
        y2 = glls * transpose(glls) * x
        isapprox(y1, y2) && eltype(y1) === Float32 && x == x2 && glls == glls2
    end
    @test let glls = rand(d1, d2), x = rand(d1), y1 = zero(x)
        glls2 = copy(glls)
        x2 = copy(x)
        mul!(y1, LyceumAI.FVP(glls2), x2)
        y2 = glls * transpose(glls) * x
        isapprox(y1, y2) && x == x2 && glls == glls2
    end
    @test let glls = rand(d1, d2), x = rand(d1), y1 = zero(x)
        glls2 = copy(glls)
        x2 = copy(x)
        mul!(y1, LyceumAI.FVP(glls2, true), x2)
        y2 = 1/d2 * glls * transpose(glls) * x
        isapprox(y1, y2) && x == x2 && glls == glls2
    end
end


text:
using BoardGamesStrategies
using Test

include("TicTacToe.jl")

b = playgame(TicTacToeBoard(), MCTSStrategy{TicTacToe}(1000, √2), MCTSStrategy{TicTacToe}(1000, √2));

b = TicTacToeBoard()
b = play(b, 1)
b = play(b, 2)

function heuristic(b::TicTacToeBoard, player::Int64)
    b.winner * -((player - 2) * 2 + 1)
end

@testset "BoardGamesStrategies.jl" begin
    #a bit random but player 1 should win
    @test 0 <= winner(playgame(b, RandomStrategy{TicTacToe}(),
                                  RandomStrategy{TicTacToe}())) <= 2

    @test winner(playgame(b, MCTSStrategy{TicTacToe}(10000, √2),
                             MCTSStrategy{TicTacToe}(10000, √2))) == 1

    @test winner(playgame(b, Minimax{TicTacToe}(10,heuristic),
                             Minimax{TicTacToe}(10,heuristic))) == 1

    @test winner(playgame(b, MinimaxAlphaBeta{TicTacToe}(10, heuristic),
                             MinimaxAlphaBeta{TicTacToe}(10, heuristic))) == 1

end


text:
function metric_callback(m::AbstractMetricParams{T}) where {T}
    min_r = inner_radius(m)
    DiscreteCallback((u, λ, integrator) -> u[6] ≤ min_r * 1.1 || u[6] > 1200.0, terminate!)
end

function create_callback_set(m::AbstractMetricParams{T}, cb::Nothing) where {T}
    mcb = metric_callback(m)
    mcb isa Tuple ? CallbackSet(mcb...) : mcb
end

function create_callback_set(
    m::AbstractMetricParams{T},
    cb::NTuple{N,SciMLBase.DECallback}
) where {T,N}
    mcb = metric_callback(m)
    mcb isa Tuple ? CallbackSet(mcb..., cb...) : CallbackSet(mcb, cb...)
end

function create_callback_set(m::AbstractMetricParams{T}, cb::SciMLBase.DECallback) where {T}
    mcb = metric_callback(m)
    mcb isa Tuple ? CallbackSet(mcb..., cb) : CallbackSet(mcb, cb)
end


text:
#! /usr/bin/julia

# Rosetta Code, SHA-1

using Nettle

function sha1sum(s::String)
    bytes2hex(sha1_hash(s))
end

println("Testing SHA-1 function against FIPS 180-1")

s = sha1sum("abc")
t = lowercase("A9993E364706816ABA3E25717850C26C9CD0D89D")
print("    \"abc\" should yield ", t)
if s == t
    println(", and it does.")
else
    println(", but it yields ", s)
end

s = sha1sum("abcdbcdecdefdefgefghfghighijhijkijkljklmklmnlmnomnopnopq")
t = lowercase("84983E441C3BD26EBAAE4AA1F95129E5E54670F1")
print("    \"abcdbc...\" should yield ", t)
if s == t
    println(", and it does.")
else
    println(", but it yields ", s)
end

s = sha1sum("a"^10^6)
t = lowercase("34AA973CD4C4DAA4F61EEB2BDBAD27316534016F")
print("    a million \"a\"s should yield ", t)
if s == t
    println(", and it does.")
else
    println(", but it yields ", s)
end

println("\nAlso")
s = "Rosetta Code"
h = sha1sum(s)
println("    ", s, " => ", h)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
function writegexf(f::IO, g::AGraphOrDiGraph)
    xdoc = XMLDocument()
    xroot = setroot!(xdoc, ElementNode("gexf"))
    xroot["xmlns"] = "http://www.gexf.net/1.2draft"
    xroot["version"] = "1.2"
    xroot["xmlns:xsi"] = "http://www.w3.org/2001/XMLSchema-instance"
    xroot["xsi:schemaLocation"] = "http://www.gexf.net/1.2draft/gexf.xsd"

    xmeta = addelement!(xroot, "meta")
    xg = addelement!(xroot, "graph")
    xg["defaultedgetype"] = is_directed(g) ? "directed" : "undirected"
    xnodes = addelement!(xg, "nodes")
    for i in 1:nv(g)
        xv = addelement!(xnodes, "node")
        xv["id"] = "$(i-1)"
    end

    xedges = addelement!(xg, "edges")
    m = 0
    for e in edges(g)
        xe = addelement!(xedges, "edge")
        xe["id"] = "$m"
        xe["source"] = "$(src(e)-1)"
        xe["target"] = "$(dst(e)-1)"
        m += 1
    end

    prettyprint(f, xdoc)
    return 1
end

function gexf_read_one_graph!(el::EzXML.Node, ::Type{G}) where G
    elnodes = getchild(el, "nodes")
    nodes = Dict{String,Int}()
    for (i,f) in enumerate(eachelement(elnodes))
        nodes[f["id"]] = i
    end

    g = G(length(nodes))

    eledges = getchild(el, "edges")
    for f in eachelement(eledges)
        n1 = f["source"]
        n2 = f["target"]
        add_edge!(g, nodes[n1], nodes[n2])
    end

    return g
end

function readgexf(io::IO, ::Type{G}) where G<:AGraphOrDiGraph
    xdoc = parsexml(read(io, String))
    xroot = root(xdoc)  # an instance of XMLElement
    nodename(xroot) == "gexf" || error("Not a Gexf file")
    xg = getchild(xroot, "graph")
    isdir = false
    if haskey(xg, "defaultedgetype")
        isdir = xg["defaultedgetype"] == "directed"  ? true  : false
    end
    H = isdir ? digraphtype(G) : graphtype(G)
    return gexf_read_one_graph!(xg, H)
end

function readnetgexf(io::IO, ::Type{G}) where G<:AGraphOrDiGraph
    xdoc = parsexml(read(io, String))
    xroot = root(xdoc)  # an instance of XMLElement
    nodename(xroot) == "gexf" || error("Not a Gexf file")
    xg = getchild(xroot, "graph")
    isdir = false
    if haskey(xg, "defaultedgetype")
        isdir = xg["defaultedgetype"] == "directed"  ? true  : false
    end

    gpropkeys=Dict{String, Tuple{String, DataType}}();
    vpropkeys=Dict{String, Tuple{String, DataType}}();
    epropkeys=Dict{String, Tuple{String, DataType}}();

    # Reading attributes
    # Note: if no default value is specified for an attribute
    # each node/edge has to implement it
    for xattr in elements(xg)
        nodename(xattr) != "attributes" && continue
        pk =  xattr["class"] == "node" ? vpropkeys :
              xattr["class"] == "edge" ? epropkeys : error("attr")

        for el in elements(xattr)
            pk[el["id"]] = (el["title"],  gexf_types_rev[el["type"]])
        end
    end

    H = isdir ? digraphtype(G) : graphtype(G)
    return gexf_read_one_net!(xg, H, gpropkeys, vpropkeys, epropkeys)
end

const gexf_types = Dict( Int32 => "integer",
                            Int64 => "integer",
                            Bool => "boolean",
                            Float32 => "float",
                            Float64 => "double",
                            String => "string",
                            )

const gexf_types_rev = Dict("integer"  =>  Int,
                            "boolean" =>  Bool,
                            "float"   =>  Float32,
                            "double"  =>  Float64,
                            "string"  =>  String,
                        )

gexfstring(x) = string(x)
gexfstring(v::Vector) = join((Printf.@sprintf("%.10g",x) for x in v), ", ")

gexfparse(T, x::String) = parse(T, x)
gexfparse(::Type{String}, x::String) = x

function gexf_read_one_net!(xg::EzXML.Node, ::Type{G},
                        gpropkeys, vpropkeys, epropkeys) where G
    nodes = Dict{String,Int}()
    nodeid = 1
    # traverse the tree to map id to 1:n
    for el in eachelement(getchild(xg,"nodes"))
        nodename(el) != "node" && continue #TODO can be removed
        nodes[el["id"]] = nodeid
        nodeid += 1
    end

    g = G(length(nodes))
    for (pname,T) in values(vpropkeys); vprop!(g, pname, T); end
    for (pname,T) in values(epropkeys); eprop!(g, pname, T); end

    for f in eachelement(getchild(xg,"nodes"))
        @assert nodename(f) == "node"
        v = nodes[f["id"]]
        for xattr in eachattribute(f)
            nattr = nodename(xattr)
            nattr == "id" && continue
            !has_vprop(g, nattr) && vprop!(g, nattr, String)
            vprop(g, nattr)[v] = f[nattr]
        end
        for el in eachelement(f)
            if nodename(el) == "attvalues"
                for xattr in eachelement(el)
                    pname, T = vpropkeys[xattr["for"]]
                    vprop(g, pname)[v] = gexfparse(T, xattr["value"])
                end
            else
                pname = nodename(el)
                if endswith(namespace(el), "viz")
                    pname = "viz:"*pname
                end
                attr = nodename.(attributes(el))
                if "value" in attr && length(attr) == 1
                    !has_vprop(g, pname) && vprop!(g, pname, String)
                    vprop(g, pname)[v] = el["value"]
                else
                    !has_vprop(g, pname) && vprop!(g, pname, Dict{String,String})
                    vprop(g, pname)[v] = Dict{String,String}(a=>el[a] for a in attr)
                end
            end
        end
    end

    for f in eachelement(getchild(xg,"edges"))
        @assert nodename(f) == "edge"
        u = nodes[f["source"]]
        v = nodes[f["target"]]
        ok, e = add_edge!(g, u, v)
        for xattr in eachattribute(f)
            nattr = nodename(xattr)
            nattr ∈ ["id","source","target"] && continue
            Tattr = nattr == "weight" ? Float64 : String
            !has_eprop(g, nattr) && eprop!(g, nattr, Tattr)
            eprop(g, nattr)[e] = Tattr == String ? f[nattr] : parse(Tattr, f[nattr])
        end
        for el in eachelement(f)
            if nodename(el) == "attvalues"
                for xattr in eachelement(el)
                    pname, T = epropkeys[xattr["for"]]
                    eprop(g, pname)[e] = gexfparse(T, xattr["value"])
                end
            else
                pname = nodename(el)
                attr = nodename.(attributes(el))
                if "value" in attr && length(attr) == 1
                    !has_eprop(g, pname) && eprop!(g, pname, String)
                    eprop(g, pname)[e] = el["value"]
                else
                    !has_eprop(g, pname) && eprop!(g, pname, Dict{String,String})
                    eprop(g, pname)[e] = Dict{String,String}(a=>el[a] for a in attr)
                end
            end
        end
        # for el in eachelement(f)
        #     pname, T = epropkeys[el["key"]]
        #     m = eprop(g, pname)
        #     m[e] = gexfparse(T, nodecontent(el))
        # end
    end
    return g
end

function writenetgexf(f::IO, g::ANetOrDiNet)
    xdoc = XMLDocument()
    xroot = setroot!(xdoc, ElementNode("gexf"))
    xroot["version"]="1.3"
    xroot["xmlns"] = "http://www.gexf.net/1.3"
    xroot["xmlns:viz"]= "http://www.gexf.net/1.3/viz"
    xroot["xmlns:xsi"] = "http://www.w3.org/2001/XMLSchema-instance"
    xroot["xsi:schemaLocation"]="http://www.gexf.net/1.3/gexf.xsd"

    xmeta = addelement!(xroot, "meta")
    xmeta["lastmodifieddate"] = string(Dates.today())
    xg = addelement!(xroot, "graph")
    xg["defaultedgetype"] = is_directed(g) ? "directed" : "undirected"
    xg["mode"] = "static"

    vattrnames = String[]
    if length(vprop(g)) > 0
        xa = addelement!(xg, "attributes")
        xa["class"] = "node"
        xa["mode"] = "static"
        for (pname, m) in vprop(g)
            pname == "label" && continue
            startswith(pname, "viz:") && continue
            T = valtype(m)
            if T ∈ keys(gexf_types)
                xaa = addelement!(xa, "attribute")
                push!(vattrnames, pname)
                xaa["id"] = pname
                xaa["title"] = pname
                xaa["type"] = gexf_types[T]
            end
        end
    end

    eattrnames = String[]
    if length(eprop(g)) > 0
        xa = addelement!(xg, "attributes")
        xa["class"] = "edge"
        xa["mode"] = "static"
        for (pname, m) in eprop(g)
            pname == "weight" && continue
            startswith(pname, "viz:") && continue
            T = valtype(m)
            xaa = addelement!(xa, "attribute")
            push!(eattrnames, pname)
            xaa["id"] = pname
            xaa["title"] = pname
            xaa["type"] = gexf_types[T]
        end
    end

    xnodes = addelement!(xg, "nodes")
    xnodes["count"] = nv(g)
    for i in 1:nv(g)
        xv = addelement!(xnodes, "node")
        xv["id"] = "$(i-1)"

        for (name, val) in vprop(g, i)
            if name == "label"
                xv["label"] = val
            elseif name ∈ vattrnames
                xvals = getchild!(xv, "attvalues")
                x = addelement!(xvals, "attvalue")
                x["for"] = name
                x["value"] = gexfstring(val)
            elseif typeof(val) <: Dict
                xval = addelement!(xv, name)
                for (k, v) in val
                    xval[k] = v
                end
            else
                xval = addelement!(xv, name)
                xval["value"] = gexfstring(val)
            end
        end
    end

    xedges = addelement!(xg, "edges")
    xedges["count"] = ne(g)
    m = 0
    for e in edges(g)
        xe = addelement!(xedges, "edge")
        xe["id"] = "$m"
        xe["source"] = "$(src(e)-1)"
        xe["target"] = "$(dst(e)-1)"
        if has_eprop(g, "weight", e)
            xe["weight"] = eprop(g,"weight")[e]
        end
        for (name, val) in eprop(g, e)
            if name == "weight"
                xe["weight"] = val
            elseif name ∈ eattrnames
                xvals = getchild!(xe, "attvalues")
                x = addelement!(xvals, "attvalue")
                x["for"] = name
                x["value"] = gexfstring(val)
            elseif typeof(val) <: Dict
                xval = addelement!(xe, name)
                for (k, v) in val
                    xval[k] = v
                end
            else
                xval = addelement!(xe, name)
                xval["value"] = gexfstring(val)
            end
        end
        m += 1
    end

    prettyprint(f, xdoc)
    return 1
end

filemap[:gexf] = (readgexf, writegexf, readnetgexf, writenetgexf)


text:
using Test, DiscreteExteriorCalculus
const DEC = DiscreteExteriorCalculus
using LinearAlgebra: norm
using Statistics: mean

@testset "Point and Simplex" begin
    p = Point(1,2,3)
    @test typeof(p) == Point{3}
    s = Simplex(Point(1,2,3), Point(4,5,6), Point(1,4,7))
    @test typeof(s) == Simplex{3, 3}
end

@testset "Barycentric and barycentric_projection" begin
    m = Metric(3)
    s = Simplex(Point(0,0,0), Point(2,0,0), Point(1,.5,0))
    for p in [s.points; Point(3,5,7); Point(-12,25,-100); Point(24.3,-3.5,-2.2)]
        b = Barycentric(m, s, p)
        @test norm(Point(b).coords - [p.coords[1:end-1]; 0]) < 1e-14
    end
end

@testset "distance" begin
    m = Metric(3)
    s = Simplex(Point(0,0,0), Point(2,0,0), Point(1.3,.7,0), Point(0,0,1))
    ps = [s.points; Point(3,5,7); Point(-12,25,-100); Point(24.3,-3.5,-2.2)]
    for p1 in ps
        for p2 in ps
            d1 = sqrt(abs(distance_square(m, Barycentric(m, s, p1), p2)))
            d2 = norm(m, p1.coords - p2.coords)
            d3 = norm(p1.coords - p2.coords) # using Euclidean metric
            @test d2 ≈ d3
            @test abs(d1 - d2) < 1e-5
        end
    end
end

@testset "circumsphere" begin
    function is_circumcenter(m::Metric{N}, s::Simplex{N, K}, p::Point{N}) where {N, K}
        if K == 1
            return s.points[1] == p
        else
            dists = [norm(m, p.coords - q.coords) for q in s.points]
            return all([d ≈ dists[1] for d in dists])
        end
    end
    function circumradius_square(m::Metric{N}, s::Simplex{N}) where N
        center, _ = circumsphere(m, s)
        return mean([norm_square(m, center.coords - v.coords) for v in s.points])
    end
    m = Metric(3)
    s = Simplex(Point(1,2,3), Point(-12,3,5), Point(8,3,4), Point(9,25,-2))
    for k in 1:4
        for sub_s in subsimplices(s, k)
            center, radius_square = circumsphere(m, sub_s)
            @test is_circumcenter(m, sub_s, center)
            @test radius_square ≈ circumradius_square(m, sub_s)
        end
    end
end

@testset "barycentric_subspace" begin
    m = Metric(2)
    s = Simplex(Point(0,0), Point(1,0), Point(0,1))
    f = Simplex(s.points[1:2]...)
    p = Point(.5, 3)
    M, v = DEC.barycentric_subspace(m, s, f, p)
    @test M * Barycentric(m, s, p).coords == v
    @test M * [.5, .5, 0] == v
    @test M * Barycentric(m, s, Point(.5, -2)).coords == v
end

@testset "rotate_about_face" begin
    m = Metric(3)
    s = Simplex(Point(-1,0,0), Point(0,0,0), Point(0,1,0))
    r, θ = 3, π/6
    p = Point(r * cos(θ), 0, r * sin(θ))
    b = DEC.rotate_about_face(m, s, p, 1)
    @test norm(Point(b).coords - [r, 0, 0]) < 1e-15

    m = Metric(2)
    s = Simplex(Point(0,0),Point(1,0))
    p = Point(1,1)
    b = DEC.rotate_about_face(m, s, p, 2)
    @test Point(b).coords ≈ [-sqrt(2), 0]
    b = DEC.rotate_about_face(m, s, p, 1)
    @test Point(b).coords ≈ [2, 0]
end

@testset "pairwise_delaunay" begin
    m = Metric(3)
    s1 = Simplex(Point(-.1,.5,0), Point(0,0,0), Point(0,1,0))
    center, _ = circumsphere(m, s1)
    θ = π/6
    for r in [.1, 1, 5, 10]
        p = Point(r * cos(θ), 0, r * sin(θ))
        proj_p = Point(r, 0, 0)
        # using Euclidean metric
        delaunay = norm(center.coords - proj_p.coords)^2 - norm(center.coords - s1.points[1].coords)^2
        s2 = Simplex(s1.points[2:end]..., p)
        @test pairwise_delaunay(m, s1, s2) ≈ delaunay
    end
end

# TODO test with other metrics


text:
module Manager

using Base.Test: @test
using HttpServer
import Requests
import StaffJoy
import JSON
import Logging

export run_server, unit_test, functional_test

if "ENV" in keys(ENV) && (
    ENV["ENV"] == "stage" || ENV["ENV"] == "prod"
    )
    env = ENV["ENV"]
    Logging.configure(level=Logging.INFO)
else
    Logging.configure(level=Logging.DEBUG)
end

# Boot Script
include("Manager/env.jl")
if "ENV" in keys(ENV)
    env = ENV["ENV"]
    if !(env in keys(env_config))
        Logging.critical("Environment $env is invalid")
    end
else
    env = "dev"
end

# We start demand at 4AM
TIME_OFFSET = 3

Logging.info("Loading environment $env")
config = env_config[env]

# Setup
include("Manager/tasking.jl")
include("Manager/server.jl")
include("Manager/calculation.jl")
include("Manager/fetching.jl")
include("Manager/setter.jl")
include("Manager/helpers.jl")
include("Manager/filters.jl")
include("Manager/messages.jl")

# Unit Tests
include("Manager/unit_tests/helpers.jl")
include("Manager/unit_tests/filters.jl")
include("Manager/unit_tests/test.jl")

# Functional Tests
include("Manager/functional_tests/test.jl")

end#module

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module KeplerTools

using StaticArrays
using Rotations
using LinearAlgebra
using Optim
using PlotlyJS

export Orbit
export CelestialBody, Star, SpaceObject
export period, orbital_angle, orbital_distance, flight_path_angle, orbital_direction
export orbital_position, orbital_velocity, state_vector
export time_orbital_position, time_orbital_velocity, time_state_vector
export kepler, time_to_mean, mean_to_time, time_to_true, true_to_time
export basis_rotation, orientation_rotation
export angle_in_plane

export p_lambert, plane_change_p_lambert
export departure_orbit, arrival_orbit, fast_departure_orbit, fast_arrival_orbit
export ejection_angle, ejection_time, ejection_position, ejection_velocity, ejection_state_vector
export insertion_angle, insertion_time, insertion_position, insertion_velocity, insertion_state_vector
export propagate_across_soi
export Burn, orientation_components, apply_burn
export Transfer, fastTransfer, match_transfer_patch_times, match_patch_positions
export Porkchop, fastPorkchop

export draw_orbit, draw_central_body, draw_orbiting_body, draw_system
export draw_burn_arrow, draw_orbit_position, draw_angle_in_plane
export draw_transfer, draw_ejection, draw_insertion
export draw_porkchop

include("orbit.jl")
include("body.jl")
include("kepler.jl")
include("anomaly.jl")
include("basis.jl")
include("angles.jl")
include("utils.jl")

include("transfer/lambert.jl")
include("transfer/departarrive.jl")
include("transfer/soipropagate.jl")
include("transfer/burn.jl")
include("transfer/transfer.jl")
include("transfer/porkchop.jl")

include("draw/draworbit.jl")
include("draw/drawobject.jl")
include("draw/drawsystem.jl")
include("draw/drawmarkers.jl")
include("draw/drawtransfer.jl")
include("draw/drawporkchop.jl")

end


text:
VML = Vectorize.VML

println("===== Testing Intel's VML =====")
# Length of test array
N = 1000
Random.seed!(13)

# Test basic utility functionality
@testset "VML: Basic Functionality" begin
    mode = VML.getmode()
    @test mode == VML.setmode(VML.VML_HA)
    @test typeof(VML.setmode(VML.VML_HA)) <: Integer
end

# Real Arithmetic
for T in [Float32, Float64]
    @testset "VML: Basic Arithmetic::$T" begin
        X = convert(Vector{T}, abs.(randn(N)))
        Y = convert(Vector{T}, abs.(randn(N)))
        @testset "Testing $f::$T" for (f, fbs) in [(:add, +), (:sub, -), (:div, /), (:mul, *),
                                                  (:hypot, :hypot), (:atan, :atan), (:pow, ^),
                                                  (:max, :max), (:min, :min)]
            @eval fbase = $fbs
            @eval fvml = Vectorize.VML.$f
            @test fbase.(X, Y) ≈ fvml(X, Y)
        end
    end
end

# Complex Arithmetic
for T in [Complex{Float32}, Complex{Float64}]
    @testset "VML: Complex Arithmetic::$T" begin
        X = convert(Vector{T}, randn(N)) .+ convert(Vector{T}, randn(N))*im
        Y = convert(Vector{T}, randn(N)) .+ convert(Vector{T}, randn(N))*im
        @testset "Testing $f::$T" for (f, fbs) in [(:add, +), (:sub, -), (:div, /), (:mul, *),
                                                  (:pow, ^)]
            @eval fbase = $fbs
            @eval fvml = Vectorize.VML.$f
            @test fbase.(X, Y) ≈ fvml(X, Y)
        end

        @testset "Testing $f::$T" for f in [:dot]
            @eval fvml = Vectorize.VML.$f
            @test dot.(X, Y) ≈ fvml(X, Y)
        end
    end
end

# Scalar Power
for T in [Float32, Float64, Complex{Float32}, Complex{Float64}]
    @testset "VML: Complex Arithmetic::$T" begin
        X = rand(T, N)
        a = rand(T)
        @testset "Testing $f::$T" for (f, fbs) in ((:pow, ^),)
            @eval fbase = $fbs
            @eval fvml = Vectorize.VML.$f
            @test fbase.(X, a) ≈ fvml(X, a)
        end
    end
end

# Real Vector Math
for T in [Float32, Float64]
    @testset "VML: Real Vector Math::$T" begin
        X = convert(Vector{T}, randn(N)) # unrestricted
        Y = convert(Vector{T}, abs.(randn(N))) .+ 1 # greater than 1
        Z = clamp.(convert(Vector{T}, randn(N)), -1,  1) # [-1, 1]
        W = clamp.(convert(Vector{T}, randn(N)), 0,  2) # [0, 2]
        @testset "Testing $f::$T" for f in [:exp, :abs, :ceil, :floor, :round,
                                            :trunc, :cos, :sin, :tan, :cosh, :sinh, :tanh,
                                            :cbrt, :erf, :erfc, :lgamma, :gamma,
                                            :sinpi, :cospi, :sind, :cosd, :tand,
                                            :exp2, :exp10, :expm1, :cis]
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(X) ≈ fvml(X)
        end

        @testset "Testing $f::$T" for f in [:sqrt, :log10, :acosh, :asinh, :log,
                                            :log2, :log1p]
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(Y) ≈ fvml(Y)
        end

        @testset "Testing $f::$T" for f in [:erfcinv]
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(W) ≈ fvml(W)
        end

        @testset "Testing $f::$T" for f in [:acos, :asin, :atan, :atanh, :erfinv]
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(Z) ≈ fvml(Z)
        end
    end
end

# Complex Vector Math
for T in [Complex{Float32}, Complex{Float64}]
    @testset "VML: Complex Vector Functions::$T" begin
        X = convert(Vector{T}, complex(randn(N)))
        @testset "Testing $f::$T" for f in [:exp, :abs, :cos, :sin, :tan, :cosh, :sinh, :tanh,
                                            :sqrt, :log10, :acosh, :asinh, :log, :acos, :asin,
                                            :atan, :atanh, :conj]
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(X) ≈ fvml(X)
        end

        @testset "Testing $f::$T" for f in []
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(Y) ≈ fvml(Y)
        end

        @testset "Testing $f::$T" for f in []
            @eval fb = $f
            @eval fvml = Vectorize.VML.$f
            @test fb.(Z) ≈ fvml(Z)
        end
    end
end

println("===== VML Tests Successful =====\n\n")


text:
"""
Module MSM implements Calvet and Fisher's Markov Switching Multifractal model.

# Example usge:
```jldocs
julia> m = MSMmodel()
MSMmodel
  k: Int64 4
  b: Float64 5.0
  m₀: Float64 1.7
  γₖ: Float64 0.9
  σ₀: Float64 0.01
# Fit model to data x
julia> fit!(m, x);
# Predict volatility of x
julia> s = predict(m);
julia> histogram(s)
# Simulate a model
julia> simulate(m, nsims = 10000)
```
"""
module MSM

export MSMmodel, fit!, fitglobal!, predict, simulate, tune

using Parameters, DataFrames, Statistics, Random, Optim, Kronecker, LinearAlgebra, Distributions, BlackBoxOptim

const PPI = 1.0 / sqrt(2π)

"""
    MSMmodel(k, b, m₀, γₖ, σ₀)

Constructor for Markov Switching Multifractal (MSM) model.

# Example:
```jldocs
julia> m = MSMmodel()
MSMmodel
  k: Int64 4
  b: Float64 5.0
  m₀: Float64 1.7
  γₖ: Float64 0.9
  σ₀: Float64 0.01
```
"""
@with_kw mutable struct MSMmodel
    @deftype Float64
    k::Int = 4
    @assert k >= 1
    b = 5.0
    @assert b > 1.0
    m₀ = 1.7
    @assert m₀ > 1.0 && m₀ < 2.0
    γₖ = 0.90
    @assert γₖ > 0.0 && γₖ < 1.0
    σ₀ = 0.01
    @assert σ₀ > 0.0
end

# transition probability matrix
function Φ(b, γₖ, k::Integer)
    γ = [1 - (1 - γₖ)^(b^(j - k)) for j ∈ 1:k] .* 0.5
    return Matrix(reduce(⊗, [[1.0-u u; u 1.0-u] for u ∈ γ]))
end

# state dependent volatility
Σ(σ₀, m₀, k::Integer, M::Vector{Int}) = [σ₀ * sqrt(((2.0 - m₀)^u) * (m₀^(k - u))) for u ∈ M]

Σ(σ₀, m₀, k::Integer, m::Int) = σ₀ * sqrt(((2.0 - m₀)^m) * (m₀^(k - m)))

"""
    fit!(model::MSMmodel, x::Vector{Float64}; time_limit=60., solver=NelderMead(), g_tol=1.0e-6)

Fit MSM model to a vector of real-valued, de-meaned data using Optim's optimization methods. Recommended method is the default method.

# Example:
```jldocs
julia> using MSM;
julia> m = MSMmodel(); fit!(m, x)
```
"""
function fit!(model::MSMmodel, x::Vector{Float64}; time_limit = 60.0, solver = NelderMead(), g_tol = 1.0e-6)
    @unpack k, b, m₀, γₖ, σ₀ = model
    p0 = [b, m₀, γₖ, σ₀]
    l = [1.001, 1.001, 0.001, 0.0001]
    h = [50.0, 1.999, 0.999, 5.0]
    M = sum.([digits(q, base = 2, pad = k) for q = 0:2^k-1])
    N = length(M)
    Π = fill(1.0 / N, N)
    A = Φ(b, γₖ, k)
    res = optimize(p -> nll(p, k, M, A, Π, x), l, h, p0, Fminbox(solver), Optim.Options(time_limit = time_limit, g_tol = g_tol))
    b, m₀, γₖ, σ₀ = res.minimizer
    @pack! model = b, m₀, γₖ, σ₀
    return res
end

"""
    fitglobal!(model::MSMmodel, x::Vector{Float64}; method = :adaptive_de_rand_1_bin_radiuslimited,  maxsteps = 10000, tracemode = :silent)

Fit MSM model to a vector of real-valued, de-meaned data using BlackBoxOptim's global optimization methods. Recommended method is the default method.

# Example:
```jldocs
julia> using MSM;
julia> m = MSMmodel(); fitglobal!(m, x)
```

tracemode = :silent suppresses optimizer trace messages.
"""
function fitglobal!(model::MSMmodel, x::Vector{Float64}; method = :adaptive_de_rand_1_bin_radiuslimited, maxsteps = 10000, tracemode = :silent)
    @unpack k, b, m₀, γₖ, σ₀ = model
    l = [1.001, 1.001, 0.001, 0.0001]
    h = [50.0, 1.999, 0.999, 5.0]
    M = sum.([digits(q, base = 2, pad = k) for q = 0:2^k-1])
    N = length(M)
    Π = fill(1.0 / N, N)
    A = Φ(b, γₖ, k)
    res = bboptimize(p -> nll(p, k, M, A, Π, x), SearchRange = collect(zip(l, h)), Method = :adaptive_de_rand_1_bin_radiuslimited, MaxSteps = maxsteps, TraceMode = tracemode)
    b, m₀, γₖ, σ₀ = best_candidate(res)
    @pack! model = b, m₀, γₖ, σ₀
    return res
end

"""
    nll(p::Vector{<:Real}, k::Integer, M::Vector{Int}, A::Matrix{Float64}, Π::Vector{Float64}, rt::Vector{Float64})

compute negative log-likelihood of MSM.
"""
function nll(p::Vector{<:Real}, k::Integer, M::Vector{Int}, A::Matrix{Float64}, Π::Vector{Float64}, rt::Vector{Float64})
    b, m₀, γₖ, σ₀ = p
    A .= Φ(b, γₖ, k)
    σ = Σ(σ₀, m₀, k, M)
    ll = 0.0
    for r in rt
        ωr = [PPI * exp(-0.5 * r * r / (s * s)) / s for s in σ]
        w = ωr .* (A * Π)
        sw = sum(w)
        ll += log(sw)
        Π = w ./ sw
    end
    return -ll
end

"""
    predict(model::MSMmodel, paths, pathlen)

Predict volatility implied by an MSM model through simulation. `pathlen` is the length of observations over which volatility (standard deviation) is computed and `paths` is the number of paths simulated. Returns a vector of length nsamples containing volatilities.
"""
function predict(model::MSMmodel, paths, pathlen)
    s = fill(0.0, paths)
    @inbounds for i = 1:paths
        r = simulate(model, pathlen)
        s[i] = std(r)
    end
    return s
end

"""
    simulate(model::MSMmodel, pathlen)

Simulate an MSM model's trajectory of length pathlen. Returns a vector.
"""
function simulate(model::MSMmodel, pathlen)
    @unpack k, b, m₀, γₖ, σ₀ = model
    A = Φ(b, γₖ, k)
    M = sum.([digits(q, base = 2, pad = k) for q = 0:2^k-1])
    N = length(M)
    u = rand(Categorical(N))
    r = Vector{Float64}(undef, pathlen)
    r[1] = Σ(σ₀, m₀, k, M[u]) * randn()
    @inbounds for j = 2:pathlen
        u = rand(Categorical(A[u, :]))
        r[j] = Σ(σ₀, m₀, k, M[u]) * randn()
    end
    return r
end

"""
    tune(x::Vector{<:Real}; ks=2:12)

Fit MSM models for different k's and output results. k > 10 may take more than a minute to fit. Returns a dataframe.

"""
function tune(x::Vector{<:Real}; ks = 2:12)
    res = DataFrame()
    for j in ks
        model = MSM(k = j, σ₀ = std(x))
        out = @timed fit!(model, x)
        @unpack k, b, m₀, γₖ, σ₀ = model
        push!(res, (; ([:exec_time, :loglik, :k, :b, :m₀, :γₖ, :σ₀] .=> [out.time, -out.value.minimum, k, b, m₀, γₖ, σ₀])...))
        println("k = $(k) took $(out.time) seconds.")
    end
    return res
end

end


text:
using DataStructures,SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
import Base.cat
import Base.∘
using ViewerGL
GL = ViewerGL

# Graphic text font implementation
# Vector definition of printable ASCII codes as 1D LAR models.
# Font design from [Geometric Programming for Computer-Aided Design](https://www.amazon.com/Geometric-Programming-Computer-Aided-Design/dp/0471899429?_encoding=UTF8&%2AVersion%2A=1&%2Aentries%2A=0), Wiley, 2003.
"""
	hpcs::Array

"""
const hpcs = Tuple{Array{Float64,N} where N,Array{Array{Int64,1},1}}[([0.0, 0.0], [[1]]), ([1.75 1.75 2.0 2.0 1.5 1.5; 1.75 5.5 0.5 1.0 1.0 0.5], [[1, 2], [3, 4], [4, 5], [5, 6], [6, 3]]), ([1.0 2.0 2.0 2.0 1.5 1.5 2.0 3.0 3.0 3.0 2.5 2.5; 4.0 5.0 5.5 6.0 6.0 5.5 4.0 5.0 5.5 6.0 6.0 5.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 9]]), ([1.0 3.0 1.0 3.0 1.25 1.75 2.25 2.75; 2.5 2.5 3.5 3.5 1.75 4.0 1.75 4.0], [[1, 2], [3, 4], [5, 6], [7, 8]]), ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 2.0 2.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 -0.5 6.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [13, 14]]), ([2.5 2.5 2.0 2.0 2.5 2.5 2.0 2.0 1.5 3.5; 0.0 0.5 0.5 0.0 5.5 6.0 6.0 5.5 5.5 11.5], [[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 5], [9, 10]]), ([4.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 5.0 4.0 2.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14]]), ([1.0 2.0 2.0 2.0 1.5 1.5; 4.0 5.0 5.5 6.0 6.0 5.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3]]), ([2.0 1.0 0.5 1.0 2.0; 0.0 1.0 3.0 5.0 6.0], [[1, 2], [2, 3], [3, 4], [4, 5]]), ([2.0 3.0 3.5 3.0 2.0; 0.0 1.0 3.0 5.0 6.0], [[1, 2], [2, 3], [3, 4], [4, 5]]), ([1.0 3.0 2.0 2.0 1.0 3.0 1.0 3.0; 3.0 3.0 2.0 4.0 2.0 4.0 4.0 2.0], [[1, 2], [3, 4], [5, 6], [7, 8]]), ([1.0 3.0 2.0 2.0; 3.0 3.0 2.0 4.0], [[1, 2], [3, 4]]), ([1.0 2.0 2.0 2.0 1.5 1.5; -1.0 0.0 0.5 1.0 1.0 0.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 3]]), ([1.0 3.0; 3.0 3.0], [[1, 2]]), ([2.0 2.0 1.5 1.5 2.0; 0.0 0.5 0.5 0.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5]]), ([1.0 3.0; 0.0 6.0], [[1, 2]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [4, 8]]), ([0.0 2.0 2.0 0.0 4.0; 4.0 6.0 0.0 0.0 0.0], [[1, 2], [2, 3], [4, 5]]), ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 4.0 5.0 6.0 6.0 5.0 4.0 0.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]]), ([0.0 4.0 2.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 2.0 1.0 0.0 0.0 1.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]]), ([4.0 0.0 4.0 4.0; 1.0 1.0 6.0 0.0], [[1, 2], [2, 3], [3, 4]]), ([4.0 0.0 0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0; 6.0 6.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10]]), ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 6.0 6.0 5.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11]]), ([0.0 0.0 4.0 0.0; 5.0 6.0 6.0 0.0], [[1, 2], [2, 3], [3, 4]]), ([1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 4.0 4.0 3.0 1.0 0.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 4.0 5.0 6.0 6.0 5.0 4.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [6, 5], [5, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14], [14, 6]]), ([0.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 5.0 6.0 6.0 5.0 3.0 2.0 2.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11]]), ([2.0 2.0 1.5 1.5 2.0 2.0 1.5 1.5; 1.0 1.5 1.5 1.0 3.0 3.5 3.5 3.0], [[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 5]]), ([2.0 2.0 1.5 1.5 1.0 2.0 2.0 2.0 1.5 1.5; 3.0 3.5 3.5 3.0 -0.5 0.5 1.0 1.5 1.5 1.0], [[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 7]]), ([3.0 0.0 3.0; 6.0 3.0 0.0], [[1, 2], [2, 3]]), ([1.0 3.0 1.0 3.0; 2.5 2.5 3.5 3.5], [[1, 2], [3, 4]]), ([1.0 4.0 1.0; 6.0 3.0 0.0], [[1, 2], [2, 3]]), ([2.0 2.0 1.5 1.5 1.75 1.75 3.0 3.0 2.0 1.0 0.0 0.0; 0.0 0.5 0.5 0.0 1.0 2.75 4.0 5.0 6.0 6.0 5.0 4.0], [[1, 2], [2, 3], [3, 4], [4, 1], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12]]), ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 2.0 1.0 2.0 3.0 2.0; 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 1.0 2.0 3.0 2.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13]]), ([0.0 0.0 1.0 3.0 4.0 4.0 0.0 4.0; 0.0 5.0 6.0 6.0 5.0 0.0 2.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [7, 8]]), ([0.0 0.0 3.0 4.0 4.0 3.0 4.0 4.0 3.0 0.0 0.0 3.0; 0.0 6.0 6.0 5.0 4.0 3.0 2.0 1.0 0.0 0.0 3.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [11, 12]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]]), ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 1.0 0.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([4.0 0.0 0.0 4.0 0.0 3.0; 0.0 0.0 6.0 6.0 3.0 3.0], [[1, 2], [2, 3], [3, 4], [5, 6]]), ([0.0 0.0 4.0 0.0 3.0; 0.0 6.0 6.0 3.0 3.0], [[1, 2], [2, 3], [4, 5]]), ([2.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10]]), ([0.0 0.0 4.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0 3.0 3.0], [[1, 2], [3, 4], [5, 6]]), ([2.0 2.0 1.0 3.0 1.0 3.0; 0.0 6.0 0.0 0.0 6.0 6.0], [[1, 2], [3, 4], [5, 6]]), ([0.0 1.0 2.0 3.0 3.0 2.0 4.0; 1.0 0.0 0.0 1.0 6.0 6.0 6.0], [[1, 2], [2, 3], [3, 4], [4, 5], [6, 7]]), ([4.0 0.0 4.0 0.0 0.0; 6.0 3.0 0.0 0.0 6.0], [[1, 2], [2, 3], [4, 5]]), ([4.0 0.0 0.0; 0.0 0.0 6.0], [[1, 2], [2, 3]]), ([0.0 0.0 2.0 4.0 4.0; 0.0 6.0 4.0 6.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5]]), ([0.0 0.0 4.0 4.0 4.0; 0.0 6.0 2.0 0.0 6.0], [[1, 2], [2, 3], [4, 5]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]]), ([0.0 0.0 3.0 4.0 4.0 3.0 0.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 4.0; 1.0 0.0 0.0 1.0 5.0 6.0 6.0 5.0 1.0 1.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [10, 11]]), ([0.0 0.0 3.0 4.0 4.0 3.0 0.0 3.0 4.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 2.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [8, 9]]), ([0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12]]), ([2.0 2.0 0.0 4.0; 0.0 6.0 6.0 6.0], [[1, 2], [3, 4]]), ([0.0 0.0 1.0 3.0 4.0 4.0; 6.0 1.0 0.0 0.0 1.0 6.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]), ([0.0 2.0 4.0; 6.0 0.0 6.0], [[1, 2], [2, 3]]), ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 6.0 3.0 0.0 3.0 0.0 3.0 6.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([0.0 4.0 0.0 4.0; 0.0 6.0 6.0 0.0], [[1, 2], [3, 4]]), ([0.0 2.0 4.0 2.0 2.0; 6.0 2.0 6.0 2.0 0.0], [[1, 2], [2, 3], [4, 5]]), ([0.0 4.0 0.0 4.0; 6.0 6.0 0.0 0.0], [[1, 2], [2, 3], [3, 4]]), ([2.0 1.0 1.0 2.0; 0.0 0.0 6.0 6.0], [[1, 2], [2, 3], [3, 4]]), ([1.0 3.0; 6.0 0.0], [[1, 2]]), ([2.0 3.0 3.0 2.0; 0.0 0.0 6.0 6.0], [[1, 2], [2, 3], [3, 4]]), ([1.0 2.0 3.0; 5.0 6.0 5.0], [[1, 2], [2, 3]]), ([1.0 4.0; 0.0 0.0], [[1, 2]]), ([2.0 2.0 3.0 2.5 2.5 2.0; 4.5 5.0 6.0 4.0 4.5 4.0], [[1, 2], [2, 3], [4, 5], [5, 1], [1, 6], [6, 4]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0 1.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10], [10, 11]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0 3.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10], [10, 11]]), ([4.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 0.0; 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]]), ([4.0 4.0 3.0 2.0 1.0 1.0 0.0 2.0; 3.0 4.0 5.0 5.0 4.0 0.0 1.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [7, 8]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 3.0 1.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 0.0 -1.0 -1.0 0.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [1, 9], [9, 10], [10, 11], [11, 12]]), ([4.0 4.0 3.0 1.0 0.0 0.0 0.0 1.0; 0.0 2.0 3.0 3.0 2.0 0.0 5.0 5.0], [[1, 2], [2, 3], [3, 4], [4, 5], [6, 7], [7, 8]]), ([1.0 3.0 1.0 3.0 2.0 2.0 2.25 2.25 1.75 1.75; 0.0 0.0 3.0 3.0 0.0 3.0 3.75 4.25 4.25 3.75], [[1, 2], [3, 4], [5, 6], [7, 8], [8, 9], [9, 10], [10, 7]]), ([1.0 3.0 2.0 2.0 1.0 0.0 2.25 2.25 1.75 1.75; 3.0 3.0 3.0 0.0 -1.0 0.0 3.75 4.25 4.25 3.75], [[1, 2], [3, 4], [4, 5], [5, 6], [7, 8], [8, 9], [9, 10], [10, 7]]), ([0.0 1.0 1.0 0.0 4.0 2.0 1.0 3.0 4.0; 0.0 0.0 3.0 3.0 0.0 0.0 1.0 3.0 3.0], [[1, 2], [2, 3], [3, 4], [5, 6], [6, 7], [7, 8], [8, 9]]), ([2.0 2.0 1.0 1.0 3.0; 0.0 5.0 5.0 0.0 0.0], [[1, 2], [2, 3], [4, 5]]), ([4.0 4.0 2.0 2.0 2.0 0.0 0.0 0.0; 0.0 3.0 2.0 0.0 3.0 2.0 0.0 3.0], [[1, 2], [2, 3], [4, 5], [5, 6], [7, 8]]), ([3.0 3.0 1.0 1.0 1.0; 0.0 3.0 2.0 0.0 3.0], [[1, 2], [2, 3], [4, 5]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 0.0 0.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]]), ([4.0 3.0 1.0 0.0 0.0 1.0 3.0 4.0 4.0 4.0; 1.0 0.0 0.0 1.0 2.0 3.0 3.0 2.0 3.0 -1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 1], [9, 10]]), ([0.0 2.0 1.0 1.0 1.0 2.0 3.0 4.0; 0.0 0.0 0.0 3.0 2.0 3.0 3.0 2.0], [[1, 2], [3, 4], [5, 6], [6, 7], [7, 8]]), ([0.0 4.0 3.0 1.0 0.0 1.0 3.0 4.0; 0.0 0.0 1.0 1.0 2.0 3.0 3.0 2.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8]]), ([1.0 3.0 2.0 2.0 2.0 3.0; 0.0 0.0 0.0 5.0 4.0 4.0], [[1, 2], [3, 4], [5, 6]]), ([0.0 1.0 1.0 2.0 3.0 4.0 4.0; 3.0 3.0 1.0 0.0 0.0 1.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([0.0 1.0 3.0 4.0; 3.0 0.0 3.0 3.0], [[1, 2], [2, 3], [3, 4]]), ([0.0 0.0 1.0 2.0 3.0 4.0 4.0; 3.0 2.0 0.0 2.0 0.0 2.0 3.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([0.0 1.0 4.0 1.0 4.0; 3.0 3.0 0.0 0.0 3.0], [[1, 2], [2, 3], [4, 5]]), ([0.0 1.0 2.5 0.0 1.0 4.0; 3.0 3.0 1.5 0.0 0.0 3.0], [[1, 2], [2, 3], [4, 5], [5, 6]]), ([0.0 0.0 3.0 0.0 3.0 4.0; 2.0 3.0 3.0 0.0 0.0 1.0], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]), ([2.5 2.0 2.0 1.5 2.0 2.0 2.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([2.0 2.0; 0.0 5.0], [[1, 2]]), ([1.5 2.0 2.0 2.5 2.0 2.0 1.5; 6.5 6.0 3.5 3.0 2.5 0.0 -0.5], [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]), ([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5], [[1, 2], [2, 3], [3, 4]])]

"""
ascii_LAR = DataStructures.OrderedDict{Int,Lar.LAR}()

	GL.ascii2lar::{Int,Lar.LAR}

*Ordered dictionary* of printable ASCII codes as one-dimensional *LAR models* in 2D.

`Key` is the `integer` ASCII code between 32 and 126.
# Example
```
julia> GL.ascii2lar[46]
([2.0 2.0 … 1.5 2.0; 0.0 0.5 … 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])

julia> GL.ascii2lar[126]
([1.0 1.75 2.75 3.5; 5.0 5.5 5.0 5.5], Array{Int64,1}[[1, 2], [2, 3], [3, 4]])
```
"""
const ascii2lar = DataStructures.OrderedDict(zip(32:126,hpcs))

# Attributes for `text` 2D graphics primitive.
# Reference to GKS ISO Graphics Standard (ISO/IEC 7942-4:1998)

const textalignment = "centre" #default value
const textangle = pi/4 #default value
const textwidth = 1.0 #default value
const textheight = 1.0 #default value
const textspacing = 0.25 #default value
const fontwidth = 4.0 #default value
const fontheight = 8.0 #default value
const fontspacing = 1.0 #default value

"""
	charpols(mystring::String)::Array{LAR,1}

Return the `array` of `LAR` models, for chars in `mystring`.
"""
function charpols(mystring)
    return [ascii2lar[code] for code in [Int(char) for char in mystring]]
end

"""
	charseq(mystring::String)::Array{Char,1}

# Example
```
julia> GL.charseq("PLaSM")
5-element Array{Char,1}:
 'P'
 'L'
 'a'
 'S'
 'M'
```
"""
function charseq(mystring)
	return [char for char in mystring]
end

"""
	text(mystring::String)::LAR

Compute the one-dim *LAR model* drawing the contents of `mystring`

# Example
```
julia> model = GL.text("PLaSM")
# output
([0.0 0.0 3.0 4.0 4.0 3.0 0.0 9.0 5.0 5.0 14.0 13.0 11.0 10.0
10.0 11.0 13.0 14.0 14.0 14.0 15.0 16.0 18.0 19.0 19.0 18.0 16.0 15.0 15.0 16.0 18.0
19.0 20.0 20.0 22.0 24.0 24.0; 0.0 6.0 6.0 5.0 3.0 2.0 2.0 0.0 0.0 6.0 1.0 0.0 0.0 1.0
2.0 3.0 3.0 2.0 0.0 3.0 1.0 0.0 0.0 1.0 2.0 3.0 3.0 4.0 5.0 6.0 6.0 5.0 0.0 6.0 4.0
6.0 0.0],
Array{Int64,1}[[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[8,9],[9,10],[11,12],[12,13],[13,14],
[14,15],[15,16],[16,17],[17,18],[18,11],[19,20],[21,22],[22,23],[23,24],[24,25],[25,26],
[26,27],[27,28],[28,29],[29,30],[30,31],[31,32],[33,34],[34,35],[35,36],[36,37]])

julia> GL.VIEW([GL.GLLar2gl(model...)])
```
"""
function text(mystring,flag=true)
	V,EV = GL.comp([ Lar.struct2lar, Lar.Struct, cat, distr,
			GL.cons([ charpols, k(Lar.t(fontspacing+fontwidth,0)) ]),charseq ])(mystring)
	out = GL.normalize3(V,flag),EV
	return out
end

"""
	GLText(string)::GL.GLMesh

Transform a string into a mesh of lines.
To display as graphical text.
# Example
```
julia> GL.text("Plasm")
([0.0 0.0 … 0.833333 0.833333; 0.0 0.25 … 0.0 0.125; 0.0 0.0 … 0.0 0.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]
, [6, 7], [8, 9], [9, 10], [11, 12], [13, 14]  …  [25, 26], [26, 27], [27, 28], [28, 29], [29, 30], [31, 32], [32, 33], [34, 3
5], [35, 36], [37, 38]])

julia> GL.GLText("Plasm")
ViewerGL.GLMesh(1, [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0], ViewerGL.GLVertexArray(-1), ViewerGL.GLVertexBuffer(-1, Float32[0.0, 0.25, 0.0, 0.0, 0.0, 0.0, 0.125, 0.25, 0.0, 0.0  …  0.0, 0.916667, 0.125, 0.0, 0.833333, 0.125, 0.0, 0.833333, 0.0, 0.0]), ViewerGL.GLVertexBuffer(-1, Float32[]), ViewerGL.GLVertexBuffer(-1, Float32[]))
```
"""
function GLText(string,color=COLORS[12])::GL.GLMesh
	GL.GLLines(GL.text(string)...,color)
end

"""
	a2a(mat::Matrix)(models::Array{LAR,1})::Struct.body

Local function, service to `textWithAttributes` implementation.
"""
function a2a(mat)
	function a2a0(models)
		assembly = []
		for model in models
			push!( assembly, Lar.Struct([ mat,model ]) )
		end
		assembly
	end
	return a2a0
end

"""
	translate(c::Number)(lar::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function translate(c)
	function translate0(lar)
		xs = lar[1][1,:]
		width = maximum(xs) - minimum(xs)
		GL.apply(Lar.t(width/c,0))(lar)
	end
	return translate0
end

"""
	align(textalignment="centre"::String)(model::LAR)::LAR

Local function, service to `textWithAttributes` implementation.
"""
function align(textalignment="centre")
	function align1(model)
		if ( textalignment == "centre" ) out=translate(-2)(model)
		elseif ( textalignment == "right" ) out=translate(-1)(model)
		elseif ( textalignment == "left" ) out=model
		end
		return out
	end
end

"""
 	cat(args)

Redefined locally, as service to `textWithAttributes` implementation.
```
julia> cat([[1,2],[3,4],[],[5,6,7,8]])==collect(1:8)
true
```
"""
function cat(args)
	return reduce( (x,y) -> append!(x,y), args; init=[] )
end

"""
	textWithAttributes(textalignment='centre', textangle=0,
		textwidth=1.0, textheight=2.0, textspacing=0.25)(strand::String)::LAR

Partial implementation of the GKS's graphics primitive `text`.

# Example

```
GL.VIEW([
	GL.GLLines(GL.textWithAttributes("centre", pi/4)("PLaSM")...),
	GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1))
	])
```
GL.VIEW([
	GL.GLLines(GL.textWithAttributes("left", pi/4)("PLaSM")...),
	GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1))
])
"""
function textWithAttributes(textalignment="centre", textangle=0,
							textwidth=1.0, textheight=2.0, textspacing=0.25)
	function textWithAttributes(strand)
		id = x->x
		mat = Lar.s(textwidth/fontwidth,textheight/fontheight)
		comp([
		   GL.apply(Lar.r(textangle)),
		   align(textalignment),
		   Lar.struct2lar,
		   Lar.Struct,
		   cat,
		   distr,
		   cons([ a2a(mat) ∘ charpols,
				k(Lar.t(textwidth+textspacing,0)) ]),
		   charseq ])(strand)
	end
end

"""
	embed(n::Int)(model::LAR)::LAR

Embed a `LAR` `model` of dimension ``d`` in a space ``R^{d+n}``.
The embedding is done by adding ``d`` zero coordinates to each vertex.

# Example

```
julia> square = Lar.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> GL.embed(1)(square)
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0; 0.0 0.0 0.0 0.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function embed(n)
	function embed0( larmodel )
		return [larmodel[1]; zeros(n,size(larmodel[1],2))], larmodel[2]
	end
	return embed0
end

"""
	numbering(model::LARmodel)

Display a *wireframe* of `model` cells, of any dimension, with their *ordinal number*.
Different `colors` and size are used for the various dimensional cells.

# Examples

```
model = Lar.cuboidGrid([3,4,2], true)
Plasm.view(Plasm.numbering()(model))

model = Lar.cuboidGrid([10,10], true)
Plasm.view(Plasm.numbering(1.5)(model))
```
"""
function numbering(numberSizeScaling=1.)
	function numbering0(model,color=COLORS[1],alpha=0.2)
		V,cells = model
		meshes = []
		if length(cells)>2
			background = GL.GLHulls(V, cells[3], color, alpha)
		end
		if size(V,1)==2
			V = GL.embed(1)(model)[1]
		end
		wireframe = V,cells[2]
		ns = numberSizeScaling
		gcode = GL.textWithAttributes("centre", 0, 0.1ns, 0.2ns, 0.025ns)
		push!(meshes,GL.GLLines(wireframe[1],wireframe[2],color))

		colors = GL.COLORS[3], GL.COLORS[7], GL.COLORS[5], GL.COLORS[8]
		for (h,skel) in enumerate(cells)
	 	  nums = []
		  for (k,cell) in enumerate(skel)
				center = sum([V[:,v] for v in cell])/length(cell)
				code = GL.embed(1)( gcode(string(k)) )
				scaling = (0.6+0.1h,0.6+0.1h,1)
				push!(nums, Lar.struct2lar( Lar.Struct([
					Lar.t(center...), Lar.s(scaling...), code ]) ))
		  end
		  for num in nums
			mesh = GL.GLLines(num[1],num[2],colors[h])
			push!( meshes, mesh )
		  end
		end
		if length(cells)>2 push!( meshes, background ) end
		return meshes
	end
	return numbering0
end

"""
#	numbering(scaling=0.1)
#		(V::Lar.Points, copEV::Lar.ChainOp, copFE::Lar.ChainOp)::Lar.Hpc

Produce the numbered `Hpc` of `planar_arrangement()` 2D output.
Vertices in `V` are stored by row.

# Example

```julia
using LinearAlgebraicRepresentation
using Plasm, SparseArrays
Lar = LinearAlgebraicRepresentation

V,EV = Lar.randomcuboids(7, 1.0);
V = GL.normalize(V,flag=true);
W = convert(Lar.Points, V');
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
V, copEV, copFE = Lar.planar_arrangement(W, cop_EW);

Plasm.view( Plasm.numbering(0.05)((V, copEV, copFE)) )
```
"""
function numbering1(scaling=0.1)
	function numbering0(model::Tuple{Lar.Points,Lar.ChainOp,Lar.ChainOp})
		(V, copEV, copFE) = model
		VV = [[k] for k=1:size(V,1)]
		EV = [SparseArrays.findnz(copEV[h,:])[1] for h=1:size(copEV,1)]
		FV = [collect(Set(cat(EV[e] for e in SparseArrays.findnz(copFE[i,:])[1]))) for i=1:size(copFE,1)]
		FV = convert(Array{Array{Int64,1},1}, FV)
		model = (convert(Lar.Points, V'), Lar.Cells[VV,EV,FV])
		return GL.numbering(scaling)(model)
	end
	return numbering0
end

function cons(funs)
	return x -> [f(x) for f in funs]
end
function comp(funs)
	function compose(f,g)
	  return x -> f(g(x))
	end
	id = x->x
	return reduce(compose, funs; init=id)
end
function distr(args)
	list,element = args
	return [ [e,element] for e in list ]
end

"""
	k(Any)(x)

*Constant* functional of FL and PLaSM languages.

Gives a constant functional that always returns the actual parameter
when applied to another parameter.

#	Examples

```
julia> GL.k(10)(100)
10

julia> GL.k(sin)(cos)
sin
```
"""
function k(any::Any)
	x->any
end

"""
	apply(affineMatrix::Matrix)(larmodel::LAR)::LAR

Apply the `affineMatrix` parameter to the vertices of `larmodel`.

# Example

```
julia> square = LinearAlgebraicRepresentation.cuboid([1,1])
([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]])

julia> Plasm.apply(LinearAlgebraicRepresentation.t(1,2))(square)
([1.0 1.0 2.0 2.0; 2.0 3.0 2.0 3.0], Array{Int64,1}[[1, 2, 3, 4]])
```
"""
function apply(affineMatrix::Array{Float64,2})
	function apply0(larmodel)
		return Lar.apply( affineMatrix, larmodel )
	end
	return apply0
end
function apply(affineMatrix::Array{Float64,2},larmodel)
	return Lar.apply( affineMatrix, larmodel )
end

"""
	COMP(Funs)

Compose an Array of functions.
```
julia> COMP([sin,cos,tan])(pi/4)
0.5143952585235492

julia> COMP([sin,cos,tan])(pi/4) == (sin∘cos)(1)
true
```
"""
function COMP(Funs)
	function ∘(f,g)
		function h(x)
			return f(g(x))
		end
		return h
	end
	return reduce(∘,Funs)
end


text:
module Polynomials

#  using GenericLinearAlgebra ## remove for now. cf: https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/pull/71#issuecomment-743928205
using LinearAlgebra
using Intervals

if VERSION >= v"1.4.0"
    import Base: evalpoly
end

include("abstract.jl")
include("show.jl")
include("plots.jl")
include("contrib.jl")

# Interface for all AbstractPolynomials
include("common.jl")

# Polynomials
include("polynomials/standard-basis.jl")
include("polynomials/mutable-arithmetics.jl")
include("polynomials/Polynomial.jl")
include("polynomials/ImmutablePolynomial.jl")
include("polynomials/SparsePolynomial.jl")
include("polynomials/LaurentPolynomial.jl")
include("polynomials/pi_n_polynomial.jl")
include("polynomials/ngcd.jl")
include("polynomials/multroot.jl")

include("polynomials/ChebyshevT.jl")

# Rational functions
if VERSION >= v"1.2.0"
    include("rational-functions/common.jl")
    include("rational-functions/rational-function.jl")
    include("rational-functions/fit.jl")
    #include("rational-transfer-function.jl")
    include("rational-functions/plot-recipes.jl")
end

# compat; opt-in with `using Polynomials.PolyCompat`
include("polynomials/Poly.jl")

end # module


text:
import JLBoost: feature_importance

import MLJBase: CrossEntropy

CrossEntropy(x, y::CategoricalVector) = CrossEntropy(x, y.refs .- 1)


text:
# generate AMGX wrappers

using Clang

const LIBAMGX_INCLUDE = joinpath(homedir(), "Applications/AMGX/base/include")
const LIBAMGX_HEADERS = [joinpath(homedir(), "Applications/AMGX/base/include/amgx_c.h"),
                         joinpath(homedir(), "Applications/AMGX/base/include/amgx_config.h")]

wc = init(; headers = LIBAMGX_HEADERS,
            output_file = joinpath(@__DIR__, "libAMGX.jl"),
            common_file = joinpath(@__DIR__, "libAMGX_common.jl"),
            clang_includes = [LIBAMGX_INCLUDE,],
            clang_args = ["-I", joinpath(LIBAMGX_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libAMGX",
            )

function main()
    cd(joinpath(@__DIR__, "..", "lib")) do
        run(wc)
    end
end

# TODO, this seems to put things in res instead of lib
isinteractive() || main()


text:
using DynamicalSystemsToolkit
using ModelingToolkit
using NonlinearSolve

function hopf()
    @parameters t β σ
    @variables y₁(t) y₂(t)
    D = Differential(t)
    return ODESystem(
        [
            D(y₁) ~ β * y₁ - y₂ + σ * y₁ * (y₁^2 + y₂^2),
            D(y₂) ~ y₁ + β * y₂ + σ * y₂ * (y₁^2 + y₂^2),
        ],
        t,
        [y₁, y₂],
        [β, σ];
        defaults=[β => 1, σ => -1],
        name=:hopf,
    )
end

function hopf_test(; n_mesh=5, n_coll=4)
    coll = CollocationSystem(hopf(), n_mesh, n_coll)
    add_boundarycondition!(coll, (u₀, t₀, u₁, t₁, p) -> u₀ .- u₁)  # periodicity
    add_boundarycondition!(coll, (u₀, t₀, u₁, t₁, p) -> [0 ~ t₀, 0 ~ u₀[2]])  # start time and phase condition
    sol = t -> (cos(t), sin(t))
    ics = generate_initialconditions(coll, sol, (0, 2π))
    prob = NonlinearProblem(NonlinearSystem(coll), ics, [], jac=true)
    return (coll, prob)
end

function test2()
    @variables x y
    sys = NonlinearSystem([x ~ y], [x, y], []; defaults=[y => 1])
    return NonlinearProblem(sys, [x => 0])
end


text:
#=
Challenge #256 [Easy] Oblique and De-Oblique

The oblique function slices a matrix (2d array) into diagonals.
The de-oblique function takes diagonals of a matrix, and reassembles the original
rectangular one.

input for oblique

 0  1  2  3  4  5
 6  7  8  9 10 11
12 13 14 15 16 17
18 19 20 21 22 23
24 25 26 27 28 29
30 31 32 33 34 35
(and the output to de-oblique)

output for oblique

0
1 6
2 7 12
3 8 13 18
4 9 14 19 24
5 10 15 20 25 30
11 16 21 26 31
17 22 27 32
23 28 33
29 34
35
(and the input to de-oblique)

bonus deambiguated de-oblique matrices

There's only one de-oblique solution for a square matrix, but when the result is
not square, another input is needed to indicate whether the output should be tall
or wide or provide specific dimentsions of output:

rectangular oblique data input

0
1 6
2 7 12
3 8 13
4 9 14
5 10 15
11 16
17

output for (wide) deoblique (3 6, INPUT) or deoblique (WIDE, INPUT)

 0  1  2  3  4  5
 6  7  8  9 10 11
12 13 14 15 16 17

output for (tall) deoblique (6 3, INPUT) or deoblique (TALL, INPUT)

 0  1  2
 6  7  3
12  8  4
13  9  5
14 10 11
15 16 17

Note

The main use of these functions in computer science is to operate on the diagonals
of a matrix, and then revert it back to a rectangular form. Usually the rectangular
dimensions are known.

(source: https://redd.it/48a4pu)
=#

#   exit at EOF
read_matrix( data::DataType=Any, input::IO=STDIN, spacing::Char=' ' ) = readdlm(
input, spacing, data, header=false, skipblanks=true, ignore_invalid_chars=true,
quotes=false, comments=false )

function de_oblique( grid::Matrix{ Int } )
end

function oblique( grid::Matrix{ Int } )
    local len::Integer = size( grid, 1 )
    local line::Array{ Int, 1 } = Array{ Int, 1 }( )
    local tmp::Integer = 0
    local new_grid::Array{ Array{ Int, 1 }, 1 } = Array{ Array{ Int, 1 }, 1 }( )

    for i in 1:len
        tmp = i
        for j in 1: i
            push!( line, grid[ j , tmp ] )
            tmp -= 1
        end
        push!( new_grid, line )
        line = Array{ Int, 1 }( )
    end

    for i in 2:len
        tmp = i
        for j in len:-1:i
            push!( line, grid[ tmp, j ] )
            tmp += 1
        end
        push!( new_grid, line )
        line = Array{ Int, 1 }( )
    end

    return new_grid
end

function main( )
    local grid::Matrix{ Int } = read_matrix( Int )
    local new_grid::Array{ Array{ Int, 1 }, 1 } = Array{ Array{ Int, 1 }, 1 }( )
    local len::Integer = wid::Integer = 0
    ( len, wid ) = size( grid )

    if len != wid
        new_grid = de_oblique( grid )
    else
        new_grid = oblique( grid )
    end

    for i in 1:size( new_grid, 1 )
        println( new_grid[ i ] )
    end
end

main( )

max_stars_repo_name: farmfirst/dailyprogrammer

max_issues_repo_name: farmfirst/dailyprogrammer

max_forks_repo_name: farmfirst/dailyprogrammer


text:
import LinearAlgebra.generic_matmatmul!

const gemm_mt_threshold = Ref(64.0)
mt_thresholds[:gemm] = gemm_mt_threshold;

function generic_matmatmul!(C::AbstractMatrix{DoubleFloat{T}}, tA, tB, A::AbstractMatrix{DoubleFloat{T}}, B::AbstractMatrix{DoubleFloat{T}},
                            _add::MulAddMul=MulAddMul()) where T
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    mC, nC = size(C)

    if iszero(_add.alpha)
        return _rmul_or_fill!(C, _add.beta)
    end
    if mA == nA == mB == nB == mC == nC == 2
        return matmul2x2!(C, tA, tB, A, B, _add)
    end
    if mA == nA == mB == nB == mC == nC == 3
        return matmul3x3!(C, tA, tB, A, B, _add)
    end
    _xgeneric_matmatmul!(C, tA, tB, A, B, _add)
end
function _xgeneric_matmatmul!(C::StridedMatrix{DoubleFloat{T}},
                            tA::AbstractChar, tB::AbstractChar,
                            A::StridedMatrix{DoubleFloat{T}},
                            B::StridedMatrix{DoubleFloat{T}},
                            _add::MulAddMul
                            ) where {T <: AbstractFloat}

    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    if ((nthreads() > 1) && (Float64(mB)*Float64(mA) > gemm_mt_threshold[])
        && _add.alpha == 1 && _add.beta == 0)
        _mt_generic_matmatmul!(C,tA,tB,A,B)
    else
        __generic_matmatmul!(C,tA,tB,A,B,_add)
    end
    C
end

function _mt_generic_matmatmul!(C::StridedMatrix{DoubleFloat{T}},
              tA::AbstractChar, tB::AbstractChar,
              A::StridedMatrix{DoubleFloat{T}},
              B::StridedMatrix{DoubleFloat{T}}) where {T <: AbstractFloat}
    require_one_based_indexing(C, A, B)
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    if  mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA), matrix B has dimensions ($mB,$nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result matrix C dimensions $(size(C)), needs ($mA,$nB)"))
    end

    if tA == 'N'
        At = Matrix(transpose(A))
    else
        At = A
    end

    Blines = [zeros(DoubleFloat{T},mB) for id in 1:nthreads()]
    if tB == 'N'
        @threads  for j = 1:nB
            id = threadid()
            Bline = Blines[id]
            gemm_kernN(C, At, B, Bline, j, mB, mA)
        end
    else
        @threads  for j = 1:nB
            id = threadid()
            Bline = Blines[id]
            gemm_kernT(C, At, B, Bline, j, mB, mA)
        end
    end
    C
end

function __generic_matmatmul!(C::StridedMatrix{DoubleFloat{T}},
                            tA::AbstractChar, tB::AbstractChar,
                             A::StridedMatrix{DoubleFloat{T}},
                             B::StridedMatrix{DoubleFloat{T}},
                             _add::MulAddMul) where {T <: AbstractFloat}
    require_one_based_indexing(C, A, B)
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    if  mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA), matrix B has dimensions ($mB,$nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result matrix C dimensions $(size(C)), needs ($mA,$nB)"))
    end

    if iszero(_add.alpha) || isempty(A) || isempty(B)
        return _rmul_or_fill!(C, _add.beta)
    end

    if tA == 'N'
        At = Matrix(transpose(A))
    else
        At = A
    end

    Bline = zeros(DoubleFloat{T},mB)
    for j = 1:nB
        @inbounds begin
            if tB == 'N'
                @simd for k = 1:mB
                    Bline[k] = B[k,j]
                end
            else
                @simd for k = 1:mB
                    Bline[k] = B[j,k]
                end
            end
            for i=1:mA
                # benchmarking showed view to beat uview here
                asum = _dot(view(At,:,i),Bline,Vec{Npref,T})
                _modify!(_add, asum, C, (i,j))
            end
        end
    end
    C
end

# isolate the middle loops so that the inference engine has a fair chance

function gemm_kernN(C::StridedMatrix{DoubleFloat{T}}, At, B, Bline, j, mB, mA) where T
    @inbounds begin
        @simd for k = 1:mB
            Bline[k] = B[k,j]
        end
        for i=1:mA
            asum = _dot(uview(At,:,i),Bline,Vec{Npref,T})
            C[i,j] = asum
        end
    end
end

function gemm_kernT(C::StridedMatrix{DoubleFloat{T}}, At, B, Bline, j, mB, mA) where T
    @inbounds begin
        @simd for k = 1:mB
            Bline[k] = B[j,k]
        end
    end
    for i=1:mA
        asum = _dot(uview(At,:,i),Bline,Vec{Npref,T})
        C[i,j] = asum
    end
end

function _generic_matmatmul!(C::StridedMatrix{Complex{DoubleFloat{T}}},
                            tA::AbstractChar, tB::AbstractChar,
                            A::StridedMatrix{Complex{DoubleFloat{T}}},
                            B::StridedMatrix{Complex{DoubleFloat{T}}},
                            _add::MulAddMul=MulAddMul()
                            ) where {T <: AbstractFloat}
    require_one_based_indexing(C, A, B)
    mA, nA = lapack_size(tA, A)
    mB, nB = lapack_size(tB, B)
    if  mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA), matrix B has dimensions ($mB,$nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result matrix C dimensions $(size(C)), needs ($mA,$nB)"))
    end

    (tB ∈ ['N','C','T']) || throw(ArgumentError("unrecognized value of tB"))
    (tA ∈ ['N','C','T']) || throw(ArgumentError("unrecognized value of tA"))

    if iszero(_add.alpha) || isempty(A) || isempty(B)
        return _rmul_or_fill!(C, _add.beta)
    end

    reAt = Matrix{DoubleFloat{T}}(undef, nA, mA)
    imAt = Matrix{DoubleFloat{T}}(undef, nA, mA)
    if tA == 'N'
        @inbounds for j=1:nA
            @simd for i=1:mA
                z = A[i,j]
                reAt[j,i] = z.re
                imAt[j,i] = z.im
            end
        end
    elseif tA == 'C'
        @inbounds for i=1:mA
            @simd for j=1:nA
                z = A[j,i]
                reAt[j,i] = z.re
                imAt[j,i] = -z.im
            end
        end
    elseif tA == 'T'
        @inbounds for i=1:mA
            @simd for j=1:nA
                z = A[j,i]
                reAt[j,i] = z.re
                imAt[j,i] = z.im
            end
        end
    end
    if ((nthreads() > 1) && (Float64(mB)*Float64(mA) > gemm_mt_threshold[])
        && _add.alpha == 1 && _add.beta == 0)
        _mt_gemmwrap(C,mA,mB,nB,tB,reAt,imAt,B,Vec{Npref,T})
    else
        _gemmcore(C,mA,mB,nB,tB,reAt,imAt,B,Vec{Npref,T},_add)
    end
    C
end

function _gemmcore(C::AbstractMatrix{Complex{DoubleFloat{T}}},
                   mA,mB,nB,tB,reAt,imAt,B,::Type{Vec{N,T}},
                   _add::MulAddMul
                   ) where {N,T <: AbstractFloat}
    reBline = zeros(DoubleFloat{T},mB)
    imBline = zeros(DoubleFloat{T},mB)
    for j = 1:nB
        @inbounds begin
            if tB == 'N'
                @simd for k = 1:mB
                    z = B[k,j]
                    reBline[k] = z.re
                    imBline[k] = z.im
                end
            elseif tB == 'C'
                @simd for k = 1:mB
                    z = B[j,k]
                    reBline[k] = z.re
                    imBline[k] = -z.im
                end
            else
                @simd for k = 1:mB
                    z = B[j,k]
                    reBline[k] = z.re
                    imBline[k] = z.im
                end
            end
            for i=1:mA
                asum1 = _dot(uview(reAt,:,i),reBline,Vec{N,T})
                asum2 = _dot(uview(reAt,:,i),imBline,Vec{N,T})
                asum3 = _dot(uview(imAt,:,i),reBline,Vec{N,T})
                asum4 = _dot(uview(imAt,:,i),imBline,Vec{N,T})
                s = complex(asum1-asum4, asum2+asum3)
                _modify!(_add, s, C, (i,j))
            end
        end
    end
end

function _mt_gemmwrap(C::AbstractMatrix{Complex{DoubleFloat{T}}},
                      mA,mB,nB,tB,reAt,imAt,B,::Type{Vec{N,T}}
                      ) where {N, T <: AbstractFloat}
    nt = nthreads()
    reBlines = [zeros(DoubleFloat{T},mB) for id in 1:nt]
    imBlines = [zeros(DoubleFloat{T},mB) for id in 1:nt]
    if tB == 'N'
        @threads for j = 1:nB
            id = threadid()
            reBline = reBlines[id]
            imBline = imBlines[id]
            _mt_gemmcoreN(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                          Vec{N,T})
        end
    elseif tB == 'C'
        @threads for j = 1:nB
            id = threadid()
            reBline = reBlines[id]
            imBline = imBlines[id]
            _mt_gemmcoreC(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                          Vec{N,T})
        end
    elseif tB == 'T'
        @threads for j = 1:nB
            id = threadid()
            reBline = reBlines[id]
            imBline = imBlines[id]
            _mt_gemmcoreT(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                          Vec{N,T})
        end
    end
    nothing
end

function _mt_gemmcoreN(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                       ::Type{Vec{N,T}}) where {N,T}
    @inbounds begin
        @simd for k = 1:mB
            z = B[k,j]
            reBline[k] = z.re
            imBline[k] = z.im
        end
        for i=1:mA
            asum1 = _dot(uview(reAt,:,i),reBline,Vec{N,T})
            asum2 = _dot(uview(reAt,:,i),imBline,Vec{N,T})
            asum3 = _dot(uview(imAt,:,i),reBline,Vec{N,T})
            asum4 = _dot(uview(imAt,:,i),imBline,Vec{N,T})
            C[i,j] = complex(asum1-asum4, asum2+asum3)
        end
    end
    nothing
end

function _mt_gemmcoreT(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                       ::Type{Vec{N,T}}) where {N,T}
    @inbounds begin
        @simd for k = 1:mB
            z = B[j,k]
            reBline[k] = z.re
            imBline[k] = z.im
        end
        for i=1:mA
            asum1 = _dot(uview(reAt,:,i),reBline,Vec{N,T})
            asum2 = _dot(uview(reAt,:,i),imBline,Vec{N,T})
            asum3 = _dot(uview(imAt,:,i),reBline,Vec{N,T})
            asum4 = _dot(uview(imAt,:,i),imBline,Vec{N,T})
            C[i,j] = complex(asum1-asum4, asum2+asum3)
        end
    end
    nothing
end

function _mt_gemmcoreC(C, reAt, imAt, B, reBline, imBline, j, mB, mA,
                       ::Type{Vec{N,T}}) where {N,T}
    @inbounds begin
        @simd for k = 1:mB
            z = B[j,k]
            reBline[k] = z.re
            imBline[k] = -z.im
        end
        for i=1:mA
            asum1 = _dot(uview(reAt,:,i),reBline,Vec{N,T})
            asum2 = _dot(uview(reAt,:,i),imBline,Vec{N,T})
            asum3 = _dot(uview(imAt,:,i),reBline,Vec{N,T})
            asum4 = _dot(uview(imAt,:,i),imBline,Vec{N,T})
            C[i,j] = complex(asum1-asum4, asum2+asum3)
        end
    end
    nothing
end

max_stars_repo_name: RalphAS/DoubleBLAS

max_issues_repo_name: RalphAS/DoubleBLAS

max_forks_repo_name: RalphAS/DoubleBLAS