function [template,I_syn, LogZ]=sparseFRAMElearnGibbs_multipleSelection(filters,nScaleGabor, nScaleDoG, nOrient, allSymbol, half, rHat,sx,sy,halfFilterSize,locationShiftLimit, nTileRow,nTileCol,lambdaLearningRate,numFilter,numWavelet,interval,Corr,threshold_corrBB, c_val_list, lower_bound_rand, upper_bound_rand,  nPartCol, nPartRow, part_sx, part_sy, gradient_threshold_scale, savingFolder)


Synfolder=[savingFolder 'Synthesis/'];
if exist(Synfolder)
    rmdir(Synfolder,'s'); 
    mkdir(Synfolder);
else
    mkdir(Synfolder);
end

templateFolder=[savingFolder 'sketchTemplate/'];
if exist(templateFolder)
    rmdir(templateFolder,'s'); 
    mkdir(templateFolder);
else
    mkdir(templateFolder);
end


I_syn = single(randn(sx*nTileRow,sy*nTileCol));  % synthesized image
SSD=[];
numGabor=nScaleGabor*nOrient;

initialLogZ =0 ;
logZRatioSeries = [];

% initialize the template
filterSelected=[];
xSelected=[];
ySelected=[];
selectedLambdas=[];

template=struct('selectedFilter',single(filterSelected),'selectedRow',single(xSelected),'selectedCol',single(ySelected),...
    'selectedLambdas',single(selectedLambdas),'numSelected',single(0));

countMap = cell(numFilter,1); % this map separate sin part and cos part, and count them separately.
for iFilter = 1:numFilter
    countMap{iFilter}=zeros(sx,sy);
end
countTotal=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
countMap_sketch = cell(numFilter,1);  % for drawing sketches, this map combine sin part and cos part, and only count them once.
for iFilter = 1:numFilter
    countMap_sketch{iFilter}=zeros(sx,sy);
end
filterList_sketch = [];
rowList_sketch=[];
colList_sketch=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I_squared_sum=zeros(1, nTileCol*nTileRow,'single');
index=0;
for j=1:nTileCol
    for i=1:nTileRow
       index=index+1;
       I_crop=I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy);  
       I_squared_sum(index) = single(sum(sum(I_crop.^2)));  % I.^2
    end
end
rHat_selected=[];
response_IB=[];

%maximumWaveletSelected = ceil((numWavelet * 2) / numWaveletSelectedPerIter) * numWaveletSelectedPerIter; 
%corrBB=zeros(maximumWaveletSelected, maximumWaveletSelected, 'single');   % the corrBB map should be initialized with the capacity of maximum possible number of selected wavelets 
corrBB=zeros(numWavelet*2, numWavelet*2, 'single');   % the corrBB map should be initialized with the capacity of maximum possible number of selected wavelets 
gradientF_selected=[];

