function varargout = sodaPlotSeq(sodaPar,sequences,metropolisRejects,varargin)


xLim = [1,max(max(sequences(:,sodaPar.iterCol,:)))];
xScale = repmat({'lin'},[sodaPar.nOptPars,1]);
yScale = repmat({'lin'},[sodaPar.nOptPars,1]);


plotMode = 'subaxes';
printToFile = false;
printToPNG = false;
printToEPS = false;
checkered = false;
showRejects = true;

authorizedOptions={'plotMode',...
                   'xScale',...
                   'yScale',...
                   'printToFile',...
                   'printToPNG',...
                   'printToEPS',...
                   'checkered',...
                   'showRejects',...
                   'xLim'};
                 
% overrule default options with user-specified options:                 
sodaParsePairs
for r=1:numel(sodaPar.parCols)
    xScale2{r,1}=r;
    xScale2{r,2}=cellfind(xScale,r);
    yScale2{r,1}=r;
    yScale2{r,2}=cellfind(yScale,r);
end
xScale = xScale2;
yScale = yScale2;
clear xScale2
clear yScale2


nPars = numel(sodaPar.parCols);
plotHandlesSequences = repmat(NaN,sodaPar.nCompl,nPars);

colorList = sodaMakeColors(sodaPar);

iGeneration = (size(sequences,1)-sodaPar.nSamplesPerCompl)/sodaPar.nOffspringPerCompl;
calcConvFromIter = (sodaPar.nSamplesPerCompl +...
                   floor(iGeneration*(1-sodaPar.convUseLastFraction))*...
                   sodaPar.nOffspringPerCompl)*sodaPar.nCompl;


for p=1:nPars
    
    switch plotMode
        case 'figures'
            figure('name',[mfilename,':',sodaPar.parNames{p}],'numbertitle','off')
        case 'subaxes'
            set(gcf,'name',mfilename,'numbertitle','off')
            subaxes(nPars,1,p,'borderBottom',0.02,...
                              'spacingRight',0.01,...
                              'spacingLeft',0.02)
        case 'subplots'
            set(gcf,'name',mfilename,'numbertitle','off')
            subplot(nPars,1,p)
    end

    
    if checkered
        
        nStratRand = floor(sodaPar.nSamples^(1/sodaPar.nOptPars));
        yv = linspace(sodaPar.parSpaceLoBound(p),sodaPar.parSpaceHiBound(p),nStratRand+1);
        nBlocksX = (nStratRand^sodaPar.nOptPars)/(nStratRand^(sodaPar.nOptPars-p));
        
        for ix = 1:nBlocksX% block counter x dir
            
            XW = nStratRand^(sodaPar.nOptPars-p);
            XLB = (ix-1)*XW + 0.5;
            XRB = ix*XW + 0.5;
            XRT = XRB;
            XLT = XLB;
            
            for iy = 1:nStratRand

                YLB = yv(iy);
                YRB = YLB;
                YRT = yv(iy+1);
                YLT = YRT;

                hp=patch([XLB,XRB,XRT,XLT,XLB],[YLB,YRB,YRT,YLT,YLB],'k');
                if (isodd(ix)&isodd(iy)) | (~isodd(ix)&~isodd(iy))
                    set(hp,'edgecolor','none','facecolor',[1,1,1]*0.9)
                else
                    set(hp,'edgecolor','none','facecolor',[1,1,1])
                end
                hold on
                
            end                 
                 
        end
    end
    
    
    nEvals = max(max(shiftdim(sequences(:,1,:),2)));
    yLim = [sodaPar.parSpaceLoBound(p),sodaPar.parSpaceHiBound(p)];
    x = [sodaPar.nSamples:sodaPar.nOffspring:nEvals]+0.5;
    plotHandlesGenerations = plot([x;x],yLim,'color',[0.9,0.9,0.9]);
    hold on
    plotHandlesGenerations = [plotHandlesGenerations;...
        plot([1;1]*calcConvFromIter+0.5,yLim,'color',[0.5,0.5,0.5])];
    set(gca,'xscale',xScale{p,2},'yscale',yScale{p,2})    
    for k=1:sodaPar.nCompl
        
        if showRejects
            plotHandlesRejects(k,p) = plot(metropolisRejects(:,1,k),metropolisRejects(:,p+1,k),'s',...
                                          'markersize',3,...
                                     'markerfacecolor','w',...
                                     'markeredgecolor',colorList(k,:));
        end

        plotHandlesSequences(k,p) = plot(sequences(:,1,k),sequences(:,p+1,k),'s',...
                                      'markersize',3,...
                                 'markerfacecolor',colorList(k,:),...
                                 'markeredgecolor',colorList(k,:));
    end
    set(gca,'ylim',yLim,'xlim',xLim,'layer','top')

    
    if isfield(sodaPar,'parNamesTex')&&...
            numel(sodaPar.parNamesTex)>=(p)&&...
            ~isempty(sodaPar.parNamesTex{p})
        ylabel(sodaPar.parNamesTex{p},'interpreter','tex',...
            'rotation',90,'horizontalalignment','right')
    else
        ylabel(sodaPar.parNames{p},'interpreter','none',...
            'rotation',90,'horizontalalignment','right')
    end

    box on
    hold off
    if ~strcmp(plotMode,'figures')
        if p~=nPars(end)
            set(gca,'xticklabel',[])
        end
    end
end

xlabel('model evaluations')




if nargout==2
   varargout{1}=plotHandlesSequences;
   varargout{2}=plotHandlesGenerations;   
end

if printToPNG||printToFile
    print(mfilename,'-dpng','-r300')
end
if printToEPS||printToFile
    print(mfilename,'-depsc','-r300')
end


% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %


function scaleType=cellfind(C,subplotIx)
scaleType='linear';
for r=1:size(C,1)
    if C{r,1}==subplotIx
        scaleType = C{r,2};
        return
    end
end

