Skip to content

Commit 9c67144

Browse files
mahdiolfatmahdiolfat
committed
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
1 parent 486816b commit 9c67144

File tree

12 files changed

+440
-143
lines changed

12 files changed

+440
-143
lines changed

‎.github/workflows/ContinuousTesting.yml

Lines changed: 22 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -3,30 +3,30 @@ name: Python Continuous Integration
33
on: [push]
44

55
jobs:
6-
build:
7-
runs-on: ubuntu-latest
8-
steps:
9-
- uses: actions/checkout@v4
6+
build:
7+
runs-on: ubuntu-latest
8+
steps:
9+
- uses: actions/checkout@v4
1010

11-
- name: Set up Python
12-
uses: actions/setup-python@v4
13-
with:
14-
python-version: '3.12'
15-
cache: 'pip'
11+
- name: Set up Python
12+
uses: actions/setup-python@v4
13+
with:
14+
python-version: '3.12'
15+
cache: 'pip'
1616

17-
- name: Install Dependencies
18-
run: |
19-
python -m pip install --upgrade pip setuptools wheel
20-
pip install -r requirements.txt
17+
- name: Install Dependencies
18+
run: |
19+
python -m pip install --upgrade pip setuptools wheel
20+
pip install -r requirements.txt
2121
22-
- name: Static Lint
23-
run: |
24-
ruff . --verbose --output-format=github
22+
- name: Static Lint
23+
run: |
24+
ruff . --verbose --output-format=github
2525
26-
- name: Static Type Check
27-
run: |
28-
mypy --verbose
26+
- name: Static Type Check
27+
run: |
28+
mypy --verbose
2929
30-
- name: Test
31-
run: |
32-
pytest --verbose
30+
- name: Test
31+
run: |
32+
pytest --verbose

‎.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,5 @@ env
66
.coverage
77
.vscode
88

9+
.ipynb_checkpoints
10+

‎examples/process_arma.ipynb

Lines changed: 187 additions & 0 deletions
Large diffs are not rendered by default.

‎pyproject.toml

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ license = "LICENSE.txt"
1010

