Repository: kymatio/kymatio

Issue:
`compute_params_filterbank` should be simpler
The current code for `compute_params_filterbank` is messy and its documentation is outdated. It says it returns dictionaries where in fact it returns list.

It just came to my attention that we don't need dictionaries at all in this part of the code, and can make everything flow with 18 lines of code

```python
xi_max = compute_xi_max(Q)
sigma_max = compute_sigma_psi(xi_max, Q, r=r_psi)

if sigma_max <= sigma_min:
    xis = []
    sigmas = []
    elbow_xi = sigma_max
else:
    xis =  [xi_max]
    sigmas = [sigma_max]

    # High-frequency (constant-Q) region: geometric progression of xi
    while sigmas[-1] > (sigma_min * math.pow(2, 1/Q)):
        xis.append(xis[-1] / math.pow(2, 1/Q))
        sigmas.append(sigmas[-1] / math.pow(2, 1/Q))
    elbow_xi = xis[-1]

# Low-frequency (constant-bandwidth) region: arithmetic progression of xi
for q in range(1, Q):
    xis.append(elbow_xi - q/Q * elbow_xi)
    sigmas.append(sigma_min)

js = list(map(get_max_dyadic_subsampling, xis, sigmas, repeat(alpha)))
return xis, sigmas, js
```
    
instead of the 30 lines that we currently have. We can also remove the private function `_move_one_dyadic_step`, so this issue aligns with the overarching goal of simplifying `filter_bank.py`: #227, #513, #825, #833 ...    
   

Discussion:

Patch:
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
index 8105fe1..6c33240 100644
--- a/CODE_OF_CONDUCT.md
+++ b/CODE_OF_CONDUCT.md
@@ -55,7 +55,7 @@ further defined and clarified by project maintainers.
 ## Enforcement
 
 Instances of abusive, harassing, or otherwise unacceptable behavior may be
-reported by contacting the project team at michael.eickenberg@berkeley.edu. All
+reported by contacting the project team at meickenberg@flatironinstitute.org. All
 complaints will be reviewed and investigated and will result in a response that
 is deemed necessary and appropriate to the circumstances. The project team is
 obligated to maintain confidentiality with regard to the reporter of an incident.
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 3767ca6..4b8c65c 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -4,9 +4,59 @@ Contributing code
 
 This document is a 'getting started' summary for contributing code, documentation, testing, and filing issues.
 
