Merge branch 'main' of https://github.com/luk036/ellalgo

.gitignore

@@ -173,3 +173,6 @@ MANIFEST
 .venv*/
 .conda*/
 .python-version
+
+# my stuffs
+.hypothesis

experiment/power_iteration.py

@@ -0,0 +1,121 @@
+import numpy as np
+
+"""
+Power method for finding the largest eigenvalue of a square matrix
+"""
+
+
+class Options:
+    def __init__(self, max_iters, tolerance):
+        self.max_iters = max_iters
+        self.tolerance = tolerance
+
+
+def norm_l1(x):
+    return np.sum(np.abs(x))
+
+
+def power_iteration(A, x, options):
+    """Power iteration method
+
+    x: assuming not zero.
+    """
+    x /= np.sqrt(np.sum(x**2))
+    for niter in range(options.max_iters):
+        x1 = x
+        x = A @ x1
+        x /= np.sqrt(np.sum(x**2))
+        if norm_l1(x - x1) <= options.tolerance or norm_l1(x + x1) <= options.tolerance:
+            return x, x @ (A @ x), niter
+    return x, x @ (A @ x), options.max_iters
+
+
+def power_iteration4(A, x, options):
+    """Power iteration method
+
+    x: assuming not zero.
+    """
+    x /= norm_l1(x)
+    for niter in range(options.max_iters):
+        x1 = x
+        x = A @ x1
+        x /= norm_l1(x)
+        if norm_l1(x - x1) <= options.tolerance or norm_l1(x + x1) <= options.tolerance:
+            x /= np.sqrt(np.sum(x**2))
+            return x, x @ (A @ x), niter
+    x /= np.sqrt(np.sum(x**2))
+    return x, x @ (A @ x), options.max_iters
+
+
+def power_iteration2(A, x, options):
+    """Power iteration method
+
+    x: assuming not zero.
+    """
+    x /= np.sqrt(np.sum(x**2))
+    new = A @ x
+    ld = x @ new
+    for niter in range(options.max_iters):
+        ld1 = ld
+        x[:] = new[:]
+        x /= np.sqrt(np.sum(x**2))
+        new = A @ x
+        ld = x @ new
+        if abs(ld1 - ld) <= options.tolerance:
+            return x, ld, niter
+    return x, ld, options.max_iters
+
+
+def power_iteration3(A, x, options):
+    """Power iteration method
+
+    x: assuming not zero.
+    """
+    new = A @ x
+    dot = x @ x
+    ld = (x @ new) / dot
+    for niter in range(options.max_iters):
+        ld1 = ld
+        x[:] = new[:]
+        dot = x @ x
+        if dot >= 1e150:
+            x /= np.sqrt(np.sum(x**2))
+            new = A @ x
+            ld = x @ new
+            if abs(ld1 - ld) <= options.tolerance:
+                return x, ld, niter
+        else:
+            new = A @ x
+            ld = (x @ new) / dot
+            if abs(ld1 - ld) <= options.tolerance:
+                x /= np.sqrt(np.sum(x**2))
+                return x, ld, niter
+    x /= np.sqrt(np.sum(x**2))
+    return x, ld, options.max_iters
+
+
+# Test data
+A = np.array([[3.7, -3.6, 0.7], [-3.6, 4.3, -2.8], [0.7, -2.8, 5.4]])
+options = Options(max_iters=2000, tolerance=1e-7)
+
+x = np.array([0.3, 0.5, 0.4])
+x1, ld, niter = power_iteration(A, x, options)
+print(x1)
+print(ld)
+
+x = np.array([0.3, 0.5, 0.4])
+x4, ld, niter = power_iteration4(A, x, options)
+print(x4)
+print(ld)
+
+options.tolerance = 1e-14
+
+x = np.array([0.3, 0.5, 0.4])
+x2, ld, niter = power_iteration2(A, x, options)
+print(x2)
+print(ld)
+
+x = np.array([0.3, 0.5, 0.4])
+x3, ld, niter = power_iteration3(A, x, options)
+print(x3)
+print(ld)

oryx-build-commands.txt

@@ -1,2 +1,2 @@
-PlatformWithVersion=Python
+PlatformWithVersion=Python 
 BuildCommands=conda env create --file environment.yml --prefix ./venv --quiet

setup.py

@@ -1,10 +1,10 @@
 """
-    Setup file for ellalgo.
-    Use setup.cfg to configure your project.
+Setup file for ellalgo.
+Use setup.cfg to configure your project.
 
-    This file was generated with PyScaffold 4.5.
-    PyScaffold helps you to put up the scaffold of your new Python project.
-    Learn more under: https://pyscaffold.org/
+This file was generated with PyScaffold 4.5.
+PyScaffold helps you to put up the scaffold of your new Python project.
+Learn more under: https://pyscaffold.org/
 """
 
 from setuptools import setup

src/ellalgo/oracles/lmi_oracle.py

@@ -29,7 +29,7 @@ def __init__(self, mat_f, mat_b):
         self.mat_f0 = mat_b
         self.ldlt_mgr = LDLTMgr(len(mat_b))
 
-    def assess_feas(self, x: np.ndarray) -> Optional[Cut]:
+    def assess_feas(self, xc: np.ndarray) -> Optional[Cut]:
         """
         The `assess_feas` function assesses the feasibility of a given input and returns a cut if it is not
         feasible.
@@ -41,7 +41,7 @@ def assess_feas(self, x: np.ndarray) -> Optional[Cut]:
 
         def get_elem(i, j):
             return self.mat_f0[i, j] - sum(
-                Fk[i, j] * xk for Fk, xk in zip(self.mat_f, x)
+                Fk[i, j] * xk for Fk, xk in zip(self.mat_f, xc)
             )
 
         if self.ldlt_mgr.factor(get_elem):

src/ellalgo/oracles/lowpass_oracle.py

@@ -39,7 +39,7 @@
 # *********************************************************************
 # number of FIR coefficients (including zeroth)
 class LowpassOracle(OracleOptim):
-    more_alt: bool = True
+    # more_alt: bool = True
     idx1: int = 0
 
     def __init__(
@@ -83,7 +83,7 @@ def assess_feas(self, x: Arr) -> Optional[ParallelCut]:
         Returns:
             [type]: [description]
         """
-        self.more_alt = True
+        # self.more_alt = True
 
         mdim, ndim = self.spectrum.shape
         for _ in range(self.nwpass):
@@ -126,7 +126,7 @@ def assess_feas(self, x: Arr) -> Optional[ParallelCut]:
             if v < 0:
                 return -col_k, -v  # single cut
 
-        self.more_alt = False
+        # self.more_alt = False
 
         # case 1 (unlikely)
         if x[0] < 0:

tests/conftest.py

@@ -1,10 +1,10 @@
 """
-    Dummy conftest.py for ellalgo.
+Dummy conftest.py for ellalgo.
 
-    If you don't know what this is for, just leave it empty.
-    Read more about conftest.py under:
-    - https://docs.pytest.org/en/stable/fixture.html
-    - https://docs.pytest.org/en/stable/writing_plugins.html
+If you don't know what this is for, just leave it empty.
+Read more about conftest.py under:
+- https://docs.pytest.org/en/stable/fixture.html
+- https://docs.pytest.org/en/stable/writing_plugins.html
 """
 
 # import pytest