1111
[tool.pytest.ini_options]
1212
pythonpath = "."
13+
log_cli = true
1314
addopts = [
1415
"--import-mode=importlib",
1516
"--cov=pyssp"
@@ -23,6 +24,9 @@ indent-width = 4
2324
target-version = "py312"
2425
extend-exclude = [".pyenv*"]
2526

27+
[tool.ruff.lint.pydocstyle]
28+
convention = "google"
29+
2630
[tool.ruff.lint]
2731
select = ["E",
2832
"F",
@@ -35,12 +39,10 @@ select = ["E",
3539
"NPY",
3640
"PERF",
3741
"C90",
38-
"RUF"]
42+
"RUF",
43+
"D417", "D414"]
3944
ignore = ["D213", "D401", "D211"]
4045

41-
[tool.ruff.lint.pydocstyle]
42-
convention = "google"
43-
4446
[tool.ruff.format]
4547
quote-style = "double"
4648
indent-style = "space"

‎pyssp/__init__.py

Lines changed: 0 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1 @@
11
"""The import root for the book sections."""
2-
3-
# import system as Chapter3
4-
# import modeling as Chapter4
5-
# import levinson as Chapter5
6-
# import lattice as Chapter6
7-
# import optimal as Chapter7
8-
# import spectrum as Chapter8
9-
# import adaptive as Chapter9
10-
# import state as Appendix

‎pyssp/modeling.py

Lines changed: 62 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -1,61 +1,53 @@
11
"""Chapter 4 modeling algorithm implementations."""
22

3+
import logging
34
from typing import NoReturn
45

56
import numpy as np
6-
import scipy as sp
77
from numpy.typing import ArrayLike
8+
from scipy import signal
89

910
from .state import convm
1011

12+
logger = logging.getLogger(__name__)
13+
1114

1215
def pade(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike]:
1316
"""Reference Page 138, Table 4.1.
1417
1518
The Pade approximation models a signal as the unis sample response
1619
of linear shift invariant system have p poles and q zeros.
1720
"""
18-
_x = np.array(x)
21+
_x = np.array(x).reshape(-1)
1922
if p + q > len(_x):
2023
raise ValueError(f"Model order {p + q} is too large.")
2124

2225
X = convm(_x, p + 1)
2326

2427
# Linear difference matrix spanning the number of zeros
2528
Xq = X[q + 1:q + p + 1, 1:p + 1].copy()
26-
print(Xq.shape)
2729
a = np.linalg.solve(-Xq, X[q + 1: q + p + 1, 0])
2830
# a(0) normalized to 1
2931
a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
3032
b = X[:q + 1, :p + 1] @ a
3133

32-
return (a, b)
34+
return a.ravel(), b.ravel()
3335

3436

35-
def prony(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, ArrayLike]:
37+
def prony(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, float]:
3638
"""Least square minimization of poles to get denominator coefficients.
3739
3840
Solves directly (Pade method) to get numerator coefficients.
3941
Also calculates minimum error achieved.
4042
4143
Condition to energy_match is on page 575
4244
"""
43-
x = np.array(x)
44-
if p + q > len(x):
45+
_x = np.array(x).reshape(-1, 1)
46+
N = len(_x)
47+
if p + q > N:
4548
raise ValueError(f"Model order {p + q} is too large.")
4649

47-
# copy and make given signal column array
4850
X = convm(x, p + 1)
49-
print(X.shape)
50-
# M = p + q
51-
N = len(x)
52-
print(f"{N=}")
53-
xc = x.copy().reshape(-1, 1)
54-
55-
# Xq = X[q + 1:q + p + 1, 1:p + 1].copy()
56-
# a = np.linalg.solve(-Xq, X[q + 1: q + p + 1, 0])
57-
# a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
58-
# b = X[:q + 1, :p + 1] @ a
5951

6052
# the factorization does not guarantee nonsingularity!
6153
# 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
6658
rx = Xq_H @ Xq1
6759
Xinv = np.linalg.inv(Xq_H @ Xq)
6860
a = -Xinv @ rx
69-
print(a.shape)
7061
# a(0) normalized to 1
7162
a = np.concatenate((np.ones(1), a)).reshape(-1, 1)
7263
# same as Pade method
7364
b = X[:q + 1, :p + 1] @ a
7465

7566
# the minimum squared error
76-
err = np.transpose(xc[q + 1:N]) @ X[q + 1:N, :p + 1] @ a
67+
err = np.transpose(_x[q + 1:N]) @ X[q + 1:N, :p + 1] @ a
7768

78-
return a, b, err
69+
return a.ravel(), b.ravel(), err.ravel()[0]
7970

8071

8172
def shanks(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike, float]:
8273
"""Shank's method."""
83-
x = np.array(x)
84-
N = len(x)
74+
_x = np.array(x).ravel().reshape(-1, 1)
75+
N = len(_x)
8576
if p + q >= N:
8677
raise ValueError(f"Model order {p + q} is too large.")
8778

8879
a, _, _ = prony(x, p, q)
89-
a = np.array(a)
90-
print(f"{a.transpose().ravel()=}")
80+
logger.info(f"{a=}")
9181
u = np.concatenate((np.ones(1), np.zeros(N - 1)))
92-
res = sp.signal.tf2zpk([1], a.ravel())
93-
res = sp.signal.zpk2sos(*res)
94-
res = sp.signal.sosfilt(res, x=u)
95-
G = convm(res.ravel(), q + 1)
96-
G0 = G[:N,].copy()
97-
print(f"{G0.shape=}")
98-
G0_H = np.transpose((G0.copy()).conjugate())
99-
x0 = (x.copy()).reshape(-1, 1)
100-
gx = G0_H @ x0
101-
# the factorization does not guarantee nonsingularity!
102-
# resulting matrix is positive *semi*-definite
103-
Ginv = np.linalg.inv(G0_H @ G0)
104-
print(f"{x.shape=}")
105-
print(f"{Ginv.shape=}")
106-
b = Ginv @ gx
107-
err = 1
82+
sos = signal.tf2sos([1], a)
83+
g = signal.sosfilt(sos, x=u)
84+
logger.info(f"{g=}")
85+
G = convm(g, q + 1)
86+
G0 = G[:N].copy()
87+
logger.info(f"{G0=}")
88+
b = np.linalg.lstsq(G0, _x, rcond=None)[0]
89+
err = _x.T @ _x - _x.T @ G[:N, :q + 1] @ b
10890

109-
return a, b, err
91+
return a, b.ravel(), err.ravel()[0]
11092

11193

112-
def spike(g: ArrayLike, n0: int, n: int) -> ArrayLike:
94+
def spike(g: ArrayLike, n0: int, n: int) -> tuple[ArrayLike, float]:
11395
"""Leaset Squares Inverse Filter."""
114-
g = np.array(g).reshape(-1, 1)
115-
m = len(g)
96+
_g = np.array(g).reshape(-1, 1)
97+
m = len(_g)
11698

11799
if m + n - 1 <= n0:
118100
raise ValueError(f"m + n - 1 must be less than {n0=}")
119101

120102
G = convm(g, n)
121103
d = np.zeros((m + n - 1, 1))
122104
d[n0] = 1
123-
print(f"{d.shape=}, {G.shape=}")
124-
G_H = np.transpose(G.conjugate())
105+
h = np.linalg.lstsq(G, d, rcond=None)[0]
106+
err = 1 - G[n0,] @ h
125107

126-
print(f"{G_H.shape=}, {G.shape=}")
127-
Ginv = np.linalg.inv(G_H @ G)
128-
h = Ginv @ G_H @ d
129-
130-
return h
108+
return h.ravel(), err.ravel()[0]
131109

132110

133111
def ipf(x: ArrayLike, p: int, q: int, n: None, a: ArrayLike) -> NoReturn:
134112
"""Iterative Pre-Filtering."""
135113
raise NotImplementedError()
136114

137115

138-
def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, np.ndarray]:
116+
def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, float]:
139117
"""The auto-correlation method."""
140118
x0 = np.array(x).ravel().reshape(-1, 1)
141119
N = len(x0)
@@ -151,36 +129,49 @@ def acm(x: ArrayLike, p: int) -> tuple[np.ndarray, np.ndarray]:
151129
a = np.concatenate((np.ones(1), a1)).reshape(-1, 1)
152130
err = np.abs(X[:N + p, 0].T @ X @ a)
153131

