From 54abb07f2d8273ecfd953cb0bbd3a2336a4fe5cb Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Wed, 22 May 2024 18:43:02 +0800 Subject: [PATCH 1/8] use latest .gitignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index d4a6f0e..7869e21 100644 --- a/.gitignore +++ b/.gitignore @@ -173,3 +173,6 @@ MANIFEST .venv*/ .conda*/ .python-version + +# my stuffs +.hypothesis From cce729fde6fca5030904a85105dc9fcc98008b86 Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Thu, 6 Jun 2024 02:18:18 +0000 Subject: [PATCH 2/8] apply black --- src/ellalgo/oracles/ldlt_mgr.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ellalgo/oracles/ldlt_mgr.py b/src/ellalgo/oracles/ldlt_mgr.py index bf4350a..3e52750 100644 --- a/src/ellalgo/oracles/ldlt_mgr.py +++ b/src/ellalgo/oracles/ldlt_mgr.py @@ -33,7 +33,6 @@ class LDLTMgr: __slots__ = ("pos", "wit", "_ndim", "_storage") - def __init__(self, ndim: int): """ The above function is the constructor for a LDLT Ext object, which initializes various attributes From 4beb219e26120a640e423fd6f46702279afffba0 Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Thu, 20 Jun 2024 04:16:35 +0000 Subject: [PATCH 3/8] apply ruff format --- setup.py | 10 +++++----- tests/conftest.py | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/setup.py b/setup.py index 4122877..cad0c37 100644 --- a/setup.py +++ b/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 diff --git a/tests/conftest.py b/tests/conftest.py index f68cb2b..22f7022 100644 --- a/tests/conftest.py +++ b/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 From 5cd4889f5cba57f2273de1488652e76a37d9405f Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Tue, 9 Jul 2024 08:59:32 +0000 Subject: [PATCH 4/8] add power_iteration --- experiment/power_iteration.py | 126 ++++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 experiment/power_iteration.py diff --git a/experiment/power_iteration.py b/experiment/power_iteration.py new file mode 100644 index 0000000..56b8ed7 --- /dev/null +++ b/experiment/power_iteration.py @@ -0,0 +1,126 @@ +import numpy as np +import math + + +class Options: + max_iters = 2000 + tolerance = 1e-9 + + +def power_iteration(A: np.ndarray, x: np.ndarray, options): + """Power iteration method""" + x = x / math.sqrt(x @ x) + for niter in range(options.max_iters): + x1 = A @ x + x1 = x1 / math.sqrt(x1 @ x1) + eps = np.linalg.norm(x - x1, np.inf) + if eps <= options.tolerance: + break + eps = np.linalg.norm(x + x1, np.inf) + if eps <= options.tolerance: + break + x = x1 + + ld = x1 @ A @ x1 + return x1, ld, niter, eps + + +def calc_core2(A, x): + x /= math.sqrt(x @ x) + new = A @ x + ld = x @ new + return new, ld + + +def power_iteration2(A: np.ndarray, x: np.ndarray, options): + """Power iteration method""" + new, ld = calc_core2(A, x) + for niter in range(options.max_iters): + ld1 = ld + x = new + new, ld = calc_core2(A, x) + eps = abs(ld1 - ld) + if eps <= options.tolerance: + break + return x, ld, niter, eps + + +def calc_core3(A, x): + new = A @ x + dot = x @ x + ld = (x @ new) / dot + return new, dot, ld + + +def power_iteration3(A: np.ndarray, x: np.ndarray, options): + """Power iteration method""" + + new, dot, ld = calc_core3(A, x) + for niter in range(options.max_iters): + ld1 = ld + x = new + xmax = max(x) + xmin = min(x) + if xmax >= 1e100 or xmin <= -1e100: + x /= 1e100 + new, dot, ld = calc_core3(A, x) + eps = abs(ld1 - ld) + if eps <= options.tolerance: + break + + x /= math.sqrt(dot) + return x, ld, niter, eps + + +def power_iteration4(A: np.ndarray, x: np.ndarray, options): + """Power iteration method""" + x = x / max(np.abs(x)) + for niter in range(options.max_iters): + x1 = A @ x + x1 = x1 / max(np.abs(x1)) + eps = np.linalg.norm(x - x1, np.inf) + if eps <= options.tolerance: + break + eps = np.linalg.norm(x + x1, np.inf) + if eps <= options.tolerance: + break + x = x1 + + x1 = x1 / math.sqrt(x1 @ x1) + ld = x1 @ A @ x1 + return x1, ld, niter, eps + + +if __name__ == "__main__": + A = np.array([[3.7, -4.5], [4.3, -5.9]]) + x, ld, niter, eps = power_iteration(A, np.array([1.1, 0.0]), Options()) + + print("1-----------------------------") + print(x) + print(ld) + print(niter) + print(eps) + + print("2-----------------------------") + + x, ld, niter, eps = power_iteration2(A, np.array([1.1, 0.0]), Options()) + print(x) + print(ld) + print(niter) + print(eps) + + print("3-----------------------------") + + x, ld, niter, eps = power_iteration3(A, np.array([1.1, 0.0]), Options()) + print(x) + print(ld) + print(niter) + print(eps) + + print("4-----------------------------") + + x, ld, niter, eps = power_iteration4(A, np.array([1.1, 0.0]), Options()) + print(x) + print(ld) + print(niter) + print(eps) From 5eb9dcd897ac7548c5f6bda47e49fff0b38d3b8b Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Tue, 9 Jul 2024 18:47:13 +0800 Subject: [PATCH 5/8] minor change --- src/ellalgo/oracles/lmi_oracle.py | 4 ++-- src/ellalgo/oracles/lowpass_oracle.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ellalgo/oracles/lmi_oracle.py b/src/ellalgo/oracles/lmi_oracle.py index bd0b8e4..a3787e6 100644 --- a/src/ellalgo/oracles/lmi_oracle.py +++ b/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): diff --git a/src/ellalgo/oracles/lowpass_oracle.py b/src/ellalgo/oracles/lowpass_oracle.py index fe64d3d..d99a3ff 100644 --- a/src/ellalgo/oracles/lowpass_oracle.py +++ b/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: From 9cda6349e7ec7140de56c9e342947394478fe5dc Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Wed, 10 Jul 2024 16:40:52 +0800 Subject: [PATCH 6/8] modify power_iteration --- experiment/power_iteration.py | 131 ++++++++++++++++------------------ 1 file changed, 63 insertions(+), 68 deletions(-) diff --git a/experiment/power_iteration.py b/experiment/power_iteration.py index 56b8ed7..e34fe20 100644 --- a/experiment/power_iteration.py +++ b/experiment/power_iteration.py @@ -13,114 +13,109 @@ def power_iteration(A: np.ndarray, x: np.ndarray, options): for niter in range(options.max_iters): x1 = A @ x x1 = x1 / math.sqrt(x1 @ x1) - eps = np.linalg.norm(x - x1, np.inf) - if eps <= options.tolerance: - break - eps = np.linalg.norm(x + x1, np.inf) - if eps <= options.tolerance: - break + if ( + sum(np.abs(x - x1)) <= options.tolerance + or sum(np.abs(x + x1)) <= options.tolerance + ): + ld = x1 @ A @ x1 + return x1, ld, niter x = x1 - ld = x1 @ A @ x1 - return x1, ld, niter, eps + ld = x @ A @ x + return x, ld, options.max_iters -def calc_core2(A, x): - x /= math.sqrt(x @ x) - new = A @ x - ld = x @ new - return new, ld +def power_iteration4(A: np.ndarray, x: np.ndarray, options): + """Power iteration method""" + x = x / sum(np.abs(x)) + for niter in range(options.max_iters): + x1 = A @ x + x1 = x1 / sum(np.abs(x1)) + if ( + sum(np.abs(x - x1)) <= options.tolerance + or sum(np.abs(x + x1)) <= options.tolerance + ): + x1 = x1 / math.sqrt(x1 @ x1) + ld = x1 @ A @ x1 + return x1, ld, niter + x = x1 + + x = x / math.sqrt(x @ x) + ld = x @ A @ x + return x, ld, options.max_iters def power_iteration2(A: np.ndarray, x: np.ndarray, options): """Power iteration method""" - new, ld = calc_core2(A, x) + x /= math.sqrt(x @ x) + new = A @ x + ld = x @ new for niter in range(options.max_iters): ld1 = ld x = new - new, ld = calc_core2(A, x) - eps = abs(ld1 - ld) - if eps <= options.tolerance: - break - return x, ld, niter, eps - - -def calc_core3(A, x): - new = A @ x - dot = x @ x - ld = (x @ new) / dot - return new, dot, ld + x /= math.sqrt(x @ x) + 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: np.ndarray, x: np.ndarray, options): """Power iteration method""" - new, dot, ld = calc_core3(A, x) + new = A @ x + dot = x @ x + ld = (x @ new) / dot for niter in range(options.max_iters): ld1 = ld x = new - xmax = max(x) - xmin = min(x) - if xmax >= 1e100 or xmin <= -1e100: - x /= 1e100 - new, dot, ld = calc_core3(A, x) - eps = abs(ld1 - ld) - if eps <= options.tolerance: - break - + dot = x @ x + if dot >= 1e150: + x /= math.sqrt(x @ x) + 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 /= math.sqrt(x @ x) + return x, ld, niter x /= math.sqrt(dot) - return x, ld, niter, eps - - -def power_iteration4(A: np.ndarray, x: np.ndarray, options): - """Power iteration method""" - x = x / max(np.abs(x)) - for niter in range(options.max_iters): - x1 = A @ x - x1 = x1 / max(np.abs(x1)) - eps = np.linalg.norm(x - x1, np.inf) - if eps <= options.tolerance: - break - eps = np.linalg.norm(x + x1, np.inf) - if eps <= options.tolerance: - break - x = x1 - - x1 = x1 / math.sqrt(x1 @ x1) - ld = x1 @ A @ x1 - return x1, ld, niter, eps + return x, ld, options.max_iters if __name__ == "__main__": - A = np.array([[3.7, -4.5], [4.3, -5.9]]) - x, ld, niter, eps = power_iteration(A, np.array([1.1, 0.0]), Options()) + A = np.array([[3.7, -3.6, 0.7], [-3.6, 4.3, -2.8], [0.7, -2.8, 5.4]]) + options = Options() + options.tolerance = 1e-7 + x, ld, niter = power_iteration(A, np.array([0.3, 0.5, 0.4]), options) print("1-----------------------------") print(x) print(ld) print(niter) - print(eps) - print("2-----------------------------") + print("4-----------------------------") - x, ld, niter, eps = power_iteration2(A, np.array([1.1, 0.0]), Options()) + x, ld, niter = power_iteration4(A, np.array([0.3, 0.5, 0.4]), options) print(x) print(ld) print(niter) - print(eps) - print("3-----------------------------") + options.tolerance = 1e-14 + print("2-----------------------------") - x, ld, niter, eps = power_iteration3(A, np.array([1.1, 0.0]), Options()) + x, ld, niter = power_iteration2(A, np.array([0.3, 0.5, 0.4]), options) print(x) print(ld) print(niter) - print(eps) - print("4-----------------------------") + print("3-----------------------------") - x, ld, niter, eps = power_iteration4(A, np.array([1.1, 0.0]), Options()) + x, ld, niter = power_iteration3(A, np.array([0.3, 0.5, 0.4]), options) print(x) print(ld) print(niter) - print(eps) From cdf88bedfebafeba484669563bdc2ba8b2f937fe Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Thu, 11 Jul 2024 11:10:12 +0800 Subject: [PATCH 7/8] better version of power_iteration --- experiment/power_iteration.py | 148 +++++++++++++++++----------------- 1 file changed, 74 insertions(+), 74 deletions(-) diff --git a/experiment/power_iteration.py b/experiment/power_iteration.py index e34fe20..96e67f5 100644 --- a/experiment/power_iteration.py +++ b/experiment/power_iteration.py @@ -1,59 +1,64 @@ import numpy as np -import math + +""" +Power method for finding the largest eigenvalue of a square matrix +""" class Options: - max_iters = 2000 - tolerance = 1e-9 + def __init__(self, max_iters, tolerance): + self.max_iters = max_iters + self.tolerance = tolerance -def power_iteration(A: np.ndarray, x: np.ndarray, options): - """Power iteration method""" - x = x / math.sqrt(x @ x) - for niter in range(options.max_iters): - x1 = A @ x - x1 = x1 / math.sqrt(x1 @ x1) - if ( - sum(np.abs(x - x1)) <= options.tolerance - or sum(np.abs(x + x1)) <= options.tolerance - ): - ld = x1 @ A @ x1 - return x1, ld, niter - x = x1 - - ld = x @ A @ x - return x, ld, options.max_iters +def norm_l1(x): + return np.sum(np.abs(x)) -def power_iteration4(A: np.ndarray, x: np.ndarray, options): - """Power iteration method""" - x = x / 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 = A @ x - x1 = x1 / sum(np.abs(x1)) - if ( - sum(np.abs(x - x1)) <= options.tolerance - or sum(np.abs(x + x1)) <= options.tolerance - ): - x1 = x1 / math.sqrt(x1 @ x1) - ld = x1 @ A @ x1 - return x1, ld, niter - x = x1 - - x = x / math.sqrt(x @ x) - ld = x @ A @ x - return x, ld, 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_iteration2(A: np.ndarray, x: np.ndarray, options): - """Power iteration method""" - x /= math.sqrt(x @ x) +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 /= math.sqrt(x @ x) + x[:] = new[:] + x /= np.sqrt(np.sum(x**2)) new = A @ x ld = x @ new if abs(ld1 - ld) <= options.tolerance: @@ -61,18 +66,20 @@ def power_iteration2(A: np.ndarray, x: np.ndarray, options): return x, ld, options.max_iters -def power_iteration3(A: np.ndarray, x: np.ndarray, options): - """Power iteration method""" +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 + x[:] = new[:] dot = x @ x if dot >= 1e150: - x /= math.sqrt(x @ x) + x /= np.sqrt(np.sum(x**2)) new = A @ x ld = x @ new if abs(ld1 - ld) <= options.tolerance: @@ -81,41 +88,34 @@ def power_iteration3(A: np.ndarray, x: np.ndarray, options): new = A @ x ld = (x @ new) / dot if abs(ld1 - ld) <= options.tolerance: - x /= math.sqrt(x @ x) + x /= np.sqrt(np.sum(x**2)) return x, ld, niter - x /= math.sqrt(dot) + x /= np.sqrt(np.sum(x**2)) return x, ld, options.max_iters -if __name__ == "__main__": - A = np.array([[3.7, -3.6, 0.7], [-3.6, 4.3, -2.8], [0.7, -2.8, 5.4]]) - options = Options() - options.tolerance = 1e-7 - x, ld, niter = power_iteration(A, np.array([0.3, 0.5, 0.4]), options) - - print("1-----------------------------") - print(x) - print(ld) - print(niter) - - print("4-----------------------------") +# 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, ld, niter = power_iteration4(A, np.array([0.3, 0.5, 0.4]), options) - print(x) - print(ld) - print(niter) +x = np.array([0.3, 0.5, 0.4]) +x1, ld, niter = power_iteration(A, x, options) +print(x1) +print(ld) - options.tolerance = 1e-14 - print("2-----------------------------") +x = np.array([0.3, 0.5, 0.4]) +x4, ld, niter = power_iteration4(A, x, options) +print(x4) +print(ld) - x, ld, niter = power_iteration2(A, np.array([0.3, 0.5, 0.4]), options) - print(x) - print(ld) - print(niter) +options.tolerance = 1e-14 - print("3-----------------------------") +x = np.array([0.3, 0.5, 0.4]) +x2, ld, niter = power_iteration2(A, x, options) +print(x2) +print(ld) - x, ld, niter = power_iteration3(A, np.array([0.3, 0.5, 0.4]), options) - print(x) - print(ld) - print(niter) +x = np.array([0.3, 0.5, 0.4]) +x3, ld, niter = power_iteration3(A, x, options) +print(x3) +print(ld) From afd97da790bba2156e08435ec8a802ab82343af4 Mon Sep 17 00:00:00 2001 From: Wai-Shing Luk Date: Thu, 18 Jul 2024 14:38:29 +0000 Subject: [PATCH 8/8] delete whitespace --- oryx-build-commands.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/oryx-build-commands.txt b/oryx-build-commands.txt index e7abcde..d647bdf 100644 --- a/oryx-build-commands.txt +++ b/oryx-build-commands.txt @@ -1,2 +1,2 @@ -PlatformWithVersion=Python +PlatformWithVersion=Python BuildCommands=conda env create --file environment.yml --prefix ./venv --quiet