function plotipls=iplsplot(Model,labeltype,optimal_lv_global,max_yaxis,plottype)
%
%  iplsplot plots results from iPLS analysis in a bar plot
%
%  Input:
%  Model (the output from ipls.m)
%  labeltype designates whether you want:
%    interval number ('intlabel'),
%  	 variable number ('varlabel')
%    wavelength number ('wavlabel')
%  optimal_lv_global: optional, the number of PLS components chosen for full spectrum model
%     if not given or given by [] the first RMSECV/RMSEP minimum is chosen
%  max_yaxis is optional (can be used to control scaling of the iPLS plot)
%  plottype is optional:
%    'Cum' (default)
%    'Cum2' (with the RMSE values of the two preceeding PLS components)
%    'n' (where is a locked number of PLS components)
%
%  Copyright, Chemometrics Group - KVL, Copenhagen, Denmark
%  Lars Nørgaard, December 2004
%
%  iplsplot(Model,labeltype,optimal_lv_global,max_yaxis,plottype);
%   iplsplot(Model,''intlabel'',[],[],''Cum2'')
%   iplsplot(Model,''intlabel'')
%
%%
%   Function to calculates the sub_pls_val for PLS modelling with selected validation method
%   PLSmodel = sub_pls_val(X,Y,no_of_lv,prepro_method,val_method,segments)
%
%   --------------------------------INPUT-------------------------------------------------
%   X               Matrix (n,p), predictor matrix (assumed to be center)
%   Y               Matrix (n,m), predictand (Y is autoscaled)
%   no_of_lv        the maximum number of PCA or PLS components
%   prepro_method   'mean', 'auto', 'mscmean' or 'mscauto' (msc is performed in each interval)
%   interval        number of intervals. If intervals is a row vector divisions are made based on the elements
%   xaxislabels     (self explainable)
%   val_method      'test', 'full', 'syst111', 'syst123', 'random', or 'manual'; the last five are cross validation based methods
%   segments        number of segments in cross validation
%                   - if val_method is 'test' then segments should be a column vector with test set indices
%                   - if val_method is 'manual' then segments should be a cell array
%
%   --------------------------------OUTPUT------------------------------------------------
%   PLSmodel        a structured array containing all model and validation information
%
%   --------------------------------Reference---------------------------------------------
%   This rotine was adapted from the code from Chemometrics Group 
%   KVL, Copenhagen, Denmark (Lars Nørgaard, July 2004)
%
%   This is a part of the GNAT
%   Copyright  2021  <Mathias Nilsson>%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
%
%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%
%   Hugo da Silva Rocha, PhD Student
%   School of Chemistry, University of Manchester,
%   hugo.rocha@postgrad.manchester.ac.uk
%%

% Error checks
if ~isstruct(Model)
    uiwait(msgbox('Model input should be a structure array (output from ipls.m)'))
    return
end

if ~strcmp(Model.type,'iPLS')
    uiwait(msgbox('The model input is not from an iPLS model'))
    return
end

plottype_cell=cell(2+Model.no_of_lv,1);
for i=1:Model.no_of_lv
    plottype_cell{i+2}=num2str(i);
end
plottype_cell{1}='Cum';
plottype_cell{2}='Cum2';

if nargin<5
    plottype='Cum';    
end

if nargin==5 && ~ismember(plottype,plottype_cell)
    uiwait(msgbox('Not legal plottype, use ''Cum'', ''Cum2'' or different number of PLS components'))
    return
end

if ~ismember(labeltype,{'intlabel','varlabel','wavlabel'})
    uiwait(msgbox('Not legal labeltype, use ''intlabel'',''varlabel'', or ''wavlabel'''))
    return
end

if Model.intervalsequi==0 && strcmp(labeltype,'intlabel')
   disp(' ')
   disp(' Manually chosen intervals are not correctly plotted with ''intlabel''')
   disp(' so please use ''varlabel'' or ''wavlabel'' as labeltype')
   disp(' ')
end
% End error checks

Xmean = mean(Model.rawX); % Mean spectrum
if min(Xmean)<0
   Xmean=Xmean+(-min(Xmean)); % To make all intensities positive
end
[~,m] = size(Model.rawX);

No_Int = Model.intervals;
if strcmp(labeltype,'intlabel')
   Xint =[Model.allint(1:No_Int,1)-0.5 Model.allint(1:No_Int,1)-0.5 Model.allint(1:No_Int,1)+0.5 Model.allint(1:No_Int,1)+0.5]';
   NumberofTicksandWhere=mean(Xint(2:3,:));
