basic testing equilibrate

project_euromir/__init__.py

@@ -44,7 +44,7 @@
     'Linux': '.so',
     'Darwin': '.dylib',
     'Windows': '.dll',
-    }
+}
 
 for _fname in _pathlib.Path(__file__).parent.iterdir():
     if _fname.suffix == _EXTS[_platform.system()]:

project_euromir/tests/__main__.py

@@ -30,6 +30,8 @@
 
 from unittest import main
 
+from .test_equilibrate import TestEquilibrate
+from .test_lbfgs import TestLBFGS
 from .test_linear_algebra import TestLinearAlgebra
 
 if __name__ == '__main__':

project_euromir/tests/equilibrate.py

@@ -36,15 +36,17 @@
 
 logger = logging.getLogger()
 
+
 def _cones_separation_matrix(zero, nonneg, second_order):
     """Sparse matrix that maps entries into which cone they belong to."""
     return sp.block_diag(
         [sp.eye(zero+nonneg)] + [np.ones((1, el)) for el in second_order]
-        + [1.]) # we add the one for use below
+        + [1.])  # we add the one for use below
+
 
-def hsde_ruiz_equilibration( # pylint: disable=too-many-arguments
-        matrix, b, c, dimensions, d=None, e=None, rho = 1., sigma = 1.,
-        eps_rows = 1E-4, eps_cols = 1E-4, max_iters=25):
+def hsde_ruiz_equilibration(  # pylint: disable=too-many-arguments
+        matrix, b, c, dimensions, d=None, e=None, rho=1., sigma=1.,
+        eps_rows=1E-4, eps_cols=1E-4, max_iters=25):
     """Ruiz equilibration of problem matrices for the HSDE system.
 
     :param matrix: Problem matrix.
@@ -104,7 +106,7 @@ def hsde_ruiz_equilibration( # pylint: disable=too-many-arguments
     e_and_sigma[-1] = sigma
 
     work_matrix = sp.diags(d_and_rho[:-1]
-        ) @ matrix @ sp.diags(e_and_sigma[:-1])
+                           ) @ matrix @ sp.diags(e_and_sigma[:-1])
     work_b = e_and_sigma[-1] * sp.diags(d_and_rho[:-1]) @ b
     work_c = d_and_rho[-1] * sp.diags(e_and_sigma[:-1]) @ c
 
@@ -128,9 +130,9 @@ def hsde_ruiz_equilibration( # pylint: disable=too-many-arguments
         norm_cols_and_b[-1] = np.linalg.norm(work_b)
 
         r1 = max(norm_rows_and_c[norm_rows_and_c > 0]
-            ) / min(norm_rows_and_c[norm_rows_and_c > 0])
+                 ) / min(norm_rows_and_c[norm_rows_and_c > 0])
         r2 = max(norm_cols_and_b[norm_cols_and_b > 0]
-            ) / min(norm_cols_and_b[norm_cols_and_b > 0])
+                 ) / min(norm_cols_and_b[norm_cols_and_b > 0])
 
         logger.info('Equilibration iter %s: r1=%s, r2=%s', i, r1, r2)
         if (r1-1 < eps_rows) and (r2-1 < eps_cols):
@@ -144,7 +146,7 @@ def hsde_ruiz_equilibration( # pylint: disable=too-many-arguments
                 norm_cols_and_b > 0]**(-0.5)
 
         work_matrix = sp.diags(d_and_rho[:-1]
-            ) @ matrix @ sp.diags(e_and_sigma[:-1])
+                               ) @ matrix @ sp.diags(e_and_sigma[:-1])
         work_b = e_and_sigma[-1] * sp.diags(d_and_rho[:-1]) @ b
         work_c = d_and_rho[-1] * sp.diags(e_and_sigma[:-1]) @ c
 

project_euromir/tests/lbfgs.py

@@ -1,3 +1,32 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
 """Pure Python implementation of L-BFGS for testing.
  
 We only need the multiplication of the gradient at the current point by the