154-
return a, err
132+
return a, err.ravel()[0]
155133

156134

157135
def covm(x: ArrayLike, p: int)-> tuple[ArrayLike, float]:
158136
"""Solve the complete Prony normal equations."""
159-
x0 = np.array(x).ravel().reshape(-1, 1)
160-
N = len(x0)
161-
if p >= len(x0):
137+
_x = np.array(x).ravel().reshape(-1, 1)
138+
N = len(_x)
139+
if p >= N:
162140
raise ValueError(f"{p=} all-pole model too large")
163141

164-
X = convm(x0, p + 1)
142+
X = convm(_x, p + 1)
165143
Xq = X[p - 1:N - 1, :p].copy()
166-
cx = X[p:N, 0].copy()
167-
Xq_H = Xq.copy().conjugate().transpose()
168-
print(f"{Xq=}")
169-
Xinv = np.linalg.inv(Xq_H @ Xq)
170-
a1 = -Xinv @ Xq_H @ cx
171-
a = np.concatenate((np.ones(1), a1)).reshape(-1, 1)
172-
err = np.abs(cx.transpose() @ X[p:N,] @ a)
173-
return a, err
144+
Xsol = np.linalg.lstsq(-Xq, X[p:N, 0], rcond=None)[0]
145+
logger.warning(f"{Xsol=}")
146+
a = np.hstack(([1], Xsol))
147+
err = np.abs(X[p:N,0] @ X[p:N,] @ a)
148+
return a, err.ravel()[0]
174149

175150

176151
def durbin(x: ArrayLike, p: int, q: int) -> tuple[ArrayLike, ArrayLike]:
177152
"""The durbin method."""
178-
x0 = np.array(x).ravel().reshape(-1, 1)
179-
# N = len(x0)
180-
if p >= len(x0):
153+
_x = np.array(x).ravel().reshape(-1, 1)
154+
N = len(_x)
155+
if p >= N:
181156
raise ValueError("p (all-pole model) too large")
182157

183-
a, eps = acm(x, p)
158+
a, eps = acm(_x, p)
184159
b, eps = acm(a / np.sqrt(eps), q)
185-
b /= np.sqrt(eps)
160+
b = b / np.sqrt(eps)
186161
return a, b
162+
163+
164+
def mywe(x: ArrayLike, p: int, q: int) -> NoReturn:
165+
"""Modified Yuler-Walker Systems of Equations.
166+
167+
Page 190.
168+
"""
169+
raise NotImplementedError()
170+
171+
172+
def eywe(x: ArrayLike, p: int, q: int) -> NoReturn:
173+
"""Extended Yuler-Walker Systems of Equations.
174+
175+
Page 193.
176+
"""
177+
raise NotImplementedError()

‎pyssp/optimal.py

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,13 @@
1-
"""Optimal Wiener Filters."""
1+
"""Optimal Filters, Chapter 7."""
22

3+
from typing import NoReturn
34

45
import numpy as np
5-
import numpy.typing as npt
6+
from numpy.typing import ArrayLike
67

78

8-
def kalman(y: list[float], A: npt.ArrayLike, C: npt.ArrayLike, sigmaw: list[float],
9-
sigmav: list[float]) -> tuple[npt.NDArray, ...]:
9+
def kalman(y: list[float], A: ArrayLike, C: ArrayLike, sigmaw: list[float],
10+
sigmav: list[float]) -> tuple[np.ndarray, ...]:
1011
"""Kalman Filter.
1112
1213
y: vector of observations N x q, n time steps, q sensors
@@ -73,3 +74,8 @@ def wiener_denoise() -> None:
7374
def wiener_systemid() -> None:
7475
"""Systemid based on FIR wiener filters."""
7576
raise NotImplementedError()
77+
78+
79+
def wiener_hopf(x: ArrayLike, p: int, q: int) -> NoReturn:
80+
"""Wiener-Hopf Systems of Equations."""
81+
raise NotImplementedError()

‎pyssp/spectrum.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
"""Power spectrum and Frequency Estimators."""
1+
"""Power spectrum and Frequency Estimators. Chapter 8."""
22

33
from typing import NoReturn
44

0 commit comments

Comments
 (0)