% ========= Scaterplot GCT x iPDC 1 a 5 ==========
% str_varX = 'pValuesGct';          % y-axis p-value holding variable name
% str_measureX = 'GCT';             % Compared measured ID on x-axis
% str_varY = 'pValuesiPdc';         % y-axis p-value holding variable name
% str_measureY = 'iPDC';            % Compared measured ID on y-axis
% str_model = 'Baccala2001a ex3: '; % Simulated model identification
%
% str_suptitle = 'MVGC'; % not used

if ~isempty(strfind(lower(str_varY),lower('iPDC'))) | ...
      ~isempty(strfind(lower(str_varY),lower('iDTF'))),
   str_varY = [str_varY 'min'];
end;

if ~isempty(strfind(lower(str_varX),lower('iPDC'))) | ...
      ~isempty(strfind(lower(str_varX),lower('iDTF'))),
   str_varX = [str_varX 'min'];
end;
%% For small p-values, let us say p = 0, attribute p = eps ~ 2.2204e-16



for kFile = 1:M,
   skF = int2str(kFile);
   for i=1:N,
      for r=1:mm,
         for s=1:mm,
            eval(['tmp1=' str_varX skF '(i,r,s);'])
            eval(['tmp2=' str_varY skF '(i,r,s);'])
            if tmp1 == 0
               eval([str_varX skF '(i,r,s) = eps;'])
            else
               eval([str_varX skF '(i,r,s)=tmp1;'])
            end;
            if tmp2 == 0
               eval([str_varY skF '(i,r,s) = eps;'])
            else
               eval([str_varY skF '(i,r,s)=tmp2;'])
            end;
         end;
      end;
   end;
   
   eval(['rho_' str_measureX '_x_' str_measureY skF '=NaN*zeros(mm,mm);'])
   eval(['pSpearman_' str_measureX '_x_' str_measureY skF '=NaN*zeros(mm,mm);'])
   eval(['angCoeff_' str_measureX '_x_' str_measureY skF '=NaN*zeros(mm,mm);'])
   eval(['angCoeffm_' str_measureX '_x_' str_measureY skF '=NaN*zeros(mm,mm);'])
end;

disp('----------------------------')
disp([str_model ': ' str_measureX ' x '  str_measureY])

for kFile = 1:M,
   skF = int2str(kFile);
   h = figure; %('position',[200 200 500 500]);
   
   set(h,'PaperUnits','inches');
   set(h,'PaperSize',[5 5]);
   
   set(h,'NumberTitle','off','MenuBar','none', ...
      'Name',[str_model ' K=' K{kFile} ': ' str_measureX ' x ' ...
      str_measureY]);
   axis equal
   %suptitle(['n_{s} = ' K{kFile}])
   
   alpha=0.01;
   if kFile == 1,
      hxp=zeros(mm,mm);
   end;
   angCoeff = NaN*zeros(mm,mm); rhos = angCoeff;
   for r=1:mm,
      for s=1:mm,
         subplot2(mm,mm,s+mm*(r-1));
         if r==s,
            text(0.5,0.5,int2str(r),'FontSize',20,  ...
               'HorizontalAlignment','center', ...
               'VerticalAlignment','middle')
            axis off
            set(gca,'xtick',[],'ytick',[],'layer','bottom','box','off')
         else
            eval(['xxx=-log10(' str_varX skF '(:,r,s));'])
            
            eval(['yyy=-log10(' str_varY skF '(:,r,s));'])
            
            xxx=xxx(:);
            yyy=yyy(:);
            inxXxx = xxx > -log10(alpha);
            inxYyy = yyy > -log10(alpha);
            xxx(isinf(xxx)) = 15.6; % ~ -log10(eps)
            yyy(isinf(yyy)) = 15.6; % ~ -log10(eps)
            
            inxXsat = xxx >= 15.59; inxYsat = yyy >= 15.59;
            inxNsat = ~inxXsat & ~inxYsat;
            
            if sum(inxNsat) > 10,
                [rho,pSpearman] = corr(xxx(inxNsat),yyy(inxNsat), ...
                                                       'type','Spearman');
                rhos(r,s) = rho;
                eval(['rho_' str_measureX '_x_' str_measureY skF ...
                                                            '(r,s)=rho;'])
                eval(['pSpearman_' str_measureX '_x_' str_measureY skF ...
                                                     '(r,s)=pSpearman;'])
                x = xxx(inxNsat); y = yyy(inxNsat);
                
                angCoeff(r,s) = sum((x-mean(x)).*(y-mean(y))) ... 
                                                     /sum((x-mean(x)).^2); 
                xxxMedian = median(xxx(inxNsat)); 
                yyyMedian = median(yyy(inxNsat));
                eval(['angCoeff_' str_measureX '_x_' str_measureY skF ...
                    '(r,s) = angCoeff(r,s);'])
                eval(['angCoeffm_' str_measureX '_x_' str_measureY skF ...
                    '(r,s) = mean(yyy(inxNsat))/mean(xxx(inxNsat));'])
            end;
            
            inx1 = inxYyy & ~inxXxx;  xxx1=xxx(inx1); yyy1=yyy(inx1);
            inx2 = ~inxYyy & inxXxx; xxx2=xxx(inx2); yyy2=yyy(inx2);
            inx3 = ~inxYyy & ~inxXxx;  xxx3=xxx(inx3); yyy3=yyy(inx3);
            inx4 = inxYyy & inxXxx; xxx4=xxx(inx4); yyy4=yyy(inx4);

            plot(xxx4,yyy4,'.r'); hold('on'); % Both axes significant
            plot(xxx1,yyy1,'.','Color',[0 .4 0])% Only significant for y-axis
            plot(xxx2,yyy2,'.b', ...          % Only significant for x-axis    
                 xxx3,yyy3,'.k');             % Nonsifnicant for both axes
            set(gca, 'FontName', 'Helvetica','FontSize',7)
            
            axis square
            xymax=max(max(xxx), max(yyy));
            xymax = 18*(xymax>=15)+xymax*(xymax<15 & xymax>=4)+4*(xymax<4);
%             if xymax >= 15,
%                xymax = 18;
%             elseif xymax >= 4,
%                %nop
%             else
%                xymax = 4;
%             end;
            
            xyinf = -0.2*xymax/4;
            if xymax < 15,
               xymax = 1.065*xymax;
            end;
            
            axis([xyinf xymax xyinf xymax])
            % Threshold lines for x and y-axes.
            plot([xyinf xymax],[2 2],'--b', 'LineWidth',[1.011])
            plot([2 2],[xyinf xymax],'--b', 'LineWidth',[1.011])
            
            
            if isempty(strfind(lower(str_varY),lower('iDTF'))), 
               %notDTF
               if isempty(strfind(lower(str_varX),lower('iDTF'))), 
                  %notDTF
                  if trueStruct(r,s) == 1,
                     eval(['iFNy = FN' str_FNY '(r,s,kFile);'])
                     if flgAdj
                        eval(['iFNyadj = FN' str_FNY 'adj(r,s,kFile);'])
                        ht1=text(0.02*xymax,0.85*xymax,['FN=' ...
                           int2str(iFNy) '/' int2str(iFNyadj)]);
                     else
                        ht1=text(0.02*xymax,0.85*xymax,['FN=' int2str(iFNy)]);
                     end;
                     eval(['iFNx = FN' str_FNX '(r,s,kFile);'])
                     ht2=text(0.37*xymax,0.05*xymax,['FN=' int2str(iFNx)]);
                  else
                     eval(['iFPy = FP' str_FNY '(r,s,kFile);'])
                     if flgAdj
                        eval(['iFPyadj = FP' str_FNY 'adj(r,s,kFile);'])
                        
                        ht1=text(0.02*xymax,0.85*xymax,['FP=' int2str(iFPy) ...
                           '/' int2str(iFPyadj)]);
                     else
                        ht1=text(0.02*xymax,0.85*xymax,['FP=' int2str(iFPy)]);
                     end;
                     eval(['iFPx = FP' str_FNX '(r,s,kFile);'])
                     ht2=text(0.37*xymax,0.05*xymax,['FP=' int2str(iFPx)]);
                  end;
               else % x-axis is iDTF
                  % But y-axis is not DTF
                  if trueStruct(r,s) == 1,
                     eval(['iFNy = FN' str_FNY '(r,s,kFile);'])
                     ht2=text(0.02*xymax,0.85*xymax,['FN=' int2str(iFNy)]);
                  else
                     eval(['iFPy = FP' str_FNY '(r,s,kFile);'])
                     ht2=text(0.02*xymax,0.85*xymax,['FP=' int2str(iFPy)]);
                  end;
                  % Indedd x-axis is DTF
                  if trueReachab(r,s) == 1,
                     eval(['iFNx = FN' str_FNX '(r,s,kFile);'])
                     if flgAdj
                        eval(['iFNxadj = FN' str_FNX 'adj(r,s,kFile);'])
                        ht1=text(0.37*xymax,0.05*xymax,['FN=' ...
                           int2str(iFNx) '/' int2str(iFNxadj)]);
                     else
                        ht1=text(0.37*xymax,0.05*xymax,['FN=' int2str(iFNx)]);
                     end;
                  else
                     eval(['iFPx = FP' str_FNX '(r,s,kFile);'])
                     if flgAdj
                        eval(['iFPxadj = FP' str_FNX 'adj(r,s,kFile);'])
                        
                        ht1=text(0.37*xymax,0.05*xymax,['FP=' int2str(iFPx) ...
                           '/' int2str(iFPxadj)]);
                     else
                        ht1=text(0.37*xymax,0.05*xymax,['FP=' int2str(iFPx)]);
                     end;
                  end;
               end;
            else  % Y-axis is iDTF.
               % First the x-axis that is not DTF
               if trueStruct(r,s) == 1,
                  eval(['iFNx = FN' str_FNX '(r,s,kFile);'])
                  ht2=text(0.37*xymax,0.05*xymax,['FN=' int2str(iFNx)]);
               else
                  eval(['iFPx = FP' str_FNX '(r,s,kFile);'])
                  ht2=text(0.37*xymax,0.05*xymax,['FP=' int2str(iFPx)]);
               end;
               % Now let's deal with y-axis is DTF, and requires reachability info
               if trueReachab(r,s) == 1,
                  eval(['iFNy = FN' str_FNY '(r,s,kFile);'])
                  if flgAdj
                     eval(['iFNyadj = FN' str_FNY 'adj(r,s,kFile);'])
                     ht1=text(0.02*xymax,0.85*xymax,['FN=' ...
                        int2str(iFNy) '/' int2str(iFNyadj)]);
                  else
                     ht1=text(0.02*xymax,0.85*xymax,['FN=' int2str(iFNy)]);
                  end;
               else
                  eval(['iFPy = FP' str_FNY '(r,s,kFile);'])
                  if flgAdj
                     eval(['iFPyadj = FP' str_FNY 'adj(r,s,kFile);'])
                     
                     ht1=text(0.02*xymax,0.85*xymax,['FP=' int2str(iFPy) ...
                        '/' int2str(iFPyadj)]);
                  else
                     ht1=text(0.02*xymax,0.85*xymax,['FP=' int2str(iFPy)]);
                  end;
               end;
            end;
            
            set(ht1,'FontName','Helvetica','FontSize',7, ...
               'FontWeight','bold','Color',[0 0 .5])
            set(ht2,'FontName','Helvetica','FontSize',7, ...
               'FontWeight','bold','Color',[0 0 .5])
            grid
         end;
      end;
      % Extract axes handles of all subplots from the figure
      axesHandles = findobj(get(h,'Children'), 'flat','Type','axes');
      % Set the axis property to square
      axis(axesHandles,'square')
   end;
   pos = get(gcf,'Position');
   dtmp=abs(pos(3)-pos(4));
   pos(2)=pos(2)-dtmp;
   tmp= max(pos(3:4)); pos(3)=tmp; pos(4)=tmp;
   set(gcf,'Position',pos, 'PaperPosition', [0.25 2.5 6 5.7468], ...
      'PaperSize', [6 5.7468]);
   
   supAxes=[.12 .10 .84 .84];
   [ax,h1]=suplabel(['-log_{10} (p_{' str_measureX '})'],'x',supAxes);
   set(h1, 'FontSize',[14])
   [ax,h2]=suplabel(['-log_{10} (p_{' str_measureY '})'],'y',supAxes);
   set(h2, 'FontSize',[14])
   
   if flgReport,
      str_SupTitle = [str_model ' K=' K{kFile} ': ' str_measureX ' x ' ...
         str_measureY];
      %       h3 = suptitle(str_SupTitle);
      %       set(h3,'FontSize',6,'FontWeight','normal')
      
      [ax,h3]=suplabel(str_SupTitle,'t');
      set(h3,'FontSize',10)
      
      format compact
      disp('----------------------------')
      disp(['K = ' K{kFile} ':'])
      disp('b = ')
      disp(angCoeff)
      disp(['rho = '])
      disp(rhos)
      
         angCoeff = NaN*zeros(mm,mm); rho = angCoeff;
   end;
end;
clear x xxx xxxMedian xxx1 xxx2 xxx3 xxx4 inxXxx y yyy yyyMedian ...
      yyy1 yyy2 yyy3 yyy4 inxYyy inxXsat inxYsat inxNsat inx1 inx2 ...
      inx3 inx4 xymax xyinf
