def _estimate_perturbation_params(
    forecast_err,
    forecast_gen,
    errdist_window_radius,
    acf_window_radius,
    interp_window_radius,
    measure_time,
    num_workers,
    use_multiprocessing,
):
    """
    Estimate perturbation generator parameters from forecast errors."""
    pert_gen = {}
    pert_gen["m"] = forecast_err.shape[0]
    pert_gen["n"] = forecast_err.shape[1]

    feature_coords = forecast_gen["feature_coords"]

    print("Estimating perturbation parameters... ", end="", flush=True)

    if measure_time:
        starttime = time.time()

    mask_finite = np.isfinite(forecast_err)

    forecast_err = forecast_err.copy()
    forecast_err[~mask_finite] = 1.0

    weights_dist = _compute_window_weights(
        feature_coords,
        forecast_err.shape[0],
        forecast_err.shape[1],
        errdist_window_radius,
    )

    acf_winfunc = _window_tukey if feature_coords.shape[0] > 1 else _window_uniform

    def worker(i):
        weights_acf = acf_winfunc(
            forecast_err.shape[0],
            forecast_err.shape[1],
            feature_coords[i, 0],
            feature_coords[i, 1],
            acf_window_radius,
            acf_window_radius,
        )

        mask = np.logical_and(mask_finite, weights_dist[i] > 0.1)
        if np.sum(mask) > 10 and np.sum(np.abs(forecast_err[mask] - 1.0) >= 1e-3) > 10:
            distpar = _fit_dist(forecast_err, stats.lognorm, weights_dist[i], mask)
            inv_acf_mapping = _compute_inverse_acf_mapping(stats.lognorm, distpar)
            mask_acf = weights_acf > 1e-4
            std = _weighted_std(forecast_err[mask_acf], weights_dist[i][mask_acf])
            if np.isfinite(std):
                acf = inv_acf_mapping(
                    _compute_sample_acf(weights_acf * (forecast_err - 1.0) / std)
                )
                acf = _fit_acf(acf)
            else:
                distpar = None
                std = None
                acf = None
        else:
            distpar = None
            std = None
            acf = None

        return distpar, std, np.sqrt(np.abs(np.fft.rfft2(acf)))

    dist_params = []
    stds = []
    acf_fft_ampl = []

    if DASK_IMPORTED and num_workers > 1:
        res = []
        for i in range(feature_coords.shape[0]):
            res.append(dask.delayed(worker)(i))
        scheduler = "threads" if not use_multiprocessing else "multiprocessing"
        res = dask.compute(*res, num_workers=num_workers, scheduler=scheduler)
        for r in res:
            dist_params.append(r[0])
            stds.append(r[1])
            acf_fft_ampl.append(r[2])
    else:
        for i in range(feature_coords.shape[0]):
            r = worker(i)
            dist_params.append(r[0])
            stds.append(r[1])
            acf_fft_ampl.append(r[2])

    pert_gen["dist_param"] = dist_params
    pert_gen["std"] = stds
    pert_gen["acf_fft_ampl"] = acf_fft_ampl

    weights = _compute_window_weights(
        feature_coords,
        forecast_err.shape[0],
        forecast_err.shape[1],
        interp_window_radius,
    )
    pert_gen["weights"] = weights / np.sum(weights, axis=0)

    if measure_time:
        print(f"{time.time() - starttime:.2f} seconds.")
    else:
        print("done.")

    return pert_gen
Estimate perturbation generator parameters from forecast errors.
pysteps/nowcasts/linda.py