def _linda_deterministic_init(
    precip,
    velocity,
    feature_method,
    max_num_features,
    feature_kwargs,
    ari_order,
    kernel_type,
    localization_window_radius,
    extrap_method,
    extrap_kwargs,
    add_perturbations,
    num_workers,
    measure_time,
):
    """Initialize the deterministic LINDA nowcast model."""
    forecast_gen = {}
    forecast_gen["velocity"] = velocity
    forecast_gen["extrap_method"] = extrap_method
    forecast_gen["ari_order"] = ari_order
    forecast_gen["add_perturbations"] = add_perturbations
    forecast_gen["num_workers"] = num_workers
    forecast_gen["measure_time"] = measure_time

    precip = precip[-(ari_order + 2) :]
    input_length = precip.shape[0]

    starttime_init = time.time()

    extrapolator = extrapolation.get_method(extrap_method)
    extrap_kwargs = extrap_kwargs.copy()
    extrap_kwargs["allow_nonfinite_values"] = True
    forecast_gen["extrapolator"] = extrapolator
    forecast_gen["extrap_kwargs"] = extrap_kwargs

    # detect features from the most recent input field
    if feature_method in {"blob", "shitomasi"}:
        precip_ = precip[-1].copy()
        precip_[~np.isfinite(precip_)] = 0.0
        feature_detector = feature.get_method(feature_method)

        if measure_time:
            starttime = time.time()

        feature_kwargs = feature_kwargs.copy()
        feature_kwargs["max_num_features"] = max_num_features

        feature_coords = np.fliplr(feature_detector(precip_, **feature_kwargs)[:, :2])

        feature_type = "blobs" if feature_method == "blob" else "corners"
        print("")
        print("Detecting features... ", end="", flush=True)
        if measure_time:
            print(
                f"found {feature_coords.shape[0]} {feature_type} in {time.time() - starttime:.2f} seconds."
            )
        else:
            print(f"found {feature_coords.shape[0]} {feature_type}.")

        if len(feature_coords) == 0:
            raise ValueError(
                "no features found, check input data and feature detector configuration"
            )
    elif feature_method == "domain":
        feature_coords = np.zeros((1, 2), dtype=int)
    else:
        raise NotImplementedError(
            "feature detector '%s' not implemented" % feature_method
        )
    forecast_gen["feature_coords"] = feature_coords

    # compute interpolation weights
    interp_weights = _compute_window_weights(
        feature_coords,
        precip.shape[1],
        precip.shape[2],
        localization_window_radius,
    )
    interp_weights /= np.sum(interp_weights, axis=0)
    forecast_gen["interp_weights"] = interp_weights

    # transform the input fields to the Lagrangian coordinates
    precip_lagr = np.empty(precip.shape)

    def worker(i):
        precip_lagr[i, :] = extrapolator(
            precip[i, :],
            velocity,
            input_length - 1 - i,
            "min",
            **extrap_kwargs,
        )[-1]

    if DASK_IMPORTED and num_workers > 1:
        res = []

    print("Transforming to Lagrangian coordinates... ", end="", flush=True)

    if measure_time:
        starttime = time.time()

    for i in range(precip.shape[0] - 1):
        if DASK_IMPORTED and num_workers > 1:
            res.append(dask.delayed(worker)(i))
        else:
            worker(i)

    if DASK_IMPORTED and num_workers > 1:
        dask.compute(*res, num_workers=min(num_workers, len(res)), scheduler="threads")
    precip_lagr[-1] = precip[-1]

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

    # compute advection mask and set nan to pixels, where one or more of the
    # advected input fields has a nan value
    mask_adv = np.all(np.isfinite(precip_lagr), axis=0)
    forecast_gen["mask_adv"] = mask_adv
    for i in range(precip_lagr.shape[0]):
        precip_lagr[i, ~mask_adv] = np.nan

    # compute differenced input fields in the Lagrangian coordinates
    precip_lagr_diff = np.diff(precip_lagr, axis=0)

    # estimate parameters of the deterministic model (i.e. the convolution and
    # the ARI process)

    print("Estimating the first convolution kernel... ", end="", flush=True)

    if measure_time:
        starttime = time.time()

    # estimate convolution kernel for the differenced component
    convol_weights = _compute_window_weights(
        feature_coords,
        precip.shape[1],
        precip.shape[2],
        localization_window_radius,
    )

    kernels_1 = _estimate_convol_params(
        precip_lagr_diff[-2],
        precip_lagr_diff[-1],
        convol_weights,
        mask_adv,
        kernel_type=kernel_type,
        num_workers=num_workers,
    )
    forecast_gen["kernels_1"] = kernels_1

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

    # compute convolved difference fields
    precip_lagr_diff_c = precip_lagr_diff[:-1].copy()
    for i in range(precip_lagr_diff_c.shape[0]):
        for _ in range(ari_order - i):
            precip_lagr_diff_c[i] = _composite_convolution(
                precip_lagr_diff_c[i],
                kernels_1,
                interp_weights,
            )

    print("Estimating the ARI(p,1) parameters... ", end="", flush=True)

    if measure_time:
        starttime = time.time()

    # estimate ARI(p,1) parameters
    weights = _compute_window_weights(
        feature_coords,
        precip.shape[1],
        precip.shape[2],
        localization_window_radius,
    )

    if ari_order == 1:
        psi = _estimate_ar1_params(
            precip_lagr_diff_c[-1],
            precip_lagr_diff[-1],
            weights,
            interp_weights,
            num_workers=num_workers,
        )
    else:
        psi = _estimate_ar2_params(
            precip_lagr_diff_c[-2:],
            precip_lagr_diff[-1],
            weights,
            interp_weights,
            num_workers=num_workers,
        )
    forecast_gen["psi"] = psi

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

    # apply the ARI(p,1) model and integrate the differences
    precip_lagr_diff_c = _iterate_ar_model(precip_lagr_diff_c, psi)
    precip_forecast = precip_lagr[-2] + precip_lagr_diff_c[-1]
    precip_forecast[precip_forecast < 0.0] = 0.0

    print("Estimating the second convolution kernel... ", end="", flush=True)

    if measure_time:
        starttime = time.time()

    # estimate the second convolution kernels based on the forecast field
    # computed above
    kernels_2 = _estimate_convol_params(
        precip_forecast,
        precip[-1],
        convol_weights,
        mask_adv,
        kernel_type=kernel_type,
        num_workers=num_workers,
    )
    forecast_gen["kernels_2"] = kernels_2

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

    if measure_time:
        return forecast_gen, precip_lagr_diff, time.time() - starttime_init
    else:
        return forecast_gen, precip_lagr_diff
Initialize the deterministic LINDA nowcast model.
pysteps/nowcasts/linda.py