text:
module KineticSolverModule

import IterativeSolvers

using ProgressMeter, MyIOModule, DelimitedFiles, LinearAlgebra, Dates, Revise, PyPlot
using Einsum, Arpack,LinearMaps, SparseArrays
using DelimitedFiles,SphericalHarmonics, TypedPolynomials,GSL

export KineticSolver, computeMoments,corysolve!, solve!,sweeping
export setTestCaseInitialCondition!,computeHenyeyGreensteinScatteringMatrix,solveWithSweeping!,solveRSNWithSweeping!,sweepingWithKernels,vis,methods!
export  kerneldecomposition, metadecomposition,kerneldecomposition
export computeScatteringKernelHenyeyGreenstein, computeScatteringKernelConvolution
export checkexpansionHG, solveWithSweeping!, computeScatteringKernelConvolution

include("quadratures/Quadrature.jl")
include("problems/getProblemSpecificParameters.jl")

mutable struct KineticSolver
    # domain is given by [x0,x1] x [y0,y1]
    x0::Float64;x1::Float64; y0::Float64; y1::Float64
    dx::Float64 # corresponding equidistant grid spacing
    dy::Float64 # same for y
    nx::Int64 # number of cells +2 ghost cells
    ny::Int64 # same for y

    norder::Int64 # Order of the quadrature
    nquadpoints::Int64 # Resulting number of quadraturepoints
    #for SN: nquadpoints = norder^2, but for other quadratures this is not true

    # index array for the specified variables
    # later we just want to write for i=rangex instead of for i=1:N
    rangex::Array{Int64,1}
    rangey::Array{Int64,1}
    rangequad::Array{Int64,1}

    ThisQuadratureType::String # Name of the Quadrature being used
    Q::Quadrature              # Quadrature struct

    rotationmagnitude::Float64 # if rotation is chosen
    convolutionmagnitude::Float64 # convolution strength
    convolutionwidth::Float64 # controls how isotropically scattering happens
    convolutiontoscattering::Float64 # controls convolution strength in regions with high scattering

    cfl::Float64 # cfl number
    dt::Float64 # time step size
    tEnd::Float64  # time until we simulate
    nt::Int64 # number of time steps

    tmp::Array{Float64,1}
    periodicXflag::Bool # periodic x boundary
    periodicYflag::Bool
    dim1flag::Bool # perform 1d simulation
    rank::Int64

    SigmaS::Array{Float64,2} # Matrix of Scattering crossection
    SigmaT::Array{Float64,2} # Matrix of Total crossection

    ictype::String ## Name of Initial Condition
    Source::Array{Float64,3} # Tensor for time independent source
    psi0::Array{Float64,3} # Tensor for the phi(t=0,x,y,q)

    outputfolderprefix::String # folder name prefix for storing the data
    staticCounter::Int

    nSweepingSteps::Int
    epsSweeping::Float64

    DecompO::Array{Float64,2}
    DecompSigma::Array{Float64,2}
    DecompM::Array{Float64,2}
    OSigma::Array{Float64,2}
    lastpsi::Array{Float64,3}
    DecompNumberElements::Int64

    ScatteringKernelConvolution::Array{Float64,2}
    OutscatteringVector::Array{Float64,1}

    contractionEstimate::Float64
    kernelTimesWeightsSparse::SparseMatrixCSC{Float64,Int64}
    psiC::Array{Float64,3} # Tensor for the phi(t=0,x,y,q) for one shot

    function KineticSolver(configfilename::String)
        # read the following words from text file
        allowed = ["nx", "ny", "cfl", "quadraturetype", "quadratureorder",
                   "convolutionmagnitude", "convolutionwidth", "convolutiontoscattering",
                   "testcaseid", "rotationmagnitude",
                   "periodicXflag", "periodicYflag", "dim1flag", "whichrank"]
        D = readtodict(configfilename,allowed)

        #				 read data from config file start					  #
        nx = Int64(D["nx"])
        ny = Int64(D["ny"])
        norder = Int64(D["quadratureorder"])
        Q = Quadrature(norder,Int64(D["quadraturetype"]))
        rotationmagnitude = D["rotationmagnitude"]
        rank = Int64(D["whichrank"])
        cfl = D["cfl"]
        PossibleQuadratureTypes = ["tensorized", "octa", "ico", "dim2", "fromfile"]
        PossibleICTypes = ["linesource", "checkerboard", "dim1", "beam"]
        ThisQuadratureType = PossibleQuadratureTypes[Int64(D["quadraturetype"])]
        periodicXflag = Int64(D["periodicXflag"]) == 1 ? true : false
        periodicYflag = Int64(D["periodicYflag"]) == 1 ? true : false
        dim1flag = Int64(D["dim1flag"]) == 1 ? true : false
        ThisICType = PossibleICTypes[Int64(D["testcaseid"])]
        if dim1flag == true ## 1D FLAG FORCES THESE THINGS:
            ny = 3; periodicYflag = true;    ThisICType = "dim1"
            text = "You chose a 1d simulation, this forces: ny = 1, periodicYflag = true, ThisICType = 3, ThisQuadratureType = 4"
            printstyled(text,color = :red)
        end
        #some variables can be deduced immediately though not being part of the config
        nquadpoints = Q.nquadpoints
        rangex = collect(1:nx).+2 #range for x and y considering ghost cells
        rangey = collect(1:ny).+2
        rangequad = collect(1:nquadpoints);
        convolutionmagnitude = Float64(D["convolutionmagnitude"]);
        convolutionwidth = Float64(D["convolutionwidth"]);
        convolutiontoscattering = Float64(D["convolutiontoscattering"]);
        #				 read data from config file stop					  #

        ######## Get problem specific parameters start  ########
        x0,x1,y0,y1,tEnd,SigmaS,SigmaT,Source,psi0 = getProblemSpecificParameters(nx,ny,nquadpoints,ThisICType,Q.pointsxyz)
        # and deduce other variables from that
        dx = (x1 - x0) / nx
        dy = (y1 - y0) / ny
        dt = cfl / 2 * (dx*dy) / (dx + dy)
        nt = ceil(tEnd / dt);
        tmp = zeros(Float64,nquadpoints) # tmp array for in place mult
        ######## Get problem specific parameters stop  ########

        identifier = string(Dates.now(),"/")
        outputfolderprefix = string("../out/",identifier)
        mkdir(outputfolderprefix)
        mkdir(string(outputfolderprefix,"data/"))
        cp("config.txt",string(outputfolderprefix,"config.txt"))
        nicewrite(string(outputfolderprefix,"data/sigmaS.txt"),SigmaS)
        nicewrite(string(outputfolderprefix,"data/sigmaT.txt"),SigmaT)
        #nicewrite(string(outputfolderprefix,"data/sourceAngular.txt"),Source)
        #nicewrite(string(outputfolderprefix,"data/phi0.txt"),phi0)
        nicewrite(string(outputfolderprefix,"data/quadpoints.txt"),Q.pointsxyz)
        nicewrite(string(outputfolderprefix,"data/quadweights.txt"),Q.weights)

        # incomplete initialization:  weightMatrix will be computed at run time
        new(x0,x1,y0,y1,dx,dy,nx + 4,ny + 4,# spatial
            norder,nquadpoints,             # quadrature
            rangex,rangey,rangequad,        # range
            ThisQuadratureType, Q,          # quadrature
            rotationmagnitude,              # how much rotation
            convolutionmagnitude,           # how much convolution
            convolutionwidth,
            convolutiontoscattering,        # how much convolution in scattering/absorption regions
            cfl,dt,tEnd,nt,                 # time
            tmp,                            # tmp array for in place mult
            periodicXflag, periodicYflag, # are boundaries periodic?
            dim1flag,                     # 1d
            rank,                         # when performing low rank computations
            SigmaS,SigmaT,                # Sigma
            ThisICType,Source,psi0,       # IC/Q/PHI0
            outputfolderprefix,           # folder
            0
            )
    end # constructor
end # struct

include("vis.jl")
include("methods.jl")
include("sweeping.jl")
include("kerneldecomposition.jl")
include("solves.jl")
include("computeKernels.jl")

end # module

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
# This file is part of Kpax3. License is MIT.

function loss_binder(R::Vector{Int},
                     P::Matrix{Float64})
  n = length(R)

  loss = 0.0
  for j in 1:(n - 1), i in (j + 1):n
    loss += (R[i] == R[j]) ? 1 - P[i, j] : P[i, j]
  end

  loss
end

ext: jl


text:
using CautiousEngine

function test_intersect_2d()
    shape1 = Dict("vertices" => [[-1. -1.], [-1. 1.], [1. 1.], [1. -1.]])
    shape2 = Dict("vertices" => [[-.5 -1.], [-.5 1.], [-4. 1.], [-4. -1.]])
    axes = [1. 0.;0. 1.]
    res = CautiousEngine.do_shapes_intersect(shape1, shape2, axes=axes)
end

function test_intersect()
    shape1 = Dict("vertices" => [[-1. -1. -1.], [-1. 1. -1.], [1. 1. -1.], [1. -1. -1.],
                                 [-1. -1. 1.], [-1. 1. 1.], [1. 1. 1.], [1. -1. 1.]])
    shape2 = Dict("vertices" => [[-5. -2. -2.], [-5. 2. -2.], [-4. 2. -2.], [-4. -2. -2.],
                                 [-5. -2. 2.], [-5. 2. 2.], [-4. 2. 2.], [-4. -2. 2.]])
    axes = [1. 0. 0.;0. 1. 0.; 0. 0. 1]
    res = CautiousEngine.do_shapes_intersect(shape1, shape2, axes=axes)
return res
end

function test_containment()

    shape1 = Dict("vertices" => [[-1. -1. -1.], [-1. 1. -1.], [1. 1. -1.], [1. -1. -1.],
                                 [-1. -1. 1.], [-1. 1. 1.], [1. 1. 1.], [1. -1. 1.]])
    shape2 = Dict("vertices" => [[-2. -2. -2.], [-2. 2. -2.], [2. 2. -2.], [2. -2. -2.],
                                 [-2. -2. 2.], [-2. 2. 2.], [2. 2. 2.], [2. -2. 2.]])
    # axes = [1. 0. 0.;0. 1. 0.; 0. 0. 1]
    res = CautiousEngine.does_shape2_contain_shape1(shape1, shape2, axes=nothing)
    return res
end

@testset "region intersection tests" begin
    @test test_intersect_2d()[1]
	@test !test_intersect()[1]
    @test test_containment()[1]
end


text:
using .IntervalBasedForest
const IBF = IntervalBasedForest

"""
TimeSeriesForestClassifier(; kwargs...)

Overview: Input n series length m
for each tree
    sample sqrt(m) intervals
    find mean, sd and slope for each interval, concatenate to form new
    data set
    build decision tree on new data set
ensemble the trees with averaged probability estimates

## Hyperparameters

- `n_trees=200`            ensemble size(number of trees)

- `random_state`           seed for random numbers

- `min_interval=3`         minimum width of an interval

- `max_depth=-1`:          max depth of the decision tree (-1=any)

- `min_samples_leaf=1`:    max number of samples each leaf needs to have

- `min_samples_split=2`:   min number of samples needed for a split

- `min_purity_increase=0`: min purity needed for a split

- `n_subfeatures=0`: number of features to select at random (0 for all,
  -1 for square root of number of features)

- `post_prune=false`:      set to `true` for post-fit pruning

- `merge_purity_threshold=1.0`:  (post-pruning) merge leaves having `>=thresh`
                           combined purity

- `pdf_smoothing=0.0`:     threshold for smoothing the predicted scores

- `display_depth=5`:       max depth to show when displaying the tree

"""
@mlj_model mutable struct TimeSeriesForestClassifier<: MMI.Probabilistic # have another abstracttype?
    n_trees::Int                     =     200::(_  ≥ 1)
    random_state                     =  nothing
    min_interval::Int                =       3::(_  ≥ 1)
    max_depth::Int                   =  (-)(1)::(_ ≥ -1)
    min_samples_leaf::Int            =       1::(_ ≥ 0)
    min_samples_split::Int           =       2::(_ ≥ 2)
    min_purity_increase::Float64     =     0.0::(_ ≥ 0)
    n_subfeatures::Int               =       0::(_ ≥ -1)
    post_prune::Bool                 =   false
    merge_purity_threshold::Float64  =     1.0::(_ ≤ 1)
    pdf_smoothing::Float64           =     0.0::(0 ≤ _ ≤ 1)
    display_depth::Int               =       5::(_ ≥ 1)
