function [x2, fv, info] = RFedAGS(problem, x, options)

    global ybases;

    stepsize = options.stepsize;
    fv  = problem.cost(x);
    gf  = problem.grad(x);
    ngf = problem.M.norm(x, gf);    
    fvs  = fv;
    ngfs = ngf;
    x1 = x; x2 = x;
    

    if ~isfield(options, 'stepsize_type')
        options.stepsize_type = 'fix';
    elseif ~isfield(options, 'lambda')
        options.lambda = 1;
    end
    
    if options.verbosity > 0
%     fprintf('iter:%6d,fv:%.6e(%.6e),ngf:%.6e,ngf/ngf0:%.6e,stepsize:%.6e\n', 0, fvs(end),fvs(end)/fvs(1),ngfs(end),ngfs(end)/ngfs(1),options.stepsize);
    fprintf('iter:%6d|fv:%.6e(%.6e)|ngf:%.6e|ngf/ngf0:%.6e|stepsize:%.6e\n', 0, fvs(end),fvs(end)/fvs(1),ngfs(end),ngfs(end)/ngfs(1),options.stepsize);
    end
    t = 0; ct = 0;   
    if options.checkperiod > 0
        xtmp = cell(options.maxiter / options.checkperiod, 1);    
    else
        xtmp = cell(options.maxiter, 1);
    end

    if options.option == 1
        num_subagent = problem.nagents;         
    elseif options.option == 2
        num_subagent = options.s;
    end  

    localzeta  = cell(num_subagent, 1);        
    num_gfs = 0; CPUtime = 0; times = 0; cc = 1;
    while t < options.maxiter
        t1 = tic;

        St = randperm(problem.nagents, num_subagent);

        if (strcmp(options.stepsize_type, 'decay')) && (mod(t, options.decaysteps) == 0) && (t>0)                                        
            ct = ct + 1;
            options.stepsize = stepsize / (options.lambda + ct);
        end             
                        
        time = toc(t1);
        % active agents locally update
        local_times = zeros(length(St), 1);        
%         tmp = zeros(1, problem.nagents);
        ybases = [];
        for i = reshape(St, 1, length(St)) % must be a row array
            % single agent 
            t2 = tic;
            prob = problem.agent{i};
            opts = options.agent{i};
            opts.stepsize_type = 'fix';
            opts.stepsize_init = options.stepsize;
            opts.option = options.option;            
            [localzeta{i}, num_gf] = localRSGD(prob, x1, opts);            
%             tmp(i) = norm(x1'*localzeta{i} + localzeta{i}'*x1,'fro');

            local_times(i) = toc(t2);
        end
        time = time + max(local_times);
        num_gfs(end+1) = num_gfs(end) + num_gf(end);
%         fprintf('-------------------------\n')
%         disp(tmp)
        
        % the server aggratates
        t3 = tic;
        x2 = Aggregation(x1, localzeta, St, prob.M);                
        t = t+1;
        x1 = x2;
        time = time + toc(t3);
        times = times + time;
        if mod(t, options.checkperiod) == 0            
            xtmp{cc} = x2; cc = cc+1;
            fv   = problem.cost(x2);
            gf   = problem.grad(x2);        
            ngf  = problem.M.norm(x2, gf);
            fvs  = [fvs;fv];
            ngfs = [ngfs;ngf];
            CPUtime =[CPUtime; CPUtime(end) + times];
            times = 0;
            if options.verbosity > 0
            fprintf('iter:%6d|fv:%.6e(%.6e)|ngf:%.6e|ngf/ngf0:%.6e|stepsize:%.6e\n', t, fvs(end),fvs(end)/fvs(1),ngfs(end),ngfs(end)/ngfs(1),options.stepsize);        
            end
        end                 
    end 

    info.fvs  = fvs;    
    info.ngfs = ngfs;
    info.xs   = xtmp;
    info.CPUtime = CPUtime;
    info.num_gfs = reshape(num_gfs, length(num_gfs),1);
    clear ybases;   
end

function [zeta, ngfs] = localRSGD(prob, x, opts)
    maxiter = opts.maxiter;    
    stepsize = opts.stepsize_init;
    x1 = x;
    zeta = zeros(size(x1));
    iter = 0;
    ngfs = 0;
    while iter < maxiter                
        idx_batch = randperm(prob.ncostterms, opts.batchsize);        
        eta = - stepsize * prob.partialgrad(x1, idx_batch); ngf = opts.batchsize;

        % update local iterate       
        x2 = prob.M.retr(x1, eta);        

        % gradient stream    
        if iter == 0
            zeta = eta;
        else
            zeta = zeta + prob.M.transp(x1, x, eta);
        end
        x1 = x2;
        iter = iter + 1;      
        ngfs(end+1) = ngfs(end) + ngf;
    end
end

function output = Aggregation(x, localzeta, St, manifold)
    N = length(St);    
    tmp = localzeta{St(1)};
    for i = 2:N
        tmp = tmp + localzeta{St(i)};
%         tmp = manifold.lincomb(x, 1, tmp, 1, localzeta{i});
    end

    output = manifold.retr(x, tmp, 1/N);
end


