function gridOfInnColorTwoSpace(sodaPar,mStatesInn,parVec,varargin)


load('./../eliminate.mat')


yUnitStr = 'units';
modeStr = 'innovations';
figNumber = 67445;
bgColor = 'auto';
showParameters = true;
showLegend = true;
showLegendAxes = true;
yLimMin = NaN;
yLimMax = NaN;
xLabel = NaN;
yLabel = NaN;
linewidth = 1.2

authorizedOptions = {'modeStr','figNumber','yUnitStr',...
    'showParameters','showLegend','bgColor','yLimMin','yLimMax',...
    'xLabel','yLabel','linewidth','lineColor1','lineColor2'};
parsepairs


useDefaultLineColor = ~exist('lineColor1','var')|~exist('lineColor2','var');

if all(isnan(mStatesInn(:)))
    disp('Only NaNs in innovations array. Aborting.')
    return
end

parValueStr = '';
nPars = sodaPar.nOptPars;
for iPar=1:nPars
    parValueStr = [parValueStr,sodaPar.parNamesTex{iPar},...
        ' = ',num2str(parVec(iPar)),'; ',char(10)];
end

% restoredefaultpath
% 
% clear
% close all
% clc
% 
% pwdStr = pwd;
% Ix = findstr(pwdStr,filesep);
% thisDisk = pwdStr(1:Ix(1));
% 
% addpath(fullfile(thisDisk,'horizon','matlab-tools'))
% 
% clear Ix 
% clear pwdStr
% 
% load('./../job.mat')
% load('./../innovations.mat')



[nTimes,nStatesKF,nMembers] = size(mStatesInn);

canvas.l = 0.12;
canvas.b = 0.06;
canvas.w = 0.53;
canvas.h = 0.85;


ax.l = 0.05;
ax.b = 0.05;
ax.w = 0.9;
ax.h = 0.9;

nRows = 15;
nCols = 7;
xRes = 3.0;
yRes = 3.0;

colWidth = canvas.w/nCols;
rowHeight = canvas.h/nRows;

axHandles = repmat(NaN,[nRows,nCols]);


timeDim = 1;
stateDim = 2;
memberDim = 3;
meanInn = repmat(0,[nTimes,nRows*nCols])
meanInn(1:nTimes,obsKFIx) = mean(mStatesInn,memberDim);
cumPosInn2D = reshape(sum((meanInn(2:nTimes,1:nRows*nCols)>0).*meanInn(2:nTimes,1:nRows*nCols),1),[nRows,nCols]);
cumNegInn2D = reshape(sum((meanInn(2:nTimes,1:nRows*nCols)<0).*meanInn(2:nTimes,1:nRows*nCols),1),[nRows,nCols]);

% cumMeanInn = reshape(sum(meanInn(2:nTimes,1:nStatesKF),timeDim),[nRows,nCols]);

colorMaxAbs = max([max(cumPosInn2D(:)),-min(cumNegInn2D(:))]);

if colorMaxAbs==0
    colorMax = 1.0;
    colorMin = -1.0;
else
    colorMax = colorMaxAbs;
    colorMin = -colorMaxAbs;
end

switch modeStr
    case 'error'
        colorLB = shiftdim([1,0,1],-1);   % magenta
        colorLT = shiftdim([1,1,1],-1);   % white
        %colorRB = shiftdim([0,0,1],-1);   % blue
        colorRB = shiftdim([0,0,0]+0.4,-1);   % black
        colorRT = shiftdim([0,1,1],-1);   % cyan
        %colorMM = shiftdim([0,0,1],-1);   % blue
        colorMM = shiftdim([0,0,0]+0.4,-1);   % black        
        colorLM = mean([colorLB,colorLT],2);
        colorRM = colorMM;%mean([colorRB,colorRT],2);
        colorMB = colorMM;%mean([colorLB,colorRB],2);
        colorMT = mean([colorLT,colorRT],2);
        colorsAll = [[colorLB,colorMB,colorRB];[colorLM,colorMM,colorRM];[colorLT,colorMT,colorRT]];
        if useDefaultLineColor
            lineColor1 = [0,0,0];
            lineColor2 = [0,0,1];
        end
    case {'innovations','state updates'}
%         colorLB = shiftdim([0,0.5,0],-1); % olive        
%         colorLT = shiftdim([1,1,1],-1);   % white
%         %colorRB = shiftdim([1,1,0],-1);   % yelllow
%         colorRB = shiftdim([0,0,0]+0.4,-1);   % black
%         colorRT = shiftdim([1,0,0],-1);   % red        
%         %colorMM = shiftdim([1,1,0],-1);   % yellow
% %         colorMM = shiftdim([0,0,0]+0.4,-1);   % black
% %         colorLM = mean([colorLB,colorLT],2);
% %         colorRM = colorMM;%mean([colorRB,colorRT],2);
% %         colorMB = colorMM;%mean([colorLB,colorRB],2);
% %         colorMT = mean([colorLT,colorRT],2);
% %         colorsAll = [[colorLB,colorMB,colorRB];[colorLM,colorMM,colorRM];[colorLT,colorMT,colorRT]];
%         colorsAll = [[colorLB,colorRB];[colorLT,colorRT]];
% 
        if useDefaultLineColor
            lineColor1 = [0,0,0];
            lineColor2 = [0,0,1];
        end
    case 'states'
        colorLB = shiftdim([0.5,0.25,0],-1); % brown
        colorLT = shiftdim([1,1,1],-1);   % white
        colorRB = shiftdim([0,0,0]+0.4,-1);   % black
        colorRT = shiftdim([1,1,0],-1);   % yellow%aquamarine
        colorMM = shiftdim([0,0,0]+0.4,-1);   % black
        colorLM = mean([colorLB,colorLT],2);
        colorRM = colorMM;%mean([colorRB,colorRT],2);
        colorMB = colorMM;%mean([colorLB,colorRB],2);
        colorMT = mean([colorLT,colorRT],2);
        colorsAll = [[colorLB,colorMB,colorRB];[colorLM,colorMM,colorRM];[colorLT,colorMT,colorRT]];
        if useDefaultLineColor
            lineColor1 = [0,0.5,1];
            lineColor2 = [0,0,1];
        end
end
% [xColorSpace,yColorSpace]=meshgrid(linspace(0,colorMax,64),linspace(colorMin,0,64));
% for iColorPage = 1:3
%     cMap2D(:,:,iColorPage) = interp2([0,colorMax],[colorMin,0],colorsAll(:,:,iColorPage),xColorSpace,yColorSpace);
% end

[xColorSpace,yColorSpace]=meshgrid(linspace(0,colorMax,64),linspace(colorMin,0,64));
for iColorPage = 1:3
%     cMap2D(:,:,iColorPage) = interp2([0,0.33,1]*colorMax,[1,0.33,0]*colorMin,colorsAll(:,:,iColorPage),xColorSpace,yColorSpace);
%     cMap2D(:,:,iColorPage) = interp2([0,1]*colorMax,[1,0]*colorMin,colorsAll(:,:,iColorPage),xColorSpace,yColorSpace);
end


if ismember(figNumber,get(0,'Children'))
    close(figNumber)
end
subplotScreen(1,2,1,'figureNumber',figNumber)
axes('position',[canvas.l,canvas.b,canvas.w,canvas.h])
set(gca,'color','none','xtick',[],'ytick',[],'box','off',...
    'xcolor',[1,1,1]*0.8,'ycolor',[1,1,1]*0.8,'xaxislocation','top',...
    'yaxislocation','right')
xlabel(['\bfx: along-slope distance [m] \rightarrow','  ',char(10)],'color','k','fontsize',16)
ylabel([char([10]),'  ',char([10]),'  ','                              \bfy: upslope distance [m] \rightarrow'],'color','k','fontsize',16)
xlim = [min(sodaPar.simPoints{1,2}(2:end)),max(sodaPar.simPoints{1,2}(2:end))];
if colorMaxAbs==0
    ylim = [-1,1];
else
    ylim = [min(meanInn(:)),max(meanInn(:))];
end

if ~isnan(yLimMin)
    ylim(1)=yLimMin;
end
if ~isnan(yLimMax)
    ylim(2)=yLimMax;
end

    

for iRow=1:nRows
    for iCol=1:nCols
        
        L = canvas.l + (iCol-1)*colWidth + ax.l*colWidth;
        B = canvas.b + (iRow-1) * rowHeight + ax.b*rowHeight;
        W = ax.w*colWidth;
        H = ax.h*rowHeight;
        axHandles(iRow,iCol)=axes('position',[L,B,W,H]);

        stateIndex = (iCol-1)*nRows+iRow;
        if ~ismember(stateIndex,obsKFIx)
            
            set(gca,'xcolor',0.5*[1,1,1],'ycolor',0.5*[1,1,1],...
                'xtick',[],'ytick',[],'box','on')
            myColor = 0.8*[1,1,1];
        else
            plot(xlim,[0,0],'linestyle','--','color',[1,1,1]*0.4,'linewidth',linewidth)
            hold on
            plot(sodaPar.simPoints{1,2}(2:end),meanInn(2:end,stateIndex),...
                'linestyle','-','color',lineColor1,'linewidth',linewidth);
            myColor = interp1([max(cumNegInn2D(obsKFIx));colorMin],[1,1,1;0,0.5,0],...
                                        cumNegInn2D(iRow,iCol));
            
        end
        
        if ~(iCol==1&iRow==1)
            set(axHandles(iRow,iCol),'xticklabel','')
            set(axHandles(iRow,iCol),'yticklabel','')
        else
            if isnan(xLabel)
                xlabel(xLabel,'fontsize',12)
            else
                xlabel(xLabel,'fontsize',12)
            end
            if isnan(yLabel)
                ylabel(modeStr,'fontsize',12)
            else
                ylabel(yLabel,'fontsize',12,'color',lineColor1)
                set(gca,'ycolor',lineColor1)
            end
        end
        if iCol==nCols
            ylabel([num2str((iRow-1)*xRes,'%.1f')],'rotation',0,'horizontalalignment','left')
            set(gca,'yaxislocation','right')
        end
        if iRow==nRows
            xlabel([num2str((iCol-1)*yRes,'%.1f')])
            set(gca,'xaxislocation','top')        
        end
        set(get(gca,'xlabel'),'color',[0,0,0])
        set(get(gca,'ylabel'),'color',[0,0,0])
        
        % determine color:
%         if ischar(bgColor)&&strcmp(bgColor,'auto')
%             for iColor = 1:3
%                 myColor(1,iColor) = interp2([0,0.33,1]*colorMax,...
%                                             [1,0.33,0]*colorMin,...
%                                             colorsAll(:,:,iColor),...
%                                             cumPosInn2D(iRow,iCol),...
%                                            cumNegInn2D(iRow,iCol));
% %             end
%         else
%             if numel(bgColor)==1
%                 myColor = bgColor;
%             else
%                 myColor = bgColor(stateIndex,1:3);
%             end
%         end

        set(axHandles(iRow,iCol),...
            'xlim',xlim,...
            'ylim',ylim,...
            'xtick',xlim,...
            'ytick',ylim,...
            'color',myColor)
        
    end
end


L = canvas.l + canvas.w + (2.25)*colWidth;
B = canvas.b + (0.5) * rowHeight + ax.b*rowHeight;
W = 1*colWidth+ax.w*colWidth;
H = rowHeight*1 + ax.h*rowHeight;

if showLegendAxes
    legendAxes=axes('position',[L-1.25*colWidth,B-0.5*rowHeight,W+1.5*colWidth,H+3*rowHeight]);
    set(legendAxes,'xtick',[],'ytick',[],'color',[0.9,0.9,0.9],'box','on')
    text(0.1,0.9,'\bfLEGEND','fontsize',14)
end
if showLegend
axes('pos',[0.83,0.12,0.05,0.15])

%     colorbarAxes=axes('position',[L,B,W,H]);
%     image([0,colorMax],[colorMin,0],cMap2D)
    myColorMap = interp1([max(cumNegInn2D(obsKFIx));colorMin],[1,1,1;0,0.5,0],...
                                linspace(max(cumNegInn2D(obsKFIx)),colorMin,64));
    image(0,[max(cumNegInn2D(obsKFIx)),colorMin],permute(myColorMap,[1,3,2]))
    set(gca,'ydir','normal','xaxislocation','bottom','xtick',[],'ytick',[-2,-1,max(cumNegInn2D(obsKFIx))])
%     xlabel(['per-node,',char(10),'time-accumulated,',char(10),'positive ',modeStr,' [',yUnitStr,']'],'fontsize',8)
    xlabel(['per-node, time-accumulated,',char(10),'negative ',modeStr,' [',yUnitStr,']'],'fontsize',8)
end

if showParameters
    L = canvas.l + (-2-1)*colWidth + ax.l*colWidth;
    B = canvas.b + (11-1) * rowHeight + ax.b*rowHeight;
    W = ax.w*colWidth;
    H = ax.h*rowHeight;
    axes('position',[L,B,W,H]);
    text(0,0,parValueStr,'horizontalalignment','center')
    axis off
end
% cBarValues = linspace(max(cumMeanInn(:)),min(cumMeanInn(:)),64)';
% imagesc(1,cBarValues,cBarValues)
% caxis([-1,1]*colormax)
% set(gca,'ydir','normal','xtick',[])
% title(['cumulative',char(10),'mean innovation [m]'])
% colormap(cMap)

% 
% set(gcf,'inverthardcopy','off','paperpositionmode','auto')
% 
% print('grid-of-innovations.png','-dpng','-r300')
% print('grid-of-innovations.eps','-depsc2','-r300','-loose')