end

function MMI.fit(m::TimeSeriesForestClassifier, verbosity::Int, Xmatrix, y)
    yplain = MMI.int(y)
    classes_seen  = filter(in(unique(y)), MMI.classes(y[1]))
    integers_seen = MMI.int(classes_seen)
    forest, intervals = IBF.TimeSeriesForestClassifier(m, Xmatrix, yplain)
    fitresult = (forest, intervals, classes_seen, integers_seen)
    cache  = nothing
    return fitresult, nothing, nothing
end

function MMI.predict(m::TimeSeriesForestClassifier, fitresult, Xmatrix)
    forest, intervals, classes_seen, integers_seen = fitresult
    scores = IBF.predict_new(Xmatrix, forest, intervals, integers_seen)
    sm_scores = smooth(scores)
    return [MMI.UnivariateFinite(classes_seen, view(sm_scores, i, :))
                   for i in 1:size(sm_scores, 1)]
end

function smooth(X)
    a, b = size(X)
    sm_scores = zeros(a, b)
    indices = map(i -> CartesianIndex(i, findmax(view(X, i, :))[2]), [1:a;])
    sm_scores[indices] .= 1
    return sm_scores
end

MMI.fitted_params(::TimeSeriesForestClassifier, fitresult) =
    (forest=fitresult[1])

#=
 NOTE: Should we add pretty printing?
=#

"""
An adapted version of the NearestNeighbors knn to work with
time series data.

## Hyperparameters

 -  `n_neighbors`            Int, set k for knn (default = 5)
 -  `weights`                mechanism for weighting a vote 'uniform', 'distance'
 -  `search_algorithm`       search method for neighbours default = "select_sort"
 -  `metric`                 distance measure for time series default = "dtw_distance"
 -  `metric_params`          dictionary for metric parameters

"""
@mlj_model mutable struct TimeSeriesKNNClassifier <: MMI.Probabilistic
    n_neighbors::Int                      =   5::(_  ≥ 1)
    weights::Symbol                       =   :uniform
    search_algorithm::Symbol              =   :select_sort
    metric::Symbol                        =   :dtw_distance
    metric_params::Array                  =   [-1.0]
end

function MMI.fit(m::TimeSeriesKNNClassifier, verbosity::Int, Xmatrix, y)
    yplain = MMI.int(y)
    classes_seen  = filter(in(unique(y)), MMI.classes(y[1]))
    integers_seen = MMI.int(classes_seen)
    fitresult = (Xmatrix, yplain, classes_seen, integers_seen)
    cache  = nothing
    return fitresult, nothing, nothing
end

function MMI.predict(m::TimeSeriesKNNClassifier, fitresult, Xmatrix_new)
    Xmatrix, yplain, classes_seen, integers_seen = fitresult
    y_pred, DistanceMatrix = Predict_new(m, Xmatrix, Xmatrix_new, yplain)
    a, b = length(y_pred), length(integers_seen)
    probas = zeros(a, b)
    @inbounds for i=1:a
        @inbounds for j=1:b
            probas[i,j] = Int(y_pred[i]) == j ? 1 : 0
        end
    end
    return  [MMI.UnivariateFinite(classes_seen, view(probas, i, :))
                for i in 1:size(probas, 1)]
end

# MMI.fitted_params(::TimeSeriesKNNClassifier, fitresult)


text:
using Revise
using JuMP
using GLPK
using ModelGraphs
using .DualDecompositionSolver

const DDSolver = DualDecompositionSolver

m1 = Model(with_optimizer(GLPK.Optimizer))
@variable(m1, xm[i in 1:2],Bin)
@constraint(m1, xm[1] + xm[2] <= 1)
@objective(m1, Max, 16xm[1] + 10xm[2])

m2 = Model(with_optimizer(GLPK.Optimizer))
@variable(m2, xs[i in 1:2], Bin)
@variable(m2, y[i in 1:2], Bin)
@constraint(m2, y[1] + y[2] <= 1)
@constraint(m2, 8xs[1] + 2xs[2] + y[1] + 4y[2] <= 10)
@objective(m2, Max, 4y[2])

m3 = Model(with_optimizer(GLPK.Optimizer))
@variable(m3, x3[i in 1:2], Bin)
@variable(m3, y[i in 1:2], Bin)
@constraint(m3, y[1] + y[2] <= 1)
@constraint(m3, 8x3[1] + 2x3[2] + y[1] + 4y[2] <= 10)
@objective(m3, Max, y[2])

## Model Graph
graph = ModelGraph()

#Add nodes
n1 = add_node!(graph,m1)
n2 = add_node!(graph,m2)
n3 = add_node!(graph,m3)

#link constraints between models
@linkconstraint(graph, [i in 1:2], n1[:xm][i] == n2[:xs][i])
@linkconstraint(graph, n3[:x3][1] + n1[:xm][1] + n2[:xs][1] <= 2)
@linkconstraint(graph, n1[:xm][2] <= n3[:y][1])

dd_model = DDModel(graph)
DDSolver.optimize!(dd_model)

#print solution
println("Dual decomposition solution: ")
solution = Float64[]
for subproblem in dd_model.subproblems
    append!(solution,value.(all_variables(subproblem)))
end
println(solution)

# #verify with full problem
glpk = with_optimizer(GLPK.Optimizer)
m,ref_map = aggregate(graph)
optimize!(m,glpk)

println()
println("Pure glpk solution:")
println(JuMP.value.(all_variables(m)))


text:
# Scanning

#################################################
#
#   Scanning along a vector
#
#################################################

macro code_scan(AN)

    h = codegen_helper(AN)
    aparams = h.contiguous_aparams
    args = h.args
    t1 = h.term(1)
    ti = h.term(:i)

    quote
        # generic scanning functions

        global _scan!
        function _scan!(s, r::AbstractArray, op::Functor{2}, $(aparams...))
            i = 1
            @inbounds r[1] = s

            for i = 2 : length(r)
                @inbounds vi = $(ti)
                s = evaluate(op, s, vi)
                @inbounds r[i] = s
            end
        end

        global scan!
        function scan!(r::AbstractArray, op::Functor{2}, $(aparams...))
            n = $(h.inputlen)
            if n > 0
                @inbounds s = $(t1)
                _scan!(s, r, op, $(args...))
            end
            return r
        end

        global scan
        function scan(op::Functor{2}, $(aparams...))
            shp = $(h.inputsize)
            n = prod(shp)
            if n > 0
                @inbounds s = $(t1)
                r = Array(typeof(s), shp)
                _scan!(s, r, op, $(args...))
                return r
            else
                Any[]
            end
        end

        # specific scanning functions
        global cumsum!, cummax!, cummin!
        cumsum!(r::AbstractArray, $(aparams...)) = scan!(r, Add(), $(args...))
        cummax!(r::AbstractArray, $(aparams...)) = scan!(r, MaxFun(), $(args...))
        cummin!(r::AbstractArray, $(aparams...)) = scan!(r, MinFun(), $(args...))

        global cumsum, cummax, cummin
        cumsum($(aparams...)) = scan(Add(), $(args...))
        cummax($(aparams...)) = scan(MaxFun(), $(args...))
        cummin($(aparams...)) = scan(MinFun(), $(args...))
    end
end

@code_scan 0
@code_scan 1
@code_scan 2
@code_scan 3

# inplace scanning

scan!(op::Functor{2}, r::ContiguousArray) = scan!(r, op, r)

cumsum!(r::ContiguousNumericArray) = scan!(Add(), r)
cummax!(r::ContiguousNumericArray) = scan!(MaxFun(), r)
cummin!(r::ContiguousNumericArray) = scan!(MinFun(), r)

#################################################
#
#   Scanning along a specific dimension
#
#################################################

macro code_scandim(AN)

    h = codegen_helper(AN)
    aparams = h.contiguous_aparams
    args = h.args
    tidx = h.term(:idx)

    if AN == 0
        offset_args = [:(offset_view(a, ao, m, n))]
    else
        offset_args = [:fun, [:(offset_view($a, ao, m, n)) for a in args[2:end]]...]
    end

    quote
        # generic scanning functions

        global _scan_eachcol!
        function _scan_eachcol!(m::Int, n::Int, r::AbstractArray, op::Functor{2}, $(aparams...))
            o = 0
            for j = 1 : n
                idx = o + 1
                @inbounds s = $(tidx)
                @inbounds r[idx] = s

                for i = 2 : m
                    idx = o + i
                    @inbounds vi = $(tidx)
                    s = evaluate(op, s, vi)
                    @inbounds r[idx] = s
                end
                o += m
            end
        end

        global _scan_eachrow!
        function _scan_eachrow!(m::Int, n::Int, r::AbstractArray, op::Functor{2}, $(aparams...))
            o = 0
            for idx = 1 : m
                @inbounds r[idx] = $(tidx)
            end

            o = m
            for j = 2 : n
                for i = 1 : m
                    idx = o + i
                    @inbounds s = r[idx - m]
                    @inbounds vi = $(tidx)
                    s = evaluate(op, s, vi)
                    @inbounds r[idx] = s
                end
                o += m
            end
        end

        global _scan!
        function _scan!(r::ContiguousArray, op::Functor{2}, $(aparams...), dim::Int)
            if !isempty(r)
                shp = size(r)
                if dim == 1
                    m = shp[1]
                    n = succ_length(shp, 1)
                    _scan_eachcol!(m, n, r, op, $(args...))

                else
                    m = prec_length(shp, dim)
                    n = shp[dim]
                    k = succ_length(shp, dim)

                    if k == 1
                        _scan_eachrow!(m, n, r, op, $(args...))
                    else
                        mn = m * n
                        ro = 0
                        ao = 0
                        for l = 1 : k
                            _scan_eachrow!(m, n, offset_view(r, ro, m, n), op, $(offset_args...))
                            ro += mn
                            ao += mn
                        end
                    end
                end
            end
            return r
        end

        global scan!
        function scan!(r::ContiguousArray, op::Functor{2}, $(aparams...), dim::Int)
            shp = $(h.inputsize)
            size(r) == shp || error("Invalid argument dimensions.")
            _scan!(r, op, $(args...), dim)
        end

        global scan
        function scan(op::Functor{2}, $(aparams...), dim::Int)
            shp = $(h.inputsize)
            tt = $(h.termtype)
            rt = result_type(op, tt, tt)
            _scan!(Array(rt, shp), op, $(args...), dim)
        end

        # specific scanning functions
        global cumsum!, cummax!, cummin!
        cumsum!(r::AbstractArray, $(aparams...), dim::Int) = scan!(r, Add(), $(args...), dim)
        cummax!(r::AbstractArray, $(aparams...), dim::Int) = scan!(r, MaxFun(), $(args...), dim)
        cummin!(r::AbstractArray, $(aparams...), dim::Int) = scan!(r, MinFun(), $(args...), dim)

        global cumsum, cummax, cummin
        cumsum($(aparams...), dim::Int) = scan(Add(), $(args...), dim)
        cummax($(aparams...), dim::Int) = scan(MaxFun(), $(args...), dim)
        cummin($(aparams...), dim::Int) = scan(MinFun(), $(args...), dim)
    end
