% Usage: 
% E = ItrLLSimpute(M)
% 
% M = the missing matrix where missing values are represented by the value 1e99
% E = the imputed dataset (matrix) with all missing values estimated

function [E] = ItrLLSimpute(M)

global m n miss_matrix gene_include m_gene_include oldE big all_k feasible write_fid k_range;

fid = 1;
fprintf(fid,'\n\n');
fprintf(fid,'----------------------------------------------\n');
fprintf(fid, 'Iterative LLSimpute/L2 with ratio training\n');
fprintf(fid,'----------------------------------------------\n');
fprintf(fid,'\n\n');

big=1e99;

miss_matrix = M; % struct array that contains miss_matrix
E = [];
k_range = [inf -inf];
all_k = [];

[m n] = size(miss_matrix);
minexp = 2;

gene_include=[]; oldE=[];

% determine which genes to include in the similar-gene pool
if (m - length(find(miss_matrix==big))) < 400 
    [oldE]=impute_rowavg(miss_matrix,minexp); 
   % consider all genes as potential similar genes
    gene_include=1:m;
    m_gene_include=m;
else    
    oldE=miss_matrix;
    [miss_gene,miss_exp]=find(miss_matrix == big);      % miss_gene is the set of row numbers to each missing value
                                                        % miss_exp is the set of column numbers to each missing value
    % consider only complete genes as potential similar genes
    gene_include = setdiff(1:m, miss_gene);
    m_gene_include = length(gene_include);
    % set minexp to 0 since we do not perform rowaverage
    minexp=0;
end

feasible = 0;
%############################################
%####### ratio estimate #####################
%############################################
min_ratio = 0.1; max_ratio = 1.0; ratio_jump = 0.1;

rho=length(find(miss_matrix==big))/m;
e=floor(rho/5);
% at least 1 position should be pretended as a missing value.
e=max([1 e]);
% the number of pretended missing values should be less than 
% 20% of non-missing positions.....
e=min([e n*0.2*(1-0.2)]);
% pretended missing positions among the non-missing positions.
p=[1:2:2*e];

% check entire row
answer=[];
ratio_count=0;
ratio = min_ratio;
while(ratio <= max_ratio)
    ratio_count=ratio_count+1; guess{ratio_count}=[];  ratio = ratio + ratio_jump;
end

