How to use the scvelo.tools.utils.test_bimodality function in scvelo

# initialize beta and gamma from extreme quantiles of s
        if True:
            weights_s = s_w >= np.percentile(s_w, perc, axis=0)
        else:
            us_norm = s_w / np.clip(np.max(s_w, axis=0), 1e-3, None) + u_w / np.clip(np.max(u_w, axis=0), 1e-3, None)
            weights_s = us_norm >= np.percentile(us_norm, perc, axis=0)

        beta, gamma = 1, linreg(convolve(u_w, weights_s), convolve(s_w, weights_s))

        u_inf, s_inf = u_w[weights_s].mean(), s_w[weights_s].mean()
        u0_, s0_ = u_inf, s_inf
        alpha = np.mean([s_inf * gamma, u_inf * beta])  # np.mean([s0_ * gamma, u0_ * beta])

        # initialize switching from u quantiles and alpha from s quantiles
        tstat_u, pval_u, means_u = test_bimodality(u_w, kde=True)
        tstat_s, pval_s, means_s = test_bimodality(s_w, kde=True)
        self.pval_steady = max(pval_u, pval_s)
        self.u_steady = means_u[1]
        self.s_steady = means_s[1]

        if self.pval_steady < .1:
            u_inf = np.mean([u_inf, self.u_steady])
            s_inf = np.mean([s_inf, self.s_steady])
            alpha = s_inf * gamma
            beta = alpha / u_inf

            weights_u = u_w >= np.percentile(u_w, perc, axis=0)
            u0_, s0_ = u_w[weights_u].mean(), s_w[weights_u].mean()

        # alpha, beta, gamma = np.array([alpha, beta, gamma]) * scaling
        t_ = tau_inv(u0_, s0_, 0, 0, alpha, beta, gamma)

elif gamma > 1.5 / scaling:
            gamma /= 1.2

        u_inf, s_inf = u_w[weights_u | weights_s].mean(), s_w[weights_s].mean()
        u0_, s0_ = u_inf, s_inf
        alpha = u_inf * beta
        # np.mean([s_inf * gamma, u_inf * beta])  # np.mean([s0_ * gamma, u0_ * beta])

        if self.init_vals is not None:  # just to test different EM start
            self.init_vals = np.array(self.init_vals)
            alpha, beta, gamma = np.array([alpha, beta, gamma]) * self.init_vals

        # initialize switching from u quantiles and alpha from s quantiles
        try:
            tstat_u, pval_u, means_u = test_bimodality(u_w, kde=True)
            tstat_s, pval_s, means_s = test_bimodality(s_w, kde=True)
        except:
            logg.warn("skipping bimodality check for", self.gene)
            tstat_u, tstat_s, pval_u, pval_s = 0, 0, 1, 1
            means_u, means_s = [0, 0], [0, 0]

        self.pval_steady = max(pval_u, pval_s)
        self.steady_u = means_u[1]
        self.steady_s = means_s[1]

        if self.pval_steady < 1e-3:
            u_inf = np.mean([u_inf, self.steady_u])
            alpha = gamma * s_inf
            beta = alpha / u_inf
            u0_, s0_ = u_inf, s_inf

        # alpha, beta, gamma = np.array([alpha, beta, gamma]) * scaling

gamma *= 1.2
        elif gamma > 1.5 / scaling:
            gamma /= 1.2

        u_inf, s_inf = u_w[weights_u | weights_s].mean(), s_w[weights_s].mean()
        u0_, s0_ = u_inf, s_inf
        alpha = u_inf * beta
        # np.mean([s_inf * gamma, u_inf * beta])  # np.mean([s0_ * gamma, u0_ * beta])

        if self.init_vals is not None:  # just to test different EM start
            self.init_vals = np.array(self.init_vals)
            alpha, beta, gamma = np.array([alpha, beta, gamma]) * self.init_vals

        # initialize switching from u quantiles and alpha from s quantiles
        try:
            tstat_u, pval_u, means_u = test_bimodality(u_w, kde=True)
            tstat_s, pval_s, means_s = test_bimodality(s_w, kde=True)
        except:
            logg.warn("skipping bimodality check for", self.gene)
            tstat_u, tstat_s, pval_u, pval_s = 0, 0, 1, 1
            means_u, means_s = [0, 0], [0, 0]

        self.pval_steady = max(pval_u, pval_s)
        self.steady_u = means_u[1]
        self.steady_s = means_s[1]

        if self.pval_steady < 1e-3:
            u_inf = np.mean([u_inf, self.steady_u])
            alpha = gamma * s_inf
            beta = alpha / u_inf
            u0_, s0_ = u_inf, s_inf

# initialize beta and gamma from extreme quantiles of s
        if True:
            weights_s = s_w >= np.percentile(s_w, perc, axis=0)
        else:
            us_norm = s_w / np.clip(np.max(s_w, axis=0), 1e-3, None) + u_w / np.clip(np.max(u_w, axis=0), 1e-3, None)
            weights_s = us_norm >= np.percentile(us_norm, perc, axis=0)

        beta, gamma = 1, linreg(convolve(u_w, weights_s), convolve(s_w, weights_s))

        u_inf, s_inf = u_w[weights_s].mean(), s_w[weights_s].mean()
        u0_, s0_ = u_inf, s_inf
        alpha = np.mean([s_inf * gamma, u_inf * beta])  # np.mean([s0_ * gamma, u0_ * beta])

        # initialize switching from u quantiles and alpha from s quantiles
        tstat_u, pval_u, means_u = test_bimodality(u_w, kde=True)
        tstat_s, pval_s, means_s = test_bimodality(s_w, kde=True)
        self.pval_steady = max(pval_u, pval_s)
        self.u_steady = means_u[1]
        self.s_steady = means_s[1]

        if self.pval_steady < .1:
            u_inf = np.mean([u_inf, self.u_steady])
            s_inf = np.mean([s_inf, self.s_steady])
            alpha = s_inf * gamma
            beta = alpha / u_inf

            weights_u = u_w >= np.percentile(u_w, perc, axis=0)
            u0_, s0_ = u_w[weights_u].mean(), s_w[weights_u].mean()

        # alpha, beta, gamma = np.array([alpha, beta, gamma]) * scaling
        t_ = tau_inv(u0_, s0_, 0, 0, alpha, beta, gamma)