skew-Hadamard matrices and related srg's

[dimpase]

implementing few basic constructions of skew-Hadamard matrices and 3 families of assoc. srg's.

[dimpase]

a draft implementation of the corresponding SRGs

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3453e16

merged...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​bf8dcc1

added Goethals-Seidel construction; this gives orders 36 and 52

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​9e0ff0f

added Goethals-Seidel construction; this gives orders 36, 52, and 92

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f4e3a06

graph constructions added (prelim version, refs etc missing)

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​4aa73f7

graph constructions added, with skew-normalizing Hadamard mats

[ncohen]

Hello Dima,

I just merged #19309 with the latest beta. Could you do *rebase* this branch on top of it (not merge) so that your commits appear on top of that other branch and are easier to review ?

Nathann

[ncohen]

(right now the branch consists of 13 commits, 10 of which are merges with tickets that belong to the latest release)

[ncohen]

I just did it: the rebased branch corresponding to this ticket can be found at public/19418.

Nathann

[dimpase]

New commits:

​8607575	trac #19309: Polhill strongly regular graphs
​2cd7910	trac #19309: Note for later
​7bbefd1	trac #19309: Merged with 6.9.rc1
​acba28b	trac #19309: Broken doctest
​732f243	trac #19309: Merged with 6.10.beta0
​e65da7e	trac #19309: Merged with 6.10.beta1
​c5a7eee	trac #19418: skew-Hadamard matrices and related srg's

[dimpase]

I assume you meant switching to this branch by 'work' :-) By the way, I've left a function _L_g_n_params in the source (it's actually not used). If you think it should be removed, I will do so.

[ncohen]

Helloooooo Dima,

First-pass review:

• Why don't you call normalise_hadamard instead of adding a keyword inside of hadamard_matrix_paleyI (without changing hadamard_matrix_paleyII)? Those two functions are not directly available to the users, who are meant to call hadamard_matrix directly.

• is_hadamard_martrix -- could you move 'skew=False' to before 'verbose=False' ? Usually the verbosity/check flags are the last ones to appear as they do not change the behaviour of the functon.

• Call is_skew_symmetric instead of doing it yourself. In theory you shouldn't have to create two copies of the matrix in memory in order to check that.

• zero_position=1 -- you have to write documentation for private functions too. Does [the paraeter that you add] appear in the lemma cited in the docstring?

• zero_position=1 -- if you meant True, write True.

• _circulant_matrix -- if it does not exist yet, this constructor should be added in the matrix/ code and be available through matrix.<tab>

• mod(j-i,n) is (j-i)%n. One day you will have to accept that you write Python code.

• Instead of a _GS_skew_hadamard function that encodes 4 matrices, why don't you reserve cases `if n == 36' (and others) in the main function skew_hadamard_matrix? That would also solve the problem that this function has no documentation whatsoever.

