function [xCS,ints,ind,target] = icoshift(xT,xP,inter,n,options,Scal)
% interval Correlation Optimized shifting
%
% [xCS,ints,ind,target] = icoshift(xT,xP,inter[,n[,options[,Scal]]])
%
% Splits a spectral database into "inter" intervals and coshift each vector
% left-right to get the maximum correlation toward a reference or toward an
% average spectrum in that interval. Missing parts on the edges after
% shifting are filled with "closest" value or with "NaNs".
%
% INPUT
% xT (1 × mT)    : target vector.
%                  Use 'average' if you want to use the average spectrum as a reference
%                  Use 'median' if you want to use the median spectrum as a reference
%                  Use 'max' if you want to use for each segment the corresponding actual spectrum having max features as a reference
%                  Use 'average2' for using the average of the average multiplied for a requested number (default=3) as a reference
%                 
% xP (nP × mP)   : Matrix of sample vectors to be aligned as a sample-set
%                  towards common target xT
% inter          : definition of alignment mode
%                  'whole'         : it works on the whole spectra (no intervals).
%                  nint            : (numeric) number of many intervals.
%                  'ndata'         : (string) length of regular intervals
%                                    (remainders attached to the last).
%                  [I1s I1e,I2s...]: interval definition. ('I(n)s' interval
%                                    n start, 'I(n)e' interval n end).
%                  (refs:refe)     : shift the whole spectra according to a
%                                    reference signal(s) in the region
%                                    refs:refe (in sampling points)
%                  'refs-refe'     : shift the whole spectra according to a
%                                    reference signal(s) in the region
%                                    refs-refe (in Scal units)
% n (1 × 1)      : (optional)
%                  n = integer n.: maximum shift correction in data
%                                  points/Scal units (cf. options(5))
%                                  in x/rows. It must be >0
%                  n = 'b' (best): the algorithm search for the best n
%                                  for each interval (it can be time consuming!)
%                  n = 'f' (fast): fast search for the best n for each interval (default)
%                  A warning is displayed for each interval if "n" appears too small
% options (1 × 5): (optional)
%                  (1) triggers plots & warnings:
%                      0 : no on-screen output
%                      1 : only warnings (default)
%                      2 : warnings and plots
%                  (2) selects filling mode
%                      0 : using not a number
%                      1 : using previous point (default)
%                  (3) turns on Co-shift preprocessing
%                      0 : no Co-shift preprocessing (default)
%                      1 : Executes a Co-shift step before carrying out iCOshift
%                  (4) max allowed shift for the Co-shift preprocessing (default = equal to n if not specified)
%                      it has to be given in Scal units if option(5)=1
%                  (5) 0 : intervals are given in No. of datapoints  (deafult)
%                      1 : intervals are given in ppm --> use Scal for inter and n
% Scal           : vector of scalars used as axis for plot (optional)
%
% OUTPUT
% xCS  (nP × mT): shift corrected vector or matrix
% ints (nI × 4) : defined intervals (Int. No., starting point, ending point, size)
% ind  (nP × nI): matrix of indexes reporting how many points each spectrum
%                 has been shifted for each interval (+ left, - right)
% target (1 x mP): actual target used for the final alignment
%
% Authors:
% Francesco Savorani - Department of Food Science
%                      Quality & Technology - Spectroscopy and Chemometrics group
%                      Faculty of Sciences
%                      University of Copenhagen - Denmark
% email: frsa@life.ku.dk - www.models.life.ku.dk
%
% Giorgio Tomasi -     Department of Basic Science and Environment
%                      Soil and Environmental Chemistry group
%                      Faculty of Life Sciences
%                      University of Copenhagen - Denmark
% email: giorgio.tomasi@ec.europa.eu - www.igm.life.ku.dk

% 170508 (FrSa) first working code
% 211008 (FrSa) improvements and bugs correction
% 111108 (Frsa) Splitting into regular intervals (number of intervals or wideness in datapoints) implemented
% 141108 (GT)   FFT alignment function implemented
% 171108 (FrSa) options implemented
% 241108 (FrSa) Automatic search for the best or the fastest n for each interval implemented
% 261108 (FrSa) Plots improved
% 021208 (FrSa) 'whole' case added & robustness improved
% 050309 (GT)   Implentation of interpolation modes (NaN); Cosmetics; Graphics
% 240309 (GT)   Fixed bug in handling missing values
% 060709 (FrSa) 'max' target and output 'target' added. Some speed, plot and robustness improvements
% 241109 (GT)   Interval and band definition in units (added options(5))
% 021209 (GT)   Minor debugging for the handling of options(5)
% 151209 (FrSa) Cosmetics and minor debugging for the handling of options(5)
% 151110 (FrSa) Option 'Max' works now also for alignment towards a reference signal
% 310311 (FrSa) Bugfix for the 'whole' case when mP < 101
% 030712 (FrSa) Introducing the 'average2' xT (target) for a better automatic target definition. Plots updated to include also this case
%% ERRORS CHECK
% Constant
BLOCKSIZE = 2^25; % To avoid out of memory errors when 'whole', the job is divided in blocks of 32MB

% Number of input arguments
if ~nargin
   
   fprintf('\n\n icoshift(target,dataset,alignment mode,max allowed shift,options,scale);')
   fprintf('\n or')
   fprintf('\n icoshift(target,dataset,alignment mode);\n')
   fprintf('\n Examples:')
   fprintf('\n icoshift(''median'',xP,20);\n')
   fprintf('splits in 20 regular intervals using the default settings and align them towards the median spectrum\n')
   fprintf('\n icoshift(''average'',xP,''200'',''40'');\n')
   fprintf('splits in regular intervals made of 200 data points using the default settings\n')
   fprintf('and align them towards the average spectrum allowing a 40 points max shift\n')
   fprintf('\n icoshift(target,xP,(1000:2000),''b'',[2 1 1]);\n')
   fprintf('aligns the whole spectra according to the reference signal in target(1000:2000)\n')
   fprintf('searching for the best max allowed shift for each interval and performing a pre-coshift step.\n')
   fprintf('Plots the achieved results\n')
   fprintf('\n icoshift(''average'',xP,user_ints,''f'',[2 0],ppm_scale);\n')
   fprintf('splits the dataset into the intervals defined by the user (user_ints) and aligns them\n')
   fprintf('towards the average spectrum using the highest allowed shift (fastest method) for each interval\n')
   fprintf('and filling the interval edges with NaNs. Plots the achieved results using the ppm_scale provided\n')
   fprintf('\n For further specifications type ''help icoshift''\n')
   return
   
end
if nargin < 3, error('Not enough input arguments'); end
if nargin < 4 || isempty(n)
   n = 'f'; % Change here the default assigned "n" value
   fprintf('\n Fast search for the best "n" set by default (initial value 50)\n')
end
if nargin < 5, options = []; end
if nargin < 6, Scal = 1:size(xP,2); end

%% Check input
% Scal #GT To handle options(5) properly
if numel(Scal) ~= length(Scal), error('Scal must be a vector'), end
if length(Scal) ~= size(xP,2), error('X and Scal are not of compatible length'), end
dScal   = diff(Scal);
incScal = Scal(1) - Scal(2);
if incScal == 0 || ~all(sign(dScal) == -sign(incScal)), error('Scal must be strictly monotonic'); end
flag_ScalDir = incScal < 0;
flag_DiScal  = any(abs(dScal) > 2 * min(abs(dScal))); % #GT Check for discontinuous Scal vector

% xT and xP
max_flag = 0;
avg2_flag = 0;
if  strcmpi(xT,'average')
    xT   = MeanwNaN(xP);
    note = 'average';
elseif strcmpi(xT,'median')
    xT = nanmedian(xP);
    note = 'median';
