diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..c9ebf2d --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:system" +} \ No newline at end of file diff --git a/generovani.py b/generovani.py deleted file mode 100644 index 8fe78c4..0000000 --- a/generovani.py +++ /dev/null @@ -1,109 +0,0 @@ -import marimo - -__generated_with = "0.19.11" -app = marimo.App(width="medium") - - -@app.cell -def _(): - import marimo as mo - - return (mo,) - - -@app.cell -def _(mo): - mo.md(""" - ### Struktura vstupu náhodných polí (xarray) - - Data jsou reprezentována pomocí datové struktury `xarray.Dataset`. - - **Definované dimenze a souřadnice (coords):** - - `i_point`: Index konkrétního bodu v nepravidelné síti. - - `i_sample`: Index vzorku (realizace) náhodného pole. - - `i_dim`: Označení prostorové osy (např. 'x', 'y'). - - **Datové proměnné (data_vars):** - - `X`: Matice prostorových souřadnic bodů. Má tvar `[i_dim, i_point]`. - - `QA`: První vygenerované náhodné pole. Má tvar `[i_point, i_sample]`. - - `QB`: Druhé vygenerované náhodné pole. Má tvar `[i_point, i_sample]`. - """) - return - - -@app.cell -def _(): - import xarray as xr - import numpy as np - import matplotlib.pyplot as plt - - n_points = 300 - n_samples = 100 - n_dim = 2 - - X_data = np.random.rand(n_dim, n_points) - QA_data = np.random.rand(n_points, n_samples) - QB_data = np.random.normal(loc=0.5, scale=0.15, size=(n_points, n_samples)) - - ds = xr.Dataset( - data_vars={ - "X": (("i_dim", "i_point"), X_data), - "QA": (("i_point", "i_sample"), QA_data), - "QB": (("i_point", "i_sample"), QB_data), - }, - coords={ - "i_point": np.arange(n_points), - "i_sample": np.arange(n_samples), - "i_dim": ["x", "y"] - } - ) - - # Přidání atributů s pojmenováním porovnávaných veličin - ds["QA"].attrs["long_name"] = "Porovnávaná veličina A (např. propustnost)" - ds["QA"].attrs["units"] = "-" - ds["QB"].attrs["long_name"] = "Porovnávaná veličina B (např. hustota puklin)" - ds["QB"].attrs["units"] = "-" - ds.attrs["description"] = "Testovací data náhodných polí ze samplů na nepravidelné síti" - - stats = xr.Dataset({ - "mean_QA": ds["QA"].mean(dim="i_sample"), - "var_QA": ds["QA"].var(dim="i_sample"), - "mean_QB": ds["QB"].mean(dim="i_sample"), - "var_QB": ds["QB"].var(dim="i_sample") - }) - - x_coords = ds["X"].sel(i_dim="x") - y_coords = ds["X"].sel(i_dim="y") - - fig, axes = plt.subplots(2, 3, figsize=(15, 9)) - fig.suptitle("Porovnání polí QA a QB na neregulární síti", fontsize=16) - - sc1 = axes[0, 0].scatter(x_coords, y_coords, c=stats["mean_QA"], cmap='viridis', s=20) - axes[0, 0].set_title("QA: Mapa průměru") - plt.colorbar(sc1, ax=axes[0, 0]) - - sc2 = axes[0, 1].scatter(x_coords, y_coords, c=stats["var_QA"], cmap='magma', s=20) - axes[0, 1].set_title("QA: Mapa rozptylu") - plt.colorbar(sc2, ax=axes[0, 1]) - - axes[0, 2].hist(ds["QA"].values.flatten(), bins=30, color='skyblue', edgecolor='black') - axes[0, 2].set_title("QA: Histogram všech hodnot") - - sc3 = axes[1, 0].scatter(x_coords, y_coords, c=stats["mean_QB"], cmap='viridis', s=20) - axes[1, 0].set_title("QB: Mapa průměru") - plt.colorbar(sc3, ax=axes[1, 0]) - - sc4 = axes[1, 1].scatter(x_coords, y_coords, c=stats["var_QB"], cmap='magma', s=20) - axes[1, 1].set_title("QB: Mapa rozptylu") - plt.colorbar(sc4, ax=axes[1, 1]) - - axes[1, 2].hist(ds["QB"].values.flatten(), bins=30, color='salmon', edgecolor='black') - axes[1, 2].set_title("QB: Histogram všech hodnot") - - plt.tight_layout() - plt.show() - return - - -if __name__ == "__main__": - app.run() diff --git a/vizualizace-xarray/analyza.py b/vizualizace-xarray/analyza.py new file mode 100644 index 0000000..95559c2 --- /dev/null +++ b/vizualizace-xarray/analyza.py @@ -0,0 +1,135 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import binned_statistic_2d, ttest_ind +from scipy.interpolate import griddata +import pywt +import warnings +warnings.filterwarnings('ignore') + +def create_binned_grid(x, y, values, num_bins=20): + bins = (num_bins, num_bins) + binned = binned_statistic_2d(x, y, values, statistic='mean', bins=bins) + return binned.statistic, binned.x_edge, binned.y_edge + +def bin_single_field(X_coords, data_var, num_bins=20): + x = X_coords.sel(i_dim="x").values + y = X_coords.sel(i_dim="y").values + grid, x_edges, y_edges = create_binned_grid(x, y, data_var.values, num_bins) + return grid, x_edges, y_edges + +def bin_all_samples(X_coords, data_var, num_bins=20): + x = X_coords.sel(i_dim="x").values + y = X_coords.sel(i_dim="y").values + n_samples = data_var.shape[1] + + grid_3d = np.zeros((num_bins, num_bins, n_samples)) + for i in range(n_samples): + grid_3d[:, :, i], x_edges, y_edges = create_binned_grid(x, y, data_var.values[:, i], num_bins) + + return grid_3d, x_edges, y_edges + +def plot_means_two_columns(x_edges, y_edges, grid_mean_A, grid_mean_B): + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + + im1 = axes[0].pcolormesh(x_edges, y_edges, grid_mean_A.T, cmap='viridis', shading='flat') + axes[0].set_title("Průměr pole A") + plt.colorbar(im1, ax=axes[0]) + + im2 = axes[1].pcolormesh(x_edges, y_edges, grid_mean_B.T, cmap='viridis', shading='flat') + axes[1].set_title("Průměr pole B") + plt.colorbar(im2, ax=axes[1]) + + plt.tight_layout() + return fig + +def perform_ttest(grid_A, grid_B): + t_stat, p_value = ttest_ind(grid_A, grid_B, axis=-1, equal_var=False, alternative='two-sided') + return t_stat, p_value + +def plot_ttest_results(x_edges, y_edges, t_stat, p_value): + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + + im1 = axes[0].pcolormesh(x_edges, y_edges, t_stat.T, cmap='coolwarm', shading='flat') + axes[0].set_title("T-statistika") + plt.colorbar(im1, ax=axes[0]) + + im2 = axes[1].pcolormesh(x_edges, y_edges, p_value.T, cmap='RdYlGn', shading='flat', vmin=0, vmax=0.1) + axes[1].set_title("p-hodnota") + plt.colorbar(im2, ax=axes[1]) + + plt.tight_layout() + return fig + +def multiscale_analysis(x, y, q_a_vals, q_b_vals): + grid_sizes = [2, 4, 8] + p_values_all = {'arith': [], 'geom': [], 'harm': []} + + fig1, axes1 = plt.subplots(len(grid_sizes), 3, figsize=(15, 5 * len(grid_sizes))) + + for idx, gs in enumerate(grid_sizes): + dx = 1.0 / gs + dy = 1.0 / gs + + p_map_arith = np.full((gs, gs), np.nan) + p_map_geom = np.full((gs, gs), np.nan) + p_map_harm = np.full((gs, gs), np.nan) + + for i in range(gs): + for j in range(gs): + mask = (x >= i*dx) & (x < (i+1)*dx) & (y >= j*dy) & (y < (j+1)*dy) + + if np.sum(mask) > 0: + qa_w = q_a_vals[mask, :] + qb_w = q_b_vals[mask, :] + + qa_arith = np.mean(qa_w, axis=0) + qb_arith = np.mean(qb_w, axis=0) + + qa_w_safe = np.clip(qa_w, 1e-10, None) + qb_w_safe = np.clip(qb_w, 1e-10, None) + + qa_geom = np.exp(np.mean(np.log(qa_w_safe), axis=0)) + qb_geom = np.exp(np.mean(np.log(qb_w_safe), axis=0)) + + qa_harm = 1.0 / np.mean(1.0 / qa_w_safe, axis=0) + qb_harm = 1.0 / np.mean(1.0 / qb_w_safe, axis=0) + + _, p_a = ttest_ind(qa_arith, qb_arith, equal_var=False) + _, p_g = ttest_ind(qa_geom, qb_geom, equal_var=False) + _, p_h = ttest_ind(qa_harm, qb_harm, equal_var=False) + + p_map_arith[j, i] = p_a + p_map_geom[j, i] = p_g + p_map_harm[j, i] = p_h + + p_values_all['arith'].append(np.nanmean(p_map_arith)) + p_values_all['geom'].append(np.nanmean(p_map_geom)) + p_values_all['harm'].append(np.nanmean(p_map_harm)) + + im_a = axes1[idx, 0].imshow(p_map_arith, origin='lower', extent=[0, 1, 0, 1], cmap='RdYlGn', vmin=0, vmax=1) + axes1[idx, 0].set_ylabel(f"Mřížka {gs}x{gs}") + if idx == 0: axes1[idx, 0].set_title("Aritmetický") + plt.colorbar(im_a, ax=axes1[idx, 0]) + + im_g = axes1[idx, 1].imshow(p_map_geom, origin='lower', extent=[0, 1, 0, 1], cmap='RdYlGn', vmin=0, vmax=1) + if idx == 0: axes1[idx, 1].set_title("Geometrický") + plt.colorbar(im_g, ax=axes1[idx, 1]) + + im_h = axes1[idx, 2].imshow(p_map_harm, origin='lower', extent=[0, 1, 0, 1], cmap='RdYlGn', vmin=0, vmax=1) + if idx == 0: axes1[idx, 2].set_title("Harmonický") + plt.colorbar(im_h, ax=axes1[idx, 2]) + + plt.tight_layout() + + fig2, ax2 = plt.subplots(figsize=(8, 5)) + ax2.plot(grid_sizes, p_values_all['arith'], marker='o', label='Aritmetický') + ax2.plot(grid_sizes, p_values_all['geom'], marker='s', label='Geometrický') + ax2.plot(grid_sizes, p_values_all['harm'], marker='^', label='Harmonický') + ax2.set_xticks(grid_sizes) + ax2.set_xlabel("Velikost okna (počet dělení)") + ax2.set_ylabel("Průměrná p-hodnota") + ax2.set_title("Závislost p-hodnoty na velikosti okna") + plt.legend() + + return fig1, fig2 + diff --git a/vizualizace-xarray/generovani.py b/vizualizace-xarray/generovani.py index 8fe78c4..eb846cf 100644 --- a/vizualizace-xarray/generovani.py +++ b/vizualizace-xarray/generovani.py @@ -1,17 +1,17 @@ import marimo -__generated_with = "0.19.11" +__generated_with = "0.20.2" app = marimo.App(width="medium") -@app.cell +@app.cell(hide_code=True) def _(): import marimo as mo return (mo,) -@app.cell +@app.cell(hide_code=True) def _(mo): mo.md(""" ### Struktura vstupu náhodných polí (xarray) @@ -22,86 +22,99 @@ def _(mo): - `i_point`: Index konkrétního bodu v nepravidelné síti. - `i_sample`: Index vzorku (realizace) náhodného pole. - `i_dim`: Označení prostorové osy (např. 'x', 'y'). - - **Datové proměnné (data_vars):** - - `X`: Matice prostorových souřadnic bodů. Má tvar `[i_dim, i_point]`. - - `QA`: První vygenerované náhodné pole. Má tvar `[i_point, i_sample]`. - - `QB`: Druhé vygenerované náhodné pole. Má tvar `[i_point, i_sample]`. """) return -@app.cell +@app.cell(hide_code=True) def _(): import xarray as xr import numpy as np - import matplotlib.pyplot as plt + import analyza - n_points = 300 - n_samples = 100 + n_points = 1000 + n_samples_A = 10 + n_samples_B = 20 n_dim = 2 X_data = np.random.rand(n_dim, n_points) - QA_data = np.random.rand(n_points, n_samples) - QB_data = np.random.normal(loc=0.5, scale=0.15, size=(n_points, n_samples)) + + # Testovací pole = 10**( np.random.normal(N, loc = -10, scale=3)) + QA_data = 10**(np.random.normal(loc=-10, scale=3, size=(n_points, n_samples_A))) + + # případně zvětšit scale pro větší rozdíl průměrů polí + QB_data = 10**(np.random.normal(loc=-10, scale=5, size=(n_points, n_samples_B))) ds = xr.Dataset( data_vars={ "X": (("i_dim", "i_point"), X_data), - "QA": (("i_point", "i_sample"), QA_data), - "QB": (("i_point", "i_sample"), QB_data), + "QA": (("i_point", "i_sample_A"), QA_data), + "QB": (("i_point", "i_sample_B"), QB_data), }, coords={ "i_point": np.arange(n_points), - "i_sample": np.arange(n_samples), + "i_sample_A": np.arange(n_samples_A), + "i_sample_B": np.arange(n_samples_B), "i_dim": ["x", "y"] } ) + return analyza, ds - # Přidání atributů s pojmenováním porovnávaných veličin - ds["QA"].attrs["long_name"] = "Porovnávaná veličina A (např. propustnost)" - ds["QA"].attrs["units"] = "-" - ds["QB"].attrs["long_name"] = "Porovnávaná veličina B (např. hustota puklin)" - ds["QB"].attrs["units"] = "-" - ds.attrs["description"] = "Testovací data náhodných polí ze samplů na nepravidelné síti" - stats = xr.Dataset({ - "mean_QA": ds["QA"].mean(dim="i_sample"), - "var_QA": ds["QA"].var(dim="i_sample"), - "mean_QB": ds["QB"].mean(dim="i_sample"), - "var_QB": ds["QB"].var(dim="i_sample") - }) +@app.cell(hide_code=True) +def _(ds): + # Průměry pro každé ze dvou polí. Průměry do samostatné buňky. + mean_QA = ds["QA"].mean(dim="i_sample_A") + mean_QB = ds["QB"].mean(dim="i_sample_B") + return mean_QA, mean_QB - x_coords = ds["X"].sel(i_dim="x") - y_coords = ds["X"].sel(i_dim="y") - fig, axes = plt.subplots(2, 3, figsize=(15, 9)) - fig.suptitle("Porovnání polí QA a QB na neregulární síti", fontsize=16) +@app.cell(hide_code=True) +def _(analyza, ds, mean_QA, mean_QB): + grid_mean_A, x_edges, y_edges = analyza.bin_single_field(ds["X"], mean_QA) + grid_mean_B, _, _ = analyza.bin_single_field(ds["X"], mean_QB) - sc1 = axes[0, 0].scatter(x_coords, y_coords, c=stats["mean_QA"], cmap='viridis', s=20) - axes[0, 0].set_title("QA: Mapa průměru") - plt.colorbar(sc1, ax=axes[0, 0]) + fig_means = analyza.plot_means_two_columns(x_edges, y_edges, grid_mean_A, grid_mean_B) + fig_means + return x_edges, y_edges - sc2 = axes[0, 1].scatter(x_coords, y_coords, c=stats["var_QA"], cmap='magma', s=20) - axes[0, 1].set_title("QA: Mapa rozptylu") - plt.colorbar(sc2, ax=axes[0, 1]) - axes[0, 2].hist(ds["QA"].values.flatten(), bins=30, color='skyblue', edgecolor='black') - axes[0, 2].set_title("QA: Histogram všech hodnot") +@app.cell(hide_code=True) +def _(analyza, ds): + grid_A, _, _ = analyza.bin_all_samples(ds["X"], ds["QA"]) + grid_B, _, _ = analyza.bin_all_samples(ds["X"], ds["QB"]) - sc3 = axes[1, 0].scatter(x_coords, y_coords, c=stats["mean_QB"], cmap='viridis', s=20) - axes[1, 0].set_title("QB: Mapa průměru") - plt.colorbar(sc3, ax=axes[1, 0]) + t_stat, p_value = analyza.perform_ttest(grid_A, grid_B) + return p_value, t_stat + + +@app.cell(hide_code=True) +def _(analyza, p_value, t_stat, x_edges, y_edges): + fig_test = analyza.plot_ttest_results(x_edges, y_edges, t_stat, p_value) + fig_test + return - sc4 = axes[1, 1].scatter(x_coords, y_coords, c=stats["var_QB"], cmap='magma', s=20) - axes[1, 1].set_title("QB: Mapa rozptylu") - plt.colorbar(sc4, ax=axes[1, 1]) - axes[1, 2].hist(ds["QB"].values.flatten(), bins=30, color='salmon', edgecolor='black') - axes[1, 2].set_title("QB: Histogram všech hodnot") +@app.cell(hide_code=True) +def _(analyza, ds): + x_vals = ds["X"].sel(i_dim="x").values + y_vals = ds["X"].sel(i_dim="y").values + q_a_vals = ds["QA"].values + q_b_vals = ds["QB"].values + + fig_multiscale_maps, fig_multiscale_line = analyza.multiscale_analysis(x_vals, y_vals, q_a_vals, q_b_vals) + return fig_multiscale_line, fig_multiscale_maps + + +@app.cell(hide_code=True) +def _(fig_multiscale_maps): + fig_multiscale_maps + return + - plt.tight_layout() - plt.show() +@app.cell(hide_code=True) +def _(fig_multiscale_line): + fig_multiscale_line return diff --git a/vizualizace-xarray/project/generovani.py b/vizualizace-xarray/project/generovani.py deleted file mode 100644 index 8d9f7dc..0000000 --- a/vizualizace-xarray/project/generovani.py +++ /dev/null @@ -1,76 +0,0 @@ -import xarray as xr -import numpy as np -import matplotlib.pyplot as plt - -# --- Konfigurační parametry --- -n_points = 300 # Počet bodů v prostoru -n_samples = 100 # Počet realizací (vzorků) -n_dim = 2 # Rozměrnost prostoru (x, y) - -# --- Generování syntetických dat --- -# Náhodné souřadnice bodů v jednotkovém čtverci [0, 1] x [0, 1] -X_data = np.random.rand(n_dim, n_points) - -# Pole QA: generováno z rovnoměrného rozdělení [0, 1] -QA_data = np.random.rand(n_points, n_samples) - -# Pole QB: generováno z normálního rozdělení (μ=0.5, σ=0.15) -QB_data = np.random.normal(loc=0.5, scale=0.15, size=(n_points, n_samples)) - -# --- Vytvoření struktury xarray Dataset --- -ds = xr.Dataset( - data_vars={ - "X": (("i_dim", "i_point"), X_data), - "QA": (("i_point", "i_sample"), QA_data), - "QB": (("i_point", "i_sample"), QB_data), - }, - coords={ - "i_point": np.arange(n_points), - "i_sample": np.arange(n_samples), - "i_dim": ["x", "y"] - } -) - -# --- Statistická analýza --- -# Výpočet průměru a rozptylu podél dimenze realizací (i_sample) -stats = xr.Dataset({ - "mean_QA": ds["QA"].mean(dim="i_sample"), - "var_QA": ds["QA"].var(dim="i_sample"), - "mean_QB": ds["QB"].mean(dim="i_sample"), - "var_QB": ds["QB"].var(dim="i_sample") -}) - -# Extraxe souřadnic pro vizualizaci -x_coords = ds["X"].sel(i_dim="x") -y_coords = ds["X"].sel(i_dim="y") - -# --- Vizualizace --- -fig, axes = plt.subplots(2, 3, figsize=(15, 9)) -fig.suptitle("Porovnání polí QA a QB na neregulární síti", fontsize=16) - -# Horní řádek: Analýza pole QA -sc1 = axes[0, 0].scatter(x_coords, y_coords, c=stats["mean_QA"], cmap='viridis', s=20) -axes[0, 0].set_title("QA: Mapa průměru") -plt.colorbar(sc1, ax=axes[0, 0]) - -sc2 = axes[0, 1].scatter(x_coords, y_coords, c=stats["var_QA"], cmap='magma', s=20) -axes[0, 1].set_title("QA: Mapa rozptylu") -plt.colorbar(sc2, ax=axes[0, 1]) - -axes[0, 2].hist(ds["QA"].values.flatten(), bins=30, color='skyblue', edgecolor='black') -axes[0, 2].set_title("QA: Histogram všech hodnot") - -# Spodní řádek: Analýza pole QB -sc3 = axes[1, 0].scatter(x_coords, y_coords, c=stats["mean_QB"], cmap='viridis', s=20) -axes[1, 0].set_title("QB: Mapa průměru") -plt.colorbar(sc3, ax=axes[1, 0]) - -sc4 = axes[1, 1].scatter(x_coords, y_coords, c=stats["var_QB"], cmap='magma', s=20) -axes[1, 1].set_title("QB: Mapa rozptylu") -plt.colorbar(sc4, ax=axes[1, 1]) - -axes[1, 2].hist(ds["QB"].values.flatten(), bins=30, color='salmon', edgecolor='black') -axes[1, 2].set_title("QB: Histogram všech hodnot") - -plt.tight_layout() -plt.show() \ No newline at end of file