text:
module TestIDRs

using IterativeSolvers
using Test
using Random
using LinearAlgebra
using SparseArrays

@testset "IDR(s)" begin

n = 10
m = 6
Random.seed!(1234567)

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n) + n * I
    b = rand(T, n)
    reltol = √eps(real(T))

    @testset "Without residual smoothing" begin
        x, history = idrs(A, b, reltol=reltol, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ reltol
    end

    # with smoothing
    @testset "With residual smoothing" begin
        x, history = idrs(A, b; smoothing=true, log=true)
        @test history.isconverged
        @test norm(A*x - b) / norm(b) ≤ reltol
    end
end

@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
    A = sprand(T, n, n, 0.5) + n * I
    b = rand(T, n)
    reltol = √eps(real(T))

    x, history = idrs(A, b, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ reltol
end

@testset "Maximum number of iterations" begin
    x, history = idrs(rand(5, 5), rand(5), log=true, maxiter=2)
    @test history.iters == 2
    @test length(history[:resnorm]) == 2
end

@testset "Near solution (#222)" begin
    x = rand(5)
    A = rand(5, 5)
    b = rand(5)

    x, history = idrs!(x, A, b, log=true)
    x_new = copy(x)
    x_new, history = idrs!(x_new, A, b, log=true)

    @test x_new ≈ x
end

@testset "Termination criterion" begin
    for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = T[ 2 -1  0
              -1  2 -1
               0 -1  2]
        n = size(A, 2)
        b = ones(T, n)
        x0 = A \ b
        perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]

        # If the initial residual is small and a small relative tolerance is used,
        # many iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = idrs!(x, A, b, log=true)
        @test 2 ≤ niters(ch) ≤ n

        # If the initial residual is small and a large absolute tolerance is used,
        # no iterations are necessary
        x = x0 + perturbation
        initial_residual = norm(A * x - b)
        x, ch = idrs!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
        @test niters(ch) == 0
    end
end

end

end # module


text:
using Reexport

include("core.jl")
@reexport using .CausalCore

include("semlang.jl")
@reexport using .SEMLang

ext: jl


text:
tests = ["Spin", "sequences"]
for t in tests
  include("$(t).jl")
end


text:
const LANGXIANG = 8

# specialize on slope type
function langxiang(polim::PolarImage, thresh, θ1::Real=0, θ2::Real=pi/2, nϕ::Integer=LANGXIANG)
    checkθ1θ2(θ1,θ2)

    ## WITHOUT SLOPE ##
    if !hasslope(polim)
        segm = segments(polim, θ1, θ2, nϕ)
        return log(mean(gapfraction.(segm, thresh))) / mean(loggapfraction.(segm, thresh))
    end

    ## WITH  SLOPE ##
    # like `segments` but also with cosτsort
    ind_first, ind_last = firstlastind(polim, θ1, θ2)
    segmvec = [eltype(polim)[] for i = 1:nϕ]
    cosτvec = [Float64[] for i = 1:nϕ]
    imgsort = polim.imgsort
    ϕsort = getϕsort(polim)
    cosτ = polim.slope.cosτsort
    adj = nϕ/2π #adjustment to segment ϕsort
    @inbounds for ind in ind_first:ind_last
        indn = ceil(Int, (ϕsort[ind] + pi) * adj)
        push!(segmvec[indn], imgsort[ind])
        push!(cosτvec[indn], cosτ[ind])
    end

    # TODO use weighted average θ instead of simple average
    θ = (θ1 + θ2) / 2
    K = [contactfreqs_iterate(segmvec[i], cosτvec[i], thresh, θ) for i = 1:nϕ]
    T = exp.( -K / cos(θ))
    for i = 1:nϕ
        if T[i] == 0.0
            T[i] = 1 / length(segmvec[i])
        end
    end
    clump_LX = log.(mean(T)) / mean(log.(T))
end

# for Chen Cihlar see 'ChenCihlar.jl'


text:
using ADSeismic
using ADCME
using PyPlot
using PyCall
using DelimitedFiles
using Statistics
use_gpu()

function mkpath_by_force(path)
  if !isdir(path)
    mkpath(path)
  end
end

output_dir = "data/acoustic"
figure_dir = "figure/acoustic/layer-model/"
result_dir = "result/acoustic/layer-model/"

mkpath_by_force(output_dir)
mkpath_by_force(figure_dir)
mkpath_by_force(result_dir)

################### Generate synthetic data #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/layer-model-anomaly.mat"; IT_DISPLAY=0)
src = load_acoustic_source("models/layer-model-anomaly.mat")
rcv_sim = load_acoustic_receiver("models/layer-model-anomaly.mat")

Rs_ = compute_forward_GPU(ap_sim, src, rcv_sim)

sess = Session(); init(sess)
Rs = run(sess, Rs_)

for i = 1:length(src)
    writedlm(joinpath(output_dir, "layermodel-r$i.txt"), Rs[i])
end

u = run(sess, ap_sim(src[div(length(src),2)+1]).u)
visualize_file(u, ap_sim.param, dirs=figure_dir)

################### Inversion using Automatic Differentiation #####################
reset_default_graph()

ap_sim = load_acoustic_model("models/layer-model-smooth.mat"; inv_vp=true, IT_DISPLAY=0)
src = load_acoustic_source("models/layer-model-smooth.mat")
rcv_sim = load_acoustic_receiver("models/layer-model-smooth.mat")

Rs = Array{Array{Float64,2}}(undef, length(src))
for i = 1:length(src)
    Rs[i] = readdlm(joinpath(output_dir, "layermodel-r$i.txt"))
end

losses, gs = compute_loss_and_grads_GPU(ap_sim, src, rcv_sim, Rs, get_collection()[1])
g = sum(gs)
loss = sum(losses)
sess = Session(); init(sess)

vp = (run(sess, ap_sim.vp.contents)).^0.5
scale = mean(vp)
vmax = maximum(vp)
vmin = minimum(vp)

function callback(x, iter, loss)
  if iter%10==0
    X = Array(reshape(x, ap_sim.param.NY+2, ap_sim.param.NX+2))
    close("all")
    pcolormesh([0:ap_sim.param.NX+1;]*ap_sim.param.DELTAX,[0:ap_sim.param.NY+1;]*ap_sim.param.DELTAY,  X*scale, vmax=vmax, vmin=vmin)
    axis("scaled")
    colorbar(shrink=0.4)
    xlabel("x (km)")
    ylabel("z (km)")
    gca().invert_yaxis()
    title("Iteration = $iter")
    savefig(joinpath(figure_dir, "inv_$(lpad(iter,5,"0")).png"))
    writedlm(joinpath(result_dir, "inv_$(lpad(iter,5,"0")).txt"), X*scale)
    open(joinpath(result_dir, "loss.txt"), "a") do io
      writedlm(io, loss)
    end
  else
    return
  end
end

# callback(run(sess, get_collection()[1])', 0, 0)
LBFGS!(sess, loss, g, get_collection()[1]; callback=callback)


text:
#!/usr/bin/env julia
#
# Called with no arguments will build the formattter system image
#   PATH_TO_CLIMATEMACHINE/.git/hooks/JuliaFormatterSysimage.so
#
# Called with a single argument the system image for the formatter will be
# placed in the path specified by the argument (relative to the callers path)

sysimage_path = abspath(
    isempty(ARGS) ?
    joinpath(@__DIR__, "..", ".git", "hooks", "JuliaFormatterSysimage.so") :
    ARGS[1],
)

@info """
Creating system image object file at:
    $(sysimage_path)
"""

using Pkg
Pkg.add("PackageCompiler")
using PackageCompiler

# If a current Manifest exist for the formatter we remove it so that we have the
# latest version
rm(joinpath(@__DIR__, "Manifest.toml"); force = true)

Pkg.activate(joinpath(@__DIR__))
PackageCompiler.create_sysimage(
    :JuliaFormatter;
    precompile_execution_file = joinpath(@__DIR__, "precompile.jl"),
    sysimage_path = sysimage_path,
)


text:
"""
Community detection using the label propagation algorithm (see [Raghavan et al.](http://arxiv.org/abs/0709.2938)).
`g`: input Graph
`maxiter`: maximum number of iterations
return : vertex assignments and the convergence history
"""
function label_propagation(g::SimpleGraph; maxiter=1000)
    n = nv(g)
    label = collect(1:n)
    active_nodes = IntSet(vertices(g))
    c = NeighComm(collect(1:n), fill(-1,n), 1)
    convergence_hist = Vector{Int}()
    random_order = Array(Int, n)
    i = 0
    while !isempty(active_nodes) && i < maxiter
        num_active = length(active_nodes)
        push!(convergence_hist, num_active)
        i += 1

        # processing nodes in random order
        for (j,node) in enumerate(active_nodes)
            random_order[j] = node
        end
        range_shuffle!(1:num_active, random_order)
        @inbounds for j=1:num_active
            u = random_order[j]
            old_comm = label[u]
            label[u] = vote!(g, label, c, u)
            if old_comm != label[u]
                for v in out_neighbors(g, u)
                    push!(active_nodes, v)
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
    fill!(c.neigh_cnt, 0)
    renumber_labels!(label, c.neigh_cnt)
    label, convergence_hist
end

"""Type to record neighbor labels and their counts."""
type NeighComm
  neigh_pos::Vector{Int}
  neigh_cnt::Vector{Int}
  neigh_last::Int
end

"""Fast shuffle Array `a` in UnitRange `r` inplace."""
function range_shuffle!(r::UnitRange, a::AbstractVector)
    (r.start > 0 && r.stop <= length(a)) || error("out of bounds")
    @inbounds for i=length(r):-1:2
        j = rand(1:i)
        ii = i + r.start - 1
        jj = j + r.start - 1
        a[ii],a[jj] = a[jj],a[ii]
    end
end

"""Return the most frequency label."""
function vote!(g::SimpleGraph, m::Vector{Int}, c::NeighComm, u::Int)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in out_neighbors(g,u)
        neigh_comm = m[neigh]
        if c.neigh_cnt[neigh_comm] < 0
            c.neigh_cnt[neigh_comm] = 0
            c.neigh_pos[c.neigh_last] = neigh_comm
            c.neigh_last += 1
        end
        c.neigh_cnt[neigh_comm] += 1
        if c.neigh_cnt[neigh_comm] > max_cnt
          max_cnt = c.neigh_cnt[neigh_comm]
        end
    end
    # ties breaking randomly
    range_shuffle!(1:c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
      if c.neigh_cnt[lbl] == max_cnt
        return lbl
      end
    end
end

function renumber_labels!(membership::Vector{Int}, label_counters::Vector{Int})
    N = length(membership)
    (maximum(membership) > N || minimum(membership) < 1) && error("Label must between 1 and |V|")
    j = 1
    @inbounds for i=1:length(membership)
        k = membership[i]
        if k >= 1
            if label_counters[k] == 0
                # We have seen this label for the first time
                label_counters[k] = j
                k = j
                j += 1
            else
                k = label_counters[k]
            end
        end
        membership[i] = k
    end
end


text:
function maximumtaylor(f, (a, b), n; absmax = false)
    x = setinterval(a, b)
    maybeabs = ifelse(absmax, abs, identity)
    PP = ArbPolyRing(parent(x), :x)

    # Compute the rest term of the Taylor expansion
    y = f(arb_series(PP([x, parent(x)(1)]), n + 1))

    restterm = (x - midpoint(x))^n*y[n]

    # If the rest term is not finite the result always be not finite
    if !isfinite(restterm)
        return restterm
    end

    # Compute the Taylor polynomial at the midpoint of x
    y = f(arb_series(PP([midpoint(x), parent(x)(1)]), n))

    # Enclose the Taylor polynomial evaluated at (a, b) - midpoint((a, b))
    mid = 0.5*(a + b)
    res = maximumpoly(y.poly, (a - mid, b - mid), absmax = absmax)
    p = res

    res += restterm

    return res
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using CairoMakie
CairoMakie.activate!() # HIDE
let
    nonStablePoincare(x,y) = Point(x*(x^2+y^2-1) - y*(x^2+y^2+1),y*(x^2+y^2-1) + x*(x^2+y^2+1))
    stableVanDerPaul(x,y) = Point(y, (1-x^2)*y -x)
    semiStable(x,y) = Point(-y+ x*(-1+x^2+y^2)^2, x+y*(-1+x^2+y^2)^2)

    fig = Figure(resolution = (900, 400), fontsize = 20)
    ax1 = fig[1, 1] = Axis(fig, xlabel = "x", ylabel = "y", title = "non-stable",
    aspect = 1, backgroundcolor = :black)
    ax2 = fig[1, 2] = Axis(fig, xlabel = "x", title = "stable",
    aspect = 1,backgroundcolor = :black)
    ax3 = fig[1, 3] = Axis(fig, xlabel = "x", title = "semi-stable",
    aspect = 1, backgroundcolor = :black)

    streamplot!(ax1, nonStablePoincare, -4..4, -4..4, colormap = Reverse(:plasma),
        gridsize= (32,32), arrow_size = 10)
    streamplot!(ax2, stableVanDerPaul, -4..4, -4..4, colormap = Reverse(:viridis),
        gridsize= (32,32), arrow_size = 10)
    streamplot!(ax3, semiStable, -4..4, -4..4, colormap = Reverse(:inferno),
        gridsize= (32,32), arrow_size = 10)
    hideydecorations!(ax2, grid = false)
    hideydecorations!(ax3, grid = false)
    limits!(ax1, -4, 4, -4, 4)
    limits!(ax2, -4, 4, -4, 4)
    limits!(ax3, -4, 4, -4, 4)
    fig
    #save("streamplots3.png", fig, px_per_unit = 2)
    save(joinpath(@__DIR__, "output", "streamplots3.png"), fig, px_per_unit = 2) # HIDE
end
using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module BenchLookup
using BenchmarkTools
using CompressHashDisplace

suite = BenchmarkGroup()

DICTIONARY = joinpath(@__DIR__, "..", "assets", "american-english.txt")

dict = Dict{String, Int}()

for (line, word) in enumerate(readlines(DICTIONARY))
    dict[word] = line
end

fd = FrozenDict(dict)
fud = FrozenUnsafeDict(dict)

word = "world"
suite["world_fd_lookup"] = @benchmarkable $fd[$word]
suite["world_fud_lookup"] = @benchmarkable $fud[$word]
suite["world_int_lookup"] = @benchmarkable $dict[$word]

word = "hello"
suite["hello_fd_lookup"] = @benchmarkable $fd[$word]
suite["hello_fud_lookup"] = @benchmarkable $fud[$word]
suite["hello_int_lookup"] = @benchmarkable $dict[$word]

word = "quintessential"
suite["quintessential_fd_lookup"] = @benchmarkable $fd[$word]
suite["quintessential_fud_lookup"] = @benchmarkable $fud[$word]
suite["quintessential_int_lookup"] = @benchmarkable $dict[$word]

end # module

BenchLookup.suite


text:
using DictionarySlicing
using Test

@testset "DictionarySlicing.jl" begin

end


text:
# ---
# title: 1146. Snapshot Array
# id: problem1146
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/snapshot-array/description/>
# hidden: true
# ---
#
# Implement a SnapshotArray that supports the following interface:
#
#   * `SnapshotArray(int length)` initializes an array-like data structure with the given length.  **Initially, each element equals 0**.
#   * `void set(index, val)` sets the element at the given `index` to be equal to `val`.
#   * `int snap()` takes a snapshot of the array and returns the `snap_id`: the total number of times we called `snap()` minus `1`.
#   * `int get(index, snap_id)` returns the value at the given `index`, at the time we took the snapshot with the given `snap_id`
#
#
#
# **Example 1:**
#
#
#
#     Input: ["SnapshotArray","set","snap","set","get"]
#     [[3],[0,5],[],[0,6],[0,0]]
#     Output: [null,null,0,null,5]
#     Explanation:
#     SnapshotArray snapshotArr = new SnapshotArray(3); // set the length to be 3
#     snapshotArr.set(0,5);  // Set array[0] = 5
#     snapshotArr.snap();  // Take a snapshot, return snap_id = 0
#     snapshotArr.set(0,6);
#     snapshotArr.get(0,0);  // Get the value of array[0] with snap_id = 0, return 5
#
#
#
# **Constraints:**
#
#   * `1 <= length <= 50000`
#   * At most `50000` calls will be made to `set`, `snap`, and `get`.
#   * `0 <= index < length`
#   * `0 <= snap_id < `(the total number of times we call `snap()`)
#   * `0 <= val <= 10^9`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
immutable Uniform <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64
    function Uniform(a::Real, b::Real)
	a < b || error("a < b required for range [a, b]")
	new(float64(a), float64(b))
    end
    Uniform() = new(0.0, 1.0)
end

## Support
@distr_support Uniform d.a d.b

## Properties
mean(d::Uniform) = (d.a + d.b) / 2.0

median(d::Uniform) = (d.a + d.b) / 2.0

mode(d::Uniform) = mean(d)
modes(d::Uniform) = Float64[]

function var(d::Uniform)
    w = d.b - d.a
    w * w / 12.0
end

skewness(d::Uniform) = 0.0
kurtosis(d::Uniform) = -1.2

entropy(d::Uniform) = log(d.b - d.a)

## Functions
pdf(d::Uniform, x::Float64) = insupport(d,x) ? 1/(d.b-d.a) : 0.0
logpdf(d::Uniform, x::Float64) = insupport(d,x) ? -log(d.b-d.a) : -Inf

function cdf(d::Uniform, q::Float64)
    if isnan(q)
        return NaN
    elseif q <= d.a
        return 0.0
    elseif q >= d.b
        return 1.0
    end
    (q-d.a)/(d.b-d.a)
end
function ccdf(d::Uniform, q::Float64)
    if isnan(q)
        return NaN
    elseif q <= d.a
        return 1.0
    elseif q >= d.b
        return 0.0
    end
    (d.b-q)/(d.b-d.a)
end

quantile(d::Uniform, p::Float64) = @checkquantile p d.a+p*(d.b-d.a)
cquantile(d::Uniform, p::Float64) = @checkquantile p d.b+p*(d.a-d.b)

function mgf(d::Uniform, t::Real)
    u = 0.5*(d.b-d.a)*t
    u == zero(u) && return one(u)
    v = 0.5*(d.a+d.b)*t
    exp(v)*(sinh(u)/u)
end
function cf(d::Uniform, t::Real)
    u = 0.5*(d.b-d.a)*t
    u == zero(u) && return complex(one(u))
    v = 0.5*(d.a+d.b)*t
    exp(im*v)*(sin(u)/u)
end

## Sampling
rand(d::Uniform) = d.a + (d.b - d.a) * rand()

## Fitting
function fit_mle{T <: Real}(::Type{Uniform}, x::Vector{T})
    if isempty(x)
        throw(ArgumentError("x cannot be empty."))
    end

    xmin = xmax = x[1]
    for i = 2:length(x)
        xi = x[i]
        if xi < xmin
            xmin = xi
        elseif xi > xmax
            xmax = xi
        end
    end

    Uniform(xmin, xmax)
end


text:
using BenchmarkTools
using UltraDark: Grids, max_phase_grad

resol = 64

grids = Grids(1.0, resol)

@benchmark max_phase_grad(grids.ψx)


text:
module QuadraticModelsCPLEX

using CPLEX
using QuadraticModels, SolverCore
using LinearAlgebra, SparseArrays

export cplex, presolve_to_file, CPXPresolveException

const cplex_statuses = Dict(1 => :acceptable,
                            2 => :unbounded,
                            3 => :infeasible,
                            4 => :infeasible,
                            10 => :max_iter,
                            11 => :max_time,
                            12 => :exception,
                            13 => :user)

function sparse_csr(I, J, V, m=maximum(I), n=maximum(J))
  csrrowptr = zeros(Int, m+1)
  # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr
  coolen = length(I)
  min(length(J), length(V)) >= coolen || throw(ArgumentError("J and V need length >= length(I) = $coolen"))
  @inbounds for k in 1:coolen
    Ik = I[k]
    if 1 > Ik || m < Ik
      throw(ArgumentError("row indices I[k] must satisfy 1 <= I[k] <= m"))
    end
    csrrowptr[Ik+1] += 1
  end

  # Compute the CSR form's rowptrs and store them shifted forward by one in csrrowptr
  countsum = 1
  csrrowptr[1] = 1
  @inbounds for i in 2:(m+1)
    overwritten = csrrowptr[i]
    csrrowptr[i] = countsum
    countsum += overwritten
  end

  # Counting-sort the column and nonzero values from J and V into csrcolval and csrnzval
  # Tracking write positions in csrrowptr corrects the row pointers
  csrcolval = zeros(Int, length(I))
  csrnzval = zeros(length(I))
  @inbounds for k in 1:coolen
    Ik, Jk = I[k], J[k]
    if 1 > Jk || n < Jk
        throw(ArgumentError("column indices J[k] must satisfy 1 <= J[k] <= n"))
    end
    csrk = csrrowptr[Ik+1]
    csrrowptr[Ik+1] = csrk + 1
    csrcolval[csrk] = Jk
    csrnzval[csrk] = V[k]
  end
  csrrowptr = csrrowptr[1:end-1]
  return csrrowptr, csrcolval, csrnzval
end

function cplex_inputQM(QM::QuadraticModel; method=4, display=1, kwargs...)

  env = CPLEX.Env()
  CPXsetintparam(env, CPXPARAM_ScreenOutput, display)   # Enable output (0=off)
  CPXsetdblparam(env, CPXPARAM_TimeLimit, 3600)  # Time limit
  CPXsetintparam(env, CPXPARAM_Threads, 1) # Single thread
  for (k, v) in kwargs
    if k==:presolve
      CPXsetintparam(env, CPXPARAM_Preprocessing_Presolve, v) # 0 = no presolve
    elseif k==:scaling
      CPXsetintparam(env, CPXPARAM_Read_Scale, -1) # -1 = no scaling
    elseif k==:crossover
      CPXsetintparam(env, CPXPARAM_SolutionType, v)  # 2 = no crossover
    end
  end
  CPXsetintparam(env, CPXPARAM_LPMethod, method)  # 4 = Use barrier
  CPXsetintparam(env, CPXPARAM_QPMethod, method) # 4 = Use barrier, 0 = automatic

  status_p = Ref{Cint}()
  lp = CPXcreateprob(env, status_p, "")
  CPXnewcols(env, lp, QM.meta.nvar, QM.data.c, QM.meta.lvar, QM.meta.uvar, C_NULL, C_NULL)
  CPXchgobjoffset(env, lp, QM.data.c0)
  if QM.meta.nnzh > 0
    Hvals = zeros(eltype(QM.data.Hvals), length(QM.data.Hvals))
    for i=1:length(QM.data.Hvals)
      if QM.data.Hrows[i] == QM.data.Hcols[i]
        Hvals[i] = QM.data.Hvals[i] / 2
      else
        Hvals[i] = QM.data.Hvals[i]
      end
    end
    Q = sparse(QM.data.Hrows, QM.data.Hcols, QM.data.Hvals, QM.meta.nvar, QM.meta.nvar)
    diag_matrix = spdiagm(0 => diag(Q))
    Q = Q + Q' - diag_matrix
    qmatcnt = zeros(Int, QM.meta.nvar)
    for k = 1:QM.meta.nvar
      qmatcnt[k] = Q.colptr[k+1] - Q.colptr[k]
    end
    CPXcopyquad(env, lp, convert(Array{Cint,1}, Q.colptr[1:end-1].-1), convert(Array{Cint,1}, qmatcnt),
                convert(Array{Cint,1}, Q.rowval.-1), Q.nzval)
  end

  Acsrrowptr, Acsrcolval, Acsrnzval = sparse_csr(QM.data.Arows, QM.data.Acols,
                                                  QM.data.Avals, QM.meta.ncon,
                                                  QM.meta.nvar)

  sense = fill(Cchar('A'), QM.meta.ncon) # lower, greater, range or equal. A is for the init
  if length(QM.meta.jinf) > 0
    error("infeasible bounds in A")
  end
  p_low, p_upp, p_rng, p_fix = 1, 1, 1, 1
  for j=1:QM.meta.ncon
    if length(QM.meta.jlow) > 0 && QM.meta.jlow[p_low] == j
      sense[j] = Cchar('G')
      if (p_low < length(QM.meta.jlow)) p_low += 1 end
    elseif length(QM.meta.jupp) > 0 && QM.meta.jupp[p_upp] == j
      sense[j] = Cchar('L')
      if (p_upp < length(QM.meta.jupp)) p_upp += 1 end
    elseif length(QM.meta.jrng) > 0 && QM.meta.jrng[p_rng] == j
      sense[j] = Cchar('R')
      if (p_rng < length(QM.meta.jrng)) p_rng += 1 end
    elseif length(QM.meta.jfix) > 0 && QM.meta.jfix[p_fix] == j
      sense[j] = Cchar('E')
      if (p_fix < length(QM.meta.jfix)) p_fix += 1 end
    else
      error("A error")
    end
  end
  rhs = zeros(QM.meta.ncon)
  drange = zeros(QM.meta.ncon)
  for j = 1:QM.meta.ncon
    if QM.meta.lcon[j] != -Inf && QM.meta.ucon[j] != Inf
      rhs[j] = QM.meta.ucon[j]
      drange[j] = QM.meta.ucon[j] - QM.meta.lcon[j]
    elseif QM.meta.lcon[j] != -Inf && QM.meta.ucon[j] == Inf
      rhs[j] = QM.meta.lcon[j]
    elseif QM.meta.lcon[j] == -Inf && QM.meta.ucon[j] != Inf
      rhs[j] = QM.meta.ucon[j]
    else
      rhs[j] = Inf
    end
  end
  CPXaddrows(env, lp, 0, QM.meta.ncon, length(Acsrcolval), rhs,
             sense, convert(Vector{Cint}, Acsrrowptr.- Cint(1)), convert(Vector{Cint}, Acsrcolval.- Cint(1)),
             Acsrnzval, C_NULL, C_NULL)

  return env, lp
end

function cplex(QM::QuadraticModel; method=4, display=1, kwargs...)

  env, lp = cplex_inputQM(QM; method=method, display=display, kwargs...)

  t = @timed begin
    if QM.meta.nnzh > 0
      CPXqpopt(env, lp)
    else
      CPXlpopt(env, lp)
    end
  end

  x = Vector{Cdouble}(undef, QM.meta.nvar)
  CPXgetx(env, lp, x, 0, QM.meta.nvar-1)
  y = Vector{Cdouble}(undef, QM.meta.ncon)
  CPXgetpi(env, lp, y, 0, QM.meta.ncon-1)
  s = Vector{Cdouble}(undef, QM.meta.nvar)
  CPXgetdj(env, lp, s, 0, QM.meta.nvar-1)
  primal_feas = Vector{Cdouble}(undef, 1)
  CPXgetdblquality(env, lp, primal_feas, CPX_MAX_PRIMAL_RESIDUAL)
  dual_feas = Vector{Cdouble}(undef, 1)
  CPXgetdblquality(env, lp, dual_feas, CPX_MAX_DUAL_RESIDUAL)
  objval_p = Vector{Cdouble}(undef, 1)
  CPXgetobjval(env, lp, objval_p)

  stats = GenericExecutionStats(get(cplex_statuses, CPXgetstat(env, lp), :unknown),
                                QM, solution = x,
                                objective = objval_p[1],
                                primal_feas = primal_feas[1],
                                dual_feas = dual_feas[1],
                                iter = Int64(CPXgetbaritcnt(env, lp)),
                                multipliers = y,
                                elapsed_time = t[2])
  return stats
end

mutable struct CPXPresolveException <: Exception
  msg::AbstractString
end

function presolve_to_file(QM::QuadraticModel; path::String = "presolved.mps", method=4, display=1, kwargs...)

  env, lp = cplex_inputQM(QM; method=method, display=display, kwargs...)
  objoff = Ref{Cdouble}()
  str_p =  string(path[1:end-3], "pre")
  status = CPXpreslvwrite(env, lp, str_p, objoff)
  status == 1103 && throw(CPXPresolveException("No presolve operation performed"))

  status_p = Ref{Cint}()
  lp_p = CPXcreateprob(env, status_p, "")

  CPXreadcopyprob(env, lp_p, str_p, C_NULL)
  status = CPXwriteprob(env, lp_p, path, C_NULL)
end

end


text:
using LinearAlgebra, DiffEqBase, ForwardDiff, StaticArrays
using SparseArrays

export dimension, get_state, DynamicalSystem
export ContinuousDynamicalSystem, DiscreteDynamicalSystem
export set_parameter!, step!
export trajectory, jacobian
export integrator, tangent_integrator, parallel_integrator
export set_state!, get_state, get_deviations, set_deviations!

#####################################################################################
#                                  DynamicalSystem                                  #
#####################################################################################
"""
    DynamicalSystem

The central structure of **DynamicalSystems.jl**. All functions of the suite that
can use known equations of motion expect an instance of this type.

## Constructing a `DynamicalSystem`
```julia
DiscreteDynamicalSystem(eom, state, p [, jacobian [, J0]]; t0::Int = 0)
ContinuousDynamicalSystem(eom, state, p [, jacobian [, J0]]; t0 = 0.0)
```
with `eom` the equations of motion function (see below).
`p` is a parameter container, which we highly suggest to use a mutable object like
`Array`, [`LMArray`](https://github.com/JuliaDiffEq/LabelledArrays.jl) or
a dictionary. Pass `nothing` in the place of `p` if your system does not have
parameters.

`t0`, `J0` allow you to choose the initial time and provide
an initialized Jacobian matrix. See `CDS_KWARGS` for the
default options used to evolve continuous systems (through `OrdinaryDiffEq`).

### Equations of motion
The are two "versions" for `DynamicalSystem`, depending on whether the
equations of motion (`eom`) are in-place (iip) or out-of-place (oop).
Here is how to define them:

* **oop** : The `eom` **must** be in the form `eom(x, p, t) -> SVector`
  which means that given a state `x::SVector` and some parameter container
  `p` it returns an [`SVector`](http://juliaarrays.github.io/StaticArrays.jl/stable/pages/api.html#SVector-1)
  (from the [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) module)
  containing the next state.
* **iip** : The `eom` **must** be in the form `eom!(xnew, x, p, t)`
  which means that given a state `x::Vector` and some parameter container `p`,
  it writes in-place the new state in `xnew`.

`t` stands for time (integer for discrete systems).
iip is suggested for big systems, whereas oop is suggested for small systems.
The break-even point at around 100 dimensions, and for using functions that use the
tangent space (like e.g. `lyapunovs` or `gali`), the break-even
point is at around 10 dimensions.

The constructor deduces automatically whether `eom` is iip or oop. It is not
possible however to deduce whether the system is continuous or discrete just from the
equations of motion, hence the 2 constructors.

### Jacobian
The optional argument `jacobian` for the constructors
is a *function* and (if given) must also be of the same form as the `eom`,
`jacobian(x, p, n) -> SMatrix`
for the out-of-place version and `jacobian!(xnew, x, p, n)` for the in-place version.

If `jacobian` is not given, it is constructed automatically using
the module [`ForwardDiff`](http://www.juliadiff.org/ForwardDiff.jl/stable/).
Even though `ForwardDiff` is very fast, depending on your exact system you might
gain significant speed-up by providing a hand-coded Jacobian and so we recommend it.

### Interface to DifferentialEquations.jl
Continuous systems are solved using
[**DifferentialEquations.jl**](http://docs.juliadiffeq.org/latest/).
The following two interfaces are provided:
```
ContinuousDynamicalSystem(prob::ODEProblem [, jacobian [, J0]])
ODEProblem(continuous_dynamical_system, tspan, args...)
```
where in the second case `args` stands for the
standard extra arguments of `ODEProblem`: `callback, mass_matrix`.

If you want to use callbacks with [`tangent_integrator`](@ref) or
[`parallel_integrator`](@ref), then invoke them with extra arguments
as shown in the [Advanced Documentation](https://juliadynamics.github.io/DynamicalSystems.jl/latest/advanced/).

## Relevant Functions
[`trajectory`](@ref), [`set_parameter!`](@ref).
"""
abstract type DynamicalSystem{
        IIP,     # is in place , for dispatch purposes and clarity
        S,       # state type
        D,       # dimension
        F,       # equations of motion
        P,       # parameters
        JAC,     # jacobian
        JM,      # jacobian matrix
        IAD}     # is auto-differentiated
    # one-liner: {IIP, S, D, F, P, JAC, JM, IAD}
    # Subtypes of DynamicalSystem have fields:
    # 1. f
    # 2. u0
    # 3. p
    # 4. t0
    # 5. jacobian (function)
    # 6. J (matrix)
end
const DS = DynamicalSystem

"""
    ContinuousDynamicalSystem(eom, state, p [, jacobian [, J]]; t0 = 0.0)
    ContinuousDynamicalSystem(prob::ODEProblem [, jacobian [, J0]])
A `DynamicalSystem` restricted to continuous-time systems (also called *ODEs*).
"""
struct ContinuousDynamicalSystem{IIP, S, D, F, P, JAC, JM, IAD} <:
                 DynamicalSystem{IIP, S, D, F, P, JAC, JM, IAD}
    f::F
    u0::S
    p::P
    t0::eltype(S)
    jacobian::JAC
    J::JM
end
const CDS = ContinuousDynamicalSystem
systemtype(::CDS) = "continuous"

"""
    DiscreteDynamicalSystem(eom, state, p [, jacobian [, J]]; t0::Int = 0)
A `DynamicalSystem` restricted to discrete-time systems (also called *maps*).
"""
struct DiscreteDynamicalSystem{IIP, S, D, F, P, JAC, JM, IAD} <:
               DynamicalSystem{IIP, S, D, F, P, JAC, JM, IAD}
    f::F
    u0::S
    p::P
    t0::Int
    jacobian::JAC
    J::JM
end
const DDS = DiscreteDynamicalSystem
systemtype(::DDS) = "discrete"

DiffEqBase.isinplace(::DS{IIP}) where {IIP} = IIP
statetype(::DS{IIP, S}) where {IIP, S} = S
stateeltype(::DS{IIP, S}) where {IIP, S} = eltype(S)
isautodiff(::DS{IIP, S, D, F, P, JAC, JM, IAD}) where
{IIP, S, D, F, P, JAC, JM, IAD} = IAD

get_state(ds::DS) = ds.u0

"""
    dimension(thing) -> D
Return the dimension of the `thing`, in the sense of state-space dimensionality.
"""
dimension(ds::DS{IIP, S, D}) where {IIP, S, D} = D

"""
    set_parameter!(ds::DynamicalSystem, index, value)
    set_parameter!(ds::DynamicalSystem, values)
Change one or many parameters of the system
by setting `p[index] = value` in the first case
and `p .= values` in the second.

The same function also works for any integrator.
"""
set_parameter!(ds, index, value) = (ds.p[index] = value)
set_parameter!(ds, values) = (ds.p .= values)

#####################################################################################
#                           State types enforcing                                   #
#####################################################################################
safe_state_type(::Val{true}, u0) = Vector(u0)
safe_state_type(::Val{false}, u0) = SVector{length(u0)}(u0...)
safe_state_type(::Val{false}, u0::SVector) = u0
safe_state_type(::Val{false}, u0::Number) = u0

safe_matrix_type(::Val{true}, Q::Matrix) = Q
safe_matrix_type(::Val{true}, Q::AbstractMatrix) = Matrix(Q)
function safe_matrix_type(::Val{false}, Q::AbstractMatrix)
    A, B = size(Q)
    SMatrix{A, B}(Q)
end
save_matrix_type(::Val{false}, Q::SMatrix) = Q
safe_matrix_type(_, a::Number) = a

#####################################################################################
#                                    Constructors                                   #
#####################################################################################
for ds in (:ContinuousDynamicalSystem, :DiscreteDynamicalSystem)

    if ds == :ContinuousDynamicalSystem
        @eval timetype(::Type{$(ds)}, s) = eltype(s)
    else
        @eval timetype(::Type{$(ds)}, s) = Int
    end

    # Main Constructor (with Jacobian and j0)
    @eval function $(ds)(
        eom::F, u0, p::P, j::JAC, j0 = nothing;
        t0 = 0.0, IAD = false, IIP = isinplace(eom, 4)) where {F, P, JAC}

        if !(typeof(u0) <: Union{AbstractVector, Number})
            throw(ArgumentError(
            "The state of a dynamical system *must* be <: AbstractVector/Number!"))
        end

        s = safe_state_type(Val{IIP}(), u0)
        S = typeof(s)
        D = length(s)

        IIP || typeof(eom(s, p, t0)) <: Union{SVector, Number} ||
        throw(ArgumentError("Equations of motion must return an `SVector` "*
        "or number for out-of-place form!"))

        J = j0 != nothing ? j0 : get_J(j, s, p, timetype($(ds), s)(t0), IIP)
        JM = typeof(J)

        return $(ds){
            IIP, S, D, F, P, JAC, JM, IAD}(eom, u0, p, t0, j, J)
    end

    # Without Jacobian
    @eval function $(ds)(eom, u0, p; t0 = 0.0)

        if !(typeof(u0) <: Union{AbstractVector, Number})
            throw(ArgumentError(
            "The state of a dynamical system *must* be <: AbstractVector/Number!"))
        end

        IIP = isinplace(eom, 4)
        s = safe_state_type(Val{IIP}(), u0)
        D = length(s)

        t = timetype($(ds), s)(t0)
        j = create_jacobian(eom, Val{IIP}(), s, p, t, Val{D}())
        J = get_J(j, s, p, t, IIP)

        return $(ds)(eom, s, p, j, J; t0 = t0, IAD = true, IIP = IIP)
    end
end

#####################################################################################
#                                Pretty-Printing                                    #
#####################################################################################
Base.summary(ds::DS) =
"$(dimension(ds))-dimensional "*systemtype(ds)*" dynamical system"

jacobianstring(ds::DS) = isautodiff(ds) ? "ForwardDiff" : "$(eomstring(ds.jacobian))"
eomstring(f::Function) = nameof(f)
eomstring(f) = nameof(typeof(f))

paramname(p::AbstractArray) = string(p)
paramname(p::Nothing) = repr(p)
paramname(p) = nameof(typeof(p))

# Credit to Sebastian Pfitzner
function printlimited(io, x; Δx = 0, Δy = 0)
    sz = displaysize(io)
    io2 = IOBuffer(); ctx = IOContext(io2, :limit => true, :compact => true,
    :displaysize => (sz[1]-Δy, sz[2]-Δx))
    Base.print_array(ctx, x)
    s = String(take!(io2))
    s = replace(s[2:end], "  " => ", ")
    Base.print(io, "["*s*"]")
end

printlimited(io, x::Number; Δx = 0, Δy = 0) = print(io, x)

function Base.show(io::IO, ds::DS)
    ps = 14
    text = summary(ds)
    u0 = get_state(ds)'

    ctx = IOContext(io, :limit => true, :compact => true, :displaysize => (10,50))

    println(io, text)
    prefix = rpad(" state: ", ps)
    print(io, prefix); printlimited(io, u0, Δx = length(prefix)); print(io, "\n")
    println(io,  rpad(" e.o.m.: ", ps),     eomstring(ds.f))
    println(io,  rpad(" in-place? ", ps),   isinplace(ds))
    println(io,  rpad(" jacobian: ", ps),   jacobianstring(ds)),
    print(io,    rpad(" parameters: ", ps), paramname(ds.p))
end

#######################################################################################
#                                    Jacobians                                        #
#######################################################################################
function create_jacobian(
    @nospecialize(f::F), ::Val{IIP}, s::S, p::P, t::T, ::Val{D}) where {F, IIP, S, P, T, D}
    if IIP
        dum = deepcopy(s)
        cfg = ForwardDiff.JacobianConfig(
        (y, x) -> f(y, x, p, t), dum, s)
        jac = (J, u, p, t) ->
        ForwardDiff.jacobian!(J, (y, x) -> f(y, x, p, t),
        dum, u, cfg, Val{false}())
        return jac
    else
        if D == 1
            return jac = (u, p, t) -> ForwardDiff.derivative((x) -> f(x, p, t), u)
        else
            # SVector methods do *not* use the config
            # cfg = ForwardDiff.JacobianConfig(
            #     (x) -> prob.f(x, prob.p, prob.tspan[1]), prob.u0)
            return jac = (u, p, t) ->
            ForwardDiff.jacobian((x) -> f(x, p, t), u #=, cfg=#)
        end
    end
end

function get_J(jacob, u0, p, t0, iip) where {JAC}
    D = length(u0)
    if iip
        J = similar(u0, (D,D))
        jacob(J, u0, p, t0)
    else
        J = jacob(u0, p, t0)
    end
    return J
end

"""
    jacobian(ds::DynamicalSystem, u = ds.u0, t = ds.t0)
Return the jacobian of the system at `u`, at `t`.
"""
function jacobian(ds::DS{true}, u = ds.u0, t = ds.t0)
    J = similar(ds.J)
    ds.jacobian(J, u, ds.p, t)
    return J
end
jacobian(ds::DS{false}, u = ds.u0, t = ds.t0) =
ds.jacobian(u, ds.p, t)

#######################################################################################
#                        Tangent Dynamics (aka linearized dynamics)                   #
#######################################################################################

# IIP Tangent Space dynamics
function create_tangent(@nospecialize(f::F), @nospecialize(jacobian::JAC), J::JM,
    ::Val{true}, ::Val{k}) where {F, JAC, JM, k}
    J = deepcopy(J)
    tangentf = (du, u, p, t) -> begin
        uv = @view u[:, 1]
        f(view(du, :, 1), uv, p, t)
        jacobian(J, uv, p, t)
        mul!((@view du[:, 2:(k+1)]), J, (@view u[:, 2:(k+1)]))
        nothing
    end
    return tangentf
end
# for the case of autodiffed systems, a specialized version is created
# so that f! is not called twice in ForwardDiff
function create_tangent_iad(f::F, J::JM, u, p, t, ::Val{k}) where {F, JM, k}
    let
        J = deepcopy(J)
        cfg = ForwardDiff.JacobianConfig(
            (du, u) -> f(du, u, p, p), deepcopy(u), deepcopy(u)
        )
        tangentf = (du, u, p, t) -> begin
            uv = @view u[:, 1]
            ForwardDiff.jacobian!(
                J, (du, u) -> f(du, u, p, t), view(du, :, 1), uv, cfg, Val{false}()
            )
            mul!((@view du[:, 2:k+1]), J, (@view u[:, 2:k+1]))
            nothing
        end
        return tangentf
    end
end

# OOP Tangent Space dynamics
function create_tangent(f::F, jacobian::JAC, J::JM,
    ::Val{false}, ::Val{k}) where {F, JAC, JM, k}

    ws_index = SVector{k, Int}(2:(k+1)...)
    tangentf = TangentOOP{F, JAC, k}(f, jacobian, ws_index)
    return tangentf
end
struct TangentOOP{F, JAC, k} <: Function
    f::F
    jacobian::JAC
    ws::SVector{k, Int}
end
function (tan::TangentOOP)(u, p, t)
    @inbounds s = u[:, 1]
    du = tan.f(s, p, t)
    J = tan.jacobian(s, p, t)
    @inbounds dW = J*u[:, tan.ws]
    return hcat(du, dW)
end

#######################################################################################
#                                Parallel Dynamics                                    #
#######################################################################################
# Create equations of motion of evolving states in parallel
function create_parallel(ds::DS{true}, states)
    st = [Vector(s) for s in states]
    L = length(st)
    paralleleom = (du, u, p, t) -> begin
        for i in 1:L
            @inbounds ds.f(du[i], u[i], p, t)
        end
    end
    return paralleleom, st
end

function create_parallel(ds::DS{false}, states)
    D = dimension(ds)
    st = [SVector{D}(s) for s in states]
    L = length(st)
    # The following may be inneficient
    paralleleom = (du, u, p, t) -> begin
        for i in 1:L
            @inbounds du[i] = ds.f(u[i], p, t)
        end
    end
    return paralleleom, st
end

#######################################################################################
#                                     Docstrings                                      #
#######################################################################################
"""
    integrator(ds::DynamicalSystem [, u0]; diffeq...) -> integ
Return an integrator object that can be used to evolve a system interactively
using `step!(integ [, Δt])`. Optionally specify an initial state `u0`.

The state of this integrator is a vector.

* `diffeq...` are keyword arguments propagated into `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.
"""
function integrator end

"""
    get_state(ds::DynamicalSystem)
Return the state of `ds`.

    get_state(integ [, i::Int = 1])
Return the state of the integrator, in the sense of the state of the dynamical system.

If the integrator is a [`parallel_integrator`](@ref), passing `i` will return
the `i`-th state. The function also correctly returns the true state of the system
for tangent integrators.
"""
get_state(integ) = integ.u

"""
    set_state!(integ, u [, i::Int = 1])
Set the state of the integrator to `u`, in the sense of the state of the
dynamical system. Works for any integrator (normal, tangent, parallel).

For parallel integrator, you can choose which state to set (using `i`).

Automatically does `u_modified!(integ, true)`.
"""
set_state!(integ, u) = (integ.u = u; u_modified!(integ, true))

"""
    tangent_integrator(ds::DynamicalSystem, Q0 | k::Int; kwargs...)
Return an integrator object that evolves in parallel both the system as well
as deviation vectors living on the tangent space, also called linearized space.

`Q0` is a *matrix* whose columns are initial values for deviation vectors. If
instead of a matrix `Q0` an integer `k` is given, then `k` random orthonormal
vectors are choosen as initial conditions.

## Keyword Arguments
* `u0` : Optional different initial state.
* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.
  These keywords can also include `callback` for [event handling](http://docs.juliadiffeq.org/latest/features/callback_functions.html).

It is *heavily* advised to use the functions [`get_state`](@ref), [`get_deviations`](@ref),
[`set_state!`](@ref), [`set_deviations!`](@ref) to manipulate the integrator.

## Description

If ``J`` is the jacobian of the system then the *tangent dynamics* are
the equations that evolve in parallel the system as well as
a deviation vector (or matrix) ``w``:
```math
\\begin{aligned}
\\dot{u} &= f(u, p, t) \\\\
\\dot{w} &= J(u, p, t) \\times w
\\end{aligned}
```
with ``f`` being the equations of motion and ``u`` the system state.
Similar equations hold for the discrete case.
"""
function tangent_integrator end

"""
    get_deviations(tang_integ)
Return the deviation vectors of the [`tangent_integrator`](@ref) in a form
of a matrix with columns the vectors.
"""
function get_deviations end

"""
    set_deviations!(tang_integ, Q)
Set the deviation vectors of the [`tangent_integrator`](@ref) to `Q`, which must
be a matrix with each column being a deviation vector.

Automatically does `u_modified!(tang_integ, true)`.
"""
function set_deviations! end

"""
    parallel_integrator(ds::DynamicalSystem, states; kwargs...)
Return an integrator object that can be used to evolve many `states` of
a system in parallel at the *exact same times*, using `step!(integ [, Δt])`.

`states` are expected as vectors of vectors.

## Keyword Arguments
* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.
  These keywords can also include `callback` for [event handling](http://docs.juliadiffeq.org/latest/features/callback_functions.html).

It is *heavily* advised to use the functions [`get_state`](@ref) and
[`set_state!`](@ref) to manipulate the integrator. Provide `i` as a second
argument to change the `i`-th state.
"""
function parallel_integrator end

"""
    trajectory(ds::DynamicalSystem, T [, u]; kwargs...) -> dataset

Return a dataset that will contain the trajectory of the system,
after evolving it for total time `T`, optionally starting from state `u`.
See [`Dataset`](@ref) for info on how to use this object.

A `W×D` dataset is returned, with `W = length(t0:dt:T)` with
`t0:dt:T` representing the time vector (*not* returned) and `D` the system dimension.
For discrete systems both `T` and `dt` must be integers.

## Keyword Arguments
* `dt` :  Time step of value output during the solving
  of the continuous system. For discrete systems it must be an integer. Defaults
  to `0.01` for continuous and `1` for discrete.
* `Ttr` : Transient time to evolve the initial state before starting saving states.
* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
  For example `abstol = 1e-9`.  Only valid for continuous systems.
  If you want to specify a solver, do so by using the name `alg`, e.g.:
  `alg = Tsit5(), maxiters = 1000`. This requires you to have been first
  `using OrdinaryDiffEq` to access the solvers. See
  `DynamicalSystemsBase.CDS_KWARGS` for default values.
  These keywords can also include `callback` for [event handling](http://docs.juliadiffeq.org/latest/features/callback_functions.html).
  Using a `SavingCallback` with `trajectory` will lead to unexpected behavior!

"""
function trajectory end


text:
function FESpace(t::Triangulation,reffes;trian=nothing, kwargs...)
  @assert trian === nothing
  model = get_active_model(t)
  FESpace(model,reffes;trian=t,kwargs...)
end

function FESpace(
  model::DiscreteModel,
  cell_fe::CellFE;
  trian = Triangulation(model),
  labels = get_face_labeling(model),
  dirichlet_tags=Int[],
  dirichlet_masks=nothing,
  constraint=nothing,
  vector_type=nothing)

  @assert num_cells(cell_fe) == num_cells(model) """\n
  The number of cells provided in the `cell_fe` argument ($(cell_fe.num_cells) cells)
  does not match the number of cells ($(num_cells(model)) cells) in the provided DiscreteModel.
  """
  _vector_type = _get_vector_type(vector_type,cell_fe,trian)
  F = _ConformingFESpace(
      _vector_type,
      model,
      labels,
      cell_fe,
      dirichlet_tags,
      dirichlet_masks,
      trian)
  V = _add_constraint(F,cell_fe.max_order,constraint)
  V
end

function _get_vector_type(vector_type,cell_fe,trian)
  if vector_type == nothing
    cell_shapefuns, cell_dof_basis = compute_cell_space(cell_fe,trian)
    T = get_dof_value_type(cell_shapefuns,cell_dof_basis)
    _vector_type = Vector{T}
  else
    @assert vector_type <: AbstractVector """\n
    The (optional) argument vector_type has to be <: AbstractVector.
    """
    _vector_type = vector_type
  end
  _vector_type
end

function _add_constraint(F,order,constraint)
  if constraint == nothing
    V = F
  elseif constraint == :zeromean
    trian = get_triangulation(F)
    dΩ = Measure(trian,order)
    V = ZeroMeanFESpace(F,dΩ)
  else
    @unreachable """\n
    The passed option constraint=$constraint is not valid.
    Valid values for constraint: nothing, :zeromean
    """
  end
  V
end

function FESpace(
  t::Triangulation,
  cell_reffe::AbstractArray{<:ReferenceFE};
  trian=nothing,
  kwargs...)
  @assert trian === nothing
  # TODO for L2 conformity and no dirichlet conditions
  # no needed to build the active model
  model = get_active_model(t)
  FESpace(model,cell_reffe;trian=t,kwargs...)
end

function FESpace(
  model::DiscreteModel,
  cell_reffe::AbstractArray{<:ReferenceFE};
  conformity=nothing,
  trian = Triangulation(model),
  labels = get_face_labeling(model),
  dirichlet_tags=Int[],
  dirichlet_masks=nothing,
  constraint=nothing,
  vector_type=nothing)

  conf = Conformity(testitem(cell_reffe),conformity)

  if _use_clagrangian(trian,cell_reffe,conf) &&
    constraint === nothing &&
    num_vertices(model) == num_nodes(model)

    V = _unsafe_clagrangian(
      cell_reffe,
      Triangulation(model),
      labels,
      vector_type,
      dirichlet_tags,
      dirichlet_masks,
      trian)
    return V
  end

  cell_fe = CellFE(model,cell_reffe,conf)
  _vector_type = _get_vector_type(vector_type,cell_fe,trian)
  if conformity in (L2Conformity(),:L2) && dirichlet_tags == Int[]
    F = _DiscontinuousFESpace(_vector_type,trian,cell_fe)
    V = _add_constraint(F,cell_fe.max_order,constraint)
  else
    V = FESpace(model,cell_fe;
      trian=trian,
      labels=labels,
      dirichlet_tags=dirichlet_tags,
      dirichlet_masks=dirichlet_masks,
      constraint=constraint,
      vector_type=_vector_type)
  end
  return V
end

function FESpace(model::DiscreteModel, reffe::Tuple{<:ReferenceFEName,Any,Any}; kwargs...)
  basis, reffe_args,reffe_kwargs = reffe
  cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
  FESpace(model,cell_reffe;kwargs...)
end

function FESpace(model::DiscreteModel, reffe::ReferenceFE; kwargs...)
  cell_reffe = Fill(reffe,num_cells(model))
  FESpace(model,cell_reffe;kwargs...)
end

"""
"""
function TestFESpace(args...;kwargs...)
  FESpace(args...;kwargs...)
end


text:
using Images

"""
Creates a raster image for the recorded spikes.

xwith: binning in time units
ywidth: binning of neuron IDs
"""
function raster_spikes(spikes::RecordedSpikes, xwidth, ywidth=1;
                       t0=first(spikes.steps),
                       t1=last(spikes.steps))
    s0 = searchsortedfirst(spikes.ts, t0)
    s1 = searchsortedlast(spikes.ts, t1)
    xsize = ceil(Int, (t1-t0) / xwidth)
    ysize = ceil(Int, length(spikes.instance) / ywidth)
    canvas = fill(0.0, ysize, xsize)
    for s in s0:s1
        t = spikes.ts[s]
        i = spikes.id[s]
        x = 1 + floor(Int, (t-t0) / xwidth)
        y = 1 + floor(Int, (i-1) / ywidth)
        if x > xsize || y > ysize
            continue
        end
        canvas[y,x] += 1.0
    end
    canvas
end

function density_plot!(canvas, xs, ys, xrange=extrema(xs), yrange=extrema(ys))
    h, w = size(canvas)
    xmin, xmax = extrema(xrange)
    ymin, ymax = extrema(yrange)
    xwidth = w/(xmax-xmin)
    ywidth = h/(ymax-ymin)
    for (x,y) in zip(xs,ys)
        i = 1 + floor(Int, (x-xmin) * xwidth)
        i > w && continue
        i < 1 && continue
        j = 1 + floor(Int, (y-ymin) * ywidth)
        j > h && continue
        j < 1 && continue
        canvas[j,i] += 1.0
    end
    canvas
end
export density_plot!

function normed(x::AbstractArray{T,3}, each=true) where T
    xnorm = if each
        maximum(abs.(x), (2,3))
    else
        maximum(abs.(x))
    end
    y = x./2xnorm
    y += 0.5
    y
end
export normed

function imgrid(rfs::AbstractArray{T,3}, nrows=1;
                pad=1, padval=0.0f0) where T
    ncols = ceil(Int, size(rfs,1)/nrows)
    h = size(rfs,2)+pad
    w = size(rfs,3)+pad
    canvas = zeros(T, (nrows*h+pad, ncols*w+pad))
    fill!(canvas, padval)
    for (i, loc) in enumerate(CartesianRange((nrows, ncols)))
        i > size(rfs,1) && break
        y, x = loc.I
        canvas[(y-1)*h+1+pad:y*h, (x-1)*w+1+pad:x*w] .= rfs[i,:,:]
    end
    Gray.(canvas)
end
export imgrid


text:
function is_jll_name(name::AbstractString)::Bool
    return endswith(name, "_jll")
end

function _get_all_dependencies_nonrecursive(working_directory::AbstractString,
                                            pkg,
                                            version)
    all_dependencies = String[]
    deps = Pkg.TOML.parsefile(joinpath(working_directory, uppercase(pkg[1:1]), pkg, "Deps.toml"))
    for version_range in keys(deps)
        if version in Pkg.Types.VersionRange(version_range)
            for name in keys(deps[version_range])
                push!(all_dependencies, name)
            end
        end
    end
    unique!(all_dependencies)
    return all_dependencies
end

const guideline_allowed_jll_nonrecursive_dependencies =
    Guideline("If this is a JLL package, only deps are Pkg, Libdl, and other JLL packages",
              data -> meets_allowed_jll_nonrecursive_dependencies(data.registry_head,
                                                                  data.pkg,
                                                                  data.version))

function meets_allowed_jll_nonrecursive_dependencies(working_directory::AbstractString,
                                                     pkg,
                                                     version)
    # If you are a JLL package, you are only allowed to have five kinds of dependencies:
    # 1. Pkg
    # 2. Libdl
    # 3. Artifacts
    # 4. JLLWrappers
    # 5. LazyArtifacts
    # 6. other JLL packages
    all_dependencies = _get_all_dependencies_nonrecursive(working_directory,
                                                          pkg,
                                                          version)
    allowed_dependencies = ("Pkg", "Libdl", "Artifacts", "JLLWrappers", "LazyArtifacts")
    for dep in all_dependencies
        if dep ∉ allowed_dependencies && !is_jll_name(dep)
            return false, "JLL packages are only allowed to depend on $(join(allowed_dependencies, ", ")) and other JLL packages"
        end
    end
    return true, ""
end


text:
using CoordinateTransformations: IdentityTransformation

@testset "Fractional coordinates to Cartesian coordinates" begin
    lattice = Lattice([
        1/2 0 0
        0 1/2 0
        0 0 1
    ])
    @test CartesianFromFractional(lattice)([2, 3, 1]) == [1, 3 / 2, 1]
    # Compared with "SeeK-path" example oC1 (SiTi)
    @testset "Base centered orthorhombic" begin
        a, b, c = 3.67939899, 4.15357986, 3.8236792350
        lattice = Lattice([a, -b, 0] / 2, [a, b, 0] / 2, [0, 0, c])
        f2c = CartesianFromFractional(lattice)
        @test f2c([1 / 2, 1 / 2, 1 / 2]) == [1.8396994950, 0, 1.9118396175]  # Si
        @test f2c([0, 0, 0]) == [0, 0, 0]  # Ti
        c2f = inv(f2c)
        @test c2f([1.8396994950, 0, 1.9118396175]) ≈ [1 / 2, 1 / 2, 1 / 2]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
end

# Example from http://ww1.iucr.org/iucr-top/comm/cteach/pamphlets/22/node34.html#SECTION00073300000000000000
@testset "Refer ZrSiO₄ primitive cell with origin at 2/m" begin
    𝐚ₚ, 𝐛ₚ, 𝐜ₚ = [1, 0, 0], [0, 1, 0], [1, 1, 1] // 2
    P = PrimitiveFromStandardized(hcat(𝐚ₚ, 𝐛ₚ, 𝐜ₚ))
    @test P([0, 1 // 4, 7 // 8]) == [-7, -5, 14] // 8  # The new coordinates of the first Zr atom
    @test P([0, 1 // 4, 3 // 8]) == [-3, -1, 6] // 8  # The new coordinates of the first Si atom
    @test P([0, 0.45, 0.215]) ≈ [-0.215, 0.235, 0.43]  # The new coordinates of the first O atom
    P⁻¹ = inv(P)
    @test P⁻¹([-7, -5, 14] // 8) == [0, 1 // 4, 7 // 8]
    @test P⁻¹([-3, -1, 6] // 8) == [0, 1 // 4, 3 // 8]
    @test P⁻¹([-0.215, 0.235, 0.43]) ≈ [0, 0.45, 0.215]
    @test P ∘ P⁻¹ == P⁻¹ ∘ P == IdentityTransformation()
end

@testset "Special points in the Brillouin zone" begin  # Reciprocal space transformation
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/orthorhombic.php
    @testset "Simple orthorhombic Brillouin zone" begin
        a, b, c = 2, 3, 5
        reci_lattice = inv(Lattice(a, b, c, 90, 90, 90))
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([1 / 2, 0, 0]) == [1 / a, 0, 0] / 2
        @test f2c([0, 1 / 2, 0]) == [0, 1 / b, 0] / 2
        @test f2c([0, 0, 1 / 2]) == [0, 0, 1 / c] / 2
        @test f2c([0, 1 / 2, 1 / 2]) == [0, 1 / b, 1 / c] / 2
        @test f2c([1 / 2, 1 / 2, 0]) == [1 / a, 1 / b, 0] / 2
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([1 / a, 0, 0] / 2) == [1 / 2, 0, 0]
        @test c2f([0, 1 / b, 0] / 2) == [0, 1 / 2, 0]
        @test c2f([0, 0, 1 / c] / 2) == [0, 0, 1 / 2]
        @test c2f([0, 1 / b, 1 / c] / 2) == [0, 1 / 2, 1 / 2]
        @test c2f([1 / a, 1 / b, 0] / 2) == [1 / 2, 1 / 2, 0]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # Compared with "SeeK-path" example oC1 (SiTi)
    # The results on http://lampx.tugraz.at/~hadley/ss1/bzones/orthorhombic_bc.php are wrong by a factor of 2.
    @testset "Base centered orthorhombic Brillouin zone" begin
        a, b, c = 3.67939899, 4.15357986, 3.8236792350
        reci_lattice = inv(Lattice([a, -b, 0] / 2, [a, b, 0] / 2, [0, 0, c]))
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([-1 / 2, 1 / 2, 0]) ≈ [0, 1 / b, 0] ≈ [0, 1.5127156619, 0] / 2pi  # Y'
        @test f2c([0, 0, 1 / 2]) ≈ [0, 0, 1 / 2 / c] ≈ [0, 0, 0.8216151148] / 2pi  # Z
        @test f2c([-1 / 2, 1 / 2, 1 / 2]) ≈
              [0, 1 / b, 1 / 2 / c] ≈
              [0, 1.5127156619, 0.8216151148] / 2pi  # T'
        @test f2c([0, 1 / 2, 0]) ≈
              [1 / 2 / a, 1 / 2 / b, 0] ≈
              [0.8538331021, 0.7563578310, 0] / 2pi  # S
        @test f2c([0, 1 / 2, 1 / 2]) ≈
              [1 / a, 1 / b, 1 / c] / 2 ≈
              [0.8538331021, 0.7563578310, 0.8216151148] / 2pi # R
        @test f2c([-0.4461772527, 0.5538227473, 1 / 2]) ≈
              [0.1838225731, 1.5127156619, 0.8216151148] / 2pi  # E0
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 1 / b, 0]) ≈ [-1 / 2, 1 / 2, 0]
        @test c2f([0, 0, 1 / 2 / c]) ≈ [0, 0, 1 / 2]
        @test c2f([0, 1 / b, 1 / 2 / c]) ≈ [-1 / 2, 1 / 2, 1 / 2]
        @test c2f([1 / 2 / a, 1 / 2 / b, 0]) ≈ [0, 1 / 2, 0]
        @test c2f([1 / a, 1 / b, 1 / c] / 2) ≈ [0, 1 / 2, 1 / 2]
        @test c2f([0.1838225731, 1.5127156619, 0.8216151148] / 2pi) ≈
              [-0.4461772527, 0.5538227473, 1 / 2]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/tetbc.php
    @testset "Body centered tetragonal Brillouin zone" begin
        a, c = 4, 6
        reci_lattice = inv(Lattice([a, a, -c] / 2, [a, -a, c] / 2, [-a, a, c] / 2))
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([0, 0, 0]) == [0, 0, 0]
        @test f2c([1 / 2, 0, 0]) == [1 / a, 1 / a, 0] / 2
        @test f2c([1 / 2, 1 / 2, -1 / 2]) == [1 / a, 0, 0]
        @test f2c([0, 1 / 2, 0]) == [1 / a, 0, 1 / c] / 2
        @test f2c([1, 1, 1] / 4) == [1 / a, 1 / a, 1 / c] / 2
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 0, 0]) == [0, 0, 0]
        @test c2f([1 / a, 1 / a, 0] / 2) == [1 / 2, 0, 0]
        @test c2f([1 / a, 0, 0]) == [1 / 2, 1 / 2, -1 / 2]
        @test c2f([1 / a, 0, 1 / c] / 2) ≈ [0, 1 / 2, 0]
        @test c2f([1 / a, 1 / a, 1 / c] / 2) == [1, 1, 1] / 4
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/hexagonal.php
    @testset "Simple hexagonal Brillouin zone" begin
        a, c = 2, 3.2
        reci_lattice = inv(Lattice([a, 0, 0], [a / 2, sqrt(3) / 2 * a, 0], [0, 0, c]))
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([0, 0, 0]) == [0, 0, 0]  # Γ
        @test f2c([1 / 2, 0, 0]) ≈ [1 / a, -1 / sqrt(3) / a, 0] / 2  # M
        @test f2c([0, 0, 1 / 2]) ≈ [0, 0, 1 / c] / 2  # A
        @test f2c([2 / 3, 1 / 3, 0]) ≈ [2 / 3 / a, 0, 0]  # K
        @test f2c([2 / 3, 1 / 3, 1 / 2]) ≈ [2 / 3 / a, 0, 1 / 2 / c]  # H
        @test f2c([1 / 2, 0, 1 / 2]) ≈ [1 / a, -1 / sqrt(3) / a, 1 / c] / 2  # L
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 0, 0]) == [0, 0, 0]
        @test c2f([1 / a, -1 / sqrt(3) / a, 0] / 2) ≈ [1 / 2, 0, 0]
        @test c2f([0, 0, 1 / c] / 2) == [0, 0, 1 / 2]
        @test c2f([2 / 3 / a, 0, 0]) == [2 / 3, 1 / 3, 0]
        @test c2f([2 / 3 / a, 0, 1 / 2 / c]) == [2 / 3, 1 / 3, 1 / 2]
        @test c2f([1 / a, -1 / sqrt(3) / a, 1 / c] / 2) ≈ [1 / 2, 0, 1 / 2]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/sc.php
    @testset "Simple cubic Brillouin zone" begin
        a = 4
        reci_lattice = inv(Lattice(a, a, a, 90, 90, 90))
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([0, 0, 0]) == [0, 0, 0]
        @test f2c([1 / 2, 1 / 2, 0]) == [1 / a, 1 / a, 0] / 2
        @test f2c([1 / 2, 1 / 2, 1 / 2]) == [1 / a, 1 / a, 1 / a] / 2
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 0, 0]) == [0, 0, 0]
        @test c2f([1 / a, 1 / a, 0] / 2) == [1 / 2, 1 / 2, 0]
        @test c2f([1 / a, 1 / a, 1 / a] / 2) == [1 / 2, 1 / 2, 1 / 2]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/fcc.php
    @testset "Face centered cubic Brillouin zone" begin
        a = 4
        reci_lattice = inv(Lattice([
            1 1 0
            0 1 1
            1 0 1
        ]) * a / 2)
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([0, 0, 0]) == [0, 0, 0]
        @test f2c([0, 1 / 2, 1 / 2]) == [0, 1 / a, 0]
        @test f2c([1 / 2, 1 / 2, 1 / 2]) == [1 / a, 1 / a, 1 / a] / 2
        @test f2c([1 / 4, 3 / 4, 1 / 2]) == [1 / 2 / a, 1 / a, 0]
        @test f2c([1 / 4, 5 / 8, 5 / 8]) == [1 / 2 / a, 2 / a, 1 / 2 / a] / 2
        @test f2c([3 / 8, 3 / 4, 3 / 8]) == [3 / 2 / a, 3 / 2 / a, 0] / 2
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 0, 0]) == [0, 0, 0]
        @test f2c([0, 1 / 2, 1 / 2]) == [0, 1 / a, 0]
        @test f2c([1 / 2, 1 / 2, 1 / 2]) == [1 / a, 1 / a, 1 / a] / 2
        @test c2f([1 / 2 / a, 1 / a, 0]) == [1 / 4, 3 / 4, 1 / 2]
        @test c2f([1 / 2 / a, 2 / a, 1 / 2 / a] / 2) == [1 / 4, 5 / 8, 5 / 8]
        @test c2f([3 / 2 / a, 3 / 2 / a, 0] / 2) == [3 / 8, 3 / 4, 3 / 8]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
    # See http://lampx.tugraz.at/~hadley/ss1/bzones/bcc.php
    @testset "Body centered cubic Brillouin zone" begin
        a = 4
        reci_lattice = inv(Lattice([
            1 -1 1
            1 1 -1
            -1 1 1
        ]) * a / 2)
        f2c = CartesianFromFractional(reci_lattice)
        @test f2c([0, 0, 0]) == [0, 0, 0]
        @test f2c([-1 / 2, 1 / 2, 1 / 2]) == [0, 0, 1 / a]
        @test f2c([1 / 4, 1 / 4, 1 / 4]) == [1 / a, 1 / a, 1 / a] / 2
        @test f2c([0, 1 / 2, 0]) == [0, 1 / a, 1 / a] / 2
        c2f = FractionalFromCartesian(reci_lattice)
        @test c2f([0, 0, 0]) == [0, 0, 0]
        @test c2f([0, 0, 1 / a]) == [-1 / 2, 1 / 2, 1 / 2]
        @test c2f([1 / a, 1 / a, 1 / a] / 2) == [1 / 4, 1 / 4, 1 / 4]
        @test c2f([0, 1 / a, 1 / a] / 2) == [0, 1 / 2, 0]
        @test c2f ∘ f2c == f2c ∘ c2f == IdentityTransformation()
    end
