So, we have P'*G**P = D as usual.
Let L be the lattice ZZ^n and L_P the matrix spanned by the columns of P.
Since P is invertible, there exists some big enough integer M so that M*L is a subset of L_P, right?
Since the elements P are orthogonal, it's way easier to reason about this lattice I think, and I'm wondering if it's a good idea to do the computations on L_P and translate back as necessary… I guess the problem is that there will be many points of L_P that, when divided by M, aren't integer: superfluous points.

I'm not sure what I'm trying to say, but I'd like to write the routine that generates all roots now, before doing the fundamental cone computation, and I'm quite confused by how it all goes: The thesis of Guglielmetti does things differently than B&P, and they have a hard to read implementation which isn't very well synchronized with their paper

More questions: (no need to answer quickly, I'm writing them down for the record)

1. Say the transformation P has columns v0,…,v_n with v0 of negative norm.
   Then in B&P, they talk about the orthogonal complement of v0: if we were doing vector spaces, we could juste take v1,…,v_n as a basis and be happy, but in our case v1,…,v_n now span a sub-lattice of the orthogonal complement of v0:
   Is there a known process to modify those vectors to get a basis for the complement itself?
   In the notation of B&P, this would be a basis of the space L'.
2. Another not-so-related one: When they build the fundamental cone P_0, they test for each candidate root whether the hyperplane is actually necessary by testing if it intersects the previous cone non-trivially.
   To do so, they call some sage functions which do the check for them.
   Is it obvious how to do this check by hand?
3. What's a good reference book/resource talking about what I need to know for these lattice/quad forms business?

More output should be human readable: e.g. one should be able to extract the diagram / Coxeter matrix from Vinberg_Algorithm. It should return some kind of an encapsulated objects that has:

1. the quad form we started with
2. v0
3. fund cone
4. the whole set of roots (if finite volume)
5. the diagram in some format that is suitable later on for drawing, e.g. by applying Draw_Diagram(x), where x = Vinberg_Algorithm(y), y = matrix of the quad form.

I tried to generate random roots, but that seems rather inefficient actually, since I guess a "random" vector won't be a root.

I would like to try a strategy consisting in generating roots in increasing distance, but with no care for whether the closest ones to v₀ are chosen.
The idea would be that, perhaps, random roots will generally be sufficient to cover a finite polyhedron if such exist.
Maybe that could even be made into a formal probabilistic argument: What's the probability that a set of n roots define a finite polyhedron, if the lattice is reflective?

Opinions, @sashakolpakov ?

Probably in drop_redundant_hyperplanes one could warmstart the solver so that any insight gained at last iteration could be reused.

Everything is in the title

More precisely, we start with the integral lattice Z^d, and an integral bilinear form on it, defined by a symmetric integral matrix A.
The diagonalization process yields an invertible matrix P and a diagonal matrix D satisfying P'AP = D with P' the transpose of P (all have integral coefficients).
Now, if we were talking about vector spaces (say over R), it would be clear that the columns of P can be viewed as a basis of R^d because P is invertible.
If I'm recalling correctly, a set of d vectors in Z^d spans Z^d as a lattice iff the determinant of the matrix they define is 1.
So the question is, is the determinant always ±1.

On include("src/vinberg.jl") always gives "WARNING: redefinition of constant deg_seq_a1. This may fail, cause incorrect answers, or produce other errors."

WARNING: using AbstractAlgebra.zero in module Main conflicts with an existing identifier.

This input

[ 4 -2 0 0 0 0 0 0 0 0 0 0 0 0 0;
-2 4 -2 2 0 0 0 0 0 0 0 0 -1 0 0;
0 -2 4 0 0 2 0 0 0 0 0 0 2 0 0;
0 2 0 4 2 2 0 0 0 0 0 0 0 0 0;
0 0 0 2 4 2 0 0 2 1 0 0 0 0 0;
0 0 2 2 2 4 2 2 1 2 0 0 1 0 0;
0 0 0 0 0 2 4 2 0 2 0 0 0 0 0;
0 0 0 0 0 2 2 4 0 2 0 0 1 0 0;
0 0 0 0 2 1 0 0 4 2 0 0 0 0 0;
0 0 0 0 1 2 2 2 2 4 2 2 1 0 0;
0 0 0 0 0 0 0 0 0 2 4 2 2 0 0;
0 0 0 0 0 0 0 0 0 2 2 4 1 0 0;
0 -1 2 0 0 1 0 1 0 1 2 1 4 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 2;
0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 ]

produces the following output

[ Info: Initialized root iterator.
[ Info: Got all roots at distance zero.
[ Info: Starting with 28 candidate roots (at distance zero) for fundamental cone.
[ Info: Returning 14 cone roots.
[ Info: Got the roots of a fundamental cone.
[ Info: : [-1, -2, -4, 0, -2, 4, -1, -2, 0, 0, -1, 0, 2, 0, 0].
[ Info: : [-1, -2, -3, 0, -2, 4, -1, -2, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [-1, -2, -3, 0, -2, 4, 0, 0, 2, -4, 1, 2, 0, 0, 0].
[ Info: : [-1, -2, -2, 1, -2, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [-1, -2, -2, 1, -2, 2, 0, 0, 2, -2, 1, 0, 0, 0, 0].
[ Info: : [-1, -2, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [0, 0, -1, -1, 0, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [0, 0, -1, -1, 0, 2, 0, 0, 0, -2, 1, 0, 0, 0, 0].
[ Info: : [0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0].
[ Info: : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0].
[ Info: : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1].
[ Info: Built the corresponding Coxeter diagram.
[ Info: Found new satisfying root: [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [0, 0, 1, 1, 0, -2, 1, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [0, 0, 1, 1, 0, -2, 0, 0, 0, 2, -1, 0, 0, 1, 0].
[ Info: Found new satisfying root: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0].
[ Info: Found new satisfying root: [1, 2, 4, 0, 2, -4, 1, 2, 0, 0, 1, 0, -2, 1, 0].
[ Info: Found new satisfying root: [1, 2, 3, 0, 2, -4, 0, 0, -2, 4, -1, -2, 0, 1, 0].
[ Info: Found new satisfying root: [1, 2, 2, -1, 2, -2, 0, 0, -2, 2, -1, 0, 0, 1, 0].
[ Info: Found new satisfying root: [1, 2, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0].
[ Info: Found new satisfying root: [1, 2, 2, -1, 2, -2, 1, 0, 0, 0, 0, 0, 0, 1, 0].

and then the process breaks off.

Could squeeze some performance.

Maybe this can be used to compute in a smaller lattice, and then extend somehow?

I used this lattice: (call it BP_D4D4-1.lat).

2 0 1 0 0 0 0 0 0
0 2 -1 0 0 0 0 0 0
1 -1 2 -1 0 0 0 0 0
0 0 -1 2 0 0 0 0 0
0 0 0 0 2 0 1 0 0
0 0 0 0 0 2 -1 0 0
0 0 0 0 1 -1 2 -1 0
0 0 0 0 0 0 -1 2 0
0 0 0 0 0 0 0 0 -1

After I added it too known_values.json, it cast several time "ERROR on roots".

Roots it found first time (as written in my json file): SArray{Tuple{9},Int64,1,9}[[-1, 2, 2, 1, 0, 0, 0, 0, 0], [0, -1, 0, 1, 0, 0, 0, 0, 0], [0, 0, -1, 0, 0, 0, 0, 0, 0], [0, 0, 0, -1, 0, 0, 0, 0, 0], [0, 0, 0, 0, -1, 2, 2, 1, 0], [0, 0, 0, 0, 0, -1, 0, 1, 0], [0, 0, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, 0, 0, 0, 0, -1, 0], [1, 0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 1, 0, 0, 0, 1], [2, -1, -2, -1, 2, -1, -2, -1, 2]]

Second time:
[[-1, 2, 2, 1, 0, 0, 0, 0, 0], [0, -1, 0, 1, 0, 0, 0, 0, 0], [0, 0, -1, 0, 0, 0, 0, 0, 0], [0, 0, 0, -1, 0, 0, 0, 0, 0], [0, 0, 0, 0, -1, 2, 2, 1, 0], [0, 0, 0, 0, 0, -1, 0, 1, 0], [0, 0, 0, 0, 0, 0, -1, 0, 0], [0, 0, 0, 0, 0, 0, 0, -1, 0], [2, -1, -2, -1, 2, -1, -2, -1, 2], [1, 0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 1, 0, 0, 0, 1]]

Looks like they coincide as sets. Can you verify how the roots are arranged by distance? Why would they be permuted? Sometimes they come out exactly as stored, and sometimes permuted.

E.g. when enumerating roots, the norm of a root returned by roots_decomposed_into is known, so that we could ship the norm with the root already, so that we don't have to compute it later on, but is it actually useful?
Probably other such double computations could be killed at other places.

I see a horrible 10x performance with the commit 513abb39fcf58d87d54f260a316568759c82d3d8 and what I had before (say 10 commits before).
Need to understand what the reason is, but I don't see it, since no big code change happened that I can see…
Suspects are:

• Using channels (but the commit introducing them said the impact was not big either way)
• Making the channel's buffer bigger (but the profiler doesn't say the problem lies here as far as I can see)
• Changes in qsolve (the profiler seems to say that qsolve_iterative is very resource hungry, but that's to be expected, and I changed almost nothing there!)
• Setting warmstart = true for the cbc optimizer (but the profiler doesn't seem to say it's that)

Because as it is, we very often have only one, and this probably hides buggy behaviour in the code.

For instance, if a collection of halfspaces is not minimal, is it true that the polyhedron they define is finite iff the diagram associated to the larger-than-needed collection is finite?

We have a bijection between hyperplanes in R^{n,1} intersecting H^n and hyperplanes in H^n.
How careful do we need to be with our roots to ensure that they define proper hyperplanes in H^n?

Start with the following:

Let's say it's distance comparison.
You need to check smth like: if H_1, H_2 are hyperplanes such that v_0 belongs to their negative half-spaces and, moreover, H_1 is contained in the negative subspace of H_2, then the distance from v_0 to H_1 is less than the distance from v_0 to H_2 (for any choice of v_0 as described).
I think this will be the 'formal" argument, but somehow it seems very intuitive.

Mark has some comments about some decompositions being impossible due to some relatively easy(?) divisibility conditions. Should see if we can use that easily?

Presumably the roots in the fundamental cone are simply those which define hyperplanes at distance zero, so aren't really different from the others.
It might be simpler to merge the searches for roots in just one loop.

For the diagonal quadratic form [1]*14+[-2], I get the following:

[ Info: Initialized root iterator.
[ Info: Got all roots at distance zero.
[ Info: Creating a channel for root enumeration (buffer size 16)
[ Info: Starting the root enumeration process.
[ Info: Found candidate roots for fundamental cone (numbering 392).
[ Info: Returning cone roots, numbering 14.
[ Info: Got the roots of a fundamental cone.
[ Info: Built the corresponding Coxeter diagram.
[ Info: Found new satisfying root: [2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1].
[ Info: Found new satisfying root: [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1].
[ Info: Found new satisfying root: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2].

Julia has exited.

It finds new roots extremely slowly compared to Guglielmetti's code, and also exits unexpectedly at some point. I have no idea what causes this issue.

Presumably better?

Citing him, this could save a lot of time for the finite volume check on some graphs.
See screenshot:

The main difficulties arise at two places:

• The root enumeration procedure must be overhauled, since the current one is tailor-made for ZZ.
• The hyperplane necessity check works only in ZZ.

Maybe instead of adapting the root enumeration procedure, we could just enumerate roots in any order: we can probably still see if some finite volume polyhedron appears, but it won't be the minimal one yet…

For the hyperplane necessity, maybe we can actually do precise enough "float" arithmetic to get a first guess, and then take the "witness" the solver gives us and try to verify that it or some close enough element is a valid formal witness?

Say we build a graph all of whose vertices are the roots at distance zero from the basepoint, and with an edge between two roots if the angle defined by their hyperplanes is acute.

Is there a relation between cliques in this graph and choices of fundamental cones?

Ideally have cache file that can be loaded on each run and updated as needed.

See the implementation: here.

What they seem to do is:

1. Take all roots in v₀^\perp = V₁.
2. Add them to the cone as necessary.

I have a problem with step 1.
The idea is that any element of the lattice can be represented as a₀v₀ + v₁ + w, where v₁ is in V₁ and w is a coset representative for the sublattice v₀ (+) V₁ of our big lattice.
Then, the procedure s.RootsDecomposedInto(s.V([0]*s.n), k) returns elements where both w and a₀ are zero, so that what remains is the v₁ component, which is indeed an element of V₁.
But it could also happen that an element of V₁ cannot be written as 0*v₀ + v₁ + 0, and thus, not come from the RootsDecomposedInto procedure.
The condition for an element of the form a₀v₀ + v₁ + w to be in V₁ is that (a₀v₀ + v₁ + w | v₀ ) == 0, which means only that a₀ == -(v₀|w) / (v₀|v₀).

Another thing that I don't understand:
Afaik, the condition that roots defining the polytope have pairwise non-positive product should hold for all of them, "cone" roots included: but B&P do not check this condition, and adding it actually changes the results! what is the reason?

Am I mistaken? pinging @sashakolpakov

To get better performance, e.g. in #33, it could be a good idea to implement a specific root enumerator since this case is much simpler, probably.

Memory usage shouldn't be too high, yet it quickly gets in the order of gigabytes…