Merge branch 'dev/normalization' into dev/general

Author: JorisVincent

docs/normalization/explore_FLODOG_parameters.py

@@ -1,3 +1,20 @@
+# -*- coding: utf-8 -*-
+# ---
+# jupyter:
+#   jupytext:
+#     cell_metadata_filter: sdmix,incorrectly_encoded_metadata,spatial_window_scalar,size,title,-all
+#     custom_cell_magics: kql
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: multyscale
+#     language: python
+#     name: python3
+# ---
+
 # %% [markdown]
 # # Exploring FLODOG normalization parameter
 # This guide describes how to change the normalization parameters
@@ -86,7 +103,7 @@
 plt.show()
 
 # %% Output default model
-normalized_4_05 = FLODOG.normalize_outputs(filters_output)
+normalized_4_05 = FLODOG.normalize_outputs(filters_output, eps=1e-6)
 output_4_05 = np.sum(normalized_4_05, axis=(0, 1))
 
 # %% [markdown]
@@ -127,7 +144,7 @@
 # FLODOG.normalization_weights = normalization_weights
 
 # %% Determine normalizing coefficients
-normalizing_coefficients_3 = multyscale.normalization.normalizers(
+normalizing_coefficients_3 = multyscale.normalization.norm_coeffs(
     filters_output, normalization_weights
 )
 
@@ -153,11 +170,11 @@
 # (the $\sigma$ of the 2D Gaussian)
 # to the spatial scale of the filter being normalized.
 # This is controlled by a parameter `spatial_window_scalar`:
-# $\sigma = \mathrm{spatial_window_scalar} \times s$.
+# $\sigma = \texttt{spatial window scalar} \times s$.
 # To explore how, here we apply various parameterizations of this normalization
 # to the same stimulus image and corresponding (weighted) filter outputs.
 
-# %% FLODOG with spatial_window_scalar = 2
+# %% FLODOG with spatial_window_scalar=2
 # Set parameter
 spatial_window_scalar = 2
 # FLODOG.spatial_window_scalar = spaial_window_scalar
@@ -169,7 +186,7 @@
 FLODOG.window_sigmas = window_sigmas
 
 # Apply spatial averaging windows to normalizing coefficients
-energies_2_3 = FLODOG.normalizers_to_RMS(normalizing_coefficients_3)
+energies_2_3 = FLODOG.norm_energies(normalizing_coefficients_3, eps=1e-6)
 
 # Visualize each energy estimate
 fig, axs = plt.subplots(*energies_2_3.shape[:2], sharex="all", sharey="all")
@@ -186,7 +203,7 @@
 plt.show()
 
 # %% Output
-normalized_2_3 = filters_output / energies_2_3
+normalized_2_3 = filters_output / (energies_2_3 + 1e-6)
 output_2_3 = np.sum(normalized_2_3, axis=(0, 1))
 
 # %% Comparing FLODOG outputs
@@ -233,7 +250,7 @@
 # to normalize the filter outputs accordingly.
 
 # %% Normalize
-normalized_4_3 = FLODOG_4_3.normalize_outputs(filters_output)
+normalized_4_3 = FLODOG_4_3.normalize_outputs(filters_output, eps=1e-6)
 
 # %% [markdown]
 # To then readout the final model prediction,
@@ -291,11 +308,11 @@
 # )
 
 # %% Outputs
-output_2_05 = FLODOG_2_05.normalize_outputs(filters_output).sum((0, 1))
+output_2_05 = FLODOG_2_05.normalize_outputs(filters_output, eps=1e-6).sum((0, 1))
 
 # %% LODOG, $\sigma=4$
 LODOG = multyscale.models.LODOG_RHS2007(shape=stimulus.shape, visextent=visextent, window_sigma=4)
-output_LODOG_4 = LODOG.normalize_outputs(filters_output).sum((0, 1))
+output_LODOG_4 = LODOG.normalize_outputs(filters_output, eps=1e-6).sum((0, 1))
 
 
 # %% [markdown]
@@ -312,7 +329,7 @@
 # - subtract the two means
 # - divide the whole model output by this difference
 
-# %% Standardize such that effect size = 1
+# %% Standardize such that effect size=1
 mask = np.load("example_stimulus_mask.npy")
 
 

docs/normalization/explore_LODOG_parameter.py

@@ -109,7 +109,7 @@
 # and thus also only need to be calculated once for all the parameterizations explored here.
 
 # %%
-normalizing_coefficients = LODOG.normalizers(filters_output)
+normalizing_coefficients = LODOG.norm_coeffs(filters_output)
 
 # %% Gaussian spatial averaging window [markdown]
 # The spatial averaging window in the LODOG model is
@@ -128,9 +128,11 @@
 
 assert np.all(window_sigmas == window_sigma)
 
-spatial_filters = multyscale.normalization.spatial_avg_windows_gaussian(
-    LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas
-)
+spatial_filters = np.ndarray(filters_output.shape)
+for o, s in np.ndindex(LODOG.window_sigmas.shape[:2]):
+    spatial_filters[o, s, :] = multyscale.normalization.spatial_kernel_gaussian(
+        LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas[o, s, :]
+    )
 plt.subplot(1, 2, 1)
 plt.imshow(spatial_filters[0, 0, ...], cmap="coolwarm", extent=visextent)
 plt.colorbar()
@@ -147,13 +149,13 @@
 plt.show()
 
 # %% [markdown]
-# The model method `.normalizers_to_RMS()` automatically generates and applies
+# The model method `.norm_energies()` automatically generates and applies
 # these spatial averaging windows to the proved normalizing coeffiencts.
 # This produces an $O \times S$ set of locally ($Y \times X$) calculated energies,
 # one for each normalizaing coefficient.
 
 # %%
-energies_4 = LODOG.normalizers_to_RMS(normalizing_coefficients)
+energies_4 = LODOG.norm_energies(normalizing_coefficients, eps=1e-6)
 
 # Visualize each local energy
 fig, axs = plt.subplots(*energies_4.shape[:2], sharex="all", sharey="all")
@@ -203,10 +205,14 @@
 assert np.all(LODOG.window_sigmas == LODOG.window_sigma)
 
 # %%
-spatial_filters = multyscale.normalization.spatial_avg_windows_gaussian(
-    LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas
-)
-plt.subplot(1, 2, 1)
+spatial_filters = np.ndarray(filters_output.shape)
+for o, s in np.ndindex(LODOG.window_sigmas.shape[:2]):
+    spatial_filters[o, s, :] = multyscale.normalization.spatial_kernel_gaussian(
+        LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas[o, s, :]
+    )
+
+# %% 
+plt.subplot(1, 2, 1) 
 plt.imshow(spatial_filters[0, 0, ...], cmap="coolwarm", extent=visextent)
 plt.colorbar()
 plt.subplot(1, 2, 2)
@@ -226,7 +232,7 @@
 # more locally restricted (smaller spatial averaging window) energies
 
 # %%
-energies_1 = LODOG.normalizers_to_RMS(normalizing_coefficients)
+energies_1 = LODOG.norm_energies(normalizing_coefficients, eps=1e-6)
 
 # Visualize each local energy
 fig, axs = plt.subplots(*energies_4.shape[:2], sharex="all", sharey="all")
@@ -242,7 +248,7 @@
 # as would be expected.
 
 # %%
-normalized_outputs_1 = filters_output / energies_1
+normalized_outputs_1 = filters_output / (energies_1 + 1e-6)
 
 # Visualize each normalized output
 vmin = min(np.min(normalized_outputs_4), np.min(normalized_outputs_1))

docs/normalization/normalization_FLODOG.py

@@ -165,8 +165,8 @@
 # the filter outputs at these spatial scales.
 
 # %% Normalizing coefficients
-normalizing_coefficients_LODOG = LODOG.normalizers(filters_output)
-normalizing_coefficients_FLODOG = FLODOG.normalizers(filters_output)
+normalizing_coefficients_LODOG = LODOG.norm_coeffs(filters_output)
+normalizing_coefficients_FLODOG = FLODOG.norm_coeffs(filters_output)
 
 # Visualize each norm. coeff.
 vmin = min(np.min(normalizing_coefficients_LODOG), np.min(normalizing_coefficients_FLODOG))
@@ -215,12 +215,17 @@
 # both in space, and in _spatial scale_ / _**F**requency_.
 
 # %% Spatial averaging window
-spatial_windows_LODOG = multyscale.normalization.spatial_avg_windows_gaussian(
-    LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas
-)
-spatial_windows_FLODOG = multyscale.normalization.spatial_avg_windows_gaussian(
-    FLODOG.bank.x, FLODOG.bank.y, FLODOG.window_sigmas
-)
+spatial_windows_LODOG = np.ndarray(filters_output.shape)
+for o, s in np.ndindex(LODOG.window_sigmas.shape[:2]):
+    spatial_windows_LODOG[o, s, :] = multyscale.normalization.spatial_kernel_gaussian(
+        LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas[o, s, :]
+    )
+
+spatial_windows_FLODOG = np.ndarray(filters_output.shape)
+for o, s in np.ndindex(FLODOG.window_sigmas.shape[:2]):
+    spatial_windows_FLODOG[o, s, :] = multyscale.normalization.spatial_kernel_gaussian(
+        FLODOG.bank.x, FLODOG.bank.y, FLODOG.window_sigmas[o, s, :]
+    )
 
 # Visualize each spatial avg. window
 fig, axs = plt.subplots(2, spatial_windows_LODOG.shape[1], sharex="all", sharey="all")
@@ -244,8 +249,8 @@
 # that form the denominators of the normalization
 
 # %% Energy estimates
-energies_LODOG = LODOG.normalizers_to_RMS(normalizing_coefficients_LODOG)
-energies_FLODOG = FLODOG.normalizers_to_RMS(normalizing_coefficients_FLODOG)
+energies_LODOG = LODOG.norm_energies(normalizing_coefficients_LODOG, eps=1e-6)
+energies_FLODOG = FLODOG.norm_energies(normalizing_coefficients_FLODOG, eps=1e-6)
 
 # Visualize
 vmin = min(np.min(energies_LODOG), np.min(energies_FLODOG))

docs/normalization/normalization_LODOG.py

@@ -172,7 +172,7 @@
 filters_output = ODOG.bank.apply(stimulus)
 weighted_outputs = ODOG.weight_outputs(filters_output)
 
-norm_outputs = LODOG.normalize_outputs(weighted_outputs)
+norm_outputs = LODOG.normalize_outputs(weighted_outputs, eps=1e-6)
 output_LODOG = norm_outputs.sum(axis=(0, 1))
 
 # %% Extract target prediction
@@ -261,10 +261,10 @@
 assert np.array_equal(norm_weights, LODOG.normalization_weights)
 
 # %% Normalizing images as weighted combination (tensor dot-product) of filter outputs
-normalizing_coefficients = multyscale.normalization.normalizers(weighted_outputs, norm_weights)
+normalizing_coefficients = multyscale.normalization.norm_coeffs(weighted_outputs, norm_weights)
+assert np.array_equal(normalizing_coefficients, LODOG.norm_coeffs(weighted_outputs))
 
-# Visualize each normalizing coefficient n_{o,s}, i.e.
-# the normalizer image for each individual filter f_{o,s}
+# Visualize each normalizing coefficient n_{o,s}, i.e. for each individual filter f_{o,s}
 fig, axs = plt.subplots(*normalizing_coefficients.shape[:2], sharex="all", sharey="all")
 for o, s in np.ndindex(normalizing_coefficients.shape[:2]):
     axs[o, s].imshow(normalizing_coefficients[o, s], cmap="coolwarm", extent=visextent)
@@ -346,46 +346,53 @@
 plt.show()
 
 # %% [markdown]
-# The function `multyscale.normalization.spatial_avg_windows_gaussian()`
+# The function `multyscale.normalization.spatial_kernel_gaussian()`
 # generates a ( $O\times S$ set of) Gaussian filters $G$,
-# where each Gaussian $G_{o',s'}$ is used to locally average filter $f_{o',s'}$.
+# where each Gaussian spatial averaging window
+# $G_{o',s'}$ is used to locally average filter $f_{o',s'}$.
 # In the LODOG model, all $G$ are identical.
 #
 # The parameter $\sigma$ controls the spatial size of each $G(\sigma)$ Gaussian filter;
 # thus this function takes in $O \times S$ $\sigma_{o', s'}$.
 # In the LODOG model, all $\sigma_{o', s'}$ are identical.
 
-# %% Spatial Gaussian
-sigmas = LODOG.window_sigmas
-G = multyscale.normalization.spatial_avg_windows_gaussian(LODOG.bank.x, LODOG.bank.y, sigmas)
+# %% All sigmas identical
+print(LODOG.window_sigmas)
 
-assert np.array_equal(G[0, 0, :], window)
+# %% Generate Gaussian filters
+spatial_windows = np.ndarray(filters_output.shape)
+for o, s in np.ndindex(LODOG.window_sigmas.shape[:2]):
+    spatial_windows[o, s, :] = multyscale.normalization.spatial_kernel_gaussian(
+        LODOG.bank.x, LODOG.bank.y, LODOG.window_sigmas[o, s, :]
+    )
+
+assert np.array_equal(spatial_windows[0, 0, :], window)
+assert np.array_equal(spatial_windows, LODOG.spatial_kernels())
 
 idx = (2, 3)
 
-plt.imshow(G[*idx], cmap="coolwarm", extent=visextent)
+plt.imshow(spatial_windows[*idx], cmap="coolwarm", extent=visextent)
 plt.show()
 
 # %% [markdown]
 # Applying this Gaussian window gives the _local_ (estimate of) of energy
 
-# %% Local RMS
-normalization_local_RMS = np.square(normalizing_coefficients.copy())
+# %% Local energy estimates
+normalization_local_energies = np.ndarray(normalizing_coefficients.shape)
 for o, s in np.ndindex(normalizing_coefficients.shape[:2]):
-    normalization_local_RMS[o, s] = multyscale.filters.apply(
-        window, normalization_local_RMS[o, s], padval=0
+    coeff = normalizing_coefficients[o, s] ** 2
+    energy = multyscale.filters.apply(
+        spatial_windows[o, s], coeff, padval=0
     )
+    energy = np.sqrt(energy + 1e-6) # minor offset to avoid negatives/0's
+    normalization_local_energies[o, s, :] = energy
 
-normalization_local_RMS = (
-    np.sqrt(normalization_local_RMS + 1e-6) + 1e-6
-)  # minor offset to avoid negatives/0's
-
-assert np.array_equal(normalization_local_RMS, LODOG.normalizers_to_RMS(normalizing_coefficients))
+assert np.allclose(normalization_local_energies, LODOG.norm_energies(normalizing_coefficients, eps=1e-6))
 
 # Visualize each local RMS
-fig, axs = plt.subplots(*normalization_local_RMS.shape[:2], sharex="all", sharey="all")
-for o, s in np.ndindex(normalization_local_RMS.shape[:2]):
-    axs[o, s].imshow(normalization_local_RMS[o, s], cmap="coolwarm", extent=visextent)
+fig, axs = plt.subplots(*normalization_local_energies.shape[:2], sharex="all", sharey="all")
+for o, s in np.ndindex(normalization_local_energies.shape[:2]):
+    axs[o, s].imshow(normalization_local_energies[o, s], cmap="coolwarm", extent=visextent)
 fig.supxlabel("Spatial scale/freq. $s'$")
 fig.supylabel("Orientation $o'$")
 plt.show()
@@ -420,14 +427,14 @@
 # get normalized quite differently.
 #
 # We now divide each filter(output) $f_{o',s'}$
-# by the spatial RMS of the normalizer coefficient $n_{o',s'}$:
+# by the spatial RMS of the normalizing coefficient $n_{o',s'}$:
 # $$f'_{o',s'} = \frac{f_{o',s'}}{RMS(n_{o',s'})}$$
 
 # %% Divisive normalization
 # Since the local RMSs tensor is the same (O, S, X, Y) shape as the filter outputs
 # we can simply divide
-normalized_outputs = weighted_outputs / normalization_local_RMS
-assert np.array_equal(normalized_outputs, norm_outputs)
+normalized_outputs = weighted_outputs / (normalization_local_energies + 1e-6)
+assert np.allclose(normalized_outputs, norm_outputs)
 
 # Visualize each normalized f'_{o',s'}
 fig, axs = plt.subplots(*normalized_outputs.shape[:2], sharex="all", sharey="all")
@@ -456,7 +463,7 @@
 
 # %% Recombine
 recombined_outputs = np.sum(normalized_outputs, axis=(0, 1))
-assert np.array_equal(recombined_outputs, output_LODOG)
+assert np.allclose(recombined_outputs, output_LODOG)
 
 plt.subplot(2, 2, 1)
 plt.imshow(recombined_outputs, cmap="coolwarm", extent=visextent)

docs/normalization/normalization_ODOG.py

@@ -244,7 +244,7 @@
             ]  # add this filter to normalizing coefficient
 
 # Plot each normalizing coefficient n_{o,s},
-# i.e., the normalizer image for each individual filter f_{o,s}
+# i.e., for each individual filter f_{o,s}
 fig, axs = plt.subplots(*norm_coeffs.shape[:2], sharex="all", sharey="all")
 for o, s in np.ndindex(norm_coeffs.shape[:2]):
     axs[o, s].imshow(
@@ -314,8 +314,8 @@
     normalization_weights, weighted_outputs, axes=([0, 1], [0, 1])
 )
 
-# Visualize each normalizing coefficient n_{o,s}, i.e.
-# the normalizer image for each individual filter f_{o,s}
+# Visualize each normalizing coefficient n_{o,s},
+# i.e., for each individual filter f_{o,s}
 fig, axs = plt.subplots(*normalizing_coefficients.shape[:2], sharex="all", sharey="all")
 for o, s in np.ndindex(normalizing_coefficients.shape[:2]):
     axs[o, s].imshow(normalizing_coefficients[o, s], cmap="coolwarm", extent=visextent)
@@ -404,15 +404,15 @@
 # for constructing the normalizing coefficients
 # from such weights, and the filter outputs to be normalized.
 #
-# The function `multyscale.normalization.normalizers()` creates these coefficients,
+# The function `multyscale.normalization.norm_coeffs()` creates these coefficients,
 # and note that these are identical to the $N$ constructed above.
 
 # %% Identical normalizing coefficients
-norm_coeffs = multyscale.normalization.normalizers(weighted_outputs, normalization_weights)
+norm_coeffs = multyscale.normalization.norm_coeffs(weighted_outputs, normalization_weights)
 assert np.allclose(norm_coeffs, normalizing_coefficients)
 
-# Visualize each normalizing coefficient n_{o,s}, i.e.
-# the normalizer image for each individual filter f_{o,s}
+# Visualize each normalizing coefficient n_{o,s},
+# i.e. for each individual filter f_{o,s}
 fig, axs = plt.subplots(*normalizing_coefficients.shape[:2], sharex="all", sharey="all")
 for o, s in np.ndindex(normalizing_coefficients.shape[:2]):
     axs[o, s].imshow(normalizing_coefficients[o, s], cmap="coolwarm", extent=visextent)
@@ -459,7 +459,7 @@
 #
 # Thus, the normalized filter output $f'_{o', s'}$
 # is calculating by dividing each filter(output) $f_{o',s'}$
-# by the energy of the normalizer coefficient $n_{o',s'}$:
+# by the energy of the normalizing coefficient $n_{o',s'}$:
 # $$f'_{o',s'} = \frac{f_{o',s'}}{RMS(n_{o',s'})}$$
 
 # %% Divisive normalization