
function [ optima, yoptima, plateaux] = calcOptima( x, y , methode, ydeb, seuil,  nf, couleur, hp2)
%calcOptima Calcule les optima de la réponse y d'une fonction=f(x)
% On retourne les x où la derivée de y par rapport à x change de signe:
% puisqu'on recherche les optima, il s'agit des x où elle passe de positif
% à négatif. On détecte aussi le dernier point des plateaux
% INPUT
%    x, y: abscisses et réponses (vecteur de même longueur)
%    methode: la méthode utilisée (string: 'Tornsey' ou 'Sobol' ou 'Sivip', voir NOTE)
%    ydeb: valeur de l'ordonnee du point precedant x
%    seuil: valeur a partir de laquelle un pic est oté (utile pour les
%    méthodes Sobol et sivip)
%    nf: nombre max d'optima voulu(utile pour la
%    méthode Tornsey)
%    couleur: triplet RGB de la couleur du tracé.
%           facultatif: Si cet argument est présent, on trace les y en
%           fonction des x en localisant les optima.
%           Si cet argument est manquant, pas de tracé graphique
%   hp2: handle de la figure parent du graphe
% OUTPUT
%    optima: les abscisses où la réponse atteint un optimum
%    yoptima: la valeur de y correspondante, si méthode est 'Sivip', ignoré
%    sinon
%    plateaux:  eventuels plateaux pour Tornsey, ignoré sinon:
%      Structure dont chaque element est un quadruplet:
%      abscisse de début du plateau, abscisse de fin du plateau, hauteur
%      moyenne, longueur
% NOTE
%    - si la méthode est 'Tornsey', on ne garde que les nf plus grands optima
%    - si la méthode est 'Sobol', on ne garde que les optima >= seuil
%    - si la méthode est 'Sivip', on ne garde que le 1ier optimum >= seuil
% =========================================================================
optima=[];
yoptima={}; % pour sivip
plateaux={}; %structure contenant les abscisses des eventuels plateaux pour Tornsey

%% On rajoute un 1ier point et un dernier point a zero pour amorcer la recherche des optima
% d'abord, on dimensionne les entrées sous forme d'un vecteur ligne
x=reshape(x, 1, max(size(x)));
y=reshape(y, 1, max(size(y)));
x=[-1 x]; %1ier point d'abscisse -1
y=[ydeb y]; %1ier point d'ordonnée ydeb
x(end+1)= max(x)+1; %dernier point
y(end+1)=-1; %dernier point d'ordonnée -1

% derivée de y par rapport à x
deriv_y = diff(y)./diff(x);

xd = x(2:length(x)); % Backward difference x values

delay_mult = deriv_y(1:length(deriv_y)-1) .*deriv_y(2:length(deriv_y));
% icritical: indices dans x où la derivée change de signe
icritical = find(delay_mult < 0);
% critical: abscisse où la derivée change de signe
d1= deriv_y(icritical);
d2= deriv_y(icritical+1);
% ioptpic: indices dans x où la derivée passe de positif à négatif
ioptpic= icritical(find(d1>d2));

% 3/6/2013: pour detecter les plateaux, on arrondit les petites differences
% des derivees pour qu'elles soient assimilées a des zeros
eps=nf/1.0e4;
deriv_y(abs(deriv_y) < eps)=0.0;

