    def __blend_cascades(self, t_sub, j, worker_state):
        """
        Blend extrapolated, NWP model, and noise cascades using predefined weights.
        Computes both full and model-only blends while also blending means and
        standard deviations across scales.
        """
        worker_state.subtimestep_index = np.where(
            np.array(worker_state.subtimesteps) == t_sub
        )[0][0]
        # First concatenate the cascades and the means and sigmas
        # precip_models = [n_models,timesteps,n_cascade_levels,m,n]
        if self.__config.blend_nwp_members:
            cascade_stack_all_components = np.concatenate(
                (
                    worker_state.precip_extrapolated_decomp[
                        None, worker_state.subtimestep_index
                    ],
                    worker_state.precip_models_cascades_timestep,
                    worker_state.noise_extrapolated_decomp[
                        None, worker_state.subtimestep_index
                    ],
                ),
                axis=0,
            )  # [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
            means_stacked = np.concatenate(
                (
                    worker_state.mean_extrapolation[None, :],
                    worker_state.mean_models_timestep,
                ),
                axis=0,
            )
            sigmas_stacked = np.concatenate(
                (
                    worker_state.std_extrapolation[None, :],
                    worker_state.std_models_timestep,
                ),
                axis=0,
            )
        else:
            cascade_stack_all_components = np.concatenate(
                (
                    worker_state.precip_extrapolated_decomp[
                        None, worker_state.subtimestep_index
                    ],
                    worker_state.precip_models_cascades_timestep[None, j],
                    worker_state.noise_extrapolated_decomp[
                        None, worker_state.subtimestep_index
                    ],
                ),
                axis=0,
            )  # [(extr_field, n_model_fields, noise), n_cascade_levels, ...]
            means_stacked = np.concatenate(
                (
                    worker_state.mean_extrapolation[None, :],
                    worker_state.mean_models_timestep[None, j],
                ),
                axis=0,
            )
            sigmas_stacked = np.concatenate(
                (
                    worker_state.std_extrapolation[None, :],
                    worker_state.std_models_timestep[None, j],
                ),
                axis=0,
            )

        # First determine the blending weights if method is spn. The
        # weights for method bps have already been determined.

        if self.__config.weights_method == "spn":
            worker_state.weights = np.zeros(
                (
                    cascade_stack_all_components.shape[0],
                    self.__config.n_cascade_levels,
                )
            )
            for i in range(self.__config.n_cascade_levels):
                # Determine the normalized covariance matrix (containing)
                # the cross-correlations between the models
                cascade_stack_all_components_temp = np.stack(
                    [
                        cascade_stack_all_components[n_model, i, :, :].flatten()
                        for n_model in range(cascade_stack_all_components.shape[0] - 1)
                    ]
                )  # -1 to exclude the noise component
                covariance_nwp_models = np.ma.corrcoef(
                    np.ma.masked_invalid(cascade_stack_all_components_temp)
                )
                # Determine the weights for this cascade level
                worker_state.weights[:, i] = calculate_weights_spn(
                    correlations=worker_state.rho_final_blended_forecast[:, i],
                    covariance=covariance_nwp_models,
                )

        # Blend the extrapolation, (NWP) model(s) and noise cascades
        worker_state.final_blended_forecast_cascades = blending.utils.blend_cascades(
            cascades_norm=cascade_stack_all_components, weights=worker_state.weights
        )

        # Also blend the cascade without the extrapolation component
        worker_state.final_blended_forecast_cascades_mod_only = (
            blending.utils.blend_cascades(
                cascades_norm=cascade_stack_all_components[1:, :],
                weights=worker_state.weights_model_only,
            )
        )

        # Blend the means and standard deviations
        # Input is array of shape [number_components, scale_level, ...]
        (
            worker_state.final_blended_forecast_means,
            worker_state.final_blended_forecast_stds,
        ) = blend_means_sigmas(
            means=means_stacked, sigmas=sigmas_stacked, weights=worker_state.weights
        )
        # Also blend the means and sigmas for the cascade without extrapolation

        (
            worker_state.final_blended_forecast_means_mod_only,
            worker_state.final_blended_forecast_stds_mod_only,
        ) = blend_means_sigmas(
            means=means_stacked[1:, :],
            sigmas=sigmas_stacked[1:, :],
            weights=worker_state.weights_model_only,
        )
Blend extrapolated, NWP model, and noise cascades using predefined weights.
Computes both full and model-only blends while also blending means and
standard deviations across scales.
pysteps/blending/steps.py