def make_feasible(x_0, c):
    """ Returns a feasible x with respect to the constraints g_3, g_4, g_5, and g_6.

    :param np.array x_0: (n,) Initial, infeasible point.
    :param dict c: Constants defining the constraints.
    :return np.array: (n,) Initial, feasible point.

    Notes:
        - x_0 should follow: x0 = [200000, 355, 1, 1, 1, 1] for one backstress
        - x_0 should follow: x0 = [200000, 355, 1, 1, 1, rho_gamma_12_inf, 1, 1] for two backstresses
    """
    n_backstresses = int((len(x_0) - 4) // 2)
    # Get the average of the bounds
    rho_yield_avg = 0.5 * (c['rho_yield_inf'] + c['rho_yield_sup'])
    rho_iso_avg = 0.5 * (c['rho_iso_inf'] + c['rho_iso_sup'])
    rho_gamma_avg = 0.5 * (c['rho_gamma_inf'] + c['rho_gamma_sup'])
    # Calculate the feasible set
    gamma2_0 = 1. / rho_gamma_avg * x_0[5]
    c1_0 = -x_0[5] * (-1. + rho_iso_avg) * (-1 + rho_yield_avg) * x_0[1]
    q_inf_0 = -c1_0 * rho_iso_avg / (x_0[5] * (-1. + rho_iso_avg))
    b_0 = 1. / rho_gamma_avg * x_0[5]
    if n_backstresses == 2:
        x_0[[2, 3, 4, 7]] = [q_inf_0, b_0, c1_0, gamma2_0]
    elif n_backstresses == 1:
        x_0[[2, 3, 4]] = [q_inf_0, b_0, c1_0]
    return x_0