function [X]=TFBSImpute(D,W,res,max_iter,p,u0)
%D=randn(r,c,num_position);% TF binding  tensor, r is the number of TFs, c is the number of cells;
%W=zeros(r,c,num_position);% weight matrix, the same size with D; 1 represents observed, 0 represents missing;

% parameter settings
%p% 1<p<2;
%u0;%u>0 
u=u0;
%res=0.05*sum(sum(sum(W.*D)))/sum(sum(sum(W)));% s.t. |Xi(k,h)-Di(k,h)|<=res, for all observed data.
a=1; % regularization papameter 

A=zeros(size(W));
X=zeros(size(W));
Z=zeros(size(W)); 

[~,~,num_position]=size(W);
% get the index of observed samples in D(:,:,i);
Wi=W(:,:,1);num_o_i=sum(sum(Wi));index_o=zeros(num_o_i,2);[num_x,num_y]=size(Wi);k=1;
for i=1:num_x
    for j=1:num_y
        if Wi(i,j)==1
            index_o(k,1)=i; index_o(k,2)=j;
            k=k+1;
        end
    end
end
temp_res =sum(sum(sum(abs(W.*(X-D)))));

% main 
for iter=1:max_iter% converge
    % step 1: update Z by solving problem (14);
    Xi=X(:,:,1);
    Ai=A(:,:,1); 
    for i=1:num_position
        Ai=A(:,:,i);Xi=X(:,:,i);
        [U,S,V]=svd(Xi-(1/u)*(Ai),'econ');%Mi=Xi-(1/u)*(Ai);[U,S,V]=svd(Mi,'econ');
        S=diag(max(0,diag(S)-1/u));%soft thresholding SVD
        Z(:,:,i)=U*S*V'; %Zi=U*S*V'; 
    end
    
    %step 2: update X
    % update N
    Q=zeros(size(Xi));
    Di=D(:,:,1);
    N=Z+(1/u)*A;
    for i=1:num_position
        Di=D(:,:,i);%??
        if i==1
           Q=(2/(2*a+u))*(a*X(:,:,2)+(u/2)*N(:,:,1));
        else if i==num_position
            Q=(2/(2*a+u))*(a*X(:,:,num_position-1)+(u/2)*N(:,:,num_position));
            else
                 Q=2/(4*a+u)*(a*X(:,:,i+1)+a*X(:,:,i-1)+(u/2)*N(:,:,i));
            end
        end
        
        Wi=W(:,:,i);
        w1=Wi & ((Q-Di)>res);
        w2=Wi & ((Q-Di)<-res);
        Xi=w1.*(Di+res)+w2.*(Di-res);
        w3=1-(w1 | w2);
        Xi=Xi+w3.*Q;

        X(:,:,i)=Xi;
   end
   
    %step 3: updata Ai,u
    A=A+u*(Z-X);
    u=p*u;
    
    if temp_res<res
         break;
    end
end