function  [IST, errIST, optima, optimaP,  titregraphe, soustitre, indice_panel] = calcStrategieOptSobolDP(N, tmax, pas, nf, graine, binf, bsup, ...
    modelepop,  nomparam,nomparamP, nomparamS, uniteT_libelle, ...
    seuil, meth, hpanel_work, obj, evnt)
%calcStrategieOpt Simulation d'une stratégie optimale par la methode
%SOBOL-SALTELLI
% INPUT
%    N: nombre de modalités (nombre de jeux de paramètres dans chaque
%    matrice)
%    tmax: temps maximal
%    nf: nombre de paramètres
%    graine aléatoire
%    binf, bsup: plages de variation des paramètres (nf)
%    modelepop: le programme de calcul du modèle
%    nomparam: le nom de tous les paramètres (nf)
%    nomparamP: le nom des parametres actifs
%    nomparamS: le nom des parametres fixes
%    seuil: valeur en % a partir de laquelle, un optima est rejeté
%    meth: 1 si on genere le lhd en triant des pseudorandom integer values
%          2 si on genere le lhd en utilisant sobolset
%   hpanel_work: handle de la figure parent du graphe
%   obj: handle du uicontrol appelant
%    evnt: 1 quand on fait passer des tests
% OUTPUT
%    IST: indices de sensibilite totaux (tmax+1, nf)
%    errIST: message d'erreur ou ''
%    optima: abscisses où les IST sont optimum (vecteur)
%    optimaP: les optima pour chaque facteur (:, nf)
%    titregraphe: titre du graphe
%    soustitre: libellé contenant les parametres de l'étude
%    indice_panel: indice du dernier panel écrit
%% ========================================================================
%% Initialisation des sorties, au cas ou l'utilisateur arrete avant leur calcul
IST=[];
errIST='';
optima=[]; % longueur=nombre d'optima, inconnu a priori
optimaP=[]; % les optima pour chaque parametre (sert pour sauvegarde/restauration)
titregraphe='';
soustitre='';
indice_panel=0;

%% Creation d'un bouton 'Wait'
hbar = waitbar(0,'1','Name','Please wait...', ...
    'CreateCancelBtn', ...
    'setappdata(gcbf,''canceling'',1)');% User can click Cancel or close button
% nbre d'etapes de calcul
nsteps=nf+5;
istep=0; % indice de l'étape courante
setappdata(hbar,'canceling',0) % état= pas cancelé


%% Initialisation de la graine
try
    rng(str2num(graine));
catch
    rng('shuffle');  % seeds the random number generator based on the current time
end


%% Generer les matrices de parametres A et B

% Rque:  taille (N x nf)
istep=istep+1;
waitbar(istep/nsteps,hbar, 'Generation of the parameter matrices A and B');
[matA, matB]= generABDP( N, nf,  binf, bsup, meth);
% matA, matB: matrices double precision (N,nf)
% Check for Cancel button press
if isempty(evnt)
    if getappdata(hbar,'canceling') || get(obj,'userdata')==0
        delete(hbar);
        hm=msgbox(['Step ', num2str(istep), ' interrupted']);
        waitfor(hm)
        return
    end
end



%% Generer les matrices C
%  tableau N x nf x nf.
istep=istep+1;
waitbar(istep/nsteps,hbar, 'Création des matrices C')
% Check for Cancel button press
if isempty(evnt)
    if getappdata(hbar,'canceling')|| get(obj,'userdata')==0
        delete(hbar);
        hm=msgbox(['Step ', num2str(istep), ' interrupted']);
        waitfor(hm)
        return
    end
end


%% Calculer le modele  sur tous ces jeux de paramètres
% Oter les eventuels  blancs dans le nom du modele
lemodele = strrep(modelepop,' ', '');
% Oter les eventuels  tirets dans le nom du modele
lemodele = strrep(lemodele,'-', '');
% Les bornes admissibles de la réponse sont determinées par les
% plages de variation de certains parametres, selon le modele
[valmin, valmax] = init_modeles_repvalides(lemodele, nomparam, binf, bsup);

%% Les temps de calcul
%% On fait des pas de temps pour que les matrices ne soient pas
% trop grosses, quand nmoda*ntemps est grand
temps=0:pas:tmax;
ntemps=length(temps);


%% Calcul pour la matrice A
istep=istep+1; % étape de calcul pour la barre de progression
waitbar(istep/nsteps,hbar, ['Compute the model ', modelepop, ' for the matrix A' ])
YA=zeros(N, ntemps);

% Appel de la version C du modele
 eval(['C', lemodele, '(double(matA), double(temps), double(valmin), double(valmax), double(YA) )']);
 
% Check for Cancel button press
if isempty(evnt)
    if getappdata(hbar,'canceling') ||  get(obj,'userdata')==0
        delete(hbar);
        hm=msgbox(['Step ', num2str(istep), ' stopped']);
        waitfor(hm)
        return
    end
end

%% Ramener a des valeurs valides
	  YA(YA<0) =0.0;


%% Calcul du denominateur
f0= sum(YA,1);
f0= (f0 .* f0)/N;
den= calcSobolIST(YA, YA) - f0;
clear YA;


%% Calcul pour la matrice B
istep=istep+1; % étape de calcul pour la barre de progression
waitbar(istep/nsteps,hbar, ['Compute the model ', modelepop, ' for the matrix B' ])
YB=zeros(N, ntemps);

% Appel de la version C du modele
eval(['C', lemodele,  '(double(matB), double(temps), double(valmin), double(valmax),double(YB) )']);


