Repository: simple-crypto/SCALib (language: python)

Issue:
LDA: handle case where we don't want to do dimensionality reduction
Currently setting  `p == ns` crashes.

When `p==ns`, no dimensionality reduction is needed, so the whole eigenvectors part can be skipped, and the projection matrix can be set to the identity matrix. https://github.com/simple-crypto/SCALib/blob/7cbe95216664d6678a52d9aaa4a1f8217b1e628d/src/scalib_ext/scalib/src/lda.rs#L219
LDA: handle case where we don't want to do dimensionality reduction
Currently setting  `p == ns` crashes.

When `p==ns`, no dimensionality reduction is needed, so the whole eigenvectors part can be skipped, and the projection matrix can be set to the identity matrix. https://github.com/simple-crypto/SCALib/blob/7cbe95216664d6678a52d9aaa4a1f8217b1e628d/src/scalib_ext/scalib/src/lda.rs#L219


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 6d3e848f..33bf15d7 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -5,7 +5,9 @@ Changelog
 Not released
 ------------
 
-* Fix freeze when there is an error in large SNR computations (i.e., when n_vars*n_samples*n_traces > 2**33).
+* Fix deadlock when there is an error in large SNR computations (i.e., when n_vars*n_samples*n_traces > 2**33).
+* Allow LDA to behave like simple pooled gaussian templates (#22)
+* Refresh build system (Tox version 4, improved CI).
 
 v0.4.3 (2022/10/27)
 -------------------
diff --git a/README.rst b/README.rst
index d1a74e46..b3f9200b 100644
--- a/README.rst
+++ b/README.rst
@@ -29,6 +29,19 @@ See install_ for details.
  
 .. _install: https://scalib.readthedocs.io/en/stable/index.html#install
 
+Dependencies upgrade policy
+===========================
+
+- For python and python packages, we follow NEP29_.
+- For rust: latest stable version.
+- OS support:
+
+  * Pre-built wheels for manylinux_ (supporting last two Ubuntu LTS) on x86_64.
+  * Pre-built wheels for Windows 10 on x86_64.
+  * Other: build yourself (CI configuration welcome).
+
+.. _NEP29: https://numpy.org/neps/nep-0029-deprecation_policy.html
+.. _manylinux: https://github.com/pypa/manylinux
 
 About us
 ========
@@ -41,6 +54,7 @@ maintained in that context.
 Contributions and Issues
 ========================
 
+Contributions welcome !
 See contribution_ and `DEVELOP.rst <DEVELOP.rst>`_.
 
 .. _contribution: https://scalib.readthedocs.io/en/stable/index.html#contributions-and-issues
diff --git a/docs/index.rst b/docs/index.rst
index 64ae9d83..7fee137c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -143,7 +143,7 @@ Where is SCALib used ?
 
 We strongly appreciate if you could mention to us your usage SCALib 
 for concrete projects so that we can add you to the lists below. Please send 
-an email to Olivier Bronchain.
+an email to Olivier Bronchain and Gaëtan Cassiers, or directly open a pull request.
 
 Scientific publications
 -----------------------
@@ -194,14 +194,18 @@ of SCALib starting from the first version. The SCALib project is part of `SIMPLE
 
 Contributions and issues
 ========================
+
 We are happy to take any suggestion for features would be useful for
-side-channel evaluators.  If you want to contribute code to the project, please
-visit `DEVELOP
-<https://github.com/simple-crypto/SCALib/blob/main/DEVELOP.rst>`_ for relevant
-information as well as the Contributor License Agreement (`CLA
-<https://www.simple-crypto.dev/organization>`_) of SIMPLE-Crypto. Please open
-an `issue <https://github.com/simple-crypto/SCALib/issues>`_ on the GitHub repository for further questions, bug report or feature
-requests. 
+side-channel evaluators.
+We very much welcome code, documentation, examples, etc. contributions.
+`DEVELOP <https://github.com/simple-crypto/SCALib/blob/main/DEVELOP.rst>`_
+contains a few informations and guides for working with SCALib code.
+Please open an `issue <https://github.com/simple-crypto/SCALib/issues>`_ on the
+GitHub repository for further questions, bug report or feature requests. 
+
+All code contributions are subject to the Contributor License Agreement (`CLA
+<https://www.simple-crypto.dev/organization>`_) of SIMPLE-Crypto, which ensures
+a thriving future for open-source hardware security.
 
 License
 =======
diff --git a/src/scalib_ext/scalib/src/lda.rs b/src/scalib_ext/scalib/src/lda.rs
index 4ca02d90..3f470138 100644
--- a/src/scalib_ext/scalib/src/lda.rs
+++ b/src/scalib_ext/scalib/src/lda.rs
@@ -211,13 +211,22 @@ impl LDA {
     ) -> Self {
         let ns = sw.shape()[0];
         let nc = means_ns.shape()[0];
-        // ---- Step 1
-        // compute the projection
-        // Partial generalized eigenvalue decomposition for sb and sw.
-        let solver =
-            geigen::GEigenSolverP::new(&sb.view(), &sw.view(), p).expect("failed to solve");
-        let projection = solver.vecs().into_owned();
-        assert_eq!(projection.dim(), (ns, p));
+        let projection = if p == ns {
+            // no LDA, simple pooled gaussian templates
+            // This is suboptimal since we then have to pay multiplication with the identity, but
+            // it should not matter much (we are probably in a cheap case anyway).
+            ndarray::Array2::eye(ns)
+        } else {
+            // LDA here
+            // ---- Step 1
+            // compute the projection
+            // Partial generalized eigenvalue decomposition for sb and sw.
+            let solver =
+                geigen::GEigenSolverP::new(&sb.view(), &sw.view(), p).expect("failed to solve");
+            let projection = solver.vecs().into_owned();
+            assert_eq!(projection.dim(), (ns, p));
+            projection
+        };
 
         // ---- Step 2
         // means per class within the subspace by projecting means_ns


Test patch:
diff --git a/tests/test_lda.py b/tests/test_lda.py
index aee35741..bc82cd91 100644
--- a/tests/test_lda.py
+++ b/tests/test_lda.py
@@ -40,7 +40,7 @@ def test_lda_pickle():
     dumped_lda = pickle.dumps(lda)
     lda = pickle.loads(dumped_lda)
 
-    lda_ref = LDA_sklearn(solver="eigen")
+    lda_ref = LDA_sklearn(solver="eigen", n_components=n_components)
     lda_ref.fit(traces, labels)
 
     lda_projection = lda.lda.get_projection()
@@ -132,31 +132,42 @@ def test_lda():
     traces_t = traces_t - means_check[labels, :]
     cov_check = np.cov(traces_t.T)
 
-    # traces_t = (lda_projection[:, :].T @ traces.T).T
-    # means_check_rust = np.zeros((nc, n_components))
-    # for i in range(nc):
-    #    I = np.where(labels == i)[0]
-    #    means_check_rust[i, :] = np.mean(traces_t[I, :], axis=0)
-    # traces_t = traces_t - means_check_rust[labels, :]
-    # cov_check_rust = np.cov(traces_t.T)
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    traces += m[labels]
+
+    prs = lda.predict_proba(traces)
+    prs_ref = lda_ref.predict_proba(traces)
+
+    assert np.allclose(prs, prs_ref, rtol=1e-2)
 
-    # assert np.allclose(means_check_rust, lda.lda.get_means().T )
-    #    assert(np.allclose(cov_check,lda.lda.get_cov()))
+
+def test_lda_noproj():
+    ns = 3
+    n_components = 3
+    nc = 4
+    n = 500
+
+    m = np.random.randint(0, 100, (nc, ns))
+    traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
+    labels = np.random.randint(0, nc, n, dtype=np.uint16)
+    traces += m[labels]
+
+    lda = LDAClassifier(nc, n_components, ns)
+    lda.fit_u(traces, labels, 1)
+    lda.solve()
+
+    lda_ref = LDA_sklearn(solver="eigen", n_components=n_components)
+    lda_ref.fit(traces, labels)
 
     traces = np.random.randint(0, 10, (n, ns), dtype=np.int16)
     labels = np.random.randint(0, nc, n, dtype=np.uint16)
     traces += m[labels]
 
     prs = lda.predict_proba(traces)
-    traces_t = (lda_ref.scalings_[:, :n_components].T @ traces.T).T
-    prs_ref = np.zeros((len(traces), nc))
-    for x in range(nc):
-        prs_ref[:, x] = scipy.stats.multivariate_normal.pdf(
-            traces_t, means_check[x], cov_check
-        )
-    prs_ref = (prs_ref.T / np.sum(prs_ref, axis=1)).T
+    prs_ref = lda_ref.predict_proba(traces)
 
-    assert np.allclose(prs, prs_ref, rtol=1e-2)
+    assert np.allclose(prs, prs_ref, rtol=1e-2, atol=1e-5)
 
 
 if __name__ == "__main__":