+
+
+Feature requests/bugs: filing issues
+------------------------------------
+We use Github issues to track all bugs and feature requests. Feel free to
+open an issue if you have found a bug or want to implement feature.
+
+It is recommended to check that your issue complies with the
+following rules before submitting:
+
+-  Verify that your issue is not being currently addressed by other
+   [issues](https://github.com/kymatio/kymatio/issues?)
+   or [pull requests](https://github.com/kymatio/kymatio/pulls).
+
+-  Please ensure all code snippets and error messages are formatted in
+   appropriate code blocks.
+   See [Creating and highlighting code blocks](https://help.github.com/articles/creating-and-highlighting-code-blocks).
+
+-  Please include your operating system type and version number, as well
+   as your Python, NumPy, SciPy, PyTorch, and Kymatio versions. This information
+   can be found by running the following code snippet:
+
+  ```python
+  import platform; print(platform.platform())
+  import sys; print("Python", sys.version)
+  import numpy; print("NumPy", numpy.__version__)
+  import scipy; print("SciPy", scipy.__version__)
+  import torch; print("PyTorch", torch.__version__)
+  import kymatio; print("Kymatio", kymatio.__version__)
+  ```
+
+- If your issue is related to GPU acceleration, please copy-paste the output
+  from the [environment collection script](https://raw.githubusercontent.com/pytorch/pytorch/master/torch/utils/collect_env.py)
+  of PyTorch.
+  You can get the script and run it with:
+  ```python
+  wget https://raw.githubusercontent.com/pytorch/pytorch/master/torch/utils/collect_env.py
+  # For security purposes, please check the contents of collect_env.py before running it.
+  python collect_env.py
+  ```
+
+- If discussing a feature request, the community should collectively decide
+whether the contribution is within the scope of the project and feasible. 
+- The goal is to reach a consensus on the roadmap for the feature and its 
+implementation.
+- The community will decide on a series of pull requests (PRs) towards merging the 
+feature. In order to make the reviewing process manageable, a decision will be 
+made on how to split a large feature into multiple smaller PRs.
+- Once a consensus is reached, the implementation can begin.
+
+
 How to contribute
 -----------------
-
 The preferred way to contribute to Kymatio is to fork the
 [main repository](http://github.com/kymatio/kymatio/) on
 GitHub:
@@ -65,47 +115,6 @@ up its [documentation](http://git-scm.com/documentation).
 
 
 
-Filing bugs
------------
-We use Github issues to track all bugs and feature requests. Feel free to
-open an issue if you have found a bug or wish to see a feature implemented.
-
-It is recommended to check that your issue complies with the
-following rules before submitting:
-
--  Verify that your issue is not being currently addressed by other
-   [issues](https://github.com/kymatio/kymatio/issues?)
-   or [pull requests](https://github.com/kymatio/kymatio/pulls).
-
--  Please ensure all code snippets and error messages are formatted in
-   appropriate code blocks.
-   See [Creating and highlighting code blocks](https://help.github.com/articles/creating-and-highlighting-code-blocks).
-
--  Please include your operating system type and version number, as well
-   as your Python, NumPy, SciPy, PyTorch, and Kymatio versions. This information
-   can be found by running the following code snippet:
-
-  ```python
-  import platform; print(platform.platform())
-  import sys; print("Python", sys.version)
-  import numpy; print("NumPy", numpy.__version__)
-  import scipy; print("SciPy", scipy.__version__)
-  import torch; print("PyTorch", torch.__version__)
-  import kymatio; print("Kymatio", kymatio.__version__)
-  ```
-
-- If your issue is related to GPU acceleration, please copy-paste the output
-  from the [environment collection script](https://raw.githubusercontent.com/pytorch/pytorch/master/torch/utils/collect_env.py)
-  of PyTorch.
-  You can get the script and run it with:
-  ```python
-  wget https://raw.githubusercontent.com/pytorch/pytorch/master/torch/utils/collect_env.py
-  # For security purposes, please check the contents of collect_env.py before running it.
-  python collect_env.py
-  ```
-
-
-
 Documentation
 -------------
 
diff --git a/examples/1d/plot_filters.py b/examples/1d/plot_filters.py
index db0ced4..06dcb13 100644
--- a/examples/1d/plot_filters.py
+++ b/examples/1d/plot_filters.py
@@ -69,8 +69,7 @@ def scattering1D_filter_factory():
   # first-order wavelet filters (`psi1_f`), and the second-order filters
   # (`psi2_f`).
   
-  phi_f, psi1_f, psi2_f, _ = scattering_filter_factory(np.log2(N), J, Q, T)
-  
+  phi_f, psi1_f, psi2_f = scattering_filter_factory(np.log2(N), J, Q, T)
   ###############################################################################
   # The `phi_f` output is a dictionary where each integer key corresponds points
   # to the instantiation of the filter at a certain resolution. In other words,
@@ -151,8 +150,8 @@ def scattering1D_filter_factory_T():
               default_str = ' (default)'
           else: 
               default_str = ''
-          phi_f, psi1_f, psi2_f, _ = scattering_filter_factory(np.log2(N), J, Q, T)
           msg = 'J=' + str(J) + ', T=' +str(T) + default_str + ': LP-filter width $\sigma_{low}$=' + str(phi_f['sigma']) + ' (=0.1/T)'                
+          phi_f, psi1_f, psi2_f = scattering_filter_factory(np.log2(N), J, Q, T)
           print(msg)
           assert(phi_f['sigma']==0.1/T)
 
diff --git a/kymatio/scattering1d/filter_bank.py b/kymatio/scattering1d/filter_bank.py
index cc5ed86..7f3bfe0 100644
--- a/kymatio/scattering1d/filter_bank.py
+++ b/kymatio/scattering1d/filter_bank.py
@@ -347,50 +347,6 @@ def get_max_dyadic_subsampling(xi, sigma, alpha):
     return j
 
 
-def _move_one_dyadic_step(cv, Q, alpha):
-    """
-    Computes the parameters of the next wavelet on the low frequency side,
-    based on the parameters of the current wavelet.
-
-    This function is used in the loop defining all the filters, starting
-    at the wavelet frequency and then going to the low frequencies by
-    dyadic steps. This makes the loop in compute_params_filterbank much
-    simpler to read.
-
-    The steps are defined as:
-    xi_{n+1} = 2^{-1/Q} xi_n
-    sigma_{n+1} = 2^{-1/Q} sigma_n
-
-    Parameters
-    ----------
-    cv : dictionary
-        stands for current_value. Is a dictionary with keys:
-        *'key': a tuple (j, n) where n is a counter and j is the maximal
-            dyadic subsampling accepted by this wavelet.
-        *'xi': central frequency of the wavelet
-        *'sigma': width of the wavelet
-    Q : int
-        number of wavelets per octave. Controls the relationship between
-        the frequency and width of the current wavelet and the next wavelet.
-    alpha : float, optional
-        tolerance parameter for the aliasing. The larger alpha,
-        the more conservative the algorithm is.
-
-    Returns
-    -------
-    new_cv : dictionary
-        a dictionary with the same keys as the ones listed for cv,
-        whose values are updated
-    """
-    factor = 1. / math.pow(2., 1. / Q)
-    n = cv['key']
-    new_cv = {'xi': cv['xi'] * factor, 'sigma': cv['sigma'] * factor}
-    # compute the new j
-    new_cv['j'] = get_max_dyadic_subsampling(new_cv['xi'], new_cv['sigma'], alpha=alpha)
-    new_cv['key'] = n + 1
-    return new_cv
-
-
 def compute_xi_max(Q):
     """
     Computes the maximal xi to use for the Morlet family, depending on Q.
@@ -418,11 +374,6 @@ def compute_params_filterbank(sigma_min, Q, alpha, r_psi=math.sqrt(0.5)):
     are translated. sigma_min limits the smallest frequential width
     among all filters, while preserving the coverage of the whole frequency
     axis.
-
-    The keys of the dictionaries are tuples of integers (j, n) where n is a
-    counter (starting at 0 for the highest frequency filter) and j is the
-    maximal dyadic subsampling accepted by this filter.
-
     Parameters
     ----------
     sigma_min : float
@@ -439,52 +390,42 @@ def compute_params_filterbank(sigma_min, Q, alpha, r_psi=math.sqrt(0.5)):
         Should be >0 and <1. Controls the redundancy of the filters
         (the larger r_psi, the larger the overlap between adjacent wavelets).
         Defaults to sqrt(0.5).
-
     Returns
     -------
-    xi : dictionary
-        dictionary containing the central frequencies of the wavelets.
-    sigma : dictionary
-        dictionary containing the frequential widths of the wavelets.
-
-    Refs
-    ----
-    Convolutional operators in the time-frequency domain, 2.1.3, V. Lostanlen,
-    PhD Thesis, 2017
-    https://tel.archives-ouvertes.fr/tel-01559667
+    xis : list
+        central frequencies of the filters.
+    sigmas : list
+        bandwidths of the filters.
+    js : list
+        maximal dyadic subsampling accepted by the filters.
+        j=0 stands for no subsampling, j=1 stands for half subsampling, etc.
     """
     xi_max = compute_xi_max(Q)
     sigma_max = compute_sigma_psi(xi_max, Q, r=r_psi)
 
-    xi = []
-    sigma = []
-    j = []
-
     if sigma_max <= sigma_min:
-        # in this exceptional case, we will not go through the loop, so
-        # we directly assign
-        last_xi = sigma_max
+        xis = []
+        sigmas = []
+        elbow_xi = sigma_max
     else:
-        # fill all the dyadic wavelets as long as possible
-        current = {'key': 0, 'j': 0, 'xi': xi_max, 'sigma': sigma_max}
-        while current['sigma'] > sigma_min:  # while we can attribute something
-            xi.append(current['xi'])
-            sigma.append(current['sigma'])
-            j.append(current['j'])
-            current = _move_one_dyadic_step(current, Q, alpha=alpha)
-        # get the last key
-        last_xi = xi[-1]
-    # fill num_interm wavelets between last_xi and 0, both excluded
-    num_intermediate = Q - 1
-    for q in range(1, num_intermediate + 1):
-        factor = (num_intermediate + 1. - q) / (num_intermediate + 1.)
-        new_xi = factor * last_xi
-        new_sigma = sigma_min
-        xi.append(new_xi)
-        sigma.append(new_sigma)
-        j.append(get_max_dyadic_subsampling(new_xi, new_sigma, alpha=alpha))
-    # return results
-    return xi, sigma, j
+        xis =  [xi_max]
+        sigmas = [sigma_max]
+
+        # High-frequency (constant-Q) region: geometric progression of xi
+        while sigmas[-1] > (sigma_min * math.pow(2, 1/Q)):
+            xis.append(xis[-1] / math.pow(2, 1/Q))
+            sigmas.append(sigmas[-1] / math.pow(2, 1/Q))
+        elbow_xi = xis[-1]
+
+    # Low-frequency (constant-bandwidth) region: arithmetic progression of xi
+    for q in range(1, Q):
+        xis.append(elbow_xi - q/Q * elbow_xi)
+        sigmas.append(sigma_min)
+
+    js = [
+        get_max_dyadic_subsampling(xi, sigma, alpha) for xi, sigma in zip(xis, sigmas)
+    ]
+    return xis, sigmas, js
 
 
 def calibrate_scattering_filters(J, Q, T, alpha, r_psi=math.sqrt(0.5), sigma0=0.1):
@@ -556,7 +497,7 @@ def calibrate_scattering_filters(J, Q, T, alpha, r_psi=math.sqrt(0.5), sigma0=0.
 
 
 def scattering_filter_factory(J_support, J_scattering, Q, T, r_psi=math.sqrt(0.5),
-                              criterion_amplitude=1e-3, normalize='l1',
+                              normalize='l1',
                               max_subsampling=None, sigma0=0.1, alpha=5.,
                               P_max=5, eps=1e-7, **kwargs):
     """
@@ -589,9 +530,6 @@ def scattering_filter_factory(J_support, J_scattering, Q, T, r_psi=math.sqrt(0.5
         Should be >0 and <1. Controls the redundancy of the filters
         (the larger r_psi, the larger the overlap between adjacent wavelets).
         Defaults to sqrt(0.5).
-    criterion_amplitude : float, optional
-        Represents the numerical error which is allowed to be lost after
-        convolution and padding. Defaults to 1e-3.
     normalize : string, optional
         Normalization convention for the filters (in the
         temporal domain). Supported values include 'l1' and 'l2'; a ValueError
@@ -643,10 +581,6 @@ def scattering_filter_factory(J_support, J_scattering, Q, T, r_psi=math.sqrt(0.5
         The keys of this dictionary are of th etype (j, n) where n is an
         integer counting the filters and j is the maximal dyadic subsampling
         which can be performed on top of this filter without aliasing.
-    t_max_phi : int
-        temporal size to use to pad the signal on the right and on the
-        left by making at most criterion_amplitude error. Assumes that the
-        temporal support of the low-pass filter is larger than all filters.
 
     Refs
     ----
@@ -733,10 +667,5 @@ def scattering_filter_factory(J_support, J_scattering, Q, T, r_psi=math.sqrt(0.5
     phi_f['sigma'] = sigma_low
     phi_f['j'] = log2_T
 
-    # compute the support size allowing to pad without boundary errors
-    # at the finest resolution  # TODO
-    t_max_phi = compute_temporal_support(
-        phi_f[0].reshape(1, -1), criterion_amplitude=criterion_amplitude)
-
     # return results
-    return phi_f, psi1_f, psi2_f, t_max_phi
+    return phi_f, psi1_f, psi2_f
diff --git a/kymatio/scattering1d/frontend/base_frontend.py b/kymatio/scattering1d/frontend/base_frontend.py
index aa7bc4a..cba2cce 100644
--- a/kymatio/scattering1d/frontend/base_frontend.py
+++ b/kymatio/scattering1d/frontend/base_frontend.py
@@ -4,8 +4,8 @@ import numbers
 
 import numpy as np
 
-from ..filter_bank import scattering_filter_factory
-from ..utils import (compute_border_indices, compute_padding, compute_minimum_support_to_pad,
+from ..filter_bank import compute_temporal_support, gauss_1d, scattering_filter_factory
+from ..utils import (compute_border_indices, compute_padding,
 compute_meta_scattering, precompute_size_scattering)
 
 
@@ -38,7 +38,6 @@ class ScatteringBase1D(ScatteringBase):
         self.alpha = 5.
         self.P_max = 5
         self.eps = 1e-7
-        self.criterion_amplitude = 1e-3
         self.normalize = 'l1'
 
         # check the shape
@@ -62,11 +61,11 @@ class ScatteringBase1D(ScatteringBase):
         self.log2_T = math.floor(math.log2(self.T))
 
         # Compute the minimum support to pad (ideally)
-        min_to_pad = compute_minimum_support_to_pad(
-            self.N, self.J, self.Q, self.T, r_psi=self.r_psi, sigma0=self.sigma0,
-            alpha=self.alpha, P_max=self.P_max, eps=self.eps,
-            criterion_amplitude=self.criterion_amplitude,
-            normalize=self.normalize)
+        phi_f = gauss_1d(
+            self.N, self.sigma0/self.T, self.normalize, self.P_max, self.eps)
+        min_to_pad = 3 * compute_temporal_support(
+            phi_f.reshape(1, -1), criterion_amplitude=1e-3)
+
         # to avoid padding more than N - 1 on the left and on the right,
         # since otherwise torch sends nans
         J_max_support = int(np.floor(np.log2(3 * self.N - 2)))
@@ -80,9 +79,8 @@ class ScatteringBase1D(ScatteringBase):
 
     def create_filters(self):
         # Create the filters
-        self.phi_f, self.psi1_f, self.psi2_f, _ = scattering_filter_factory(
+        self.phi_f, self.psi1_f, self.psi2_f = scattering_filter_factory(
             self.J_pad, self.J, self.Q, self.T, normalize=self.normalize,
-            criterion_amplitude=self.criterion_amplitude,
             r_psi=self.r_psi, sigma0=self.sigma0, alpha=self.alpha,
             P_max=self.P_max, eps=self.eps)
 
diff --git a/kymatio/scattering1d/utils.py b/kymatio/scattering1d/utils.py
index 4d5516f..55c44de 100644
--- a/kymatio/scattering1d/utils.py
+++ b/kymatio/scattering1d/utils.py
@@ -1,6 +1,6 @@
 import numpy as np
 import math
-from .filter_bank import scattering_filter_factory, calibrate_scattering_filters
+from .filter_bank import calibrate_scattering_filters
 
 def compute_border_indices(log2_T, J, i0, i1):
     """
@@ -69,74 +69,6 @@ def compute_padding(J_pad, N):
         raise ValueError('Too large padding value, will lead to NaN errors')
     return pad_left, pad_right
 
-def compute_minimum_support_to_pad(N, J, Q, T, criterion_amplitude=1e-3,
-                                       normalize='l1', r_psi=math.sqrt(0.5),
-                                       sigma0=1e-1, alpha=5., P_max=5, eps=1e-7):
-
-
-    """
-    Computes the support to pad given the input size and the parameters of the
-    scattering transform.
-
-    Parameters
-    ----------
-    N : int
-        temporal size of the input signal
-    J : int
-        scale of the scattering
-    Q : int
-        number of wavelets per octave
-    T : int
-        temporal support of low-pass filter, controlling amount of imposed
-        time-shift invariance and maximum subsampling
-    normalize : string, optional
-        normalization type for the wavelets.
-        Only `'l2'` or `'l1'` normalizations are supported.
-        Defaults to `'l1'`
-    criterion_amplitude: float `>0` and `<1`, optional
-        Represents the numerical error which is allowed to be lost after
-        convolution and padding.
-        The larger criterion_amplitude, the smaller the padding size is.
-        Defaults to `1e-3`
-    r_psi : float, optional
-        Should be `>0` and `<1`. Controls the redundancy of the filters
-        (the larger r_psi, the larger the overlap between adjacent
-        wavelets).
-        Defaults to `sqrt(0.5)`.
-    sigma0 : float, optional
-        parameter controlling the frequential width of the
-        low-pass filter at J_scattering=0; at a an absolute J_scattering,
-        it is equal to :math:`\\frac{\\sigma_0}{2^J}`.
-        Defaults to `1e-1`.
-    alpha : float, optional
-        tolerance factor for the aliasing after subsampling.
-        The larger the alpha, the more conservative the value of maximal
-        subsampling is.
-        Defaults to `5`.
-    P_max : int, optional
-        maximal number of periods to use to make sure that the Fourier
-        transform of the filters is periodic.
-        `P_max = 5` is more than enough for double precision.
-        Defaults to `5`.
-    eps : float, optional
-        required machine precision for the periodization (single
-        floating point is enough for deep learning applications).
-        Defaults to `1e-7`.
-
-    Returns
-    -------
-    min_to_pad: int
-        minimal value to pad the signal on one size to avoid any
-        boundary error.
-    """
-    J_tentative = int(np.ceil(np.log2(N)))
-    _, _, _, t_max_phi = scattering_filter_factory(
-        J_tentative, J, Q, T, normalize=normalize, to_torch=False,
-        max_subsampling=0, criterion_amplitude=criterion_amplitude,
-        r_psi=r_psi, sigma0=sigma0, alpha=alpha, P_max=P_max, eps=eps)
-    min_to_pad = 3 * t_max_phi
-    return min_to_pad
-
 
 def precompute_size_scattering(J, Q, T, max_order=2, detail=False):
     """Get size of the scattering transform


Test patch:
diff --git a/tests/scattering1d/test_numpy_scattering1d.py b/tests/scattering1d/test_numpy_scattering1d.py
index b375ed4..c9da21d 100644
--- a/tests/scattering1d/test_numpy_scattering1d.py
+++ b/tests/scattering1d/test_numpy_scattering1d.py
@@ -96,5 +96,5 @@ class TestScattering1DNumpy:
                     default_str = ' (default)'
                 else:
                     default_str = ''
-                phi_f, psi1_f, psi2_f, _ = scattering_filter_factory(np.log2(N), J, Q, T)
+                phi_f, psi1_f, psi2_f = scattering_filter_factory(np.log2(N), J, Q, T)
                 assert(phi_f['sigma']==0.1/T)