end

@code_scandim 0
@code_scandim 1
@code_scandim 2
@code_scandim 3

# inplace scanning

scan!(op::Functor{2}, r::ContiguousArray, dim::Int) = scan!(r, op, r, dim)

cumsum!{T<:Number}(r::ContiguousArray{T}, dim::Int) = scan!(Add(), r, dim)
cummax!{T<:Real}(r::ContiguousArray{T}, dim::Int) = scan!(MaxFun(), r, dim)
cummin!{T<:Real}(r::ContiguousArray{T}, dim::Int) = scan!(MinFun(), r, dim)


text:
# * Binary utilities

function showbin(io::IO, c::BitVector)
    # Printing order = LSB ... MSB
    for (i,e) in enumerate(c)
        write(io, e ? "1" : ".")
        i%4 == 0 && write(io, " ")
    end
end
showbin(c::BitVector) = showbin(stdout, c)

showbin(io::IO, i::Integer) = showbin(io, BitVector(digits(i, base=2)))

struct BitPattern{O}
    obj::O
end

Base.show(io::IO, b::BitPattern) = showbin(io, b.obj)

macro showbin(var)
    name = string(var)
    quote
        println($name, " = ", BitPattern($(esc(var))))
    end
end

# * Orbitals

struct Orbitals{O,Overlaps,HasOverlap,NonOrthogonalities}
    orbitals::Vector{O}
    overlaps::Overlaps
    has_overlap::HasOverlap
    non_orthogonalities::NonOrthogonalities
    function Orbitals(orbitals::Vector{O}, overlaps::AbstractVector{<:OrbitalOverlap}) where O
        orb_map = Dict(o => i for (i,o) in enumerate(orbitals))
        n = length(orbitals)
        S = spdiagm(n,n,0 => fill(one(NBodyTerm), n))
        for oo in overlaps
            i = orb_map[oo.a]
            j = orb_map[oo.b]
            S[i,j] = OrbitalOverlap(oo.a, oo.b)
            S[j,i] = OrbitalOverlap(oo.b, oo.a)
        end
        N = falses(n,n)
        for oo in overlaps
            i = orb_map[oo.a]
            j = orb_map[oo.b]
            N[i,j] = true
            N[j,i] = true
        end
        has_overlap =  N .| Matrix(I, size(N))
        non = map(|, vec(reduce(|, N, dims=1)), vec(reduce(|, N, dims=2)))
        new{O,typeof(S),BitMatrix,BitVector}(orbitals, S, has_overlap, non)
    end
    function Orbitals(orbitals::Vector{O}) where O
        no = length(orbitals)
        has_overlap = Matrix(I,no,no)
        new{O,UniformScaling{Bool},BitMatrix,Nothing}(orbitals, I, has_overlap, nothing)
    end
end

Base.length(orbitals::Orbitals) = length(orbitals.orbitals)
Base.eachindex(orbitals::Orbitals) = eachindex(orbitals.orbitals)
Base.getindex(orbitals::Orbitals, i) = orbitals.orbitals[i]

# * BitConfiguration

struct BitConfiguration{Orbitals,Occupations}
    orbitals::Orbitals
    occupations::Occupations
end

Base.show(io::IO, c::BitConfiguration, sel=eachindex(c.orbitals)) =
    write(io, join(string.(c.orbitals[filter(∈(sel), orbital_numbers(c))]), " "))

showbin(io::IO, c::BitConfiguration) = showbin(io, c.occupations)
showbin(c::BitConfiguration) = showbin(c.occupations)

Base.:(==)(a::BitConfiguration, b::BitConfiguration) =
    a.orbitals == b.orbitals &&
    a.occupations == b.occupations

num_particles(c::BitConfiguration) = count(c.occupations)
orbital_numbers(c::BitConfiguration) = findall(c.occupations)
get_orbitals(c::BitConfiguration) = c.orbitals[orbital_numbers(c)]

function configuration_orbital_numbers(c::BitConfiguration, mask::BitVector)
    occ = copy(c.occupations)
    n = 0
    ons = Vector{Int}()
    while !iszero(occ)
        occ[1] && (n += 1)
        mask[1] && push!(ons, n)
        occ <<= 1
        mask <<= 1
    end
    ons
end

for f in (:(&), :(|), :xor)
    @eval begin
        Base.$(f)(a::BitConfiguration, b::BitConfiguration) = map($f, a.occupations, b.occupations)
        Base.$(f)(a::BitConfiguration, b::BitVector) = map($f, a.occupations, b)
        Base.$(f)(a::BitVector, b::BitConfiguration) = map($f, a, b.occupations)
    end
end

#=
Adapted from

- Scemama, A., & Giner, E. (2013). An efficient implementation of
  Slater–Condon rules. CoRR, arXiv:1311.6244.

=#

num_excitations(a::BitConfiguration, b::BitConfiguration) = count(a ⊻ b)

holes_mask(a::BitConfiguration, b::BitConfiguration) = (a ⊻ b) & a
holes(a::BitConfiguration, b::BitConfiguration) = findall(holes_mask(a, b))

particles_mask(a::BitConfiguration, b::BitConfiguration) = (a ⊻ b) & b
particles(a::BitConfiguration, b::BitConfiguration) = findall(particles_mask(a, b))

function holes_particles_mask(a::BitConfiguration, b::BitConfiguration)
    x = a ⊻ b
    x, x & a, x & b
end

function holes_particles(a::BitConfiguration, b::BitConfiguration)
    x, hm, pm = holes_particles_mask(a, b)
    findall(hm), findall(pm)
end

function relative_sign(a::BitConfiguration, b::BitConfiguration,
                       holes, particles; verbosity=0)
    isempty(holes) && return 1

    perm = 0
    av = copy(a.occupations)
    n = length(av)

    if verbosity > 0
        print("   ")
        @showbin a
        print("   ")
        @showbin b
    end

    @inbounds for (h,p) in zip(holes, particles)
        lo,hi = minmax(h,p)
        mask = vcat(falses(lo), trues(hi-lo-1), falses(n-hi+1))
        perm += count(map(&, mask, av))
        if verbosity > 0
            @show h,p,lo,hi
            ma = map(&, mask, av)
            print("  ")
            @showbin av
            @showbin mask
            print("  ")
            @showbin ma
        end
        av[h] = false
        av[p] = true
    end
    @assert av == b.occupations

    iseven(perm) ? 1 : -1
end

# * BitConfigurations

struct BitConfigurations{Orbitals}
    orbitals::Orbitals
    configurations::BitMatrix
end

get_orbitals(c::Configuration) = c.orbitals
get_orbitals(v::AbstractVector) = v

function BitConfigurations(configurations::AbstractVector,
                           overlaps=OrbitalOverlap[])
    orbitals = unique(reduce(vcat, map(get_orbitals, configurations)))
    orb_map = Dict(o => i for (i,o) in enumerate(orbitals))
    C = falses(length(orbitals), length(configurations))
    p = Progress(prod(size(C)))
    for (j,c) in enumerate(configurations)
        for o in get_orbitals(c)
            C[orb_map[o],j] = true
            ProgressMeter.next!(p)
        end
    end
    BitConfigurations(Orbitals(orbitals, overlaps), C)
end

Base.getindex(m::BitConfigurations, i::Integer) = BitConfiguration(m.orbitals, m.configurations[:,i])
Base.getindex(m::BitConfigurations, i) = BitConfigurations(m.orbitals, m.configurations[:,i])
Base.length(m::BitConfigurations) = size(m.configurations, 2)
Base.isempty(m::BitConfigurations) = length(m) == 0

function AtomicLevels.core(m::BitConfigurations)
    isempty(m) && return 1:0
    # Find all initial orbitals which are occupied in all
    # configurations.
    i = findfirst(.!vec(reduce(&, m.configurations, dims=2))) - 1
    sel = 1:i
    if i > 0 && !isnothing(m.orbitals.overlaps)
        # We require that the core orbitals are all canonical ones,
        # i.e. no non-orthogonalities are allowed. This simplifies the
        # energy expression, since the core can then only contribute
        # through Slater–Condon rules.
        non = m.orbitals.non_orthogonalities[sel]
        imax = sel[end]
        sel = 1:min(imax, something(findfirst(non), imax+1))
    end
    sel
end

function Base.show(io::IO, m::BitConfigurations)
    norb = length(m.orbitals)
    ncfg = length(m)
    write(io, "$(norb)-orbital $(ncfg)-configuration BitConfigurations")
end

function Base.show(io::IO, ::MIME"text/plain", m::BitConfigurations)
    show(io, m)
    println(io)
    c = core(m)
    !isempty(c) && write(io, "Common core: $(c)")
    isempty(m) && return
    println(io)
    ncfg = length(m)
    fmt = FormatExpr("{1:$(length(digits(ncfg)))d}: ")
    printfmt(io, fmt, 1)
    ref = m[1]
    sel = (isempty(c) ? 0 : c[end])+1:length(m.orbitals)
    show(io, ref, sel)
    show_conf = i -> begin
        println(io)
        exc = m[i]
        h,p = holes_particles(ref, exc)
        printfmt(io, fmt, i)
        write(io, join(string.(m.orbitals[h]), " "))
        write(io, " -> ")
        write(io, join(string.(m.orbitals[p]), " "))
    end
    screen_rows = fld(max(3,first(displaysize(io))-6), 2)
    ncfg = length(m)
    first_block = 2:min(ncfg, screen_rows)
    second_block = max(1,ncfg - screen_rows):ncfg
    if isempty((2:first_block[end]+2) ∩ second_block)
        foreach(show_conf, first_block)
        print(io, "\n⋮")
        foreach(show_conf, second_block)
    else
        foreach(show_conf, 2:second_block[end])
    end
end

function orbital_overlap_matrix(sd::BitConfigurations, i, j)
    a = sd[i]
    b = sd[j]

    sd.orbitals.overlaps[orbital_numbers(a), orbital_numbers(b)]
end

