Refactor code.
alberto
1 month ago
parent
7948c92c74
commit
1cba3734ec
@ -0,0 +1,177 @@
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.ticker as ticker
|
||||
import numpy as np
|
||||
import numpy.linalg as LA
|
||||
import scienceplots # noqa: F401
|
||||
import scipy
|
||||
|
||||
# Requires latex installed
|
||||
from pygsp import graphs
|
||||
|
||||
from afgl.util.build_T_matrix import build_T_matrix
|
||||
from afgl.util.lanczos import lanczos
|
||||
from afgl.util.plot import latex_log_formatter
|
||||
|
||||
|
||||
def plot_graphs(G_ER, G_Sensor, s: np.ndarray, N: int, p: float) -> None:
|
||||
fig, axs = plt.subplots(2, 2, figsize=(6.6, 5))
|
||||
|
||||
# Set coordinates
|
||||
G_ER.set_coordinates()
|
||||
G_Sensor.set_coordinates()
|
||||
|
||||
signal_ER = filter_signal_with_fourier(G_ER, s)
|
||||
signal_S = filter_signal_with_fourier(G_Sensor, s)
|
||||
|
||||
# TOP LEFT
|
||||
G_ER.plot(s, ax=axs[0, 0], vertex_size=15, edge_width=0.5, edge_color="gray")
|
||||
axs[0, 0].set_title(rf"Erdős-Rényi Graph $(N = {N}, p = {p})$", pad=20)
|
||||
axs[0, 0].set_axis_off()
|
||||
|
||||
# BOTTOM LEFT
|
||||
G_ER.plot(
|
||||
signal_ER, ax=axs[1, 0], vertex_size=15, edge_width=0.5, edge_color="gray"
|
||||
)
|
||||
axs[1, 0].set_title("", pad=20)
|
||||
axs[1, 0].set_axis_off()
|
||||
|
||||
# TOP RIGHT
|
||||
G_Sensor.plot(s, ax=axs[0, 1], vertex_size=15, edge_width=0.5, edge_color="gray")
|
||||
axs[0, 1].set_title(rf"Sensor Network $(N = {N})$", pad=20)
|
||||
axs[0, 1].set_axis_off()
|
||||
|
||||
# BOTTOM RIGHT
|
||||
G_Sensor.plot(
|
||||
signal_S, ax=axs[1, 1], vertex_size=15, edge_width=0.5, edge_color="gray"
|
||||
)
|
||||
axs[1, 1].set_title("", pad=20)
|
||||
axs[1, 1].set_axis_off()
|
||||
|
||||
# Prevent label/title overlap
|
||||
plt.savefig("./out/printed_graphs.pdf", bbox_inches="tight")
|
||||
|
||||
|
||||
def g_extended(t: np.ndarray) -> np.ndarray:
|
||||
return np.sin(0.5 * np.pi * np.cos(np.pi * t) ** 2)
|
||||
|
||||
|
||||
"""
|
||||
Evaluates the function sin(0.5*pi*cos(pi*t)^2)chi_[-1/2,1/2] where chi_I is the
|
||||
characteristic function of I, as defined in example 1 (see [1]).
|
||||
"""
|
||||
|
||||
|
||||
def g(T: np.ndarray) -> np.ndarray:
|
||||
if scipy.sparse.issparse(T):
|
||||
# Operator & not supporting sparse matrix
|
||||
T = T.toarray()
|
||||
Chi = ((T >= -1 / 2) & (T <= 1 / 2)).astype(int)
|
||||
# Apply g_extended where Chi is True, else output 0
|
||||
return np.where(Chi, g_extended(T), 0)
|
||||
|
||||
|
||||
"""
|
||||
Computes the approximation g_M (see [1]) using Lanczos
|
||||
"""
|
||||
|
||||
|
||||
def compute_g_M(
|
||||
V: np.ndarray, alp: np.ndarray, beta: np.ndarray, s: np.ndarray
|
||||
) -> np.ndarray:
|
||||
M = len(alp)
|
||||
e_1 = np.zeros(M)
|
||||
e_1[0] = 1
|
||||
T = build_T_matrix(alp, beta)
|
||||
y = LA.norm(s) * (g(T) @ e_1)
|
||||
return V @ y
|
||||
|
||||
|
||||
def filter_signal_with_fourier(G, s: np.ndarray) -> np.ndarray:
|
||||
G.compute_fourier_basis()
|
||||
U = G.U
|
||||
return (U @ np.diag(g(G.e)) @ U.T) @ s
|
||||
|
||||
|
||||
def plot_error_comparison(
|
||||
l_err_ER: np.ndarray, t_err_ER: np.ndarray, l_err_S: np.ndarray, t_err_S: np.ndarray
|
||||
) -> None:
|
||||
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6.6, 2.5))
|
||||
|
||||
# Left plot (Erdos-Renyi)
|
||||
ax1.plot(l_err_ER, label=r"$\left\lVert g_{M+3} - g_M \right\rVert_2$")
|
||||
ax1.plot(t_err_ER, label=r"$\left\lVert e_M \right\rVert_2$")
|
||||
ax1.set_title("Erdős-Rényi graph")
|
||||
|
||||
# Right plot (Sensor)
|
||||
ax2.plot(l_err_S, label=r"$\left\lVert g_{M+3} - g_M \right\rVert_2$")
|
||||
ax2.plot(t_err_S, label=r"$\left\lVert e_M \right\rVert_2$")
|
||||
ax2.set_title("Sensor graph")
|
||||
|
||||
# Apply identical formatting to both subplots
|
||||
for ax in (ax1, ax2):
|
||||
ax.xaxis.set_major_locator(ticker.MultipleLocator(50))
|
||||
ax.set_yscale("log")
|
||||
ax.yaxis.set_major_formatter(ticker.FuncFormatter(latex_log_formatter))
|
||||
ax.legend()
|
||||
|
||||
# Prevents overlapping of labels between the subplots
|
||||
plt.tight_layout()
|
||||
|
||||
plt.savefig("./out/ex1_estimate.pdf", bbox_inches="tight")
|
||||
|
||||
|
||||
def run_comparison_1_for_graph(
|
||||
G, s: np.ndarray, M_MAX: int
|
||||
) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""Compares the error with error generated by Lanczos.
|
||||
|
||||
Args:
|
||||
G: Graph
|
||||
s: Signal vector
|
||||
|
||||
Returns:
|
||||
[lanczos_err, true_err]: Vectors of errors norm(g_{M+3} - g_M) and norm(e_M)
|
||||
as defined in [1]
|
||||
"""
|
||||
G.compute_laplacian("combinatorial")
|
||||
L = G.L
|
||||
|
||||
j = 3
|
||||
V, alp, beta = lanczos(L, s, M_MAX + j)
|
||||
|
||||
lanczos_err = np.zeros(M_MAX + j)
|
||||
true_err = np.zeros(M_MAX + j)
|
||||
|
||||
GLs = filter_signal_with_fourier(G, s)
|
||||
|
||||
for M in range(2, M_MAX + j):
|
||||
g_M = compute_g_M(V[:, 0:M], alp[0:M], beta[0 : M - 1], s)
|
||||
g_Mj = compute_g_M(V[:, 0 : M + j], alp[0 : M + j], beta[0 : M + j - 1], s)
|
||||
|
||||
lanczos_err[M - 1] = LA.norm(g_Mj - g_M)
|
||||
true_err[M - 1] = LA.norm(GLs - g_M)
|
||||
|
||||
return lanczos_err, true_err
|
||||
|
||||
|
||||
def run() -> None:
|
||||
"""Ripete il test corrispondente ad Example 1 dell'articolo limitandosi al
|
||||
metodo di Lanczos (no Chebyshev) e utilizzando come funzione g(t) = sin(0.5π
|
||||
cos(πt)2) * \chi_{[-0.5, 0.5]}.
|
||||
"""
|
||||
N = 500
|
||||
M_MAX = 200
|
||||
p = 0.04
|
||||
|
||||
s = np.random.randint(1, 10000, N).astype(float)
|
||||
# Normalize s as in request
|
||||
s /= LA.norm(s)
|
||||
|
||||
G_ER = graphs.ErdosRenyi(N, p)
|
||||
G_S = graphs.Sensor(N)
|
||||
|
||||
l_err_ER, t_err_ER = run_comparison_1_for_graph(G_ER, s, M_MAX)
|
||||
l_err_S, t_err_S = run_comparison_1_for_graph(G_S, s, M_MAX)
|
||||
|
||||
plot_error_comparison(l_err_ER, t_err_ER, l_err_S, t_err_S)
|
||||
plot_graphs(G_ER, G_S, s, N, p)
|
||||
@ -0,0 +1,13 @@
|
||||
def run() -> None:
|
||||
"""Genera i grafi di di Erdos-Reny di grandezza crescente (ad esempio 250,
|
||||
500,1000, 2000, 4000) e parametro p = 0.04 e misura il tempo
|
||||
computazionale del metodo di Lanczos utilizzando come soglia per il criterio
|
||||
d'arresto epsilon = 10^-2 (o una soglia a scelta).
|
||||
|
||||
Args:
|
||||
None
|
||||
|
||||
Returns:
|
||||
None
|
||||
"""
|
||||
pass
|
||||
@ -0,0 +1,28 @@
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import scienceplots # noqa: F401
|
||||
from pygsp import plotting
|
||||
|
||||
|
||||
def latex_sci(val: float, decimals: int = 2) -> str:
|
||||
"""Converts a value to LaTeX scientific notation A x 10^{B}."""
|
||||
if val == 0:
|
||||
return "0"
|
||||
exponent = int(np.floor(np.log10(abs(val))))
|
||||
mantissa = val / 10**exponent
|
||||
return rf"{mantissa:.{decimals}f} \times 10^{{{exponent}}}"
|
||||
|
||||
|
||||
def latex_log_formatter(y: float, pos: int) -> str:
|
||||
"""Custom formatter to render tick labels as LaTeX 10^{n}."""
|
||||
if y <= 0:
|
||||
return ""
|
||||
# Extract the exponent using log10
|
||||
n = int(np.round(np.log10(y)))
|
||||
return f"$10^{{{n}}}$"
|
||||
|
||||
|
||||
def plot_setup() -> None:
|
||||
plotting.BACKEND = "matplotlib"
|
||||
plt.style.use(["science"])
|
||||
# TODO match font with document
|
||||
Loading…
Reference in New Issue