How to use the scvelo.tools.velocity.Velocity function in scvelo

for cat in categories:
            groups = cat if cat is not None else groups

            cell_subset = groups_to_bool(adata, groups, groupby)
            _adata = adata if groups is None else adata[cell_subset]
            velo = Velocity(_adata, groups_for_fit=groups_for_fit, groupby=groupby, constrain_ratio=constrain_ratio,
                            min_r2=min_r2, r2_adjusted=r2_adjusted, use_raw=use_raw)
            velo.compute_deterministic(fit_offset=fit_offset, perc=perc)

            if mode == 'stochastic':
                if filter_genes and len(set(velo._velocity_genes)) > 1:
                    adata._inplace_subset_var(velo._velocity_genes)
                    residual = velo._residual[:, velo._velocity_genes]
                    _adata = adata if groups is None else adata[cell_subset]
                    velo = Velocity(_adata, residual=residual, groups_for_fit=groups_for_fit, groupby=groupby, constrain_ratio=constrain_ratio)
                velo.compute_stochastic(fit_offset, fit_offset2, mode, perc=perc)

            write_residuals(adata, vkey, velo._residual, cell_subset)
            write_residuals(adata, f'variance_{vkey}', velo._residual2, cell_subset)
            write_pars(adata, vkey, velo.get_pars(), velo.get_pars_names(), add_key=cat)

            if filter_genes and len(set(velo._velocity_genes)) > 1:
                adata._inplace_subset_var(velo._velocity_genes)

    else:
        raise ValueError('Mode can only be one of these: steady_state, deterministic, stochastic or dynamical.')

    if f'{vkey}_genes' in adata.var.keys() and np.sum(adata.var[f'{vkey}_genes']) < 10:
        logg.warn('Too few genes are selected as velocity genes. '
                  'Consider setting a lower threshold for min_r2 or min_likelihood.')

adata.layers[f'{vkey}_u'] = np.ones(adata.shape) * np.nan
        adata.layers[f'{vkey}_u'][:, gene_subset] = wt

        if filter_genes and len(set(vgenes)) > 1:
            adata._inplace_subset_var(vgenes)

    elif mode in {'steady_state', 'deterministic', 'stochastic'}:
        categories = adata.obs[groupby].cat.categories \
            if groupby is not None and groups is None and groups_for_fit is None else [None]

        for cat in categories:
            groups = cat if cat is not None else groups

            cell_subset = groups_to_bool(adata, groups, groupby)
            _adata = adata if groups is None else adata[cell_subset]
            velo = Velocity(_adata, groups_for_fit=groups_for_fit, groupby=groupby, constrain_ratio=constrain_ratio,
                            min_r2=min_r2, r2_adjusted=r2_adjusted, use_raw=use_raw)
            velo.compute_deterministic(fit_offset=fit_offset, perc=perc)

            if mode == 'stochastic':
                if filter_genes and len(set(velo._velocity_genes)) > 1:
                    adata._inplace_subset_var(velo._velocity_genes)
                    residual = velo._residual[:, velo._velocity_genes]
                    _adata = adata if groups is None else adata[cell_subset]
                    velo = Velocity(_adata, residual=residual, groups_for_fit=groups_for_fit, groupby=groupby, constrain_ratio=constrain_ratio)
                velo.compute_stochastic(fit_offset, fit_offset2, mode, perc=perc)

            write_residuals(adata, vkey, velo._residual, cell_subset)
            write_residuals(adata, f'variance_{vkey}', velo._residual2, cell_subset)
            write_pars(adata, vkey, velo.get_pars(), velo.get_pars_names(), add_key=cat)

            if filter_genes and len(set(velo._velocity_genes)) > 1:

