Resolve #18: Levinson (chapter 5) implementation coverage testing

removed prints from implementations and tests, added logging instead
mahdiolfat · Dec 11, 2023 · f6a6ebd · f6a6ebd
1 parent 3b9d8a1
commit f6a6ebd
Show file tree

Hide file tree

Showing 16 changed files with 579 additions and 195 deletions.
diff --git a/.github/workflows/ContinuousTesting.yml b/.github/workflows/ContinuousTesting.yml
@@ -3,30 +3,30 @@ name: Python Continuous Integration
 on: [push]
 
 jobs:
-    build:
-        runs-on: ubuntu-latest
-        steps:
-            - uses: actions/checkout@v4
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
 
-            - name: Set up Python
-              uses: actions/setup-python@v4
-              with:
-                python-version: '3.12'
-                cache: 'pip'
+      - name: Set up Python
+        uses: actions/setup-python@v4
+        with:
+          python-version: '3.12'
+          cache: 'pip'
 
-            - name: Install Dependencies
-              run: |
-                python -m pip install --upgrade pip setuptools wheel
-                pip install -r requirements.txt
+      - name: Install Dependencies
+        run: |
+          python -m pip install --upgrade pip setuptools wheel
+          pip install -r requirements.txt
 
-            - name: Static Lint
-              run: |
-                ruff . --verbose --output-format=github 
+      - name: Static Lint
+        run: |
+          ruff . --verbose --output-format=github 
 
-            - name: Static Type Check
-              run: |
-                mypy --verbose
+      - name: Static Type Check
+        run: |
+          mypy --verbose
 
-            - name: Test
-              run: |
-                pytest --verbose
+      - name: Test
+        run: |
+          pytest --verbose
diff --git a/.gitignore b/.gitignore
@@ -6,3 +6,5 @@ env
 .coverage
 .vscode
 
+.ipynb_checkpoints
+
diff --git a/examples/process_arma.ipynb b/examples/process_arma.ipynb
diff --git a/pyproject.toml b/pyproject.toml
@@ -10,6 +10,7 @@ license = "LICENSE.txt"
 
 [tool.pytest.ini_options]
 pythonpath = "."
