Fix #464: crash of Pauli-flow finding algorithm due to incorrect numba signatures.

[matulni]

This commit proposes a fix to issue #464. The example triggering the crash is added to the test suit.

Detailed explanation

To improve the efficiency of just-in-time compiled functions in the graphix._linalg module, we indicate the memory layout of the array parameters. In particular, we write u8[:,::1] to indicate that we expect C_CONTIGUOUS arrays which is numpy's default. If we pass an F_CONTIGUOUS array to these functions, numba will crash as in the reported bug.

Computing the transpose of a C_CONTIGUOUS array with numpy.transpose only changes the metadata flags of the array to avoid unnecessary copying:

>>> import numpy as np
>>> a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> print(a.flags)
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
>>> print(a.transpose().flags)
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

This transformation occurred inadvertently in :func: graphix.flow._find_gpflow._compute_correction_matrix_general_case:

# Steps 5, 6 and 7
ker_flow_demand_matrix = flow_demand_matrix.null_space().transpose()

Why did all the tests pass ?

Very often, the kernel of the flow demand matrix is a one-row or one-column matrix. In this case, the corresponding numpy array is both C_CONTIGUOUS and F_CONTIGUOUS, so computing the transpose does not invalidate the numba decorator signature:

>>> import numpy as np
>>> a = np.array([[1, 2, 3]])
>>> print(a.flags)
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
>>> print(a.transpose().flags)
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

What is the fix?

The function numpy.ascontiguousarray allows to reexpress an F_CONTIGUOUS array as C_CONTIGUOUS. Note that if the original array is C_CONTIGUOUS in the first place, calling numpy.ascontiguousarray does not incur in an extra coyping:

import numpy as np
import timeit

a = np.eye(100)
at = a.transpose()

assert a.flags["C_CONTIGUOUS"]
assert at.flags["F_CONTIGUOUS"]

# Original matrix
assert np.ascontiguousarray(a).flags["C_CONTIGUOUS"]
assert a is np.ascontiguousarray(a)
t = timeit.Timer(lambda: np.ascontiguousarray(a))  
print(t.timeit(1_000_000))  # 0.06663083399871539

# Transpose 
assert np.ascontiguousarray(at).flags["C_CONTIGUOUS"]
assert a is not np.ascontiguousarray(at)
t = timeit.Timer(lambda: np.ascontiguousarray(at))  
print(t.timeit(1_000_000))  # 6.021290332999342

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 88.95%. Comparing base (d5343cf) to head (3649af5).
⚠️ Report is 1 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #465   +/-   ##
=======================================
  Coverage   88.95%   88.95%           
=======================================
  Files          44       44           
  Lines        6530     6530           
=======================================
  Hits         5809     5809           
  Misses        721      721           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[thierry-martinez]

It is a bit disappointing that functions in _linalg.py make assumptions about array layouts, especially since we don't document this. Are you sure these assumptions impact performance? Do benchmarks degrade if we use the signature @nb.njit("uint8[:,:](uint8[:,:], uint8[:,:])", parallel=True) instead?

[matulni]

Commit 5cf616c ensures that all numba-jitted functions called from graphix._linalg.py receive c-contiguous functions. New tests at the level of _linalg.py are added.

Further, initialisation of graphix._linalg.MatGF2 preserves the memory layout of the original data now, instead of generating C_CONTIGUOUS arrays which was an unreported bug.

[thierry-martinez]

Very nice! Just a small suggestion, but I approve already.

[thierry-martinez]

LGTM! Thanks!