% Check for Cancel button press
if isempty(evnt)
    if getappdata(hbar,'canceling') ||  get(obj,'userdata')==0
        delete(hbar);
        hm=msgbox(['Step ', num2str(istep), ' stopped']);
        waitfor(hm)
        return
    end
end



%% Ramener a des valeurs valides
	  YB(YB<0) =0.0;

IST=zeros(ntemps, nf);
YC=zeros(N, ntemps);
arrayC=zeros(N,nf);

%% Calcul pour les nf matrices C et des IST
for j=1:nf
    istep=istep+1;
    waitbar(istep/nsteps,hbar,['Compute the model for the ', num2str(j) , '-st matrix C']);
    arrayC(:,:) = matB;
    arrayC(:,j) = matA(:,j);
    eval(['C', lemodele,  '(double(arrayC), double(temps),double(valmin), double(valmax), double(YC) )']);
   
    % Check for Cancel button press
    if isempty(evnt)
        if getappdata(hbar,'canceling') ||  get(obj,'userdata')==0
            delete(hbar);
            hm=msgbox(['Step ', num2str(istep), ' interrupted']);
            waitfor(hm)
            return
        end
    end



%% Ramener a des valeurs valides
	  YC(YC<0) =0.0;


	sumBC= calcSobolIST(YB, YC);
    IST(:, j) = 1 - ((sumBC -f0) ./ (den));
    
end % fin nf
clear matA matB arrayC;
clear  YB YC
delete(hbar); % a faire avant les graphiques

% VERIF: a un temps donné, il faut qu'il y ait au moins 1 IST >0
indices=find(max(IST')<=0);
if ~isempty(indices)
    %ne pas afficher ce message par errordlg, a cause des graphiques qui suivent: le plot sera vide
    
    errIST= sprintf('%s\n%s\n%s\n', ...
        ' ERROR: At any time, there should be at least one positive IST. ', ...
        'There are no positive IST at the times: ', num2str(indices));
    hm=msgbox(errIST);
    waitfor(hm);
    return
end

%% Enlever les IST qui correspondent a des parametres fixes
nfS=length(nomparamS);

for i=1:nfS
    indice=find(strcmpi(nomparam, nomparamS{i}));
    IST(:,indice)=[];
    nomparam(indice)=[];
end
nf=nf-nfS;

%% Les IST <0 correspondent à des valeurs des paramètres incompatibles
IST(IST<0)=0;

%% Les transformer en % de leur somme, en chaque temps
a= sum(IST,2);
% On a deja vérifié que a>0 en chaque temps
a=repmat(a,1, nf);
IST= (100 * IST) ./ a;

%% Les IST=NaN correspondent au cas ou aucun IST>0 a l'un des temps
IST(isnan(IST))=0;


%% Pour le trace graphique
couleur={'k', 'r','g','b','k', 'm', 'c', 'k'};
if nf>length(couleur)
    errordlg('calcStrategie Opt::Internal error: pas assez de couleurs');
    return
end

%% Calcul des optima et trace graphique
%disp('*** Calcul des optima et trace graphique')
if ~ishandle(hpanel_work) %interruption via quitter
    return
end

indice_panel=1;
tag=['panel' num2str(indice_panel) '_tag'];
hp2=uipanel('parent',hpanel_work,...
    'BorderType', 'none', ...
    'tag', tag, ...
    'BackgroundColor',[1 1 1],...
    'Position',[0 0.12 1 0.88]);


h=subplot(1,1,1, 'parent',hp2);


ydeb=-1; % pour amorcer la recherche d'optimum
for p=1:nf
    e=calcOptima(temps, IST(:,p), 'Sobol', ydeb, seuil, 0, couleur(p), h);
    optima(length(optima)+1:length(e)+length(optima))= e;
    optimaP(1:length(e),p)=e;
    hold on
end
grid on
%% Trier les abscisses des optima en éliminant les doubles
optima=unique(sort(optima));
% la partie du titre qui sera reprise dans le tableau final des resultats
soustitre= ['Model: ', modelepop, ...
    ';  Sobol-Saltelli method'];
%	       '; LHS:', num2str(N, '%3.1e'), ...
%	       '; threshold:', num2str(seuil), '; seed: ', num2str(graine, '%10.0f'), ...
%    errIST];


titregraphe = [ '  {\bf{Total Sensitivity Index and maxima}} - ' , soustitre];
ylab='TSI (%)';
xlab=['times (' uniteT_libelle ')'];
annotGraphe(optima, nomparamP, titregraphe,xlab, ylab, h, hpanel_work);
% Mettre dans une structure 'filtrexplot' situé dans 'base'
% ce qui est nécessaire aux graphiques de façon à pouvoir
% les retracer si clic gauche dessus
sauve_hp('filtrexplot', 1,  'temps',  temps, ...
    'optimaP', optimaP, 'optima', optima, ...
    'nomparam', nomparamP, ...
    'titregraphe', titregraphe, 'IST', IST, ...
    'couleur', couleur);

set(gca,'buttondownfcn',@dessin_simul_sobol);
end %fin calcStrategieOptSobol
%% ========================================================================
function [unIST ] = calcSobolIST(Ya, Yb)
%calcSobolIST fonction intermediaire 
% INPUT
%    matrices (nmoda x ntemps)

% OUTPUT
%    unIST: indices de sensibilité totaux (scalaire)
%% ========================================================================
res = Ya .* Yb;
unIST =sum(res,1);

end % fin calcSobolIST




