function code_PROTDPAPNV
clc;
close all;
clear all;
format short

%% random seed
rng(213)



%% the setup of test image
N=256;
image_type=4;

% 1  Test image: cameraman
% 2  Test image: barbara
% 3  Test image: phantom
% 4  Test image: RSCB (random signed cameraman + i random signed barbara)
% 5  Test image: RPP (phantom with random phases)

%% the number of random illumination & oversampling ratio of Fourier intensity data
if (image_type<3.5)
    number_illumination=1;
    oversampling_ratio=4;
else
    number_illumination=2;
    oversampling_ratio=4;
end

if (image_type==1)
cameraman=double(imread('image_cameraman256.png'));
test_image=cameraman(1:N,1:N);
elseif (image_type==2)
barbara=double(imread('image_barbara260.png'));
test_image=barbara(1:N,1:N);
elseif (image_type==3)
test_image=double(imread('phantom.png')); 
elseif (image_type==4)
cameraman=double(imread('image_cameraman256.png'));
barbara=double(imread('image_barbara260.png')); 
random_sign_matrix_real=sign(round(rand(N))-0.5);
random_sign_matrix_imag=sign(round(rand(N))-0.5);
test_image=random_sign_matrix_real.*cameraman(1:N,1:N)+1i*random_sign_matrix_imag.*barbara(1:N,1:N);
elseif (image_type==5)
phantom=double(imread('phantom.png'));
test_image=exp(1i*2*pi*rand(N,N)).*phantom(1:N,1:N); 
end

%% methods for generating initial guesses  (1 : Execute  0 : Idle)
initial_selection=3;

if (initial_selection==3)
spectral_method=0;
truncated_spectral_method=0;
null_vector_method=1;
elseif (initial_selection==2)
spectral_method=0;
truncated_spectral_method=1;
null_vector_method=0;
elseif (initial_selection==1)
spectral_method=1;
truncated_spectral_method=0;
null_vector_method=0;
end

%% algorithms
max_iterate=300;
max_iterate_powermethod=1000;
serialAP=1;
parallelAP=1;
WF=1;

%% zeropadding for oversampling the fourier magnitude data 
zero_padding_size=ceil(N*((oversampling_ratio)^0.5-1)/2);
z=zero_padding_size;

filter=zeros(N+2*z);
filter(z+1:z+N,z+1:z+N)=1;
    
extension=zeros(N+2*z); 
extension(z+1:z+N,z+1:z+N)=test_image;
test_image=extension;




if (image_type<3.5)
figure;
imshow(uint8(test_image(z+1:z+N,z+1:z+N)));
title ('test image');
elseif (image_type==4)   
figure;
imshow(uint8(abs(real(test_image(z+1:z+N,z+1:z+N)))));
title ('real part of the test image');
figure;
imshow(uint8(abs(imag(test_image(z+1:z+N,z+1:z+N)))));
title ('imaginary part of the test image');
elseif (image_type==5)
figure;
imshow(uint8(abs(test_image(z+1:z+N,z+1:z+N))));
title ('absolute value of the test image');
end






%% random phase masks
mask_box=zeros(N+2*z,N+2*z,number_illumination);
for mask_index=1:number_illumination
    if (mask_index==1)
    mask_box(:,:,mask_index)=exp(1i*2*pi.*rand(N+2*z)); 
    elseif (mask_index==2)
    mask_box(:,:,mask_index)=exp(1i*2*pi.*rand(N+2*z)); 
    end
end


 
 %% data (absolute value)
 mask_test_image=zeros(N+2*z,N+2*z,number_illumination);
 data_box=zeros(N+2*z,N+2*z,number_illumination);
for mask_index=1:number_illumination
 mask_test_image(:,:,mask_index)=mask_box(:,:,mask_index).*test_image(:,:); 
 data_box(:,:,mask_index)=abs(fft2(mask_test_image(:,:,mask_index))); 
end


square_norm_x=0;
for mask_index=1:number_illumination
   square_norm_x=square_norm_x+norm(data_box(:,:,mask_index),'fro')^2/(N+2*z)^2;
end   
square_norm_x=square_norm_x/number_illumination;
   