% detecter les plateaux, i.e là ou la derivee est nulle
icritical = find(deriv_y(1:length(deriv_y)-1)==0);
d1= deriv_y(icritical);
d2= deriv_y(icritical+1);
% ioptplat: indices dans x où la derivée passe de nul à négatif (fin de
% plateau, suivie d'une decroissance)
ioptplat= icritical(find(d1>d2));
% Mettre les 2 resultats en vecteur
if ~isempty(ioptpic)
    ioptpicv = reshape(ioptpic, 1, max(size(ioptpic)));
else
    ioptpicv=[];
end
if ~isempty(ioptplat)
    ioptplatv= reshape(ioptplat, 1, max(size(ioptplat)));
else
    ioptplatv=[];
end
iopt=[ioptpicv ioptplatv]; % les indices des optima dans xd (décalés de 1)
iopt=iopt+1; % les indices des optima dans x et y

switch methode
    case 'Tornsey'
        %% Ne garder que les nf plus grands (pour Tornsey seumement)
        if length(iopt) > nf
            maxy=sort(y(iopt), 'descend');
            maxy=maxy(1:(nf));
            [Lia,iopt] = ismember(maxy, y); % iopt: indices correspondants dans x et y
        end
        
        optima=x(iopt); % les abscisses correspondant a un des nf plus grands optima
        
        %Rajouter le point 0
        optima(end+1)=0;
        % eliminer les doubles (le 0 et le dernier point peuvent etre en
        % double) et trier
        optima=sort(unique(optima));
        % Détecter les eventuels plateaux s'ils sont de
        % longueur suffisante
        % seuil: temps maximal en heure de detection des plateau (minimal time step)
        if seuil >0 && ~isempty(icritical)
            % detecter les points consécutifs des plateaux
            ideb=1;
            while ideb < length(icritical)
                i=ideb+1;
                while i <= length(icritical)
                    if icritical(i) ~=icritical(i-1)+1
                        break
                    end
                    i=i+1;
                end % fin while i
                ifin=i-1;
                xdeb=x(icritical(ideb)+1);
                xfin=x(icritical(ifin)+1);
                lgplat= xfin - xdeb;
                if lgplat > seuil
                    plateaux{end+1}= [xdeb xfin mean(y((icritical(ideb)+1):(icritical(ifin)+1))) lgplat];
                end
                ideb=i;
            end % fin while ideb
        end %fin seuil
        
    case 'Sobol'
        %% Eliminer les optima correspondant a un IST <seuil
        optima=x(iopt);
        optima(find(y(iopt)<seuil)) =NaN;
        optima=optima(~isnan(optima));
        optima=sort(optima);
        %% On rajoute le dernier temps (pour les graphiques, il faut au moins un
        % optima, sinon la légende va etre décalée)
        % indice=end-1 car on a rajouté un point bidon
        if length(optima)==0 || optima(end) ~= x(end-1)
            optima(end+1) = x(end-1);
        end
    case 'Sivip'
        %% SIVIP
        sol=2; % on fixe le mode d'élimination des petits pics
        if sol==1
            % version ou on ote les petits pics relatifs
            while (~isempty(iopt) &&  y(iopt(1)) < seuil)
                % eliminer les optima < seuil
                disp(['The time ' num2str(x(iopt(1))) ' is ignored because its optima, ' num2str(y(iopt(1))) ...
                    ', is less than the threshold']);
                iopt(1)=[];
                
            end % fin while
            % eliminer les optima qui dépassent de moins du seuil du point
            % précédent et du point suivant
            ntemps=length(x)-1;
            while  (~isempty(iopt)) && ...
                    (iopt(1) >1) && ( iopt(1)< ntemps) &&...
                    (y(iopt(1)) - y(iopt(1)-1) < seuil) && ...
                    (y(iopt(1)) - y(iopt(1)+1) < seuil)
                disp(['Le temps ' num2str(x(iopt(1))) ' est ote car son optima, ' num2str(y(iopt(1))) ...
                    ' dépasse de moins du seuil du point précédent ', num2str(y(iopt(1)-1)), ...
                    ' et du point suivant ', num2str(y(iopt(1)+1))])
                iopt(1)=[];
                
            end % fin while
            if ~isempty(iopt)
                optima=x(iopt(1));
                yoptima=y(iopt(1));
            else
                % si pas d'optimum sur la suite des temps: le dernier point
                optima=x(end-1);
                yoptima=y(end-1);
            end
        end %%fin sol==1
        
        if sol==2
            
            if ~isempty(iopt)
                optima=x(iopt);
                yoptima=y(iopt);
                while length(optima) >0 && yoptima(1) < seuil
                    disp(['The time ' num2str(optima(1)) ' is removed because its optima, ' num2str(yoptima(1)) ...
                        ' is less than the threshold']);
                    yoptima(1)=[];
                    optima(1)=[];
                end
                if length(optima) >0
                    optima=optima(1);
                    yoptima=yoptima(1);
                else
                    % pas d'optimum>=seuil sur la suite des temps: le dernier point
                    optima=x(end-1);
                    yoptima=y(end-1);
                end
            else
                % si pas d'optimum sur la suite des temps: le dernier point
                optima=x(end-1);
                yoptima=y(end-1);
            end
        end % fin sol=2
        
        
        
    otherwise
        errordlg(['Internal error calcOptima: method not recognized reconnue ' methode]);
        return
end % fin switch



%% Trace graphique
if (nargin > 6)
    % Restaurer x et y: oter le 1ier et dernier éléments bidons
    x(1,1)=NaN; x(1,end)=NaN;
    x=x(~isnan(x));
    y(1,1)=NaN; y(1,end)=NaN;
    y=y(~isnan(y));
    plot(hp2, x,y,  couleur{1})
    xlim([min(x) max(x)]);
    hold on
    [Lia,Locb] = ismember(optima,x);
    plot(hp2, optima ,y(Locb), ['s' couleur{1}]);
    hold on
    
end


end % fin calcOptima
%% ========================================================================
