function [template,currSample, LogZ]=sparseFRAMElearn(template, nIteration,filters,rHat,sx,sy,halfFilterSize,locationShiftLimit, nTileRow,nTileCol,epsilon,L,lambdaLearningRate,numSample,isSaved,savingFolder)



nValidFeature= sum((sx-2*halfFilterSize-2*locationShiftLimit).*(sy-2*halfFilterSize-2*locationShiftLimit));
SSD=zeros(nIteration,1);
%initialLogZ = log((2*pi))*(nValidFeature/2);
initialLogZ =0 ;
logZRatioSeries = zeros(nIteration,1);


% the field of enlargeTemplate stores a big tiled template for multiple chains HMC sampling 
template.enlargedTemplate = enlargingTemplate(nTileRow, nTileCol, sx, sy, template);
prevSample = single(randn(sx*nTileRow,sy*nTileCol));  

rHat_selected = zeros(1,template.numSelected);
for i=1:template.numSelected
   rHat_selected(i)=rHat{template.selectedFilter(i)}(template.selectedRow(i),template.selectedCol(i));
end

gradientF_selected=zeros(1,template.numSelected);
for iter = 1:nIteration
    disp( [ 'iteration: ' num2str(iter)]);
    tic  
  
     [rModel_selected, currSample]=multiChainHMC_sparse(template, filters, prevSample, epsilon, L, numSample, nTileRow, nTileCol);  
    
     % compute z ratio
     logZRatio = computeLogZRatio_sparse(prevSample,currSample,template,filters,gradientF_selected,lambdaLearningRate, nTileRow, nTileCol);
     logZRatioSeries(iter)=logZRatio;
     %
     
     disp(['one iteration learning time: ' num2str(toc) ' seconds']);

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

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

%%%save the multiple chains of (projected) synthesized images separately
%% 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
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*(inv(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
curr_projected=zeros(nTileRow*sx,nTileCol*sy);
currSample2=zeros(nTileRow*sx,nTileCol*sy);
index=0;
for j=1:nTileCol
   for i=1:nTileRow
     index=index+1;
     crop=currSample(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']);
     currSample2(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'));
     curr_projected(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy)=img_projected;
     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']);
     curr_projected(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy)=img_projected;
   end
end
imwrite(currSample2,[EachChainSynFolder 'raw.png']);
imwrite(curr_projected,[EachChainReconstructionFolder 'reconstruct.png']);

%%%%%%



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

  disp([' Final LogZ: ' num2str(LogZ)]);
  bar(1:nIteration,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;
end



