def _get_efron_values_batch(self, X, T, E, weights, beta):
        n, d = X.shape
        hessian = np.zeros((d, d))
        gradient = np.zeros((d,))
        log_lik = 0
        risk_phi, tie_phi = 0, 0
        risk_phi_x, tie_phi_x = np.zeros((d,)), np.zeros((d,))
        risk_phi_x_x, tie_phi_x_x = np.zeros((d, d)), np.zeros((d, d))
        _, counts = np.unique(-T, return_counts=True)
        scores = weights * np.exp(np.dot(X, beta))
        pos = n
        for count_of_removals in counts:
            slice_ = slice(pos - count_of_removals, pos)
            X_at_t = X[slice_]
            weights_at_t = weights[slice_]
            phi_i = scores[slice_, None]
            phi_x_i = phi_i * X_at_t
            phi_x_x_i = np.dot(X_at_t.T, phi_x_i)
            risk_phi = risk_phi + array_sum_to_scalar(phi_i)
            risk_phi_x = risk_phi_x + matrix_axis_0_sum_to_array(phi_x_i)
            risk_phi_x_x = risk_phi_x_x + phi_x_x_i
            deaths = E[slice_]
            tied_death_counts = array_sum_to_scalar(deaths.astype(int))
            if tied_death_counts == 0:
                pos -= count_of_removals
                continue
            xi_deaths = X_at_t[deaths]
            weights_deaths = weights_at_t[deaths]
            x_death_sum = matrix_axis_0_sum_to_array(weights_deaths[:, None] * xi_deaths)
            weight_count = array_sum_to_scalar(weights_deaths)
            weighted_average = weight_count / tied_death_counts
            if tied_death_counts > 1:
                tie_phi = array_sum_to_scalar(phi_i[deaths])
                tie_phi_x = matrix_axis_0_sum_to_array(phi_x_i[deaths])
                tie_phi_x_x = np.dot(xi_deaths.T, phi_i[deaths] * xi_deaths)
                increasing_proportion = np.arange(tied_death_counts) / tied_death_counts
                denom = 1.0 / (risk_phi - increasing_proportion * tie_phi)
                numer = risk_phi_x - np.outer(increasing_proportion, tie_phi_x)
                a1 = np.einsum("ab,i->ab", risk_phi_x_x, denom) - np.einsum(
                    "ab,i->ab", tie_phi_x_x, increasing_proportion * denom
                )
            else:
                denom = 1.0 / np.array([risk_phi])
                numer = risk_phi_x
                a1 = risk_phi_x_x * denom
            summand = numer * denom[:, None]
            a2 = summand.T.dot(summand)
            gradient = gradient + x_death_sum - weighted_average * summand.sum(0)
            log_lik = log_lik + np.dot(x_death_sum, beta) + weighted_average * np.log(denom).sum()
            hessian = hessian + weighted_average * (a2 - a1)
            pos -= count_of_removals
        return hessian, gradient, log_lik