elseif strcmpi(xT,'average2')   %AVG2
    xT = MeanwNaN(xP);           %AVG2
    avg2_flag = 1;                        %AVG2
    if options(1) ~= 0;           %AVG2
        avg_power = input('which multiplier for the 2nd average? [default is 3] ');      %AVG2
    else avg_power = 3;      %AVG2
    end
    if isempty(avg_power); avg_power = 3;end                %AVG2
    note = 'average2';                 %AVG2
elseif strcmpi(xT,'max')        % MAX
    xT = zeros(1,size(xP,2));    % MAX
    max_flag = 1;                % MAX
    note = 'max';                % MAX
else note = 0;
end
[nT,mT] = size(xT);
[nP,mP] = size(xP);
if (mT ~= mP), error('Target "xT" and sample "xP" must have the same number of columns'); end
if isnumeric(inter) %FRSA bugfix for 'whole' case
    if inter > mP
        error('ERROR: number of intervals "inter" must be smaller than number of variables in xP');
    end
end

% Options
% options_co   = [1 1 0 0 nargin > 5];
options_co   = [1 1 0 0 0]; % FrSa: default options(5) has to be = 0
if length(options) < length(options_co), options(end + 1:length(options_co)) = options_co(length(options) + 1:end); end
if options(1) < 0 || options(1) > 2,   error('options(1) must be 0, 1 or 2'); end
if options(5), 
   prec = abs(min(unique(dScal))); 
   if flag_DiScal, 
      warning('iCOshift:discontinuous_Scal','Scal vector is not continuous, the defined intervals might not reflect actual ranges'); 
   end
end

% Execute preliminary coshift if so required
flag_coshift = ~strcmpi(inter,'whole') && options(3);
if flag_coshift
   
   if options(4) == 0
      n_co = n;
   else
      n_co = options(4);
   if nargin >= 6 && options(5), n_co = dscal2dpts(n_co,Scal,prec); end % GT To handle scalar input
   end
   if max_flag             % MAX
      xT = MeanwNaN(xP);   % MAX
   end                     % MAX
   xPo = xP;    %FRSA save the original xP for plotting purpose
   [xP,nil,wint] = icoshift(xT,xP,'whole',n_co,[0 1 0]); %#ok<*ASGLU> %FRSA preprocessing COshift doesn't use NaNs
   title_suff    = ' (after coshift)';
   if  strcmpi(note,'average')
      xT = nanmean(xP); %FRSA uses the new CO-shifted average as a target
   end
   if  strcmpi(xT,'median')
      xT = nanmedian(xP); %FRSA uses the new CO-shifted median as a target
   end
   if  strcmpi(note,'average2')
      xT = nanmean(xP); %FRSA uses the new CO-shifted average as a first iteration average target
   end
   
else
   title_suff = '';
end

% "Inter" % GT: handle scalars automatically
flag_allsegs = false;
whole        = false;
flag2        = false;
if ischar(inter)
   
   if  strcmpi(inter,'whole')
      inter = [1,mP];
      whole = true;
   else
      
      if ~any(inter == '-')
      
         interv = str2double(inter);
         if nargin < 6 || ~options(5), interv = round(interv);
         else interv = dscal2dpts(interv,Scal,prec);
         end
         [inter]      = defints(xP,interv,options(1));
         flag_allsegs = true;
      
      else
         
         interv = regexp(inter,'(-{0,1}\d*\.{0,1}\d*)-(-{0,1}\d*\.{0,1}\d*)','tokens');
         interv = sort(scal2pts(str2double(cat(1,interv{:})),Scal,prec));
         if numel(interv) ~= 2, error('Invalid range for reference signal'), end
         inter = interv(1):interv(2);
         flag2 = true;
         
      end
      
   end
   
elseif isa(inter,'numeric')
   
   if numel(inter) == 1
      
      if fix(inter) == inter
         
         flag_allsegs = true;
         remain       = mod(mP,inter);
         N            = fix(mP/inter);
         % Distributes vars_left_over in the first "vars_left_over" intervals
         startint     = [(1:(N+1):(remain - 1)*(N+1)+1)'; ((remain - 1) * (N + 1) + 1 + 1 + N:N:mP)'];
         % endint=[startint(2:inter)-1; mP]; %ints with adjacent edges
         endint       = [startint(2:inter); mP]; % ints with common edges
         inter        = [startint endint]';
         inter        = inter(:)';
         
      else error('The number of intervals must be an integer');
      end
      % else intdif = diff(inter);
      %     if any(intdif(2:2:length(intdif)) < 0)
      %         uiwait(msgbox('The user-defined intervals are overlapping: is it intentional?','Warning','warn'));
      %     end
      
   else      
      flag2 = isequal(fix(inter),inter) && length(inter) > 1 && ...
         isequal([1,max(inter) - min(inter) + 1],size(inter)) && isequal(unique(diff(inter,1,2)),1);
      if ~flag2 && options(5), 
         inter = scal2pts(inter,Scal,prec);
         if any(inter(1:2:end) > inter(2:2:end)) && ~flag_ScalDir 
            inter = flipud(reshape(inter,2,length(inter)/2));
            inter = inter(:)';
         end
      end % Converts scalars to points 
   end 
   
end

[nint,mint] = size(inter);
wmsg = [];
scfl = isequal(fix(Scal),Scal) & ~options(5);
if ischar(n) && (numel(n) > 1 || ~any(ismember('bf',lower(n))))
   error('"n" must be a scalar/''f''/''b''')
elseif isa(n,'numeric')
   
   if any(n <= 0),  error('ERROR: shift(s) "n" must be larger than zero'); end
   if numel(n) > 1, wmsg = sprintf('"n" must be a scalar/character; first element (i.e. %i) used',round(n)); end
   if scfl && n ~= fix(n),  wmsg = sprintf('"n" must be an integer if Scal is ignorde; first element (i.e. %i) used',round(n));
   elseif options(5), n = dscal2dpts(n,Scal,prec);
   end
   if ~flag2 && any(diff((reshape(inter,2,mint/2)),1,1) < n)
      error('ERROR: shift "n" must be not larger than the size of the smallest interval');
   end
   n = round(n(1));
   if ~isempty(wmsg)
      warning('iCoShift:Input',wmsg)
      fprintf('\n press a key to continue...')
      pause
   end
   
end