@@ -20,9 +49,11 @@
     https://web.archive.org/web/20231002054213/https://aria42.com/blog/2014/12/understanding-lbfgs
 """
 
-from future import __typing__
+from __future__ import annotations
+
 import numpy as np
 
+
 def lbfgs_multiply(
     current_gradient: np.array,
     past_steps: np.array,
@@ -78,7 +109,7 @@ def lbfgs_multiply(
     for i in range(memory):
         betas[i] = rhos[i] * np.dot(past_grad_diffs[i], r)
         r += (alphas[i] - betas[i]) * past_steps[i]
-    
+
     return r
 
 
@@ -102,5 +133,5 @@ def _lbfgs_multiply_dense(
             past_grad_diffs[i], past_steps[i])
         right = left.T
         H = left @ H @ right + rho * np.outer(past_steps[i], past_steps[i])
-    
+
     return H @ current_gradient

project_euromir/tests/test_equilibrate.py

@@ -0,0 +1,114 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+import logging
+import time  # on msw timing function calls <1s does not work
+from unittest import TestCase
+
+import numpy as np
+import scipy as sp
+
+import project_euromir as lib
+
+from . import equilibrate
+
+logging.basicConfig(level='INFO')
+
+
+def _make_Q(matrix, b, c):
+    m = len(b)
+    n = len(c)
+    return sp.sparse.bmat([
+        [None, matrix.T, c.reshape(n, 1)],
+        [-matrix, None, b.reshape(m, 1)],
+        [-c.reshape(1, n), -b.reshape(1, m), None]
+    ], format='csc')
+
+
+class TestEquilibrate(TestCase):
+    """Test equilibration functions."""
+
+    def _spectrum_Q_more_stable(
+            self, matrix, b, c, matrix_transf, b_transf, c_transf):
+        """Check that the spectrum of Q after transf has fewer changes."""
+        Q_orig = _make_Q(matrix, b, c)
+        Q_transf = _make_Q(matrix_transf, b_transf, c_transf)
+
+        spectrum_orig = np.linalg.eigh(Q_orig.toarray())[0]
+        spectrum_transf = np.linalg.eigh(Q_transf.toarray())[0]
+
+        diff_spectrum_orig = np.diff(spectrum_orig)
+        diff_spectrum_transf = np.diff(spectrum_transf)
+
+        # import matplotlib.pyplot as plt
+        # plt.plot(spectrum_orig)
+        # plt.plot(spectrum_transf)
+        # plt.show()
+
+        for thres in np.logspace(-8, 0, 9):
+            print('thres', thres)
+            self.assertLessEqual(
+                np.sum(diff_spectrum_orig < thres),
+                np.sum(diff_spectrum_transf < thres),
+            )
+
+    def _generate(self, m, n, density=.01, seed=0):
+        """Generate problem data."""
+        np.random.seed(seed)
+        matrix = sp.sparse.random(
+            m=m, n=n, format='csc', density=density,
+            # this gives very bad unscaled matrices
+            data_rvs=sp.stats.poisson(25, loc=10,).rvs)
+        b = np.random.randn(m)
+        c = np.random.randn(n)
+        return matrix, b, c
+
+    def test_base(self):
+        """Test Python implementation."""
+        m = 50
+        n = 20
+
+        for i, density in enumerate(np.linspace(0.01, 1, 10)):
+            matrix, b, c = self._generate(m, n, density=density, seed=i)
+
+            d, e, sigma, rho, matrix_transf, b_transf, c_transf = \
+                equilibrate.hsde_ruiz_equilibration(
+                    matrix, b, c, dimensions={
+                        'zero': m//2, 'nonneg': m//2, 'second_order': ()})
+            # breakpoint()
+
+            print(np.linalg.norm(c_transf))
+            print(np.linalg.norm(b_transf))
+            print(np.linalg.norm(matrix_transf.toarray(), axis=1))
+            print(np.linalg.norm(matrix_transf.toarray(), axis=0))
+
+            self._spectrum_Q_more_stable(
+                matrix, b, c, matrix_transf, b_transf, c_transf)
+
+        # breakpoint()

project_euromir/tests/test_lbfgs.py

@@ -0,0 +1,42 @@
+# BSD 3-Clause License
+
+# Copyright (c) 2024-, Enzo Busseti
+
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+
+# 1. Redistributions of source code must retain the above copyright notice, this
+#    list of conditions and the following disclaimer.
+
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+#    this list of conditions and the following disclaimer in the documentation
+#    and/or other materials provided with the distribution.
+
+# 3. Neither the name of the copyright holder nor the names of its
+#    contributors may be used to endorse or promote products derived from
+#    this software without specific prior written permission.
+
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+import time  # on msw timing function calls <1s does not work
+from unittest import TestCase
+
+import numpy as np
+import scipy.sparse as sp
+
+import project_euromir as lib
+
+from . import lbfgs
+
+
+class TestLBFGS(TestCase):
+    """Test L-BFGS functions."""

project_euromir/tests/test_linear_algebra.py

@@ -27,7 +27,7 @@
 # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-import time # on msw timing function calls <1s does not work
+import time  # on msw timing function calls <1s does not work
 from unittest import TestCase
 
 import numpy as np