function xout = andersen(x0, F, param)

% This is the AA-I safe algorithm adapted slightly

% Original: https://github.com/cvxgrp/nonexp_global_aa1 by 
% Paper  https://stanford.edu/~boyd/papers/pdf/scs_2.0_v_global.pdf 
% Zhang, Donoghue, Boyd 2018


beta     = 0.5;                              % biased towards more progressive iterations
itermax  = param.itermax;
xout     = x0;
rec      = struct();
tol      = param.tol; 
step     = 1/2; 

    mem_size = param.mem_size;
    theta    = param.theta;
    tau      = param.tau;
    D        = param.D;
    epsilon  = param.epsilon;
    
    Fx0      = F(x0);
    g0       = x0 - Fx0;
    
    Ubar     = norm(g0);
    Shat_mem = [];
    H_vecs1  = [];
    H_vecs2  = [];
    count    = 0;
    m        = 0;

    rec.restart = [];
    rec.safeguard = [];
    
    for i = 1 : itermax
        m = m + 1; % default increase of memory
        
        if mod(i, 1) == 1
            fprintf('iteration = %d\n', i);
        end
        
        if i == 1
        
            x1 = x0 + step*(Fx0 - x0);
            
        else
            
            x1 = x0 + step*(x0 - H_AA(g0, H_vecs1, H_vecs2) - x0);
        
        end
        
        % update using the AA trial update
        
        s0         = x1 - x0;
        [Fx1, err] = F(x1);
        g1         = x1 - Fx1;
        y0         = g1 - g0;
        
        %%% Safeguard checking
        
        Ubar0 = norm(g0);
        
        if Ubar0 <= D * Ubar * (count + 1)^(-1-epsilon) || i == 1
            
            count = count + 1;
            x0    = x1;
            Fx0   = Fx1;
            
        else
            
            rec.safeguard = [rec.safeguard, i];
            
            x1         = beta * x0 + (1 - beta) * Fx0;
            x0         = x1;
            [Fx0, err] = F(x0);
            
        end
        
        xout  = x0;
        g_1   = g0;                 % maintain g_{k-1} for Powell's trick
        g0    = x0 - Fx0;
        
 
        if err < tol, break, end
        
        %%% Restart checking
        if m <= mem_size
            if ~isempty(Shat_mem) % only do when nonempty memory
                s0hat = s0 - ((Shat_mem * s0)' * Shat_mem)';
            else
                %fprintf('iter=%d: no orthogonalization\n', i);
                s0hat = s0;
            end
            if norm(s0hat) < tau * norm(s0) % restart
                rec.restart = [rec.restart, i];
                s0hat = s0;
                m = 1;
                Shat_mem = [];
                H_vecs1 = [];
                H_vecs2 = [];
            end
        else % memory exceeds
            s0hat = s0;
            m = 1;
            Shat_mem = [];
            H_vecs1 = [];
            H_vecs2 = [];
        end
        
        Shat_mem = [Shat_mem; s0hat' / norm(s0hat)];
        
        %%% Powell regularization
        gamma0 = s0hat' * H_AA(y0, H_vecs1, H_vecs2) / norm(s0hat)^2;
        theta0 = phi(gamma0, theta); 
        y0tilde = theta0 * y0 - (1-theta0) * g_1;
        
        %%% Update H_vecs
        Hytilde = H_AA(y0tilde, H_vecs1, H_vecs2);
        hvec1   = s0 - Hytilde;
        hvec2   = H_AAt(s0hat, H_vecs1, H_vecs2) / (s0hat' * Hytilde);
        H_vecs1 = [H_vecs1, hvec1];
        H_vecs2 = [H_vecs2, hvec2];
        
    end
    
