def line_search(oracle, x_k, t, d_k, alpha_0=1., c1=1e-4):
    """
    Returns
    -------
    alpha : float or None if failure
        Chosen step size
    """

    phi = lambda alpha: oracle.f_t_directional(x_k, t, d_k, alpha)
    derphi = lambda alpha: oracle.grad_directional(x_k, t, d_k, alpha)
    phi_0 = phi(0.)
    derphi_0 = derphi(0.)
    alpha = alpha_0

    while phi(alpha) > phi_0 + c1 * alpha * derphi_0:
        alpha = alpha / 2.

    return alpha


def newton_barrier(oracle, common_x_0, t, tolerance=1e-5,
                   max_iter=100, theta=0.99, c1=1e-4):

    common_x_k = common_x_0
    start_grad_norm = norm(oracle.common_grad(common_x_0, t))
    x_k, u_k = np.split(common_x_k, 2)
    prev_x_k, prev_u_k = x_k, u_k

    for iter_num in range(max_iter + 1):

        if not np.any(np.isfinite(common_x_k)):
            return prev_x_k, prev_u_k

        try:
            common_grad = oracle.common_grad(common_x_k, t)
            common_grad_norm = norm(common_grad)
            if common_grad_norm ** 2 < tolerance * start_grad_norm ** 2:
                break

            common_hess = oracle.common_hess(common_x_k, t)

        except ValueError:
            return prev_x_k, prev_u_k

        try:
            if sp.issparse(common_hess):
                common_d_k = spla.spsolve(common_hess, -common_grad)
            else:
                c, lower = cho_factor(common_hess)
                common_d_k = cho_solve((c, lower), -common_grad)
        except LinAlgError:
            return prev_x_k, prev_u_k

        d_x, d_u = np.split(common_d_k, 2)
        x_k, u_k = np.split(common_x_k, 2)
        prev_x_k = x_k
        prev_u_k = u_k

        dx_du = d_x - d_u
        _dx_du = -d_x - d_u

        pos_mask = dx_du > 0
        neg_mask = _dx_du > 0

        alphas_first = theta * (u_k - x_k)[pos_mask] / dx_du[pos_mask]
        alpha_second = theta * (u_k + x_k)[neg_mask] / _dx_du[neg_mask]

        alpha_max_to_use = np.concatenate(([1.], alphas_first, alpha_second))
        alpha_max_to_use = alpha_max_to_use.min()

        alpha = line_search(oracle=oracle, x_k=common_x_k, t=t,
                            d_k=common_d_k, alpha_0=alpha_max_to_use, c1=c1)

        common_x_k = common_x_k + alpha * common_d_k

    x_star, u_star = np.split(common_x_k, 2)
    return x_star, u_star


def barrier_method_lasso(A, b, reg_coef, x_0, u_0, tolerance=1e-5,
                         tolerance_inner=1e-8, max_iter=100,
                         max_iter_inner=20, t_0=1, gamma=10,
                         c1=1e-4, lasso_duality_gap=None,
                         trace=False, display=False):
    """
    Log-barrier method for solving the problem:
        minimize    f(x, u) := 1/2 * ||Ax - b||_2^2 + reg_coef * \sum_i u_i
        subject to  -u_i <= x_i <= u_i.

    The method constructs the following barrier-approximation of the problem:
        phi_t(x, u) := t * f(x, u) - sum_i( log(u_i + x_i) + log(u_i - x_i) )
    and minimize it as unconstrained problem by Newton's method.

    In the outer loop `t` is increased and we have a sequence of approximations
        { phi_t(x, u) } and solutions { (x_t, u_t)^{*} } which converges in `t`
    to the solution of the original problem.

    Parameters
    ----------
    A : np.array
        Feature matrix for the regression problem.
    b : np.array
        Given vector of responses.
    reg_coef : float
        Regularization coefficient.
    x_0 : np.array
        Starting value for x in optimization algorithm.
    u_0 : np.array
        Starting value for u in optimization algorithm.
    tolerance : float
        Epsilon value for stopping criterion.
    max_iter : int
        Maximum number of iterations for interior point method.
    max_iter_inner : int
        Maximum number of iterations for inner Newton's method.
    t_0 : float
        Starting value for `t`.
    gamma : float
        Multiplier for changing `t` during the iterations:
        t_{k + 1} = gamma * t_k.
    c1 : float
        Armijo's constant for line search in Newton's method.
    lasso_duality_gap : callable object or None.
        If calable the signature is lasso_duality_gap(x, Ax_b, ATAx_b, b, regcoef)
        Returns duality gap value for esimating the progress of method.
    trace : bool
        If True, the progress information is appended into history dictionary
        during training. Otherwise None is returned instead of history.
    display : bool
        If True, debug information is displayed during optimization.
        Printing format is up to a student and is not checked in any way.

    Returns
    -------
    (x_star, u_star) : tuple of np.array
        The point found by the optimization procedure.
    message : string
        "success" or the description of error:
            - 'iterations_exceeded': if after max_iter iterations of the method x_k still doesn't satisfy
                the stopping criterion.
            - 'computational_error': in case of getting Infinity or None value during the computations.
    history : dictionary of lists or None
        Dictionary containing the progress information or None if trace=False.
        Dictionary has to be organized as follows:
            - history['time'] : list of floats, containing time in seconds passed from the start of the method
            - history['func'] : list of function values f(x_k) on every step of the algorithm
            - history['duality_gap'] : list of duality gaps
            - history['x'] : list of np.arrays, containing the trajectory of the algorithm. ONLY STORE IF x.size <= 2
    """
    history = defaultdict(list) if trace else None
    converged = False
    start_time = time()
    x_k = x_0
    u_k = u_0
    t_k = t_0
    theta = 0.99

    oracle = oracles.LassoBarrierOracul(A, b, reg_coef)

    for iter_num in range(max_iter + 1):

        gap = lasso_duality_gap(x_k, oracle.Ax_b(x_k),
                                oracle.ATAx_b(x_k), b, reg_coef)

        f_k = oracle.f(x_k, u_k)

        if trace:
            history['time'].append(time() - start_time)
            history['func'].append(f_k)
            history['duality_gap'].append(gap)
            if x_k.size <= 2:
                history['x'].append(x_k)

        if display:
            print("f_k = {}".format(f_k))

        if gap < tolerance:
            converged = True
            break

        common_x_k = np.concatenate((x_k, u_k), axis=0)
        x_k, u_k = newton_barrier(oracle, common_x_k, t_k,
                                  tolerance=tolerance_inner,
                                  max_iter=max_iter_inner, theta=theta, c1=c1)
        t_k *= gamma

    if converged:
        return (x_k, u_k), 'success', history
    else:
        return (x_k, u_k), 'iterations_exceeded', history