Nearest-Neighbour-Search-Neural-Networks-Product-Quantization/MNIST/Product Quantization/PQNNSwithADC.m

function [Performances, normCCost, iXMinDist] = PQNNSwithADC( kk, mm, nRecall, flagVec)

    %clear all
    close all
    clc
    
    addpath('../Functions');
    addpath('./Functions');
    addpath('./Functions/Yael Library');

    % kk=100;                   % Number of Centroids for each set
    % mm=16;                    % Number of Centroid Sets to take into account
    % kCent=kk^mm;               % Total cardinality of search space
    % nRecall=1;               % Number of Recalls
    
    imgSz=784;                  % Image Descriptor size
    subImgSz=imgSz/mm;         % Size of sub images to analyze
    %Nlearn=6e4;                 % Number of pictures learn set
    Ntrain=6e4;                 % Number of pictures train set
    Ntests=1e4;                 % Number of pictures test set
    PQCentCompon=flagVec(1);    % flag to activate Centroids Computation 
    Quantizationon=flagVec(2);  % flag to activate Data Sub Quantization Indices
    ADCCompon=flagVec(3);       % flag to activate ADC Computations
    dispPQCenton=flagVec(4);    % flag to activate Centroids Displayt
    justGetResultson=flagVec(5);% flag to just get results from stored data


    if(~justGetResultson)

    
    %% Computing/Loading PQ SubCentroids 

        ndots=10;
        str2print='';
        ll=10;
        while mm/ll>1
            ndots=ndots-1;
            ll=ll*10;
        end
        ll=10;
        while kk/ll>1
            ndots=ndots-1;
            ll=ll*10;
        end
        for ll=0:ndots
            str2print=strcat(str2print,'.');
        end

                %% Reading Binary Learning Data Set

            fprintf('Loading Learning Data Set ...................... ');
