    def compute_forecast(self):
        """
        Generate a blended nowcast ensemble by using the Short-Term Ensemble
        Prediction System (STEPS) method.

        Parameters
        ----------
        precip: array-like
          Array of shape (ar_order+1,m,n) containing the input precipitation fields
          ordered by timestamp from oldest to newest. The time steps between the
          inputs are assumed to be regular.
        precip_models: array-like
          Either raw (NWP) model forecast data or decomposed (NWP) model forecast data.
          If you supply decomposed data, it needs to be an array of shape
          (n_models,timesteps+1) containing, per timestep (t=0 to lead time here) and
          per (NWP) model or model ensemble member, a dictionary with a list of cascades
          obtained by calling a method implemented in :py:mod:`pysteps.cascade.decomposition`.
          If you supply the original (NWP) model forecast data, it needs to be an array of shape
          (n_models,timestep+1,m,n) containing precipitation (or other) fields, which will
          then be decomposed in this function.

          Depending on your use case it can be advantageous to decompose the model
          forecasts outside beforehand, as this slightly reduces calculation times.
          This is possible with :py:func:`pysteps.blending.utils.decompose_NWP`,
          :py:func:`pysteps.blending.utils.compute_store_nwp_motion`, and
          :py:func:`pysteps.blending.utils.load_NWP`. However, if you have a lot of (NWP) model
          members (e.g. 1 model member per nowcast member), this can lead to excessive memory
          usage.

          To further reduce memory usage, both this array and the ``velocity_models`` array
          can be given as float32. They will then be converted to float64 before computations
          to minimize loss in precision.

          In case of one (deterministic) model as input, add an extra dimension to make sure
          precip_models is four dimensional prior to calling this function.
        velocity: array-like
          Array of shape (2,m,n) containing the x- and y-components of the advection
          field. The velocities are assumed to represent one time step between the
          inputs. All values are required to be finite.
        velocity_models: array-like
          Array of shape (n_models,timestep,2,m,n) containing the x- and y-components
          of the advection field for the (NWP) model field per forecast lead time.
          All values are required to be finite.

          To reduce memory usage, this array
          can be given as float32. They will then be converted to float64 before computations
          to minimize loss in precision.
        time_steps: int or list of floats
          Number of time steps to forecast or a list of time steps for which the
          forecasts are computed (relative to the input time step). The elements of
          the list are required to be in ascending order.
        issue_time: datetime
          is issued.
        config: StepsBlendingConfig
            Provides a set of configuration parameters for the nowcast ensemble generation.

        Returns
        -------
        out: ndarray
          If return_output is True, a four-dimensional array of shape
          (n_ens_members,num_timesteps,m,n) containing a time series of forecast
          precipitation fields for each ensemble member. Otherwise, a None value
          is returned. The time series starts from t0+timestep, where timestep is
          taken from the input precipitation fields precip. If measure_time is True, the
          return value is a three-element tuple containing the nowcast array, the
          initialization time of the nowcast generator and the time used in the
          main loop (seconds).

        See also
        --------
        :py:mod:`pysteps.extrapolation.interface`, :py:mod:`pysteps.cascade.interface`,
        :py:mod:`pysteps.noise.interface`, :py:func:`pysteps.noise.utils.compute_noise_stddev_adjs`

        References
        ----------
        :cite:`Seed2003`, :cite:`BPS2004`, :cite:`BPS2006`, :cite:`SPN2013`, :cite:`PCH2019b`

        Notes
        -----
        1. The blending currently does not blend the beta-parameters in the parametric
        noise method. It is recommended to use the non-parameteric noise method.

        2. If blend_nwp_members is True, the BPS2006 method for the weights is
        suboptimal. It is recommended to use the SPN2013 method instead.

        3. Not yet implemented (and neither in the steps nowcasting module): The regression
        of the lag-1 and lag-2 parameters to their climatological values. See also eq.
        12 - 19 in :cite: `BPS2004`. By doing so, the Phi parameters change over time,
        which enhances the AR process. This can become a future development if this
        turns out to be a warranted functionality.
        """

        self.__check_inputs()
        self.__print_forecast_info()
        # Measure time for initialization
        if self.__config.measure_time:
            self.__start_time_init = time.time()

        # Slice the precipitation field to only use the last ar_order + 1 fields
        self.__precip = self.__precip[-(self.__config.ar_order + 1) :, :, :].copy()
        self.__initialize_nowcast_components()
        self.__prepare_radar_and_NWP_fields()

        # Determine if rain is present in both radar and NWP fields
        if self.__params.zero_precip_radar and self.__params.zero_precip_model_fields:
            return self.__zero_precipitation_forecast()
        else:
            # Prepare the data for the zero precipitation radar case and initialize the noise correctly
            if self.__params.zero_precip_radar:
                self.__prepare_nowcast_for_zero_radar()
            else:
                self.__state.precip_noise_input = self.__precip.copy()
            self.__initialize_noise()
            self.__estimate_ar_parameters_radar()
            self.__multiply_precip_cascade_to_match_ensemble_members()
            self.__initialize_random_generators()
            self.__prepare_forecast_loop()
            self.__initialize_noise_cascades()
            if self.__config.measure_time:
                self.__init_time = self.__measure_time(
                    "initialization", self.__start_time_init
                )

            self.__blended_nowcast_main_loop()
            # Stack and return the forecast output
            if self.__config.return_output:
                self.__state.final_blended_forecast = np.stack(
                    [
                        np.stack(self.__state.final_blended_forecast[j])
                        for j in range(self.__config.n_ens_members)
                    ]
                )
                if self.__config.measure_time:
                    return (
                        self.__state.final_blended_forecast,
                        self.__init_time,
                        self.__mainloop_time,
                    )
                else:
                    return self.__state.final_blended_forecast
            else:
                return None