numWaveletSelectedPerIter=nPartCol*nPartRow;
nStroke=0;
iter=0; % iteration counter
while (template.numSelected <= numWavelet*2 - numWaveletSelectedPerIter)   % 

     iter=iter+1;
     disp( [ 'iteration: ' num2str(iter)]);      
     if mod(iter-1,interval)==0     % we add one basis every few steps. The interval is set by users.  
         
        toDrawSketchTemplate=0;
        %%% forward selection (multiple selection version)
        %[x_added, y_added, filter_added]=forwardSelection(filters, halfFilterSize, locationShiftLimit, rHat, I_syn, nTileRow, nTileCol, sx, sy);    
        [x_added_list, y_added_list, filter_added_list, gradient_added_list, numAdded]=forwardSelection_multiple(filters, halfFilterSize, locationShiftLimit, ...
            rHat, I_syn, nTileRow, nTileCol, sx, sy, nPartCol, nPartRow, part_sx, part_sy, gradient_threshold_scale);   % multiple selection
       
       disp( [ num2str(numAdded) ' filters are selected!']); 
       for iAdd=1:numAdded
        
        x_added = x_added_list(iAdd);
        y_added = y_added_list(iAdd);
        filter_added = filter_added_list(iAdd);
        gradient_added = gradient_added_list(iAdd);
        
        if countMap{filter_added}(x_added, y_added)==0
                       
            % record the selected wavelet
            countMap{filter_added}(x_added, y_added)=1;
            countTotal = countTotal+1;
            
            % add the information to template
            template.selectedFilter = [template.selectedFilter filter_added];
            template.selectedRow = [template.selectedRow x_added];
            template.selectedCol = [template.selectedCol y_added];
            template.selectedLambdas = [template.selectedLambdas 0];           
            template.numSelected = template.numSelected +1;    
            
            gradientF_selected = [gradientF_selected 0];        
            
            % for drawing sketch template, we create a template combining sin part and cos part together
            if filter_added <= numGabor; %Gabor
                 Scale = floor( ( filter_added-1)/ (2*nOrient)) + 1;
                 Orient = mod(filter_added - 1, nOrient) + 1;
                 filter_added_ForSketch = Orient + (Scale-1) * nOrient;
            else  % DoG
                 filter_added_ForSketch = filter_added - nScaleGabor*nOrient;                
            end
            if countMap_sketch{filter_added_ForSketch}(x_added, y_added)==0
                    countMap_sketch{filter_added_ForSketch}(x_added, y_added)=1;
                    filterList_sketch = [filterList_sketch filter_added_ForSketch];
                    rowList_sketch=[rowList_sketch x_added];
                    colList_sketch=[colList_sketch y_added]; 
                    toDrawSketchTemplate=1;   % at least one new sketch is selected, we then draw sketch template 
                    nStroke= nStroke+1;
            end

            rHat_selected=[rHat_selected rHat{filter_added}(x_added,y_added)];            
            
            mexc_updateCorrelationTable(corrBB, Corr, single(halfFilterSize), template);     % we update the correlation map instead of recomputing it
                    
            nonzero_corrBB_List=cell(1,template.numSelected);  % neighborhood system (i-th cell contains the indices of all bases whose correlations with B_i are not equal to 0.)
            corrBB2=corrBB(1:template.numSelected, 1:template.numSelected); % crop the valid part from corrBB

            for i=1:template.numSelected
                  nonzero_corrBB_List{i}=single(find(corrBB2(i,:)>threshold_corrBB | corrBB2(i,:)<-threshold_corrBB));
            end
            %%%%
            
            %response_IB=zeros(nTileCol*nTileRow,template.numSelected,'single');
            response_added=zeros(nTileCol*nTileRow, 1,'single');
            index=0;
            for j=1:nTileCol
                for i=1:nTileRow
                    index=index+1;
                    I_crop=I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy);                     
                    response_added(index) = sparseFiltering(I_crop, filters, filter_added, x_added, y_added);                          
                end
            end
            response_IB=[response_IB, response_added];
             
            disp(['filter: ' num2str(filter_added), '; x: ' num2str(x_added), '; y: ', num2str(y_added), ' with gradient = ', num2str(gradient_added)]);            
            
        end   
        
       end
       
       if toDrawSketchTemplate   % draw sketch template only if we have new sketch added
             sketchTemplate = zeros(sx, sy); 
             drawSketch( nOrient * nScaleGabor+nScaleDoG, half, sx, sy, allSymbol(1, :),  sketchTemplate , ...
                double(filterList_sketch), double(rowList_sketch), double(colList_sketch), length(filterList_sketch));
             sketchTemplate = -sketchTemplate;
             gLow = min(sketchTemplate(:));
             gHigh = max(sketchTemplate(:));
             sketchTemplate = (sketchTemplate-gLow)/(gHigh-gLow);              
             imwrite(sketchTemplate ,[templateFolder 'sketch' num2str(nStroke,'%04d') '.png']); 
       end             
     end  
    
     old_I_syn=I_syn;  % before updating the I_syn, we save the previous syn for computing Z.
     
     tic      
     % multi-chain sampling
     index=0;     
     for j=1:nTileCol
        for i=1:nTileRow
            index=index+1;
            I_syn_single_chain=I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy);             
            response_IB_single_chain = response_IB(index,:);
            I_squared_sum_single_chain = I_squared_sum(index);
            
            % gibbs sampling
            [I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy), I_squared_sum(index), response_IB(index,:)]= ...
                Gibbs_sparseV5(template, filters, I_syn_single_chain, corrBB(1:template.numSelected, 1:template.numSelected), single(c_val_list), ...
                I_squared_sum_single_chain, response_IB_single_chain, nonzero_corrBB_List, lower_bound_rand, upper_bound_rand);        
        end
     end     
     
     % compute z ratio
     logZRatio = computeLogZRatio_sparse_gibbs(old_I_syn,I_syn,template,filters,gradientF_selected,lambdaLearningRate, nTileRow, nTileCol, sx, sy);      
     logZRatioSeries =[logZRatioSeries, logZRatio];
     %
       
     rModel_selected=mean(abs(response_IB),1);   % |<I, B>|        
     
     disp(['one iteration learning time: ' num2str(toc) ' seconds']);

     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5   
     % compute gradient and do graidnet ascent       
     gradientF_selected = rHat_selected-rModel_selected;
    
     SSD = [SSD, mean(abs(gradientF_selected))];
            
     template.selectedLambdas = template.selectedLambdas + lambdaLearningRate * gradientF_selected;     
                
     if mod(iter-1,interval)==0
       % save synthesied image
       gLow = min(I_syn(:));
       gHigh = max(I_syn(:));
       disp([ 'min: ' num2str(gLow) ' max: ' num2str(gHigh)]);
       disp([ 'SSD: ' num2str(SSD(iter))]);
       img = (I_syn-gLow)/(gHigh-gLow);
       imwrite(img,[Synfolder num2str(iter,'%04d') '.png']);
     end
     