+log_cli = true
 addopts = [
     "--import-mode=importlib",
     "--cov=pyssp"
@@ -23,6 +24,9 @@ indent-width = 4
 target-version = "py312"
 extend-exclude = [".pyenv*"]
 
+[tool.ruff.lint.pydocstyle]
+convention = "google"
+
 [tool.ruff.lint]
 select = ["E",
           "F",
@@ -35,12 +39,10 @@ select = ["E",
           "NPY",
           "PERF",
           "C90",
-          "RUF"]
+          "RUF",
+          "D417", "D414"]
 ignore = ["D213", "D401", "D211"]
 
-[tool.ruff.lint.pydocstyle]
-convention = "google"
-
 [tool.ruff.format]
 quote-style = "double"
 indent-style = "space"

diff --git a/pyssp/__init__.py b/pyssp/__init__.py
@@ -1,10 +1 @@
 """The import root for the book sections."""
-
-# import system as Chapter3
-# import modeling as Chapter4
-# import levinson as Chapter5
-# import lattice as Chapter6
-# import optimal as Chapter7
-# import spectrum as Chapter8
-# import adaptive as Chapter9
-# import state as Appendix
diff --git a/pyssp/lattice.py b/pyssp/lattice.py
@@ -1,11 +1,14 @@
 """Implementation of algorithm from Chapter 6."""
 
+import logging
 from typing import NoReturn
 
 import numpy as np
 import scipy as sp
 from numpy.typing import ArrayLike
 
+logger = logging.getLogger(__name__)
+
 
 def fcov(x: ArrayLike, p: int) -> tuple[ArrayLike, ArrayLike]:
     """Figure 6.15, Page 310.
@@ -26,19 +29,18 @@ def fcov(x: ArrayLike, p: int) -> tuple[ArrayLike, ArrayLike]:
     err = np.empty((p, 1))
 
     for j in range(p):
-        print(j)
+        logger.info(j)
         N = N - 1
-        # print(f"{eplus=}, {eplus.shape=}")
-        # print(f"{eminus=}, {eminus.shape=}")
+        logger.info(f"{eplus=}, {eplus.shape=}")
+        logger.info(f"{eminus=}, {eminus.shape=}")
         gamma[j] = (np.transpose(-eminus) @ eplus) / (np.transpose(eminus) @ eminus)
         temp1 = eplus + gamma[j] * eminus
         temp2 = eminus + np.conjugate(gamma[j]) * eplus
         err[j] = np.transpose(temp1) @ temp1
         eplus = temp1[1:N]
         eminus = temp2[:N - 1]
-        print(gamma)
-        print(err)
-        print()
+        logger.info(gamma)
+        logger.info(err)
 
     return gamma, err
 
@@ -60,10 +62,10 @@ def burg(x: ArrayLike, p: int) -> tuple[ArrayLike, ArrayLike]:
     err = np.empty((p, 1))
 
     for j in range(p):
-        print(j)
+        logger.info(j)
         N = N - 1
-        # print(f"{eplus=}, {eplus.shape=}")
-        # print(f"{eminus=}, {eminus.shape=}")
+        logger.info(f"{eplus=}, {eplus.shape=}")
+        logger.info(f"{eminus=}, {eminus.shape=}")
         eplusmag = np.transpose(eplus) @ eplus
         eminusmag = np.transpose(eplus) @ eplus
         gamma[j] = (np.transpose(-2 * eminus) @ eplus) / (eplusmag + eminusmag)
@@ -72,7 +74,6 @@ def burg(x: ArrayLike, p: int) -> tuple[ArrayLike, ArrayLike]:
         err[j] = np.transpose(temp1) @ temp1 + np.transpose(temp2) @ temp2
         eplus = temp1[1:N]
         eminus = temp2[:N - 1]
-        print()
 
     return gamma, err
 
@@ -105,7 +106,6 @@ def mcov(x: ArrayLike, p: int) -> tuple[ArrayLike, ArrayLike]:
     b = b1 + b2
     a = sp.linalg.solve_toeplitz(Rx[:, 1], b)
     a = np.concatenate(([1], a))
-    # print(a.shape)
     err = np.dot(R[0], a) + np.dot(np.flip(R[p]), a)
 
     return a, err
diff --git a/pyssp/levinson.py b/pyssp/levinson.py
@@ -1,9 +1,13 @@
 """Chapter 5 algorithm implementations."""
 
+import logging
+from typing import NoReturn
 
 import numpy as np
 from numpy.typing import ArrayLike
 
+logger = logging.getLogger()
+
 
 def rtoa(r: ArrayLike) -> tuple[np.ndarray, float]:
     """Recursively map a set of autocorrelations to a set of model parameters.
@@ -22,8 +26,8 @@ def rtoa(r: ArrayLike) -> tuple[np.ndarray, float]:
         anT = np.conjugate(np.flipud(a))
         a = an + gamma * np.concatenate((np.zeros(1), anT.ravel())).reshape(-1, 1)
         epsilon = epsilon * (1 - np.abs(gamma)**2)
-        print(f"{gamma=},\n{a=},\n{epsilon=}\n")
-    return a, epsilon
+        logger.info(f"{gamma=},\n{a=},\n{epsilon=}\n")
+    return a.ravel(), epsilon.ravel()[0]
 
 
 def gtoa(gamma: ArrayLike) -> ArrayLike:
@@ -40,13 +44,13 @@ def gtoa(gamma: ArrayLike) -> ArrayLike:
     a = np.ones((1, 1))
     _g = np.array(gamma)
     p = len(_g)
-    for j in range(1, p):
+    for j in range(1, p + 1):
         a = np.concatenate((a, np.zeros((1, 1))))
         _a = a.copy()
         af = np.conjugate(np.flipud(_a))
-        a = a + _g[j] * af
+        a = a + _g[j - 1] * af
 
-    return a
+    return a.ravel()
 
 
 def atog(a: ArrayLike) -> ArrayLike:
@@ -70,23 +74,25 @@ def atog(a: ArrayLike) -> ArrayLike:
     gamma[p - 2] = _a[p - 2]
 
     for j in range(p - 2, 0, -1):
-        # print(f"{gamma=}, {_a=}")
+        # logger.info(f"{gamma=}, {_a=}")
         ai1 = _a[:j].copy()
         ai2 = _a[:j].copy()
         af = np.flipud(np.conjugate(ai1))
-        # print(f"{ai1=}, {ai2=}, {af=}")
+        # logger.info(f"{ai1=}, {ai2=}, {af=}")
         s1 = ai2 - gamma[j] * af
         s2 = 1 - np.abs(gamma[j])**2
         _a = np.divide(s1, s2)
-        # print(f"{s1=}, {s2=}, {_a=}")
+        # logger.info(f"{s1=}, {s2=}, {_a=}")
         gamma[j - 1] = _a[j - 1]
 
-    return gamma
+    return gamma.ravel()
 
 
 def gtor(gamma: ArrayLike, epsilon: None | float = None) -> ArrayLike:
     """Find the autocorrelation sequence from the reflection coefficients and the modeling error.
 
+    Also called the Inverse Levinson-Durbin Recursion.
+
     Page 241, Figure 5.9.
     """
     _g = np.array(gamma)
@@ -99,19 +105,18 @@ def gtor(gamma: ArrayLike, epsilon: None | float = None) -> ArrayLike:
         aa0 = np.concatenate((aa, np.zeros((1, 1)))).reshape(-1, 1)
         aaf = np.conjugate(np.flipud(aa1))
         aa = aa0 + _g[j] * aaf
-        print(aa)
+        logger.info(aa)
         rf = -np.fliplr(r) @ aa
-        print(rf)
-        print(rf.shape)
-        print(r)
-        print(r.shape)
-        print()
+        logger.info(rf)
+        logger.info(rf.shape)
+        logger.info(r)
+        logger.info(r.shape)
         r = np.concatenate((r[0], rf[0])).reshape(1, -1)
 
     if epsilon is not None:
         r = r * epsilon / np.prod(1 - np.abs(gamma)**2)
 
-    return r
+    return r.ravel()
 
 
 def ator(a: ArrayLike, b: float) -> ArrayLike:
@@ -147,27 +152,40 @@ def glev(r: ArrayLike, b: ArrayLike) -> ArrayLike:
     x = np.array([_b[0] / _r[0]]).reshape(-1, 1)
     epsilon = _r[0]
     for j in range(1, p):
-        print(j)
-        print(f"{_r=}, {_r.shape=}")
+        logger.info(j)
+        logger.info(f"{_r=}, {_r.shape=}")
         _r1 = np.transpose(np.array(_r[1:j + 1]))
-        print(f"{_r1=}, {_r1.shape=}")
-        print(f"{x=}, {x.shape=}")
-        print(f"{a=}, {a.shape=}")
+        logger.info(f"{_r1=}, {_r1.shape=}")
+        logger.info(f"{x=}, {x.shape=}")
+        logger.info(f"{a=}, {a.shape=}")
         g = _r1 @ np.flipud(a)
-        print(f"{g=}, {g.shape=}")
-        gamma = -g / epsilon
-        print(f"{gamma=}, {gamma.shape=}")
+        logger.info(f"{g=}, {g.shape=}")
+        gamma = g * -1 / epsilon
+        logger.info(f"{gamma=}, {gamma.shape=}")
         _a0 = np.concatenate([a, [[0]]])
         _af = np.conjugate(np.flipud(_a0))
-        print(f"{_a0=}, {_a0.shape=}")
-        print(f"{_af=}, {_af.shape=}")
+        logger.info(f"{_a0=}, {_a0.shape=}")
+        logger.info(f"{_af=}, {_af.shape=}")
         a = _a0 + gamma * _af
         epsilon = epsilon * (1 - np.abs(gamma)**2)
-        print(f"{epsilon=}")
+        logger.info(f"{epsilon=}")
         delta = _r1 @ np.flipud(x)
         q = (_b[j] - delta[0, 0]) / epsilon
         x = np.concatenate([x, [[0]]])
         x = x + q * np.conjugate(np.flipud(a))
-        print()
 
-    return x
+
+    return x.ravel()
+
+
+def shur_cohn_test() -> NoReturn:
+    """Check stability of any rational filter."""
+    raise NotImplementedError()
+
+
+def splitlev(r: ArrayLike, b: ArrayLike) -> ArrayLike:
+    """Implement the split levinson recursion.
+
+    Table 5.7, Page 274.
+    """
+    raise NotImplementedError()
Original file line number	Diff line number	Diff line change
Expand Up		@@ -6,3 +6,5 @@ env
		.coverage
		.vscode

		.ipynb_checkpoints