function  NSIMCA_QTplots(XfitResults,XvalResults,XtestResults, Method, numF, classid, classts)
% this function reads output from NSIMCA and:
% 
% - make H/Q plot with different limits for a model with numF dimensionality 
%% INPUT
%% XfitResults, XvalResults,XtestResults structured array(factors,
%% class) from NSIMCA
%% Method ('par' or 'tuc')
%% numF =
%% [numFc1; numFc2; ....] for a PARAFAC model
%% [numFc1m1 numFc1m2 numFc1m3...;....] for a Tucker  model wher c is for
%% class and m is for mode
%% classid = vector containing class assignation as numeric
%% classts = vector containing class assignation as numeric for the test
%% set samples
%
%** NSIMCA_QTplots(XfitResults,XvalResults,XtestResults, Method, numF, classid, classts)

%% OUTPUT
%% figures
%%
%%  1. H_fit/Q_fit 
%%  2. H_CV/Q_CV


%get number of factors and number of classes
[nf nc]= size(XvalResults);

% loop over classes
for indc=1:nc
    fmin =XfitResults(1,indc).finit;
    Fmax =XfitResults(1,indc).ftot;
    DimXcal=XfitResults(1,indc).DimXcal;
    ord=size(DimXcal,2);
    

    if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
        indf=numF(indc);
    else
        for i=fmin:Fmax
            PossibleCombs(i,:)=XfitResults(i,indc).PossibleCombs;
        end
        [nn indf]=intersect(PossibleCombs(:,1:ord),numF(indc,1:ord),'rows');
    end    
    
    H_fit=[XfitResults(indf,indc).hT_fit];    
    Q_fit=[XfitResults(indf,indc).Qfit];  
    Qlimit_fit=[XfitResults(indf,indc).Qlimit_fit];  
    
    Qlimit_fitAP=[XfitResults(indf,indc).Qlimit_fitAP];  
    Hlim_fitAP=[XfitResults(indf,indc).Hlim_fitAP];  
    limF_hTfit=[XfitResults(indf,indc).limF_hTfit];  

    
    %% CrossValidated
    H_CV=[XvalResults(indf,indc).hT_CV];    
    Q_CV=[XvalResults(indf,indc).Q_CV];
          
    limAP_QCV=[XvalResults(indf,indc).limAP_QCV];
    limAP_hT_CV=[XvalResults(indf,indc).limAP_hT_CV];  
    %
    Hrlim_CV = [XvalResults(indf,indc).hTrlim_CV];
    Qrlim_CV = [XvalResults(indf,indc).Qrlim_CV]; % added Feb 07
    
    
    %% other  classes
    if nc > 1
        H_others=[XtestResults(indf,indc).hT_test];
        Q_others=[XtestResults(indf,indc).Qtest];
    else
        H_others=[];
        Q_others=[];
    end
    %% external test set 
    if nargin > 6
        H_test_est=[XtestResults(indf,indc).hT_test_test];
        Q_test_est=[XtestResults(indf,indc).Qtest_test];
    else
        H_test_est=[];
        Q_test_est=[];
    end

    %
    
    %%
    
    if ~strcmp(lower(Method(1:3)),'tuc') %% Method PARAFAC
        titlef=strcat('Model of class= ',int2str(indc),'Factors= ',int2str(indf));
    else
        titlef=strcat('Model of class= ',int2str(indc),'Factors= ',int2str(nn));
    end
    % H/Q fit/CV
    if ~exist('classts')
        classts=[];
    end
    makeplot(H_fit,Q_fit,H_CV,Q_CV,H_others,Q_others,H_test_est,Q_test_est,[limF_hTfit Hlim_fitAP Hrlim_CV limAP_hT_CV],[Qlimit_fit Qlimit_fitAP Qrlim_CV limAP_QCV],titlef, classid, indc, 'H', classts)
    fnam=strcat('c_',int2str(indc),'_fitCV')
    hgsave(fnam)   

   
end


    %%
 

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

function makeplot(x,y,xcv, ycv, xo,yo,xt,yt,limx,limy,titlef,ci,i, xlab, cts)
    
ntr=length(x);
if ~isempty(xo); no=length(xo);end
if ~isempty(xt);nts=length(xt);end
nc=nonzeros(unique(ci));
ci_c=find(ci);
classid=ci(ci_c); 
nc=setdiff(nc,i);
indtr=find(classid==i);
ind_o=setdiff([1:length(classid)],indtr);
class_o=classid(ind_o);

s=['.r';'.b'; '*g'; 'xm'; '+c'; '^y'; 'pk';'hr'];
scv=['or';'ob'; 'og'; 'om'; 'oc'; 'oy'; 'ok';'<r'];
st=['dr';'db';'dg';'dm';'dc';'dy';'dk';'>r'];
sl=['-k';'-g';'-r';'-m'; '-c'];
figure;
plot(x,y,s(i,:))
hold on;plot(xcv,ycv,scv(i,:))
title(titlef);
if ~isempty(xt);hold on; plot(xt(find(cts==i)),yt(find(cts==i)),st(i,:));end
t{1}=strcat('train',int2str(i));
t{2}='CV';
if ~isempty(cts)
    t{3}=strcat('test',int2str(i));
    tind=3;
else
    tind=2;
end
xlabel(xlab);
ylabel ('Q');
%%leg=3
%ncp=nonzeros(unique(cts));
for j=1:length(nc)
    ii=find(class_o==nc(j));
    if ~isempty (cts); its=find(cts==nc(j));end
    if ~isempty(xo);hold on; plot(xo(ii),yo(ii),s(nc(j),:)); t{tind+2*(j-1)+1}=strcat('train',int2str(nc(j))); end
    if ~isempty(xt);hold on; plot(xt(its),yt(its),st(nc(j),:)); t{tind+2*j}=strcat('test',int2str(nc(j))); end
end
lt=length(t);
% limts in fit and CV
if xlab=='H'
    t{lt+1}='H_l_i_mfitF/ Q_l_i_mfit' ;
    t{lt+2}='H_l_i_mfitAP/ Q_l_i_mfitAP' ;
    t{lt+3}='H_l_i_mCVper/ Q_l_i_mCVper'; 
    t{lt+4}='H_l_i_mCVAP/ Q_l_i_mCVAP';
end

for il=1:length(limy)
    ly=limy(il)
    hline(ly,sl(il,:));
end
for il=1:length(limx)
    lx=limx(il)
    vline(lx,sl(il,:));
end
axis tight
legend(t);