"""
    non_zero_cofactors(sd, N, i, j)

Find all non-zero cofactors of the orbital overlap matrix between the
Slater determinants `i` & `j` of `sd`, generated when striking out `N`
rows & columns. This routine is tailored towards the case when few
non-orthogonalities are present, e.g. approximately proportional to
the number of orbitals.

Non-orthogonality between spin-orbitals is handled by dividing the
them into two subspaces:

1. The orthogonal spin-orbitals that are common to both Slater
   determinants (core orbitals),

2. All non-orthogonal orbitals, and the orbitals which differ between
   the Slater determinants (i.e. holes of `i` and particles of `j`).

The relative phase between the Slater determinants is determined by
group `2` alone, by permuting the particles to the positions of the
holes, we find this phase. We can then formally permute them together
to a diagonal block at lower-right corner of the orbital overlap
matrix without incurring a phase change, since we need to permute the
same number of rows and columns. We thus get this structure:

                     ╷       ╷
                     │ 1 │   │
     det(Sᵢⱼ) = (-)ᵏ │───┼───│
                     │   │ 2 │
                     ╵       ╵

where `k` is decided by the permutation necessary to put the particles
in the positions of the holes.

Obviously, the determinant of the orbital matrix is now given by
`det(Sᵢⱼ) = (-)ᵏ*det(2)`, since we trivially have `det(1)==1`.

Depending on the rank of `2` (determined by the number of
hole–particle pairs and which spin-orbitals are non-orthogonal), we
need to strike out at least `size(2,1)-rank(2)` rows/columns from `2`,
and at most `min(N,size(2,1))`, i.e. for each value of `n ∈
size(2,1)-rank(2):min(N,size(2,1))`, we need to additionally strike
out `m = N - n` rows from `1`, but since the determinant of subspace
`1` is unity, regardless of how many rows/columns we've stricken out,
this is a trivial excercise. Of course, we also require that `m ≤
size(1,1)`.
"""
function non_zero_cofactors(sd::BitConfigurations, N, i, j; verbosity=0)
    verbosity > 0 && @show N
    a = sd[i]
    b = sd[j]
    np = num_particles(a)
    np == num_particles(b) ||
        throw(ArgumentError("Configurations do not have the same number of particles"))

    ks = Vector{Vector{Int}}()
    ls = Vector{Vector{Int}}()
    Ds = Vector{NBodyMatrixElement}()

    # First deal with trivial special cases.
    N > np && return ks, ls, Ds

    aon = orbital_numbers(a)
    bon = orbital_numbers(b)

    if N == np
        push!(ks, aon)
        push!(ls, bon)
        push!(Ds, one(NBodyMatrixElement))
        return ks, ls, Ds
    end

    isnothing(sd.orbitals.non_orthogonalities) &&
        error("Straight Slater–Condon not yet implemented")

    if verbosity > 0
        S = orbital_overlap_matrix(sd, i, j)
        display(S)
        i == j && @show det(S)
    end

    v = sd.orbitals.non_orthogonalities
    excitations, hmask, pmask = holes_particles_mask(a, b)

    common = a & b
    # For these, we may employ Slater–Condon rules
    canonical_orbitals = findall(map(&, common, .~v))

    # For these, we have to apply the Löwdin rules, but we can
    # precompute the determinant of this subspace once, for all those
    # determinant minors which do no strike out rows/columns in this
    # subspace.
    x = map(|, v, excitations)
    anon = findall(a & x)
    bnon = findall(b & x)

    h = findall(hmask)
    p = findall(pmask)

    S2 = sd.orbitals.overlaps[anon, bnon]
    non2 = sd.orbitals.has_overlap[anon, bnon]
    if verbosity > 0
        display(S2)
        display(sparse(non2))
        isempty(S2) && @info "We're essentially in Slater–Condon land"
    end

    # Minimum order of Γ⁽ᵖ⁾
    prmin = count(iszero, reduce(|, non2, dims=2))
    pcmin = count(iszero, reduce(|, non2, dims=1))
    pmin = max(prmin, pcmin)

    if N < pmin
        verbosity > 0 && @info "Too many vanishing rows/columns"
        return ks, ls, Ds
    end

    verbosity > 0 && @info "We need to strike out $(pmin)..$(min(N, length(anon))) rows/columns from S2"

    rs = relative_sign(a, b, h, p, verbosity=verbosity)
    verbosity > 0 && @show rs

    if N == 0
        # As a special case, for the zero-body operator, return the
        # properly phased determinant of S2.
        push!(Ds, rs*det(S2))
        return ks, ls, Ds
    end

    for n = pmin:min(N, length(anon))
        # n is many rows/columns we're striking out from S2;
        # consequently m is how many rows/columns we have to strike
        # out from S1
        m = N - n
        verbosity > 0 && @show n, m
        kk,ll = nonzero_minors(n, S2)
        for (k,l) in zip(kk,ll)
            verbosity > 0 && @show k,l
            Dkl = cofactor(k, l, S2)
            iszero(Dkl) && continue
            verbosity > 0 && @show Dkl
            # We now employ the Slater–Condon rules for S1
            for co in combinations(canonical_orbitals, m)
                ko = vcat(anon[k], co)
                lo = vcat(bnon[l], co)
                push!(ks, ko)
                push!(ls, lo)
                push!(Ds, rs*Dkl)
            end
        end
    end

    ks, ls, Ds
end

function NBodyMatrixElement(bcs::BitConfigurations, op::NBodyOperator{N}, i, j) where N
    orbitals = bcs.orbitals.orbitals
    ks,ls,Ds = non_zero_cofactors(bcs, N, i, j)
    NBodyMatrixElement(orbitals, op, orbitals,
                       zip(ks, ls, Ds))
end

function NBodyMatrixElement(bcs::BitConfigurations, op::LinearCombinationOperator, i, j)
    terms = NBodyTerm[]
    for (o,coeff) in op.operators
        iszero(coeff) && continue
        nbme = coeff*NBodyMatrixElement(bcs, o, i, j)
        !iszero(nbme) && append!(terms, nbme.terms)
    end
    NBodyMatrixElement(terms)
end

function Base.Matrix(bcs::BitConfigurations, rows, op::QuantumOperator, cols; verbosity=0)
    m,n = length(rows), length(cols)

    I = [Int[] for i in 1:Threads.nthreads()]
    J = [Int[] for i in 1:Threads.nthreads()]
    V = [NBodyMatrixElement[] for i in 1:Threads.nthreads()]

    p = if verbosity > 0
        @info "Generating energy expression"
        Progress(m*n)
    end

    Threads.@threads for i in eachindex(rows)
        r = rows[i]
        tid = Threads.threadid()
        for (j,c) in enumerate(cols)
            me = NBodyMatrixElement(bcs, op, r, c)

            !isnothing(p) && ProgressMeter.next!(p)

            iszero(me) && continue

            push!(I[tid], i)
            push!(J[tid], j)
            push!(V[tid], me)
        end
    end

    sparse(reduce(vcat, I),
           reduce(vcat, J),
           reduce(vcat, V), m, n)
end

function Base.Matrix(bcs::BitConfigurations, op::QuantumOperator; kwargs...)
    m = length(bcs)
    Matrix(bcs, 1:m, op, 1:m; kwargs...)
end

export BitConfigurations


text:
# using Revise
using Test, BifurcationKit, LinearAlgebra, SparseArrays, Arpack
const BK = BifurcationKit
####################################################################################################
# test the type BorderedArray and the different methods associated to it
z_pred = BorderedArray(rand(10),1.0)
tau_pred = BorderedArray(rand(10),2.0)
BK.minus!(z_pred, tau_pred)
BK.eltype(z_pred)

axpy!(2. /3, tau_pred, z_pred)
axpby!(2. /3, tau_pred, 1.0, z_pred)
dot(z_pred, tau_pred)

dottheta = BK.DotTheta((x, y) -> dot(x, y) / length(x))

dottheta(z_pred, 0.1)
dottheta(z_pred, tau_pred, 0.1)
dottheta(z_pred.u, tau_pred.u, 1.0, 1.0, 0.1)

z = BorderedArray(z_pred, rand(10))
z2 = BorderedArray(z_pred, rand(10))
zero(z2);zero(z_pred)
@test length(z_pred) == 11

copyto!(z,z2)
BK.minus(z.u,z2.u);BK.minus!(z.u,z2.u)
BK.minus(1.,2.);BK.minus!(1.,2.)
rmul!(z_pred, 1.0)
rmul!(z_pred, true)
mul!(z_pred, tau_pred, 1.0)

z_predC = BorderedArray(ComplexF64.(z_pred.u), ComplexF64.(z_pred.u))
z3 = similar(z_predC, ComplexF64)
mul!(z3, z3, 1.0)

z_sim = BorderedArray(rand(3), rand(2))
z_sim2 = similar(z_sim)
typeof(z_sim) == typeof(z_sim2)
####################################################################################################
# test of the linear  solvers
J0 = rand(100,100) * 0.1 - I
rhs = rand(100)
sol_explicit = J0 \ rhs

ls = DefaultLS()
_sol, = ls(J0, rhs)
@test _sol == sol_explicit

ls = GMRESIterativeSolvers(N=100, reltol = 1e-16)
_sol, = ls(J0, rhs)
@test _sol ≈ sol_explicit

ls = GMRESKrylovKit(rtol = 1e-16)
_sol, = ls(J0, rhs)
@test _sol ≈ sol_explicit

# Case with a shift
sol_explicit = (0.9J0 + 0.1I) \ rhs

ls = DefaultLS()
_sol, = ls(J0, rhs; a₀ = 0.1, a₁ = 0.9)
@test _sol == sol_explicit

ls = GMRESIterativeSolvers(N=100, reltol = 1e-16)
_sol, = ls(J0, rhs; a₀ = 0.1, a₁ = 0.9)
@test _sol ≈ sol_explicit

ls = GMRESKrylovKit(rtol = 1e-16)
_sol, = ls(J0, rhs; a₀ = 0.1, a₁ = 0.9)
@test _sol ≈ sol_explicit
####################################################################################################
# test the bordered linear solvers
J0 = rand(100,100) * 0.9 - I
rhs = rand(100)
sol_explicit = (J0 + 0.2spdiagm(0 => vcat(ones(99),0))) \ rhs

linBdsolver = BK.BorderingBLS(solver = DefaultLS(), checkPrecision=true)
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd1u
@test sol_explicit[end] ≈ sol_bd1p

ls = GMRESIterativeSolvers(reltol = 1e-11, N = length(rhs)-1)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

ls = GMRESKrylovKit(dim = length(rhs) - 1, rtol = 1e-11)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

linBdsolver = BK.MatrixBLS(ls)
sol_bd3u, sol_bd3p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd3u
@test sol_explicit[end] ≈ sol_bd3p

BK.MatrixFreeBLS(nothing)
linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd3u, sol_bd3p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd3u
@test sol_explicit[end] ≈ sol_bd3p rtol = 1e-6

linBdsolver = BK.MatrixFreeBLS(GMRESIterativeSolvers(reltol = 1e-9, N = size(J0, 1)))
sol_bd4u, sol_bd4p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd4u
@test sol_explicit[end] ≈ sol_bd4p

# test the bordered linear solvers as used in newtonPALC
xiu = rand()
xip = rand()

linBdsolver = BK.BorderingBLS(DefaultLS())
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end], xiu, xip)

linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J0[1:end-1,1:end-1], J0[1:end-1,end], J0[end,1:end-1], J0[end,end], rhs[1:end-1], rhs[end], xiu, xip)

