function [rModel_selected, pSample]=multiChainLangevin_sparse(template, filters, pCurrSample,stepsize,L, nRuns, nRow, nCol)
% function to sample from multiple parallel chains, with gaussian prior
% nRow and nCol defines the size of the chain grid.
% This is a CPU version of sparse multiple chain HMC sampling 
% we assume the pCurrSample is already a large canvas containing a grid of nRow by nCol sampled images

%% 0. prepare the enlarged template
pTemplate = template.enlargedTemplate;   % p for paralle
pTemplate.selectedLambdas = repmat(template.selectedLambdas,1,nRow*nCol);
numSelected = template.numSelected;


%% 1. correct the bounary condition, and run HMC
pSample = pCurrSample;
for iRun= 1:nRuns
    pSample = multi_Langevin(pTemplate,filters,pSample,stepsize,L);
end
%% 2 compute the rModel_selected

[pRSample] = sparseFeatureExtraction(single(pSample),filters, pTemplate,true);
rModel_selected = mean(reshape(pRSample,numSelected,[]),2);
rModel_selected = rModel_selected';

end % end of main function


function I=multi_Langevin(pTemplate, filters, pCurrSample, stepsize, L)
% The main routine for sampling using HMC.
%
%% The traditional leapFrog steps


I=pCurrSample;
for i=1:L
   noise = normrnd(0,1,size(pCurrSample));
   I = I- (stepsize^2)*0.5 * grad_U_sparse(I, filters, pTemplate) + noise;    
end

end % gpuMultiHMC
        




