text:
function plmdca_asym(Z::Array{T,2},W::Vector{Float64};
                decimation::Bool=false,
                fracmax::Real=0.3,
                fracdec::Real=0.1,
                remove_dups::Bool=true,
                min_separation::Int=1,
                theta=:auto,
                lambdaJ::Real=0.01,
                lambdaH::Real=0.01,
                epsconv::Real=1.0e-5,
                maxit::Int=1000,
                verbose::Bool=true,
                method::Symbol=:LD_LBFGS) where T <: Integer

    all(x -> x > 0, W) || throw(DomainError("vector W should normalized and with all positive elements"))
    isapprox(sum(W), 1) || throw(DomainError("sum(W) ≠ 1. Consider normalizing the vector W"))
    N, M = size(Z)
    M = length(W)
    q = Int(maximum(Z))

    plmalg = PlmAlg(method, verbose, epsconv, maxit)
    plmvar = PlmVar(N, M, q, q * q, lambdaJ, lambdaH, Z, W)
    Jmat, pslike = if decimation  == false
        MinimizePLAsym(plmalg, plmvar)
    else
        decvar = DecVar{2}(fracdec, fracmax, ones(Bool, (N - 1) * q * q, N))
        DecimateAsym!(plmvar, plmalg, decvar)
    end
    score, FN, Jtensor, htensor, Jt1,Jt2 =  ComputeScore(Jmat, plmvar, min_separation)
    return PlmOut(sdata(pslike), Jtensor, htensor, score,Jmat,Jt1,Jt2)

end
plmdca(Z,W;kwds...) = plmdca_asym(Z, W;kwds...)

function plmdca_asym(filename::String;
                theta::Union{Symbol,Real}=:auto,
                max_gap_fraction::Real=0.9,
                remove_dups::Bool=true,
                kwds...)
    time = @elapsed W, Z, N, M, q = ReadFasta(filename, max_gap_fraction, theta, remove_dups)
    println("preprocessing took $time seconds")
    plmdca_asym(Z, W; kwds...)
end

plmdca(filename::String; kwds...) = plmdca_asym(filename; kwds...)

function MinimizePLAsym(alg::PlmAlg, var::PlmVar)

    LL = (var.N - 1) * var.q2 + var.q
    x0 = zeros(Float64, LL)
    vecps = SharedArray{Float64}(var.N)
    Jmat = zeros(LL,var.N) |> SharedArray
    Threads.@threads for site = 1:var.N # 1:12
        opt = Opt(alg.method, length(x0))
        ftol_abs!(opt, alg.epsconv)
        xtol_rel!(opt, alg.epsconv)
        xtol_abs!(opt, alg.epsconv)
        ftol_rel!(opt, alg.epsconv)
        maxeval!(opt, alg.maxit)
        min_objective!(opt, (x, g) -> optimfunwrapper(x, g, site, var))
        elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
        alg.verbose && @printf("site = %d\t pl = %.4f\t time = %.4f\t", site, minf, elapstime)
        alg.verbose && println("exit status = $ret")
        vecps[site] = minf
        Jmat[:,site] .= minx
    end
    # Jmat = @distributed hcat for site = 1:var.N # 1:12
    #     opt = Opt(alg.method, length(x0))
    #     ftol_abs!(opt, alg.epsconv)
    #     xtol_rel!(opt, alg.epsconv)
    #     xtol_abs!(opt, alg.epsconv)
    #     ftol_rel!(opt, alg.epsconv)
    #     maxeval!(opt, alg.maxit)
    #     min_objective!(opt, (x, g) -> optimfunwrapper(x, g, site, var))
    #     elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
    #     alg.verbose && @printf("site = %d\t pl = %.4f\t time = %.4f\t", site, minf, elapstime)
    #     alg.verbose && println("exit status = $ret")
    #     vecps[site] = minf
    #     minx
    # end
    return sdata(Jmat), vecps
end

function ComputeUL(alg::PlmAlg, var::PlmVar, site::Int, LL::Int)

    N  = var.N
    q2 = var.q2
    lb = -Inf * ones(Float64, LL)
    ub =  Inf * ones(Float64, LL)
    tiny::Float64 = 1.0e-6
    offset::Int = 0

    for i = 1:site - 1
        for s = 1:q2
            lb[offset + s] = -tiny
            ub[offset + s] =  tiny
        end
    end
    offset += q2

    for i = site + 1:N
        for s = 1:q2
            lb[offset + s] = -tiny
            ub[offset + s] =  tiny
        end
        offset += q2
    end
    return lb, ub
end

function PLsiteAndGrad!(x::Vector{Float64}, grad::Vector{Float64}, site::Int, plmvar::PlmVar)

    LL = length(x)
    q2 = plmvar.q2
    q = plmvar.q
    N = plmvar.N
    M = plmvar.M
    Z = sdata(plmvar.Z)
    W = sdata(plmvar.W)
    IdxZ = sdata(plmvar.IdxZ)

    for i = 1:LL - q
	    grad[i] = 2.0 * plmvar.lambdaJ  * x[i]
	end
	for i = (LL - q + 1):LL
	    grad[i] = 4.0 * plmvar.lambdaH * x[i]
	end

	pseudolike = 0.0
	vecene = zeros(Float64, q)
	lnorm = 0.0
	expvecenesumnorm = zeros(Float64, q)

    @inbounds for m = 1:M
        izm = view(IdxZ, :, m)
        zsm = Z[site,m]
		fillvecene!(vecene, x, site, izm, q, N)
	    lnorm = logsumexp(vecene)
	    expvecenesumnorm .= @. exp(vecene - lnorm)
	    pseudolike -= W[m] * (vecene[ zsm ] - lnorm)
        @avx for i = 1:site - 1
	        for s = 1:q
	            grad[ izm[i] + s ] += W[m] * expvecenesumnorm[s]
	        end
	        grad[ izm[i] + zsm ] -= W[m]
	    end
        @avx for i = site + 1:N
	        for s = 1:q
	            grad[ izm[i] - q2 + s ] += W[m] *  expvecenesumnorm[s]
	        end
	        grad[ izm[i] - q2 + zsm ] -= W[m]
	    end
        @avx for s = 1:q
	        grad[ (N - 1) * q2 + s ] += W[m] * expvecenesumnorm[s]
	    end
		grad[ (N - 1) * q2 + zsm ] -= W[m]
	end
	pseudolike += L2norm_asym(x, plmvar)
end

# Energy filling

function fillvecene!(vecene::Vector{Float64}, x::Vector{Float64}, site::Int, IdxSeq::AbstractArray{Int,1}, q::Int, N::Int)

	q2 = q * q
    @inbounds for l = 1:q
        scra::Float64 = 0.0
        @avx for i = 1:site - 1 # Begin sum_i \neq site J
            scra += x[IdxSeq[i] + l]
        end
        # skipping sum over residue site
        @avx for i = site + 1:N
            scra +=  x[IdxSeq[i] - q2 + l]
        end # End sum_i \neq site J
        scra +=  x[(N - 1) * q2 + l] # sum H
        vecene[l] = scra
    end

end

function logsumexp(X::Vector{Float64})
    u = maximum(X)
    isfinite(u) || return float(u)
    return u + log(sum(x -> exp(x - u), X))
end

function L2norm_asym(vec::Array{Float64,1}, plmvar::PlmVar)
    q = plmvar.q
    N = plmvar.N
    lambdaJ = plmvar.lambdaJ
    lambdaH = plmvar.lambdaH

    LL = length(vec)

    mysum1 = 0.0
    @inbounds @avx for i = 1:(LL - q)
        mysum1 += vec[i] * vec[i]
    end
    mysum1 *= lambdaJ

    mysum2 = 0.0
    @inbounds @avx for i = (LL - q + 1):LL
        mysum2 += vec[i] * vec[i]
    end
    mysum2 *= 2lambdaH

    return mysum1 + mysum2
end

function ComputeScore(Jmat::Array{Float64,2}, var::PlmVar, min_separation::Int)

    q = var.q
    N = var.N
    JJ = reshape(Jmat[1:end - q,:], q, q, N - 1, N)
    Jtemp1 = zeros(q, q, Int(N * (N - 1) / 2))
    Jtemp2 = zeros(q, q, Int(N * (N - 1) / 2))
    l = 1
    for i = 1:(N - 1)
        for j = (i + 1):N
            Jtemp1[:,:,l] = JJ[:,:,j - 1,i] # J_ij as estimated from from g_i.
            Jtemp2[:,:,l] = JJ[:,:,i,j]' # J_ij as estimated from from g_j.
            l = l + 1
        end
    end

    hplm = fill(0.0, q, N)
    for i in 1:N
        hplm[:,i] .= Jmat[end - q + 1:end,i]
    end

    Jtensor1 = inflate_matrix(Jtemp1, N)
    Jtensor2 = inflate_matrix(Jtemp2, N)
    Jplm = (Jtensor1 + Jtensor2) / 2 # for the energy I do not want to gauge

    ctr = 0
    for i in 1:N - 1
        for j in i + 1:N
            ctr += 1
            Jtensor1[:,:,i,j] = Jtemp1[:,:,ctr] - repeat(mean(Jtemp1[:,:,ctr], dims=1), q, 1) - repeat(mean(Jtemp1[:,:,ctr], dims=2), 1, q) .+ mean(Jtemp1[:,:,ctr])
            Jtensor1[:,:,j,i] = Jtensor1[:,:,i,j]'
            Jtensor2[:,:,i,j] = Jtemp2[:,:,ctr] - repeat(mean(Jtemp2[:,:,ctr], dims=1), q, 1) - repeat(mean(Jtemp2[:,:,ctr], dims=2), 1, q) .+ mean(Jtemp2[:,:,ctr])
            Jtensor2[:,:,j,i] = Jtensor2[:,:,i,j]'
        end
    end # zerosumgauge the different tensors

    Jtensor = (Jtensor1 + Jtensor2) / 2

    FN = compute_APC(Jtensor, N, q)
    score = compute_ranking(FN, min_separation)
    return score, FN, Jplm, hplm, Jtensor1,Jtensor2
end

max_stars_repo_name: giovact/PlmDCA

max_issues_repo_name: giovact/PlmDCA

max_forks_repo_name: giovact/PlmDCA


text:
using DifferentialEvolution
using Test

@testset "DifferentialEvolution.jl" begin
    include("nsga_ii_test.jl")
    include("constraint_test.jl")
end

ext: jl


text:
@testset "CookieCutter" begin
  problem = SimulationProblem(CartesianGrid(100,100), (:facies => Int, :property => Float64), 3)

  solver = CookieCutter(DummySimSolver(:facies=>NamedTuple()), Dict(0=>DummySimSolver(), 1=>DummySimSolver()))

  @test sprint(show, solver) == "CookieCutter"
  @test sprint(show, MIME"text/plain"(), solver) == "CookieCutter\n  └─facies ⇨ DummySimSolver\n    └─0 ⇨ DummySimSolver\n    └─1 ⇨ DummySimSolver"

  Random.seed!(1234)
  solution = solve(problem, solver)

  if visualtests
    @test_reference "data/COOKIE.png" plot(solution)
  end
end


text:
using Base.Test
import NeuralNetworks
nn = NeuralNetworks

m = nn.SpatialConv(Float32,3,16,5,5,1,1)
input = rand(Float32, 3, 32, 32)
output = nn.forward!(m, input)

@test 1 == 1


text:
"""
        critic(decisionMat, fns)

Apply CRITIC (Combined Compromise Solution) method for a given matrix and criteria types.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
critic() applies the CRITIC method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::CRITICResult`: CRITICResult object that holds multiple outputs including weighting and best index.

# Examples
```julia-repl

julia> decmat
3×4 Array{Float64,2}:
 12.9918  0.7264  -1.1009  1.59814
  4.1201  5.8824   3.4483  1.02156
  4.1039  0.0     -0.5076  0.984469

julia> df = makeDecisionMatrix(decmat)

3×4 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4
     │ Float64  Float64  Float64  Float64
─────┼─────────────────────────────────────
   1 │ 12.9918   0.7264  -1.1009  1.59814
   2 │  4.1201   5.8824   3.4483  1.02156
   3 │  4.1039   0.0     -0.5076  0.984469

julia> fns = [maximum, maximum, minimum, maximum];

julia> result = critic(df, fns);

julia> result.w
4-element Array{Float64,1}:
 0.16883905506169491
 0.41844653698732126
 0.24912338769165807
 0.16359102025932576

julia> result.bestIndex
2
```
# References

Diakoulaki, D., Mavrotas, G., & Papayannakis, L. (1995). Determining objective weights in multiple criteria problems: The critic method. Computers & Operations Research, 22(7), 763–770. doi:10.1016/0305-0548(94)00059-h
Akçakanat, Ö., Aksoy, E., Teker, T. (2018). CRITIC ve MDL Temelli EDAS Yöntemi ile TR-61 Bölgesi Bankalarının Performans Değerlendirmesi. Süleyman Demirel Üniversitesi Sosyal Bilimler Enstitüsü Dergisi, 1 (32), 1-24.

"""
function critic(decisionMat::DataFrame, fns::Array{Function,1})::CRITICResult

    row, col = size(decisionMat)
    colMax = colmaxs(decisionMat)
    colMin = colmins(decisionMat)

    A = similar(decisionMat)

    for i in 1:row
        for j in 1:col
            if fns[j] == maximum
                @inbounds A[i, j] = (decisionMat[i, j] - colMin[j]) / (colMax[j] - colMin[j])
            elseif fns[j] == minimum
                @inbounds A[i, j] = (colMax[j] - decisionMat[i, j]) / (colMax[j] - colMin[j])
            end
        end
    end

    # normalizedMat = convert(Matrix, A)
    normalizedMat = Matrix(A)

    corMat = 1 .- cor(normalizedMat)

    scores = zeros(Float64, col)
    for i in 1:col
        scores[i] = sum(corMat[:, i]) .* std(normalizedMat[:, i])
    end

    w = zeros(Float64, col)

    for i in 1:col
        w[i] = scores[i] ./ sum(scores)
    end

    rankings = sortperm(w)

    bestIndex = rankings |> last

    result = CRITICResult(
        decisionMat,
        w,
        rankings,
        bestIndex
    )

    return result