@test sol_bd1u ≈ sol_bd2u
@test sol_bd1p ≈ sol_bd2p
####################################################################################################
# test the bordered linear solvers, Complex case
# we test the linear system with MA formulation of Hopf bifurcation
J = [0 1 0; -1 0 0; 0 0 1.]
V, Ve = eigen(J)
λ = V[2]
v = Ve[:,2]
V, Ve = eigen(J')
λ2 = V[1]
@assert λ ≈ conj(λ2)
w = Ve[:,1]
@test det(J - λ*I(3)) ≈ 0

J0 = [[J - λ*I(3) w];[(v') 0]]
@test ~(det(J0) ≈ 0)
@test det([[J - λ*I(3) w];[conj(v') 0]]) ≈ 0

n = size(J, 1)
rhs = zeros(ComplexF64, n+1)
sol_explicit = J0 \ rhs

# extract the border
J11 = J0[1:end-1,1:end-1]
J12 = J0[1:end-1,end]
J21 = J0[end,1:end-1] |> conj
J22 = J0[end,end]

@test norm([[J11 J12]; [J21' J22]] - J0) ≈ 0

J0_b = hcat(J11, J12); J0_b = vcat(J0_b, hcat(J21' , J22))
    @test norm(J0_b - J0) ≈ 0

linBdsolver = BK.MatrixBLS()
    sol_bd3u, sol_bd3p, _, _ = @time linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
    @test sol_explicit[1:end-1] ≈ sol_bd3u
    @test sol_explicit[end] ≈ sol_bd3p

# without the shift the frist block is singular
linBdsolver = BK.BorderingBLS(solver = DefaultLS(), checkPrecision=true)
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd1u
@test sol_explicit[end] ≈ sol_bd1p

ls = GMRESIterativeSolvers(reltol = 1e-11, N = length(rhs)-1)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

ls = GMRESKrylovKit(dim = length(rhs) - 1, rtol = 1e-11)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

BK.MatrixFreeBLS(nothing)
linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd3u, sol_bd3p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd3u
@test sol_explicit[end] ≈ sol_bd3p rtol = 1e-6

linBdsolver = BK.MatrixFreeBLS(GMRESIterativeSolvers(reltol = 1e-9, N = size(J0, 1)))
sol_bd4u, sol_bd4p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end]; shift = 0.2)
@test sol_explicit[1:end-1] ≈ sol_bd4u
@test sol_explicit[end] ≈ sol_bd4p

# test the bordered linear solvers as used in newtonPALC
ξu = 1.
ξp = rand()

linBdsolver = BK.BorderingBLS(DefaultLS())
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end], ξu, ξp; shift = 0.2)

linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end], ξu, ξp; shift = 0.2)

@test sol_bd1u ≈ sol_bd2u
@test sol_bd1p ≈ sol_bd2p
####################################################################################################
# test the bordered linear solvers, Complex case
# random case
n = 5
J0 = rand(n,n) * 0.9 - I
J0 = hcat(J0, rand(ComplexF64, n))
J0 = vcat(J0, rand(ComplexF64, n+1)')
rhs = rand(ComplexF64, size(J0, 1))
sol_explicit = J0 \ rhs

# extract the border
J11 = real.(J0[1:end-1,1:end-1])
J12 = J0[1:end-1,end]
J21 = J0[end,1:end-1] |> conj
J22 = J0[end,end]

@test norm([[J11 J12]; [J21' J22]] - J0) ≈ 0

J0_b = hcat(J11, J12); J0_b = vcat(J0_b, hcat(J21' , J22))
    @test norm(J0_b - J0) ≈ 0

linBdsolver = BK.MatrixBLS()
sol_bd3u, sol_bd3p, _, _ = @time linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd3u
@test sol_explicit[end] ≈ sol_bd3p

linBdsolver = BK.BorderingBLS(solver = DefaultLS(), checkPrecision=true)
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd1u
@test sol_explicit[end] ≈ sol_bd1p

ls = GMRESIterativeSolvers(reltol = 1e-11, N = length(rhs)-1)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

ls = GMRESKrylovKit(dim = length(rhs) - 1, rtol = 1e-11)
linBdsolver = BK.BorderingBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd2u
@test sol_explicit[end] ≈ sol_bd2p

BK.MatrixFreeBLS(nothing)
linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd3u, sol_bd3p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd3u
@test sol_explicit[end] ≈ sol_bd3p rtol = 1e-6

linBdsolver = BK.MatrixFreeBLS(GMRESIterativeSolvers(reltol = 1e-9, N = size(J0, 1)))
sol_bd4u, sol_bd4p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end])
@test sol_explicit[1:end-1] ≈ sol_bd4u
@test sol_explicit[end] ≈ sol_bd4p

# test the bordered linear solvers as used in newtonPALC
xiu = rand()
xip = rand()

linBdsolver = BK.BorderingBLS(DefaultLS())
sol_bd1u, sol_bd1p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end], xiu, xip)

linBdsolver = BK.MatrixFreeBLS(ls)
sol_bd2u, sol_bd2p, _, _ = linBdsolver(J11, J12, J21, J22, rhs[1:end-1], rhs[end], xiu, xip)

@test sol_bd1u ≈ sol_bd2u
@test sol_bd1p ≈ sol_bd2p
####################################################################################################
# test the linear solvers for matrix free formulations
J0 = I + sprand(100,100,0.1)
Jmf = x -> J0 * x
x0 = rand(100)
ls = DefaultLS()
out = ls(J0, x0)

# speciric linear solver built from a bordered linear solver
linSolver = BK.LSFromBLS()
sol_lin0 = linSolver(J0, x0)
@test sol_lin0[1] ≈ out[1]
sol_lin0 = linSolver(J0, x0, 2x0)
@test sol_lin0[1] ≈ out[1]
@test sol_lin0[2] ≈ 2out[1]

ls = GMRESKrylovKit(rtol = 1e-9, dim = 100)
outkk = ls(J0, x0)
@test out[1] ≈ outkk[1]
outkk = ls(Jmf, x0)
@test out[1] ≈ outkk[1]
outkk = ls(Jmf, x0; a₀ = 0., a₁ = 1.)
outkk = ls(Jmf, x0; a₀ = 0., a₁ = 1.5)
outkk = ls(Jmf, x0; a₀ = 1., a₁ = 1.)
outkk = ls(Jmf, x0; a₀ = 1., a₁ = 1.5)
outkk = ls(Jmf, x0; a₀ = 0.5, a₁ = 1.5)

# test preconditioner
Pl = lu(J0*0.9)
ls = GMRESKrylovKit(rtol = 1e-9, dim = 100, Pl = Pl)
outkk = ls(J0, x0)
@test out[1] ≈ outkk[1]
outkk = ls(Jmf, x0)
@test out[1] ≈ outkk[1]
outkk = ls(Jmf, x0; a₀ = 0.5, a₁ = 1.5)

ls = GMRESIterativeSolvers(N = 100, reltol = 1e-9)
outit = ls(J0, x0)
@test out[1] ≈ outit[1]
outkk = ls(J0, x0; a₀ = 0., a₁ = 1.)
outit = ls(J0, x0; a₀ = 0., a₁ = 1.5)
outit = ls(J0, x0; a₀ = 1., a₁ = 1.)
outit = ls(J0, x0; a₀ = 1., a₁ = 1.5)
outit = ls(J0, x0; a₀ = 0.5, a₁ = 1.5)

ls = GMRESIterativeSolvers(N = 100, reltol = 1e-9, ismutating = true)
Jom = (o,x) -> mul!(o,J0,x)
outit = ls(Jom, x0)
@test out[1] ≈ outit[1]
####################################################################################################
# test the shifted linear systems
rhs = rand(size(J0, 1))
sol0 = J0\rhs;

ls0 = GMRESIterativeSolvers(N = size(J0,1), reltol = 1e-10)
sol1, _ = ls0(J0, rhs)
@test norm(sol0 .- sol1, Inf) < 1e-8

h = 0.81
sol0 = (I - h.*J0)\rhs
sol1 = (I/h - J0)\rhs
@test norm(sol0 - sol1/h, Inf) < 1e-8

sol0,_ = ls0(I - h*J0, rhs)
sol1,_ = ls0(J0, rhs; a₀ = 1.0, a₁ = -h)
@test norm(sol0 - sol1, Inf) < 1e-8

ls0 = GMRESKrylovKit(atol = 1e-10)
sol0,_ = ls0(I - h*J0, rhs)
sol1,_ = ls0(J0, rhs; a₀ = 1.0, a₁ = -h)
@test norm(sol0 - sol1, Inf) < 1e-8

sol0,_ = ls0(I - h*J0, rhs)
sol1,_ = ls0(J0, rhs; a₀ = 1.0/h, a₁ = -1.)
@test norm(sol0 - sol1/h, Inf) < 1e-8

sol0,_ = ls0(I - h*J0, rhs)
sol1,_ = ls0(J0, rhs; a₀ = 1., a₁ = -h)
@test norm(sol0 - sol1, Inf) < 1e-8

####################################################################################################
# test the eigen solvers for matrix free formulations
# eil = BK.EigIterativeSolvers(tol = 1e-9)
out = Arpack.eigs(J0, nev = 20, which = :LR)

eil = BK.EigKrylovKit(tol = 1e-9)
outkk = eil(J0, 20)
geteigenvector(eil, outkk[2], 2)

eil = BK.EigKrylovKit(tol = 1e-9, x₀ = x0)
outkkmf = eil(Jmf, 20)
geteigenvector(eil, outkkmf[2], 2)

eil = BK.EigArpack(v0 = copy(x0))
outdefault = eil(J0, 20)
@test out[1] ≈ outdefault[1]
outdefault = eil(x ->J0*x, 20)
@test out[1] ≈ outdefault[1]

eil = BK.EigArnoldiMethod(;x₀ = x0)
outam = eil(J0, 20)
outam = eil(Jmf, 20)
geteigenvector(eil, outam[2], 2)

eil = BK.EigArnoldiMethod(;x₀ = x0, sigma = 1.)
outam = eil(J0, 20)
outam = eil(Jmf, 20)
geteigenvector(eil, outam[2], 2)


text:
# using Pkg
println("I am being built...")
# Pkg.add(PackageSpec(url="https://github.com/konkam/ExactWrightFisher.jl"))
# Pkg.add(PackageSpec(url="https://github.com/konkam/FeynmanKacParticleFilters.jl"))

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# ### Built In Types
# - [Integers and Floating Point Numbers](http://docs.julialang.org/en/latest/manual/integers-and-floating-point-numbers/)
# - [Strings and Characters](http://docs.julialang.org/en/latest/manual/strings/)
# - [Arrays](http://docs.julialang.org/en/latest/manual/arrays/)
# - [Dictionaries](http://docs.julialang.org/en/latest/stdlib/collections/#associative-collections)

1 + 2.0
typeof(3.3)

println("Hello World")
println("1 + 2 is: $(1 + 2)")

'a'
'\u2603'
int('☃')

['a' 'b' 'c']

['a' 'b' 'c';
'd' 'e' 'f';]

[1, 2, 3]

Dict(:erik => "python",
     :justin => "ruby",
     :colin => "scala",
)

# ### Variables
# - [Variables](http://docs.julialang.org/en/latest/manual/variables/)
# - [Variables and Scope](http://docs.julialang.org/en/latest/manual/variables-and-scoping/)
# - [Type Assertions](http://julia.readthedocs.org/en/latest/manual/types/#type-declarations)

a = 7
erik = :erik
erik::Symbol
erik::Array

# ### Control Flow
# - [Control Flow](http://docs.julialang.org/en/latest/manual/control-flow/)
# - [List Comprehensions](http://docs.julialang.org/en/latest/manual/arrays/#comprehensions)

a = 3
if a > 4
    println("a is too big")
else
    println("a is $a")
end

x = 1 == 2 ? "yeah!" : "nope..."

for x = 1:3:10
    println(x)
end

o_count = 0
for x = "one two three"
    if x == 'o'
        o_count += 1
    end
end
o_count

[x^2 for x = 1:5]

["$x $y" for x = 1:3, y='a':'c']

# ### Functions
# - [Functions](http://docs.julialang.org/en/latest/manual/functions/)

function minimum(l)
    min = l[1]
    for x = l[2:end]
        min = min < x ? min : x
    end
    min
end
minimum([2, -3, 6, -4, 3])

f(x) = 2x
f(3)

((x, y) -> x + y)(2, 3)

map(c -> c + 2, "Hello")

filter(x -> x % 3 == 0, 1:10)

max_stars_repo_name: Wilduck/IntroToJulia

max_issues_repo_name: swans-one/IntroToJulia

max_forks_repo_name: swans-one/IntroToJulia


text:
字符串="""
cat-kit-kitten
dog-doge-puppy
"""

emmm="iee"

println(字符串)
println(字符串[6:9])# 切片

println(emmm*"emmm")
println(emmm^6)

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
"""
    _to_union(datatype)
Make a minimal `Union` type from a collection of data types.
"""
_to_union(datatype) = Union{(isa(datatype, Type) ? [datatype] : datatype)...}

"""
    _find_rand_argmax(d::DictionaryView)
Compute `argmax` of `d` and select one element randomly.
"""
function _find_rand_argmax(d::DictionaryView)
    max = -Inf
    argmax = Vector{Int}()
    for (k, v) in pairs(d)
        if v > max
            max = v
            argmax = [k]
        elseif v == max
            push!(argmax, k)
        end
    end
    return rand(argmax)
end

abstract type FunctionContainer end
apply(fc::FC) where {FC <: FunctionContainer} = fc.f
apply(fc::FC, x, X) where {FC <: FunctionContainer} = convert(Float64, apply(fc)(x, X))
apply(fc::FC, x) where {FC <: FunctionContainer} = convert(Float64, apply(fc)(x))


text:
# This file was generated, do not modify it. # hide
ŷ = predict_mode(clf, X2)
misclassification_rate(ŷ, y)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
# This file provides glue from julia to our hnsw rust crate
#
#

using Printf

using Logging
using Base.CoreLogging

const libhnswso = "libhnsw_rs"

# to be called before anything
"""
# initialization of Julia DL_LOAD_PATH

 `function setRustlibPath(path::String)`

This function tells julia where is installed the rust dynamic library implementing the
Hnsw algorithm.
It must be called after `using Hnsw` and before any function call

The argument is the path to the rust library.
"""
function setRustlibPath(path::String)
    push!(Base.DL_LOAD_PATH, path)
end

"""
    A structure to encapsulate the Rust structure.
"""
mutable struct Hnswrs
end

logger = ConsoleLogger(stdout, CoreLogging.Debug)
global_logger(logger)

"""
# Struct Neighbour.

It contains the id of a neighbour and distance to the query point.

When searching for the neighbour of a given point, this struct is the basic block
of vector returned by search methods.
"""
struct Neighbour
    id :: UInt
    dist :: Float32
end

"""
# struct Neighbourhood.

A pointer to Neighbours
Structure returned by request searchNeighbour

"""
struct Neighbourhood
    nbgh :: Int64
    neighbours :: Ptr{Neighbour}
end

"""
# struct NeighbourhoodVect

 To retrive answer to parallel search of neigbourhood.
"""
struct NeighbourhoodVect
    nb_request :: UInt
    neighbourhoods :: Ptr{Neighbourhood}
end

"""
# implementedTypes

implementedTypes is a dictionary that has 2 functionalities
1.  It lists the types on which an instantiation is made for languages others than Rust.
2.  It enables the mapping to the correct method in rust interface by concatenating the type to the
    name of the method to be called.
"""
implementedTypes = Dict{DataType, String}()
implementedTypes[Float32] = "f32"
implementedTypes[UInt8] = "u8"
implementedTypes[UInt16] = "u16"
implementedTypes[Int32] = "i32"
implementedTypes[UInt32] = "u32"

function checkForImplementedType(type::DataType)
    if haskey(implementedTypes, type)
        return implementedTypes[type]
    else
        println("unimplemented type = ", type)
        throw("hnswrs : unimplement type")
    end
end

"""

# function hnswInit
`function hnswInit(type :: DataType, maxNbConn::Int64, efConstruction::Int64, distname::String)`

## Args
- type of data vector.
        The names of types are String and correspond to rust type names i.e
        f32, i32, u16, u8. So the type arg are "f32" , "i32" and so on.
        The subsequent request insertion or search must be made with data corresponding
        to the type used in initialization of Hnsw_api. The rust library will panic otherwise.

- maxNbConn. The maximum number of connection by node
- search parameter
- distname

 ## Return
    - A pointer to Hnswrs

"""
function hnswInit(type :: DataType, maxNbConn::Int64, efConstruction::Int64, distname::String)
    # check for type
    rust_type_name = checkForImplementedType(type)
    @eval hnsw = ccall(
            $(string("init_hnsw_", rust_type_name), libhnswso),
            Ptr{Hnswrs}, # return type
            (UInt64, UInt64, Int64, Ptr{UInt8},),
            UInt64($maxNbConn), UInt64($efConstruction), UInt64(length($distname)), pointer($distname)
        )
end

"""

# function hnswInit
`function hnswInit(type :: DataType, maxNbConn::Int64, efConstruction::Int64, f :: Ptr{Cvoid})`

## Args
- type of data vector.
    The names of types are String and correspond to rust type names i.e
    f32, i32, u16, u8. So the type arg are "f32" , "i32" and so on.
    The subsequent request insertion or search must be made with data corresponding
    to the type used in initialization of Hnsw_api. The rust library will panic otherwise.

- maxNbConn. The maximum number of connection by node
- search parameter
- a C-api function pointer to the distance to be used

 ## Return
    - A pointer to Hnswrs

"""
function hnswInit(type :: DataType, maxNbConn::Int64, efConstruction::Int64, f :: Ptr{Cvoid})
    @info "recieving function ptr : " f
    # check for type
    rust_type_name = checkForImplementedType(type)
    @eval hnsw = ccall(
        $(string("init_hnsw_ptrdist_", rust_type_name), libhnswso),
        Ptr{Hnswrs}, # return type
        (UInt64, UInt64, Ptr{Cvoid}),
        UInt64($maxNbConn), UInt64($efConstruction), $f
    )
end

###################  insert method

"""
# one_insert

` function one_insert(ptr::Ref{Hnswrs}, data::Vector{T}, id::UInt64) where {T <: Number} `

The function first checks that it is called for an implemented type.
It generates the name of the rust function to be called and
passes the call to @eval as we cannot call directly ccall with a
non constant couple (fname, library) Cf Julia manual
"""
function one_insert(ptr::Ref{Hnswrs}, data::Vector{T}, id::UInt64) where {T <: Number}
    rust_type_name = checkForImplementedType(eltype(data))
    @eval ccall(
        $(string("insert_", rust_type_name), libhnswso),
        Cvoid,
        (Ref{Hnswrs}, UInt, Ref{$T}, UInt64),
        $ptr, UInt(length($data)), $data, UInt64($id))
end

###################  parallel insert

"""
    # parallel insertion

` function parallel_insert(ptr::Ref{Hnswrs}, datas::Vector{Tuple{Vector{T}, UInt}}) where {T <: Number} `
"""
function parallel_insert(ptr::Ref{Hnswrs}, datas::Vector{Tuple{Vector{T}, UInt}}) where {T <: Number}
    # get data type of first field of tuple in datas
    d_type = eltype(fieldtype(eltype(datas),1))
    rust_type_name = checkForImplementedType(d_type)
    # split vector of tuple
    nb_vec = length(datas)
    # get length of first component of first data (i.e dim of first vector)
    dim = length(datas[1][1])
    # make a Vector{Ref{Float32}} where each ptr is a ref to datas[i] memory beginning
    vec_ref = map(x-> pointer(x[1]), datas)
    ids_ref = map(x-> x[2], datas)
    @eval neighbourhood_vec_ptr = ccall(
        $(string("parallel_insert_", rust_type_name), libhnswso),
        Cvoid,
        (Ref{Hnswrs}, UInt, UInt, Ref{Ptr{$T}}, Ref{UInt}),
        $ptr, UInt($nb_vec), UInt($dim), $vec_ref, $ids_ref
    )
end

########   search method

"""
# function one_search

 `  function one_search(ptr::Ref{Hnswrs}, vector::Vector{T}, knbn::Int64, ef_search ::Int64) where {T<:Number}`
"""
function one_search(ptr::Ref{Hnswrs}, vector::Vector{T}, knbn::Int64, ef_search ::Int64) where {T<:Number}
    rust_type_name = checkForImplementedType(eltype(vector))
    #
    @eval neighbours_ptr = ccall(
        $(string("search_neighbours_", rust_type_name), libhnswso),
        Ptr{Neighbourhood},
        (Ref{Hnswrs}, UInt, Ref{$T}, UInt, UInt),
        $ptr, UInt(length($vector)), $vector, UInt($knbn), UInt($ef_search)
    )
    # now return a Vector{Neighbourhood}
    # @debug "\n search (rust) returned pointer, will do unsafe_load"  neighbours_ptr
    neighbourhood = unsafe_load(neighbours_ptr::Ptr{Neighbourhood})
    # @debug "\n search rs returned neighbours "  neighbourhood
    neighbours = unsafe_wrap(Array{Neighbour,1}, neighbourhood.neighbours, NTuple{1,Int64}(neighbourhood.nbgh); own = true)
    # @debug "neighbours : " neighbours
    # we got Vector{Neighbour}
    return neighbours
end

# we must return a Vector{Vector{Neighbour}} , one Vector{Neighbour} per request input
"""
# parallel_search function

` function parallel_search(ptr::Ref{Hnswrs}, datas::Vector{Vector{T}}, knbn::Int64, ef_search:: Int64) where {T<:Number} `

parallel search of a Vector of Vectors with search parameters.
"""
function parallel_search(ptr::Ref{Hnswrs}, datas::Vector{Vector{T}}, knbn::Int64, ef_search:: Int64) where {T<:Number}
    d_type = eltype(eltype(datas))
    rust_type_name = checkForImplementedType(d_type)
    #
    nb_vec = length(datas)
    len = length(datas[1])
    # make a Vector{Ref{Float32}} where each ptr is a ref to datas[i] memory beginning
    vec_ref = map(x-> pointer(x), datas)
    @eval neighbourhood_vec_ptr = ccall(
        $(string("parallel_search_neighbours_", rust_type_name), libhnswso),
        Ptr{NeighbourhoodVect},
        (Ref{Hnswrs}, UInt, UInt, Ref{Ptr{$T}}, UInt, UInt),
        $ptr, UInt($nb_vec), UInt($len), $vec_ref, UInt($knbn), UInt($ef_search)
    )
    # @debug "\n parallel_search_neighbours rust returned pointer" neighbourhood_vec_ptr
    neighbourhoods_vec = unsafe_load(neighbourhood_vec_ptr::Ptr{NeighbourhoodVect})
    neighbourhoods = unsafe_wrap(Array{Neighbourhood,1}, neighbourhoods_vec.neighbourhoods, NTuple{1,Int64}(neighbourhoods_vec.nb_request); own = true)
    neighbourhoods_answer = Vector{Vector{Neighbour}}(undef, nb_vec)
    # now we must unwrap each neighbourhood
    for i in 1:nb_vec
        # @debug "\n unwraping neighbourhood of request " i
        neighbourhoods_answer[i] = unsafe_wrap(Array{Neighbour,1}, neighbourhoods[i].neighbours, NTuple{1,Int64}(neighbourhoods[i].nbgh); own = true)
    end
    #
    return neighbourhoods_answer
end

function filedump(ptr::Ref{Hnswrs}, d_type :: DataType, filename::String)
    rust_type_name = checkForImplementedType(d_type)
    resdump = @eval ccall(
        $(string("file_dump_", rust_type_name), libhnswso),
        Clonglong,
        (Ref{Hnswrs}, UInt64, Ptr{UInt8},),
        $ptr, UInt64(length($filename)), pointer($filename)
    )
    resdump
end

struct LoadHnswDescription
    dumpmode :: UInt8
    #max number of connections in layers != 0
    max_nb_connection :: UInt8
    # number of observed layers
    nb_layer :: UInt8
    # search parameter
    ef :: UInt64
    # total number of points
    nb_point:: UInt64
    # dimesion of data vector
    data_dimension :: UInt64
    # length and pointer on dist name
    distname_len :: UInt64
    distname :: Ptr{UInt8}
    # T typename
    t_name_len :: UInt64
    t_name :: Ptr{UInt8}
end

"""
# struct HnswDescription

When trying to reload a Hnsw structure from a previous dump
it is necessary to known some characteristics of the data dumped.
So first a call to getDescription is made, then wth info returnded
it is possible to call loadHnsw with the adequate parameters.
The parameters necessay to call loadHnsw are:

-   DataType of vectors stored in structure (it also help not inserting
        different size of data)
-   The distance name (just now the case of Ptr distance is not treated, The name of distance
        is set to DistPtr)
    The others descriptors are just given for information but not necessary.
"""
struct HnswDescription
    maxNbConn::Int64
    # number of observed layers
    nb_layer :: Int64
    # search parameter
    ef :: Int64
    # total number of points
    nb_point:: Int64
    # dimesion of data vector
    data_dimension :: Int64
    # type of vector
    type :: DataType
    # name of distance
    distname::String
    # pointer on distance function
    distfunctPtr :: Union{Some{Ptr{Cvoid}} , Nothing}
end

#  to load a description of graph stored in a file
#  filename is the base of the filename (without suffixes "hnsw.graph" of hnsw.data"

"""
# `function getDescription(filename::String `

This function returns a description of graph stored such as type of data stored (Float32,....)
    type of distance and search parameters.
    Information in data type is necessary to be able to instantiate the rust library and is used
    by function loadHnsw(filename :: String, ....)
"""
function getDescription(filename :: String)
    #
    description_ptr = ccall(
        (:load_hnsw_description, libhnswso),
        Ptr{LoadHnswDescription},
        (UInt64, Ptr{UInt8},),
        UInt64(length(filename)), pointer(filename)
    )
    #
    if description_ptr == C_NULL
        println("call to getDescription failed , filename was :", filename)
        return nothing
    end
    #
    ffiDescription  = unsafe_load(description_ptr::Ptr{LoadHnswDescription})
    @info " data dimension : " ffiDescription.data_dimension
    typename_u = unsafe_wrap(Array{UInt8,1}, ffiDescription.t_name, NTuple{1,UInt64}(ffiDescription.t_name_len); own = true)
    typename = String(typename_u)
    @info "getDescription got typename : " typename
    # get key for typename
    allkeys = collect(keys(implementedTypes))
    keyindex = findfirst(x-> implementedTypes[x] == typename , allkeys)
    if keyindex === nothing
        # this should not occur
        println("type not implemented : ", typename)
        return nothing
    end
    keytype = allkeys[keyindex]
    @info " keytype : " keytype
    # now we have rust type name and we know it is implemented
    # we must check if distname is "DistPtr"
    distname_u = unsafe_wrap(Array{UInt8,1}, ffiDescription.distname, NTuple{1,UInt64}(ffiDescription.distname_len); own = true)
    distname = String(distname_u)
    # dump info
    println("Description :")
    println("distance name : ", distname)
    println("data type : ", typename)
    #
    HnswDescription(Int64(ffiDescription.max_nb_connection),
                    Int64(ffiDescription.nb_layer),
                    Int64(ffiDescription.ef),
                    Int64(ffiDescription.nb_point),
                    Int64(ffiDescription.data_dimension),
                    keytype,
                    distname,
                    nothing
    )
end

"""
# `function loadHnsw(filename :: String, type :: DataType, distname :: String)`

This function reloads data from the 2 files created when dumping data
the graph and data files.
The filename sent as arg is the base of the names used to dump files in.
It does not have the suffixes ".hnsw.graph" and ".hnsw.data"

This function returns a couple (HnswDescription, )
"""
function loadHnsw(filename :: String, type :: DataType, distname :: String)
    # append hnsw.graph and load description
    graphFileName  = filename*".hnsw.graph"
    # get a HnswDescription
    description = getDescription(graphFileName)
    # check for types
    if description === nothing
        @warn "some error occurred, could not load a coherent description"
        return nothing
    end
    println("dimension of data %d", description.data_dimension)
    rust_type_name = checkForImplementedType(type)
    # call rust stub
    @eval hnsw = ccall(
            $(string("load_hnswdump_", rust_type_name), libhnswso),
            Ptr{Hnswrs}, # return type
            (UInt64, Ptr{UInt8}),
            UInt64(length($filename)), pointer($filename)
        )
    #
    maxNbConn = description.maxNbConn
    # we do not know how rust constructed it...
    efConstruction = 0
    distname_load = description.distname
    # coherence check
    findres = findfirst(distname, distname_load)
    if findres == nothing
        @warn "some error occurred, distances do not match, expected %s, got %s", distname, distname_load
        return nothing
    end
    hnswapi = HnswApi(hnsw, type, maxNbConn, efConstruction, distname, nothing)
    (description, hnswapi)
end


text:
using GLMakie, GeometryBasics, RPRMakie, RadeonProRender
using Colors
using Colors: N0f8
RPR = RadeonProRender

f = (u, v) -> cos(v)*(6 - (5/4 +sin(3*u))*sin(u-3*v))
g = (u, v) -> sin(v)*(6 - (5/4 +sin(3*u))*sin(u-3*v))
h = (u, v) -> -cos(u-3*v) * (5/4 +sin(3*u));
u = range(0, stop=2π, length=150)
v = range(0, stop=2π, length=150)

context = RPR.Context()
matsys = RPR.MaterialSystem(context, 0)
material = RPR.UberMaterial(matsys)

# for (k, v) in pairs(RPR.defaults(RPR.UberMaterial))
#     setproperty!(material, k, v)
# end

fig, ax, pl = surface(f.(u,v'),
             g.(u,v'),
             h.(u,v'),
             ambient=Vec3f(0.5),
             diffuse=Vec3f(1),
             specular=0.5,
             colormap=:balance, figure=(resolution=(1500, 1000),),
             show_axis=false, material=material)

function Input(fig, val::RGB)
    hue = Slider(fig, range = 1:380, width=200)
    lightness = Slider(fig, range = LinRange(0, 1, 100), width=200)
    labels = [Label(fig, "hue"; halign = :left), Label(fig, "light"; halign = :left)]
    layout = grid!(hcat(labels, [hue, lightness]))
    hsl = HSL(val)
    set_close_to!(hue, hsl.h)
    set_close_to!(lightness, hsl.l)
    color = map((h, l)-> RGB(HSL(h, 0.9, l)), hue.value, lightness.value)
    return color, layout
end

function Input(fig, val::Vec4)
    s = Slider(fig, range = LinRange(0, 1, 100), width=200)
    set_close_to!(s, first(val))
    return map(x-> Vec4f(x), s.value), s
end

function Input(fig, val::Bool)
    toggle = Toggle(fig, active = val)
    return toggle.active, toggle
end

sliders = (
    reflection_color = Input(fig, RGB(0, 0, 0)),
    reflection_weight = Input(fig, Vec4(0)),
    reflection_roughness = Input(fig, Vec4(0)),
    reflection_anisotropy = Input(fig, Vec4(0)),
    reflection_anisotropy_rotation = Input(fig, Vec4(0)),
    reflection_mode = Input(fig, Vec4(0)),
    reflection_ior = Input(fig, Vec4(0)),
    reflection_metalness = Input(fig, Vec4(0)),

    refraction_color = Input(fig, RGB(0, 0, 0)),
    refraction_weight = Input(fig, Vec4(0)),
    refraction_roughness = Input(fig, Vec4(0)),
    refraction_ior = Input(fig, Vec4(0)),
    refraction_absorption_color = Input(fig, RGB(0, 0, 0)),
    refraction_absorption_distance = Input(fig, Vec4(0)),
    refraction_caustics = Input(fig, true),

    sss_scatter_color = Input(fig, RGB(0, 0, 0)),
    sss_scatter_distance = Input(fig, Vec4(0)),
    sss_scatter_direction = Input(fig, Vec4(0)),
    sss_weight = Input(fig, Vec4(0)),
    sss_multiscatter = Input(fig, false),
    backscatter_weight = Input(fig, Vec4(0)),
    backscatter_color = Input(fig, RGB(0, 0, 0)),
)

labels = []
inputs = []
refresh = Observable(nothing)
for (key, (obs, input)) in pairs(sliders)
    push!(labels, Label(fig, string(key); align=:left))
    push!(inputs, input)
    on(obs) do value
        setproperty!(material, key, value)
        notify(refresh)
    end
end

fig[1, 2] = grid!(hcat(labels, inputs), width=500)
GLMakie.inline!(false)
display(fig)
context, task = RPRMakie.replace_scene_rpr!(ax.scene, context, matsys; refresh=refresh)


text:
using Flux
using Random

mutable struct ImageDQNAgent_old{M, TN, O, LU, AP<:AbstractValuePolicy, Φ, ER<:AbstractImageReplay} <: AbstractAgent
    model::M
    target_network::TN
    opt::O
    lu::LU
    ap::AP
    er::ER
    batch_size::Int
    tn_counter_init::Int
    target_network_counter::Int
    wait_time::Int
    wait_time_counter::Int
    exp_wait_size::Int
    action::Int
    prev_s::Φ
end

ImageDQNAgent_old(model,
              target_network,
              image_replay,
              opt,
              lu,
              ap,
              size_buffer,
              batch_size,
              tn_counter_init,
              wait_time,
              exp_wait_size) =
    ImageDQNAgent_old(model,
                  target_network,
                  opt,
                  lu,
                  ap,
                  image_replay,
                  batch_size,
                  tn_counter_init,
                  tn_counter_init,
                  wait_time,
                  0,
                  exp_wait_size,
                  0,
                  zeros(Int, image_replay.hist))

function get_action(agent::ImageDQNAgent_old, s)

end

function _start!(agent::ImageDQNAgent_old,
                       env_s_tp1,
                       rng::AbstractRNG;)

    # Start an Episode
    agent.prev_s .= add!(agent.er, env_s_tp1)

    prev_s =
        cat(getindex(agent.er.image_buffer, agent.prev_s)./256f0;
            dims=4) |> gpu

    agent.action = sample(agent.ap,
                          cpu(agent.model(prev_s)),
                          rng)

    return agent.action
end

function _step!(agent::ImageDQNAgent_old,
                      env_s_tp1,
                      r,
                      terminal,
                      rng::AbstractRNG;)

    cur_s = add!(agent.er,
                 (agent.prev_s,
                  findfirst((a)->a==agent.action,
                            agent.ap.action_set),
                  env_s_tp1,
                  r,
                  terminal))

    update_params!(agent, rng)

    agent.prev_s .= cur_s

    prev_s =
        cat(getindex(agent.er.image_buffer, agent.prev_s)./256f0;
            dims=4) |> gpu

    agent.action = sample(agent.ap,
                          cpu(agent.model(prev_s)),
                          rng)

    return agent.action
end

function update_params!(agent::ImageDQNAgent_old, rng)

    if size(agent.er)[1] > agent.exp_wait_size
        agent.wait_time_counter -= 1
        if agent.wait_time_counter <= 0

            e = sample(agent.er,
                       agent.batch_size;
                       rng=rng)
            s = gpu(e.s)
            r = gpu(e.r)
            t = gpu(e.t)
            sp = gpu(e.sp)

            update!(agent.model,
                    agent.lu,
                    agent.opt,
                    s,
                    e.a,
                    sp,
                    r,
                    t,
                    agent.target_network)

            agent.wait_time_counter = agent.wait_time
        end
    end

    # Target network updates
    if !(agent.target_network isa Nothing)
        if agent.target_network_counter == 1
            agent.target_network_counter = agent.tn_counter_init
            for ps ∈ zip(params(agent.model),
                         params(agent.target_network))
                ps[2] .= ps[1]
            end
        else
            agent.target_network_counter -= 1
        end
    end

    return nothing
end


text:
using Base.Test
using Quaternions

qx = qrotation([1,0,0], pi/4)
@test_approx_eq qx*qx qrotation([1,0,0], pi/2)
@test_approx_eq qx^2 qrotation([1,0,0], pi/2)
theta = pi/8
qx = qrotation([1,0,0], theta)
c = cos(theta); s = sin(theta)
Rx = [1 0 0; 0 c -s; 0 s c]
@test_approx_eq rotationmatrix(qx) Rx
theta = pi/6
qy = qrotation([0,1,0], theta)
c = cos(theta); s = sin(theta)
Ry = [c 0 s; 0 1 0; -s 0 c]
@test_approx_eq rotationmatrix(qy) Ry
theta = 4pi/3
qz = qrotation([0,0,1], theta)
c = cos(theta); s = sin(theta)
Rz = [c -s 0; s c 0; 0 0 1]
@test_approx_eq rotationmatrix(qz) Rz

@test_approx_eq rotationmatrix(qx*qy*qz) Rx*Ry*Rz
@test_approx_eq rotationmatrix(qy*qx*qz) Ry*Rx*Rz
@test_approx_eq rotationmatrix(qz*qx*qy) Rz*Rx*Ry

a, b = qrotation([0,0,1], deg2rad(0)), qrotation([0,0,1], deg2rad(180))
@test_approx_eq slerp(a,b,0.0) a
@test_approx_eq slerp(a,b,1.0) b
@test_approx_eq slerp(a,b,0.5) qrotation([0,0,1], deg2rad(90))

@test_approx_eq angle(qrotation([1,0,0], 0)) 0
@test_approx_eq angle(qrotation([0,1,0], pi/4)) pi/4
@test_approx_eq angle(qrotation([0,0,1], pi/2)) pi/2


text:
using Compat.LinearAlgebra: chkstride1, checksquare, LAPACK.chklapackerror,
    BlasInt, BLAS.@blasfunc

const liblapack = Base.liblapack_name

hermfact!(w::Vector,Z::Matrix,H::Hermitian) = syevr!(w,Z,H.uplo,H.data,0.0,0.0,0,0,-1.0)

for (syevr1, syevr2, elty, relty) in
    ((:dsyevr_,:zheevr_,:ComplexF64,:Float64),
    (:ssyevr_,:cheevr_,:ComplexF32,:Float32))
    @static if VERSION < v"0.7.0-"
        @eval begin
            # Real symmetric
            function syevr!(w::Vector{$relty}, Z::Matrix{$relty}, uplo::Char, A::StridedMatrix{$relty},
                vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
                chkstride1(A)
                n = checksquare(A)
                length(w) == n || throw(ArgumentError("vector w should have length $n."))
                size(Z) == (n,n) || throw(ArgumentError("matrix Z should have size $n by $n."))
                jobz = 'V'
                range = 'A'
                lda = stride(A,2)
                m = Ref{Int}()
                ldz = n
                isuppz = Vector{Int}(2n)
                work   = Vector{$relty}(1)
                lwork  = BlasInt(-1)
                iwork  = Vector{BlasInt}(1)
                liwork = BlasInt(-1)
                info   = Ref{BlasInt}()
                for i = 1:2
                    ccall((@blasfunc($syevr1), liblapack), Void,
                    (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty},
                    Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
                    Ptr{BlasInt}),
                    &jobz, &range, &uplo, &n,
                    A, &max(1,lda), &vl, &vu,
                    &il, &iu, &abstol, m,
                    w, Z, &max(1,ldz), isuppz,
                    work, &lwork, iwork, &liwork,
                    info)
                    chklapackerror(info[])
                    if i == 1
                        lwork = BlasInt(real(work[1]))
                        work = Vector{$relty}(lwork)
                        liwork = iwork[1]
                        iwork = Vector{BlasInt}(liwork)
                    end
                end
                return w, Z
            end
            # Hermitian
            function syevr!(w::Vector{$relty}, Z::Matrix{$elty}, uplo::Char, A::StridedMatrix{$elty},
                vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
                chkstride1(A)
                n = checksquare(A)
                length(w) == n || throw(ArgumentError("vector w should have length $n."))
                size(Z) == (n,n) || throw(ArgumentError("matrix Z should have size $n by $n."))
                jobz = 'V'
                range = 'A'
                lda = max(1,stride(A,2))
                m = Ref{Int}()
                ldz = n
                isuppz = Vector{Int}(2n)
                work   = Vector{$elty}(1)
                lwork  = BlasInt(-1)
                rwork  = Vector{$relty}(1)
                lrwork = BlasInt(-1)
                iwork  = Vector{BlasInt}(1)
                liwork = BlasInt(-1)
                info   = Ref{BlasInt}()
                for i = 1:2
                    ccall((@blasfunc($syevr2), liblapack), Void,
                    (Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt},
                    Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
                    Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt},
                    Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt},
                    Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                    &jobz, &range, &uplo, &n,
                    A, &lda, &vl, &vu,
                    &il, &iu, &abstol, m,
                    w, Z, &ldz, isuppz,
                    work, &lwork, rwork, &lrwork,
                    iwork, &liwork, info)
                    chklapackerror(info[])
                    if i == 1
                        lwork = BlasInt(real(work[1]))
                        work = Vector{$elty}(lwork)
                        lrwork = BlasInt(rwork[1])
                        rwork = Vector{$relty}(lrwork)
                        liwork = iwork[1]
                        iwork = Vector{BlasInt}(liwork)
                    end
                end
                return w, Z
            end
        end
    else
        @eval begin
            # Real symmetric
            function syevr!(w::Vector{$relty}, Z::Matrix{$relty}, uplo::AbstractChar, A::AbstractMatrix{$relty},
                vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
                chkstride1(A)
                n = checksquare(A)
                length(w) == n || throw(ArgumentError("vector w should have length $n."))
                size(Z) == (n,n) || throw(ArgumentError("matrix Z should have size $n by $n."))
                jobz = 'V'
                range = 'A'
                lda = stride(A,2)
                m = Ref{BlasInt}()
                ldz = n
                isuppz = similar(A, BlasInt, 2*n)
                work   = Vector{$relty}(undef, 1)
                lwork  = BlasInt(-1)
                iwork  = Vector{BlasInt}(undef, 1)
                liwork = BlasInt(-1)
                info   = Ref{BlasInt}()
                for i = 1:2  # first call returns lwork as work[1] and liwork as iwork[1]
                    ccall((@blasfunc($syevr1), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                    Ptr{$relty}, Ref{BlasInt}, Ref{$relty}, Ref{$relty},
                    Ref{BlasInt}, Ref{BlasInt}, Ref{$relty}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{$relty}, Ref{BlasInt}, Ptr{BlasInt},
                    Ptr{$relty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{BlasInt}),
                    jobz, range, uplo, n,
                    A, max(1,lda), vl, vu,
                    il, iu, abstol, m,
                    w, Z, max(1,ldz), isuppz,
                    work, lwork, iwork, liwork,
                    info)
                    chklapackerror(info[])
                    if i == 1
                        lwork = BlasInt(real(work[1]))
                        resize!(work, lwork)
                        liwork = iwork[1]
                        resize!(iwork, liwork)
                    end
                end
                return w, Z
            end
            # Hermitian
            function syevr!(w::Vector{$relty}, Z::Matrix{$elty}, uplo::AbstractChar, A::AbstractMatrix{$elty},
                vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat)
                chkstride1(A)
                n = checksquare(A)
                length(w) == n || throw(ArgumentError("vector w should have length $n."))
                size(Z) == (n,n) || throw(ArgumentError("matrix Z should have size $n by $n."))
                jobz = 'V'
                range = 'A'
                lda = max(1,stride(A,2))
                m = Ref{BlasInt}()
                ldz = n
                isuppz = similar(A, BlasInt, 2*n)
                work   = Vector{$elty}(undef, 1)
                lwork  = BlasInt(-1)
                rwork  = Vector{$relty}(undef, 1)
                lrwork = BlasInt(-1)
                iwork  = Vector{BlasInt}(undef, 1)
                liwork = BlasInt(-1)
                info   = Ref{BlasInt}()
                for i = 1:2  # first call returns lwork as work[1], lrwork as rwork[1] and liwork as iwork[1]
                    ccall((@blasfunc($syevr2), liblapack), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ref{$elty},
                    Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{BlasInt},
                    Ptr{$relty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{$relty}, Ref{BlasInt},
                    Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}),
                    jobz, range, uplo, n,
                    A, lda, vl, vu,
                    il, iu, abstol, m,
                    w, Z, ldz, isuppz,
                    work, lwork, rwork, lrwork,
                    iwork, liwork, info)
                    chklapackerror(info[])
                    if i == 1
                        lwork = BlasInt(real(work[1]))
                        resize!(work, lwork)
                        lrwork = BlasInt(rwork[1])
                        resize!(rwork, lrwork)
                        liwork = iwork[1]
                        resize!(iwork, liwork)
                    end
                end
                return w, Z
            end
        end
    end
end


text:
module WorldVocoder

using WORLD_jll: libworld

const DEFAULT_FRAME_PERIOD = 5.0

include("dio.jl")
export dio

include("harvest.jl")
export harvest

include("cheaptrick.jl")
export cheaptrick

include("stonemask.jl")
export stonemask

include("d4c.jl")
export d4c

include("synthesis.jl")
export synthesis

export world

struct Harvest end
struct Dio end

"""
    world(x::Vector{Float64}, fs::Int)

Returns `(f0 = refined_f0, sp = spectrogram, ap = aperiodicity)` after running the full WORLD pipeline
consisting of

- `dio` (f0 extraction)
- `stonemask` (f0 refinement)
- `cheaptrick` (spectrogram extraction)
- `d4c` (aperiodicity extraction)
"""
function world(x::Vector{Float64}, fs::Int, algorithm = Harvest())
    f0, timestamps = extract_f0(x, fs, algorithm)
    spectrogram = cheaptrick(x, fs, timestamps, f0)
    fft_size = (size(spectrogram, 1) - 1) * 2
    aperiodicity = d4c(x, fs, timestamps, f0, fft_size)
    (f0 = f0, sp = spectrogram, ap = aperiodicity)
end

function world(x::AbstractVector{<:Real}, fs, algorithm = Harvest())
    world(convert(Vector{Float64}, x), convert(Int, fs), algorithm)
end

function world(x::AbstractMatrix{<:Real}, fs, algorithm = Harvest())
    if any(==(1), size(x))
        x = vec(x)
    else
        error("Only signals that are vector-like (one dimension length 1) are allowed.")
    end
    world(convert(Vector{Float64}, x), convert(Int, fs), algorithm)
end

function extract_f0(x, fs, ::Harvest)
    refined_f0, timestamps = harvest(x, fs)
end

function extract_f0(x, fs, ::Dio)
    f0, timestamps = dio(x, fs)
    refined_f0 = stonemask(x, fs, timestamps, f0)
    refined_f0, timestamps
end
end


text:
using Test
using MultivariateMoments

@testset "SymMatrix" begin
    Q = SymMatrix([1, 2, 3], 2)
    @test Q[1:2, 1:2] isa Matrix{Int}
    @test Q[1:2, 1:2] == [1 2; 2 3]
    @test square_getindex(Q, 1:2) isa SymMatrix{Int}
    @test square_getindex(Q, 1:2).Q == [1, 2, 3]
    @test square_getindex(Q, 2:-1:1) isa SymMatrix{Int}
    @test square_getindex(Q, 2:-1:1).Q == [3, 2, 1]
    @test square_getindex(Q, 1:1) isa SymMatrix{Int}
    @test square_getindex(Q, 1:1).Q == [1]
    @test square_getindex(Q, 2:2) isa SymMatrix{Int}
    @test square_getindex(Q, 2:2).Q == [3]
    @test similar(Q, Float64, (Base.OneTo(2), Base.OneTo(2))) isa Matrix{Float64}
    @test eltype(Q) == Int
    symmetric_setindex!(Q, 4, 1, 1)
    @test Q.Q == [4, 2, 3]
    symmetric_setindex!(Q, 5, 1, 2)
    @test Q.Q == [4, 5, 3]
    P = copy(Q)
    @test P.n == 2
    symmetric_setindex!(P, 6, 2, 2)
    @test Q.Q == [4, 5, 3]
    @test P.n == 2
    @test P.Q == [4, 5, 6]
    R = map(i -> i - 1, P)
    @test R.n == 2
    @test R.Q == [3, 4, 5]
    for S in [similar(Q, 1, 1), similar(Q, (1, 1)),
              similar(typeof(Q), 1, 1), similar(typeof(Q), (1, 1))]
        @test typeof(S) == typeof(Q)
        symmetric_setindex!(S, 2, 1, 1)
        @test S[1, 1] == 2
    end
end