% Run code for NFR

clear all
clc

randn('seed',5);
rand('seed',5);

fid = fopen('fused_run.txt','w');

repeats = 10;
sigma = 0.1; % noise level
r = 6;
p = 0.5;
lambda_init = 0.1;

fopen(fid);

for index = 4:4:20
    iter_SD = 0; time_SD = 0; Fval_SD = 0;
    iter_SM = 0; time_SM = 0; Fval_SM = 0;
    iter_SM_1 = 0; time_SM_1 = 0; Fval_SM_1 = 0;
    
    n = 500*index;
    c2 = sigma*sqrt(n)/40;
    c1 = c2;
    
    % Generate signal
    length = 50*index;
    fprintf('n = %g, Start Now!\n',n);
    for repeat = 1: repeats
        J = randperm(n/length);
        I = sort(J(1:r),'ascend');
        x = zeros(n,1);
        for i = 1:r
            if randn > 0
                x(length*I(i) - 30*index - randi(3):length*I(i) ) = randi(3);
            else
                x(length*I(i) - 30*index - randi(3):length*I(i) ) = -randi(3);
            end
        end
        b = x + sigma*randn(n,1);
        
        t_sd = tic;
        fprintf(' ********** NFR by SDCAM **********\n')
        clear opts
        opts.x0 = ones(n,1);opts.tol = 1e-5;opts.freq = 100000;
        opts.maxiter = 10000;opts.disp = 1;opts.nonM = 5;lambda = lambda_init;
        % Initial function value
        Dx = diff(opts.x0);
        prox = lp_prox_mapping(Dx, c2*lambda, p);
        fval_1 = 0.5*norm(opts.x0 - b)^2 + c1*norm(opts.x0,1);
        fval_2 = norm(prox - Dx)^2/(2*lambda) + c2*norm(prox,p)^p;
        Fval_x0 = fval_1 + fval_2;
        % Initialization and solve by SDCAM
        xfeas = opts.x0;
        Fval_feas = Fval_x0;
        iter_1 = 0;
        while 1 == 1
            if Fval_feas < Fval_x0
                fprintf('swap initial point, lambda = %g\n', lambda)
                opts.x0 = xfeas;
            end
            [x,fval_1,iter,time] = SDCAM(b,c1,c2,p,lambda,opts);
            iter_1 = iter_1 + iter; 
            fprintf(' One inner loop done, lambda = %g\n', lambda)
            lambda = lambda/10;
            opts.x0 = x;
            opts.tol = max(opts.tol/1.5, 1e-6);
            if lambda < 1e-9 % note the difference
                break
            end
            Dx = diff(opts.x0);
            prox = lp_prox_mapping(Dx, c2*lambda, p);
            fval_2 = norm(prox - Dx)^2/(2*lambda) + c2*norm(prox,p)^p;
            Fval_x0 = fval_1 + fval_2;
        end
        time_1 = toc(t_sd);
        Fval_1 = fval_1 + c2*norm(diff(x),p)^p;
        
        iter_SD = iter_SD + iter_1;
        time_SD = time_SD + time_1;
        Fval_SD = Fval_SD + Fval_1;
        
        t_sm_1 = tic;
        fprintf(' ********** NFR by NPG_Smoothing with 1e-7**********\n')
        clear opts
        opts.x0 = ones(n,1);opts.tol = 1e-5;opts.freq = 100000;
        opts.maxiter = 10000;opts.disp = 1;opts.nonM = 5;lambda = lambda_init;
        iter_2 = 0; 
        while 1 == 1
            [x_1,fval_sm,iter,time] = NPG_sm(b,c1,c2,p,lambda,opts);
            iter_2 = iter_2 + iter; 
            fprintf(' One inner loop done, lambda = %g\n', lambda)
            lambda = lambda/10;
            opts.x0 = x_1;
            opts.tol = max(opts.tol/1.5, 1e-6);
            if lambda < 1e-7
                break
            end
        end
        time_2 = toc(t_sm_1);
        Fval_2 = fval_sm + c2*norm(diff(x_1),p)^p;
        
        iter_SM = iter_SM + iter_2;
        time_SM = time_SM + time_2;
        Fval_SM = Fval_SM + Fval_2;      
            
        t_sm_2 = tic;
        fprintf(' ********** NFR by NPG_Smoothing with 1e-8 **********\n')
        clear opts
        opts.x0 = ones(n,1);opts.tol = 1e-5;opts.freq = 100000;
        opts.maxiter = 10000;opts.disp = 1;opts.nonM = 5;lambda = lambda_init;
        iter_3 = 0; 
        while 1 == 1
            [x_2,fval_sm,iter,time] = NPG_sm(b,c1,c2,p,lambda,opts);
            iter_3 = iter_3 + iter; 
            fprintf(' One inner loop done, lambda = %g\n', lambda)
            lambda = lambda/10;
            opts.x0 = x_2;
            opts.tol = max(opts.tol/1.5, 1e-6);
            if lambda < 1e-8
                break
            end
        end
        time_3 = toc(t_sm_2);
        Fval_3 = fval_sm + c2*norm(diff(x_2),p)^p;
        
        iter_SM_1 = iter_SM_1 + iter_3;
        time_SM_1 = time_SM_1 + time_3;
        Fval_SM_1 = Fval_SM_1 + Fval_3;      
    end
    
    fprintf(fid, ' %7.0f & %7.0f & %7.0f & %7.0f & %5.1f & %5.1f & %5.1f & %6.5e & %6.5e & %6.5e\n', n, iter_SD/repeats, iter_SM/repeats, iter_SM_1/repeats, ...
        time_SD/repeats, time_SM/repeats, time_SM_1/repeats, Fval_SD/repeats, Fval_SM/repeats, Fval_SM_1/repeats);   
end

fclose(fid);


