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 diff --git a/experiment/power_iteration.py b/experiment/power_iteration.py new file mode 100644 index 0000000..96e67f5 --- /dev/null +++ b/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) 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 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/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: 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