• 'and' (see wikipedia: ​https://en.wikipedia.org/wiki/Bitwise_operation#AND) is a useful operator in computer science. E.g.:

+    if M is None: # try Williamson-Goethals-Seidel construction
+        if _GS_skew_hadamard(n, existence=True):

becomes

+    if M is None and_GS_skew_hadamard(n, existence=True):

• switch_skewhad_pow2 -- please respect the current standard for the names in graphs.<tab>. Upper case, no underscore, ends with Graph. Actually, it could be named 'SwitchSkewhadGraph? and only mention the pow2` part in the docstring.

• graphs.Pseudo_L_2n_4n_m_1? Seriously? If you cannot give it a clear and meaningful name easily, then take it as a sign that this function should not be exposed to the users directly.

• _L_g_n_params -- yes, please remove this function if you do not use it.

Nathann

[dimpase]

Replying to ncohen:

First-pass review:

thanks, and apologies for sloppy code (back to Oxford with huge jetlag)

• Why don't you call normalise_hadamard instead of adding a keyword inside of hadamard_matrix_paleyI (without changing hadamard_matrix_paleyII)? Those two functions are not directly available to the users, who are meant to call hadamard_matrix directly.

normalize_hadamard does not do what I need; I need a different type of normalization, in fact, where the matrix has the 1st row consisting of 1s, and still H+H.T==2*I, i.e. the 1st column of H is all -1 (with exception of the top left entry)

Should I perhaps introduce skew_normalize_hadamard, that can only be applied to skew Hadamard matrices? (I'd have uses for it in graph constructions, too)

• is_hadamard_martrix -- could you move 'skew=False' to before 'verbose=False' ? Usually the verbosity/check flags are the last ones to appear as they do not change the behaviour of the functon.

right. I forgot that keyword args can be positional, too...

• Call is_skew_symmetric instead of doing it yourself. In theory you shouldn't have to create two copies of the matrix in memory in order to check that.

• zero_position=1 -- you have to write documentation for private functions too. Does [the paraeter that you add] appear in the lemma cited in the docstring?

• zero_position=1 -- if you meant True, write True.

• _circulant_matrix -- if it does not exist yet, this constructor should be added in the matrix/ code and be available through matrix.<tab>

OK

• mod(j-i,n) is (j-i)%n. One day you will have to accept that you write Python code.

• Instead of a _GS_skew_hadamard function that encodes 4 matrices, why don't you reserve cases `if n == 36' (and others) in the main function skew_hadamard_matrix? That would also solve the problem that this function has no documentation whatsoever.

I'd rather add documentation, thanks for reminding. I find skew_hadamard_matrix already hard to read, and it would positively turn ugly if I start adding that explicit data for matrices there.

• 'and' (see wikipedia: ​https://en.wikipedia.org/wiki/Bitwise_operation#AND) is a useful operator in computer science. E.g.:

+    if M is None: # try Williamson-Goethals-Seidel construction
+        if _GS_skew_hadamard(n, existence=True):

becomes

+    if M is None and_GS_skew_hadamard(n, existence=True):

oops - there is an extra if there that should go...

• switch_skewhad_pow2 -- please respect the current standard for the names in graphs.<tab>. Upper case, no underscore, ends with Graph. Actually, it could be named 'SwitchSkewhadGraph? and only mention the pow2` part in the docstring.

• graphs.Pseudo_L_2n_4n_m_1? Seriously? If you cannot give it a clear and meaningful name easily, then take it as a sign that this function should not be exposed to the users directly.

How about we have a PseudoLatinSquaresGraph(m,n) or PseudoOrthogonalArrayBlockGraph(m,n) ? (I'd prefer the former name; or maybe we can have both?)

It would include Pseudo_L_2n_4n_m_1 and OrthogonalArrayBlockGraph as subcases. (Perhaps there are more subcases, in fact, but we don't have to worry about it now).

• _L_g_n_params -- yes, please remove this function if you do not use it.

sure.

[ncohen]

Helloooooooooooooooooooo,

(back to Oxford with huge jetlag)

Wow. Common timezone again. Cool ;-)

normalize_hadamard does not do what I need; I need a different type of normalization, in fact, where the matrix has the 1st row consisting of 1s, and still H+H.T==2*I, i.e. the 1st column of H is all -1 (with exception of the top left entry)

Should I perhaps introduce skew_normalize_hadamard, that can only be applied to skew Hadamard matrices? (I'd have uses for it in graph constructions, too)

Sounds right. Then you can have a proper documentation and stuff. If you plan to use it elsewhere that's even better.

I'd rather add documentation, thanks for reminding. I find skew_hadamard_matrix already hard to read, and it would positively turn ugly if I start adding that explicit data for matrices there.

Okayokay.

How about we have a PseudoLatinSquaresGraph(m,n) or PseudoOrthogonalArrayBlockGraph(m,n) ? (I'd prefer the former name; or maybe we can have both?)

+1 to the former. With a description of what it is in its doc, if possible ^^;

Nathann

[dimpase]

Replying to ncohen:

(back to Oxford with huge jetlag)

Wow. Common timezone again. Cool ;-)

an hour difference, still :-)

normalize_hadamard does not do what I need; I need a different type of normalization, in fact, where the matrix has the 1st row consisting of 1s, and still H+H.T==2*I, i.e. the 1st column of H is all -1 (with exception of the top left entry)

Should I perhaps introduce skew_normalize_hadamard, that can only be applied to skew Hadamard matrices? (I'd have uses for it in graph constructions, too)

Sounds right. Then you can have a proper documentation and stuff. If you plan to use it elsewhere that's even better.

oops, this won't fly. The problem is that once a skew-Hadamard matrix was de-skewed (e.g. by the usual normalization procedure, producing the 1st row and column of 1s), there is no general way known to make it skew again. (and in fact not every Hadamard matrix can be made skew, AFAIK).

That is, I have to disable normalization right at the point it is done in hadamard_matrix_paleyI. (writing an ad hoc code to undo this specific normalization seems silly).

[dimpase]

Replying to ncohen:

Helloooooo Dima,

First-pass review:

• Why don't you call normalise_hadamard instead of adding a keyword inside of hadamard_matrix_paleyI (without changing hadamard_matrix_paleyII)? Those two functions are not directly available to the users, who are meant to call hadamard_matrix directly.

• is_hadamard_martrix -- could you move 'skew=False' to before 'verbose=False' ? Usually the verbosity/check flags are the last ones to appear as they do not change the behaviour of the functon.

• Call is_skew_symmetric instead of doing it yourself. In theory you shouldn't have to create two copies of the matrix in memory in order to check that.

this won't fly either, as a skew Hadamard matrix H is not skew-symmetric. What is skew-symmetric is the matrix H-I, so it has to be created anyway

(by the way, this explains why making a Hadamard matrix skew is not too obvious; indeed, the trick mentioned in the docs of is_skew_symmetrizable would only work if H had constant diagonal.)

[ncohen]

Yo,

oops, this won't fly. The problem is that once a skew-Hadamard matrix was de-skewed (e.g. by the usual normalization procedure, producing the 1st row and column of 1s), there is no general way known to make it skew again. (and in fact not every Hadamard matrix can be made skew, AFAIK).

That is, I have to disable normalization right at the point it is done in hadamard_matrix_paleyI. (writing an ad hoc code to undo this specific normalization seems silly).

I cannot say that I understand what you have in mind, but as long as what happens is made clear in the documentation I guess that it is okay.

this won't fly either, as a skew Hadamard matrix H is not skew-symmetric. What is skew-symmetric is the matrix H-I, so it has to be created anyway

Then create H-I, but don't create -S nor S.T.

Nathann

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0a5e36f	adding missing docs and tests, part I
​c284eda	moved circulant to matrix.<tab>

[ncohen]

a couple of remarks about the last commit:

• except without the name of a specific exception makes Jeroen scream.

• It is a weird that the value of 'sparse' is ignored when 'v.is_sparse' says otherwise. The expected behaviour would be sparse=None by default (auto detect) and use the user-provided value otherwise.

[dimpase]

Replying to ncohen:

a couple of remarks about the last commit:

• except without the name of a specific exception makes Jeroen scream.

OK, sure, I will make it catch NoneType

• It is a weird that the value of 'sparse' is ignored when 'v.is_sparse' says otherwise. The expected behaviour would be sparse=None by default (auto detect) and use the user-provided value otherwise.

I want to implement something like this:

    def f(vector v):
         stuff...

    def f(list v, sparse=False):
         stuff...

and this is what my code does.

The user-provided value of sparse for vector v is present in v already. If the user wants to make it explicit, he can call the function with the parameter vector(v, sparse=whatever) instead of just v

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6f54887

function renaming, documenting, and doctesting, and removing extra if

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6b80f00

added AttributeError to check for in except:

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​135c0ea	for a skew Hadamard matrix, check that the diagonal entries are all 1
​c04916b	return skew-normalized Hadamard mats by default

[dimpase]

looks like I have addressed all the points raised; (I didn't mention in commit messages that the code now does checking of skew-Hadamardness inplace). ready for another round, hopefully...

---

New commits:

​135c0ea	for a skew Hadamard matrix, check that the diagonal entries are all 1
​c04916b	return skew-normalized Hadamard mats by default

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​535e109

putting reference in the right file

[dimpase]

now make doc-clean && make passes.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b3fb30b

Merge branch 'develop' of git://trac.sagemath.org/sage into skewhad

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​01763fd

remove #random from paleyI/II and fix the test

[dimpase]

Replying to git:

Branch pushed to git repo; I updated commit sha1. New commits:

​01763fd

remove #random from paleyI/II and fix the test

I guess these # random were leftovers from old code. Let's get rid of them.

[ncohen]

I want to implement something like this:

    def f(vector v):
         stuff...

    def f(list v, sparse=False):
         stuff...

and this is what my code does.

The user-provided value of sparse for vector v is present in v already. If the user wants to make it explicit, he can call the function with the parameter vector(v, sparse=whatever) instead of just v

Usually the explicit value overrides the implicit value:

sage: Graph(graphs.PetersenGraph())._backend
<type 'sage.graphs.base.sparse_graph.SparseGraphBackend'>
sage: Graph(graphs.PetersenGraph(),sparse=False)._backend
<type 'sage.graphs.base.dense_graph.DenseGraphBackend'>

It is not the case in your code.

Nathann

[ncohen]

could be something like that

def thing(v,sparse=None):
    if sparse is None:
        try:
            sparse = v.is_sparse()
        except AttributeError:
            sparse = False

[dimpase]

Replying to ncohen:

could be something like that

def thing(v,sparse=None):
    if sparse is None:
        try:
            sparse = v.is_sparse()
        except AttributeError:
            sparse = False

I'll be happy to fix it as you prefer. (You could also add a commit, of course :))

[ncohen]

Okayokay. Well I'm reviewing the code at the moment.

By the way on my way to here I finished a wonderful script that extracts tables from a pdf file. Some french guys need that to turn into open data what the government releases to feel open (i.e. ugly pdf files).

​http://www.performance-publique.budget.gouv.fr/sites/performance_publique/files/farandole/ressources/2016/pap/pdf/jaunes/jaune2016_commissions.pdf

But we'll show 'em.

Nathann

[ncohen]

Helloooooooooo Dima,

Here is another pass:

• I don't think that you need to import functions at the top of every single '::' block. If you imported them in a previous doctest of the same function they are usually around. I say 'usually' because I don't think I read that anywhere, but it just works. If the order in which doctests are run can be 'random', I think that it just shuffles the functions and not the doctests inside of a function.

• docstring of zero_position: bad alignment of 'place it first'. Misses one space.

• Same line: the sentence is not complete -- what about adding 'When it' at its beginning?

• You did not doctest your new flag. Also, the code *looks* wrong.

• Docstring of williamson_goethals_seidel_skew_hadamard_matrix. Things seem to miss from the first sentence of the paragraph. Or it's me being sick again. P.S.: syntactically, a ')' is missing at least.

• Same function. The following doc should be in the function that actually returns those matrices:

  +    Matrices for `n=36` and `52` are given in [GS70s]_. Matrices for `n=92` are given
  +    in [Wall71]_.
  

• _GS_skew_hadamard -- the name of that function is a bit too vague. What about GS_skew_hadamard_smallcases? You don't necessarily need to make it private, by the way. It lives in a module and does not appear in the global namespace, so only the interested guys will find it. And it's not technically complicated to use, the input/output is well defined and reliable. It can't be misused.

• private functions need doctests too. 'cdef' is the only situation in which you don't need them.

• def pmtoZ(s): could be defined in the only 'if' that requires it (not very important).

• docstring of skew_hadamard_matrix -- in the block of doc about 'existence', one reads 'but that design does not exist'. Could you replace 'design' by 'matrix'?

• I don't see why the following doctest has a 'long' flag

  sage: skew_hadamard_matrix(784,existence=True) # long time
  

  If you actually built the matrix, however, that would be correct.

• I don't think that check=False is useful when existence=True

  skew_hadamard_matrix(n//2,existence=True, check=False)
  

• In the n%8==0 case you implement the product construction of skew hadamard matrices with a 2x2 matrix only. Is there a reason to that?

• if _GS_skew_hadamard(n, existence=True) -- it hardly matters, but it would be 'cleaner' to call this (quick) function before checking the more complicated (recursive) constructions.

• Shouldn't you also check that the matrix is 'skew-normalized', as we do with the regular normalization?

  +    if skew_normalize:
  +        dd = diagonal_matrix(M[0])
  +        M = dd*M*dd
  +    if check:
  +        assert is_hadamard_matrix(M, normalized=False, skew=True)
  

• If you are too humble to do it yourself, I will change that to 'Return the Pasechnik graph of order n' :-)

  +    Pseudo-`OA(2n-1,4n-1)`-graph from a skew Hadamard matrix of order `4n`
  

• Please add your new graphs to the graphs.<tab> constructor.

Nathann

[dimpase]

Replying to ncohen:

• In the n%8==0 case you implement the product construction of skew hadamard matrices with a 2x2 matrix only. Is there a reason to that?

huh? I have there

    elif n % 8 == 0:
        if skew_hadamard_matrix(n//2,existence=True, check=False):
            if existence:
                return true()
            H = skew_hadamard_matrix(n//2,check=False)
            M = block_matrix([[H,H], [-H.T,H.T]])

        else: # try Williamson construction
            for d in divisors(n)[2:-2]: # skip 1, 2, n/2, and n
             ... more stuff

so I first check if I can multiply with 2x2 matrix, and if not, I try Williamson construction.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​fd7894a	check that skew_normalize worked
​3503d62	rest of fixes for hadamard_matrix.py

[dimpase]

Replying to ncohen:

• I don't think that you need to import functions at the top of every single '::' block. If you imported them in a previous doctest of the same function they are usually around. I say 'usually' because I don't think I read that anywhere, but it just works. If the order in which doctests are run can be 'random', I think that it just shuffles the functions and not the doctests inside of a function.

OK. Removed some of these...

• docstring of zero_position: bad alignment of 'place it first'. Misses one space.

• Same line: the sentence is not complete -- what about adding 'When it' at its beginning?

• You did not doctest your new flag. Also, the code *looks* wrong.

well, actually in this case the code is not doing what is advertised; without my flag the matrix for n=4k-1 is not symmetric w.r.t. the anti-diagonal; the submatrix obtained by removal of the 1st row and the 1st column is symmetric w.r.t. the anti-diagonal. I changed the doc accordingly. And added a doctest.

• Docstring of williamson_goethals_seidel_skew_hadamard_matrix. Things seem to miss from the first sentence of the paragraph. Or it's me being sick again. P.S.: syntactically, a ')' is missing at least.

fixed.

• Same function. The following doc should be in the function that actually returns those matrices:

  +    Matrices for `n=36` and `52` are given in [GS70s]_. Matrices for `n=92` are given
  +    in [Wall71]_.
  

moved

• _GS_skew_hadamard -- the name of that function is a bit too vague. What about GS_skew_hadamard_smallcases? You don't necessarily need to make it private, by the way. It lives in a module and does not appear in the global namespace, so only the interested guys will find it. And it's not technically complicated to use, the input/output is well defined and reliable. It can't be misused.

renamed as you proposed.

• private functions need doctests too. 'cdef' is the only situation in which you don't need them.

OK, OK... added

• def pmtoZ(s): could be defined in the only 'if' that requires it (not very important).

well, if we added more data in +- - format, it would be at the right place as it is.

• docstring of skew_hadamard_matrix -- in the block of doc about 'existence', one reads 'but that design does not exist'. Could you replace 'design' by 'matrix'?

done

• I don't see why the following doctest has a 'long' flag

  sage: skew_hadamard_matrix(784,existence=True) # long time
  

  If you actually built the matrix, however, that would be correct.

good catch - fixed

• I don't think that check=False is useful when existence=True

  skew_hadamard_matrix(n//2,existence=True, check=False)
  

uh oh... fixed.

• if _GS_skew_hadamard(n, existence=True) -- it hardly matters, but it would be 'cleaner' to call this (quick) function before checking the more complicated (recursive) constructions.

Oh, well, if you feel like changing this, go ahead...

• Shouldn't you also check that the matrix is 'skew-normalized', as we do with the regular normalization?

  +    if skew_normalize:
  +        dd = diagonal_matrix(M[0])
  +        M = dd*M*dd
  +    if check:
  +        assert is_hadamard_matrix(M, normalized=False, skew=True)
  

added a check.

• If you are too humble to do it yourself, I will change that to 'Return the Pasechnik graph of order n' :-)

  +    Pseudo-`OA(2n-1,4n-1)`-graph from a skew Hadamard matrix of order `4n`
  

if you must...

• Please add your new graphs to the graphs.<tab> constructor.

actually, I removed them, but I'll add them again...

---

New commits:

​fd7894a	check that skew_normalize worked
​3503d62	rest of fixes for hadamard_matrix.py

---

New commits:

​fd7894a	check that skew_normalize worked
​3503d62	rest of fixes for hadamard_matrix.py

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​595f7e9

added graphs to graphs.<tab>, and doc(test)s additions and fixes

[ncohen]

sage: skew_hadamard_matrix(816)
...
TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 12 by 12 dense matrices over Integer Ring' and 'Full MatrixSpace of 68 by 68 dense matrices over Integer Ring'

[ncohen]

Can you tell me why you apparently refuse my request about the behaviour of 'sparse' in 'circulant'?

Nathann

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​06fd9b4

fix the bug, and add docs

[dimpase]

Replying to ncohen:

Can you tell me why you apparently refuse my request about the behaviour of 'sparse' in 'circulant'?

oops - forgot about it. Tonight I'll have more time to look at it.

The last commit fixes the n=816 bug (--facepalm--) - thanks for the example. I checked now the values of n up to 1200, they appears to work fine...

[ncohen]

Okayyy. Thanks for the fix ;-)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c50c78d

fixing circulant's behaviour

[ncohen]

Okayyyyyyyyyyyy... Thanks for this branch ! ;-)

Nathann

[dimpase]

OK, with all the changes, patchbots (most of them) seem to be happy. And the reviewer?

[ncohen]

He's happy too.

[dimpase]

Thanks! By the way, _check_database() now says for me:

...
Sage can build a (1024, 462, 206, 210), Brouwer's database cannot
Sage can build a (1024, 561, 308, 306), Brouwer's database cannot

In Andries Brouwer's database:
- 448 impossible entries
- 2950 undecided entries
- 1140 realizable entries (Sage misses 36 of them)

On the other hand the graphs above can be found in ​http://www.win.tue.nl/~aeb/graphs/srg/srgtab1001-1050.html Does it mean that Sage's copy of the DB needs an update?

OK, so should I next take out Mathon's v=45?

[ncohen]

Does it mean that Sage's copy of the DB needs an update?

Yes. We gave him info to update his database but we didn't update our copy. By the way, the difference between our copy and his can be seen as 'one of the side-effects' of the work we have been doing. We may need that info.

OK, so should I next take out Mathon's v=45?

Take as many as you can out. I don't think that there is any left that I can implement: those that I thought manageable were apparently wrong constructions, and the authors 'play dead'.

Nathann