%% initialization
if (image_type<3.5)
random_matrix_1=rand(N+2*z);
random_initialization=(square_norm_x^0.5)*random_matrix_1/norm(random_matrix_1,'fro');
else    
random_matrix_2=rand(N+2*z).*exp(1i*2*pi*rand(N+2*z));
random_initialization=(square_norm_x^0.5)*random_matrix_2/norm(random_matrix_2,'fro');
end


%% spectral method 
if (spectral_method==1)
initialization_spectral=random_initialization;
relative_residual_spectral_initial=zeros(1,max_iterate_powermethod);
    for tt = 1 : max_iterate_powermethod 
    temp1=zeros(N+2*z);
      for mask_index=1:number_illumination
      temp2 =(data_box(:,:,mask_index).^2).*fft2(mask_box(:,:,mask_index).*initialization_spectral(:,:)); 
      IFFT2_temp2=ifft2(temp2);
      temp1=temp1+IFFT2_temp2(:,:)./mask_box(:,:,mask_index)/number_illumination;
      end
    
      if (image_type==1 || image_type==2 || image_type==3)
      temp1=max(real(temp1),0);
      end
    
    initialization_spectral= sqrt(square_norm_x)*temp1.*filter/norm(temp1.*filter,'fro');
    %%%%%%%%%%% compute the relative residual %%%%%%%%%%
       residual=0;
       for mask_index=1:number_illumination
       mask_initialization_spectral= mask_box(:,:,mask_index).*initialization_spectral;  
       FFT_mask_initialization_spectral=fft2(mask_initialization_spectral);
       abs_FFT_mask_initialization_spectral=abs(FFT_mask_initialization_spectral);
       residual=residual+norm(data_box(:,:,mask_index)-abs_FFT_mask_initialization_spectral(:,:),'fro')^2/norm(data_box(:,:,mask_index),'fro')^2; 
       end
       relative_residual_spectral_initial(tt)=sqrt(residual/number_illumination);
       %fprintf('relative residual of the initial guess via the spectral method  %d=   %d; \n', tt, relative_residual_spectral_initial(tt))
    end 

      figure;
      plot(1:max_iterate_powermethod,relative_residual_spectral_initial,'r')
      xlabel('iteration','FontSize',12,'FontWeight','bold')
      ylabel('relative residual for the spectral method','FontSize',12,'FontWeight','bold')
      axis on
    