elseif strcmp(labeltype,'wavlabel')
   a = Model.allint(1:No_Int,2);
   b = Model.allint(1:No_Int,3);

   % To reverse wavenumber axis before plotting; will be reversed back when the final plot is made
   NewAxisLabels=Model.xaxislabels; % Important; original axislabels are used in the last three lines of the program
   if NewAxisLabels(1)>NewAxisLabels(2)
       if size(NewAxisLabels,1)==1
           NewAxisLabels=fliplr(NewAxisLabels);
       elseif size(NewAxisLabels,2)==1
           NewAxisLabels=flipud(NewAxisLabels);
       end
   end
   
   Xint = [NewAxisLabels(a)' NewAxisLabels(a)' NewAxisLabels(b)' NewAxisLabels(b)']';
   NumberofTicksandWhere=[Xint(2,:) Xint(3,end)];

elseif strcmp(labeltype,'varlabel')
   Xtext = sprintf('Variable number');
   Xint  = [Model.allint(1:No_Int,2) Model.allint(1:No_Int,2) Model.allint(1:No_Int,3) Model.allint(1:No_Int,3)]';
   NumberofTicksandWhere=[Xint(2,:) Xint(3,end)];
end

for i=1:(Model.intervals+1)
   RMSE(i,:) = Model.PLSmodel{i}.RMSE;
end