end

"""
        critic(setting)

Apply CRITIC (Combined Compromise Solution) method for a given matrix and criteria types.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

# Description
critic() applies the CRITIC method to rank n alterntives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::CRITICResult`: CRITICResult object that holds multiple outputs including weighting and best index.
"""
function critic(setting::MCDMSetting)::CRITICResult
    critic(
        setting.df,
        setting.fns
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
using GLVisualize, GLWindow, Colors, IndirectArrays

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
Iso surface volume Plot.
"""

function volume_data(N)
    id = zeros(UInt32, N, N, N)
    id[:, :, 1:N÷2] = 1
    id[:, :, N÷2+1:N] = 2
    vol = IndirectArray(id, [RGBA{N0f8}(1,0,0,1), RGBA{N0f8}(0,0,1,1)])
end

vol = visualize(volume_data(128), :absorption)
_view(vol, window)

if !isdefined(:runtests)
    renderloop(window)
end


text:
export BinaryDistillationColumn

"""
    BinaryDistillationColumn

Type representing the IFAC Binary Distillation Column benchmark problem.

From E. J. Davison, ‘Benchmark Problems for Control System Design: Report of the IFAC Theory Committee’, IFAC, 1990.
"""
struct BinaryDistillationColumn <: ControlBenchmarkProblem end

@addbenchmark BinaryDistillationColumn "ifac01" "Binary distillation column" "A distillation column to separate two compounds" Set([:linear, :stable, :continuous, :process, :ifac])

"""
    controlbenchmark(::BinaryDistillationColumn) -> (sys, input, disturbance)

Generate the components of the binary distillation column benchmark problem.

`sys` is the linear statespace system describing the distillation column.

`input` is the set of `Bounds` representing the constraints on the input actuators.

`disturbance` is the bound representing the constraint on the disturbance input (the input feed).
"""
function controlbenchmark( ::BinaryDistillationColumn )
    # The state transition matrix
    #           R   C    Value
    Aparts = [  1   1  -1.4e-2;
                1   2   4.3e-3;
                2   1   9.5e-3;
                2   2  -1.38e-2;
                2   3   4.6e-3;
                2  11   5.0e-4;
                3   2   9.5e-3;
                3   3  -1.41e-2;
                3   4   6.3e-3;
                3  11   2.0e-4;
                4   3   9.5e-3;
                4   4  -1.58e-2;
                4   5   1.1e-2;
                5   4   9.5e-3;
                5   5  -3.12e-2;
                5   6   1.5e-2;
                6   5   2.02e-4;
                6   6  -3.52e-2;
                6   7   2.2e-2;
                7   6   2.02e-2;
                7   7  -4.22e-2
                7   8   2.8e-2;
                8   7   2.02e-2;
                8   8  -4.82e-2;
                8   9   3.7e-2;
                8  11   2.0e-4;
                9   8   2.02e-2;
                9   9  -5.7e-2;
                9  10   4.2e-2;
                9  11   5.0e-4;
               10   9   2.02e-2;
               10  10  -4.83e-2;
               10  11   5.0e-4;
               11   1   2.55e-2;
               11  10   2.55e-2;
               11  11  -1.85e-2 ]

    A = sparse( Aparts[:, 1], Aparts[:, 2], Aparts[:, 3], 11, 11 )

    # The input matrix
    Bdense = [  0.0      0.0     0.0;
                5.0e-6  -4.0e-5  2.5e-3;
                2.0e-6  -2.0e-5  5.0e-3;
                1.0e-6  -1.0e-5  5.0e-3;
                0.0      0.0     5.0e-3;
                0.0      0.0     5.0e-3;
               -5.0e-6   1.0e-5  5.0e-3;
               -1.0e-5   3.0e-5  5.0e-3;
               -4.0e-5   5.0e-6  2.5e-3;
               -2.0e-5   2.0e-6  2.5e-3;
                4.6e-4   4.6e-4  0.0 ]

    B = sparse( Bdense )

    # The output matrix
    C = spzeros( 3, 11 )
    C[1, 10] = 1.0
    C[2,  1] = 1.0
    C[3, 11] = 1.0

    # The feed-forward matrix
    D = spzeros( 3, 3 )

    # The disturbance matrix
    # TODO: Extend ss to take this disturbance matrix
    E = spzeros( 11, 1 )
    E[5, 1] = 1.0e-2

    # Input constraints
    inputBound = Bounds( [ Bound( -2.5,  2.5 );
                           Bound( -2.5,  2.5 );
                           Bound( -0.30, 0.30) ] )

    # Disturbance bound for FL
    disturbanceBound = Bound( -1.0, 1.0 )

    # Create the constrained system
    return (; :sys => HeteroStateSpace( A, B, C, D, Continuous() ),
              :inputBounds => inputBound,
              :disturbanceBounds => disturbanceBound )
end


text:
# Autogenerated wrapper script for XGBoost_jll for i686-linux-gnu-cxx03
export libxgboost, xgboost

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("XGBoost")
JLLWrappers.@declare_library_product(libxgboost, "libxgboost.so")
JLLWrappers.@declare_executable_product(xgboost)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libxgboost,
        "lib/libxgboost.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        xgboost,
        "bin/xgboost",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export Solver, SolverParameters, SolverState, LRPolicy, MomPolicy, solve, SGD,
    get_learning_rate, get_momentum, update
abstract Solver

type SolverState
    iter :: Int
    obj_val :: Float32
    learning_rate :: Float32
    momentum :: Float32
    accuracy_log :: IOStream
    loss_log :: IOStream
    SolverState(iter::Int, obj_val::Float32, learning_rate::Float32, momentum::Float32) = new(iter, obj_val, learning_rate, momentum)
end

abstract LearningRatePolicy
module LRPolicy
using ..Latte.LearningRatePolicy
type Fixed <: LearningRatePolicy
    base_lr :: Float32
end

# base_lr * gamma ^ (floor(iter / stepsize))
type Step <: LearningRatePolicy
    base_lr  :: Float32
    gamma    :: Float32
    stepsize :: Int
end

function Step(;base_lr=0.01f0, gamma=0.1f0, stepsize=100000)
    Step(base_lr, gamma, stepsize)
end

# base_lr * gamma ^ iter
type Exp <: LearningRatePolicy
    base_lr :: Float32
    gamma   :: Float32
end

type Inv <: LearningRatePolicy
    base_lr :: Float32
    gamma   :: Float32
    power   :: Float32
end

type Decay <: LearningRatePolicy
    base_lr :: Float32
    decay   :: Float32
end

type Poly <: LearningRatePolicy
    base_lr  :: Float32
    max_iter :: Float32
    power    :: Float32
end
end # module LRPolicy

get_learning_rate(policy::LRPolicy.Fixed, state::SolverState) = policy.base_lr
get_learning_rate(policy::LRPolicy.Step, state::SolverState) =
    policy.base_lr * policy.gamma ^ (floor(state.iter / policy.stepsize))
get_learning_rate(policy::LRPolicy.Exp, state::SolverState) =
    policy.base_lr * policy.gamma ^ state.iter
get_learning_rate(policy::LRPolicy.Inv, state::SolverState) =
    policy.base_lr * (1 + policy.gamma * state.iter) ^ (-policy.power)
get_learning_rate(policy::LRPolicy.Decay, state::SolverState) =
    policy.base_lr / (1 + state.iter * policy.decay)
get_learning_rate(policy::LRPolicy.Poly, state::SolverState) =
    policy.base_lr * (1 - state.iter / policy.max_iter) ^ (policy.power)

abstract MomentumPolicy
module MomPolicy
using ..Latte.MomentumPolicy
type Fixed <: MomentumPolicy
  base_mom :: Float32
end

# min(base_mom * gamma ^ (floor(iter / stepsize)), max_mom)
type Step <: MomentumPolicy
    base_mom :: Float32
    gamma    :: Float32
    stepsize :: Int
    max_mom  :: Float32
end

type Linear <: MomentumPolicy
    base_mom :: Float32
    gamma    :: Float32
    stepsize :: Int
    max_mom  :: Float32
end
end # module MomPolicy

get_momentum(policy::MomPolicy.Fixed, state::SolverState) = policy.base_mom
get_momentum(policy::MomPolicy.Step, state::SolverState) =
    min(policy.base_mom * policy.gamma ^ (floor(state.iter / policy.stepsize)), policy.max_mom)
get_momentum(policy::MomPolicy.Linear, state::SolverState) =
    min(policy.base_mom + floor(state.iter / policy.stepsize) * policy.gamma, policy.max_mom)

type SolverParameters
    lr_policy::LearningRatePolicy
    mom_policy::MomentumPolicy
    max_epoch::Int
    regu_coef::Float32
    snapshot_dir::AbstractString
end

function SolverParameters(;
        lr_policy=LR.Decay(.01f0, 5.0f-7),
        mom_policy=MomPolicy.Fixed(0.9),
        max_epoch=300,
        regu_coef=.0005,
        snapshot_dir="")
    if snapshot_dir == ""
        snapshot_dir = string(now())
    end
    SolverParameters(lr_policy, mom_policy, max_epoch, regu_coef, snapshot_dir)
end

type SGD <: Solver
    params :: SolverParameters
    state  :: SolverState
    SGD(params::SolverParameters) = new(params, SolverState(0, 0.0f0, 0.0f0, 0.0f0))
end

function update(solver::Solver, net::Net, param_id::UInt64)
    for param in net.params
        if object_id(param) == param_id
            update(solver, param)
            break
        end
    end
end

function update(sgd::SGD, net::Net)
    for param in net.params
        l2_regularization(sgd.params.regu_coef * param.regu_coef, param.value, param.gradient)
        sgd_update(sgd.state.learning_rate * param.learning_rate,
                   sgd.state.momentum, param.value, param.gradient, param.hist)
    end
end

function update(sgd::SGD, param::Param)
    @latte_mpi(@eval(ccall((:wait, $libComm), Void, (Cint,), $(param.request))))
    if !LOSSY_GRADIENTS
        gradient = sum(param.gradient, ndims(param.gradient))
    else
        gradient = param.gradient
    end
    l2_regularization(sgd.params.regu_coef * param.regu_coef, param.value, gradient)
    sgd_update(sgd.state.learning_rate * param.learning_rate,
               sgd.state.momentum, param.value, gradient, param.hist)
end

function sgd_update{T}(learning_rate::Float32, momentum::Float32,
                       param::Array{T}, gradient::Array{T}, hist::Array{T})
    momentum = convert(T, momentum)
    learning_rate = convert(T, learning_rate)
    # We use BLAS here because the array notation in the comments does not work
    # with current julia, it creates a copy of param and does not increment the
    # array in place

    # hist *= momentum
    BLAS.scal!(length(hist), momentum, pointer(hist), 1)
    BLAS.axpy!(length(hist), learning_rate, pointer(gradient), 1,
               pointer(hist), 1)
    # param -= hist
    BLAS.axpy!(length(hist), convert(T, -1), pointer(hist), 1, pointer(param), 1)
end

function clip_gradients(sgd::SGD, net::Net)
    for param in net.params
        if param.clip_gradients < 0.0f0
            return
        end
    end
    sumsq = 0.0
    for param in net.params
        sumsq += dot(param.gradient[:], param.gradient[:])
    end
    l2norm_diff = sqrt(sumsq)
    if l2norm_diff > param.clip_gradients
        scale = param.clip_gradients / l2norm_diff
        for param in net.params
            BLAS.scal!(length(param.gradient), scale, pointer(param.gradient), 1)
        end
    end
end

function regularize(sgd::SGD, net::Net)
    for param in net.params
        l2_regularization(sgd.params.regu_coef * param.regu_coef, param.value, param.gradient)
    end
end

function l2_regularization(regu_coef::Float32, param::Array, gradient::Array)
    BLAS.axpy!(length(param), convert(eltype(param), 2.0 * regu_coef), pointer(param), 1, pointer(gradient), 1)
end

function cleanup(solver::Solver)
    close(solver.state.accuracy_log)
    close(solver.state.loss_log)
end

function solve(solver::Solver, net::Net)
    init(net)

    solver.state.learning_rate = get_learning_rate(solver.params.lr_policy,
                                                   solver.state)
    solver.state.momentum = get_momentum(solver.params.mom_policy,
                                         solver.state)
    @latte_mpi broadcast_initial_params(net)

    @latte_mpi(if get_inter_rank(net) == 0 && get_net_subrank(net) + 1 == net.num_subgroups
        if isdir(solver.params.snapshot_dir)
            log_info("Snapshot directory exists, overwriting.")
            rm(solver.params.snapshot_dir; recursive=true)
        end
        mkdir(solver.params.snapshot_dir)
        solver.state.accuracy_log = open("$(solver.params.snapshot_dir)/accuracy.csv", "w")
        solver.state.loss_log = open("$(solver.params.snapshot_dir)/loss.csv", "w")
        atexit(() -> cleanup(solver))
    end
    , begin
        if isdir(solver.params.snapshot_dir)
            log_info("Snapshot directory exists, overwriting.")
            rm(solver.params.snapshot_dir; recursive=true)
        end
        mkdir(solver.params.snapshot_dir)
        solver.state.accuracy_log = open("$(solver.params.snapshot_dir)/accuracy.csv", "w")
        solver.state.loss_log = open("$(solver.params.snapshot_dir)/loss.csv", "w")
    end)
    log_info("Entering solve loop")
    curr_train_epoch = net.train_epoch
    while curr_train_epoch < solver.params.max_epoch
        # if LATTE_BATCH_DROPOUT
            rand_values(net)
        # end
        solver.state.iter += 1
        forward(net; solver=solver)
        clear_∇(net)
        backward(net)

        solver.state.obj_val = get_loss(net)
        solver.state.learning_rate = get_learning_rate(solver.params.lr_policy, solver.state)
        solver.state.momentum = get_momentum(solver.params.mom_policy, solver.state)

        clear_values(net)
        if solver.state.iter % 20 == 0
            @latte_mpi(if get_net_subrank(net) + 1 == net.num_subgroups
                log_info("Iter $(solver.state.iter) - Loss: $(solver.state.obj_val)")
                if get_inter_rank(net) == 0 && get_net_subrank(net) + 1 == net.num_subgroups
                    write(solver.state.loss_log, "$(solver.state.iter),$(solver.state.obj_val)\n")
                end
            end, begin
                log_info("Iter $(solver.state.iter) - Loss: $(solver.state.obj_val)")
                write(solver.state.loss_log, "$(solver.state.iter),$(solver.state.obj_val)\n")
            end)
        end
        @latte_mpi if net.num_subgroups > 1
            sync_intra_train_epoch(net)
        end
        if curr_train_epoch != net.train_epoch
            log_info("Epoch $(curr_train_epoch) - Testing...")
            acc = test(net)
            @latte_mpi(if get_net_subrank(net) + 1 == net.num_subgroups
                total_acc = @eval ccall((:reduce_accuracy, $libComm), Cfloat, (Cfloat,), $acc)
                if total_acc >= 0.0f0 && get_inter_rank(net) == 0
                    log_info("Epoch $(curr_train_epoch) - Test Result: $total_acc%")
                    write(solver.state.accuracy_log, "$(curr_train_epoch),$total_acc\n")
                end
            end, begin
                write(solver.state.accuracy_log, "$(curr_train_epoch),$acc\n")
                log_info("Epoch $(curr_train_epoch) - Test Result: $acc%")
            end)
            curr_train_epoch = net.train_epoch
        end
    end
end


text:
using Flux
using Flux: onehot, chunk, batchseq, throttle, crossentropy
using StatsBase: wsample
using Base.Iterators: partition
using CuArrays
using CUDAnative: device!
using Random

ϵ = 1.0f-32

epochs = 2
batch_size = 50
sequence = 50
gpu_device = 0

device!(gpu_device)
CuArrays.allowscalar(false)

input_file = joinpath(dirname(@__FILE__),"input.txt")

isfile(input_file) ||
    download("https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt",
             input_file)

text = collect(String(read(input_file)))
alphabet = [unique(text)...,'_']
text = map(ch -> Float32.(onehot(ch,alphabet)),text)
stop = Float32.(onehot('_',alphabet))

N = length(alphabet)
seqlen = sequence
nbatch = batch_size

Xs = collect(partition(batchseq(chunk(text, nbatch), stop), seqlen))
txt = circshift(text,-1)
txt[end] = stop
Ys = collect(partition(batchseq(chunk(txt, nbatch), stop), seqlen))

model = Chain(
  LSTM(N, 128),
  LSTM(128, 256),
  LSTM(256, 128),
  Dense(128, N),
  softmax)
  m = model |>gpu

opt = ADAM(0.01)
tx, ty = (Xs[5]|>gpu, Ys[5]|>gpu)

function loss(xx, yy)
  out = 0.0f0
  for (idx, x) in enumerate(xx)
    out += crossentropy(m(x) .+ ϵ, yy[idx])
  end
  Flux.reset!(m)
  out
end

idxs = length(Xs)
for epoch_idx in 1:epochs
  for (idx,(xs,ys)) in enumerate(zip(Xs, Ys))
    Flux.train!(loss, params(m), [(xs|>gpu,ys|>gpu)], opt)
    lss = loss(tx,ty)
    if idx % 10 == 0
      @info "epoch# $(epoch_idx)/$(epochs)-$(idx)/$(idxs) loss = $(lss)"
    end
  end
end

# Sampling
function sample(m, alphabet, len)
  m = cpu(m)
  Flux.reset!(m)
  buf = IOBuffer()
  c = rand(alphabet)
  for i = 1:len
    write(buf, c)
    c = wsample(alphabet, m(onehot(c, alphabet)))
  end
  return String(take!(buf))
end
@info sample(m, alphabet, 1000)

max_stars_repo_name: mrchaos/model-zoo

max_issues_repo_name: mrchaos/model-zoo

max_forks_repo_name: mrchaos/model-zoo


text:
function cuMPS(psi::MPS)
    phi = copy(psi)
    for site in 1:length(psi)
        phi.data[site] = cuITensor(psi.data[site])
    end
    return phi
end

cu(ψ::MPS) = cuMPS(ψ)

cuMPS() = MPS()
function cuMPS(::Type{T}, sites::Vector{<:Index}; linkdims::Integer=1) where {T<:Number}
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  if N == 1
    v[1] = emptyITensor(T, sites[1])
    return MPS(v)
  end

  space = if hasqns(sites)
    [QN() => linkdims]
  else
    linkdims
  end

  l = [Index(1, "Link,l=$ii") for ii=1:N-1]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = cuITensor(l[ii], s)
    elseif ii == N
      v[ii] = cuITensor(l[ii-1], s)
    else
      v[ii] = cuITensor(l[ii-1],s,l[ii])
    end
  end
  return MPS(v,0,N+1)
end
cuMPS(sites::Vector{<:Index}, args...; kwargs...) = cuMPS(Float64, sites, args...; kwargs...)
function cuMPS(N::Int; ortho_lims::UnitRange=1:N)
    return MPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims)
end

function randomCuMPS(sites)
  M = cuMPS(sites)
  for i in eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  M.llim = 1
  M.rlim = length(M)
  return M
end
function randomCuMPS(N::Int; ortho_lims::UnitRange=1:N)
  return randomCuMPS(Vector{ITensor}(undef, N); ortho_lims=ortho_lims)
end

const productCuMPS = cuMPS

function cuMPS(::Type{T}, ivals::Vector{<:Pair{<:Index}}) where {T<:Number}
  N     = length(ivals)
  As    = Vector{ITensor}(undef,N)
  links = Vector{Index}(undef,N)
  for n=1:N
    s = ITensors.ind(ivals[n])
    links[n] = Index(1,"Link,l=$n")
    if n == 1
      A = ITensor(T, s, links[n])
      A[ivals[n],links[n](1)] = 1.0
    elseif n == N
      A = ITensor(T, links[n-1], s)
      A[links[n-1](1),ivals[n]] = 1.0
    else
      A = ITensor(T, links[n-1], s, links[n])
      A[links[n-1](1),ivals[n],links[n](1)] = 1.0
    end
    As[n] = cuITensor(A)
  end
  return MPS(As,0,2)
end
cuMPS(ivals::Vector{<:IndexVal}) = cuMPS(Float64, ivals)

function cuMPS(::Type{T}, sites::Vector{<:Index}, states_) where {T<:Number}
  if length(sites) != length(states_)
    throw(DimensionMismatch("Number of sites and and initial vals don't match"))
  end
  N = length(states_)
  M = cuMPS(N)

  if N == 1
    M[1] = state(sites[1], states_[1])
    return M
  end

  states = [state(sites[j], states_[j]) for j in 1:N]

  if hasqns(states[1])
    lflux = QN()
    for j in 1:(N - 1)
      lflux += flux(states[j])
    end
    links = Vector{QNIndex}(undef, N - 1)
    for j in (N - 1):-1:1
      links[j] = Index(lflux => 1; tags="Link,l=$j", dir=In)
      lflux -= flux(states[j])
    end
  else
    links = [Index(1; tags="Link,l=$n") for n in 1:N]
  end

  M[1] = cuITensor(T, sites[1], links[1])
  M[1] += cuITensor(states[1] * state(links[1], 1))
  for n in 2:(N - 1)
    M[n] = cuITensor(T, dag(links[n - 1]), sites[n], links[n])
    M[n] += cuITensor(state(dag(links[n - 1]), 1) * states[n] * state(links[n], 1))
  end
  M[N] = cuITensor(T, dag(links[N - 1]), sites[N])
  M[N] += cuITensor(state(dag(links[N - 1]), 1) * states[N])

  return M
end

function cuMPS(
  ::Type{T}, sites::Vector{<:Index}, state::Union{String,Integer}
) where {T<:Number}
  return cuMPS(T, sites, fill(state, length(sites)))
end

function cuMPS(::Type{T}, sites::Vector{<:Index}, states::Function) where {T<:Number}
  states_vec = [states(n) for n in 1:length(sites)]
  return cuMPS(T, sites, states_vec)
end

cuMPS(sites::Vector{<:Index}, states) = cuMPS(Float64, sites, states)


text:
@testset "$TEST $G" begin

triangle = random_regular_graph(3,2)
quadrangle = random_regular_graph(4,2)
pentagon = random_regular_graph(5,2)

d, λ = edit_distance(triangle, quadrangle, subst_cost=MinkowskiCost(1:3,1:4))
@test d == 1.0
@test λ == Tuple[(1,1),(2,2),(3,3),(0,4)]

d, λ = edit_distance(quadrangle, triangle, subst_cost=MinkowskiCost(1:4,1:3))
@test d == 1.0
@test λ == Tuple[(1,1),(2,2),(3,3),(4,0)]

d, λ = edit_distance(triangle, pentagon, subst_cost=MinkowskiCost(1:3,1:5))
@test d == 2.0
@test λ == Tuple[(1,1),(2,2),(3,3),(0,4),(0,5)]

d, λ = edit_distance(pentagon, triangle, subst_cost=MinkowskiCost(1:5,1:3))
@test d == 2.0
@test λ == Tuple[(1,1),(2,2),(3,3),(4,0),(5,0)]

cost = MinkowskiCost(1:3,1:3)
bcost = BoundedMinkowskiCost(1:3,1:3)
for i=1:3
  @test cost(i,i) == 0.
  @test bcost(i,i) == 2/3
end

end # testset


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Subnet <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function Subnet(;id=nothing)
        o = new()
        validate_property(Subnet, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type Subnet

const _property_map_Subnet = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_Subnet = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ Subnet }) = collect(keys(_property_map_Subnet))
Swagger.property_type(::Type{ Subnet }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_Subnet[name]))}
Swagger.field_name(::Type{ Subnet }, property_name::Symbol) =  _property_map_Subnet[property_name]

function check_required(o::Subnet)
    true
end

function validate_property(::Type{ Subnet }, name::Symbol, val)
end


text:
"""
    federated_averaging(round_weights)

Implement the federated averaging algorithm descibed in:
- [Communication-Efficient Learning of Deep Networks from Decentralized Data](https://arxiv.org/abs/1602.05629)

`round_weights` must be an array containing the weights from each clients
participating in the round.
"""
function federated_averaging(round_weights)
    mapreduce(w -> w[:], +, round_weights) / length(round_weights)
end


text:
module PointGroups

using GenericSVD: svd
using NearestNeighbors: KDTree, nn
using StaticArrays: SArray, SVector, cross, det, norm, qr, svd
using Suppressor: @suppress

export cyclic_group, rotoreflection_group, pinwheel_group, pyramid_group,
    dihedral_group, prismatic_group, antiprismatic_group,
    chiral_tetrahedral_group, full_tetrahedral_group, pyritohedral_group,
    chiral_octahedral_group, full_octahedral_group,
    chiral_icosahedral_group, full_icosahedral_group,
    tetrahedron_vertices, tetrahedron_edge_centers, tetrahedron_face_centers,
    tetrahedron_rotoinversion_centers, pyritohedron_vertices,
    octahedron_vertices, octahedron_edge_centers, octahedron_face_centers,
    icosahedron_vertices, icosahedron_edge_centers, icosahedron_face_centers,
    CHIRAL_POLYHEDRAL_POINT_GROUPS, POLYHEDRAL_POINT_GROUPS,
    multiplication_table, count_central_elements, degenerate_orbits,
    isometric, isometries, identify_point_group

######################################################### 3D SYMMETRY OPERATIONS

function rotation_matrix_z(::Type{T}, i::Int, n::Int) where {T}
    _zero = zero(T)
    _one = one(T)
    if i % n == 0
        return SArray{Tuple{3,3},T,2,9}(
            +_one, _zero, _zero,
            _zero, +_one, _zero,
            _zero, _zero, +_one)
    elseif iseven(n) && (i % n) == (n >> 1)
        return SArray{Tuple{3,3},T,2,9}(
            -_one, _zero, _zero,
            _zero, -_one, _zero,
            _zero, _zero, +_one)
    end
    c = cospi(T(2*i) / T(n))
    s = sinpi(T(2*i) / T(n))
    return SArray{Tuple{3,3},T,2,9}(
        c, s, _zero, -s, c, _zero, _zero, _zero, _one)
end

function rotoreflection_matrix_z(::Type{T}, i::Int, n::Int) where {T}
    _zero = zero(T)
    _one = one(T)
    if i % n == 0
        return SArray{Tuple{3,3},T,2,9}(
            +_one, _zero, _zero,
            _zero, +_one, _zero,
            _zero, _zero, -_one)
    elseif iseven(n) && (i % n) == (n >> 1)
        return SArray{Tuple{3,3},T,2,9}(
            -_one, _zero, _zero,
            _zero, -_one, _zero,
            _zero, _zero, -_one)
    end
    c = cospi(T(2*i) / T(n))
    s = sinpi(T(2*i) / T(n))
    return SArray{Tuple{3,3},T,2,9}(
        c, s, _zero, -s, c, _zero, _zero, _zero, -_one)
end

function rotation_matrix_to_x(v::SVector{3,T}) where {T}
    _one = one(T)
    x, y, z = v
    if signbit(x)
        x, y, z = -x, -y, -z
    end
    xp1 = x + _one
    y2 = y * y
    yz = y * z
    z2 = z * z
    return SArray{Tuple{3,3},T,2,9}(
        x, -y, -z,
        y, _one - y2 / xp1, -yz / xp1,
        z, -yz / xp1, _one - z2 / xp1)
end

function rotation_matrix_to_y(v::SVector{3,T}) where {T}
    _one = one(T)
    x, y, z = v
    if signbit(y)
        x, y, z = -x, -y, -z
    end
    yp1 = y + _one
    x2 = x * x
    xz = x * z
    z2 = z * z
    return SArray{Tuple{3,3},T,2,9}(
        _one - x2 / yp1, x, -xz / yp1,
        -x, y, -z,
        -xz / yp1, z, _one - z2 / yp1)
end

function rotation_matrix_to_z(v::SVector{3,T}) where {T}
    _one = one(T)
    x, y, z = v
    if signbit(z)
        x, y, z = -x, -y, -z
    end
    zp1 = z + _one
    x2 = x * x
    xy = x * y
    y2 = y * y
    return SArray{Tuple{3,3},T,2,9}(
        _one - x2 / zp1, -xy / zp1, x,
        -xy / zp1, _one - y2 / zp1, y,
        -x, -y, z)
end

function rotation_matrix_x_pi(::Type{T}) where {T}
    _zero = zero(T)
    _one = one(T)
    return SArray{Tuple{3,3},T,2,9}(
        +_one, _zero, _zero,
        _zero, -_one, _zero,
        _zero, _zero, -_one)
end

function reflection_matrix_y(::Type{T}) where {T}
    _zero = zero(T)
    _one = one(T)
    return SArray{Tuple{3,3},T,2,9}(
        +_one, _zero, _zero,
        _zero, -_one, _zero,
        _zero, _zero, +_one)
end

############################################################# AXIAL POINT GROUPS

cyclic_group(::Type{T}, n::Int) where {T} =
    [rotation_matrix_z(T, i, n) for i = 0 : n-1]

function rotoreflection_group(::Type{T}, n::Int) where {T}
    @assert iseven(n)
    return [(iseven(i)
        ? rotation_matrix_z(T, i, n)
        : rotoreflection_matrix_z(T, i, n))
        for i = 0 : n-1]
end

pinwheel_group(::Type{T}, n::Int) where {T} = vcat(
    [rotation_matrix_z(T, i, n) for i = 0 : n-1],
    [rotoreflection_matrix_z(T, i, n) for i = 0 : n-1])

pyramid_group(::Type{T}, n::Int) where {T} = vcat(
    [rotation_matrix_z(T, i, n) for i = 0 : n-1],
    [rotation_matrix_z(T, i, n) * reflection_matrix_y(T) for i = 0 : n-1])

dihedral_group(::Type{T}, n::Int) where {T} = vcat(
    [rotation_matrix_z(T, i, n) for i = 0 : n-1],
    [rotation_matrix_z(T, i, n) * rotation_matrix_x_pi(T) for i = 0 : n-1])

prismatic_group(::Type{T}, n::Int) where {T} = vcat(
    [rotation_matrix_z(T, i, n) for i = 0 : n-1],
    [rotation_matrix_z(T, i, n) * rotation_matrix_x_pi(T) for i = 0 : n-1],
    [rotoreflection_matrix_z(T, i, n) for i = 0 : n-1],
    [rotation_matrix_z(T, i, n) * reflection_matrix_y(T) for i = 0 : n-1])

antiprismatic_group(::Type{T}, n::Int) where {T} = vcat(
    [rotation_matrix_z(T, i, n) for i = 0 : n-1],
    [rotoreflection_matrix_z(T, i, 2*n) for i = 1 : 2 : 2*n],
    [rotation_matrix_z(T, i, n) * reflection_matrix_y(T) for i = 0 : n-1],
    [rotoreflection_matrix_z(T, i, 2*n) * reflection_matrix_y(T)
     for i = 1 : 2 : 2*n])

######################################################## POLYHEDRAL POINT GROUPS

function chiral_tetrahedral_group(::Type{T}) where {T}
    M = SArray{Tuple{3,3},T,2,9}
    _zero = zero(T)
    _one = one(T)

    return [
        M(+_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, +_one),
        M(+_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, +_one),

        M(_zero, _zero, +_one, +_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, -_one, +_one, _zero, _zero, _zero, -_one, _zero),
        M(_zero, _zero, -_one, -_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, +_one, -_one, _zero, _zero, _zero, -_one, _zero),

        M(_zero, +_one, _zero, _zero, _zero, +_one, +_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, -_one, +_one, _zero, _zero),
        M(_zero, +_one, _zero, _zero, _zero, -_one, -_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, +_one, -_one, _zero, _zero),
    ]
end

function full_tetrahedral_group(::Type{T}) where {T}
    M = SArray{Tuple{3,3},T,2,9}
    _zero = zero(T)
    _one = one(T)

    return [
        M(+_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, +_one),
        M(+_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, +_one),

        M(+_one, _zero, _zero, _zero, _zero, +_one, _zero, +_one, _zero),
        M(+_one, _zero, _zero, _zero, _zero, -_one, _zero, -_one, _zero),
        M(-_one, _zero, _zero, _zero, _zero, +_one, _zero, -_one, _zero),
        M(-_one, _zero, _zero, _zero, _zero, -_one, _zero, +_one, _zero),

        M(_zero, _zero, +_one, +_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, -_one, +_one, _zero, _zero, _zero, -_one, _zero),
        M(_zero, _zero, -_one, -_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, +_one, -_one, _zero, _zero, _zero, -_one, _zero),

        M(_zero, +_one, _zero, +_one, _zero, _zero, _zero, _zero, +_one),
        M(_zero, -_one, _zero, +_one, _zero, _zero, _zero, _zero, -_one),
        M(_zero, -_one, _zero, -_one, _zero, _zero, _zero, _zero, +_one),
        M(_zero, +_one, _zero, -_one, _zero, _zero, _zero, _zero, -_one),

        M(_zero, +_one, _zero, _zero, _zero, +_one, +_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, -_one, +_one, _zero, _zero),
        M(_zero, +_one, _zero, _zero, _zero, -_one, -_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, +_one, -_one, _zero, _zero),

        M(_zero, _zero, +_one, _zero, +_one, _zero, +_one, _zero, _zero),
        M(_zero, _zero, -_one, _zero, -_one, _zero, +_one, _zero, _zero),
        M(_zero, _zero, +_one, _zero, -_one, _zero, -_one, _zero, _zero),
        M(_zero, _zero, -_one, _zero, +_one, _zero, -_one, _zero, _zero),
    ]
end

pyritohedral_group(::Type{T}) where {T} =
    vcat(+chiral_tetrahedral_group(T),
         -chiral_tetrahedral_group(T))

function chiral_octahedral_group(::Type{T}) where {T}
    M = SArray{Tuple{3,3},T,2,9}
    _zero = zero(T)
    _one = one(T)

    return [
        M(+_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, +_one),
        M(+_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, +_one),

        M(_zero, _zero, +_one, +_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, -_one, +_one, _zero, _zero, _zero, -_one, _zero),
        M(_zero, _zero, -_one, -_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, +_one, -_one, _zero, _zero, _zero, -_one, _zero),

        M(_zero, +_one, _zero, _zero, _zero, +_one, +_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, -_one, +_one, _zero, _zero),
        M(_zero, +_one, _zero, _zero, _zero, -_one, -_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, +_one, -_one, _zero, _zero),

        M(+_one, _zero, _zero, _zero, _zero, +_one, _zero, -_one, _zero),
        M(+_one, _zero, _zero, _zero, _zero, -_one, _zero, +_one, _zero),
        M(-_one, _zero, _zero, _zero, _zero, +_one, _zero, +_one, _zero),
        M(-_one, _zero, _zero, _zero, _zero, -_one, _zero, -_one, _zero),

        M(_zero, -_one, _zero, +_one, _zero, _zero, _zero, _zero, +_one),
        M(_zero, +_one, _zero, +_one, _zero, _zero, _zero, _zero, -_one),
        M(_zero, +_one, _zero, -_one, _zero, _zero, _zero, _zero, +_one),
        M(_zero, -_one, _zero, -_one, _zero, _zero, _zero, _zero, -_one),

        M(_zero, _zero, +_one, _zero, -_one, _zero, +_one, _zero, _zero),
        M(_zero, _zero, -_one, _zero, +_one, _zero, +_one, _zero, _zero),
        M(_zero, _zero, +_one, _zero, +_one, _zero, -_one, _zero, _zero),
        M(_zero, _zero, -_one, _zero, -_one, _zero, -_one, _zero, _zero),
    ]
end

full_octahedral_group(::Type{T}) where {T} =
    vcat(+chiral_octahedral_group(T),
         -chiral_octahedral_group(T))

function chiral_icosahedral_group(::Type{T}) where {T}
    M = SArray{Tuple{3,3},T,2,9}
    _zero = zero(T)
    _one = one(T)
    two = _one + _one
    four = two + two
    five = four + _one
    half = inv(two)
    quarter = inv(four)
    hphi = quarter * (sqrt(five) + _one)
    hpsi = quarter * (sqrt(five) - _one)

    return [
        # Identity matrix
        M(+_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, +_one),

        # 180-degree rotations about coordinate axes
        M(+_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, +_one, _zero, _zero, _zero, -_one),
        M(-_one, _zero, _zero, _zero, -_one, _zero, _zero, _zero, +_one),

        # 120-degree rotations about [±1, ±1, ±1]
        M(_zero, _zero, +_one, +_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, -_one, +_one, _zero, _zero, _zero, -_one, _zero),
        M(_zero, _zero, -_one, -_one, _zero, _zero, _zero, +_one, _zero),
        M(_zero, _zero, +_one, -_one, _zero, _zero, _zero, -_one, _zero),
        M(_zero, +_one, _zero, _zero, _zero, +_one, +_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, -_one, +_one, _zero, _zero),
        M(_zero, +_one, _zero, _zero, _zero, -_one, -_one, _zero, _zero),
        M(_zero, -_one, _zero, _zero, _zero, +_one, -_one, _zero, _zero),

        M(+half, +hphi, +hpsi, +hphi, -hpsi, -half, -hpsi, +half, -hphi),
        M(+half, +hphi, +hpsi, -hphi, +hpsi, +half, +hpsi, -half, +hphi),
        M(+half, +hphi, -hpsi, +hphi, -hpsi, +half, +hpsi, -half, -hphi),
        M(+half, +hphi, -hpsi, -hphi, +hpsi, -half, -hpsi, +half, +hphi),
        M(+half, -hphi, +hpsi, +hphi, +hpsi, -half, +hpsi, +half, +hphi),
        M(+half, -hphi, +hpsi, -hphi, -hpsi, +half, -hpsi, -half, -hphi),
        M(+half, -hphi, -hpsi, +hphi, +hpsi, +half, -hpsi, -half, +hphi),
        M(+half, -hphi, -hpsi, -hphi, -hpsi, -half, +hpsi, +half, -hphi),

        M(-half, +hphi, +hpsi, +hphi, +hpsi, +half, +hpsi, +half, -hphi),
        M(-half, +hphi, +hpsi, -hphi, -hpsi, -half, -hpsi, -half, +hphi),
        M(-half, +hphi, -hpsi, +hphi, +hpsi, -half, -hpsi, -half, -hphi),
        M(-half, +hphi, -hpsi, -hphi, -hpsi, +half, +hpsi, +half, +hphi),
        M(-half, -hphi, +hpsi, +hphi, -hpsi, +half, -hpsi, +half, +hphi),
        M(-half, -hphi, +hpsi, -hphi, +hpsi, -half, +hpsi, -half, -hphi),
        M(-half, -hphi, -hpsi, +hphi, -hpsi, -half, +hpsi, -half, +hphi),
        M(-half, -hphi, -hpsi, -hphi, +hpsi, +half, -hpsi, +half, -hphi),

        M(+hphi, +hpsi, +half, +hpsi, +half, -hphi, -half, +hphi, +hpsi),
        M(+hphi, +hpsi, +half, -hpsi, -half, +hphi, +half, -hphi, -hpsi),
        M(+hphi, +hpsi, -half, +hpsi, +half, +hphi, +half, -hphi, +hpsi),
        M(+hphi, +hpsi, -half, -hpsi, -half, -hphi, -half, +hphi, -hpsi),
        M(+hphi, -hpsi, +half, +hpsi, -half, -hphi, +half, +hphi, -hpsi),
        M(+hphi, -hpsi, +half, -hpsi, +half, +hphi, -half, -hphi, +hpsi),
        M(+hphi, -hpsi, -half, +hpsi, -half, +hphi, -half, -hphi, -hpsi),
        M(+hphi, -hpsi, -half, -hpsi, +half, -hphi, +half, +hphi, +hpsi),

        M(-hphi, +hpsi, +half, +hpsi, -half, +hphi, +half, +hphi, +hpsi),
        M(-hphi, +hpsi, +half, -hpsi, +half, -hphi, -half, -hphi, -hpsi),
        M(-hphi, +hpsi, -half, +hpsi, -half, -hphi, -half, -hphi, +hpsi),
        M(-hphi, +hpsi, -half, -hpsi, +half, +hphi, +half, +hphi, -hpsi),
        M(-hphi, -hpsi, +half, +hpsi, +half, +hphi, -half, +hphi, -hpsi),
        M(-hphi, -hpsi, +half, -hpsi, -half, -hphi, +half, -hphi, +hpsi),
        M(-hphi, -hpsi, -half, +hpsi, +half, -hphi, +half, -hphi, -hpsi),
        M(-hphi, -hpsi, -half, -hpsi, -half, +hphi, -half, +hphi, +hpsi),

        M(+hpsi, +half, +hphi, +half, -hphi, +hpsi, +hphi, +hpsi, -half),
        M(+hpsi, +half, +hphi, -half, +hphi, -hpsi, -hphi, -hpsi, +half),
        M(+hpsi, +half, -hphi, +half, -hphi, -hpsi, -hphi, -hpsi, -half),
        M(+hpsi, +half, -hphi, -half, +hphi, +hpsi, +hphi, +hpsi, +half),
        M(+hpsi, -half, +hphi, +half, +hphi, +hpsi, -hphi, +hpsi, +half),
        M(+hpsi, -half, +hphi, -half, -hphi, -hpsi, +hphi, -hpsi, -half),
        M(+hpsi, -half, -hphi, +half, +hphi, -hpsi, +hphi, -hpsi, +half),
        M(+hpsi, -half, -hphi, -half, -hphi, +hpsi, -hphi, +hpsi, -half),

        M(-hpsi, +half, +hphi, +half, +hphi, -hpsi, -hphi, +hpsi, -half),
        M(-hpsi, +half, +hphi, -half, -hphi, +hpsi, +hphi, -hpsi, +half),
        M(-hpsi, +half, -hphi, +half, +hphi, +hpsi, +hphi, -hpsi, -half),
        M(-hpsi, +half, -hphi, -half, -hphi, -hpsi, -hphi, +hpsi, +half),
        M(-hpsi, -half, +hphi, +half, -hphi, -hpsi, +hphi, +hpsi, +half),
        M(-hpsi, -half, +hphi, -half, +hphi, +hpsi, -hphi, -hpsi, -half),
        M(-hpsi, -half, -hphi, +half, -hphi, +hpsi, -hphi, -hpsi, +half),
        M(-hpsi, -half, -hphi, -half, +hphi, -hpsi, +hphi, +hpsi, -half),
    ]
end

full_icosahedral_group(::Type{T}) where {T} =
    vcat(+chiral_icosahedral_group(T),
         -chiral_icosahedral_group(T))

################################################### DEGENERATE POLYHEDRAL ORBITS

function tetrahedron_vertices(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _one = one(T)
    irt3 = inv(sqrt(_one + _one + _one))
    return [
        V(+irt3, +irt3, +irt3),
        V(+irt3, -irt3, -irt3),
        V(-irt3, +irt3, -irt3),
        V(-irt3, -irt3, +irt3),
    ]
end

function tetrahedron_edge_centers(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _one = one(T)
    _zero = zero(T)
    return [
        V(+_one, _zero, _zero),
        V(-_one, _zero, _zero),
        V(_zero, +_one, _zero),
        V(_zero, -_one, _zero),
        V(_zero, _zero, +_one),
        V(_zero, _zero, -_one),
    ]
end

function tetrahedron_face_centers(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _one = one(T)
    irt3 = inv(sqrt(_one + _one + _one))
    return [
        V(-irt3, -irt3, -irt3),
        V(-irt3, +irt3, +irt3),
        V(+irt3, -irt3, +irt3),
        V(+irt3, +irt3, -irt3),
    ]
end

function tetrahedron_rotoinversion_centers(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _zero = zero(T)
    _one = one(T)
    irt2 = inv(sqrt(_one + _one))
    return [
        V(_zero, +irt2, +irt2),
        V(_zero, +irt2, -irt2),
        V(_zero, -irt2, +irt2),
        V(_zero, -irt2, -irt2),
        V(+irt2, _zero, +irt2),
        V(-irt2, _zero, +irt2),
        V(+irt2, _zero, -irt2),
        V(-irt2, _zero, -irt2),
        V(+irt2, +irt2, _zero),
        V(+irt2, -irt2, _zero),
        V(-irt2, +irt2, _zero),
        V(-irt2, -irt2, _zero),
    ]
end

function pyritohedron_vertices(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _one = one(T)
    irt3 = inv(sqrt(_one + _one + _one))
    return [
        V(+irt3, +irt3, +irt3),
        V(+irt3, +irt3, -irt3),
        V(+irt3, -irt3, +irt3),
        V(+irt3, -irt3, -irt3),
        V(-irt3, +irt3, +irt3),
        V(-irt3, +irt3, -irt3),
        V(-irt3, -irt3, +irt3),
        V(-irt3, -irt3, -irt3),
    ]
end

octahedron_vertices(::Type{T}) where {T} =
    tetrahedron_edge_centers(T)

octahedron_edge_centers(::Type{T}) where {T} =
    tetrahedron_rotoinversion_centers(T)

octahedron_face_centers(::Type{T}) where {T} =
    pyritohedron_vertices(T)

function icosahedron_vertices(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _zero = zero(T)
    _one = one(T)
    two = _one + _one
    five = two + two + _one
    half = inv(two)
    sqrt_five = sqrt(five)
    phi = half * (sqrt_five + _one)
    psi = half * (sqrt_five - _one)
    qdrt_five = sqrt(sqrt_five)
    a = inv(qdrt_five * sqrt(phi))
    b = inv(qdrt_five * sqrt(psi))
    return [
        V(_zero, +b, +a),
        V(_zero, -b, +a),
        V(_zero, +b, -a),
        V(_zero, -b, -a),
        V(+a, _zero, +b),
        V(+a, _zero, -b),
        V(-a, _zero, +b),
        V(-a, _zero, -b),
        V(+b, +a, _zero),
        V(-b, +a, _zero),
        V(+b, -a, _zero),
        V(-b, -a, _zero),
    ]
end

function icosahedron_edge_centers(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _zero = zero(T)
    _one = one(T)
    two = _one + _one
    four = two + two
    five = four + _one
    half = inv(two)
    quarter = inv(four)
    hphi = quarter * (sqrt(five) + _one)
    hpsi = quarter * (sqrt(five) - _one)
    return [
        V(+_one, _zero, _zero),
        V(-_one, _zero, _zero),
        V(_zero, +_one, _zero),
        V(_zero, -_one, _zero),
        V(_zero, _zero, +_one),
        V(_zero, _zero, -_one),
        V(+half, +hphi, +hpsi),
        V(+half, +hphi, -hpsi),
        V(+half, -hphi, +hpsi),
        V(+half, -hphi, -hpsi),
        V(-half, +hphi, +hpsi),
        V(-half, +hphi, -hpsi),
        V(-half, -hphi, +hpsi),
        V(-half, -hphi, -hpsi),
        V(+hpsi, +half, +hphi),
        V(-hpsi, +half, +hphi),
        V(+hpsi, +half, -hphi),
        V(-hpsi, +half, -hphi),
        V(+hpsi, -half, +hphi),
        V(-hpsi, -half, +hphi),
        V(+hpsi, -half, -hphi),
        V(-hpsi, -half, -hphi),
        V(+hphi, +hpsi, +half),
        V(+hphi, -hpsi, +half),
        V(-hphi, +hpsi, +half),
        V(-hphi, -hpsi, +half),
        V(+hphi, +hpsi, -half),
        V(+hphi, -hpsi, -half),
        V(-hphi, +hpsi, -half),
        V(-hphi, -hpsi, -half),
    ]
end

function icosahedron_face_centers(::Type{T}) where {T}
    V = SArray{Tuple{3},T,1,3}
    _zero = zero(T)
    _one = one(T)
    three = _one + _one + _one
    five = three + _one + _one
    six = three + three
    irt3 = inv(sqrt(three))
    a = sqrt((three + sqrt(five)) / six)
    b = sqrt((three - sqrt(five)) / six)
    return [
        V(_zero, +b, +a),
        V(_zero, -b, +a),
        V(_zero, +b, -a),
        V(_zero, -b, -a),
        V(+a, _zero, +b),
        V(+a, _zero, -b),
        V(-a, _zero, +b),
        V(-a, _zero, -b),
        V(+b, +a, _zero),
        V(-b, +a, _zero),
        V(+b, -a, _zero),
        V(-b, -a, _zero),
        V(+irt3, +irt3, +irt3),
        V(+irt3, +irt3, -irt3),
        V(+irt3, -irt3, +irt3),
        V(+irt3, -irt3, -irt3),
        V(-irt3, +irt3, +irt3),
        V(-irt3, +irt3, -irt3),
        V(-irt3, -irt3, +irt3),
        V(-irt3, -irt3, -irt3),
    ]
end

const CHIRAL_POLYHEDRAL_POINT_GROUPS = [
    (chiral_tetrahedral_group, [tetrahedron_vertices,
        tetrahedron_edge_centers, tetrahedron_face_centers]),
    (chiral_octahedral_group, [octahedron_vertices,
        octahedron_edge_centers, octahedron_face_centers]),
    (chiral_icosahedral_group, [icosahedron_vertices,
        icosahedron_edge_centers, icosahedron_face_centers]),
]

const POLYHEDRAL_POINT_GROUPS = [
    (chiral_tetrahedral_group, [tetrahedron_vertices,
        tetrahedron_edge_centers, tetrahedron_face_centers]),
    (full_tetrahedral_group, [
        tetrahedron_vertices, tetrahedron_edge_centers,
        tetrahedron_face_centers, tetrahedron_rotoinversion_centers]),
    (pyritohedral_group, [
        pyritohedron_vertices, tetrahedron_edge_centers]),
    (chiral_octahedral_group, [octahedron_vertices,
        octahedron_edge_centers, octahedron_face_centers]),
    (full_octahedral_group, [octahedron_vertices,
        octahedron_edge_centers, octahedron_face_centers]),
    (chiral_icosahedral_group, [icosahedron_vertices,
        icosahedron_edge_centers, icosahedron_face_centers]),
    (full_icosahedral_group, [icosahedron_vertices,
        icosahedron_edge_centers, icosahedron_face_centers]),
]

####################################################### ABSTRACT GROUP STRUCTURE

function inf_norm(xs::AbstractArray{T}) where {T}
    result = zero(T)
    @simd ivdep for x in xs
        result = max(result, abs(x))
    end
    return result
end

function multiplication_table(
        group::Vector{SArray{Tuple{N,N},T,2,M}}) where {N,T,M}
    n = length(group)
    indices, dists = nn(
        KDTree(SVector{M,T}.(group)),
        [SVector{M,T}(g * h) for h in group for g in group])
    return (reshape(indices, n, n), inf_norm(dists))
end

function count_central_elements(mul_table::Matrix{Int})
    m, n = size(mul_table)
    @assert m == n
    return count(all(
        @inbounds mul_table[i,j] == mul_table[j,i]
        for j = 1 : n) for i = 1 : n)
end

################################################################ ORBIT STRUCTURE

function rotation_axis(mat::SArray{Tuple{3,3},T,2,9}, epsilon) where {T}
    @assert inf_norm(mat' * mat - one(mat)) <= epsilon
    u, s, v = @suppress svd(mat - sign(det(mat)) * one(mat))
    @assert s[1] > epsilon
    @assert s[2] > epsilon
    @assert zero(T) <= s[3] <= epsilon
    return v[:,3]
end

function connected_components(adjacency_lists::Dict{V,Vector{V}}) where {V}
    visited = Dict{V,Bool}()
    for (v, l) in adjacency_lists
        visited[v] = false
    end
    components = Vector{V}[]
    for (v, l) in adjacency_lists
        if !visited[v]
            visited[v] = true
            current_component = [v]
            to_visit = copy(l)
            while !isempty(to_visit)
                w = pop!(to_visit)
                if !visited[w]
                    visited[w] = true
                    push!(current_component, w)
                    append!(to_visit, adjacency_lists[w])
                end
            end
            push!(components, current_component)
        end
    end
    @assert allunique(vcat(components...))
    return components
end

function degenerate_orbits(group::Vector{SArray{Tuple{3,3},T,2,9}},
                           epsilon) where {T}
    Point = SArray{Tuple{3},T,1,3}
    points = Vector{Point}()
    for (i, g) in enumerate(group)
        try
            axis = rotation_axis(g, epsilon)
            push!(points, +axis)
            push!(points, -axis)
        catch e
            if !(e isa AssertionError)
                rethrow(e)
            end
        end
    end
    clusters = Vector{Tuple{Point,Vector{Point}}}()
    for point in points
        found = false
        for (center, cluster) in clusters
            if norm(point - center) <= epsilon
                found = true
                push!(cluster, point)
                break
            end
        end
        if !found
            push!(clusters, (point, [point]))
        end
    end
    @assert all(
        norm(p - q) <= epsilon
        for (_, cluster) in clusters
        for p in cluster for q in cluster)
    n = length(clusters)
    @assert all(
        @inbounds norm(clusters[i][1] - clusters[j][1]) > epsilon
        for i = 1 : n-1 for j = i+1 : n)
    adjacency_lists = Dict(i => Int[] for i = 1 : n)
    for (i, (center, _)) in enumerate(clusters)
        for g in group
            p = g * center
            dist, j = minimum(
                (norm(p - q), k)
                for (k, (q, _)) in enumerate(clusters))
            @assert dist <= epsilon
            if i != j
                @inbounds push!(adjacency_lists[i], j)
            end
        end
    end
    return [[@inbounds clusters[i][1] for i in comp]
            for comp in connected_components(adjacency_lists)]
end

############################################################# FINDING ISOMETRIES

function labeled_distances(points::Vector{SVector{N,T}}) where {T,N}
    return [(sqrt(sum((points[i] - points[j]).^2)), i, j)
            for j = 2 : length(points)
            for i = 1 : j-1]
end

function bucket_by_first(items::Vector{T}, epsilon) where {T}
    result = Vector{T}[]
    if length(items) == 0
        return result
    end
    push!(result, [items[1]])
    for i = 2 : length(items)
        if abs(items[i][1] - result[end][end][1]) <= epsilon
            push!(result[end], items[i])
        else
            push!(result, [items[i]])
        end
    end
    for bucket in result
        @assert abs(bucket[end][1] - bucket[1][1]) <=epsilon
    end
    return result
end

middle(x::AbstractVector) = x[(length(x) + 1) >> 1]

function candidate_isometries(
        a_points::Vector{SVector{3,T}},
        b_points::Vector{SVector{3,T}}, epsilon) where {T}
    _one = one(T)
    two = _one + _one
    four = two + two
    one_fourth = inv(four)
    seven_fourths = (four + two + _one) * one_fourth
    result = SArray{Tuple{3,3},T,2,9}[]
    a_buckets = bucket_by_first(sort!(labeled_distances(a_points)), epsilon)
    b_buckets = bucket_by_first(sort!(labeled_distances(b_points)), epsilon)
    if length.(a_buckets) != length.(b_buckets)
        return result
    end
    a_midpoints = [middle(bucket)[1] for bucket in a_buckets]
    b_midpoints = [middle(bucket)[1] for bucket in b_buckets]
    for (a, b) in zip(a_midpoints, b_midpoints)
        if !(abs(a - b) <= epsilon)
            return result
        end
    end
    selected_bucket = minimum(
        (length(bucket), abs(one(T) - middle(bucket)[1]), i)
        for (i, bucket) in enumerate(a_buckets)
        if one_fourth <= middle(bucket)[1] <= seven_fourths)[3]
    _, i, j = middle(a_buckets[selected_bucket])
    u1, v1 = a_points[i], a_points[j]
    pos_mat = inv(hcat(u1, v1, cross(u1, v1)))
    neg_mat = inv(hcat(u1, v1, cross(v1, u1)))
    for (_, k, l) in b_buckets[selected_bucket]
        u2, v2 = b_points[k], b_points[l]
        fwd_mat = hcat(u2, v2, cross(u2, v2))
        rev_mat = hcat(v2, u2, cross(v2, u2))
        push!(result, fwd_mat * pos_mat)
        push!(result, fwd_mat * neg_mat)
        push!(result, rev_mat * pos_mat)
        push!(result, rev_mat * neg_mat)
    end
    for mat in result
        @assert inf_norm(mat' * mat - one(mat)) <= epsilon
    end
    return result
end

function isometric(
        a_points::Vector{SVector{3,T}},
        b_points::Vector{SVector{3,T}}, epsilon) where {T}
    b_tree = KDTree(b_points)
    for mat in candidate_isometries(a_points, b_points, epsilon)
        indices, distances = nn(b_tree, [mat * p for p in a_points])
        if allunique(indices) && (inf_norm(distances) <= epsilon)
            return true
        end
    end
    return false
end

function isometries(
        a_points::Vector{SVector{3,T}},
        b_points::Vector{SVector{3,T}}, epsilon) where {T}
    result = SArray{Tuple{3,3},T,2,9}[]
    b_tree = KDTree(b_points)
    for mat in candidate_isometries(a_points, b_points, epsilon)
        indices, distances = nn(b_tree, [mat * p for p in a_points])
        if allunique(indices) && (inf_norm(distances) <= epsilon)
            push!(result, mat)
        end
    end
    return result
end

isometries(points::Vector{SVector{3,T}}, epsilon) where {T} =
    isometries(points, points, epsilon)

##################################################### POINT GROUP IDENTIFICATION

@enum MatrixType begin
    IdentityMatrix
    PureInversion
    PureRotation # 180-degree rotation
    PureReflection
    GeneralRotation
    GeneralRotoinversion
end

function identify_point_group(group::Vector{SArray{Tuple{3,3},T,2,9}},
                              epsilon) where {T}
    n = length(group)
    mul_table, dist = multiplication_table(group)
    @assert dist <= epsilon
    @assert (n, n) == size(mul_table)

    identity_matrix = one(SArray{Tuple{3,3},T,2,9})
    if n == 1
        return ("C_1", identity_matrix)
    end

    id_dist, id = minimum(
        (inf_norm(mat - identity_matrix), i)
        for (i, mat) in enumerate(group))
    @assert id_dist <= epsilon
    @assert issorted(view(mul_table, id, :))
    @assert issorted(view(mul_table, :, id))

    is_abelian = (mul_table == mul_table')

    order_table = zeros(Int, n)
    for g = 1 : n
        acc = g
        ord = 1
        while acc != id
            acc = mul_table[acc,g]
            ord += 1
        end
        order_table[g] = ord
    end
    @assert minimum(order_table) == 1
    max_order = maximum(order_table)
    is_cyclic = (max_order == n)

    _zero = zero(T)
    _one = one(T)
    two = _one + _one
    matrix_types = Vector{MatrixType}(undef, n)
    axes = zeros(SVector{3,T}, n)

    for (i, mat) in enumerate(group)
        u, s, v = @suppress svd(mat - one(mat); full=true)
        if inf_norm(s - [_zero, _zero, _zero]) <= epsilon
            @assert id == i
            @assert !signbit(det(mat))
            matrix_types[i] = IdentityMatrix
        elseif inf_norm(s - [two, two, two]) <= epsilon
            @assert signbit(det(mat))
            matrix_types[i] = PureInversion
        elseif inf_norm(s - [two, two, _zero]) <= epsilon
            @assert !signbit(det(mat))
            matrix_types[i] = PureRotation
            axes[i] = v[:,3]
        elseif inf_norm(s - [two, _zero, _zero]) <= epsilon
            @assert signbit(det(mat))
            matrix_types[i] = PureReflection
            axes[i] = v[:,1]
        elseif abs(s[1] - s[2]) <= epsilon
            @assert !signbit(det(mat))
            @assert s[3] <= epsilon
            matrix_types[i] = GeneralRotation
            axes[i] = v[:,3]
        else
            @assert signbit(det(mat))
            @assert abs(s[1] - two) <= epsilon
            @assert abs(s[2] - s[3]) <= epsilon
            matrix_types[i] = GeneralRotoinversion
            axes[i] = v[:,1]
        end
    end

    has_inversion = (PureInversion in matrix_types)
    has_pure_reflection = (PureReflection in matrix_types)
    is_chiral = (!has_inversion && !has_pure_reflection
                                && !(GeneralRotoinversion in matrix_types))

    principal_axes = [axes[i]
        for i = 1 : n
        if order_table[i] == max_order && (
            matrix_types[i] == PureRotation ||
            matrix_types[i] == GeneralRotation ||
            matrix_types[i] == GeneralRotoinversion)]
    if isempty(principal_axes)
        principal_axis = zero(SVector{3,T})
        axis_rotation_matrix = identity_matrix
    else
        principal_axis = first(principal_axes)
        axis_rotation_matrix = rotation_matrix_to_z(principal_axis)
    end

    if n == 2
        if is_chiral
            return ("C_2", axis_rotation_matrix)
        elseif PureInversion in matrix_types
            return ("C_i", identity_matrix)
        else
            return ("C_s", rotation_matrix_to_z(axes[3 - id]))
        end
    end

    orthogonal_rotation_indices = [i
        for i = 1 : n
        if (matrix_types[i] == PureRotation &&
            abs(principal_axis' * axes[i]) <= epsilon)]
    if !isempty(orthogonal_rotation_indices)
        orthogonal_axis =
            axis_rotation_matrix * axes[first(orthogonal_rotation_indices)]
        @assert abs(orthogonal_axis[3]) <= epsilon
        axis_rotation_matrix =
            rotation_matrix_to_x(orthogonal_axis) * axis_rotation_matrix
    end

    orthogonal_reflection_indices = [i
        for i = 1 : n
        if (matrix_types[i] == PureReflection &&
            abs(principal_axis' * axes[i]) <= epsilon)]
    if !isempty(orthogonal_reflection_indices)
        secondary_axis =
            axis_rotation_matrix * axes[first(orthogonal_reflection_indices)]
        @assert abs(secondary_axis[3]) <= epsilon
        secondary_rotation_matrix =
            rotation_matrix_to_y(secondary_axis) * axis_rotation_matrix
    end

    if n == 4
        if is_cyclic
            return (is_chiral ? "C_4" : "S_4", axis_rotation_matrix)
        elseif is_chiral
            return ("D_2", axis_rotation_matrix)
        elseif has_inversion
            return ("C_2h", axis_rotation_matrix)
        else
            return ("C_2v", secondary_rotation_matrix)
        end
    end

    has_max_order_negative_element = any(
        order_table[i] == max_order && (
            matrix_types[i] == PureInversion ||
            matrix_types[i] == PureReflection ||
            matrix_types[i] == GeneralRotoinversion)
        for i = 1 : n)

    if n == 8
        if is_cyclic
            return (is_chiral ? "C_8" : "S_8", axis_rotation_matrix)
        elseif is_chiral
            return ("D_4", axis_rotation_matrix)
        elseif is_abelian
            if max_order == 2
                return ("D_2h", axis_rotation_matrix)
            else
                @assert max_order == 4
                return ("C_4h", axis_rotation_matrix)
            end
        elseif has_max_order_negative_element
            return ("D_2d", secondary_rotation_matrix)
        else
            return ("C_4v", secondary_rotation_matrix)
        end
    end

    if is_cyclic
        @assert is_abelian
        if is_chiral
            return ("C_$n", axis_rotation_matrix)
        elseif PureReflection in matrix_types
            @assert iseven(n)
            @assert isodd(n >> 1)
            return ("C_$(n >> 1)h", axis_rotation_matrix)
        else
            @assert iseven(n)
            return ("S_$n", axis_rotation_matrix)
        end
    end

    if is_abelian
        @assert iseven(n)
        @assert iseven(n >> 1)
        return ("C_$(n >> 1)h", axis_rotation_matrix)
    end

    is_axial = (max_order >= 3) && all(
        abs(sqrt(abs(axes[i]' * axes[j])) - _one) <= epsilon
        for i = 1 : n, j = 1 : n
        if (order_table[i] == max_order &&
            order_table[j] == max_order))

    if is_axial
        @assert iseven(n)
        if is_chiral
            return ("D_$(n >> 1)", axis_rotation_matrix)
        elseif !has_max_order_negative_element
            refl_index = first(i for i = 1 : n
                               if matrix_types[i] == PureReflection)
            refl_axis = axis_rotation_matrix * axes[refl_index]
            return ("C_$(n >> 1)v",
                    rotation_matrix_to_y(refl_axis) * axis_rotation_matrix)
        else
            @assert iseven(n >> 1)
            if xor(iseven(n >> 2), PureInversion in matrix_types)
                return ("D_$(n >> 2)d", secondary_rotation_matrix)
            else
                return ("D_$(n >> 2)h", axis_rotation_matrix)
            end
        end
    else
        if n == 12
            return ("T", axis_rotation_matrix)
        elseif n == 48
            return ("O_h", axis_rotation_matrix)
        elseif n == 60
            return ("I", axis_rotation_matrix)
        elseif n == 120
            return ("I_h", axis_rotation_matrix)
        else
            @assert n == 24
            if is_chiral
                return ("O", axis_rotation_matrix)
            else
                return (PureInversion in matrix_types ? "T_h" : "T_d",
                        axis_rotation_matrix)
            end
        end
    end
end

###################################################### TEST POINT CONFIGURATIONS

function equatorial_points(height::T, offset::T, n::Int) where {T}
    _one = one(T)
    two = _one + _one
    a = sqrt(_one + height * height)
    return [SVector{3,T}(
            cospi(two * (T(i) + offset) / n),
            sinpi(two * (T(i) + offset) / n),
            height) / a
        for i = 0 : n-1]
end

function alternating_equatorial_points(height::T, offset::T, n::Int) where {T}
    @assert iseven(n)
    _one = one(T)
    two = _one + _one
    a = sqrt(_one + height * height)
    return [SVector{3,T}(
            cospi(two * (T(i) + offset) / n),
            sinpi(two * (T(i) + offset) / n),
            iseven(i) ? height : -height) / a
        for i = 0 : n-1]
end

cyclic_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(-inv(T(5)), inv(T(20)), n),
    equatorial_points(inv(T(10)), T(3)/T(20), n))

rotoreflection_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(zero(T), inv(T(5)), n),
    alternating_equatorial_points(inv(T(10)), inv(T(5)), n))

pinwheel_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(+inv(T(5)), inv(T(10)), n),
    equatorial_points(-inv(T(5)), inv(T(10)), n))

pyramid_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(inv(T(5)), zero(T), n),
    equatorial_points(zero(T), +inv(T(10)), n),
    equatorial_points(zero(T), -inv(T(10)), n))

dihedral_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(zero(T), +inv(T(10)), n),
    equatorial_points(zero(T), -inv(T(10)), n),
    equatorial_points(+inv(T(5)), +inv(T(10)), n),
    equatorial_points(-inv(T(5)), -inv(T(10)), n))

prismatic_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(zero(T), +inv(T(10)), n),
    equatorial_points(zero(T), -inv(T(10)), n),
    equatorial_points(+inv(T(5)), zero(T), n),
    equatorial_points(-inv(T(5)), zero(T), n))

antiprismatic_point_configuration(::Type{T}, n::Int) where {T} = vcat(
    equatorial_points(zero(T), zero(T), n),
    equatorial_points(zero(T), +inv(T(10)), n),
    equatorial_points(zero(T), -inv(T(10)), n),
    equatorial_points(+inv(T(5)), zero(T), n),
    equatorial_points(zero(T), inv(T(2)), n),
    equatorial_points(zero(T), T(2)/T(5), n),
    equatorial_points(zero(T), T(3)/T(5), n),
    equatorial_points(-inv(T(5)), inv(T(2)), n))

################################################ TEST POINT GROUP IDENTIFICATION

function matching_distance(a::Vector{SArray{X,T,R,N}},
                           b::Vector{SArray{X,T,R,N}}) where {X,T,R,N}
    if length(a) != length(b)
        return typemax(T)
    end
    a_points = SVector{N,T}.(a)
    b_points = SVector{N,T}.(b)
    b_tree = KDTree(b_points)
    idxs, dists = nn(b_tree, a_points)
    if !allunique(idxs)
        return typemax(T)
    else
        return maximum(dists)
    end
end

random_rotation(::Type{T}, ::Val{N}) where {T,N} =
    qr(SArray{Tuple{N,N},T,2,N*N}(randn(N, N))).Q

function test_alignment(
        points::Vector{SVector{3,T}},
        expected_group::Vector{SArray{Tuple{3,3},T,2,9}},
        epsilon) where {T}
    rotation = random_rotation(T, Val(3))
    rotated_points = [rotation * p for p in points]
    name, align = identify_point_group(
        isometries(rotated_points, epsilon),
        epsilon)
    @assert abs(det(align) - one(T)) <= epsilon
    dist = matching_distance(
        isometries([align * p for p in rotated_points], epsilon),
        expected_group)
    @assert dist <= epsilon
    return name
end

function test_axial_group_alignment(::Type{T}, n_max::Int,
                                    epsilon) where {T}
    cyclic_groups = [
        test_alignment(cyclic_point_configuration(T, n),
                       cyclic_group(T, n), epsilon)
        for n = 1 : n_max]
    @assert cyclic_groups == ["C_$n" for n = 1 : n_max]

    rotoreflection_groups = [
        test_alignment(rotoreflection_point_configuration(T, n),
                       rotoreflection_group(T, n), epsilon)
        for n = 2 : 2 : n_max]
    @assert rotoreflection_groups == [n == 2 ? "C_i" : "S_$n"
                                      for n = 2 : 2 : n_max]

    pinwheel_groups = [
        test_alignment(pinwheel_point_configuration(T, n),
                       pinwheel_group(T, n), epsilon)
        for n = 1 : n_max]
    @assert pinwheel_groups == [n == 1 ? "C_s" : "C_$(n)h"
                                for n = 1 : n_max]

    pyramid_groups = [
        test_alignment(pyramid_point_configuration(T, n),
                       pyramid_group(T, n), epsilon)
        for n = 2 : n_max]
    @assert pyramid_groups == ["C_$(n)v" for n = 2 : n_max]

    dihedral_groups = [
        test_alignment(dihedral_point_configuration(T, n),
                       dihedral_group(T, n), epsilon)
        for n = 2 : n_max]
    @assert dihedral_groups == ["D_$n" for n = 2 : n_max]

    prismatic_groups = [
        test_alignment(prismatic_point_configuration(T, n),
                       prismatic_group(T, n), epsilon)
        for n = 2 : n_max]
    @assert prismatic_groups == ["D_$(n)h" for n = 2 : n_max]

    antiprismatic_groups = [
        test_alignment(antiprismatic_point_configuration(T, n),
                       antiprismatic_group(T, n), epsilon)
        for n = 2 : n_max]
    @assert antiprismatic_groups == ["D_$(n)d" for n = 2 : n_max]
end

################################################################################

end # module PointGroups


text:
@testset "runtestTrimLefToFQuote" begin
    println("Testing runtestTrimLefToFQuote")
    include("../src/bfgaTrimLefToFQuote.jl")
    using .bfgaTrimLefToFQuote

    include("../src/Bf.jl")
    using .BfInterpreter

    @time model = test_serial()
    #@time model = test_parallel()
    prems = model.refSet[length(model.refSet)]
    println("Generation : $(model.params.currentGeneration)")
    println("& fitness : $(model.specific_fitness.fitness(prems, model.instructionsSet))")
    #bfcode :: Array{Float64,1}  = prems.dna
    println("code :  $(prems.program)") #$(join(genesToBfInstructions(bfcode), "") ) ")

    for i in 1:length(words)
        input = words[i]
        goal = _results[i]
        n = length(input)
        output1, _ = execute(prems.program, input, model.instructionsSet)
        @test length(output1) >= n
        @test collect(input1[1:n1]) == output1[1:n]
    end

    # do the test_parallel after ...
    #test_parallel()
end


text:
# load a BioNetGen defined BCR network, solve the ODEs
# plot a specific observable and compare to the BioNetGen solution

using DiffEqBase, Catalyst, Plots, OrdinaryDiffEq, Sundials, DataFrames, CSVFiles, LinearAlgebra, SparseArrays
using ReactionNetworkImporters
using TimerOutputs

# parameters
doplot = false
tf = 10000.
build_jac = false
sparse_jac = false
zeroout_jac = false

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/bcr")
fname    = joinpath(datadir, "bcr.net")
gdatfile = joinpath(datadir, "bcr.gdat")
print("getting gdat file...")
gdatdf = DataFrame(load(File(format"CSV", gdatfile), header_exists=true, spacedelim=true) )
println("done")

# we'll time the DiffEq solvers
const to = TimerOutput()
reset_timer!(to)

# BioNetGen network
@timeit to "bionetgen" prnbng = loadrxnetwork(BNGNetwork(),fname);
show(to)
rnbng = prnbng.rn
@timeit to "bODESystem" bosys = convert(ODESystem, rnbng)
show(to)
@timeit to "bODEProb" boprob = ODEProblem(bosys, Float64[], (0.,tf), Float64[])
show(to)
u = zeros(numspecies(rnbng))
p = zeros(length(parameters(rnbng)))
du = similar(u);
@timeit to "f1" boprob.f(du,u,p,0.)
@timeit to "f2" boprob.f(du,u,p,0.)
if build_jac
    J = zeros(length(u),length(u))
    #J = similar(rnbng.odefun.jac_prototype)
    @timeit to "J1" rnbng.jac(J,u,p,0.)
    @timeit to "J2" rnbng.jac(J,u,p,0.)
end
show(to)
println()

# DiffEq solver
#reset_timer!(to);
#@timeit to "BNG_CVODE_BDF-LU-1" begin bsol = solve(boprob, CVODE_BDF(), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-LU-2" begin bsol = solve(boprob, CVODE_BDF(), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-GMRES-1" begin bsol = solve(boprob, CVODE_BDF(linear_solver=:GMRES), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "BNG_CVODE_BDF-GMRES-2" begin bsol = solve(boprob, CVODE_BDF(linear_solver=:GMRES), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
# #reset_timer!(to); @timeit to "BNG_RODAS5_BDF" begin bsol2 = solve(boprob, rodas5(autodiff=false), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "KenCarp4-1" begin sol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
#@timeit to "KenCarp4-2" begin bsol = solve(boprob,KenCarp4(autodiff=false,linsolve=LinSolveFactorize(lu)), saveat=1., abstol=1e-8, reltol=1e-8); end; show(to)
@timeit to "KenCarp4-1" begin sol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)
@timeit to "KenCarp4-2" begin bsol = solve(boprob,KenCarp4(autodiff=false), abstol=1e-8, reltol=1e-8, saveat=1.); end; show(to)

# Activated Syk from DiffEq
@unpack Activated_Syk = rnbng

if doplot
    plotlyjs()
    plot(gdatdf[!,:time][2:end], gdatdf[!,:Activated_Syk][2:end], xscale=:log10, label="AsykGroup", linestyle=:dot)
    plot!(bsol.t[2:end], sol[Activated_Syk][2:end], label="Activated_Syk", xscale=:log10)
end

# test the error, note may be large in abs value though small relatively
# #norm(gdatdf[:Activated_Syk] - asynbng, Inf)
# #norm(asynbng - basyk', Inf)
norm(gdatdf[!,:Activated_Syk] - sol[Activated_Syk], Inf)

#@assert all(abs.(gdatdf[!,:Activated_Syk] - basyk') .< 1e-6 * abs.(gdatdf[:Activated_Syk]))


text:
using MLStyle

"""
    callify(mycall, ast)

Replace every `f(args...; kwargs..)` with `mycall(f, args...; kwargs...)`
"""
function callify(mycall, ast)
    leaf(x) = x
    function branch(f, head, args)
        default() = Expr(head, map(f, args)...)
        head == :call || return default()

        if first(args) == :~ && length(args) == 3
            return default()
        end

        # At this point we know it's a function call
        length(args) == 1 && return Expr(:call, mycall, first(args))

        fun = args[1]
        arg2 = args[2]

        if arg2 isa Expr && arg2.head == :parameters
            # keyword arguments (try dump(:(f(x,y;a=1, b=2))) to see this)
            return Expr(:call, mycall, arg2, fun, map(f, Base.rest(args, 3))...)
        else
            return Expr(:call, mycall, map(f, args)...)
        end

    end

    foldast(leaf, branch)(ast)
end

# struct Provenance{T,S}
#     value::T
#     sources::S
# end

# getvalue(p::Provenance) = p.value
# getvalue(x) = x

# getsources(p::Provenance) = p.sources
# getsources(x) = Set()

# function trace_provenance(f, args...; kwargs...)
#     (newargs, arg_sources) = (getvalue.(args), union(getsources.(args)...))

#     k = keys(kwargs)
#     v = values(kwargs)
#     newkwargs = NamedTuple{k}(map(getvalue, v))

#     k = keys(kwargs)
#     v = values(NamedTuple(kwargs))
#     newkwargs = NamedTuple{k}(getvalue.(v))
#     kwarg_sources = union(getsources.(args)...)

#     sources = union(arg_sources, kwarg_sources)
#     Provenance(f(newargs...; newkwargs), sources)
# end

# macro call(expr)
#     callify(expr)
# end

# julia> callify(:(f(g(x,y))))
# :(call(f, call(g, x, y)))

# julia> callify(:(f(x; a=3)))
# :(call(f, x; a = 3))

# julia> callify(:(a+b))
# :(call(+, a, b))

# julia> callify(:(call(f,3)))
# :(call(f, 3))

# f(x) = x+1

# @call f(2)

# using SymbolicUtils

# @syms x::Vector{Float64} i::Int

# @call getindex(x,i)


text:
append!(truths, [
    Axiom("Definition of And", P ∧ False ⟹ False),
    Axiom("Definition of And", False ∧ P ⟹ False),
    Axiom("Definition of And", P ∧ True ⩶ P),
    Axiom("Definition of And", True ∧ P ⩶ P),
    Axiom("And is associative", @associative(∧, P, Q, R)),
    Axiom("And is commutative", @commutative(∧, P, Q)),
    Axiom("And is idempotent", @idempotent(∧, P)),
    Axiom("De Morgan's laws", P ∧ Q ⩶ ~(~P ∨ ~Q)),
    Axiom("And is absorbing over Or", P ∧ (P ∨ Q) ⟹ Q),
    Axiom("And is distributive over Or", @distributive(∧, ∨, P, Q, R)),
    Axiom("And is factorisable over Or", @factorisable(∧, ∨, P, Q, R)),

    Axiom("Law of noncontradiction", P ∧ ~P ⩶ False),
    Axiom("Law of excluded middle", P ∨ ~P ⩶ True)

    Axiom("Definition of Or", P ∨ True ⟹ True),
    Axiom("Definition of Or", True ∨ P ⟹ True),
    Axiom("Definition of Or", P ∨ False ⩶ P),
    Axiom("Definition of Or", False ∨ P ⩶ P),
    Axiom("Or is associative", @associative(∨, P, Q, R)),
    Axiom("Or is commutative", @commutative(∨, P, Q)),
    Axiom("Or is idempotent", @idempotent(∨, P)),
    Axiom("De Morgan's laws", P ∨ Q ⩶ ~(~P ∧ ~Q)),
    Axiom("Or is absorbing over And", P ∨ (P ∧ Q) ⟹ Q),
    Axiom("Or is distributive over And", @distributive(∨, ∧, P, Q, R)),
    Axiom("Or is factorisable over And", @factorisable(∨, ∧, P, Q, R)),
    Axiom("Or is the weaker And", P ∧ Q ⟹ P ∨ Q),

    Axiom("Definition of Not", ~False ⩶ True),
    Axiom("Definition of Not", ~True ⩶ False),
    # prove these
    #Axiom("De Morgan's laws", ~(P ∨ Q), ~P ∧ ~Q),
    #Axiom("De Morgan's laws", ~(P ∨ Q), ~P ∧ ~Q),
    Axiom("Not is involute", @involute(~, P)),

    Axiom("Definition of Implies", (P ⟹ Q) ⩶ ~P ∨ Q),
    # Theorem("Contrapositive", P ⟹ Q, ~Q ⟹ ~P, Proof([])),

    Axiom("Definition of Equivalent", (P ⩶ Q) ⩶ (P ⟹ Q) ∧ (Q ⟹ P)),
])


text:
using Datetime, TimeSeries, Requests, DataFrames

module Quandl

    using Datetime, TimeSeries, Requests, DataFrames

    export quandlget,
           quandl,
           quandlsearch,
           set_auth_token

    include("api.jl")
    include("timearray.jl")
    include("utilities.jl")
    include("dataframe.jl")

end


text:
type Class_stable_ls <: abstract_ls_info
    step_size_P::Float64
    step_size_D::Float64
    num_steps::Int64
    predict_red::Float64
    frac_progress::Float64
    actual_red::Float64
    cur_merit::Float64
    do_ls::Bool

    function Class_stable_ls(iter::Class_iterate, dir::Class_point, pars::Class_parameters)
        this = new()

        this.step_size_D = 0.0
        this.step_size_P = 0.0
        this.num_steps = 0
        this.predict_red = merit_function_predicted_reduction(iter, dir, 1.0);
        this.frac_progress = NaN
        this.actual_red = NaN
        this.cur_merit = eval_merit_function(iter, pars)
        this.do_ls = this.predict_red >= 0.0

        return this
    end

    function Class_stable_ls(step_size_P::Float64,step_size_D::Float64,num_steps::Int64,frac_progress::Float64, predict_red::Float64)
        this = new(step_size_P, step_size_D, num_steps, frac_progress, predict_red)
        return this
    end
end

function accept_func_stable!(accept::abstract_ls_info, iter::Class_iterate, candidate::Class_iterate, dir::Class_point, step_size::Float64, pars::Class_parameters, timer::class_advanced_timer)
  #@show comp_merit(iter)
  if accept.predict_red < 0.0
    #@show comp_merit(candidate)
    #old_merit = eval_merit_function(iter, pars)
    old_merit = accept.cur_merit
    #new_merit = eval_merit_function(candidate, pars)
    #@show comp_merit(candidate), new_merit, old_merit
    #@show new_merit - old_merit

    accept.actual_red = eval_merit_function_difference(iter, candidate, pars)

    #@show accept.predict_red, accept.actual_red
    #predict_red = merit_function_predicted_reduction(iter, dir, step_size);
    predict_red = accept.predict_red * step_size

    if accept.actual_red > 0.0
      return :negative_progress, step_size * pars.ls.backtracking_factor
    else
      accept.frac_progress = accept.actual_red / predict_red
      if accept.frac_progress > pars.ls.predict_reduction_factor
          #@show accept.frac_progress, accept.actual_red
          return :success, step_size
      else
          return :not_enough_progress, step_size * pars.ls.backtracking_factor
      end
    end
  else
      return :not_enough_progress, step_size * pars.ls.backtracking_factor
  end
end

function accept_func!(accept::Class_stable_ls, iter::Class_iterate, candidate::Class_iterate, dir::Class_point, step_size::Float64, filter::Array{Class_filter,1}, pars::Class_parameters, timer::class_advanced_timer)
    return accept_func_stable!(accept, iter, candidate, dir, step_size, pars, timer)
end


text:
module BaseImplementation
using  ...Dispatch.Structure
import Base: showerror, ==, show, display

(met::MetaMethod)(exprs...) = met.method(exprs...)

function show(io::IO, met::MetaMethod)
  args = map(met.tree.child.children) do tree
    s = sprint(show, tree)
    s = replace(s, r"pattern`", "", 1)[1:end-1]
  end
  label = met.label == unlabeled ? "[---]" : "[$(met.label)]"

  print(io, label, "<", join(args, " "), ">")
end

function show(io::IO, table::MetaMethodTable)
  if get(io, :compact, false)
    print(io, table.name)
  else
    len = length(table.methods)
    println(io, table.name, " with $len method$(len==1 ? "" : "s"):")
    for method in table.methods
      println(io, " ", method)
    end
  end
end

showerror(io::IO, err::MetaMethodError) =
  print(io, "$(err.name) has no method that matches $(err.expr).")

end


text:
include("gherkin/feature_test.jl")
include("gherkin/scenario_test.jl")
include("gherkin/scenario_outline_test.jl")
include("gherkin/tag_test.jl")

include("executor_test.jl")
include("step_def_test.jl")
include("result_accumulator_test.jl")
include("asserts_test.jl")
include("outlines_test.jl")


text:
"""
    MOSUM(x, G, var_est_method, alpha, criterion, eta, epsilon)

Runs the MOSUM procedure for the univariate data x with bandwidth G, and returns the number and position of found changepoints.

See also: [`@BS`](@ref), [`@WBS`](@ref)

# Returns
* `out::Dict{String,Union{Bool, Float64, Int64, Array}}`:
  * `out["number"]` : Integer number of detected changes
  * `out["Reject"]` : Boolean hypothesis test outcome
  * `out["threshold"]` : Float of testing threshold
  * `out["detector"]` : vector of mosum detector values
  * `out["var.estimation"]` : vector of estimated variance values

# Example
```julia-repl
# Sample Normal time series with changing mean
n = 1000
λ = 300
μ, σ = Normal(0.0, 10.0), 1.0
x, cps = @changepoint_sampler n λ Normal(μ, σ)
G = 150
# Run MOSUM procedure
MOSUM_out = MOSUM(x, G)
# Plot MOSUM detector
using Gadfly
Gadfly.plot(MOSUM_out)
```

# References
Eichinger, Birte, and Claudia Kirch. "A MOSUM procedure for the estimation of multiple random change points." Bernoulli 24.1 (2018): 526-564.
Meier, Alexander, Claudia Kirch, and Haeran Cho. "mosum: A package for moving sums in change point analysis." (2018).
"""
function MOSUM( x::Array{Float64} , G::Int64; var_est_method::String = "mosum",
     alpha::Float64 = 0.1, criterion::String = "eta", eta::Float64 = 0.4, epsilon::Float64 = 0.2)

    CP = Array{Int64}(undef,0)
    q = 0
    Reject = false
    n = length(x)

    csum = cumsum(x)
    sqsum = cumsum(x.^2)
    variance = zeros(Float64, n)
    lvar = zeros(Float64, n)
    rvar = zeros(Float64, n)
    lsqsums = zeros(Float64, n)
    rsqsums =  zeros(Float64, n)
    detector = zeros(Float64, n)
    lsums = zeros(Float64, n)
    rsums = zeros(Float64, n)

    #populate
    for t in 1:(n-G)
        rsums[t] = csum[t+G] - csum[t]
        lsums[t+G] = rsums[t]

        rsqsums[t]  = sqsum[t+G] - sqsum[t]
        lsqsums[t+G] = rsqsums[t]

        rlmean = rsums[t]/G
        rvar[t] = rsqsums[t]/G -  (1/G - 2)*(rlmean^2)
        lvar[t+G] = rvar[t]
    end

    # normalise
    for t in (G+1):(n-G)
        if var_est_method == "mosum"
            variance[t] = (rvar[t] + lvar[t])/ 2
        end
        if var_est_method == "mosum.min"
            variance[t] = min( rvar[t], lvar[t])
        end
        detector[t] = abs(rsums[t] - lsums[t])/ sqrt(2*G*variance[t])
    end

    #threshold
    a = sqrt(2*n/G); b = 2log(n/G) + log(log(n/G))/2 + log(3/2) - log(pi)/2
    c = - log(log( (1-alpha)^(-.5) ))
    D = (b+c)/a

    mosum_stat, k = findmax(detector)

    if mosum_stat > D
        got_cps = get_cps_mosum(detector, D, G, criterion, eta, epsilon)
        CP = got_cps[1]
        q = got_cps[2]
        Reject = got_cps[3]
    end

#    cost = 0.0
#    CP = sort(CP)
#    for j in 1:(length(CP)-1)
#        cost = cost + segment_cost(CP[j]+1,CP[j+1]) + pen
#    end
#    cost = cost + segment_cost(CP[end]+1,n)
    # organise output into a dictionary
    out = Dict{String, Union{Array, Bool, Int64, Float64} }()
    out["number"] = q
    out["detector"] = detector
    out["Reject"] = Reject
    out["changepoints"] = CP
    out["threshold"] = D
    out["var.estimation"] = variance
    return out

end

function get_cps_mosum(stat::Array{Float64} , D_n::Float64, G::Int64,  criterion::String = ["epsilon","eta"][1],
     eta::Float64 = 0.4, epsilon::Float64 = 0.2 )
  cps = Array{Int64}(undef,0) ##assign empty cps
  Reject = false
  if criterion == "epsilon"
      n = length(stat)
      rshift = stat[2:end]
      push!(rshift, 0.0)
      lshift = stat[1:(end-1)]
      pushfirst!(lshift, 0.0)
      over =  stat .> D_n #findall(x -> x > D_n , stat ) #indices greater than D_n
      vv = over .& (lshift .< D_n )#& lshift >0 ) #lowers
      v = findall(vv .== 1)
      ww = over .& (rshift .< D_n)
      w = findall(ww .== 1)
      nu_remove = findall(w-v .>= epsilon * G) #nu(epsilon) test for distance between
      v_nu = v[nu_remove]; w_nu = w[nu_remove] #c(w[nu_remove],n)
      sub_pairs = [v_nu w_nu]

   q = size(sub_pairs)[1]
   if  q>0
       for ii in 1:q
           interval = sub_pairs[ii,1]:sub_pairs[ii,2]
           kk = findmax(stat[interval])[2] #internal cp location
           push!(cps, kk + 1 + sub_pairs[ii,1]) #- G-p
       end
   end
 end
   if criterion == "eta"
       n = length(stat)
       window = Int(floor(eta*G))
       for t in (G+1):(n-G)
           if (stat[t] > D_n) & (stat[t] == maximum(stat[(t .- window):(t .+ window)]) )
               push!(cps, t) ##add to cps
           end
       end
   end

  q = length(cps)
  if q>0
      Reject = true
  end

  return cps, q, Reject
end


text:
using DashBootstrapComponents

navbar = dbc_navbarsimple(
    children = [
        dbc_navitem(dbc_navlink("Page 1", href = "#")),
        dbc_dropdownmenu(
            children = [
                dbc_dropdownmenuitem("More pages", header = true),
                dbc_dropdownmenuitem("Page 2", href = "#"),
                dbc_dropdownmenuitem("Page 3", href = "#"),
            ],
            nav = true,
            in_navbar = true,
            label = "More",
        ),
    ],
    brand = "NavbarSimple",
    brand_href = "#",
    color = "primary",
    dark = true,
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module midiModule

using MIDI
function getMidiDatas(filePath)

#midiファイル読み込み
println("koko")
midi = MIDI.readMIDIFile(filePath)
#1トラック目取得
notes = getnotes(midi.tracks[1], midi.tpq)
#ピッチ、発音位置、ベロシティー取得
pitchSequence = []
positionSequence = []
velocitySequence = []
unit = 120
beforeNotePosition = 0
beforeNoteDuration = 0
cnt = 0
for note in notes
    cnt += 1
    push!(pitchSequence, Int(note.pitch))
    push!(positionSequence, Int(note.position))
    # ベロシティーは、休符のところは0、音符が伸びて鳴っている箇所はそのノートが持つベロシティーを、duration=120ごとに1要素として配列で取得
    # 休符がある場合
    if note.position - beforeNotePosition - beforeNoteDuration + 1 > 0 && cnt != 1
        #その幅だけ休符を追加
        append!(velocitySequence, repeat([0], div(note.position - beforeNotePosition - beforeNoteDuration + 1, unit)))
    end
    append!(velocitySequence, repeat([Int(note.velocity)], div(Int(note.duration + 1), unit)))
    beforeNotePosition = Int(note.position)
    beforeNoteDuration = Int(note.duration)
end
return pitchSequence, positionSequence, velocitySequence
end

end


text:
"""
    HANS

HANS (Heauristic Analysis for NLI Systems) is a dataset for NLI.

It contains the set of examples used in the 2019 paper "Right for the
Wrong Reasons: Diagnosing Syntactic Heuristics in Natural Language
Inference" by R. Tom McCoy, Ellie Pavlick, and Tal Linzen.

For details, see the [2019 paper](https://www.aclweb.org/anthology/P19-1334) "Right for the Wrong Reasons:
Diagnosing Syntactic Heuristics in Natural Language Inference"
by by R. Tom McCoy, Ellie Pavlick, and Tal Linzen.

Consists of a set of examples for evaluation, provided with `test_tsv`.

```julia
HANS.test_tsv()
```
"""
module HANS

import ...register_data
using DataDeps: @datadep_str

hansfile(filename) =
    joinpath(datadep"HANS", "hans-0e52769572ab97da1b919baf222f355bdf481d88", filename)

test_tsv() = hansfile("heuristics_evaluation_set.txt")

function __init__()
    register_data(
        "HANS",
        """
        HANS (Heauristic Analysis for NLI Systems) is a dataset for natural language inference.
        It contains the set of examples used in the 2019 paper "Right for the Wrong Reasons: Diagnosing Syntactic Heuristics in Natural Language Inference" by R. Tom McCoy, Ellie Pavlick, and Tal Linzen.
        """,
        "https://github.com/tommccoy1/hans/archive/0e52769572ab97da1b919baf222f355bdf481d88.zip",
        "24316fa3d7298dfd5198059ed1c4eab0b84876cfbc8851df48e96abaa36ead8c")
end

end # module


text:
using AdventOfCode2018
using AdventOfCode2018.Day1

input = get_aoc_input(1)
frequency_shifts = parse_input(input)

println("Solution 1: ", part_1(frequency_shifts))
println("Solution 2: ", part_2(frequency_shifts))


text:
#figuring out matrix
using Random
a = [1,2,3,4,5;6,7,8,9]
for i = 1:length(a)
    print(i)
end

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie