diff --git a/Old/1_fft_basics/1.py b/Old/1_fft_basics/1.py new file mode 100644 index 0000000..666f9f4 --- /dev/null +++ b/Old/1_fft_basics/1.py @@ -0,0 +1,49 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.fft import fft, ifft, fftfreq +p = True +N = 100 +x, dx = np.linspace(0, 2*np.pi, N, endpoint=False, retstep=True) +f_x = 1 + 2*np.sin(x) + 10*np.sin(5*x) + 3*np.cos(5*x) +F_k = fft(f_x) +if p: + print("dx:", dx) + print("x:", x) + print("f(x):", f_x) + print("F_k:", F_k) + +f_x_reconstructed = ifft(F_k) + +# --- VYKRESLENÍ VÝSLEDKŮ PRO KONTROLU --- +plt.figure(figsize=(12, 8)) + +# Graf 1: Původní a složená funkce (měly by se překrývat) +plt.subplot(2, 1, 1) +plt.plot(x, f_x, label="Původní f(x)", linewidth=3) +# ifft vrací komplexní čísla, pro graf vezmeme reálnou část (.real) +plt.plot(x, f_x_reconstructed.real, '--', label="Složená (iFFT)", color='red') +plt.title("Původní signál vs. Rekonstruovaný signál (Úspěšný test!)") +plt.legend() +plt.grid() + +# Graf 2: Zobrazení frekvenčního spektra (co vlastně FFT našla) +plt.subplot(2, 1, 2) +# Spočítáme frekvence pro osu X (zajímají nás jen kladné frekvence, proto úprava do N//2) +xf = fftfreq(N, dx)[:N//2] +print("xf:", xf) +# Amplituda (síla) signálu - musíme normalizovat dělením (2/N) +amplitudy = (2.0/N) * np.abs(F_k[:N//2]) +print("Amplitudy:", amplitudy) +# Ruční úprava nulté frekvence (DC složky - to je ta "1" na začátku rovnice) +amplitudy[0] = amplitudy[0] / 2 + +# Vykreslíme to jako sloupečky (stem plot) +plt.stem(xf, amplitudy) +plt.xlim(-0.1, 2) # Přiblížíme začátek grafu, aby byly vidět špičky +plt.title("Spektrum (FFT) - Všimni si špiček na frekvencích, které odpovídají rovnici") +plt.xlabel("Frekvence (Hz)") +plt.ylabel("Amplituda") +plt.grid() + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/Old/1_fft_basics/1ext.py b/Old/1_fft_basics/1ext.py new file mode 100644 index 0000000..fc32847 --- /dev/null +++ b/Old/1_fft_basics/1ext.py @@ -0,0 +1,140 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.fft import fft, ifft, fftfreq +p = False + +# ========================================== +# 1. Definice mřížky (Grid) a signálu +# ========================================== +N = 100 +x, dx = np.linspace(0, 2*np.pi, N, endpoint=False, retstep=True) + +# Rozložení signálu na jednotlivé komponenty (abychom viděli, co se z čeho skládá) +f1 = 1 * np.ones_like(x) # DC složka (konstanta 1) +f2 = 2 * np.sin(x) # Frekvence 1 +f3 = 10 * np.sin(5*x) # Frekvence 5 (sinusovka) +f4 = 3 * np.cos(5*x) # Frekvence 5 (kosinusovka) + +# Výsledná složená funkce ze zadání +f_x = f1 + f2 + f3 + f4 + +if p: + print("dx:", dx) + print("x:", x) + print("f(x):", f_x) + +# ========================================== +# 2. Výpočet FFT a zpětné transformace +# ========================================== +F_k = fft(f_x) +f_x_reconstructed = ifft(F_k) + +# Výpočet frekvencí pro osu X +xf_vsechny = fftfreq(N, dx) + +# Pro srozumitelné spektrum nás většinou zajímá jen kladná část spektra +xf_kladne = xf_vsechny[:N//2] +F_k_kladne = F_k[:N//2] + +# Amplitudy (síla signálu) - normalizace pomocí (2/N) +amplitudy = (2.0/N) * np.abs(F_k_kladne) +amplitudy[0] = amplitudy[0] / 2 # DC složka se nedělí dvěma, vracíme ji na správnou hodnotu + +# Fáze (posun signálu v radiánech) +faze = np.angle(F_k_kladne) +# Odfiltrování šumu z fáze: tam, kde je nulová nebo minimální amplituda, nemá smysl počítat fázi +faze[amplitudy < 0.1] = 0 + +# Výkonové spektrum (Power Spectrum - energie na dané frekvenci) +vykon = np.abs(F_k_kladne)**2 + + +# ========================================== +# 3. VYKRESLENÍ - OBRÁZEK 1: Časová oblast +# ========================================== +plt.figure(figsize=(12, 10)) +plt.subplot(3, 1, 1) +plt.plot(x, f_x, label="Výsledný signál f(x)", color='black', linewidth=2) +plt.title("1. Celkový složený signál f(x) v čase") +plt.ylabel("Amplituda") +plt.grid(True); plt.legend() + +plt.subplot(3, 1, 2) +plt.plot(x, f1, label="DC (1)", linestyle='--') +plt.plot(x, f2, label="2*sin(x)", linestyle='--') +plt.plot(x, f3, label="10*sin(5x)", linestyle='--') +plt.plot(x, f4, label="3*cos(5x)", linestyle='--') +plt.title("2. Jednotlivé " + "skryté" + " složky, ze kterých se signál skládá") +plt.ylabel("Amplituda") +plt.grid(True); plt.legend() + +plt.subplot(3, 1, 3) +plt.plot(x, f_x, label="Původní f(x)", linewidth=4, alpha=0.5) +plt.plot(x, f_x_reconstructed.real, '--', label="Složená (iFFT)", color='red') +plt.title("3. Původní signál vs. Rekonstruovaný signál (Úspěšný test!)") +plt.xlabel("x (Čas / Prostor)") +plt.ylabel("Amplituda") +plt.grid(True); plt.legend() +plt.tight_layout() + + +# ========================================== +# VYKRESLENÍ - OBRÁZEK 2: Surová komplexní čísla z FFT +# ========================================== +plt.figure(figsize=(12, 8)) +# Reálná část reprezentuje podíly KOSINUSOVÝCH funkcí +plt.subplot(2, 1, 1) +plt.stem(range(N), F_k.real, basefmt=" ") +plt.title("4. Surový výstup FFT: Reálná část (Odpovídá zastoupení KOSINUSŮ ve frekvencích)") +plt.ylabel("Reálná hodnota") +plt.grid(True) + +# Imaginární část reprezentuje podíly SINUSOVÝCH funkcí +plt.subplot(2, 1, 2) +plt.stem(range(N), F_k.imag, basefmt=" ") +plt.title("5. Surový výstup FFT: Imaginární část (Odpovídá zastoupení SINUSŮ ve frekvencích)") +plt.xlabel("Index pole (frekvenční koš)") +plt.ylabel("Imaginární hodnota") +plt.grid(True) +plt.tight_layout() + + +# ========================================== +# VYKRESLENÍ - OBRÁZEK 3: Polární forma (Amplituda a Fáze) +# ========================================== +plt.figure(figsize=(12, 8)) +plt.subplot(2, 1, 1) +plt.stem(xf_kladne, amplitudy, basefmt=" ") +plt.xlim(-0.1, 1.2) # Ořízneme graf jen na zajímavou část, ať vidíme špičky +plt.title("6. Amplitudové spektrum - To nejdůležitější (Síla jednotlivých frekvencí)") +plt.ylabel("Amplituda") +plt.grid(True) + +# Přidání textu přímo nad špičky v grafu pro přehlednost +for i, amp in enumerate(amplitudy): + if amp > 0.5: + plt.text(xf_kladne[i], amp + 0.5, f"{amp:.1f}", ha='center', fontweight='bold') + +plt.subplot(2, 1, 2) +plt.stem(xf_kladne, faze, basefmt=" ") +plt.xlim(-0.1, 1.2) +plt.title("7. Fázové spektrum - Fázový posun odhalených vln (v radiánech)") +plt.xlabel("Frekvence") +plt.ylabel("Fáze [rad]") +plt.grid(True) +plt.tight_layout() + + +# ========================================== +# VYKRESLENÍ - OBRÁZEK 4: Výkonové spektrum (Power Spectrum) +# ========================================== +plt.figure(figsize=(12, 4)) +plt.stem(xf_kladne, vykon, basefmt=" ", markerfmt='ro', linefmt='r-') +plt.xlim(-0.1, 1.2) +plt.title("8. Výkonové spektrum (Power Spectral Density) - Kde je ukryta většina energie signálu") +plt.xlabel("Frekvence") +plt.ylabel("Výkon (Absolutní hodnota na druhou)") +plt.grid(True) +plt.tight_layout() + +plt.show() \ No newline at end of file diff --git a/Old/2_grf_1d/fft1d.py b/Old/2_grf_1d/fft1d.py new file mode 100644 index 0000000..cf56f23 --- /dev/null +++ b/Old/2_grf_1d/fft1d.py @@ -0,0 +1,108 @@ +""" +Generování náhodných polí pomocí spektrálních metod (FFT / NUFFT). +Úkoly 3.3: + 1) make_hermitian – reálná funkce z náhodných komplexních koeficientů + 2) generate_grf_fft – GRF na pravidelné mřížce (Gaussovská korelace) + 3) generate_grf_nufft – GRF na nepravidelné mřížce (finufft) + + literatura + +""" + +import numpy as np +from scipy.fft import fftfreq +import finufft + + +# SPEKTRÁLNÍ HUSTOTA + +def spectral_density_gaussian(omega, phi, sigma=1.0): + # C(r) = sigma^2 * exp(-r^2 / 2phi^2) => S(w) = sigma^2 * sqrt(2pi)*phi * exp(-w^2*phi^2/2) + return sigma**2 * np.sqrt(2 * np.pi) * phi * np.exp(-0.5 * (omega * phi)**2) + +def spectral_density_exponential(omega, phi, sigma=1.0): + # C(r) = sigma^2 * exp(-|r|/phi) => S(w) = sigma^2 * 2phi / (1 + (w*phi)^2) + return sigma**2 * 2 * phi / (1 + (omega * phi)**2) + + +# HERMITOVSKÁ SYMETRIE – zaručí reálný výstup IFFT + +def make_hermitian(f_hat, centered=False): + """ + Vyrobí hermitovsky symetrické koeficienty: f_hat[-k] = conj(f_hat[k]) + + centered=False – FFT pořadí: DC na indexu 0, záporné frekvence na konci + centered=True – centrované pořadí: DC uprostřed na indexu N//2 + """ + N = len(f_hat) + f = f_hat.copy() + + if centered: + zi = N // 2 # index DC (nulové frekvence) + f[zi] = f[zi].real + for i in range(1, N // 2): + f[zi - i] = np.conj(f[zi + i]) # záporné frekvence = sdružené kladných + else: + f[0] = f[0].real # DC + if N % 2 == 0: + f[N // 2] = f[N // 2].real # Nyquistova frekvence + for k in range(1, N // 2): + f[-k] = np.conj(f[k]) + + return f + + +# FFT – PRAVIDELNÁ MŘÍŽKA + +def generate_grf_fft(L, N, phi, sigma=1.0, seed=None, + spectral_density_fn=spectral_density_gaussian): + """ + Generuje 1D náhodné pole na pravidelné mřížce. + - náhodné koeficienty z N(0,1) komplexního rozdělení + - modulace filtrem sqrt(S(omega)) + - hermitovská symetrie -> reálný výstup přes IFFT + spectral_density_fn: spectral_density_gaussian | spectral_density_exponential + """ + rng = np.random.default_rng(seed) + dx = L / N + x = np.arange(N) * dx + + omega = 2 * np.pi * fftfreq(N, d=dx) # FFT pořadí + S_k = spectral_density_fn(omega, phi, sigma) + + white = (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2) + f_hat = make_hermitian(white * np.sqrt(S_k), centered=False) + + field = np.fft.ifft(f_hat).real * N / np.sqrt(L) + return x, field + + +# NUFFT – NEPRAVIDELNÁ MŘÍŽKA + +def generate_grf_nufft(L, N_freq, M_points, phi, sigma=1.0, seed=None, x_points=None, + spectral_density_fn=spectral_density_gaussian): + """ + Generuje 1D náhodné pole na nepravidelné mřížce pomocí NUFFT typu 2 + (pravidelné frekvence -> nepravidelné prostorové body). + finufft očekává souřadnice v [-pi, pi]. + """ + rng = np.random.default_rng(seed) + + # centrované pořadí: k = -N/2, ..., -1, 0, +1, ..., N/2-1 + k = np.arange(-N_freq // 2, N_freq // 2) + omega = k * (2 * np.pi / L) + S_omega = spectral_density_fn(omega, phi, sigma) + + white = (rng.standard_normal(N_freq) + 1j * rng.standard_normal(N_freq)) / np.sqrt(2) + f_hat = make_hermitian(white * np.sqrt(S_omega), centered=True) + + if x_points is None: + x_nufft = rng.uniform(-np.pi, np.pi, M_points) + else: + x_nufft = (np.asarray(x_points) / L) * 2 * np.pi - np.pi # [0,L] -> [-pi,pi] + + field = finufft.nufft1d2(x_nufft, f_hat.astype(np.complex128)) + x_final = (x_nufft + np.pi) * (L / (2 * np.pi)) # [-pi,pi] -> [0,L] + return x_final, field.real + + diff --git a/Old/2_grf_1d/main2.py b/Old/2_grf_1d/main2.py new file mode 100644 index 0000000..5e04df5 --- /dev/null +++ b/Old/2_grf_1d/main2.py @@ -0,0 +1,35 @@ +import matplotlib.pyplot as plt +import numpy as np +import fft1d as rf + +L, N, phi = 100.0, 1024, 3.0 + +fig, axes = plt.subplots(2, 2, figsize=(14, 8)) + +# --- Gaussovská --- +x, f = rf.generate_grf_fft(L, N, phi, seed=42) +axes[0][0].plot(x, f, lw=1) +axes[0][0].set_title(f"FFT Gaussovská (φ={phi})") + +x_nu, f_nu = rf.generate_grf_nufft(L, 256, 500, phi, seed=42) +idx = np.argsort(x_nu) +axes[0][1].plot(x_nu[idx], f_nu[idx], '-', alpha=0.4, color='gray', lw=0.8) +axes[0][1].scatter(x_nu, f_nu, s=12, c=f_nu, cmap='viridis', zorder=3) +axes[0][1].set_title(f"NUFFT Gaussovská (φ={phi})") + +# --- Exponenciální --- +x, f = rf.generate_grf_fft(L, N, phi, seed=42, spectral_density_fn=rf.spectral_density_exponential) +axes[1][0].plot(x, f, lw=1) +axes[1][0].set_title(f"FFT Exponenciální (φ={phi})") + +x_nu, f_nu = rf.generate_grf_nufft(L, 256, 500, phi, seed=42, spectral_density_fn=rf.spectral_density_exponential) +idx = np.argsort(x_nu) +axes[1][1].plot(x_nu[idx], f_nu[idx], '-', alpha=0.4, color='gray', lw=0.8) +axes[1][1].scatter(x_nu, f_nu, s=12, c=f_nu, cmap='viridis', zorder=3) +axes[1][1].set_title(f"NUFFT Exponenciální (φ={phi})") + +for ax in axes.flat: + ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/Old/3_grf_2d/3.2d.py b/Old/3_grf_2d/3.2d.py new file mode 100644 index 0000000..a7a38c3 --- /dev/null +++ b/Old/3_grf_2d/3.2d.py @@ -0,0 +1,216 @@ +""" +Upravit zadání + +""" + +import random +import numpy as np +import finufft +# ============================================================================= +# CORRELATION TŘÍDY +# ============================================================================= + +class GaussianCorrelation: + """ + C(r) = sigma^2 * exp(-|r|^2 / 2*phi^2) + S_nD(omega) = sigma^2 * (sqrt(2pi)*phi)^dim * exp(-|omega|^2*phi^2/2) + """ + def __init__(self, L, N_freq, phi, sigma=1.0, dim=1): + self.phi = phi + self.sigma = sigma + self.dim = dim + self.N_freq = N_freq + self.L = L + k = np.arange(-N_freq // 2, N_freq // 2) + self.f_points = k * (2 * np.pi / L) # 1D osa frekvencí [rad/m] + + def spectral_density(self, omega_mag): + # omega_mag = |omega| (skalár nebo pole) + return (self.sigma**2 + * (np.sqrt(2 * np.pi) * self.phi)**self.dim + * np.exp(-0.5 * (omega_mag * self.phi)**2)) + + +class ExponentialCorrelation: + """ + C(r) = sigma^2 * exp(-|r|/phi) + 1D: S(w) = sigma^2 * 2*phi / (1 + (w*phi)^2) + 2D: S(w) = sigma^2 * 2*pi*phi^2 / (1 + (|w|*phi)^2)^(3/2) + """ + def __init__(self, L, N_freq, phi, sigma=1.0, dim=1): + self.phi = phi + self.sigma = sigma + self.dim = dim + self.N_freq = N_freq + self.L = L + k = np.arange(-N_freq // 2, N_freq // 2) + self.f_points = k * (2 * np.pi / L) + + def spectral_density(self, omega_mag): + if self.dim == 1: + return self.sigma**2 * 2 * self.phi / (1 + (omega_mag * self.phi)**2) + elif self.dim == 2: + return self.sigma**2 * 2 * np.pi * self.phi**2 / (1 + (omega_mag * self.phi)**2)**1.5 + else: + raise NotImplementedError("Exponential: dim > 2 není implementováno") + + +# ============================================================================= +# HERMITOVSKÁ SYMETRIE – zaručí reálný výstup IFFT +# centrované pořadí: DC uprostřed na indexu N//2 (v každé dimenzi) +# podmínka: f_hat[-k] = conj(f_hat[k]) pro všechny k +# ============================================================================= + +def make_hermitian_nd(f_hat): + """ + nD hermitovská symetrie v centrovaném pořadí. + f_hat shape: (N,) pro 1D, (N, N) pro 2D, (N, N, N) pro 3D. + """ + N = f_hat.shape[0] + zi = N // 2 # index DC v každé dimenzi + f = f_hat.copy() + + # DC (všechny indexy == zi) musí být reálná + dc_idx = tuple([zi] * f.ndim) + f[dc_idx] = f[dc_idx].real + + # pro každý multi-index k: f[-k] = conj(f[k]) + # iterujeme přes všechny indexy, záporné nastavíme jako sdružené kladných + for idx in np.ndindex(*f.shape): + # přeskočit DC a body s nulovým indexem (ty jsou sdružené sami sobě) + shifted = tuple(i - zi for i in idx) + if all(s <= 0 for s in shifted): + continue + neg_idx = tuple((zi - s) % N for s in shifted) + f[neg_idx] = np.conj(f[idx]) + + return f + + +def make_white_noise(N_freq, dim=1, seed=None): + """Komplexní bílý šum tvaru (N_freq,)*dim.""" + rng = np.random.default_rng(seed) + shape = (N_freq,) * dim + size = N_freq ** dim + flat = (rng.standard_normal(size) + 1j * rng.standard_normal(size)) / np.sqrt(2) + return flat.reshape(shape) + + +# ============================================================================= +# GENERATE GRF +# ============================================================================= + +def generate_grf_nufft(x_points, corr, weights=None, seed=None): + """ + Generuje náhodné pole v dim dimenzích pomocí NUFFT typu 2. + + x_points : (M,) pro 1D + (M, 2) pro 2D + (M, 3) pro 3D + corr : GaussianCorrelation | ExponentialCorrelation (s dim nastaveným) + weights : bílý šum tvaru (N_freq,)*dim; pokud None, vygeneruje se nový + """ + dim = corr.dim + N_freq = corr.N_freq + + if weights is None: + weights = make_white_noise(N_freq, dim=dim, seed=seed) + + # frekvenční mřížka v nD: |omega| pro každý bod + if dim == 1: + omega_mag = np.abs(corr.f_points) + else: + grids = np.meshgrid(*([corr.f_points] * dim), indexing='ij') + omega_mag = np.sqrt(sum(g**2 for g in grids)) # shape (N_freq,)*dim + + S = corr.spectral_density(omega_mag) + f_hat = make_hermitian_nd(weights * np.sqrt(S)) + + # souřadnice do [-pi, pi] + x_points = np.asarray(x_points) + L = corr.L + + if dim == 1: + x_nu = (x_points / L) * 2 * np.pi - np.pi + field = finufft.nufft1d2(x_nu, f_hat.astype(np.complex128)) + + elif dim == 2: + x_nu = (x_points[:, 0] / L) * 2 * np.pi - np.pi + y_nu = (x_points[:, 1] / L) * 2 * np.pi - np.pi + field = finufft.nufft2d2(x_nu, y_nu, f_hat.astype(np.complex128)) + + elif dim == 3: + x_nu = (x_points[:, 0] / L) * 2 * np.pi - np.pi + y_nu = (x_points[:, 1] / L) * 2 * np.pi - np.pi + z_nu = (x_points[:, 2] / L) * 2 * np.pi - np.pi + field = finufft.nufft3d2(x_nu, y_nu, z_nu, f_hat.astype(np.complex128)) + + else: + raise NotImplementedError("dim > 3 není podporováno finufft") + + return field.real + + +def generate_grf_fft(x_points, corr, weights=None, seed=None): + """Alias pro pravidelnou mřížku – volá generate_grf_nufft.""" + return generate_grf_nufft(x_points, corr, weights=weights, seed=seed) + + +# ============================================================================= +# demo +# ============================================================================= + +if __name__ == "__main__": + import matplotlib.pyplot as plt + + L, N, phi = 100.0, 256, 3.0 + seed = random.randint(3, 9) + # --- 1D --- + x = np.linspace(0, L, N) + x_irr = np.sort(np.random.uniform(0, L, 300)) + white1d = make_white_noise(N, dim=1, seed=seed) + + fig, axes = plt.subplots(2, 2, figsize=(14, 8)) + for row, (label, CorrClass) in enumerate([("Gaussovská", GaussianCorrelation), + ("Exponenciální", ExponentialCorrelation)]): + corr = CorrClass(L, N, phi, dim=1) + axes[row][0].plot(x, generate_grf_fft(x, corr, weights=white1d), lw=1) + axes[row][0].set_title(f"1D FFT {label} (φ={phi})") + f_nu = generate_grf_nufft(x_irr, corr, weights=white1d) + axes[row][1].plot(x_irr, f_nu, '-', alpha=0.4, color='gray', lw=0.8) + axes[row][1].scatter(x_irr, f_nu, s=12, c=f_nu, cmap='viridis', zorder=3) + axes[row][1].set_title(f"1D NUFFT {label} (φ={phi})") + for ax in axes.flat: ax.grid(True, alpha=0.3) + plt.tight_layout(); plt.show() + + # --- 2D --- + N2 = 64 # menší kvůli N^2 bodům + white2d = make_white_noise(N2, dim=2, seed=seed) + g = np.linspace(0, L, N2) + xx, yy = np.meshgrid(g, g) + pts2d = np.column_stack([xx.ravel(), yy.ravel()]) + + # nepravidelné body pro NUFFT + pts2d_irr = np.column_stack([np.random.uniform(0, L, N2**2), + np.random.uniform(0, L, N2**2)]) + + fig, axes = plt.subplots(2, 2, figsize=(12, 10)) + for col, (label, CorrClass) in enumerate([("Gaussovská", GaussianCorrelation), + ("Exponenciální", ExponentialCorrelation)]): + corr2d = CorrClass(L, N2, phi, dim=2) + + # pravidelná mřížka (FFT) + field2d = generate_grf_nufft(pts2d, corr2d, weights=white2d).reshape(N2, N2) + im = axes[0][col].imshow(field2d, extent=[0, L, 0, L], origin='lower', cmap='viridis') + plt.colorbar(im, ax=axes[0][col]) + axes[0][col].set_title(f"2D FFT {label} (φ={phi})") + + # nepravidelná mřížka (NUFFT) – interpolováno na pravidelnou mřížku pro vizualizaci + f_irr = generate_grf_nufft(pts2d_irr, corr2d, weights=white2d) + from scipy.interpolate import griddata + field_interp = griddata(pts2d_irr, f_irr, (xx, yy), method='linear') + im2 = axes[1][col].imshow(field_interp, extent=[0, L, 0, L], origin='lower', cmap='viridis') + plt.colorbar(im2, ax=axes[1][col]) + axes[1][col].set_title(f"2D NUFFT {label} (φ={phi}, interpolováno)") + + plt.tight_layout(); plt.show() \ No newline at end of file diff --git a/Old/3_grf_2d/3.py b/Old/3_grf_2d/3.py new file mode 100644 index 0000000..7a4d587 --- /dev/null +++ b/Old/3_grf_2d/3.py @@ -0,0 +1,168 @@ +""" +Generování náhodných polí pomocí spektrálních metod (FFT / NUFFT). +""" +import numpy as np +from scipy.fft import fftfreq +import finufft +# ============================================================================= +# CORRELATION TŘÍDY +# ============================================================================= +class GaussianCorrelation: + """ + C(r) = sigma^2 * exp(-r^2 / 2*phi^2) + S(w) = sigma^2 * sqrt(2pi) * phi * exp(-w^2*phi^2 / 2) + """ + def __init__(self, L, N_freq, phi, sigma=1.0): + self.phi = phi + self.sigma = sigma + # frekvenční body v centrovaném pořadí: k = -N/2, ..., N/2-1 + k = np.arange(-N_freq // 2, N_freq // 2) + self.f_points = k * (2 * np.pi / L) # úhlové frekvence [rad/m] + + def spectral_density(self, omega): + # S(w) = sigma^2 * sqrt(2pi) * phi * exp(-w^2*phi^2/2) + return self.sigma**2 * np.sqrt(2 * np.pi) * self.phi * np.exp(-0.5 * (omega * self.phi)**2) + +class ExponentialCorrelation: + """ + C(r) = sigma^2 * exp(-|r|/phi) + S(w) = sigma^2 * 2*phi / (1 + (w*phi)^2) + """ + def __init__(self, L, N_freq, phi, sigma=1.0): + self.phi = phi + self.sigma = sigma + k = np.arange(-N_freq // 2, N_freq // 2) + self.f_points = k * (2 * np.pi / L) + + def spectral_density(self, omega): + # S(w) = sigma^2 * 2*phi / (1 + (w*phi)^2) + return self.sigma**2 * 2 * self.phi / (1 + (omega * self.phi)**2) + +# ============================================================================= +# HERMITOVSKÁ SYMETRIE – zaručí reálný výstup IFFT +# ============================================================================= +def make_hermitian(f_hat, centered=False): + """ + Vyrobí hermitovsky symetrické koeficienty: f_hat[-k] = conj(f_hat[k]) + + centered=False – FFT pořadí: DC na indexu 0 + centered=True – centrované pořadí: DC uprostřed na indexu N//2 + """ + N = len(f_hat) + f = f_hat.copy() + if centered: + zi = N // 2 # index DC + f[zi] = f[zi].real + for i in range(1, N // 2): + f[zi - i] = np.conj(f[zi + i]) # f[-k] = conj(f[k]) + else: + f[0] = f[0].real # DC + if N % 2 == 0: + f[N // 2] = f[N // 2].real # Nyquistova frekvence + for k in range(1, N // 2): + f[-k] = np.conj(f[k]) + return f + + +def make_white_noise(N, seed=None): + """Komplexní bílý šum ze standardního normálního rozdělení.""" + rng = np.random.default_rng(seed) + return (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2) + + +# ============================================================================= +# FFT – PRAVIDELNÁ MŘÍŽKA +# ============================================================================= + +def generate_grf_fft(x_points, corr, weights=None, seed=None): + """ + Generuje 1D náhodné pole na pravidelné mřížce pomocí FFT. + + x_points : ndarray (N,) – pravidelná prostorová mřížka + corr : GaussianCorrelation | ExponentialCorrelation + weights : komplexní koeficienty šumu (N_freq,); pokud None, vygenerují se nové + """ + N_freq = len(corr.f_points) + if weights is None: + weights = make_white_noise(N_freq, seed=seed) + + S = corr.spectral_density(corr.f_points) + f_hat = make_hermitian(weights * np.sqrt(S), centered=True) + + # stejná code path jako NUFFT -> identické výsledky na pravidelné mřížce + return generate_grf_nufft(x_points, corr, weights=weights) + + +# ============================================================================= +# NUFFT – NEPRAVIDELNÁ MŘÍŽKA +# ============================================================================= + +def generate_grf_nufft(x_points, corr, weights=None, seed=None): + """ + Generuje 1D náhodné pole na nepravidelné mřížce pomocí NUFFT typu 2. + + x_points : ndarray (M,) – libovolné prostorové souřadnice v [0, L] + corr : GaussianCorrelation | ExponentialCorrelation + weights : komplexní koeficienty šumu (N_freq,); pokud None, vygenerují se nové + -> stejné weights jako u FFT = srovnatelné realizace + """ + N_freq = len(corr.f_points) + if weights is None: + weights = make_white_noise(N_freq, seed=seed) + + S = corr.spectral_density(corr.f_points) + f_hat = make_hermitian(weights * np.sqrt(S), centered=True) + + # finufft očekává souřadnice v [-pi, pi] + L = x_points[-1] - x_points[0] + x_nufft = (x_points / L) * 2 * np.pi - np.pi # [0,L] -> [-pi,pi] + + field = finufft.nufft1d2(x_nufft, f_hat.astype(np.complex128)) + return field.real + + +# ============================================================================= +# demo.py +# ============================================================================= + +if __name__ == "__main__": + import matplotlib.pyplot as plt + + L, N, phi = 100.0, 1024, 3.0 + N_freq = N # stejný pro FFT i NUFFT -> srovnatelné realizace + + x_regular = np.linspace(0, L, N) + x_irregular = np.sort(np.random.uniform(0, L, 500)) + + corr_g = GaussianCorrelation(L, N_freq, phi) + corr_e = ExponentialCorrelation(L, N_freq, phi) + + white = make_white_noise(N_freq) # stejný šum pro obě metody + + fig, axes = plt.subplots(2, 2, figsize=(14, 8)) + + # Gaussovská + f_fft = generate_grf_fft(x_regular, corr_g, weights=white) + axes[0][0].plot(x_regular, f_fft, lw=1) + axes[0][0].set_title(f"FFT Gaussovská (φ={phi})") + + f_nu = generate_grf_nufft(x_irregular, corr_g, weights=white) + axes[0][1].plot(x_irregular, f_nu, "-", alpha=0.4, color="gray", lw=0.8) + axes[0][1].scatter(x_irregular, f_nu, s=12, c=f_nu, cmap="viridis", zorder=3) + axes[0][1].set_title(f"NUFFT Gaussovská (φ={phi})") + + # Exponenciální + f_fft = generate_grf_fft(x_regular, corr_e, weights=white) + axes[1][0].plot(x_regular, f_fft, lw=1) + axes[1][0].set_title(f"FFT Exponenciální (φ={phi})") + + f_nu = generate_grf_nufft(x_irregular, corr_e, weights=white) + axes[1][1].plot(x_irregular, f_nu, "-", alpha=0.4, color="gray", lw=0.8) + axes[1][1].scatter(x_irregular, f_nu, s=12, c=f_nu, cmap="viridis", zorder=3) + axes[1][1].set_title(f"NUFFT Exponenciální (φ={phi})") + + for ax in axes.flat: + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/Old/3_grf_2d/Figure_1.png b/Old/3_grf_2d/Figure_1.png new file mode 100644 index 0000000..31851d1 Binary files /dev/null and b/Old/3_grf_2d/Figure_1.png differ diff --git a/Old/3_grf_2d/gstest.py b/Old/3_grf_2d/gstest.py new file mode 100644 index 0000000..434603e --- /dev/null +++ b/Old/3_grf_2d/gstest.py @@ -0,0 +1,20 @@ +import numpy as np +import matplotlib.pyplot as plt + +import gstools as gs +x = y = np.linspace(0, 100, 100) +model = gs.TPLStable( + dim=2, # spatial dimension + var=1, # variance (C is calculated internally, so variance is actually 1) + len_low=0, # lower truncation of the power law + len_scale=10, # length scale (a.k.a. range), len_up = len_low + len_scale + nugget=0.1, # nugget + anis=0.5, # anisotropy between main direction and transversal ones + angles=np.pi / 4, # rotation angles + alpha=1.5, # shape parameter from the stable model + hurst=0.7, # hurst coefficient from the power law +) +srf = gs.SRF(model, mean=1.0, seed=19970221) +srf.structured([x, y]) +srf.plot() +plt.show() # blocks until you close the window \ No newline at end of file diff --git a/grf.py b/grf.py new file mode 100644 index 0000000..915ac9b --- /dev/null +++ b/grf.py @@ -0,0 +1,294 @@ +""" +grf.py – knihovna pro generování náhodných polí (GRF) +""" + +import numpy as np +from dataclasses import dataclass, field +import finufft +from scipy.fft import fft, ifft, fftfreq, fftshift + + +# ============================================================================= +# KORELAČNÍ TŘÍDY +# ============================================================================= +@dataclass +class GaussianCorrelation: + """Gaussian correlation model. + + Correlation function: + C(r) = sigma^2 * exp(-|r|^2 / (2 * phi^2)) + + Spectral density: + S(omega) = sigma^2 * (sqrt(2*pi) * phi)^dim * exp(-|omega|^2 * phi^2 / 2) + + Parameters + ---------- + L : float + Domain size [m]. Assumes a cubic domain [0, L]^dim. + N_freq : int + Number of frequency points per axis. Total N_freq^dim nodes in Fourier space. + phi : float + Correlation length [m]. Controls the spatial range of dependence. + sigma : float + Standard deviation of the field (square root of variance). Default 1.0. + dim : int + Spatial dimension (1 or 2). Default 1. + + Attributes + ---------- + f_points : np.ndarray, shape (N_freq,) + Centered angular frequencies [rad/m]: + omega_k = k * (2*pi / L) for k = -N_freq//2, ..., N_freq//2 - 1 + """ + L: float + N_freq: int + phi: float + sigma: float = 1.0 + dim: int = 1 + f_points: np.ndarray = field(init=False, repr=False) + + def __post_init__(self) -> None: + k = np.arange(-self.N_freq // 2, self.N_freq // 2) + self.f_points = k * (2 * np.pi / self.L) + + def spectral_density(self, omega_mag: np.ndarray) -> np.ndarray: + """Evaluate spectral density S(|omega|). + + Parameters + ---------- + omega_mag : np.ndarray, arbitrary shape + Magnitudes of angular frequencies [rad/m]. + + Returns + ------- + S : np.ndarray, same shape as omega_mag + """ + return ( + self.sigma ** 2 + * (np.sqrt(2 * np.pi) * self.phi) ** self.dim + * np.exp(-0.5 * (omega_mag * self.phi) ** 2) + ) + + +@dataclass +class ExponentialCorrelation: + """Exponential (Matérn-1/2) correlation model. + + Correlation function: + C(r) = sigma^2 * exp(-|r| / phi) + + Spectral density: + 1D: S(omega) = sigma^2 * 2*phi / (1 + (omega*phi)^2) + 2D: S(omega) = sigma^2 * 2*pi*phi^2 / (1 + (|omega|*phi)^2)^(3/2) + + Parameters + ---------- + L : float + Domain size [m]. Assumes a cubic domain [0, L]^dim. + N_freq : int + Number of frequency points per axis. Total N_freq^dim nodes in Fourier space. + phi : float + Correlation length [m]. Controls the spatial range of dependence. + sigma : float + Standard deviation of the field (square root of variance). Default 1.0. + dim : int + Spatial dimension (1 or 2). Default 1. + + Attributes + ---------- + f_points : np.ndarray, shape (N_freq,) + Centered angular frequencies [rad/m]: + omega_k = k * (2*pi / L) for k = -N_freq//2, ..., N_freq//2 - 1 + """ + + L: float + N_freq: int + phi: float + sigma: float = 1.0 + dim: int = 1 + f_points: np.ndarray = field(init=False, repr=False) + + def __post_init__(self) -> None: + k = np.arange(-self.N_freq // 2, self.N_freq // 2) + self.f_points = k * (2 * np.pi / self.L) + + def spectral_density(self, omega_mag: np.ndarray) -> np.ndarray: + """Evaluate spectral density S(|omega|). + + Parameters + ---------- + omega_mag : np.ndarray, arbitrary shape + Magnitudes of angular frequencies [rad/m]. + + Returns + ------- + S : np.ndarray, same shape as omega_mag + + Raises + ------ + NotImplementedError + For dim > 2. + """ + if self.dim == 1: + return self.sigma ** 2 * 2 * self.phi / (1 + (omega_mag * self.phi) ** 2) + elif self.dim == 2: + return ( + self.sigma ** 2 + * 2 * np.pi * self.phi ** 2 + / (1 + (omega_mag * self.phi) ** 2) ** 1.5 + ) + else: + raise NotImplementedError( + f"ExponentialCorrelation: dim={self.dim} not implemented (max dim=2)" + ) + + + + +# ============================================================================= +# HERMITOVSKÁ SYMETRIE f_hat[-k] = conj(f_hat[k]) → reálný IFFT výstup +# ============================================================================= + +def make_hermitian_nd(f_hat): + """nD hermitovská symetrie v centrovaném pořadí (DC uprostřed na indexu N//2).""" + N = f_hat.shape[0] + zi = N // 2 + f = f_hat.copy() + f[tuple([zi] * f.ndim)] = f[tuple([zi] * f.ndim)].real # DC reálná + for idx in np.ndindex(*f.shape): + shifted = tuple(i - zi for i in idx) + if all(s <= 0 for s in shifted): + continue + f[tuple((zi - s) % N for s in shifted)] = np.conj(f[idx]) + return f + + +def make_white_noise( + N_freq: int, + dim: int = 1, + seed: int | None = None, + use_gstools_rng: bool = False, +) -> np.ndarray: + """Generate complex white noise in Fourier space. + + Each component is an independent complex Gaussian variable with zero mean + and unit variance: z ~ CN(0, 1). + + Parameters + ---------- + N_freq : int + Number of frequency points per axis. + dim : int + Spatial dimension. Default 1. + seed : int or None + Random seed for reproducibility. None = random seed. + use_gstools_rng : bool + True - use gstools.random.RNG for direct comparison with GSTools. + False - use numpy.random.default_rng (default, recommended). + + Returns + ------- + noise : np.ndarray, shape (N_freq,) * dim, dtype complex128 + """ + size = N_freq ** dim + if use_gstools_rng: + from gstools.random import RNG + gs_rng = RNG(seed) + rs = gs_rng.random + flat = (rs.normal(size=size) + 1j * rs.normal(size=size)) / np.sqrt(2) + else: + rng = np.random.default_rng(seed) + flat = (rng.standard_normal(size) + 1j * rng.standard_normal(size)) / np.sqrt(2) + return flat.reshape((N_freq,) * dim) + + +# ============================================================================= +# GRF GENERATION +# ============================================================================= + +def generate_grf( + x_points: np.ndarray, + corr: GaussianCorrelation | ExponentialCorrelation, + weights: np.ndarray | None = None, + seed: int | None = None, +) -> np.ndarray: + """Generate a GRF realization at arbitrary points via NUFFT type 2. + + x_points : np.ndarray, shape (M,) for 1D or (M, 2) for 2D + corr : GaussianCorrelation | ExponentialCorrelation + weights : white noise, shape (N_freq,)^dim; None = generate from seed + + Returns field.real, shape (M,) + """ + dim = corr.dim + N_freq = corr.N_freq + L = corr.L + + if weights is None: + weights = make_white_noise(N_freq, dim=dim, seed=seed) + + # |omega| pro každý frekvenční bod + if dim == 1: + omega_mag = np.abs(corr.f_points) + else: + grids = np.meshgrid(*([corr.f_points] * dim), indexing='ij') + omega_mag = np.sqrt(sum(g**2 for g in grids)) + + S = corr.spectral_density(omega_mag) + f_hat = make_hermitian_nd(weights * np.sqrt(S)) + + x_points = np.asarray(x_points) + + if dim == 1: + x_nu = (x_points / L) * 2 * np.pi - np.pi # [0,L] -> [-pi,pi] + field = finufft.nufft1d2(x_nu, f_hat.astype(np.complex128)) + elif dim == 2: + x_nu = (x_points[:, 0] / L) * 2 * np.pi - np.pi + y_nu = (x_points[:, 1] / L) * 2 * np.pi - np.pi + field = finufft.nufft2d2(x_nu, y_nu, f_hat.astype(np.complex128)) + else: + raise NotImplementedError("dim > 2") + + # normalizace: C(0) = sigma² = (1/2π)*∫S(ω)dω ≈ (Δω/2π)*ΣS(ωk) = (1/L)*ΣS(ωk) + # -> faktor (1/L)^(dim/2) + norm = (1.0 / L) ** (dim / 2) + return (field * norm).real + +# ============================================================================= +# EMPIRICKÝ V`ARIOGRAM gamma(h) = 0.5 * E[(f(x)-f(x+h))^2] +# ============================================================================= + +def empirical_variogram(x, field, n_bins=30, n_sample=300): + """ + Estimate the empirical variogram from a single field realization. + Works for both 1D and 2D. Subsamples n_sample points to reduce O(M^2) cost. + + x : shape (M,) for 1D, (M, 2) for 2D + field : shape (M,) + """ + rng = np.random.default_rng(0) + x = np.asarray(x) + idx = rng.choice(len(field), size=min(len(field), n_sample), replace=False) + xi, fi = x[idx], field[idx] + + # vzdálenost – funguje pro 1D i 2D + if xi.ndim == 1: + dists = np.abs(xi[:, None] - xi[None, :]) # (n, n) + else: + diff = xi[:, None, :] - xi[None, :, :] # (n, n, 2) + dists = np.sqrt((diff ** 2).sum(axis=-1)) # (n, n) + + # jen horní trojúhelník (každý pár jednou) + i_upper, j_upper = np.triu_indices(len(xi), k=1) + pairs_h = dists[i_upper, j_upper] + pairs_sq = (fi[i_upper] - fi[j_upper]) ** 2 + + bins = np.linspace(0, pairs_h.max(), n_bins + 1) + h_vals, g_vals = [], [] + for k in range(n_bins): + mask = (pairs_h >= bins[k]) & (pairs_h < bins[k + 1]) + if mask.sum() > 5: + h_vals.append(0.5 * (bins[k] + bins[k + 1])) + g_vals.append(0.5 * pairs_sq[mask].mean()) + + return np.array(h_vals), np.array(g_vals) \ No newline at end of file diff --git a/gstoolsnufft.py b/gstoolsnufft.py new file mode 100644 index 0000000..d9ebf8d --- /dev/null +++ b/gstoolsnufft.py @@ -0,0 +1,71 @@ +""" +gstools_nufft3.py – Rekonstrukce GSTools pole pomocí NUFFT typ 3 + +GSTools (RandMeth) počítá: + f(x) = sqrt(1/N) * sum_k [ z1_k*cos(w_k*x) + z2_k*sin(w_k*x) ] + = Re[ sum_k c_k * exp(i*w_k*x) ] kde c_k = z1_k - i*z2_k + +finufft.nufft1d3: neuniformní frekvence s_k -> hodnoty na libovolných bodech x_j +""" + +import numpy as np +import matplotlib.pyplot as plt +import gstools as gs +import finufft + +L, phi, seed = 100.0, 5.0, 42 +N = 512 +x = np.linspace(0, L, N) + +# --- 1. Vyhodnotíme GSTools normálně --- +model = gs.Gaussian(dim=1, var=1.0, len_scale=phi) +srf = gs.SRF(model, seed=seed) +field_gs = np.asarray(srf(x)).ravel() + +# --- 2. Vytáhneme interní stav generátoru --- +gen = srf.generator +print("Atributy generátoru:", [a for a in dir(gen) if not a.startswith('__')]) + +z1 = gen._z_1 # shape (N_modes,) +z2 = gen._z_2 # shape (N_modes,) +modes = gen._cov_sample # shape (1, N_modes) pro 1D + +omega = modes[0] # frekvence w_k (1D) + +print(f"Počet módů: {len(omega)}") +print(f"Rozsah frekvencí: [{omega.min():.4f}, {omega.max():.4f}]") + +# --- 3. Rekonstrukce přes NUFFT typ 3 --- +# f(x_j) = Re[ sum_k c_k * exp(i * w_k * x_j) ] +# c_k = (z1_k - i*z2_k) / sqrt(N_modes) +N_modes = len(omega) +c = (z1 - 1j * z2) / np.sqrt(N_modes) + +# nufft1d3: vstupy jsou neuniformní frekvence a neuniformní body +# pozor: finufft očekává x v libovolných reálných hodnotách (ne nutně [-pi,pi]) +field_nufft3 = finufft.nufft1d3( + omega.astype(np.float64), # neuniformní frekvence s_k + c.astype(np.complex128), # komplexní váhy c_k + x.astype(np.float64), # body kde chceme hodnoty + eps=1e-9 +).real + +# --- 4. Porovnání --- +max_diff = np.max(np.abs(field_gs - field_nufft3)) +print(f"\nMax |GSTools - NUFFT3|: {max_diff:.2e}") + +fig, axes = plt.subplots(2, 1, figsize=(12, 7)) + +axes[0].plot(x, field_gs, lw=1.5, label="GSTools přímý výpočet") +axes[0].plot(x, field_nufft3, lw=1, label="NUFFT typ 3 rekonstrukce", alpha=0.8, ls='--') +axes[0].set_title(f"1D pole – GSTools vs NUFFT3 rekonstrukce (φ={phi})") +axes[0].set_xlabel("x [m]") +axes[0].legend(); axes[0].grid(True, alpha=0.3) + +axes[1].plot(x, field_gs - field_nufft3, lw=1, color='red') +axes[1].set_title(f"Rozdíl: max = {max_diff:.2e}") +axes[1].set_xlabel("x [m]"); axes[1].set_ylabel("GSTools − NUFFT3") +axes[1].grid(True, alpha=0.3) + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/main_fields.py b/main_fields.py new file mode 100644 index 0000000..1fc76d1 --- /dev/null +++ b/main_fields.py @@ -0,0 +1,153 @@ +""" +main_fields.py – Bod 3 + 4 + +Bod 3: GRF 1D – pravidelná a nepravidelná mřížka, body vyznačeny na ose +Bod 4: GRF 2D – pravidelná a nepravidelná mřížka +""" + +import numpy as np +import matplotlib.pyplot as plt +import finufft +from scipy.interpolate import griddata + +import grf + +L, N, phi, seed = 100.0, 512, 5.0, 42 +N_sparse = 50 # počet bodů řídké NUFFT mřížky +N_dense = 500 # počet bodů husté NUFFT mřížky +rng = np.random.default_rng(seed) + + +# ============================================================================= +# BOD 3: GRF 1D +# ============================================================================= + +x_reg = np.linspace(0, L, N) +x_irr = np.sort(rng.uniform(0, L, N // 2)) +white = grf.make_white_noise(N, dim=1, seed=seed) + +fig, axes = plt.subplots(2, 3, figsize=(18, 8)) +for row, (label, CorrClass) in enumerate([("Gaussovská", grf.GaussianCorrelation), + ("Exponenciální", grf.ExponentialCorrelation)]): + corr = CorrClass(L, N, phi, dim=1) + + f_r = grf.generate_grf(x_reg, corr, weights=white) + ax = axes[row][0] + ax.plot(x_reg, f_r, lw=1) + ax.plot(x_reg, np.full(N, f_r.min() - 0.3), + '|', color='steelblue', ms=6, label="body mřížky") + ax.set_title(f"FFT pravidelná – {label} φ={phi}") + ax.legend(fontsize=8); ax.grid(True, alpha=0.3) + + # reference: hustá pravidelná mřížka + x_ref = np.linspace(0, L, N) + f_ref = grf.generate_grf(x_ref, corr, weights=white) + + # řídká nepravidelná + x_sparse = np.sort(rng.uniform(0, L, N_sparse)) + f_sparse = grf.generate_grf(x_sparse, corr, weights=white) + # interpolace řídkých bodů na jemnou mřížku + from scipy.interpolate import interp1d + f_sparse_interp = interp1d(x_sparse, f_sparse, kind='cubic', + bounds_error=False, fill_value='extrapolate')(x_ref) + + ax = axes[row][1] + ax.plot(x_ref, f_sparse_interp, '-', color='darkorange', lw=1.5, alpha=0.4, + label=f"interp z {N_sparse} bodů") + ax.plot(x_ref, f_ref, 'k-', lw=1, label="reference") + ax.scatter(x_sparse, f_sparse, s=25, color='darkorange', zorder=3) + ax.plot(x_sparse, np.full(len(x_sparse), f_ref.min() - 0.3), + '|', color='darkorange', ms=8) + ax.set_title(f"NUFFT řídká ({N_sparse} bodů) – {label}") + ax.legend(fontsize=8); ax.grid(True, alpha=0.3) + + # hustá nepravidelná + x_dense = np.sort(rng.uniform(0, L, N_dense)) + f_dense = grf.generate_grf(x_dense, corr, weights=white) + f_dense_interp = interp1d(x_dense, f_dense, kind='cubic', + bounds_error=False, fill_value='extrapolate')(x_ref) + + ax = axes[row][2] + ax.plot(x_ref, f_dense_interp, '-', color='steelblue', lw=1.5, alpha=0.4, + label=f"interp z {N_dense} bodů") + ax.plot(x_ref, f_ref, 'k-', lw=1, label="reference") + ax.scatter(x_dense, f_dense, s=4, color='steelblue', zorder=3) + ax.plot(x_dense, np.full(len(x_dense), f_ref.min() - 0.3), + '|', color='steelblue', ms=4) + ax.set_title(f"NUFFT hustá ({N_dense} bodů) – {label}") + ax.legend(fontsize=8); ax.grid(True, alpha=0.3) + +plt.suptitle("Bod 3: GRF 1D") +plt.tight_layout(); plt.show() + + +# ============================================================================= +# BOD 4: GRF 2D +# ============================================================================= + +N2 = 64 +N2_sparse = 20 # řídká mřížka – počet bodů na osu (N2_sparse^2 celkem) +N2_dense = 64 # hustá mřížka – počet bodů na osu (N2_dense^2 celkem) + +white2d = grf.make_white_noise(N2, dim=2, seed=seed) +g = np.linspace(0, L, N2) +xx, yy = np.meshgrid(g, g) + +n_sparse = N2_sparse**2 +n_dense = N2_dense**2 +pts_sparse = np.column_stack([rng.uniform(0, L, n_sparse), rng.uniform(0, L, n_sparse)]) +pts_dense = np.column_stack([rng.uniform(0, L, n_dense), rng.uniform(0, L, n_dense)]) + +for label, CorrClass in [("Gaussovská", grf.GaussianCorrelation), + ("Exponenciální", grf.ExponentialCorrelation)]: + corr2d = CorrClass(L, N2, phi, dim=2) + + # 1. Pravidelná referenční mřížka (stejná) + pts_reg = np.column_stack([xx.ravel(), yy.ravel()]) + field_ref = grf.generate_grf(pts_reg, corr2d, weights=white2d).reshape(N2, N2) + + # 2. Generování z řídkých a hustých bodů (stejné) + f_sparse = grf.generate_grf(pts_sparse, corr2d, weights=white2d) + field_sparse_interp = griddata(pts_sparse, f_sparse, (xx, yy), method='linear') + + f_dense = grf.generate_grf(pts_dense, corr2d, weights=white2d) + field_dense_interp = griddata(pts_dense, f_dense, (xx, yy), method='linear') + + + # --- VIZUALIZACE (4 subploty, tečky všude, kam patří) --- + fig, axes = plt.subplots(1, 4, figsize=(22, 5)) # Širší plátno + vmin, vmax = field_ref.min(), field_ref.max() + + # Nastavení pro scatter-bar + sm = plt.cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(vmin=vmin, vmax=vmax)) + + # 1. Referenční pole + im0 = axes[0].imshow(field_ref, extent=[0,L,0,L], origin='lower', cmap='viridis', vmin=vmin, vmax=vmax) + axes[0].set_title(f"Referenční pravidelná ({N2}×{N2})") + plt.colorbar(im0, ax=axes[0]) + + # 2. NOVÝ: Samotná nepravidelná mřížka (Barevné tečky) + sc = axes[1].scatter(pts_sparse[:, 0], pts_sparse[:, 1], c=f_sparse, cmap='viridis', s=18, vmin=vmin, vmax=vmax) + axes[1].set_title(f"Nepravidelná mřížka\n({n_sparse} vzorků)") + axes[1].set_xlim(0, L); axes[1].set_ylim(0, L) + axes[1].set_aspect('equal') + axes[1].grid(True, alpha=0.2) # Jemná mřížka pro orientaci + plt.colorbar(sc, ax=axes[1]) + + # 3. PŮVODNÍ (s tečkami): Interpolace z řídké mřížky + im2 = axes[2].imshow(field_sparse_interp, extent=[0,L,0,L], origin='lower', cmap='viridis', vmin=vmin, vmax=vmax) + # Tady jsou ty bílé tečky navrch + axes[2].plot(pts_sparse[:, 0], pts_sparse[:, 1], '.', color='white', ms=2.5, alpha=0.5, label='vzorky') + axes[2].set_title(f"NUFFT řídká ({n_sparse} bodů) + interp") + plt.colorbar(im2, ax=axes[2]) + + # 4. PŮVODNÍ (s tečkami): Interpolace z husté mřížky + im3 = axes[3].imshow(field_dense_interp, extent=[0,L,0,L], origin='lower', cmap='viridis', vmin=vmin, vmax=vmax) + # Tady jsou ty bílé tečky navrch (menší ms= a alpha=) + axes[3].plot(pts_dense[:, 0], pts_dense[:, 1], '.', color='white', ms=1.8, alpha=0.4) + axes[3].set_title(f"NUFFT hustá ({n_dense} bodů) + interp") + plt.colorbar(im3, ax=axes[3]) + + plt.suptitle(f"Bod 4: GRF 2D – {label} φ={phi}") + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/main_gstools.py b/main_gstools.py new file mode 100644 index 0000000..f4c5bff --- /dev/null +++ b/main_gstools.py @@ -0,0 +1,212 @@ +""" +main_gstools.py – Bod 5: Porovnání naší NUFFT metody s GSTools + +Srovnání přes empirický variogram – obě metody by měly konvergovat +k teoretické křivce gamma(h) = sigma^2 * (1 - exp(-h^2 / 2*phi^2)) + +POZOR na parametrizaci GSTools Gaussian: + GSTools: C(r) = exp(-pi * r^2 / (4 * len_scale^2)) + Naše: C(r) = exp(-r^2 / (2 * phi^2)) + Převod: len_scale_gs = phi * sqrt(pi/2) ≈ phi * 1.2533 + +POZOR na spektrální rozlišení: + Frekvenční krok Δω = 2π/L – při velkém phi (phi > L/8) spadne + do spektrální hustoty jen pár frekvenčních bodů a naše NUFFT + metoda nebude správně odhadovat rozptyl. Řešení: zvýšit L. +""" + +import numpy as np +import matplotlib.pyplot as plt +import gstools as gs +from scipy.interpolate import griddata +from scipy.stats import norm as sp_norm +import grf + + +L, N, phi, seed = 100.0, 1024, 11.0, 42 + +# Převod phi -> len_scale pro GSTools Gaussian +# GSTools: C(r) = exp(-pi*r^2 / (4*ls^2)), naše: C(r) = exp(-r^2 / (2*phi^2)) +# => ls = phi * sqrt(pi/2) +len_scale_gs = phi * np.sqrt(np.pi / 2) + +# Varování při špatném poměru phi/L +delta_omega = 2 * np.pi / L +spectral_width = 1.0 / phi +n_spectral_pts = spectral_width / delta_omega +if n_spectral_pts < 3: + print(f"VAROVÁNÍ: phi={phi}, L={L} -> do spektra spadnou jen " + f"~{n_spectral_pts:.1f} frekvenční body. " + f"Zvyš L (doporučeno L > {int(8*phi)+1}) nebo sniž phi.") + +# ============================================================================= +# 1D – jedna realizace +# ============================================================================= + +x = np.linspace(0, L, N) +corr = grf.GaussianCorrelation(L, N, phi, sigma=1.0, dim=1) +model = gs.Gaussian(dim=1, var=1.0, len_scale=len_scale_gs) + +white = grf.make_white_noise(N, dim=1, seed=seed) +field_nufft = grf.generate_grf(x, corr, weights=white) +field_gs = np.asarray(gs.SRF(model, seed=seed)(x)).ravel() + +# ============================================================================= +# VARIOGRAM 1D – průměr přes N_real realizací +# ============================================================================= + +N_real = 100 +n_bins = 40 +# Oříznutí na L/3 – při větších lagy je málo párů a odhad je nespolehlivý +h_max = L / 3 +bin_edges = np.linspace(0, h_max, n_bins + 1) +bin_c = 0.5 * (bin_edges[:-1] + bin_edges[1:]) + +g_nufft = np.zeros(n_bins) +g_gs = np.zeros(n_bins) +cnt_nufft = np.zeros(n_bins, dtype=int) +cnt_gs = np.zeros(n_bins, dtype=int) + +for i in range(N_real): + w = grf.make_white_noise(N, dim=1, seed=i) + f_n = grf.generate_grf(x, corr, weights=w) + h_v, g = grf.empirical_variogram(x, f_n, n_bins=n_bins) + # empirical_variogram vrací variabilní délku – mapujeme do fixních binů + for k, hk in enumerate(h_v): + bi = np.searchsorted(bin_edges[1:], hk) + if 0 <= bi < n_bins: + g_nufft[bi] += g[k] + cnt_nufft[bi] += 1 + + f_g = np.asarray(gs.SRF(model, seed=i)(x)).ravel() + h_v, g = grf.empirical_variogram(x, f_g, n_bins=n_bins) + for k, hk in enumerate(h_v): + bi = np.searchsorted(bin_edges[1:], hk) + if 0 <= bi < n_bins: + g_gs[bi] += g[k] + cnt_gs[bi] += 1 + +g_nufft = np.divide(g_nufft, cnt_nufft, where=cnt_nufft > 0) +g_gs = np.divide(g_gs, cnt_gs, where=cnt_gs > 0) + +# Maska – zobraz jen biny kde máme data +mask_n = cnt_nufft > 0 +mask_g = cnt_gs > 0 + +# Teoretický variogram: C(r) = exp(-r^2 / 2*phi^2) => gamma(h) = 1 - C(h) +h_th = np.linspace(0, h_max, 300) +gamma_th = 1 - np.exp(-0.5 * (h_th / phi)**2) + +# ============================================================================= +# 2D – generování a variogram přímo z nepravidelných bodů +# ============================================================================= + +N2 = 128 +g2 = np.linspace(0, L, N2) +xx, yy = np.meshgrid(g2, g2) +rng_np = np.random.default_rng(seed) +pts_irr = np.column_stack([rng_np.uniform(0, L, N2**2), + rng_np.uniform(0, L, N2**2)]) + +corr2d = grf.GaussianCorrelation(L, N2, phi, sigma=1.0, dim=2) +white2d = grf.make_white_noise(N2, dim=2, seed=seed) + +# Generujeme pole v nepravidelných bodech +f_irr2d = grf.generate_grf(pts_irr, corr2d, weights=white2d) + +# ✅ Interpolace jen pro vizualizaci – variogram počítáme z původních bodů +field2d_nufft = griddata(pts_irr, f_irr2d, (xx, yy), method='linear') + +model2d = gs.Gaussian(dim=2, var=1.0, len_scale=len_scale_gs) +field2d_gs = gs.SRF(model2d, seed=seed).structured([g2, g2]) + +# ✅ Variogram 2D – počítáme z nepravidelných bodů, ne z interpolované mřížky +h_max_2d = L / 3 +n_bins_2d = 30 + +# NUFFT: přímo z pts_irr a f_irr2d +h_2d_nufft, g_2d_nufft = grf.empirical_variogram(pts_irr, f_irr2d, + n_bins=n_bins_2d) + +# GSTools: z pravidelné mřížky (jako (M,2) pole bodů) +pts_reg = np.column_stack([xx.ravel(), yy.ravel()]) +f_gs_reg = field2d_gs.ravel() +h_2d_gs, g_2d_gs = grf.empirical_variogram(pts_reg, f_gs_reg, + n_bins=n_bins_2d) + +# Oříznutí na L/3 +mask_2d_n = h_2d_nufft <= h_max_2d +mask_2d_g = h_2d_gs <= h_max_2d + +h_th_2d = np.linspace(0, h_max_2d, 300) +gamma_th_2d = 1 - np.exp(-0.5 * (h_th_2d / phi)**2) + +# ============================================================================= +# GRAFY +# ============================================================================= + +fig = plt.figure(figsize=(16, 9)) +layout = fig.add_gridspec(2, 3, hspace=0.45, wspace=0.32) + +# --- 1D realizace --- +ax0 = fig.add_subplot(layout[0, 0]) +ax0.plot(x, field_gs, lw=1.2, label="GSTools") +ax0.plot(x, field_nufft, lw=1, label="naše NUFFT", alpha=0.8) +ax0.set_title(f"1D realizace (φ={phi})") +ax0.set_xlabel("x [m]") +ax0.legend(fontsize=8); ax0.grid(True, alpha=0.3) + +# --- Variogram 1D průměr --- +ax1 = fig.add_subplot(layout[0, 1]) +ax1.plot(bin_c[mask_n], g_nufft[mask_n], 'o-', ms=3, + label=f"naše NUFFT (avg {N_real})") +ax1.plot(bin_c[mask_g], g_gs[mask_g], 's-', ms=3, + label=f"GSTools (avg {N_real})") +ax1.plot(h_th, gamma_th, 'k--', lw=2, label="teoretický") +ax1.set_title(f"Variogram 1D – průměr přes {N_real} realizací\n" + f"(zobrazeno do h = L/3 = {h_max:.0f} m)") +ax1.set_xlabel("h [m]"); ax1.set_ylabel("γ(h)") +ax1.set_xlim(0, h_max) +ax1.legend(fontsize=8); ax1.grid(True, alpha=0.3) + +# --- GSTools vario_estimate (1 realizace) --- +bin_edges_gs = np.linspace(0, h_max, n_bins + 1) +bin_c_gs, vario_gs_builtin = gs.vario_estimate([x], field_gs, + bin_edges=bin_edges_gs) +ax2 = fig.add_subplot(layout[0, 2]) +ax2.scatter(bin_c_gs, vario_gs_builtin, s=15, label="gs.vario_estimate") +ax2.plot(h_th, gamma_th, 'k--', lw=2, label="teoretický") +ax2.set_title(f"GSTools vario_estimate (1 realizace)\n" + f"(zobrazeno do h = L/3 = {h_max:.0f} m)") +ax2.set_xlabel("h [m]"); ax2.set_ylabel("γ(h)") +ax2.set_xlim(0, h_max) +ax2.legend(fontsize=8); ax2.grid(True, alpha=0.3) + +# --- 2D NUFFT (interpolovaná vizualizace) --- +ax3 = fig.add_subplot(layout[1, 0]) +im3 = ax3.imshow(field2d_nufft, extent=[0,L,0,L], origin='lower', cmap='viridis') +plt.colorbar(im3, ax=ax3) +ax3.set_title(f"2D naše NUFFT (φ={phi})") + +# --- 2D GSTools --- +ax4 = fig.add_subplot(layout[1, 1]) +im4 = ax4.imshow(field2d_gs, extent=[0,L,0,L], origin='lower', cmap='viridis') +plt.colorbar(im4, ax=ax4) +ax4.set_title(f"2D GSTools (φ={phi})") + +# --- Variogram 2D – z nepravidelných bodů --- +ax5 = fig.add_subplot(layout[1, 2]) +ax5.plot(h_2d_nufft[mask_2d_n], g_2d_nufft[mask_2d_n], 'o-', ms=3, + label="naše NUFFT (z irr. bodů)") +ax5.plot(h_2d_gs[mask_2d_g], g_2d_gs[mask_2d_g], 's-', ms=3, + label="GSTools (z reg. mřížky)") +ax5.plot(h_th_2d, gamma_th_2d, 'k--', lw=2, label="teoretický") +ax5.set_title(f"Variogram 2D (do L/3 = {h_max_2d:.0f} m)\n" + f"✅ počítáno z původních bodů, ne z interpolace") +ax5.set_xlabel("h [m]"); ax5.set_ylabel("γ(h)") +ax5.set_xlim(0, h_max_2d) +ax5.legend(fontsize=8); ax5.grid(True, alpha=0.3) + +plt.suptitle(f"Bod 5: Srovnání s GSTools – Gaussovská korelace φ={phi}", fontsize=13) +plt.tight_layout() +plt.show() \ No newline at end of file