switch plottype
    case {'Cum','Cum2'}
		% [minRMSE,min_ix]=min(RMSE'); % Global minima; could be changed using e.g. F-test or equal
		
		% First local minima is better; could be changed using e.g. F-test or equal NOT IMPLEMENTED
		% Ones appended to make the search stop if the first local miminum is the last PLSC
		RedRMSE=RMSE(:,2:end); % PLSC0 is excluded in finding the first local minimum
		SignMat=[sign(diff(RedRMSE')') ones(No_Int+1,1)];
		
		for i=1:size(RedRMSE,1)
            for j=1:size(RedRMSE,2)
                if SignMat(i,j)==1
                    min_ix(i)=j; % Note: PLSC0 is excluded
                    minRMSE(i)=RedRMSE(i,j); % Note: PLSC0 is excluded
                    break
                end
            end
		end
    case plottype_cell(3:end)
		RedRMSE=RMSE(:,2:end); % PLSC0 is excluded
		for i=1:size(RedRMSE,1)
           min_ix(i)=str2num(plottype); % Same number of components for all intervals
           minRMSE(i)=RedRMSE(i,min_ix(i)); % Note: PLSC0 is excluded
        end
end % switch
        
if nargin==2
  optimal_lv_global=min_ix(Model.intervals+1);
end

if nargin >2 && isempty(optimal_lv_global)
  optimal_lv_global=min_ix(Model.intervals+1);
end

set(0,'Units','pixels');
Scrsiz=get(0,'ScreenSize');
ScrLength=Scrsiz(3);
ScrHight=Scrsiz(4);
bdwidth=10;
% [left(->) bottom(up) width hight]
figpos=[0.1*ScrLength 0.15*ScrHight 0.85*ScrLength 0.75*ScrHight];
% figure
% set(1,'Position',figpos)
Response = [zeros(No_Int,1) minRMSE(1:No_Int)' minRMSE(1:No_Int)' zeros(No_Int,1)]';

for i=1:No_Int
    if min_ix(i)==1
        minRMSEminusOne(i)=NaN;
    else
        minRMSEminusOne(i)=RMSE(i,min_ix(i));
    end
end
ResponseMinusOnePC = [zeros(No_Int,1) minRMSEminusOne' minRMSEminusOne' zeros(No_Int,1)]';

for i=1:No_Int
    if min_ix(i)<=2
        minRMSEminusTwo(i)=NaN;
    else
        minRMSEminusTwo(i)=RMSE(i,min_ix(i)-1);
    end
end
ResponseMinusTwoPC = [zeros(No_Int,1) minRMSEminusTwo' minRMSEminusTwo' zeros(No_Int,1)]';

%% Cumulated plots
% if nargin==5 & strcmp(plottype,'Cum2')
%     % fill(Xint(:),Response(:),'c',Xint(:),ResponseMinusOnePC(:),'c',Xint(:),ResponseMinusTwoPC(:),'c') % DOESN'T WORK!!
%     if strcmp(labeltype,'wavlabel') & (Model.xaxislabels(1)>Model.xaxislabels(2))
%         fill(flipud(Xint(:)),ResponseMinusTwoPC(:),'w',flipud(Xint(:)),ResponseMinusOnePC(:),'w', flipud(Xint(:)),Response(:),[0.75 0.75 0.75])
%     else
%         fill(Xint(:),ResponseMinusTwoPC(:),'w',Xint(:),ResponseMinusOnePC(:),'w', Xint(:),Response(:),[0.75 0.75 0.75])
%     end
% elseif nargin <5 | ismember(plottype,plottype_cell([1 3:end]))
%     if strcmp(labeltype,'wavlabel') & (Model.xaxislabels(1)>Model.xaxislabels(2))
%         fill(flipud(Xint(:)),Response(:),[0.75 0.75 0.75])
%     else
%         fill(Xint(:),Response(:),[0.75 0.75 0.75]) % Note: substitute [0.75 0.75 0.75]'c' for cyan
%     end
% else
%     % Nothing
% end
% 
% if strcmp(Model.val_method,'test')
%     plottitle = sprintf('Dotted line is RMSEP (%g LV''s) for global model / Italic numbers are optimal LVs in interval model',optimal_lv_global);
%     ylabel('RMSEP','FontSize',10)
% elseif strcmp(Model.val_method,'none')
%     plottitle = sprintf('Dotted line is RMSEC (%g LV''s) for global model / Italic numbers are optimal LVs in interval model',optimal_lv_global);
%     ylabel('RMSEC','FontSize',10)
% else
%     plottitle = sprintf('Dotted line is RMSECV (%g LV''s) for global model / Italic numbers are optimal LVs in interval model',optimal_lv_global);
%     ylabel('RMSECV','FontSize',10)
% end
% title(plottitle,'FontSize',10,'FontWeight','Bold')
% % xlabel(Xtext)
% 
% hold on
%     axis tight;
%     horzline(RMSE(No_Int+1,optimal_lv_global+1),':k')
%     actualaxis=axis;
%     if nargin >= 4 & ~isempty(max_yaxis)
%         axis([actualaxis(1) actualaxis(2) actualaxis(3) max_yaxis]);
%         actualaxis(4)=max_yaxis;
%     end
%     Xaxis = linspace(actualaxis(1),actualaxis(2),m);
%     if strcmp(labeltype,'wavlabel') & (Model.xaxislabels(1)>Model.xaxislabels(2))
%         plot(fliplr(Xaxis),Xmean./max(Xmean)*actualaxis(4),'-k') % Scaled spectrum
%     else
%         plot(Xaxis,Xmean./max(Xmean)*actualaxis(4),'-k') % Scaled spectrum
%     end
%     set(gca,'XTick',NumberofTicksandWhere)
%     for i=1:Model.intervals
%       	if strcmp(labeltype,'wavlabel') & (Model.xaxislabels(1)>Model.xaxislabels(2))
%             text(mean(Xint(2:3,i)),0.03*(actualaxis(4)-actualaxis(3))+actualaxis(3),int2str(min_ix(Model.intervals-(i-1))),'Color','k','FontAngle','italic');
%         else
%             text(mean(Xint(2:3,i)),0.03*(actualaxis(4)-actualaxis(3))+actualaxis(3),int2str(min_ix(i)),'Color','k','FontAngle','italic');
%         end
%     end
% hold off
% 
% % To reverse in case of reversed wavenumber axis
% if strcmp(labeltype,'wavlabel') & (Model.xaxislabels(1)>Model.xaxislabels(2))
%     set(gca,'XDir','reverse');
% end

plotipls.Xint = Xint;
plotipls.plottype_cell = plottype_cell;
plotipls.Response = Response;
plotipls.ResponseMinusOnePC = ResponseMinusOnePC;
plotipls.ResponseMinusTwoPC = ResponseMinusTwoPC;
plotipls.optimal_lv_global = optimal_lv_global;
plotipls.RMSE = RMSE;
plotipls.max_yaxis = max_yaxis;
plotipls.NumberofTicksandWhere = NumberofTicksandWhere;
plotipls.min_ix = min_ix;
plotipls.Xmean = Xmean;



function horzline(ordinate,linetype_color)    
[~,m] = size(ordinate);    
V = axis;
if ishold
  plot(V(1:2),[1 1]*ordinate,linetype_color);
else
  hold on
  plot(V(1:2),[1 1]*ordinate,linetype_color);
  hold off
end