%             TSetMat=fvecs_read('../Data/sift/sift_learn.fvecs');
            fileOffset = 16;    % bytes of file offset
            fileID = fopen('../Data/train-images.idx3-ubyte');
            XSetVec = fread(fileID,imgSz*Ntrain+fileOffset,'uint8');
            XSetVec = XSetVec(fileOffset+1:end);
            XSetMat = zeros(imgSz,Ntrain);
            for jj=1:Ntrain
                XSetMat(:,jj)=XSetVec((jj-1)*imgSz+1:jj*imgSz);
            end
            fclose(fileID);
            fprintf('Done.\n');
            clear XSetVec

        
        if(PQCentCompon)
        %% Preparing Training Data for Computation
    
            fprintf('Computing m=%1.0f, k=%1.0f PQ subcentroids ...%s',mm,kk,str2print);        
            CSubSetMat=zeros(subImgSz,kk,mm);
            for jj=1:mm
                CSubSetMat(:,:,jj)=yael_kmeans(...
                        single(XSetMat((jj-1)*subImgSz+1:jj*subImgSz,:)),...
                        kk, 'niter', 100, 'verbose', 0);
                %[~, CCtemp]=kmeans(TSubSetMat(:,:,jj)',kk);
                %CSubSetMat(:,:,jj)=CCtemp';
                pp=double(jj/mm*100);
                if(pp<10)
                    fprintf('\b\b%2.0f%%',pp);
                else
                    fprintf('\b\b\b%2.0f%%',pp);
                end
            end
            fprintf('\b\b Done.\n');
%             clearvars TSetMat
            
            save(strcat('./Data/k',int2str(kk),'m',int2str(mm),'CSubSetMat.mat'),'CSubSetMat'); 
        else
            fprintf('Loading m=%1.0f, k=%1.0f PQ subcentroids ....%s',mm,kk,str2print);
            load(strcat('./Data/k',int2str(kk),'m',int2str(mm),'CSubSetMat.mat'));
            fprintf(' Done.\n');
        end
    
%         clearvars TSubSetMat
        
    %% Displaying PQ SubCentroidss
        if(dispPQCenton)
            fprintf('Displaying PQ subcentroids: ...    \n',mm);

            mmax=min(3,mm);
    
            for jj=1:mmax
                figure('units','normalized','outerposition',[0 0 1 1])
                suptitle('Centroid Portions of Pictures')
                imshow255_texmex(0,imgSz,[255*ones(kk,(jj-1)*subImgSz) ...
                               CSubSetMat(:,:,jj)' ...
                               255*ones(kk,(mm-jj)*subImgSz)]);
            end
        end

    
    %% Preparing Training Data for Computation

%         fprintf('Loading Training Data Set ...................... ');
%         XSetMat=fvecs_read('../Data/sift/sift_base.fvecs');
%         fprintf('Done.\n');
    
        
    %% Computing/Loading PQ SubCentroids for Training Data

        if(Quantizationon)
            nSteps=100;         % number of steps for computation, be sure it divides 1e6
            fprintf('Computing PQ subcentroids for Training Data ....    ');
            iXQ=zeros(mm,Ntrain);
            for jj=1:mm
                for ll=1:Ntrain/nSteps
                    [~,iXQ(jj,(ll-1)*nSteps+1:ll*nSteps)]=Quantization(...
                        XSetMat((jj-1)*subImgSz+1:jj*subImgSz,(ll-1)*nSteps+1:ll*nSteps),...
                        CSubSetMat(:,:,jj));
                end
                pp=double(jj/mm*100);
                if(pp<10)
                    fprintf('\b\b%2.0f%%',pp);
                else
                    fprintf('\b\b\b%2.0f%%',pp);fprintf('Computing m=%1.0f, k=%1.0f PQ subcentroids ...%s',mm,kk,str2print);        
            CSubSetMat=zeros(subImgSz,kk,mm);
            for jj=1:mm
                CSubSetMat(:,:,jj)=yael_kmeans(...
                        single(XSetMat((jj-1)*subImgSz+1:jj*subImgSz,:)),...
                        kk, 'niter', 100, 'verbose', 0);
                %[~, CCtemp]=kmeans(TSubSetMat(:,:,jj)',kk);
                %CSubSetMat(:,:,jj)=CCtemp';
                pp=double(jj/mm*100);
                if(pp<10)
                    fprintf('\b\b%2.0f%%',pp);
                else
                    fprintf('\b\b\b%2.0f%%',pp);
                end
            end
            fprintf('\b\b Done.\n');
                end
            end
            fprintf('\b\b\b\b\b\bDone.\n');
            save(strcat('./Data/ADC/ADCk',int2str(kk),'m',int2str(mm),'iXQ.mat'),'iXQ'); 
        else
            fprintf('Loading PQ subcentroids for Training Data ...... ');
            load(strcat('./Data/ADC/ADCk',int2str(kk),'m',int2str(mm),'iXQ.mat'));
            fprintf('Done.\n');
        end

        clearvars XSetMat
        
    %% Preparing Testing Data for Computation

        fprintf('Loading Testing Data Set ....................... ');
%         YSetMat=fvecs_read('../Data/sift/sift_query.fvecs');
        fileID = fopen('../Data/t10k-images.idx3-ubyte');
        YSetVec = fread(fileID,imgSz*Ntests+fileOffset,'uint8');
        YSetVec = YSetVec(fileOffset+1:end);
        YSetMat = zeros(imgSz,Ntests);
        for jj=1:Ntests
            YSetMat(:,jj)=YSetVec((jj-1)*imgSz+1:jj*imgSz);
        end
        fclose(fileID);
        fprintf('Done.\n');
        clear YSetVec
        
    %% Computing/Loading NNs by applying ADC

        if(ADCCompon)
            fprintf('Solving NNS Applying ADC ........................      ');
        
            [iXMinDist]= AsymmetricDistanceComputation(CSubSetMat, iXQ, YSetMat, nRecall);
            save(strcat('./Data/ADC/ADCk',int2str(kk),'m',int2str(mm),'iXMinDist.mat'),'iXMinDist');
            fprintf('\b\b\bDone.\n');
        else
            fprintf('Loading NNS results with ADC ...................... ');
            load(strcat('./Data/ADC/ADCk',int2str(kk),'m',int2str(mm),'iXMinDist.mat'));
            fprintf('Done.\n');
        end
    
        clearvars iXQ YSetMat CSubSetMat 
    else
        fprintf('Loading NNS results with ADC ...................... ');
        load(strcat('./Data/ADC/ADCk',int2str(kk),'m',int2str(mm),'iXMinDist.mat'));
        fprintf('Done.\n');
end
    
%% Loading results    

fprintf('Loading Results Data Set ............');
% iSet=ivecs_read('../Data/sift/sift_groundtruth.ivecs')+1;
fileOffsetL = 8;    % bytes of file offset
fileIDL = fopen('../Data/train-labels.idx1-ubyte');
XLab = fread(fileIDL,Ntrain+fileOffsetL,'uint8');
XLab = XLab(fileOffsetL+1:end);
fclose(fileIDL);
fileIDL = fopen('../Data/t10k-labels.idx1-ubyte');
YLab = fread(fileIDL,Ntests+fileOffsetL,'uint8');
YLab = YLab(fileOffsetL+1:end);
fclose(fileIDL);
fprintf('Done.\n');

%iSet=iSet(1,:);
    
%% Testing PQ Data Set

fprintf('Testing PQ Test Data Set ...    \n');

%     reqRecall=zeros(Ntests,1);
%     
%     for ii = 1:Ntests 
%         iXeq=sum(XLab(iXMinDist(ii,:))==YLab(ii));
%         if numel(iXeq)==1
%             reqRecall(ii)=iXeq;
%         else
%             reqRecall(ii)=nRecall+1; 
%         end
%     end
%     
%     reqRecall=sort(reqRecall);
%     
%     ll=1;
      recAtR=[1 2 5 10 20 50 100 200 500 1000 2000 5000 10000];
      recAtR=recAtR(recAtR<=nRecall);
%     % recAtR=nRecall;
%     Performances=zeros(1,length(recAtR));
%     
%     for ii=recAtR
%         if ii <= nRecall
%             Performances(ll) = length (find (reqRecall <= ii & reqRecall <= nRecall)) / Ntests * 100;
%             fprintf ('Recall@%3d = %.3f\n', ii, Performances(ll));
%             ll=ll+1;
%         end
%     end
%    
    Performances=zeros(1,length(recAtR));
    for rr=1:length(recAtR)
        if recAtR(rr)==1
            Performances(rr)=sum(XLab(iXMinDist(:,1))==YLab)/(Ntests);
        else
            Performances(rr)=sum(max((XLab(iXMinDist(:,1:recAtR(rr)))==repmat(YLab,1,recAtR(rr)))'))/(Ntests);
        end
        fprintf ('Recall@%3d = %.3f\n', recAtR(rr), Performances(rr)*100);
    end

%% Computational Cost Evaluation

maxCCost=Ntrain*Ntests*imgSz;
quantCCost=Ntests*kk*imgSz;
PQCCostADC=Ntrain*Ntests*mm+quantCCost;

normCCost=PQCCostADC/maxCCost;
fprintf('Computational Cost: %2.2f %% w.r.t Exhaustive NNS\n',normCCost*100);