for i=1:m
    missidxj=find(miss_matrix(i,:) == big); len_miss=length(missidxj);
    nomissidxj=find(miss_matrix(i,:) < big);  len_nomiss=length(nomissidxj);
    if (len_nomiss < minexp | len_nomiss < 2) 
    elseif (len_miss > 0)  
         
        % determine the number of artifical missing values among p
        missnum=min([len_miss length(p)]); 
        % too small nomissing entries to make artificial missing entries
        if (len_nomiss < 10)
            missnum=1;
        end
        p1=p(1:missnum);

        % generate the artificial missing values
        missidxj=nomissidxj(p1);
        nomissidxj(p1)=[];         
        answer=[answer; oldE(i,missidxj)'];
        
        % for ratio survey
        ratio_count=0;
        ratio = min_ratio;
        while(ratio <= max_ratio)
            ratio_count=ratio_count+1;
            % gather the nearest neighbors (similar genes)
            [A,B,w] = closegene(i,missidxj,nomissidxj, ratio); 
            A=A';
            %linear combination of experiments (performing least squares
            %optimization)
            X = pinv(A)*w';  
                        
            guess{ratio_count} = [guess{ratio_count}; B'*X];
            ratio = ratio + ratio_jump;
        end%while k
    end%if
    
end%i

ratio = min_ratio; ratio_count=0; xratio=[];
while(ratio <= max_ratio)
    ratio_count=ratio_count+1;
    nrmse(ratio_count) = sqrt( mean( (guess{ratio_count}-answer).^2 )) / std( answer );
    xratio=[xratio ratio];  ratio = ratio + ratio_jump;
end

% try to find local minimum
[minnrmse,min_ratioidx]=min(nrmse); ratio=xratio(min_ratioidx);

%############################################
%######### ratio estimate   #################
%############################################

fprintf(fid, 'Feasible ratio for choosing neighbors: %f \n\n', ratio);
gene_include = 1:m; m_gene_include = m;
[oldE]=impute_rowavg(miss_matrix,minexp); 
feasible = 1;
% once the feasible ratio is determined, continue to the imputation stage
E = inner_impute(ratio);
k_range
return;

% This function imputes the missing values of the given set using the given ratio
function [newE] = inner_impute(ratio)

global m n miss_matrix big oldE newE gene_include m_gene_include all_itr_k feasible write_fid;

all_itr_k = [];
fid = 1;
total=m*n; 

% the number of minimal experiments for estimating missing values by llsq
minexp=2;
 
newE = miss_matrix;  

fprintf(fid, 'miss_matrix(%d,%d) total: %d minexp: %.2f \n\n', ...
    m, n, total, minexp);
fprintf(fid, 'Estimating missing values...\n');  

% we don't know the exact answer for the missing values  
guess=[];

%############### iterative estimation method #########################

itrs = 0;
while 1
    all_itr_k = [];
    guess = [];
    for i=1:m
        
        missidxj=find(miss_matrix(i,:) == big); len_miss=length(missidxj);
        nomissidxj=find(miss_matrix(i,:) < big);  len_nomiss=length(nomissidxj);
        if (len_nomiss < minexp) | (len_nomiss < 2)
            %fprintf('%dth gene: skip due to nomiss_exp(%d)<%.2f or < 2\n', i, len_nomiss,minexp);      
        elseif (len_miss > 0)  
            % find genes that are closer to gene i determined by the threshold
            [A,B,w] = closegene(i,missidxj,nomissidxj, ratio);
            
            A=A';
            %linear combination of experiments
            X = pinv(A)*w';  
            guess = [guess; B'*X];
        end       
    end
    %store estimated values
    s=1;  
    for i=1:m
        missidxj=find(miss_matrix(i,:) == big); len_miss=length(missidxj);  
        nomissidxj=find(miss_matrix(i,:) < big);   len_nomiss=length(nomissidxj);
        if (len_nomiss < minexp) | (len_nomiss < 2)
            % skip
        elseif (len_miss > 0)  
            newE(i,missidxj) = guess(s:s+len_miss-1)'; s=s+len_miss;
        end
    end
    
    itrs = itrs + 1;
    
    if itrs == 5
        fprintf('\nTotal number of iterations: %d\n', itrs);
        return;
    else        
        oldE = newE;
    end    
end
%#######################################################################################
return;


% In this function, the value k is determined before returning from the function
function [A,B,w] = closegene(i, missidxj, nomissidxj, ratio)
global m n miss_matrix oldE gene_include m_gene_include all_k all_itr_k feasible big k_range;
% L2-norm distance calculation

mm1=1;  mm2=m_gene_include;
AA1=oldE(i,nomissidxj); BB1=oldE(gene_include,nomissidxj);   
AA2=sum(AA1.^2,2); BB2=sum(BB1.^2,2);
distance=repmat(AA2,1,mm2)+repmat(BB2',mm1,1)-2*AA1*BB1';

if any(isnan(distance))
   [find(isnan(distance)) i]
end

[distsorted,sortidx]=sort(distance); 
% gene number
gene=gene_include(sortidx);    % gene is the sorted set (based on distance) of gene_include
% remove the same gene
gene(find(gene==i))=[];
% here mink is the first few genes that satisfy the min distance to gene i
mink = min(length(find(distsorted <= (mean(distsorted)*ratio))), length(gene));
mink = max(mink, 1);

if (feasible)
    k_range = [min(k_range(1), mink) max(k_range(2), mink)];
end

% here A, B and w are according to the LLSimpute paper
A=oldE(gene(1:mink),nomissidxj);       
B=oldE(gene(1:mink),missidxj);         
w=miss_matrix(i,nomissidxj);         

return;

function [E] = impute_rowavg(miss_set,minexp)
global m n big;

% return the matrix that has the missing values estimated using row average

E=miss_set; gene0=[]; gene1=[];
for i=1:m
      missidxj=find(miss_set(i,:) == big);
      if (length(missidxj) > 0)
         nomissidxj=find(miss_set(i,:) < big);
         exp=length(nomissidxj);
         if(exp==0) 
             gene0=[gene0 i];
         else
             % row average
            avg=mean(miss_set(i,nomissidxj)); 
            E(i,missidxj)=avg;
            if(exp < minexp)
                gene1=[gene1 i];
            end
        end   
     end
end

return;


