clear all;

rng(1)

tic

%% test image

image_type=3;
n = 144; % n (number of pixels) can not be too large. Otherwise, your memory will not be enough.
% image_type 1 : white noise Type 1; choose n s.t. mod(n,8)=0
% image_type 2 : White noise Type 2; choose n s.t. mod(n,2)=0
% image_type 3 : Phantom with random phases; choose n s.t. n=N^2 for some integer N


%% the number of random illumination & oversampling ratio of Fourier intensity data


    if (image_type==1)  
    matrix1=0:n-1;
    matrix1=matrix1'/n;
    matrix2=(-n/8)+1:(n/8);
    matrix3=exp(2*pi*1i*matrix1*matrix2);
    matrix4=normrnd(0,1,[n/4,1])+1i*normrnd(0,1,[n/4,1]);
    test_image=matrix3*matrix4;
    elseif (image_type==2)
    matrix1=0:n-1;
    matrix1=matrix1'/n;
    matrix2=(-n/2+1):(n/2);
    matrix3=exp(2*pi*1i*matrix1*matrix2);
    matrix4=(1/sqrt(8))*normrnd(0,1,[n,1])+1i*(1/sqrt(8))*normrnd(0,1,[n,1]);
    test_image=matrix3*matrix4;
    elseif (image_type==3)
    N=ceil(sqrt(n));
    test_image=256*exp(1i*2*pi*rand(N,N)).*phantom(N); 
    test_image=test_image(:);
    test_image=test_image(1:n);
    end

%% normalize test image
norm_x=norm(test_image,'fro');
test_image=test_image/norm(test_image,'fro');



%% plot the reshaped test image RPP
if (image_type==3)
reshaped_test_image = norm_x*reshape(test_image,sqrt(n),sqrt(n));
figure;
imshow(uint8(abs(reshaped_test_image)));
title ('absolute value of the reshaped test image');
end



%% the oversampling ratio

min_ratio = 25;  % minimum oversampling ratio
max_ratio = 500; % maximum oversampling ratio
step_size_ratio = 25; % the step size for increasing the oversampling ratio from min_ratio to max_ratio


number_step=floor(((max_ratio-min_ratio)/step_size_ratio)+1);
max_ratio=min_ratio+(number_step-1)*step_size_ratio;


% parameter setup for the threshold used in the null vector method
a1=1;
a2=2;
alpha=a1/a2; % or 3/4;
dynamic_alpha=1; % dynamical threshold for the null vector method



operatornorm_error_null = zeros(number_step,1); % recording error
frobeniusnorm_error_null = zeros(number_step,1);

hwait=waitbar(0,'Please wait >>>>>>>>');

    

sensing_matrix=(1/sqrt(2))*normrnd(0,1,[n min_ratio*n])+(1/sqrt(2))*1i*normrnd(0,1,[n min_ratio*n]); % initial sensing matrix
data_box=abs(sensing_matrix'*test_image); % diffraction pattern
new=[];
    

for k=1:number_step
        ratio_temp=k/number_step;    
        PerStr=(100*ratio_temp); 
        str=['running ...',num2str(PerStr),'%'];
        waitbar(ratio_temp,hwait,str);
        
        
sensing_matrix=[sensing_matrix new]; %extending the size of the last sensing matrix
 
if (k>1) % New data will be generated when k>1. 
new_data_box=abs(new'*test_image);
data_box=[data_box;new_data_box];
end

new=(1/sqrt(2))*normrnd(0,1,[n step_size_ratio*n])+(1/sqrt(2))*1i*normrnd(0,1,[n step_size_ratio*n]); % new sensing matrix

    
   if (dynamic_alpha==0) % setting the threshold for the data
   threshold_data=quantile(data_box, alpha); 
   else
       value_quantile= ((min_ratio+(k-1)*step_size_ratio))^(alpha-1); %sampling_ratio
       threshold_data=quantile(data_box,value_quantile);
   end
   
       sub_column_index=find(data_box>threshold_data);
       [Q,R] = qr(sensing_matrix',0);
       inverse_R=R^(-1);
       clear R
       Q_complement=Q(sub_column_index,:);
       clear Q
       reverse=Q_complement'*Q_complement;
       clear Q_complement
       [initialization_null,EB]=eigs(reverse,1);
       clear reverse
       initialization_null=inverse_R*initialization_null;
       clear inverse_R
       initialization_null= initialization_null/norm(initialization_null,'fro');

      %%%%%%%%%%% compute paper norm %%%%%%%%%%%
       frobeniusnorm_error_null(k,1)= norm(test_image- exp(-1i*angle(trace(test_image'*initialization_null)))*initialization_null, 'fro');
       aa = frobeniusnorm_error_null(k,1);
       operatornorm_error_null(k,1)=((1-(1-aa^2/2)^2)*2)^0.5;
       %Relerr_null(k,seed)= norm(test_image- exp(-1i*angle(trace(test_image'*initialization_null)))*initialization_null, 'fro');
       %fprintf('relative error (Frobenius Norm) =   %d; \n',frobeniusnorm_error_null(k,1));
       name=['L = ',num2str(min_ratio+(k-1)*step_size_ratio), ' relative error (Frobenius Norm) =' num2str(frobeniusnorm_error_null(k,1)),'\n'];
       fprintf(name);
   end

close(hwait);


toc     
     

      figure;      
      plot(min_ratio:step_size_ratio:max_ratio,frobeniusnorm_error_null(1:number_step),'r'); 
      hold on
      plot(min_ratio:step_size_ratio:max_ratio,operatornorm_error_null(1:number_step),'b'); 
      hold off
      xlabel('oversampling ratio (L)','FontSize',12,'FontWeight','bold')
      ylabel('error','FontSize',12,'FontWeight','bold')
      axis on
      legend('frobenius norm','operator norm')