Resolve #13: implemented ARMA process and

- created an example jupytor notebook for the ARMA process
lint yaml to 2 space indentations
- modeling.pade validation tests complete
- more ruff linting updates

.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

.gitignore

@@ -6,3 +6,5 @@ env
 .coverage
 .vscode
 
+.ipynb_checkpoints
+

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"

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

pyssp/modeling.py

@@ -1,61 +1,53 @@
 """Chapter 4 modeling algorithm implementations."""
 
+import logging
 from typing import NoReturn
 
 import numpy as np
-import scipy as sp
 from numpy.typing import ArrayLike
+from scipy import signal
 
 from .state import convm
 
+logger = logging.getLogger(__name__)
+
 
 def pade(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike]:
     """Reference Page 138, Table 4.1.
 
     The Pade approximation models a signal as the unis sample response
     of linear shift invariant system have p poles and q zeros.
     """
-    _x = np.array(x)
+    _x = np.array(x).reshape(-1)
     if p + q > len(_x):
         raise ValueError(f"Model order {p + q} is too large.")
 
     X = convm(_x, p + 1)
 
     # Linear difference matrix spanning the number of zeros
     Xq = X[q + 1:q + p + 1, 1:p + 1].copy()
-    print(Xq.shape)
     a = np.linalg.solve(-Xq, X[q + 1: q + p + 1, 0])
     # a(0) normalized to 1
     a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
     b = X[:q + 1, :p + 1] @ a
 
-    return (a, b)
+    return a.ravel(), b.ravel()
 
 
-def prony(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, ArrayLike]:
+def prony(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, float]:
     """Least square minimization of poles to get denominator coefficients.
 
     Solves directly (Pade method) to get numerator coefficients.
     Also calculates minimum error achieved.
 
     Condition to energy_match is on page 575
     """
-    x = np.array(x)
-    if p + q > len(x):
+    _x = np.array(x).reshape(-1, 1)
+    N = len(_x)
+    if p + q > N:
         raise ValueError(f"Model order {p + q} is too large.")
 
-    # copy and make given signal column array
     X = convm(x, p + 1)
-    print(X.shape)
-    # M = p + q
-    N = len(x)
-    print(f"{N=}")
-    xc = x.copy().reshape(-1, 1)
-
-    # Xq = X[q + 1:q + p + 1, 1:p + 1].copy()
-    # a = np.linalg.solve(-Xq, X[q + 1: q + p + 1, 0])
-    # a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
-    # b = X[:q + 1, :p + 1] @ a
 
     # the factorization does not guarantee nonsingularity!
     # resulting matrix is positive *semi*-definite: all zeros are
@@ -66,76 +58,62 @@ def prony(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, ArrayLike
     rx = Xq_H @ Xq1
     Xinv = np.linalg.inv(Xq_H @ Xq)
     a = -Xinv @ rx
-    print(a.shape)
     # a(0) normalized to 1
     a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
     # same as Pade method
     b = X[:q + 1, :p + 1] @ a
 
     # the minimum squared error
-    err = np.transpose(xc[q + 1:N]) @ X[q + 1:N, :p + 1] @ a
+    err = np.transpose(_x[q + 1:N]) @ X[q + 1:N, :p + 1] @ a
 
-    return a, b, err
+    return a.ravel(), b.ravel(), err.ravel()[0]
 
 
 def shanks(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, float]:
     """Shank's method."""
-    x = np.array(x)
-    N = len(x)
+    _x = np.array(x).ravel().reshape(-1, 1)
+    N = len(_x)
     if p + q >= N:
         raise ValueError(f"Model order {p + q} is too large.")
 
     a, _, _ = prony(x, p, q)
-    a = np.array(a)
-    print(f"{a.transpose().ravel()=}")
+    logger.info(f"{a=}")
     u = np.concatenate((np.ones(1), np.zeros(N - 1)))
-    res = sp.signal.tf2zpk([1], a.ravel())
-    res = sp.signal.zpk2sos(*res)
-    res = sp.signal.sosfilt(res, x=u)
-    G = convm(res.ravel(), q + 1)
-    G0 = G[:N,].copy()
-    print(f"{G0.shape=}")
-    G0_H = np.transpose((G0.copy()).conjugate())
-    x0 = (x.copy()).reshape(-1, 1)
-    gx = G0_H @ x0
-    # the factorization does not guarantee nonsingularity!
-    # resulting matrix is positive *semi*-definite
-    Ginv = np.linalg.inv(G0_H @ G0)
-    print(f"{x.shape=}")
-    print(f"{Ginv.shape=}")
-    b = Ginv @ gx
-    err = 1
+    sos = signal.tf2sos([1], a)
+    g = signal.sosfilt(sos, x=u)
+    logger.info(f"{g=}")
+    G = convm(g, q + 1)
+    G0 = G[:N].copy()
+    logger.info(f"{G0=}")
+    b = np.linalg.lstsq(G0, _x, rcond=None)[0]
+    err = _x.T @ _x - _x.T @ G[:N, :q + 1] @ b
 
-    return a, b, err
+    return a, b.ravel(), err.ravel()[0]
 
 
-def spike(g: ArrayLike, n0: int, n: int) -> ArrayLike:
+def spike(g: ArrayLike, n0: int, n: int) -> tuple[ArrayLike, float]:
     """Leaset Squares Inverse Filter."""
-    g = np.array(g).reshape(-1, 1)
-    m = len(g)
+    _g = np.array(g).reshape(-1, 1)
+    m = len(_g)
 
     if m + n - 1 <= n0:
         raise ValueError(f"m + n - 1 must be less than {n0=}")
 
     G = convm(g, n)
     d = np.zeros((m + n - 1, 1))
     d[n0] = 1
-    print(f"{d.shape=}, {G.shape=}")
-    G_H = np.transpose(G.conjugate())
+    h = np.linalg.lstsq(G, d, rcond=None)[0]
+    err = 1 - G[n0,] @ h
 
-    print(f"{G_H.shape=}, {G.shape=}")
-    Ginv = np.linalg.inv(G_H @ G)
-    h = Ginv @ G_H @ d
-
-    return h
+    return h.ravel(), err.ravel()[0]
 
 
 def ipf(x: ArrayLike, p: int, q: int, n: None, a: ArrayLike) -> NoReturn:
     """Iterative Pre-Filtering."""
     raise NotImplementedError()
 
 
-def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, np.ndarray]:
+def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, float]:
     """The auto-correlation method."""
     x0 = np.array(x).ravel().reshape(-1, 1)
     N = len(x0)
