function res = probadis(x,p)
%res = PROBADIS(x,p)
%        Input   x     vector of real numbers (support points in IR)
%                p     vector of probability weights associated to x
%                      i.e. non-negative real numbers such that sum(p) == 1
%        Output  res   random number chosen from the finite discrete 
%                      distribution on x(1),...,x(n) with probability weights
%                      p(1),...,p(n) where n == length(x) == length(p)
%
% Renvoie UNE réalisation de la loi discrète à support fini sur IR dont les 
% points de support sont les composantes du vecteur x, et les poids sont les 
% composantes du vecteur p. Donc x(i) a une probabilité p(i) d'être renvoyé.
%
% La méthode consiste à écrire l'intervalle [0,1] comme une réunion disjointe
% d'intervalles de longueurs p(1),..., p(n) puis à regarder au quel appartient
% la réalisation d'une loi uniforme obtenue par un appel à rand. Les valeurs 
% renvoyées par des appels successifs sont donc pseudo-indépendantes. 
% Pour des raisons d'efficacité, les conditions nécessaires suivantes ne sont 
% pas contrôlées par cette fonction :
%  - le nombre de paramètres passés est exactement 2
%  - les deux paramètres x et p sont bien des vecteurs et sont de même longueur
%  - les composantes de p sont positives ou nulles et leur somme vaut 1.
%
% See also RDISCR.
%
% ### Copyright (C) D. Chafaï, 2003-12-06.
% ### http://www.lsp.ups-tlse.fr/Chafai/agregation.html
% ### Licence GNU General Public License http://www.gnu.org/copyleft/gpl.html
%

% On pourrait implémenter cette fonction de la façon suivante :
%
%  INDICES = find(cumsum(p) >= rand);
%  res = x(INDICES(1));
%  return;
%
% Cette méthode est correcte mais inefficace car elle ne tient pas compte de 
% la monotonie de cumsum(p), ce qui entraine des tests inutiles une fois que
% la valeur critique a été franchie. En jargon, c'est un 'firstmatch' qu'il 
% nous faut, pas un 'matchall'. Le nombre de if impliqués dans find est 
% toujours égal à la taille de ce qu'on lui passe en paramètre.
%
% Morale de l'histoire : la brièveté d'un code n'assure pas sa performance !
% Et l'abscence de if dans un code ne signifie pas qu'il ne fait pas appel
% indirectement à des if, et encore moins que cela est fait de manière 
% optimale ! La fonction find n'est qu'une boucle for contenant un if.
%
% Ci-dessous, nous utilisons une version plus rapide dans le plus pur style
% for-if avec un nombre de if optimal. On pourrait adapter l'ordre des 
% intervalles testés (et donc l'arbre associéé) aux poids p(i) de façon à
% tester d'abord les intervalles les plus probables. Est-ce vraiment mieux ?
% Exercice !

n = length(x);    % le nombre d'atomes.
r = rand;         % une réalisation de loi uniforme sur [0,1].
a = 0; b = p(1);  % [a,b] = sous-intervalle de proba p(i) pour l'uniforme.
for i = 1:n-1     % parcours de tous les sous-intervalles juxtaposés.
   if ((r >= a) & (r < b))
      res = x(i);
      return;     % on a trouvé le bon intervalle, on sort.
   end
   a = b; b = b + p(i+1); % on passe à l'intervalle suivant.
end
res = x(n); % le bon intervalle est le dernier.
return;