Generate a blended nowcast ensemble by using the Short-Term Ensemble
Prediction System (STEPS) method.

Parameters
----------
precip: array-like
Array of shape (ar_order+1,m,n) containing the input precipitation fields
ordered by timestamp from oldest to newest. The time steps between the
inputs are assumed to be regular.
precip_models: array-like
Either raw (NWP) model forecast data or decomposed (NWP) model forecast data.
If you supply decomposed data, it needs to be an array of shape
(n_models,timesteps+1) containing, per timestep (t=0 to lead time here) and
per (NWP) model or model ensemble member, a dictionary with a list of cascades
obtained by calling a method implemented in :py:mod:`pysteps.cascade.decomposition`.
If you supply the original (NWP) model forecast data, it needs to be an array of shape
(n_models,timestep+1,m,n) containing precipitation (or other) fields, which will
then be decomposed in this function.

Depending on your use case it can be advantageous to decompose the model
forecasts outside beforehand, as this slightly reduces calculation times.
This is possible with :py:func:`pysteps.blending.utils.decompose_NWP`,
:py:func:`pysteps.blending.utils.compute_store_nwp_motion`, and
:py:func:`pysteps.blending.utils.load_NWP`. However, if you have a lot of (NWP) model
members (e.g. 1 model member per nowcast member), this can lead to excessive memory
usage.

To further reduce memory usage, both this array and the ``velocity_models`` array
can be given as float32. They will then be converted to float64 before computations
to minimize loss in precision.

In case of one (deterministic) model as input, add an extra dimension to make sure
precip_models is four dimensional prior to calling this function.
velocity: array-like
Array of shape (2,m,n) containing the x- and y-components of the advection
field. The velocities are assumed to represent one time step between the
inputs. All values are required to be finite.
velocity_models: array-like
Array of shape (n_models,timestep,2,m,n) containing the x- and y-components
of the advection field for the (NWP) model field per forecast lead time.
All values are required to be finite.

To reduce memory usage, this array
can be given as float32. They will then be converted to float64 before computations
to minimize loss in precision.
time_steps: int or list of floats
Number of time steps to forecast or a list of time steps for which the
forecasts are computed (relative to the input time step). The elements of
the list are required to be in ascending order.
issue_time: datetime
is issued.
config: StepsBlendingConfig
Provides a set of configuration parameters for the nowcast ensemble generation.

Returns
-------
out: ndarray
If return_output is True, a four-dimensional array of shape
(n_ens_members,num_timesteps,m,n) containing a time series of forecast
precipitation fields for each ensemble member. Otherwise, a None value
is returned. The time series starts from t0+timestep, where timestep is
taken from the input precipitation fields precip. If measure_time is True, the
return value is a three-element tuple containing the nowcast array, the
initialization time of the nowcast generator and the time used in the
main loop (seconds).

See also
--------
:py:mod:`pysteps.extrapolation.interface`, :py:mod:`pysteps.cascade.interface`,
:py:mod:`pysteps.noise.interface`, :py:func:`pysteps.noise.utils.compute_noise_stddev_adjs`

References
----------
:cite:`Seed2003`, :cite:`BPS2004`, :cite:`BPS2006`, :cite:`SPN2013`, :cite:`PCH2019b`

Notes
-----
1. The blending currently does not blend the beta-parameters in the parametric
noise method. It is recommended to use the non-parameteric noise method.

2. If blend_nwp_members is True, the BPS2006 method for the weights is
suboptimal. It is recommended to use the SPN2013 method instead.

3. Not yet implemented (and neither in the steps nowcasting module): The regression
of the lag-1 and lag-2 parameters to their climatological values. See also eq.
12 - 19 in :cite: `BPS2004`. By doing so, the Phi parameters change over time,
which enhances the AR process. This can become a future development if this
turns out to be a warranted functionality.
pysteps/blending/steps.py