@@ -151,36 +129,49 @@ def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, np.ndarray]:
     a = np.concatenate((np.ones(1), a1)).reshape(-1, 1)
     err = np.abs(X[:N + p, 0].T @ X @ a)
 
-    return a, err
+    return a, err.ravel()[0]
 
 
 def covm(x: ArrayLike, p: int)-> tuple[ArrayLike, float]:
     """Solve the complete Prony normal equations."""
-    x0 = np.array(x).ravel().reshape(-1, 1)
-    N = len(x0)
-    if p >= len(x0):
+    _x = np.array(x).ravel().reshape(-1, 1)
+    N = len(_x)
+    if p >= N:
         raise ValueError(f"{p=} all-pole model too large")
 
-    X = convm(x0, p + 1)
+    X = convm(_x, p + 1)
     Xq = X[p - 1:N - 1, :p].copy()
-    cx = X[p:N, 0].copy()
-    Xq_H = Xq.copy().conjugate().transpose()
-    print(f"{Xq=}")
-    Xinv = np.linalg.inv(Xq_H @ Xq)
-    a1 = -Xinv @ Xq_H @ cx
-    a = np.concatenate((np.ones(1), a1)).reshape(-1, 1)
-    err = np.abs(cx.transpose() @ X[p:N,] @ a)
-    return a, err
+    Xsol = np.linalg.lstsq(-Xq, X[p:N, 0], rcond=None)[0]
+    logger.warning(f"{Xsol=}")
+    a = np.hstack(([1], Xsol))
+    err = np.abs(X[p:N,0] @ X[p:N,] @ a)
+    return a, err.ravel()[0]
 
 
 def durbin(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike]:
     """The durbin method."""
-    x0 = np.array(x).ravel().reshape(-1, 1)
-    # N = len(x0)
-    if p >= len(x0):
+    _x = np.array(x).ravel().reshape(-1, 1)
+    N = len(_x)
+    if p >= N:
         raise ValueError("p (all-pole model) too large")
 
-    a, eps = acm(x, p)
+    a, eps = acm(_x, p)
     b, eps = acm(a / np.sqrt(eps), q)
-    b /= np.sqrt(eps)
+    b = b / np.sqrt(eps)
     return a, b
+
+
+def mywe(x: ArrayLike, p: int, q: int) -> NoReturn:
+    """Modified Yuler-Walker Systems of Equations.
+
+    Page 190.
+    """
+    raise NotImplementedError()
+
+
+def eywe(x: ArrayLike, p: int, q: int) -> NoReturn:
+    """Extended Yuler-Walker Systems of Equations.
+
+    Page 193.
+    """
+    raise NotImplementedError()

pyssp/optimal.py

@@ -1,12 +1,13 @@
-"""Optimal Wiener Filters."""
+"""Optimal Filters, Chapter 7."""
 
+from typing import NoReturn
 
 import numpy as np
-import numpy.typing as npt
+from numpy.typing import ArrayLike
 
 
-def kalman(y: list[float], A: npt.ArrayLike, C: npt.ArrayLike, sigmaw: list[float],
-           sigmav: list[float]) -> tuple[npt.NDArray, ...]:
+def kalman(y: list[float], A: ArrayLike, C: ArrayLike, sigmaw: list[float],
+           sigmav: list[float]) -> tuple[np.ndarray, ...]:
     """Kalman Filter.
 
     y: vector of observations N x q,  n time steps, q sensors
@@ -73,3 +74,8 @@ def wiener_denoise() -> None:
 def wiener_systemid() -> None:
     """Systemid based on FIR wiener filters."""
     raise NotImplementedError()
+
+
+def wiener_hopf(x: ArrayLike, p: int, q: int) -> NoReturn:
+    """Wiener-Hopf Systems of Equations."""
+    raise NotImplementedError()

pyssp/spectrum.py

@@ -1,4 +1,4 @@
-"""Power spectrum and Frequency Estimators."""
+"""Power spectrum and Frequency Estimators. Chapter 8."""
 
 from typing import NoReturn