function [u,u_best] = GPS_algorithm(input)
diffpat = input.diffpat;
support = input.support; 
ntimes  = input.iterations; 
isreal  = input.isreal;

[n1, n2]      = size(diffpat)                ;

x_center = ceil(n1/2); y_center = ceil(n2/2);

index0 = support == 0;



% filter parameter
X = 1:ntimes; filtercount = 10;
sub_iter = ntimes/filtercount;
FX=(filtercount + 1-ceil(X*filtercount/ntimes))*ceil(ntimes/(1*filtercount));
FX=((FX-ceil(ntimes/filtercount))*(2*n1)/max(FX))+(2*n1/10);
[KK,JJ] = meshgrid(1:n1,1:n2);
R2 = (KK - x_center).^2 + (JJ - y_center).^2 ;

%% OSS
Z_initial = diffpat;
stopper = find(diffpat==-1);
unknown = diffpat==-1; known = ~unknown;
Z_initial (stopper) = 0;

% initialization
if isfield(input,'sigma')
    sigma = input.sigma;
else
    sigma=0.01;
end

if isfield(input, 'initial')
    z = fftshift(fft2(input.initial));
    z = z /norm(z,'fro')*norm(diffpat,'fro');
else
    b = fftshift(fft(support.*randn(size(support)) ));
    %z  = exp(1i*phase_angle) .* Z_initial;    
    z  = exp(1i*angle(b)) .* Z_initial;
end

fprintf('Done with initialization, starting loop\n')

y = zeros(n1,n2);
dt=1; ds=1;
errorF_best=1;
u_best = zeros(size(z));
%% Loop
for t = 1:ntimes
    z_old = z;
    Fy = fftshift(fft2(y));
    d = z - dt * Fy;
    d_stopper = d(stopper);
    phase = angle(d);
    if t>0.8*ntimes,  sigma=0.1; end
    z = ( exp(1i*phase) .* diffpat + sigma*d ) / (1+ sigma);
    z(stopper) = d_stopper;
    if t==1 && isfield(input, 'initial') 
        z = z_old;  
    end
    
    % update y
    %if t==0, u = input.initial;
    u = ifft2(ifftshift(z) );    
    u_hat = ifft2(ifftshift(2*z - z_old));
    y = y + ds*u_hat;
    if isreal
        y_hat = real(y); y_hat(y_hat<0 | index0) = 0;    
        y =  y - y_hat;
    else
        y(support)=0;
    end
    
    % smooth out z FX(t)^2)
    if mod(t,sub_iter)==1
        kfilter = exp( -R2 /(2* FX(t)^2 )); kfilter = kfilter/max(max(kfilter));
    end
    % type 2: Poisson noise removal on Fourier space
    %y_temp = y.*kfilter; y(index0) = y_temp(index0);
    %y = y.*kfilter;
    
    % type 1: low pass filter
    %Fy = my_fft(y) .* kfilter; y = my_ifft(Fy);
    Fy = fftshift(fft2(y)) .* kfilter; y = ifft2(ifftshift(Fy));
    
    % type 3: apply type1 and type2 alternatively
    %if mod(t,2)==0, Fy = fftshift(fft2(y)) .* kfilter; y = ifft2(ifftshift(Fy)); else y = y.*kfilter; end
    
    
    if mod(t,sub_iter)==0
        z = Z2D_best;
        y = Y2D_best;
    end
    
    if mod(t,2)==0
        %u = ifft2(ifftshift(z));
        rec = support .* (u);
        z_temp = abs( my_fft(rec));
        errorF = sum(sum(abs(z_temp(known)-diffpat(known)))) / sum(sum(diffpat(known)));
        
        if t>10 && errorF <= errorF_best
            errorF_best = errorF;
            Z2D_best = z;
            Y2D_best = y;
            u_best = u;
        end
                
    end
    
    if mod(t,100)==0
        fprintf('%d. ErrorF = %f\n', t, errorF) ;
        rec_crop = croppedOut(rec,[512,512]);
        figure(101); img( rec_crop, 'mag', 'colormap','gray'); drawnow();
    end
    
end

fprintf('All done!\n')


