9.3 zadani

Alevember · Alevember · commit fc49770a94c8 · 2026-03-24T10:44:40.000+01:00
diff --git a/vizualizace-xarray/analyza.py 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
@@ -7,6 +7,7 @@
 @app.cell(hide_code=True)
 def _():
     import marimo as mo
+
     return (mo,)
 
 
@@ -21,11 +22,6 @@ 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
 
@@ -34,132 +30,93 @@ def _(mo):
 def _():
     import xarray as xr
     import numpy as np
-    import matplotlib.pyplot as plt
-    from scipy.stats import binned_statistic_2d
+    import analyza
 
     n_points = 1000
-    n_samples = 100
+    n_samples_A = 10
+    n_samples_B = 20
     n_dim = 2
 
     X_data = np.random.rand(n_dim, n_points)
-    QA_data = np.random.uniform(0.01, 1.0, size=(n_points, n_samples))
-    QB_data = np.random.normal(loc=0.5, scale=0.15, size=(n_points, n_samples))
-    QB_data = np.clip(QB_data, 0.01, None)
+
+    # 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
+
+
+@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
+
+
+@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)
+
+    fig_means = analyza.plot_means_two_columns(x_edges, y_edges, grid_mean_A, grid_mean_B)
+    fig_means
+    return x_edges, y_edges
+
+
+@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"])
+
+    t_stat, p_value = analyza.perform_ttest(grid_A, grid_B)
+    return p_value, t_stat
 
-    ds["QA"].attrs["long_name"] = "Porovnávaná veličina A"
-    ds["QA"].attrs["units"] = "-"
-    ds["QB"].attrs["long_name"] = "Porovnávaná veličina B"
-    ds["QB"].attrs["units"] = "-"
-
-    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 = ds["X"].sel(i_dim="x").values
-    y = ds["X"].sel(i_dim="y").values
-    
-    # Pro grafy průměru a rozptylu použít imgshow nebo podobnou funkci
-    num_bins = 20
-    bins = (num_bins, num_bins)
-
-    binned_mean_QA = binned_statistic_2d(x, y, stats["mean_QA"].values, statistic='mean', bins=bins)
-    binned_var_QA = binned_statistic_2d(x, y, stats["var_QA"].values, statistic='mean', bins=bins)
-    binned_mean_QB = binned_statistic_2d(x, y, stats["mean_QB"].values, statistic='mean', bins=bins)
-    binned_var_QB = binned_statistic_2d(x, y, stats["var_QB"].values, statistic='mean', bins=bins)
-
-    x_edges = binned_mean_QA.x_edge
-    y_edges = binned_mean_QA.y_edge
-
-    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
-    
-    im1 = axes[0, 0].pcolormesh(x_edges, y_edges, binned_mean_QA.statistic.T, cmap='viridis', shading='flat')
-    axes[0, 0].set_title("QA: Mapa průměru")
-    axes[0, 0].set_aspect('equal')
-    plt.colorbar(im1, ax=axes[0, 0])
-
-    im2 = axes[0, 1].pcolormesh(x_edges, y_edges, binned_var_QA.statistic.T, cmap='magma', shading='flat')
-    axes[0, 1].set_title("QA: Mapa rozptylu")
-    axes[0, 1].set_aspect('equal')
-    plt.colorbar(im2, 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")
-
-    im3 = axes[1, 0].pcolormesh(x_edges, y_edges, binned_mean_QB.statistic.T, cmap='viridis', shading='flat')
-    axes[1, 0].set_title("QB: Mapa průměru")
-    axes[1, 0].set_aspect('equal')
-    plt.colorbar(im3, ax=axes[1, 0])
-
-    im4 = axes[1, 1].pcolormesh(x_edges, y_edges, binned_var_QB.statistic.T, cmap='magma', shading='flat')
-    axes[1, 1].set_title("QB: Mapa rozptylu")
-    axes[1, 1].set_aspect('equal')
-    plt.colorbar(im4, 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")
-
-    plt.tight_layout()
-    plt.show()
-
-    # Vypočtěte různé druhy průměrů (aritmetický, geometrický, harmonický) na podčtvercích (multi-scale) pro různé velikosti oken. Zatím pro jednu velikost okna.
-    q_values = stats["mean_QA"].values
-    grid_size = 10
-    A_arith = np.zeros((grid_size, grid_size))
-    A_geom = np.zeros((grid_size, grid_size))
-    A_harm = np.zeros((grid_size, grid_size))
-    
-    dx = 1.0 / grid_size
-    dy = 1.0 / grid_size
-    
-    for i in range(grid_size):
-        for j in range(grid_size):
-            mask = (x >= i*dx) & (x < (i+1)*dx) & (y >= j*dy) & (y < (j+1)*dy)
-            q_in_window = q_values[mask]
-            
-            if len(q_in_window) > 0:
-                A_arith[j, i] = np.mean(q_in_window)
-                A_geom[j, i] = np.exp(np.mean(np.log(q_in_window)))
-                A_harm[j, i] = 1.0 / np.mean(1.0 / q_in_window)
-            else:
-                A_arith[j, i] = np.nan
-                A_geom[j, i] = np.nan
-                A_harm[j, i] = np.nan
-
-    fig2, axes2 = plt.subplots(1, 3, figsize=(15, 5))
-    
-    im5 = axes2[0].imshow(A_arith, origin='lower', extent=[0, 1, 0, 1], cmap='viridis')
-    axes2[0].set_title("Aritmetický průměr")
-    plt.colorbar(im5, ax=axes2[0])
-    
-    im6 = axes2[1].imshow(A_geom, origin='lower', extent=[0, 1, 0, 1], cmap='viridis')
-    axes2[1].set_title("Geometrický průměr")
-    plt.colorbar(im6, ax=axes2[1])
-    
-    im7 = axes2[2].imshow(A_harm, origin='lower', extent=[0, 1, 0, 1], cmap='viridis')
-    axes2[2].set_title("Harmonický průměr")
-    plt.colorbar(im7, ax=axes2[2])
-    
-    plt.tight_layout()
-    plt.show()
-
-    return ds, fig, fig2
+
+@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
+
+
+@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
+
+
+@app.cell(hide_code=True)
+def _(fig_multiscale_line):
+    fig_multiscale_line
+    return
 
 
 if __name__ == "__main__":
-    app.run()
+    app.run()
diff --git a/vizualizace-xarray/project/generovani.py b/vizualizace-xarray/project/generovani.py