end

@testset "Test coordinate transformation with units" begin
    alat = 8.8097
    reci_lattice = inv(
        Lattice(
            [
                -2.739741805 2.739741805 2.739741805
                2.739741805 -2.739741805 2.739741805
                2.592418910 2.592418910 -2.592418910
            ] / alat * u"angstrom" .|> u"bohr",
        ),
    )
    # Compared with Quantum ESPRESSO results
    qe_crystal =
        [
            0.0 0.0 0.0
            -0.0 -0.0 0.016129
            0.0 -0.0 0.0322581
            -0.0 0.0 0.0483871
            0.0 -0.0 0.0645161
            -0.0 0.0 0.0806452
            -0.0 -0.0 0.0967742
            -0.0 0.0 0.1129032
            0.0 -0.0 0.1290323
            -0.0 0.0 0.1451613
            -0.0 0.0 0.1612903
            -0.0 -0.0 0.1774194
            0.0 -0.0 0.1935484
            -0.0 -0.0 0.2096774
            -0.0 0.0 0.2258065
            -0.0 0.0 0.2419355
            -0.0 -0.0 0.2580645
            -0.0 -0.0 0.2741935
            -0.0 0.0 0.2903226
            0.0 -0.0 0.3064516
            -0.0 0.0 0.3225806
            0.0 -0.0 0.3387097
            -0.0 -0.0 0.3548387
            -0.0 -0.0 0.3709677
            0.0 -0.0 0.3870968
            -0.0 0.0 0.4032258
            0.0 -0.0 0.4193548
            -0.0 -0.0 0.4354839
            -0.0 0.0 0.4516129
            -0.0 -0.0 0.4677419
            -0.0 -0.0 0.483871
            0.0 -0.0 0.5
            -0.016129 0.016129 0.5
            -0.0322581 0.0322581 0.5
            -0.0483871 0.0483871 0.5
            -0.0645161 0.0645161 0.5
            -0.0806452 0.0806452 0.5
            -0.0967742 0.0967742 0.5
            -0.1129032 0.1129032 0.5
            -0.1290323 0.1290323 0.5
            -0.1451613 0.1451613 0.5
            -0.1612903 0.1612903 0.5
            -0.1774194 0.1774194 0.5
            -0.1935484 0.1935484 0.5
            -0.2096774 0.2096774 0.5
            -0.2258065 0.2258065 0.5
            -0.2419355 0.2419355 0.5
            -0.2580645 0.2580645 0.5
            -0.2741935 0.2741935 0.5
            -0.2903226 0.2903226 0.5
            -0.3064516 0.3064516 0.5
            -0.3225806 0.3225806 0.5
            -0.3387097 0.3387097 0.5
            -0.3548387 0.3548387 0.5
            -0.3709677 0.3709677 0.5
            -0.3870968 0.3870968 0.5
            -0.4032258 0.4032258 0.5
            -0.4193548 0.4193548 0.5
            -0.4354839 0.4354839 0.5
            -0.4516129 0.4516129 0.5
            -0.4677419 0.4677419 0.5
            -0.483871 0.483871 0.5
            -0.5 0.5 0.5
            -0.4886364 0.4886364 0.4886364
            -0.4772727 0.4772727 0.4772727
            -0.4659091 0.4659091 0.4659091
            -0.4545455 0.4545455 0.4545455
            -0.4431818 0.4431818 0.4431818
            -0.4318182 0.4318182 0.4318182
            -0.4204545 0.4204545 0.4204545
            -0.4090909 0.4090909 0.4090909
            -0.3977273 0.3977273 0.3977273
            -0.3863636 0.3863636 0.3863636
            -0.375 0.375 0.375
            -0.3636364 0.3636364 0.3636364
            -0.3522727 0.3522727 0.3522727
            -0.3409091 0.3409091 0.3409091
            -0.3295455 0.3295455 0.3295455
            -0.3181818 0.3181818 0.3181818
            -0.3068182 0.3068182 0.3068182
            -0.2954545 0.2954545 0.2954545
            -0.2840909 0.2840909 0.2840909
            -0.2727273 0.2727273 0.2727273
            -0.2613636 0.2613636 0.2613636
            -0.25 0.25 0.25
            -0.2386364 0.2386364 0.2386364
            -0.2272727 0.2272727 0.2272727
            -0.2159091 0.2159091 0.2159091
            -0.2045455 0.2045455 0.2045455
            -0.1931818 0.1931818 0.1931818
            -0.1818182 0.1818182 0.1818182
            -0.1704545 0.1704545 0.1704545
            -0.1590909 0.1590909 0.1590909
            -0.1477273 0.1477273 0.1477273
            -0.1363636 0.1363636 0.1363636
            -0.125 0.125 0.125
            -0.1136364 0.1136364 0.1136364
            -0.1022727 0.1022727 0.1022727
            -0.0909091 0.0909091 0.0909091
            -0.0795455 0.0795455 0.0795455
            -0.0681818 0.0681818 0.0681818
            -0.0568182 0.0568182 0.0568182
            -0.0454545 0.0454545 0.0454545
            -0.0340909 0.0340909 0.0340909
            -0.0227273 0.0227273 0.0227273
            -0.0113636 0.0113636 0.0113636
            0.0 0.0 0.0
            0.010769 0.010769 -0.010769
            0.021538 0.021538 -0.021538
            0.032307 0.032307 -0.032307
            0.0430761 0.0430761 -0.0430761
            0.0538451 0.0538451 -0.0538451
            0.0646141 0.0646141 -0.0646141
            0.0753831 0.0753831 -0.0753831
            0.0861521 0.0861521 -0.0861521
            0.0969211 0.0969211 -0.0969211
            0.1076901 0.1076901 -0.1076901
            0.1184591 0.1184591 -0.1184591
            0.1292282 0.1292282 -0.1292282
            0.1399972 0.1399972 -0.1399972
            0.1507662 0.1507662 -0.1507662
            0.1615352 0.1615352 -0.1615352
            0.1723042 0.1723042 -0.1723042
            0.1830732 0.1830732 -0.1830732
            0.1938422 0.1938422 -0.1938422
            0.2046113 0.2046113 -0.2046113
            0.2153803 0.2153803 -0.2153803
            0.2261493 0.2261493 -0.2261493
            0.2369183 0.2369183 -0.2369183
            0.2476873 0.2476873 -0.2476873
            0.2584563 0.2584563 -0.2584563
            0.2692253 0.2692253 -0.2692253
            0.2799944 0.2799944 -0.2799944
            0.2907634 0.2907634 -0.2907634
            0.3015324 0.3015324 -0.3015324
            0.3123014 0.3123014 -0.3123014
            0.3230704 0.3230704 -0.3230704
            0.3338394 0.3338394 -0.3338394
            0.3446084 0.3446084 -0.3446084
            0.3553774 0.3553774 -0.3553774
            0.3661465 0.3661465 -0.3661465
            0.3769155 0.3769155 -0.3769155
            0.3876845 0.3876845 -0.3876845
            0.3984535 0.3984535 -0.3984535
            0.4092225 0.4092225 -0.4092225
            0.4199915 0.4199915 -0.4199915
            0.4307605 0.4307605 -0.4307605
            0.4415296 0.4415296 -0.4415296
            0.4522986 0.4522986 -0.4522986
            0.4630676 0.4630676 -0.4630676
            0.4738366 0.4738366 -0.4738366
            -0.4738366 0.5261634 0.4738366
            -0.5 0.5 0.5
            0.0 -0.0 0.5
            0.0108696 0.0108696 0.4891304
            0.0217391 0.0217391 0.4782609
            0.0326087 0.0326087 0.4673913
            0.0434783 0.0434783 0.4565217
            0.0543478 0.0543478 0.4456522
            0.0652174 0.0652174 0.4347826
            0.076087 0.076087 0.423913
            0.0869565 0.0869565 0.4130435
            0.0978261 0.0978261 0.4021739
            0.1086957 0.1086957 0.3913043
            0.1195652 0.1195652 0.3804348
            0.1304348 0.1304348 0.3695652
            0.1413043 0.1413043 0.3586957
            0.1521739 0.1521739 0.3478261
            0.1630435 0.1630435 0.3369565
            0.173913 0.173913 0.326087
            0.1847826 0.1847826 0.3152174
            0.1956522 0.1956522 0.3043478
            0.2065217 0.2065217 0.2934783
            0.2173913 0.2173913 0.2826087
            0.2282609 0.2282609 0.2717391
            0.2391304 0.2391304 0.2608696
            0.25 0.25 0.25
            0.2380952 0.2619048 0.2380952
            0.2261905 0.2738095 0.2261905
            0.2142857 0.2857143 0.2142857
            0.202381 0.297619 0.202381
            0.1904762 0.3095238 0.1904762
            0.1785714 0.3214286 0.1785714
            0.1666667 0.3333333 0.1666667
            0.1547619 0.3452381 0.1547619
            0.1428571 0.3571429 0.1428571
            0.1309524 0.3690476 0.1309524
            0.1190476 0.3809524 0.1190476
            0.1071429 0.3928571 0.1071429
            0.0952381 0.4047619 0.0952381
            0.0833333 0.4166667 0.0833333
            0.0714286 0.4285714 0.0714286
            0.0595238 0.4404762 0.0595238
            0.047619 0.452381 0.047619
            0.0357143 0.4642857 0.0357143
            0.0238095 0.4761905 0.0238095
            0.0119048 0.4880952 0.0119048
            -0.0 0.5 0.0
            -0.0 0.484375 0.0
            -0.0 0.46875 0.0
            0.0 0.453125 -0.0
            -0.0 0.4375 0.0
            -0.0 0.421875 0.0
            -0.0 0.40625 0.0
            0.0 0.390625 -0.0
            -0.0 0.375 0.0
            -0.0 0.359375 0.0
            -0.0 0.34375 0.0
            -0.0 0.328125 0.0
            0.0 0.3125 0.0
            -0.0 0.296875 0.0
            -0.0 0.28125 0.0
            0.0 0.265625 0.0
            -0.0 0.25 0.0
            -0.0 0.234375 0.0
            -0.0 0.21875 0.0
            -0.0 0.203125 0.0
            -0.0 0.1875 0.0
            -0.0 0.171875 0.0
            0.0 0.15625 0.0
            -0.0 0.140625 0.0
            -0.0 0.125 0.0
            -0.0 0.109375 0.0
            -0.0 0.09375 0.0
            0.0 0.078125 0.0
            -0.0 0.0625 0.0
            -0.0 0.046875 0.0
            -0.0 0.03125 0.0
            -0.0 0.015625 0.0
            0.0 0.0 0.0
            0.0 0.0 -0.016129
            -0.0 0.0 -0.0322581
            0.0 -0.0 -0.0483871
            -0.0 0.0 -0.0645161
            0.0 -0.0 -0.0806452
            0.0 0.0 -0.0967742
            0.0 -0.0 -0.1129032
            -0.0 0.0 -0.1290323
            0.0 -0.0 -0.1451613
            0.0 -0.0 -0.1612903
            0.0 0.0 -0.1774194
            -0.0 0.0 -0.1935484
            0.0 0.0 -0.2096774
            0.0 -0.0 -0.2258065
            0.0 -0.0 -0.2419355
            0.0 0.0 -0.2580645
            0.0 0.0 -0.2741935
            0.0 -0.0 -0.2903226
            -0.0 0.0 -0.3064516
            0.0 -0.0 -0.3225806
            -0.0 0.0 -0.3387097
            0.0 0.0 -0.3548387
            0.0 0.0 -0.3709677
            -0.0 0.0 -0.3870968
            0.0 -0.0 -0.4032258
            -0.0 0.0 -0.4193548
            0.0 0.0 -0.4354839
            0.0 -0.0 -0.4516129
            0.0 0.0 -0.4677419
            0.0 0.0 -0.483871
            -0.0 0.0 -0.5
            0.016129 -0.016129 -0.5
            0.0322581 -0.0322581 -0.5
            0.0483871 -0.0483871 -0.5
            0.0645161 -0.0645161 -0.5
            0.0806452 -0.0806452 -0.5
            0.0967742 -0.0967742 -0.5
            0.1129032 -0.1129032 -0.5
            0.1290323 -0.1290323 -0.5
            0.1451613 -0.1451613 -0.5
            0.1612903 -0.1612903 -0.5
            0.1774194 -0.1774194 -0.5
            0.1935484 -0.1935484 -0.5
            0.2096774 -0.2096774 -0.5
            0.2258065 -0.2258065 -0.5
            0.2419355 -0.2419355 -0.5
            0.2580645 -0.2580645 -0.5
            0.2741935 -0.2741935 -0.5
            0.2903226 -0.2903226 -0.5
            0.3064516 -0.3064516 -0.5
            0.3225806 -0.3225806 -0.5
            0.3387097 -0.3387097 -0.5
            0.3548387 -0.3548387 -0.5
            0.3709677 -0.3709677 -0.5
            0.3870968 -0.3870968 -0.5
            0.4032258 -0.4032258 -0.5
            0.4193548 -0.4193548 -0.5
            0.4354839 -0.4354839 -0.5
            0.4516129 -0.4516129 -0.5
            0.4677419 -0.4677419 -0.5
            0.483871 -0.483871 -0.5
            0.5 -0.5 -0.5
            0.4886364 -0.4886364 -0.4886364
            0.4772727 -0.4772727 -0.4772727
            0.4659091 -0.4659091 -0.4659091
            0.4545455 -0.4545455 -0.4545455
            0.4431818 -0.4431818 -0.4431818
            0.4318182 -0.4318182 -0.4318182
            0.4204545 -0.4204545 -0.4204545
            0.4090909 -0.4090909 -0.4090909
            0.3977273 -0.3977273 -0.3977273
            0.3863636 -0.3863636 -0.3863636
            0.375 -0.375 -0.375
            0.3636364 -0.3636364 -0.3636364
            0.3522727 -0.3522727 -0.3522727
            0.3409091 -0.3409091 -0.3409091
            0.3295455 -0.3295455 -0.3295455
            0.3181818 -0.3181818 -0.3181818
            0.3068182 -0.3068182 -0.3068182
            0.2954545 -0.2954545 -0.2954545
            0.2840909 -0.2840909 -0.2840909
            0.2727273 -0.2727273 -0.2727273
            0.2613636 -0.2613636 -0.2613636
            0.25 -0.25 -0.25
            0.2386364 -0.2386364 -0.2386364
            0.2272727 -0.2272727 -0.2272727
            0.2159091 -0.2159091 -0.2159091
            0.2045455 -0.2045455 -0.2045455
            0.1931818 -0.1931818 -0.1931818
            0.1818182 -0.1818182 -0.1818182
            0.1704545 -0.1704545 -0.1704545
            0.1590909 -0.1590909 -0.1590909
            0.1477273 -0.1477273 -0.1477273
            0.1363636 -0.1363636 -0.1363636
            0.125 -0.125 -0.125
            0.1136364 -0.1136364 -0.1136364
            0.1022727 -0.1022727 -0.1022727
            0.0909091 -0.0909091 -0.0909091
            0.0795455 -0.0795455 -0.0795455
            0.0681818 -0.0681818 -0.0681818
            0.0568182 -0.0568182 -0.0568182
            0.0454545 -0.0454545 -0.0454545
            0.0340909 -0.0340909 -0.0340909
            0.0227273 -0.0227273 -0.0227273
            0.0113636 -0.0113636 -0.0113636
            0.0 0.0 0.0
            -0.010769 -0.010769 0.010769
            -0.021538 -0.021538 0.021538
            -0.032307 -0.032307 0.032307
            -0.0430761 -0.0430761 0.0430761
            -0.0538451 -0.0538451 0.0538451
            -0.0646141 -0.0646141 0.0646141
            -0.0753831 -0.0753831 0.0753831
            -0.0861521 -0.0861521 0.0861521
            -0.0969211 -0.0969211 0.0969211
            -0.1076901 -0.1076901 0.1076901
            -0.1184591 -0.1184591 0.1184591
            -0.1292282 -0.1292282 0.1292282
            -0.1399972 -0.1399972 0.1399972
            -0.1507662 -0.1507662 0.1507662
            -0.1615352 -0.1615352 0.1615352
            -0.1723042 -0.1723042 0.1723042
            -0.1830732 -0.1830732 0.1830732
            -0.1938422 -0.1938422 0.1938422
            -0.2046113 -0.2046113 0.2046113
            -0.2153803 -0.2153803 0.2153803
            -0.2261493 -0.2261493 0.2261493
            -0.2369183 -0.2369183 0.2369183
            -0.2476873 -0.2476873 0.2476873
            -0.2584563 -0.2584563 0.2584563
            -0.2692253 -0.2692253 0.2692253
            -0.2799944 -0.2799944 0.2799944
            -0.2907634 -0.2907634 0.2907634
            -0.3015324 -0.3015324 0.3015324
            -0.3123014 -0.3123014 0.3123014
            -0.3230704 -0.3230704 0.3230704
            -0.3338394 -0.3338394 0.3338394
            -0.3446084 -0.3446084 0.3446084
            -0.3553774 -0.3553774 0.3553774
            -0.3661465 -0.3661465 0.3661465
            -0.3769155 -0.3769155 0.3769155
            -0.3876845 -0.3876845 0.3876845
            -0.3984535 -0.3984535 0.3984535
            -0.4092225 -0.4092225 0.4092225
            -0.4199915 -0.4199915 0.4199915
            -0.4307605 -0.4307605 0.4307605
            -0.4415296 -0.4415296 0.4415296
            -0.4522986 -0.4522986 0.4522986
            -0.4630676 -0.4630676 0.4630676
            -0.4738366 -0.4738366 0.4738366
            0.4738366 -0.5261634 -0.4738366
            0.5 -0.5 -0.5
            -0.0 0.0 -0.5
            -0.0108696 -0.0108696 -0.4891304
            -0.0217391 -0.0217391 -0.4782609
            -0.0326087 -0.0326087 -0.4673913
            -0.0434783 -0.0434783 -0.4565217
            -0.0543478 -0.0543478 -0.4456522
            -0.0652174 -0.0652174 -0.4347826
            -0.076087 -0.076087 -0.423913
            -0.0869565 -0.0869565 -0.4130435
            -0.0978261 -0.0978261 -0.4021739
            -0.1086957 -0.1086957 -0.3913043
            -0.1195652 -0.1195652 -0.3804348
            -0.1304348 -0.1304348 -0.3695652
            -0.1413043 -0.1413043 -0.3586957
            -0.1521739 -0.1521739 -0.3478261
            -0.1630435 -0.1630435 -0.3369565
            -0.173913 -0.173913 -0.326087
            -0.1847826 -0.1847826 -0.3152174
            -0.1956522 -0.1956522 -0.3043478
            -0.2065217 -0.2065217 -0.2934783
            -0.2173913 -0.2173913 -0.2826087
            -0.2282609 -0.2282609 -0.2717391
            -0.2391304 -0.2391304 -0.2608696
            -0.25 -0.25 -0.25
            -0.2380952 -0.2619048 -0.2380952
            -0.2261905 -0.2738095 -0.2261905
            -0.2142857 -0.2857143 -0.2142857
            -0.202381 -0.297619 -0.202381
            -0.1904762 -0.3095238 -0.1904762
            -0.1785714 -0.3214286 -0.1785714
            -0.1666667 -0.3333333 -0.1666667
            -0.1547619 -0.3452381 -0.1547619
            -0.1428571 -0.3571429 -0.1428571
            -0.1309524 -0.3690476 -0.1309524
            -0.1190476 -0.3809524 -0.1190476
            -0.1071429 -0.3928571 -0.1071429
            -0.0952381 -0.4047619 -0.0952381
            -0.0833333 -0.4166667 -0.0833333
            -0.0714286 -0.4285714 -0.0714286
            -0.0595238 -0.4404762 -0.0595238
            -0.047619 -0.452381 -0.047619
            -0.0357143 -0.4642857 -0.0357143
            -0.0238095 -0.4761905 -0.0238095
            -0.0119048 -0.4880952 -0.0119048
            0.0 -0.5 -0.0
            0.0 -0.484375 -0.0
            0.0 -0.46875 -0.0
            -0.0 -0.453125 0.0
            0.0 -0.4375 -0.0
            0.0 -0.421875 -0.0
            0.0 -0.40625 -0.0
            -0.0 -0.390625 0.0
            0.0 -0.375 -0.0
            0.0 -0.359375 -0.0
            0.0 -0.34375 -0.0
            0.0 -0.328125 -0.0
            -0.0 -0.3125 -0.0
            0.0 -0.296875 -0.0
            0.0 -0.28125 -0.0
            -0.0 -0.265625 -0.0
            0.0 -0.25 -0.0
            0.0 -0.234375 -0.0
            0.0 -0.21875 -0.0
            0.0 -0.203125 -0.0
            0.0 -0.1875 -0.0
            0.0 -0.171875 -0.0
            -0.0 -0.15625 -0.0
            0.0 -0.140625 -0.0
            0.0 -0.125 -0.0
            0.0 -0.109375 -0.0
            0.0 -0.09375 -0.0
            -0.0 -0.078125 -0.0
            0.0 -0.0625 -0.0
            0.0 -0.046875 -0.0
            0.0 -0.03125 -0.0
            0.0 -0.015625 -0.0
            0.0 0.0 0.0
        ]'
    qe_cartesian =
        [
            0.0 0.0 0.0
            0.0137224 0.0137224 -0.0
            0.0274447 0.0274447 -0.0
            0.0411671 0.0411671 -0.0
            0.0548894 0.0548894 -0.0
            0.0686118 0.0686118 -0.0
            0.0823341 0.0823341 -0.0
            0.0960565 0.0960565 -0.0
            0.1097788 0.1097788 -0.0
            0.1235012 0.1235012 -0.0
            0.1372235 0.1372235 -0.0
            0.1509459 0.1509459 -0.0
            0.1646682 0.1646682 -0.0
            0.1783906 0.1783906 -0.0
            0.1921129 0.1921129 -0.0
            0.2058353 0.2058353 -0.0
            0.2195576 0.2195576 -0.0
            0.23328 0.23328 -0.0
            0.2470023 0.2470023 -0.0
            0.2607247 0.2607247 -0.0
            0.2744471 0.2744471 -0.0
            0.2881694 0.2881694 -0.0
            0.3018918 0.3018918 -0.0
            0.3156141 0.3156141 -0.0
            0.3293365 0.3293365 -0.0
            0.3430588 0.3430588 -0.0
            0.3567812 0.3567812 -0.0
            0.3705035 0.3705035 -0.0
            0.3842259 0.3842259 -0.0
            0.3979482 0.3979482 -0.0
            0.4116706 0.4116706 -0.0
            0.4253929 0.4253929 -0.0
            0.4391153 0.4116706 -0.0
            0.4528376 0.3979482 -0.0
            0.46656 0.3842259 -0.0
            0.4802823 0.3705035 -0.0
            0.4940047 0.3567812 -0.0
            0.5077271 0.3430588 -0.0
            0.5214494 0.3293365 -0.0
            0.5351718 0.3156141 -0.0
            0.5488941 0.3018918 -0.0
            0.5626165 0.2881694 -0.0
            0.5763388 0.2744471 -0.0
            0.5900612 0.2607247 -0.0
            0.6037835 0.2470023 -0.0
            0.6175059 0.23328 -0.0
            0.6312282 0.2195576 -0.0
            0.6449506 0.2058353 -0.0
            0.6586729 0.1921129 -0.0
            0.6723953 0.1783906 -0.0
            0.6861176 0.1646682 -0.0
            0.69984 0.1509459 0.0
            0.7135623 0.1372235 -0.0
            0.7272847 0.1235012 -0.0
            0.741007 0.1097788 0.0
            0.7547294 0.0960565 -0.0
            0.7684518 0.0823341 -0.0
            0.7821741 0.0686118 0.0
            0.7958965 0.0548894 -0.0
            0.8096188 0.0411671 -0.0
            0.8233412 0.0274447 0.0
            0.8370635 0.0137224 -0.0
            0.8507859 0.0 -0.0
            0.8314498 -0.0 -0.0
            0.8121138 0.0 0.0
            0.7927777 -0.0 -0.0
            0.7734417 0.0 -0.0
            0.7541057 -0.0 -0.0
            0.7347696 -0.0 -0.0
            0.7154336 -0.0 -0.0
            0.6960975 0.0 0.0
            0.6767615 0.0 -0.0
            0.6574254 -0.0 -0.0
            0.6380894 0.0 -0.0
            0.6187534 -0.0 0.0
            0.5994173 0.0 0.0
            0.5800813 0.0 -0.0
            0.5607452 -0.0 -0.0
            0.5414092 -0.0 -0.0
            0.5220731 -0.0 -0.0
            0.5027371 0.0 0.0
            0.4834011 0.0 0.0
            0.464065 -0.0 0.0
            0.444729 -0.0 -0.0
            0.4253929 0.0 -0.0
            0.4060569 0.0 0.0
            0.3867209 0.0 -0.0
            0.3673848 -0.0 -0.0
            0.3480488 -0.0 -0.0
            0.3287127 -0.0 -0.0
            0.3093767 -0.0 0.0
            0.2900406 -0.0 0.0
            0.2707046 -0.0 -0.0
            0.2513686 0.0 0.0
            0.2320325 0.0 -0.0
            0.2126965 0.0 -0.0
            0.1933604 -0.0 -0.0
            0.1740244 0.0 0.0
            0.1546883 -0.0 0.0
            0.1353523 -0.0 -0.0
            0.1160163 0.0 -0.0
            0.0966802 -0.0 -0.0
            0.0773442 -0.0 -0.0
            0.0580081 0.0 -0.0
            0.0386721 0.0 -0.0
            0.019336 -0.0 -0.0
            0.0 0.0 0.0
            -0.0 -0.0 0.0193656
            -0.0 -0.0 0.0387312
            -0.0 -0.0 0.0580968
            -0.0 -0.0 0.0774623
            -0.0 -0.0 0.0968279
            -0.0 -0.0 0.1161935
            0.0 0.0 0.1355591
            -0.0 -0.0 0.1549247
            0.0 0.0 0.1742903
            -0.0 -0.0 0.1936559
            0.0 0.0 0.2130214
            -0.0 -0.0 0.232387
            -0.0 -0.0 0.2517526
            0.0 0.0 0.2711182
            0.0 0.0 0.2904838
            0.0 0.0 0.3098494
            -0.0 -0.0 0.329215
            -0.0 -0.0 0.3485805
            0.0 0.0 0.3679461
            0.0 0.0 0.3873117
            0.0 0.0 0.4066773
            -0.0 -0.0 0.4260429
            -0.0 -0.0 0.4454085
            0.0 0.0 0.4647741
            0.0 0.0 0.4841397
            0.0 0.0 0.5035052
            -0.0 -0.0 0.5228708
            -0.0 -0.0 0.5422364
            0.0 0.0 0.561602
            -0.0 -0.0 0.5809676
            0.0 0.0 0.6003332
            0.0 0.0 0.6196988
            0.0 0.0 0.6390643
            -0.0 -0.0 0.6584299
            0.0 0.0 0.6777955
            -0.0 -0.0 0.6971611
            0.0 0.0 0.7165267
            0.0 0.0 0.7358923
            -0.0 -0.0 0.7552579
            0.0 0.0 0.7746234
            -0.0 -0.0 0.793989
            0.0 0.0 0.8133546
            -0.0 -0.0 0.8327202
            -0.0 -0.0 0.8520858
            0.8507859 0.0 0.0470488
            0.8507859 0.0 -0.0
            0.4253929 0.4253929 -0.0
            0.4253929 0.4253929 0.0195464
            0.4253929 0.4253929 0.0390928
            0.4253929 0.4253929 0.0586392
            0.4253929 0.4253929 0.0781856
            0.4253929 0.4253929 0.097732
            0.4253929 0.4253929 0.1172784
            0.4253929 0.4253929 0.1368248
            0.4253929 0.4253929 0.1563712
            0.4253929 0.4253929 0.1759176
            0.4253929 0.4253929 0.1954641
            0.4253929 0.4253929 0.2150105
            0.4253929 0.4253929 0.2345569
            0.4253929 0.4253929 0.2541033
            0.4253929 0.4253929 0.2736497
            0.4253929 0.4253929 0.2931961
            0.4253929 0.4253929 0.3127425
            0.4253929 0.4253929 0.3322889
            0.4253929 0.4253929 0.3518353
            0.4253929 0.4253929 0.3713817
            0.4253929 0.4253929 0.3909281
            0.4253929 0.4253929 0.4104745
            0.4253929 0.4253929 0.4300209
            0.4253929 0.4253929 0.4495673
            0.4253929 0.4051361 0.4495673
            0.4253929 0.3848793 0.4495673
            0.4253929 0.3646225 0.4495673
            0.4253929 0.3443657 0.4495673
            0.4253929 0.3241089 0.4495673
            0.4253929 0.3038521 0.4495673
            0.4253929 0.2835953 0.4495673
            0.4253929 0.2633385 0.4495673
            0.4253929 0.2430817 0.4495673
            0.4253929 0.2228249 0.4495673
            0.4253929 0.2025681 0.4495673
            0.4253929 0.1823113 0.4495673
            0.4253929 0.1620545 0.4495673
            0.4253929 0.1417976 0.4495673
            0.4253929 0.1215408 0.4495673
            0.4253929 0.101284 0.4495673
            0.4253929 0.0810272 0.4495673
            0.4253929 0.0607704 0.4495673
            0.4253929 0.0405136 0.4495673
            0.4253929 0.0202568 0.4495673
            0.4253929 0.0 0.4495673
            0.4120994 0.0 0.4355183
            0.3988059 0.0 0.4214694
            0.3855123 0.0 0.4074204
            0.3722188 0.0 0.3933714
            0.3589253 0.0 0.3793224
            0.3456318 0.0 0.3652734
            0.3323382 0.0 0.3512245
            0.3190447 0.0 0.3371755
            0.3057512 0.0 0.3231265
            0.2924576 0.0 0.3090775
            0.2791641 0.0 0.2950286
            0.2658706 0.0 0.2809796
            0.2525771 0.0 0.2669306
            0.2392835 0.0 0.2528816
            0.22599 0.0 0.2388326
            0.2126965 0.0 0.2247837
            0.1994029 0.0 0.2107347
            0.1861094 0.0 0.1966857
            0.1728159 0.0 0.1826367
            0.1595224 0.0 0.1685877
            0.1462288 0.0 0.1545388
            0.1329353 0.0 0.1404898
            0.1196418 0.0 0.1264408
            0.1063482 0.0 0.1123918
            0.0930547 0.0 0.0983429
            0.0797612 0.0 0.0842939
            0.0664676 0.0 0.0702449
            0.0531741 0.0 0.0561959
            0.0398806 0.0 0.0421469
            0.0265871 0.0 0.028098
            0.0132935 0.0 0.014049
            0.0 0.0 0.0
            -0.0137224 -0.0137224 0.0
            -0.0274447 -0.0274447 0.0
            -0.0411671 -0.0411671 0.0
            -0.0548894 -0.0548894 0.0
            -0.0686118 -0.0686118 0.0
            -0.0823341 -0.0823341 0.0
            -0.0960565 -0.0960565 0.0
            -0.1097788 -0.1097788 0.0
            -0.1235012 -0.1235012 0.0
            -0.1372235 -0.1372235 0.0
            -0.1509459 -0.1509459 0.0
            -0.1646682 -0.1646682 0.0
            -0.1783906 -0.1783906 0.0
            -0.1921129 -0.1921129 0.0
            -0.2058353 -0.2058353 0.0
            -0.2195576 -0.2195576 0.0
            -0.23328 -0.23328 0.0
            -0.2470023 -0.2470023 0.0
            -0.2607247 -0.2607247 0.0
            -0.2744471 -0.2744471 0.0
            -0.2881694 -0.2881694 0.0
            -0.3018918 -0.3018918 0.0
            -0.3156141 -0.3156141 0.0
            -0.3293365 -0.3293365 0.0
            -0.3430588 -0.3430588 0.0
            -0.3567812 -0.3567812 0.0
            -0.3705035 -0.3705035 0.0
            -0.3842259 -0.3842259 0.0
            -0.3979482 -0.3979482 0.0
            -0.4116706 -0.4116706 0.0
            -0.4253929 -0.4253929 0.0
            -0.4391153 -0.4116706 0.0
            -0.4528376 -0.3979482 0.0
            -0.46656 -0.3842259 0.0
            -0.4802823 -0.3705035 0.0
            -0.4940047 -0.3567812 0.0
            -0.5077271 -0.3430588 0.0
            -0.5214494 -0.3293365 0.0
            -0.5351718 -0.3156141 0.0
            -0.5488941 -0.3018918 0.0
            -0.5626165 -0.2881694 0.0
            -0.5763388 -0.2744471 0.0
            -0.5900612 -0.2607247 0.0
            -0.6037835 -0.2470023 0.0
            -0.6175059 -0.23328 0.0
            -0.6312282 -0.2195576 0.0
            -0.6449506 -0.2058353 0.0
            -0.6586729 -0.1921129 0.0
            -0.6723953 -0.1783906 0.0
            -0.6861176 -0.1646682 0.0
            -0.69984 -0.1509459 -0.0
            -0.7135623 -0.1372235 0.0
            -0.7272847 -0.1235012 0.0
            -0.741007 -0.1097788 -0.0
            -0.7547294 -0.0960565 0.0
            -0.7684518 -0.0823341 0.0
            -0.7821741 -0.0686118 -0.0
            -0.7958965 -0.0548894 0.0
            -0.8096188 -0.0411671 0.0
            -0.8233412 -0.0274447 -0.0
            -0.8370635 -0.0137224 0.0
            -0.8507859 0.0 0.0
            -0.8314498 0.0 0.0
            -0.8121138 -0.0 -0.0
            -0.7927777 0.0 0.0
            -0.7734417 -0.0 0.0
            -0.7541057 0.0 0.0
            -0.7347696 0.0 0.0
            -0.7154336 0.0 0.0
            -0.6960975 -0.0 -0.0
            -0.6767615 -0.0 0.0
            -0.6574254 0.0 0.0
            -0.6380894 -0.0 0.0
            -0.6187534 0.0 -0.0
            -0.5994173 -0.0 -0.0
            -0.5800813 -0.0 0.0
            -0.5607452 0.0 0.0
            -0.5414092 0.0 0.0
            -0.5220731 0.0 0.0
            -0.5027371 -0.0 -0.0
            -0.4834011 -0.0 -0.0
            -0.464065 0.0 -0.0
            -0.444729 0.0 0.0
            -0.4253929 0.0 0.0
            -0.4060569 -0.0 -0.0
            -0.3867209 -0.0 0.0
            -0.3673848 0.0 0.0
            -0.3480488 0.0 0.0
            -0.3287127 0.0 0.0
            -0.3093767 0.0 -0.0
            -0.2900406 0.0 -0.0
            -0.2707046 0.0 0.0
            -0.2513686 -0.0 -0.0
            -0.2320325 -0.0 0.0
            -0.2126965 0.0 0.0
            -0.1933604 0.0 0.0
            -0.1740244 -0.0 -0.0
            -0.1546883 0.0 -0.0
            -0.1353523 0.0 0.0
            -0.1160163 -0.0 0.0
            -0.0966802 0.0 0.0
            -0.0773442 0.0 0.0
            -0.0580081 -0.0 0.0
            -0.0386721 -0.0 0.0
            -0.019336 0.0 0.0
            0.0 0.0 0.0
            0.0 0.0 -0.0193656
            0.0 0.0 -0.0387312
            0.0 0.0 -0.0580968
            0.0 0.0 -0.0774623
            0.0 0.0 -0.0968279
            0.0 0.0 -0.1161935
            -0.0 -0.0 -0.1355591
            0.0 0.0 -0.1549247
            -0.0 -0.0 -0.1742903
            0.0 0.0 -0.1936559
            -0.0 -0.0 -0.2130214
            0.0 0.0 -0.232387
            0.0 0.0 -0.2517526
            -0.0 -0.0 -0.2711182
            -0.0 -0.0 -0.2904838
            -0.0 -0.0 -0.3098494
            0.0 0.0 -0.329215
            0.0 0.0 -0.3485805
            -0.0 -0.0 -0.3679461
            -0.0 -0.0 -0.3873117
            -0.0 -0.0 -0.4066773
            0.0 0.0 -0.4260429
            0.0 0.0 -0.4454085
            -0.0 -0.0 -0.4647741
            -0.0 -0.0 -0.4841397
            -0.0 -0.0 -0.5035052
            0.0 0.0 -0.5228708
            0.0 0.0 -0.5422364
            -0.0 -0.0 -0.561602
            0.0 0.0 -0.5809676
            -0.0 -0.0 -0.6003332
            -0.0 -0.0 -0.6196988
            -0.0 -0.0 -0.6390643
            0.0 0.0 -0.6584299
            -0.0 -0.0 -0.6777955
            0.0 0.0 -0.6971611
            -0.0 -0.0 -0.7165267
            -0.0 -0.0 -0.7358923
            0.0 0.0 -0.7552579
            -0.0 -0.0 -0.7746234
            0.0 0.0 -0.793989
            -0.0 -0.0 -0.8133546
            0.0 0.0 -0.8327202
            0.0 0.0 -0.8520858
            -0.8507859 -0.0 -0.0470488
            -0.8507859 0.0 0.0
            -0.4253929 -0.4253929 0.0
            -0.4253929 -0.4253929 -0.0195464
            -0.4253929 -0.4253929 -0.0390928
            -0.4253929 -0.4253929 -0.0586392
            -0.4253929 -0.4253929 -0.0781856
            -0.4253929 -0.4253929 -0.097732
            -0.4253929 -0.4253929 -0.1172784
            -0.4253929 -0.4253929 -0.1368248
            -0.4253929 -0.4253929 -0.1563712
            -0.4253929 -0.4253929 -0.1759176
            -0.4253929 -0.4253929 -0.1954641
            -0.4253929 -0.4253929 -0.2150105
            -0.4253929 -0.4253929 -0.2345569
            -0.4253929 -0.4253929 -0.2541033
            -0.4253929 -0.4253929 -0.2736497
            -0.4253929 -0.4253929 -0.2931961
            -0.4253929 -0.4253929 -0.3127425
            -0.4253929 -0.4253929 -0.3322889
            -0.4253929 -0.4253929 -0.3518353
            -0.4253929 -0.4253929 -0.3713817
            -0.4253929 -0.4253929 -0.3909281
            -0.4253929 -0.4253929 -0.4104745
            -0.4253929 -0.4253929 -0.4300209
            -0.4253929 -0.4253929 -0.4495673
            -0.4253929 -0.4051361 -0.4495673
            -0.4253929 -0.3848793 -0.4495673
            -0.4253929 -0.3646225 -0.4495673
            -0.4253929 -0.3443657 -0.4495673
            -0.4253929 -0.3241089 -0.4495673
            -0.4253929 -0.3038521 -0.4495673
            -0.4253929 -0.2835953 -0.4495673
            -0.4253929 -0.2633385 -0.4495673
            -0.4253929 -0.2430817 -0.4495673
            -0.4253929 -0.2228249 -0.4495673
            -0.4253929 -0.2025681 -0.4495673
            -0.4253929 -0.1823113 -0.4495673
            -0.4253929 -0.1620545 -0.4495673
            -0.4253929 -0.1417976 -0.4495673
            -0.4253929 -0.1215408 -0.4495673
            -0.4253929 -0.101284 -0.4495673
            -0.4253929 -0.0810272 -0.4495673
            -0.4253929 -0.0607704 -0.4495673
            -0.4253929 -0.0405136 -0.4495673
            -0.4253929 -0.0202568 -0.4495673
            -0.4253929 -0.0 -0.4495673
            -0.4120994 -0.0 -0.4355183
            -0.3988059 -0.0 -0.4214694
            -0.3855123 -0.0 -0.4074204
            -0.3722188 -0.0 -0.3933714
            -0.3589253 -0.0 -0.3793224
            -0.3456318 -0.0 -0.3652734
            -0.3323382 -0.0 -0.3512245
            -0.3190447 -0.0 -0.3371755
            -0.3057512 -0.0 -0.3231265
            -0.2924576 -0.0 -0.3090775
            -0.2791641 -0.0 -0.2950286
            -0.2658706 -0.0 -0.2809796
            -0.2525771 -0.0 -0.2669306
            -0.2392835 -0.0 -0.2528816
            -0.22599 -0.0 -0.2388326
            -0.2126965 -0.0 -0.2247837
            -0.1994029 -0.0 -0.2107347
            -0.1861094 -0.0 -0.1966857
            -0.1728159 -0.0 -0.1826367
            -0.1595224 -0.0 -0.1685877
            -0.1462288 -0.0 -0.1545388
            -0.1329353 -0.0 -0.1404898
            -0.1196418 -0.0 -0.1264408
            -0.1063482 -0.0 -0.1123918
            -0.0930547 -0.0 -0.0983429
            -0.0797612 -0.0 -0.0842939
            -0.0664676 -0.0 -0.0702449
            -0.0531741 -0.0 -0.0561959
            -0.0398806 -0.0 -0.0421469
            -0.0265871 -0.0 -0.028098
            -0.0132935 -0.0 -0.014049
            0.0 0.0 0.0
        ]' * u"bohr^-1"
    f2c = CartesianFromFractional(reci_lattice)
    @test maximum(abs.(f2c(qe_crystal) - qe_cartesian)) < 1e-5u"bohr^-1"
end


text:
module Primes

export prime_factors

_primes = [2, 3]

function is_prime(n, known_primes)
    for prime in known_primes
        if n == prime || prime > n
            return true
        end
        if n % prime == 0
            return false
        end
    end
    return true
end

function next_prime(n)
    inx = findfirst(x -> x > n, _primes)
    if inx == 0
        # Compute prime
        x = n + 1
        while !is_prime(x, _primes)
            x+=1
        end
        push!(_primes, x)
        return x
    end
    return _primes[inx]
end

function prime_factors(n, result=[])

    current_prime = 2
    while (n % current_prime) != 0
        current_prime = next_prime(current_prime)
    end

    if is_prime(n, _primes)
        push!(result, n)
        return result
    end
    push!(result, current_prime)

    return prime_factors(n / current_prime, result)
end

end # End of Module Primes

max_stars_repo_name: ddugue/project-euler

max_issues_repo_name: ddugue/project-euler

max_forks_repo_name: ddugue/project-euler