Relerr_spectral_initial= norm(test_image- exp(-1i*angle(trace(test_image.*initialization_spectral))) * initialization_spectral, 'fro')/sqrt(square_norm_x);
fprintf('relative error of the initial guess via the spectral method =   %d; \n',Relerr_spectral_initial);
global_phase=exp(1i*angle(conj(trace(test_image'*initialization_spectral))));
   if (image_type<3.5)
    figure;
    imshow(uint8((initialization_spectral(z+1:z+N,z+1:z+N)))); 
    title ('initial guess via the spectral method');
   elseif (image_type==4)   
   figure;
   imshow(uint8(abs(real(global_phase*initialization_spectral(z+1:z+N,z+1:z+N)))));
   title ('the real part of initial guess via the spectral method');
   figure;
   imshow(uint8(abs(imag(global_phase*initialization_spectral(z+1:z+N,z+1:z+N)))));
   title ('the imaginary part of initial guess via the spectral method');
   elseif (image_type==5)
   figure;
   imshow(uint8(abs(initialization_spectral(z+1:z+N,z+1:z+N))));
   title ('the absolute value of initial guess via the spectral method');
   end

   
end


%% truncated spectral method
if (truncated_spectral_method==1)
alpha_square=5;
initialization_Tspectral(:,:)=random_initialization(:,:);
truncation_threshold=alpha_square*square_norm_x/(N+2*z)/(N+2*z); 


   T0=zeros(N+2*z,N+2*z,number_illumination);
   for mask_index=1:number_illumination
   T0(:,:,mask_index)=(data_box(:,:,mask_index).^2/(N+2*z)/(N+2*z)<=truncation_threshold);
   end
    
    relative_residual_Tspectral_initial=zeros(1,max_iterate_powermethod);
    for tt = 1 : max_iterate_powermethod 
    temp1=zeros(N+2*z);
      for mask_index=1:number_illumination
      temp2 =T0(:,:,mask_index).*(data_box(:,:,mask_index).^2).*fft2(mask_box(:,:,mask_index).*initialization_Tspectral(:,:)); 
      IFFT2_temp2=ifft2(temp2);
      temp1=temp1+IFFT2_temp2(:,:)./mask_box(:,:,mask_index)/number_illumination;
      end
    
      if (image_type==1 || image_type==2 || image_type==3)
      temp1=max(real(temp1),0);
      end
    
     initialization_Tspectral=sqrt(square_norm_x)*temp1.*filter/norm(temp1.*filter,'fro');
    
      %%%%%%%%%%% compute the relative residual %%%%%%%%%%
       residual=0;
       for mask_index=1:number_illumination
       mask_initialization_Tspectral= mask_box(:,:,mask_index).*initialization_Tspectral;  
       FFT_mask_initialization_Tspectral=fft2(mask_initialization_Tspectral);
       abs_FFT_mask_initialization_Tspectral=abs(FFT_mask_initialization_Tspectral);
       residual=residual+norm(data_box(:,:,mask_index)-abs_FFT_mask_initialization_Tspectral(:,:),'fro')^2/norm(data_box(:,:,mask_index),'fro')^2; 
       end
       relative_residual_Tspectral_initial(tt)=sqrt(residual/number_illumination);
       %fprintf('relative residual of the initial guess via the Tspectral method  %d=   %d; \n', tt, relative_residual_Tspectral_initial(tt))
    end 
      figure;
      plot(1:max_iterate_powermethod,relative_residual_Tspectral_initial,'r')
      xlabel('iteration','FontSize',12,'FontWeight','bold')
      ylabel('relative residual for the truncated spectral method','FontSize',12,'FontWeight','bold')
      axis on


Relerr_Tspectral_initial= norm(test_image- exp(-1i*angle(trace(test_image'*initialization_Tspectral))) * initialization_Tspectral, 'fro')/sqrt(square_norm_x);
fprintf('relative error of the initial guess via the truncated spectral method =   %d; \n',Relerr_Tspectral_initial);
    
   if (image_type==1 || image_type==2 || image_type==3)
    figure;
    imshow(uint8((initialization_Tspectral(z+1:z+N,z+1:z+N)))); 
    title ('initial guess via the truncated spectral method');
   elseif (image_type==4)   
   figure;
   imshow(uint8(abs(real(exp(-1i*angle(trace(test_image'*initialization_Tspectral))) *initialization_Tspectral(z+1:z+N,z+1:z+N)))));
   title ('the real part of initial guess via the truncated spectral method');
   figure;
   imshow(uint8(abs(imag(exp(-1i*angle(trace(test_image'*initialization_Tspectral))) *initialization_Tspectral(z+1:z+N,z+1:z+N)))));
   title ('the imaginary part of initial guess via the truncated spectral method');
   elseif (image_type==5)
   figure;
   imshow(uint8(abs(initialization_Tspectral(z+1:z+N,z+1:z+N))));
   title ('the absolute value of initial guess via the truncated spectral method');
   end

end









%% null vector method (the median parameter setup)
if (null_vector_method==1)
partition=100;
%record_of_ratio=0;
record_of_residual=500;
record_of_initial=zeros(N+2*z);
   for ii=50:50
   ratio_null=ii/partition;

   sort_data=sort(data_box(:));
   medium_data=sort_data(round((N+2*z)^2*number_illumination*ratio_null));   
   frequency_choose=zeros(N+2*z,N+2*z,number_illumination);
      for mask_index=1:number_illumination
      frequency_choose(:,:,mask_index)=data_box(:,:,mask_index)>medium_data;      
      end

   initialization_null(:,:)=random_initialization(:,:);
   relative_residual_nullvector=zeros(1,max_iterate_powermethod);
   relative_residual=zeros(1,max_iterate_powermethod);
      for tt=1:max_iterate_powermethod
      temp2_null=zeros(N+2*z);
        for mask_index=1:number_illumination
        temp1_null =frequency_choose(:,:,mask_index).*fft2(mask_box(:,:,mask_index).*initialization_null(:,:)); 
        temp2_null=temp2_null+conj(mask_box(:,:,mask_index)).*ifft2(temp1_null);
        end 

        if (image_type==1 || image_type==2 || image_type==3)
        temp2_null=max(real(temp2_null),0);
        end

        initialization_null= sqrt(square_norm_x)*(temp2_null.*filter)/norm(temp2_null.*filter,'fro');
        
       %%%%%%%%%%% compute the relative residual %%%%%%%%%%
       residual=0;
       for mask_index=1:number_illumination
       mask_nullvector= mask_box(:,:,mask_index).*initialization_null;  
       FFT_mask_nullvector=fft2(mask_nullvector);
       abs_FFT_mask_nullvector=abs(FFT_mask_nullvector);
       residual=residual+norm(data_box(:,:,mask_index)-abs_FFT_mask_nullvector(:,:),'fro')^2/norm(data_box(:,:,mask_index),'fro')^2; 
       end
       relative_residual_nullvector(tt)=sqrt(residual/number_illumination);
      end
   
      figure;
      plot(1:max_iterate_powermethod,relative_residual_nullvector,'r')
      xlabel('iteration','FontSize',12,'FontWeight','bold')
      ylabel('relative residual for the null vector method','FontSize',12,'FontWeight','bold')
      axis on

      residual=0;
      for mask_index=1:number_illumination
      mask_initial= mask_box(:,:,mask_index).*initialization_null;  
      FFT_mask_initial=fft2(mask_initial);
      abs_FFT_mask_initial=abs(FFT_mask_initial);
      
      residual=residual+sum(sum((data_box(:,:,mask_index)-abs_FFT_mask_initial(:,:)).^2))/(norm(data_box(:,:,mask_index),'fro')^2); % error calculation
      end
      relative_residual(ii)=sqrt(residual)/number_illumination;

      if (relative_residual(ii)<record_of_residual)
      record_of_ratio=ratio_null;
      record_of_residual=relative_residual(ii);
      record_of_initial=initialization_null;
      end
    
   end
   Ans_of_Gamma=1-record_of_ratio;
   fprintf('record_of_gamma =   %d; \n',Ans_of_Gamma);
   initialization_null=record_of_initial;
   %fprintf('record_of_residual =   %d; \n',record_of_residual);
   Relerr_null_vector= norm(test_image- exp(1i*angle(conj(trace(test_image'*initialization_null)))) * initialization_null, 'fro')/sqrt(square_norm_x);
   fprintf('relative error of the initial guess via the null vector method =   %d; \n',Relerr_null_vector);
   global_phase=exp(1i*angle(conj(trace(test_image'*initialization_null))));
   if (image_type<3.5)
    figure;
    imshow(uint8((global_phase*initialization_null(z+1:z+N,z+1:z+N)))); 
    title ('initial guess via the null vector method');
   elseif (image_type==4)   
   figure;
   imshow(uint8(abs(real(global_phase*initialization_null(z+1:z+N,z+1:z+N)))));
   title ('the real part of initial guess via the null vector method');
   figure;
   imshow(uint8(abs(imag(global_phase*initialization_null(z+1:z+N,z+1:z+N)))));
   title ('the imaginary part of initial guess via the null vector method');
   elseif (image_type==5)
   figure;
   imshow(uint8(abs(initialization_null(z+1:z+N,z+1:z+N))));
   title ('the absolute value of initial guess via the null vector method');
   end
end



%% the selection of initial input for SAP, PAP and WF below. 
if (initial_selection==3) 
input_SAP=initialization_null;
input_PAP=initialization_null;
input_WF=initialization_null;
elseif (initial_selection==2)
input_SAP=initialization_Tspectral;
input_PAP=initialization_Tspectral;
input_WF=initialization_Tspectral;
elseif (initial_selection==1)
input_SAP=initialization_spectral;
input_PAP=initialization_spectral;
input_WF=initialization_spectral;
end


% SPA, PAP, and WF algorithms:
%%% SAP %%%
if (serialAP==1)
relative_error_SAP=zeros(1,max_iterate+1);
relative_residual_SAP=zeros(1,max_iterate+1);  
guess=input_SAP;
Relerr_guess= norm(test_image- exp(-1i*angle(trace(test_image'*guess))) * guess, 'fro')/sqrt(square_norm_x);
relative_error_SAP(1)=Relerr_guess;

residual_vector=0;
   for index=1:number_illumination
   mask_guess= mask_box(:,:,index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);
   residual_vector=residual_vector+norm(data_box(:,:,index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,index),'fro')^2; % error calculation
   end
relative_residual_SAP(1)=sqrt(residual_vector/number_illumination);
fprintf('Initial input of SAP: (relative residual, relative error)=   (%d,%d); \n',relative_residual_SAP(1),relative_error_SAP(1))
 


for iterate=1:max_iterate
  for mask_index=1:number_illumination    
  mask(:,:)=mask_box(:,:,mask_index); 
  data(:,:)=data_box(:,:,mask_index);
 
  mask_guess=mask.*guess; 
  FFT2_mask_guess=fft2(mask_guess); 
  fourier_phase(:,:)=exp(1i*angle(FFT2_mask_guess(:,:)));
  data_fitting=data.*fourier_phase;
  IFFT_data_fitting=ifft2(data_fitting); 
  guess=IFFT_data_fitting./mask; 
 end

   if (image_type==1 || image_type==2 || image_type==3)
   guess=max(real(guess),0).*filter;
   else 
   guess=guess.*filter;
   end
   
%%%%%%%%%%% compute the relative residual %%%%%%%%%%
   residual=0;
   for mask_index=1:number_illumination
   mask_guess= mask_box(:,:,mask_index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);    
   residual=residual+norm(data_box(:,:,mask_index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,mask_index),'fro')^2; 
   end
   relative_residual_SAP(iterate+1)=sqrt(residual/number_illumination);


Relerr_guess= norm(test_image- exp(-1i*angle(trace(test_image'*guess))) * guess, 'fro')/sqrt(square_norm_x);
relative_error_SAP(iterate+1)=Relerr_guess;  
   
% global phase adjustment
guess=guess*exp(-1i*angle(trace(test_image'*guess)));
end

fprintf('Final output of SAP: (relative residual,relative error)  = (%d,%d); \n',relative_residual_SAP(iterate+1),relative_error_SAP(max_iterate+1));
end




%%% PAP %%%
if (parallelAP==1)
relative_error_PAP=zeros(1,max_iterate+1);
relative_residual_PAP=zeros(1,max_iterate+1);  
guess=input_PAP;
Relerr_guess= norm(test_image- exp(-1i*angle(trace(test_image'*guess))) * guess, 'fro')/sqrt(square_norm_x);
relative_error_PAP(1)=Relerr_guess;

   residual=0;
   for index=1:number_illumination
   mask_guess= mask_box(:,:,index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);   
   residual=residual+norm(data_box(:,:,index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,index),'fro')^2; % error calculation
   end
relative_residual_PAP(1)=sqrt(residual/number_illumination);

fprintf('Initial input of PAP: (relative residual, relative error)=   (%d,%d); \n',relative_residual_PAP(1),relative_error_PAP(1))
 
for iterate=1:max_iterate
 AA=zeros(N+2*z); 
      for index=1:number_illumination    
      mask(:,:)=mask_box(:,:,index); 
      data(:,:)=data_box(:,:,index);
      mask_guess=mask.*guess;  
      FFT2_mask_guess=fft2(mask_guess); 
      fourier_phase=exp(1i*angle(FFT2_mask_guess));
      data_fitting=data.*fourier_phase;
      IFFT_data_fitting=ifft2(data_fitting);  
      AA=AA+IFFT_data_fitting./mask/number_illumination;
      end
 guess=AA;
 %-------------
   if (image_type==1 || image_type==2 || image_type==3)
       guess=max(real(guess),0).*filter;     
   else
       guess=guess.*filter; 
   end
  %-------------------------  
   
%%%%%%%%%%% compute the relative residual %%%%%%%%%%
   residual=0;
   for index=1:number_illumination
   mask_guess= mask_box(:,:,index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);  
   residual=residual+norm(data_box(:,:,index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,index),'fro')^2; % error calculation
   end
   relative_residual_PAP(iterate+1)=sqrt(residual/number_illumination);

   Relerr_guess= norm(test_image- exp(-1i*angle(trace(test_image'*guess))) * guess, 'fro')/sqrt(square_norm_x);
   relative_error_PAP(iterate+1)=Relerr_guess;

guess=exp(-1i*angle(trace(test_image'*guess)))*guess;
end
fprintf('Final output of PAP: (relative residual,relative error) = (%d,%d); \n',relative_residual_PAP(iterate+1),relative_error_PAP(max_iterate+1));
end


%%% WF %%%
if (WF==1)
relative_error_WF=zeros(1,max_iterate+1);
relative_residual_WF=zeros(1,max_iterate+1);  
guess=input_WF;
Relerr_guess= norm(test_image- exp(-1i*angle(trace(test_image'*guess))) * guess, 'fro')/sqrt(square_norm_x);
relative_error_WF(1)=Relerr_guess;

   residual=0;
   for index=1:number_illumination
   mask_guess= mask_box(:,:,index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);   
   residual=residual+norm(data_box(:,:,index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,index),'fro')^2; % error calculation
   end
relative_residual_WF(1)=sqrt(residual/number_illumination);
fprintf('Initial input of WF:(relative residual, relative error)=   (%d,%d); \n',relative_residual_WF(1),relative_error_WF(1))
 


for gradientmethod=1:max_iterate
    AA=zeros(N+2*z);
    for  mask_index=1:number_illumination
        cc=fft2(mask_box(:,:,mask_index).*guess(:,:));  
        dd(:,:)=(abs(cc(:,:)).^2-data_box(:,:,mask_index).^2).*cc(:,:);
        conjmask(:,:)=conj(mask_box(:,:,mask_index));
        AA=AA+ifft2(dd).*conjmask;
    end
    alpha = - log(0.997) ; 
    mu = 1-exp(-alpha*gradientmethod) ; 
    parameter_mu = min(mu,0.4);
    guess=guess-parameter_mu*AA/number_illumination/(square_norm_x);
    
    if (image_type==1 || image_type==2 || image_type==3)
    guess=max(real(guess),0).*filter;
    else
    guess=guess.*filter;
    end
    global_phase=exp(1i*angle(conj(trace(test_image'*guess))));
    relative_error_WF(gradientmethod+1)=(sum(sum(abs(test_image(:,:)-global_phase*guess(:,:)).^2))/square_norm_x)^0.5;
%%%%%%%%%%% prepare to compute the residual error %%%%%%%%%%
   residual=0;
   for index=1:number_illumination
   mask_guess= mask_box(:,:,index).*guess;  
   FFT_mask_guess=fft2(mask_guess);
   abs_FFT_mask_guess=abs(FFT_mask_guess);    
   residual=residual+norm(data_box(:,:,index)-abs_FFT_mask_guess(:,:),'fro')^2/norm(data_box(:,:,index),'fro')^2; % error calculation
   end
   relative_residual_WF(gradientmethod+1)=sqrt(residual/number_illumination);
   guess=exp(-1i*angle(trace(test_image'*guess)))*guess;
end

fprintf('Final output of WF: (relative residual,relative error) = (%d,%d); \n',relative_residual_WF(gradientmethod+1),relative_error_WF(gradientmethod+1));
end






figure;
plot(1:max_iterate+1,relative_residual_SAP,'r',1:max_iterate+1,relative_residual_PAP,'g',1:max_iterate+1,relative_residual_WF,'b')
xlabel('iteration','FontSize',12,'FontWeight','bold')
ylabel('relative residual (R: SAP; G: PAP; B:WF)','FontSize',12,'FontWeight','bold')
axis on




figure;
plot(1:max_iterate+1,relative_error_SAP,'r',1:max_iterate+1,relative_error_PAP,'g',1:max_iterate+1,relative_error_WF,'b')
xlabel('iteration','FontSize',12,'FontWeight','bold')
ylabel('relative error (R: SAP; G: PAP; B:WF)','FontSize',12,'FontWeight','bold')
axis on










end