kwargs_.update(**kwargs)

        if 'residuals' in mode:
            u, s = get_reads(vdata, use_raw=adata.uns['recover_dynamics']['use_raw'])
            if kwargs_['fit_basal_transcription']:
                u, s = u - u0, s - s0
            o = vdata.layers['fit_t'] < t_
            vt = (u * beta - s * gamma)  # ds/dt
            wt = (alpha * o - beta * u) * scaling  # du/dt
        else:
            vt, wt = get_divergence(vdata, mode='velocity', **kwargs_)

        vgenes = adata.var.fit_likelihood > min_likelihood
        if min_r2 is not None:
            if 'fit_r2' not in adata.var.keys():
                velo = Velocity(adata, groups_for_fit=groups_for_fit, groupby=groupby, constrain_ratio=constrain_ratio,
                                min_r2=min_r2, use_raw=use_raw)
                velo.compute_deterministic(fit_offset=fit_offset, perc=perc)
                adata.var['fit_r2'] = velo._r2
            vgenes &= adata.var.fit_r2 > min_r2

        lb, ub = np.nanpercentile(adata.var.fit_scaling, [10, 90])
        vgenes = vgenes & (adata.var.fit_scaling > np.min([lb, .03])) & (adata.var.fit_scaling < np.max([ub, 3]))

        adata.var[f'{vkey}_genes'] = vgenes

        adata.layers[vkey] = np.ones(adata.shape) * np.nan
        adata.layers[vkey][:, gene_subset] = vt


        adata.layers[f'{vkey}_u'] = np.ones(adata.shape) * np.nan
        adata.layers[f'{vkey}_u'][:, gene_subset] = wt

adata.uns["recover_dynamics"] = {
        "fit_connected_states": fit_connected_states,
        "fit_basal_transcription": fit_basal_transcription,
        "use_raw": use_raw,
    }

    if isinstance(var_names, str) and var_names not in adata.var_names:
        if var_names in adata.var.keys():
            var_names = adata.var_names[adata.var[var_names].values]
        elif use_raw or var_names == "all":
            var_names = adata.var_names
        elif "_genes" in var_names:
            from .velocity import Velocity

            velo = Velocity(adata, use_raw=use_raw)
            velo.compute_deterministic(perc=[5, 95])
            var_names = adata.var_names[velo._velocity_genes]
            adata.var["fit_r2"] = velo._r2
        else:
            raise ValueError("Variable name not found in var keys.")
    if not isinstance(var_names, str):
        var_names = list(np.ravel(var_names))

    var_names = make_unique_list(var_names, allow_array=True)
    var_names = np.array([name for name in var_names if name in adata.var_names])
    if len(var_names) == 0:
        raise ValueError("Variable name not found in var keys.")
    if n_top_genes is not None and len(var_names) > n_top_genes:
        X = adata[:, var_names].layers[("spliced" if use_raw else "Ms")]
        var_names = var_names[np.argsort(np.sum(X, 0))[::-1][:n_top_genes]]
    if return_model is None:

-------
    Returns or updates `adata`
    """
    adata = data.copy() if copy else data

    if "Ms" not in adata.layers.keys() or "Mu" not in adata.layers.keys():
        use_raw = True
    if isinstance(var_names, str) and var_names not in adata.var_names:
        if var_names in adata.var.keys():
            var_names = adata.var_names[adata.var[var_names].values]
        elif use_raw or var_names == "all":
            var_names = adata.var_names
        elif "_genes" in var_names:
            from .velocity import Velocity

            velo = Velocity(adata, use_raw=use_raw)
            velo.compute_deterministic(perc=[5, 95])
            var_names = adata.var_names[velo._velocity_genes]
            adata.var["fit_r2"] = velo._r2
        else:
            raise ValueError("Variable name not found in var keys.")
    if not isinstance(var_names, str):
        var_names = list(np.ravel(var_names))

    var_names = make_unique_list(var_names, allow_array=True)
    var_names = [name for name in var_names if name in adata.var_names]
    if len(var_names) == 0:
        raise ValueError("Variable name not found in var keys.")
    if return_model is None:
        return_model = len(var_names) < 5

    # fit dynamical model first, if not done yet.