end

% re-estimate logz
LogZ = initialLogZ + sum(logZRatioSeries);
disp([' Final LogZ: ' num2str(LogZ)]);

%% solution of least squares min |I_syn - B_selected * C|^2
%% C = inv(B'B)B'I;  I_reconstruct = B*C = B* inv(B'B)B' I_syn, we compute and store R=B inv(B'B)B' as reconstruction matrix
tic;
B=[];
for i=1:template.numSelected    
   Bi = filterTiling(filters{template.selectedFilter(i)}, sx, sy, template.selectedRow(i), template.selectedCol(i));
   B = [B,Bi(:)];
end
R=B*(pinv(B'*B)*(B'));

%%%% save each chain of synthesized images respectively in the last iteration
EachChainSynFolder=[savingFolder 'EachChainSyn/'];
if exist(EachChainSynFolder)
    rmdir(EachChainSynFolder,'s'); 
    mkdir(EachChainSynFolder);
else
    mkdir(EachChainSynFolder);
end
%%%% reconstruct each chain of synthesized images respectively in the last iteration
EachChainReconstructionFolder=[savingFolder 'EachChainReconSyn/'];
if exist(EachChainReconstructionFolder)
    rmdir(EachChainReconstructionFolder,'s'); 
    mkdir(EachChainReconstructionFolder);
else
    mkdir(EachChainReconstructionFolder);
end
I_multi_chain_projected = zeros(nTileRow*sx, nTileCol*sy);
I_multi_chain = zeros(nTileRow*sx, nTileCol*sy);

index=0;
for j=1:nTileCol
   for i=1:nTileRow
     index=index+1;
     crop=I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy);
     gLow = min(crop(:));
     gHigh = max(crop(:));
     img = (crop-gLow)/(gHigh-gLow);       
     imwrite(img,[EachChainSynFolder 'original_chain-' num2str(index,'%04d') '.png']);
     I_multi_chain(1+(i-1)*sx: sx+(i-1)*sx, 1+(j-1)*sy: sy+(j-1)*sy) = img; 
     %% project the synthesized images onto the domain specified by the selected basis
     img_projected = reshape(R*crop(:),sx,sy);  % R = B*(inv(B'*B)*(B'));
     gLow = min(img_projected(:));
     gHigh = max(img_projected(:));
     img_projected = (img_projected-gLow)/(gHigh-gLow);       
     imwrite(img_projected,[EachChainReconstructionFolder 'reconstruct_chain-' num2str(index,'%04d') '.png']);
     I_multi_chain_projected(1+(i-1)*sx: sx+(i-1)*sx, 1+(j-1)*sy: sy+(j-1)*sy) = img_projected; 
   end
end
imwrite(I_multi_chain,[EachChainSynFolder 'raw.png']);
imwrite(I_multi_chain_projected,[EachChainReconstructionFolder 'reconstruct.png']);


% save parts
PartFolder=[Synfolder 'PartSyn/'];
if exist(PartFolder)
    rmdir(PartFolder,'s'); 
    mkdir(PartFolder);
else
    mkdir(PartFolder);
end

index=0;
for i=1:nPartRow
    for j=1:nPartCol
        index=index+1;
        img_parts=img_projected(1+(i-1)*part_sx:part_sx+(i-1)*part_sx, 1+(j-1)*part_sy:part_sy+(j-1)*part_sy);
        gLow = min(img_parts(:));
        gHigh = max(img_parts(:));
        img_parts = (img_parts-gLow)/(gHigh-gLow);  
        imwrite(img_parts,[PartFolder 'part-' num2str(index,'%04d') '.png']);
    end
end


h = figure; 
plot(1:iter,SSD);
title('Convergence');
xlabel('Iterations');
ylabel('SSD');
saveas(h,[savingFolder 'plot.png'],'png');

disp([' Final LogZ: ' num2str(LogZ)]);
bar(1:iter,logZRatioSeries);
title(['Initial Log Z: ' num2str(initialLogZ,'%e') ', final logZ: ' num2str(LogZ,'%e')]);
saveas(gcf,fullfile(savingFolder,'logZRatios.png'));
disp([' The plots saved!' ]);
close all;