%% Missing values
flag     = isnan(cat(1,xT,xP));
frag     = false;
Ref      = @(c) reshape(c,2,length(c)/2)';
vec      = @(A) A(:);
[mi,pmi] = min(inter);
[ma,pma] = max(inter);
if any(vec(flag))
   
   % If there are missing values in the target, check that the same pattern
   % of missing is present in the samples. If not an error is produced.
   if isequal(flag(ones(nP,1),:),flag(2:end,:))
      Select = @any;
   else
      Select = @all;
   end
   if flag2
      
      intern = RemoveNaN([1,pma - pmi + 1],cat(1,xT(:,inter),xP(:,inter)),Select);
      if size(intern,1) ~= 1, error('Reference region contains a pattern of missing that cannot be handled consistently')
      elseif ~isequal(intern,[1,inter(end) - inter(1) + 1]), warning('iCoShift:miss_refreg','The missing values at the boundaries of the reference region will be ignored')
      end
      intern = inter(1) + intern(1) - 1:inter(1) + intern(2) - 1;
      
   else
      [intern,flag_nan] = RemoveNaN(Ref(inter),cat(1,xT,xP),Select);
      intern            = vec(intern')';
   end
   if isempty(intern), error('iCoShift cannot handle this pattern of missing values.'), end
   if length(intern) ~= length(inter) && ~flag2
      
      if whole
         
         if length(intern) > 2
            
            if options(1) == 2
               xPBU = xP;
               xTBU = xT;
            end
            [XSeg,InOr] = ExtractSegments(cat(1,xT,xP),Ref(intern));
            inter       = [InOr(1),InOr(end)];
            InOr        = cat(2,Ref(intern),InOr);
            xP          = XSeg(2:end,:);
            xT          = XSeg(1,:);
            frag        = true;
            
         end
         
      else
         warning('iCoShift:Missing_values','\nTo handle the pattern of missing values, %i segments are created/removed',abs(length(intern) - length(inter))/2)
         inter       = intern;
         [nint,mint] = size(inter);
      end
      
   end
   
end

%% INITIALISE OUTPUT
xCS = xP;

%% ALIGN
[mi,pmi] = min(inter);
[ma,pma] = max(inter);
flag     = length(inter) > 1 && isequal([1,pma - pmi + 1],size(inter)) && isequal(unique(diff(inter,1,2)),1);
if flag % a Reference signal has been chosen
   
   if options(1)
      
      if  strcmpi(n,'b')      %     Option 'best search'
         fprintf('\n Automatic searching for the best "n" for the reference window "refW" enabled \n That can take a longer time \n')
      elseif  strcmpi(n,'f')  %     Option 'fast search'
         fprintf('\n Fast automatic searching for the best "n" for the reference window "refW" enabled \n')
      end
      
   end
   if max_flag                                  % MAX
         [amax,bmax] = max(sum(xP(:,mi:ma),2)); % MAX
         xT(mi:ma) = xP(bmax,mi:ma);            % MAX
   end
   ind                         = NaN(nP,1);
   missind                     = ~all(isnan(xP),2);
   [xCS(missind,:),ind(missind)] = coshifta(xT,xP(missind,:),inter,n,[options(1:3),BLOCKSIZE]);
   ints                        = [1 mi ma];
   
else % splits the dataset in intervals
   
   % Intervals definition
   if mint > 1
      if rem(mint,2), error('Wrong definition of intervals ("inter")'); end
      if ma>mP, error('Intervals ("inter") exceed samples matrix dimension'); end
      allint=[(1:round(mint/2))' inter(1:2:mint)' inter(2:2:mint)'];
   end
   %  Cheks interval overlapping
   sallint = sortrows(allint,2);
   sinter = reshape(sallint(:,2:3)',1,numel(sallint(:,2:3)));
   intdif = diff(sinter);
   if any(intdif(2:2:length(intdif)) < 0)
      uiwait(msgbox('The user-defined intervals are overlapping: is that intentional?','Warning','warn'));
   end
   clear sallint; clear sinter; clear intdif
   
   ints = allint;
   ints(:,4)=(ints(:,3)-ints(:,2))+1;
   ind  = zeros(nP,size(allint,1));
   % end Interval definition
   
   if options(1)
      
      if  strcmpi(n,'b')
         fprintf('\n Automatic searching for the best "n" for each interval enabled \n That can take a longer time \n')
      elseif  strcmpi(n,'f')
         fprintf('\n Fast automatic searching for the best "n" for each interval enabled \n' )
      end
      
   end
   
   % Local co-shifting
   %     top_int = NaN(1,size(allint,1));
   for i = 1:size(allint,1)
      
      if options(1) ~= 0;
         if whole, fprintf('\n Co-shifting the whole %g samples... \r',nP);
         else fprintf('\n Co-shifting interval no. %g of %g... \r',i,size(allint,1));
         end
      end
      intervalnow = xP(:,allint(i,2):allint(i,3));
      if max_flag                               % MAX
         [amax,bmax] = max(sum(intervalnow,2)); % MAX
         target  = intervalnow(bmax,:);         % MAX
         xT(allint(i,2):allint(i,3)) = target;  % MAX
      else                                      % MAX
         target  = xT(:,allint(i,2):allint(i,3)); 
      end
      missind     = ~all(isnan(intervalnow),2);
      if ~all(isnan(target)) && sum(missind) ~= 0
         [cosh_interval,loc_ind]              = coshifta(target,intervalnow(missind,:),0,n,[options(1:3),BLOCKSIZE]);
         xCS(missind,allint(i,2):allint(i,3)) = cosh_interval;
         ind(missind,i)                       = loc_ind';
      else
         xCS(:,allint(i,2):allint(i,3)) = intervalnow;
      end
      %         top_int(i) = max(max(intervalnow));
      home
      
   end
   % AVG2 secong iteration (FrSa)
   if avg2_flag
       for i = 1:size(allint,1)
           if options(1) ~= 0;
               if whole, fprintf('\n Co-shifting again the whole %g samples... \r',nP);
               else fprintf('\n Co-shifting again interval no. %g of %g... \r',i,size(allint,1));
               end
           end
%            intervalnow = xCS(:,allint(i,2):allint(i,3)); % (FrSa) further align the previously aligned matrix
           intervalnow = xP(:,allint(i,2):allint(i,3));  % align the raw matrix using the newly calculated average2 target
           target1  = mean(xCS(:,allint(i,2):allint(i,3)));
           min_interv=min(target1);
           target = (target1-min_interv).*avg_power; %AVG2 power multiplier and baseline correction
           missind     = ~all(isnan(intervalnow),2);
           if ~all(isnan(target)) && sum(missind) ~= 0
               [cosh_interval,loc_ind]              = coshifta(target,intervalnow(missind,:),0,n,[options(1:3),BLOCKSIZE]);
               xCS(missind,allint(i,2):allint(i,3)) = cosh_interval;
               xT(allint(i,2):allint(i,3)) = target;
               ind(missind,i)                       = loc_ind';
           else
               xCS(:,allint(i,2):allint(i,3)) = intervalnow;
           end
           %         top_int(i) = max(max(intervalnow));
           home
       end
   end
   
end

% If n == 'whole' and some values are missing, reconstruct signal
if frag
   
   Xn  = NaN(nP,mP);
   for i_sam = 1:nP
      
      for i_seg = 1:size(InOr,1),
         
         Xn(i_sam,InOr(i_seg,1):InOr(i_seg,2)) = xCS(i_sam,InOr(i_seg,3):InOr(i_seg,4));
         if loc_ind(i_sam) < 0
            if flag_nan(i_seg,1,i_sam), Xn(i_sam,InOr(i_seg,1):InOr(i_seg,1) - loc_ind(i_sam)) = NaN; end
         elseif loc_ind(i_sam) > 0
            if flag_nan(i_seg,2,i_sam), Xn(i_sam,InOr(i_seg,2) - loc_ind(i_sam) + 1:InOr(i_seg,2)) = NaN; end
         end
         
      end
      
   end
   if options(1) == 2 % if plots are activated
      xP  = xPBU;
      xT  = xTBU;
   end
   xCS = Xn;
   
end
target = xT;
%% PLOT COMPARISON rigid shift
if flag_coshift %Plots the actual raw data
   %     xP = xPo;
   %     clear xPo
   sp1 = 3; sp2 = 1;    %FRSA Subplot type
else sp1 = 2; sp2 = 1;   %FRSA Subplot type
end

if options(1)==2
   
   if flag_allsegs, Colour = [1 1 1];
   else Colour = [1.0000    0.8314    0.8314];
   end
   pbu = get(0,'units');
   set(0,'Units','pixels');
   Scrsiz = get(0,'ScreenSize');
   set(0,'Units',pbu);
   ScrLength = Scrsiz(3);
   ScrHight  = Scrsiz(4);
   figpos    = [0.1*ScrLength 0.15*ScrHight 0.85*ScrLength 0.75*ScrHight];
   figure('unit','pix','Position',figpos);
   st        = ceil(mP/(2^14));
   x         = 1:st:mP;
   CM        = ExtCM(size(xP,1),1 - MyCM);
   top       = max(xP(:));
   bot       = min(xP(:));
   sh        = NaN(1,2);
   
   sh(sp2)     = subplot(sp1,1,sp2,'nextplot','add','colororder',CM,'color','w');
   if flag_coshift
      ph        = plot(Scal(x),xPo(:,x)'); %line([mi ma; mi ma], [0 0; top top], 'Color', 'r', 'LineStyle', ':', 'LineWidth', 2)
   else ph        = plot(Scal(x),xP(:,x)');
   end
   ax        = axis;
   ax(4)     = ax(4)*1.1;
   if ax(3) > bot
      ax(3) = bot*2;
   end
   if round(bot) >= 0
      ax(3) = -top/15;
   end
   if flag
      plot(Scal([mi ma; mi ma]), ax(3:4),':r','LineWidth',2);
      text((Scal(mi) + Scal(ma))/2,max(max(xP(:,inter))*1.1),'Ref.','Color','k', 'Rotation', 90, 'clip','on');
   else
      
      %         xpts = [(allint(1:end,2))' (allint(1:end,3))'; (allint(1:end,2))' (allint(1:end,3))'];
      XX   = [(allint(1:end,2))'; (allint(1:end,2))';(allint(1:end,3))';(allint(1:end,3))'];
      YY   = [ax(3);ax(4);ax(4);ax(3)]*ones(1,size(allint,1));
      if ~whole
         h = fill(Scal(XX),YY,Colour);
         set(h,'facealpha',0.4,{'tag'},cellstr(num2str(allint(:,1))),'buttondown',@DispInfo,'edgecolor','b','linewidth',1)
         %          plot(xpts, ax(3:4),'b','LineWidth',1);
         text((Scal(allint(:,2)) + Scal(allint(:,3)))/2,ones(size(allint,1),1)*top,num2str(allint(:,1)), 'Rotation', 90,'clip','on');
      end
      
   end
   axis tight
   title('Raw data');
   delete(ph)
   if flag_coshift
      ph1 = plot(Scal(x),xPo(:,x)'); %line([mi ma; mi ma], [0 0; top top], 'Color', 'r', 'LineStyle', ':', 'LineWidth', 2)
      avrg = nanmean(xPo);
      clear xPo;
   else ph1 = plot(Scal(x),xP(:,x)');
   end
   set(ph1,{'tag'},cellstr(num2str((1:nP)','Sam. #%i')),'buttondown',@DispInfo)
   
   % 2nd plot
   if flag_coshift
      sp2 = sp2+1;
      sh(sp2) = subplot(sp1,1,sp2,'nextplot','add','colororder',CM,'color','w');
      ph = plot(Scal(x),xP(:,x)');%line([mi ma; mi ma], [0 0; top top], 'Color', 'g', 'LineStyle', ':', 'LineWidth', 2)
      if flag
         plot(Scal([mi ma; mi ma]), ax(3:4),':g','LineWidth',2);
         text((Scal(mi) + Scal(ma))/2,max(max(xP(:,inter))*1.1),'Ref.','Color','k', 'Rotation', 90, 'clip','on');
         title(['Reference aligned data',title_suff]);
      else
         
         if ~whole
            
            h = fill(Scal(XX),YY,Colour);
            set(h,'facealpha',0.5,{'tag'},cellstr(num2str(allint(:,1))),'buttondown',@DispInfo,'edgecolor','b','linewidth',1)
            %          plot(xpts, ax(3:4),'b','LineWidth',1);
            if flag_coshift
               title('pre-coshifted data');
            else title(['Interval-aligned data',title_suff]);
            end
            text((Scal(allint(:,2)) + Scal(allint(:,3)))/2,ones(size(allint,1),1)*top,num2str(allint(:,1)), 'Rotation', 90,'clip','on');
            
         else title('Whole spectra aligned data');
         end
         
      end
      axis tight
      delete(ph)
      ph2 = plot(Scal(x),xP(:,x)');%line([mi ma; mi ma], [0 0; top top], 'Color', 'g', 'LineStyle', ':', 'LineWidth', 2)
      set(ph2,{'tag'},cellstr(num2str((1:nP)','Sam. #%i')),'buttondown',@DispInfo)
      linkaxes(sh,'xy')
   end
   
   % 2nd or 3rd plot
   sp2 = sp2+1;
   sh(sp2) = subplot(sp1,1,sp2,'nextplot','add','colororder',CM,'color','w');
   ph = plot(Scal(x),xCS(:,x)');%line([mi ma; mi ma], [0 0; top top], 'Color', 'g', 'LineStyle', ':', 'LineWidth', 2)
   if flag
      plot(Scal([mi ma; mi ma]), ax(3:4),':g','LineWidth',2);
      text((Scal(mi) + Scal(ma))/2,max(max(xP(:,inter))*1.1),'Ref.','Color','k', 'Rotation', 90, 'clip','on');
      title(['Reference aligned data',title_suff]);
   else
      
      if ~whole
         
         h = fill(Scal(XX),YY,Colour);
         set(h,'facealpha',0.5,{'tag'},cellstr(num2str(allint(:,1))),'buttondown',@DispInfo,'edgecolor','b','linewidth',1)
         %          plot(xpts, ax(3:4),'b','LineWidth',1);
         title(['Interval-aligned data',title_suff]);
         text((Scal(allint(:,2)) + Scal(allint(:,3)))/2,ones(size(allint,1),1)*top,num2str(allint(:,1)), 'Rotation', 90,'clip','on');
         
      else title('Whole spectra aligned data');
      end
      
   end
   axis tight
   delete(ph)
   ph3 = plot(Scal(x),xCS(:,x)');%line([mi ma; mi ma], [0 0; top top], 'Color', 'g', 'LineStyle', ':', 'LineWidth', 2)
   set(ph3,{'tag'},cellstr(num2str((1:nP)','Sam. #%i')),'buttondown',@DispInfo)
   linkaxes(sh,'xy')
   
   if flag_coshift
      subplot(sp1,1,1)
      h = plot (Scal(x),avrg(x), 'Color','r', 'LineWidth', 4);
      legend (h, 'average');
      subplot(sp1,1,2)
      h = plot (Scal(x),xT(x), 'Color','r', 'LineWidth', 4);
      if avg2_flag %AVG2
          note = ['average2 * ', num2str(avg_power)];                 %AVG2
      end
      legend (h, note);
   else
      if any(strcmpi(note,{'average','median', 'max', 'average2'})) % MAX
         subplot(sp1,1,1)
         h = plot (Scal(x),xT(x), 'Color','r', 'LineWidth', 4);
         if avg2_flag %AVG2
             note = ['average2 * ', num2str(avg_power)];                 %AVG2
         end
         legend (h, note);
         a = sum(abs(ind),2);
         b = find(a == min(a));
         result = ['Sample(s) No. [' num2str(b') '] is/are the most similar to the ' note ' one'];
         msgbox(result)
      else subplot(sp1,1,1)
         h = plot (Scal(x),xT(x), 'Color','r', 'LineWidth', 4);
         legend (h, 'Target');
      end
   end
   
   if ~flag_ScalDir, set(sh,'xdir','rev'), end % Reverts X axis direction if Scal is monotonically decreasing
   for i_sam = 1:nP
      if flag_coshift
         setappdata(ph1(i_sam),'Extra',ph3(i_sam))
         setappdata(ph2(i_sam),'Extra',ph1(i_sam))
         setappdata(ph3(i_sam),'Extra',ph2(i_sam))
      else
         setappdata(ph1(i_sam),'Extra',ph3(i_sam))
         setappdata(ph3(i_sam),'Extra',ph1(i_sam))
      end
      setappdata(ph1(i_sam),'flag_AvoidEndlessRecursion',1)
      setappdata(ph3(i_sam),'flag_AvoidEndlessRecursion',1)
      if flag_coshift
         setappdata(ph2(i_sam),'flag_AvoidEndlessRecursion',1)
      end
   end
   
end
if flag_coshift,
   ind = ind + wint * ones(1,size(ind,2)); end

function [xW,ind,R] = coshifta(xT,xP,refW,n,options)
% function [xW,ind,R] = coshifta(xT,xP,refW,n,options)
% Based on Correlation optimized shifting (COshift)
% Author:
% Frans van den Berg* (FvdB) email: fb@life.ku.dk
% Improved by:
% Francesco Savorani* (FS)   email: frsa@life.ku.dk
% Giorgio Tomasi*§    (GT)   email: gto@life.ku.dk
%
% * Department of Food Science Quality & Technology - Spectroscopy and Chemometrics group - www.models.life.ku.dk
% § Department of Basic Science and Environment - Soil and Environmental Chemistry group - www.igm.life.ku.dk
% Faculty of Life Sciences - University of Copenhagen - Denmark
%
% loosely based on : V.G. van Mispelaar et.al. 'Quantitative analysis of target components by comprehensive two-dimensional gas chromatography'
%                   J. Chromatogr. A 1019(2003)15-29

% HISTORY
% 1.00.00 2008 Jun 19 First documented version
% 2.00.00 2009 Nov 24 Blocking in case of large X (32MB blocks by default - options(4))


% ERRORS CHECK
if (isempty(refW)), refW = 0; end
if all(refW >= 0), rw = length(refW);
else rw = 1;
end
options_def = [ones(1,3) 2^25];
options(end + 1:length(options_def)) = options_def(length(options) + 1:end);
if options(2) == 1, Filling = -inf;
elseif options(2) == 0, Filling = NaN;
else error('options(2) must be 0 or 1');
end

if  strcmpi(xT,'average'), xT = nanmean(xP); end
[nT,mT] = size(xT);
[nP,mP] = size(xP);
[nR,mR] = size(refW);

if (mT ~= mP), error('Target "xT" and sample "xP" must be of compatible size (same vectors, same matrices or row + matrix of rows)'); end
if any(n <= 0),error('Shift(s) "n" must be larger than zero'); end
if (nR ~= 1), error('Reference windows "refW" must be either a single vector or 0'); end
if rw > 1 && (min(refW)<1) || (max(refW)>mT), error('Reference window "refW" must be a subset of xP'); end
if (nT ~= 1), error('Target "xT" must be a single row spectrum/chromatogram'); end

auto = 0;
if  strcmpi(n,'b')
   auto     = 1;                      % switch for the best automatic serach on
   if rw ~=1
      n        = max(fix(0.05 * mR),10); % change here the first searching point for the best "n"
      src_step = fix(mR * 0.05);         % Change here to define the searching step
   else
      n        = max(fix(0.05 * mP),10); % change here the first searching point for the best "n"
      src_step = fix(mP * 0.05);         % Change here to define the searching step
   end
   try_last = 0;                          %FRSA
elseif strcmpi(n,'f')
   
   auto = 1; % switch for the fast automatic search on %FRSA modified from 1 to 0
   if rw ~=1
      n        = mR-1; %floor(mR/2);     %FRSA speed up fast option
               src_step = round(mR/2)-1;        %FRSA
   else
      n        = mP-1; %floor(mP/2);     %FRSA
               src_step = round(mP/2)-1;        %FRSA
   end
   try_last = 1;                          %FRSA speed up fast option
end
if (nT ~= 1), error('ERROR: Target "xT" must be a single row spectrum/chromatogram'); end
% End ERRORS CHECK

xW                             = NaN(nP,mP); % GT # Insert NaNs
ind                            = zeros(1,nP);
BLOCKSIZE                      = options(4);
nBlocks                        = whos('xP');
nBlocks                        = ceil(nBlocks.bytes/BLOCKSIZE);
SamxBlock                      = floor(nP/nBlocks);
indBlocks                      = SamxBlock(1,ones(1,nBlocks));
indBlocks(1:rem(nP,SamxBlock)) = SamxBlock + 1;
indBlocks                      = [0,cumsum(indBlocks)];
   
% try_last = 0;

if auto == 1
   
   while auto == 1 % Automatic search for the best "n" for each interval
      
      if Filling == -inf,    xtemp = [xP(:,ones(1,n)) xP xP(:,mP(1,ones(1,n)))];
      elseif isnan(Filling), xtemp = [NaN(nP,n),xP,NaN(nP,n)];
      end
      if rw == 1, refW = 1:mP; end
      ind = NaN(nP,1);
      R   = [];
      for i_block = 1:nBlocks
         block_indices = indBlocks(i_block) + 1:indBlocks(i_block + 1);
         [dummy,ind(block_indices),Ri] = ... %FFT Co-Shifting
            CC_FFTShift(xT(1,refW),xP(block_indices,refW),[-n n 2 1 Filling]); %#ok<ASGLU> 
            R = cat(1,R,Ri);
      end
      temp_index = -n:n;
      for i_sam = 1:nP
         index       = find(temp_index == ind(i_sam));
         xW(i_sam,:) = xtemp(i_sam,index:index + mP - 1);
      end
      if (max(abs(ind)) == n) && try_last ~= 1 %to check
         if n+src_step >= size(refW,2)
             try_last = 1; continue %FRSA to avoid n > interval size
         end
         n = n+src_step ;
         continue
      elseif (max(abs(ind)) < n) && n+src_step < size(refW,2) && try_last ~= 1
         n = n+src_step ;
         try_last = 1;
         continue
      else
         auto = 0;
         if options(1) ~= 0
            fprintf('\n Best shift allowed for this interval = %g \n',n);
         end
      end
      
   end
   
else % Uses the same defined "n" for each interval
   
   if Filling == -inf,    xtemp = [xP(:,ones(1,n)) xP xP(:,mP(1,ones(1,n)))];
   elseif isnan(Filling), xtemp = [NaN(nP,n),xP,NaN(nP,n)];
   end
   if rw == 1, refW = 1:mP; end
   ind = NaN(nP,1);
   R   = [];
   for i_block = 1:nBlocks
      block_indices = indBlocks(i_block) + 1:indBlocks(i_block + 1);
      [dummy,ind(block_indices),Ri] = ... %FFT Co-Shifting
         CC_FFTShift(xT(1,refW),xP(block_indices,refW),[-n n 2 1 Filling]); %#ok<ASGLU>
      R = cat(1,R,Ri);
   end
   temp_index = -n:n;
   for i_sam = 1:nP
      index       = find(temp_index == ind(i_sam));
      xW(i_sam,:) = xtemp(i_sam,index:index + mP - 1);
   end
   if (max(abs(ind)) == n) && options(1) ~= 0
      disp('Warning: Scrolling window size "n" may not be enough wide because extreme limit has been reached')
   end
   
end

function [inter] = defints(xP,interv,opt)
% Defines a vector "inter" [startint1 endint1 startint2 endint2...startintN endintN]
% for splitting a dataset "xP" into intervals having an imposed wideness of "interv" datapoints.
%
% SYNTAX
% [inter] = defints(xP,interv,opt)
%
% INPUT
% Xp (n x m): matrix to be split in regular intervals "interv" points wide
% interv (1 x 1): number of datapoints defining intervals' wideness
% opt (1 x 1): triggers warnings on/off
%
% OUTPUT
% inter [startint1 endint1 startint2 endint2...startintN endintN]: vector
% defining intervals boundaries
%
% Author: Francesco Savorani
%         frsa@life.ku.dk
% 11/11/08 (FrSa)
%
[nP,mP] = size(xP);
sizechk = mP/interv-round(mP/interv);
plus = (mP/interv-round(mP/interv))*interv;
mbx=['Warning: the last interval will not fulfill the selected intervals size "inter" = ', num2str(interv)];
if plus >= 0
   mbx2=['Size of the last interval = ', num2str(plus)];
else
   mbx2=['Size of the last interval = ', num2str(interv+plus)];
end
mbx3={mbx;mbx2};
if opt(1)==2 && (sizechk ~= 0)
   uiwait(msgbox(mbx3,'Warning','warn'));
elseif opt(1)~=0 && (sizechk ~= 0)
   fprintf('\r Warning: the last interval will not fulfill the selected intervals size "inter"=%g. \r Size of the last interval = %g \r',interv,plus);
end
t = cat(2,0:interv:mP,mP);
if t(end) == t(end - 1), t(end) = []; end
t     = cat(1,t(1:end - 1) + 1,t(2:end));
inter = t(:)';

function [Xwarp,Shift,Values] = CC_FFTShift(T,X,Options)
% Calculate optimal rigid shift using Fast Fourier Transform
% Signals are expected to be on the rows (horizontal slabs).
%
% SYNTAX
% [Xwarp,Shift,Values] = CC_FFTShift(T,X,Options)
%
% INPUT
% T      : target
% X      : signals
% Options: 1 x 4 vector. If element is NaN, the default value is used
%          (1) lower bound for shift
%          (2) upper bound for shift
%          (3) mode along which shift is to be determined (the last by default)
%          (4) 0 -> find maximum mean of cross-correlation for multivariate signals
%              1 -> find maximum product of cross-correlation for multivariate signals (default)
%          (5) fill in value (default: NaN), -inf uses the previous/following point
%
% OUTPUT
% Xwarp  : aligned signals
% Shifts : (i) shift for the i-th signal
% Values : (i + 1,j) cost function for the i-th signal and shift (j),
%          the first row contains the investigated shifts
%
% Author: Giorgio Tomasi
%         giorgio.tomasi@gmail.com
%
% Created      : August 2006
% Last modified: 24 March, 2009; 01:15

% HISTORY
% 1.00.00 August 2006    -> First working version
% 1.01.00 23 March, 2009 -> Fixed problem with single signals
% 2.00.00 23 March, 2009 -> Included handling of leading and trailing NaNs
% 2.00.01 23 Nov,   2009 -> Minor M-lint stuff

dimX = size(X);
dimT = size(T);
if nargin < 3, Options = []; end

% Default options, max shift is half the target's length
OptionsDefault                          = [-fix(dimT(end) * 0.5),fix(dimT(end) * 0.5) ndims(T) 1 NaN];
Options(end + 1:length(OptionsDefault)) = OptionsDefault(length(Options) + 1:length(OptionsDefault));
Options(isnan(Options))                 = OptionsDefault(isnan(Options));
if Options(1) > Options(2), error('Lower bound for shift is larger than upper bound'), end
TimeDim = Options(3);
if not(isequal(dimX([2:TimeDim - 1,TimeDim + 1:end]),dimX([2:TimeDim - 1,TimeDim + 1:end])))
   error('Target and signals do not have compatible dimensions')
end

% Normalise to length 1 along time dimension
ord        = [TimeDim,2:TimeDim - 1,TimeDim + 1:ndims(X),1]; % Leaves the sample mode as last
X_fft      = permute(X,ord);
X_fft      = reshape(X_fft,dimX(TimeDim),prod(dimX(ord(2:end))));
X_fft      = full(X_fft / sparse(1:prod(dimX(ord(2:end))),1:prod(dimX(ord(2:end))),sqrt(SumwNaN(X_fft.^2))));
T          = permute(T,ord);
T          = reshape(T,dimT(TimeDim),prod(dimT(ord(2:end))));
T          = Normalise(T);
[nP,mP]    = size(X_fft);
nT         = size(T,1);

% Remove leading NaN's
flag_miss = any(isnan(X_fft(:))) || any(isnan(T(:)));
if flag_miss
   
   if ndims(X) > 2, error('Multidimensional handling of missing not implemented, yet'), end
   MissOff = NaN(1,mP);
   for i_signal = 1:mP
      
      Limits = RemoveNaN([1,nP],X_fft(:,i_signal)',@all);
      if ~isequal(size(Limits),[1,2]), error('Missing values can be handled only if leading or trailing'); end
      if any(cat(2,Limits(1),mP - Limits(2)) > max(abs(Options(1:2)))),
         error('Missing values band larger than largest admitted shift'),
      end
      MissOff(i_signal) = Limits(1);
      if MissOff(i_signal) > 1, X_fft(1:Limits(2) - Limits(1) + 1,i_signal)  = X_fft(Limits(1):Limits(2),i_signal); end
      if Limits(2) < nP, X_fft(Limits(2) - Limits(1) + 1:nP,i_signal) = 0; end
      
   end
   Limits                            = RemoveNaN([1,nT],T',@all);
   T(1:Limits(2) - Limits(1) + 1,:)  = T(Limits(1):Limits(2),:);
   T(Limits(2) - Limits(1) + 1:nP,:) = 0;
   MissOff                           = MissOff(1:mP) - Limits(1);
   
end

% Zero pad to avoid pollution (cf. Press & Teukolsky pg. 540 and 545)
X_fft = cat(1,X_fft,zeros(max(abs(Options(1:2))),prod(dimX(ord(2:end)))));
T     = cat(1,T,zeros(max(abs(Options(1:2))),prod(dimT(ord(2:end)))));

% Calculate cross-correlation
len_fft = max(size(X_fft,1),size(T,1));
Shift   = Options(1):Options(2);
if Options(1) < 0 && Options(2) > 0
   ind = [len_fft + Options(1) + 1:len_fft,1:Options(2) + 1];
elseif Options(1) < 0 && Options(2) < 0
   ind = len_fft + Options(1) + 1:len_fft + Options(2) + 1;
elseif Options(1) < 0 && Options(2) == 0
   ind = [len_fft + Options(1) + 1:len_fft + Options(2),1];
else % Options(1) >= 0 && Options(2) > 0
   ind = Options(1) + 1:Options(2) + 1;
end
X_fft = fft(X_fft,len_fft,1);
T_fft = fft(T,len_fft,1);
T_fft = conj(T_fft);
T_fft = reshape(T_fft(:,:,ones(dimX(1),1)),size(X_fft));
cc    = ifft(X_fft .* T_fft,len_fft,1);
cc    = reshape(cc(ind,:),[Options(2) - Options(1) + 1,prod(dimX(ord(2:end - 1))),dimX(1)]);
if Options(4) == 0
   cc = squeeze(mean(cc,2))';
elseif Options(4) == 1
   cc = squeeze(prod(cc,2))';
else
   error('Invalid options for correlation of multivariate signals')
end
[nil,pos]  = max(cc,[],2);
Values     = cat(1,Shift,cc);
Shift      = Shift(pos);
if flag_miss, Shift = Shift + MissOff; end
Xwarp      = NaN([dimX(1),dimT(2:end)]);
[ind,indw] = deal(repmat({':'},ndims(X),1));
for i_X = 1:dimX(1)
   
   ind{1}  = i_X;
   indw{1} = i_X;
   if Shift(i_X) >= 0
      
      ind{TimeDim}  = Shift(i_X) + 1:dimX(TimeDim);
      indw{TimeDim} = 1:(dimX(TimeDim) - Shift(i_X));
      if Options(5) == -inf
         ind{TimeDim}  = cat(2,ind{TimeDim},dimX(TimeDim(ones(1,abs(Shift(i_X))))));
         indw{TimeDim} = cat(2,indw{TimeDim},dimX(TimeDim) - Shift(i_X) + 1:dimX(TimeDim));
      end
      
   elseif Shift(i_X) < 0
      
      ind{TimeDim}  = 1:(dimX(TimeDim) + Shift(i_X));
      indw{TimeDim} = -Shift(i_X) + 1:dimX(TimeDim);
      if Options(5) == -inf
         ind{TimeDim}  = cat(2,ones(1,-Shift(i_X)),ind{TimeDim});
         indw{TimeDim} = cat(2,1:-Shift(i_X),indw{TimeDim});
      end
      
   end
   Xwarp(indw{:}) = X(ind{:});
   
end

function DispInfo(varargin)
% Displays the 'tag' (or the 'userdata' if a cell vector of strings).
% The tag remains visible on the plot for 5 seconds (the time can be
% changed setting the DispInfoDelay property of the figure). It changes
% linewidth and color for the same amount of time as to highlight the object.
%
% DispInfo(Handle);
%
% Inputs:
% Handle: handle of the object whose tag is to be displayed (optional). Callback object is used if Handle
%         is not specified
%
% Outputs:
% None
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Created      : 05 March, 2003; 18:37
% Last modified: 09 March, 2009; 14:18

if nargin < 1
   h = gcbo;
else
   h = varargin{1};
end
s = get(h,'tag');
g = get(h,'userdata');
ex = getappdata(h,'Extra');
if ~isempty(g) && isa(g,'cell') && all(cellfun(@(x) isa(x,'char'),g))
   s = g;
end
a = get(gcf,'units');
set(gcf,'units','points')
t = get(gcf,'currentpoint');
set(gcf,'units',a)
a = axes('units','points','position',[t,50,10],'box','on','xtick',[],'ytick',[]);
b = text(.51,.5,s,'verticalalignment','middle','horizontalalignment','center','fontsize',9);
u = get(b,'units');
set(b,'units','points')
c = get(b,'extent');
set(b,'units',u);
set(a,'position',[t,c(3) * 1.05 c(4) * 1.2])
if isequal(get(h,'type'),'line')
   
   d = get(h,'linewidth');
   e = get(h,'color');
   if isequal(e,'r') || isequal(e,[1 0 0])
      col = [0 0 0.8];
   else
      col = 'r';
   end
   set(h,'linewidth',d+1,'color',col)
   
elseif isequal(get(h,'type'),'patch')
   
   d = get(h,'linewidth');
   e = get(h,'edgecolor');
   if isequal(e,'r') || isequal(e,[1 0 0])
      col = [0 0 0.8];
   else
      col = 'r';
   end
   set(h,'edgecolor',col,'linewidth',d + 1);
   
end
if isempty(getappdata(gcf,'DispInfoDelay'))
   Delay = 5;
else
   Delay = getappdata(gcf,'DispInfoDelay');
end
t = timer('StartDelay',Delay,'TimerFcn',@RemInfoAx,'userdata',{d e h,a},'stopfcn',@Delt);
start(t);
setappdata(h,'flag_AvoidEndlessRecursion',0)
if ~isempty(ex)
   a = getappdata(ex,'flag_AvoidEndlessRecursion');
   if isempty(a) || a, arrayfun(@DispInfo,ex); end
end

function RemInfoAx(varargin)
o = get(varargin{1},'userdata');
[d e h,a] = deal(o{:});
if isequal(get(h,'type'),'line')
   set(h,'linewidth',d,'color',e)
elseif isequal(get(h,'type'),'patch')
   set(h,'edgecolor',e,'linewidth',d);
end
setappdata(h,'flag_AvoidEndlessRecursion',1)
delete(a)

function Delt(varargin)
delete(varargin{1})

function [An,flag] = RemoveNaN(B, Signal, Select)
% Rearrange segments so that they do not include NaN's
%
% [Bn] = RemoveNaN(B, Signal, Select)
%
% INPUT
% B     : (p × 2) Boundary matrix (i.e. [Seg_start(1) Seg_end(1); Seg_start(2) Seg_end(2);...]
% Signal: (n × 2) Matrix of signals (with signals on rows)
% Select: (1 × 1) function handle to selecting operator
%                 e.g. @any (default) eliminate a column from signal matrix
%                                     if one or more elements are missing
%                      @all           eliminate a column from signal matrix
%                                     if all elements are missing
%
% OUTPUT
% Bn  : (q × 2)     new Boundary matrix in which NaN's are removed
% flag: (q × 2 × n) flag matrix if there are NaN before (column 1) or after (column 2)
%                   the corresponding segment in the n signals.
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Created      : 25 February, 2009
% Last modified: 23 March, 2009; 18:02

% HISTORY
% 1.00.00 09 Mar 09 -> First working version
% 2.00.00 23 Mar 09 -> Added output for adjacent NaN's in signals
% 2.01.00 23 Mar 09 -> Added Select input parameter

error(nargchk(2,3,nargin))
if nargin < 3, Select = @any; end
C      = NaN(size(B));
count  = 0;
Signal = isnan(Signal);
for i_el = 1:size(B,1)
   
   ind = B(i_el,1):B(i_el,2);
   in  = Select(Signal(:,ind),1);
   if any(in)
      
      p = diff([0 in],1,2);
      a = find(p < 0);
      b = find(p > 0) - 1;
      if ~in(1), a = cat(2,1,a); else b = b(2:end); end
      if ~in(end), b = cat(2,b,length(ind)); end
      a = unique(a);
      b = unique(b);
      C(count + 1:count + length(a),:) = ind(cat(2,a(:),b(:)));
      count = count + length(a);
      
   else
      count      = count + 1;
      C(count,:) = B(i_el,:);
   end
   
end
An = C;
if nargout > 1
   
   flag            = false(size(C,1),2,size(Signal,1));
   Cinds           = C(:,1) > 1;
   Cinde           = C(:,2) < size(Signal,2);
   flag(Cinds,1,:) = Signal(:,C(Cinds,1) - 1)';
   flag(Cinde,2,:) = Signal(:,C(Cinde,2) + 1)';
   
end

function Cols = ExtCM(n,CM)
% Extend colormap to have n colors
if nargin < 2
   CM = MyCM;
end
if isa(CM,'char')
   CM = eval(CM);
end
if isa(CM,'function_handle')
   CM = CM();
end
if ~isequal(size(CM),[64,3]) || any(CM(:) < 0 | CM(:) > 1)
   error('CM argument is not a colormap')
end
if n < 0 || ~isfinite(n)
   error('n must be positive and finite')
end
Cols           = interp1((1:64)',CM,linspace(1,64,n)');
Cols(Cols < 0) = 0;
Cols(Cols > 1) = 1;

function Colours = MyCM
% Purple to orange colour map
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
Colours(:,1) = linspace(0,1,64)';
Colours(:,2) = abs(linspace(-.5,.5,64)');
Colours(:,3) = flipud(Colours(:,1));

function [S] = SumwNaN(X, Dim, flag)
% Sum with missing values
%
% [S] = SumwNaN(X, Dim, flag)
% 
% INPUT
% X   : data array
% Dim : dimension along which the sum is to be performed
% flag: true if any missing / missingness of X
% 
% OUTPUT
% S: sum along Dim of X ignoring missing values
% 
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
% 
% Created      : 10 March, 2009; 15:59
% Last modified: 24 November, 2009; 16:03
if nargin < 1,error('not enough input arguments'); end
if nargin < 2,Dim = 1;end
if nargin < 3,flag = [];end
if isempty(flag) || (numel(flag) == 1 && flag)
   flag = isnan(X);
   N    = true;
elseif isequal(size(flag),size(X)), N = true;
elseif numel(flag) == 1 && ~flag,   N = false;
else error('Invalid flag'),
end
if N, X(flag) = 0; end
S = sum(X,Dim);

function [S] = MeanwNaN(X, Dim, flag)
% Sum with missing values
%
% [S] = MeanwNaN(X, Dim, flag)
% 
% INPUT
% X   : data array
% Dim : dimension along which the mean is calculated
% flag: true if any missing / missingness of X
% 
% OUTPUT
% S: sum along Dim of X ignoring missing values
% 
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
% 
% Created      : 24 November, 2009; 09:04
% Last modified: 24 November, 2009; 09:25

% HISTORY
% 1.00.00 2009 Nov 24 -> Dunction derived from SumwNaN
error(nargchk(1,3,nargin))
if nargin < 2,Dim = 1;end
if nargin < 3,flag = [];end
if isempty(flag) || (numel(flag) == 1 && flag)
   flag = isnan(X);
   N    = true;
elseif isequal(size(flag),size(X)), N = true;
elseif numel(flag) == 1 && ~flag,   N = false;
else error('Invalid flag'),
end
if N, 
   X(flag) = 0; 
   flag    = ~flag;
   c       = sum(flag,Dim);
   c(~c)   = NaN;
else
   c = size(X,Dim);
end
S = sum(X,Dim) ./ c;

function [Xn] = Normalise(X, Flag)
% Column-wise normalise matrix
% NaN's are ignored
%
% [Xn] = Normalise(X, Flag)
%
% INPUT
% X   : Marix
% Flag: true if any NaNs are present (optional - it saves time for large matrices)
%
% OUTPUT
% Xn: Column-wise normalised matrix
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Created      : 09 March, 2009; 13:18
% Last modified: 09 March, 2009; 13:50
if nargin < 2
   Patt = ~isnan(X);
   Flag = any(~Patt(:));
else
   if Flag, Patt = ~isnan(X); end
end
[M,N] = size(X);
Xn    = NaN(M,N);
if Flag
   
   for i_n = 1:N
      n = norm(X(Patt(:,i_n),i_n));
      if ~n, n = 1; end
      Xn(Patt(:,i_n),i_n) = X(Patt(:,i_n),i_n) / n;
   end
   
else
   
   for i_n = 1:N
      n = norm(X(:,i_n));
      if ~n, n = 1; end
      Xn(:,i_n) = X(:,i_n) / n;
   end
   
end

function [XSeg,SegNew] = ExtractSegments(X, Segments)
% Extract segments from signals
%
% [XSeg] = ExtractSegments(X, Segments)
%
% INPUT
% X       : (n × p) data matrix
% Segments: (s × 2) segment boundary matrix
%
% OUTPUT
% XSeg: (n × q) data matrix in which segments have been removed
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Created      : 23 March, 2009; 07:51
% Last modified: 23 March, 2009; 15:07

% HISTORY
% 0.00.01 23 Mar 09 -> Generated function with blank help
% 1.00.00 23 Mar 09 -> First working version

% Some initialisations
[n,p] = size(X);
Sd    = diff(Segments,1,2) + 1;
q     = sum(Sd + 1);
[s,t] = size(Segments);

% Check input
error(nargchk(2,2,nargin,'struct'))

flag_si = t ~= 2;
flag_in = any(Segments(:) ~= fix(Segments(:)));
flag_ob = any(Segments(:,1) < 1) || any(Segments(:,2) > p);
flag_ni = any(diff(Segments(:,1)) < 0) || any(diff(Segments(:,2)) < 0);
flag_ab = any(Sd < 2);
if flag_si, error('Segment boundary matrix must have two columns'), end
if flag_in, error('Segment boundaries must be integers'), end
if flag_ob, error('Segment boundaries outside of segment'), end
if flag_ni, error('Segments boundaries must be monotonically increasing'), end
if flag_ab, error('Segments must be at least two points long'), end

% Initialise output
if nargout > 0, XSeg = NaN(n,q); else return, end

% Calculate Segments
cdim = cat(1,0,cumsum(Sd));
for i_seg = 1:s
   XSeg(:,cdim(i_seg) + 1:cdim(i_seg + 1)) = X(:,Segments(i_seg,1):Segments(i_seg,2));
end
if nargout > 1, SegNew = cat(2,cdim((1:s)') + 1,cdim((2:s + 1)')); end

function pts = scal2pts(ppmi,ppm,prec)
% Transforms scalars in data points
%
% pts = scal2pts(values,scal)
%
% INPUT
% values: scalars whose position is sought
% scal  : vector scalars
% prec  : precision (optional) to handle endpoints
%
% OUTPUT
% pts   : position of the requested scalars (NaN if it is outside of 'scal')
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Created      : 12 February, 2009; 17:43
% Last modified: 11 March, 2009; 15:14

% HISTORY
% 1.00.00 12 Feb 09 -> First working version
% 1.01.00 11 Mar 09 -> Added input parameter check
if nargin < 3, prec = min(abs(unique(diff(ppm)))); end
dimppmi = size(ppmi);
ppmi    = ppmi(:);
ppm     = ppm(:);
rev     = ppm(1) > ppm(2);
if rev, ppm = ppm(end:-1:1); end
ubound  = (ppmi - ppm(end)) < prec & (ppmi - ppm(end)) > 0;
lbound  = (ppm(1) - ppmi) < prec & (ppm(1) - ppmi) > 0;
ppmi(ubound) = ppm(end);
ppmi(lbound) = ppm(1);
if nargin < 2, error('Not enough input arguments'); end
if length(ppmi) > length(ppm), warning('icoshift:scal2pts','ppm vector is shorter than the value''s'), end
[xxi,k]              = sort(ppmi(:));
[nil,j]                = sort([ppm(:);xxi(:)]);
r(j)                 = 1:length(j);
r                    = r(length(ppm) + 1:end) - (1:length(ppmi));
r(k)                 = r;
r(ppmi == ppm(end))  = length(ppm);
ind                  = find((r > 0) & (r <= length(ppm)));
ind                  = ind(:);
pts                  = Inf(size(ppmi));
pts(ind)             = r(ind);
ptsp1                = min(length(ppm),abs(pts + 1));
ptsm1                = max(1,abs(pts - 1));
ind                  = find(isfinite(pts));
dp0                  = abs(ppm(pts(ind)) - ppmi(ind));
dpp1                 = abs(ppm(ptsp1(ind)) - ppmi(ind));
dpm1                 = abs(ppm(ptsm1(ind)) - ppmi(ind));
pts(ind(dpp1 < dp0)) = pts(ind(dpp1 < dp0)) + 1;
pts(ind(dpm1 < dp0)) = pts(ind(dpm1 < dp0)) - 1;
if isempty(pts), pts = []; end
pts(not(isfinite(pts))) = NaN;
if rev, pts = length(ppm) - pts + 1; end
if ~isequal(size(pts),ppmi), pts = reshape(pts,dimppmi); end

function I = dscal2dpts(d,ppm,varargin)
% Translates an interval width from scal to the best approximation in sampling points.
%
% I = dppm2dpts(Delta,scal,prec)
%
% INPUT
% Delta: interval widths in scale units
% scal : scale
% prec : precision on the scal axes
%
% OUTPUT
% I: interval widths in sampling points
%
% Author: Giorgio Tomasi
%         Giorgio.Tomasi@gmail.com
%
% Last modified: 21st February, 2009
%
if isempty(d), I = []; return; end
if nargin < 2, error('Not enough input arguments'); end
if d <= 0; error('Delta must be positive'); end
if ppm(1) < ppm(2), I = scal2pts(ppm(1) + d,ppm,varargin{:}) - 1;
else I = length(ppm) - scal2pts(ppm(end) + d,ppm,varargin{:}) + 1;
end