% APRENDIZAJE AUTOMÁTICO: UN ENFOQUE PRÁCTICO
% Autores: Gonzalo Pajares, María Guijarro, P. Javier Herrera y Jesús M. de la Cruz
% Copyright RA-MA, 2010
% Capítulo 2: APRENDIZAJE ESTADÍSTICO

% parámetros del clasificador
  clear all;
  close all;
  
 %definición de parámetros para las ventanas modales de diálogo
 options.Resize='on';
 options.WindowStyle='modal';
 options.Interpreter='tex';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% Paso 1: carga de datos desde la Base de Conocimiento %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%cargar los datos existentes en la Base de Conocimiento
if exist('MuestrasAprendizaje.mat') ~= 0 
  load MuestrasAprendizaje
else
  X = [0 0 0];  
end

[d1,d2] = size(X);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Paso 2: si se selecciona si (s, por defecto) %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%         se pregunta que si la carga se hará desde fichero %%%%%%%
%%%%%%%%%         o por selección de datos con el ratón    %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 's';
prompt = {'s/n'};
dlg_title = '¿nuevos datos?';
num_lines = 1;
def = {'s'};
answer = inputdlg(prompt,dlg_title,num_lines,def,options);
n = answer{1};
if (strcmp(n,'s') ~= 0)

    n = 's';
    num_img = 0;
    while(strcmp(n,'s') ~= 0) % pregunta de más imágenes

      % numero de imágenes
      num_img = num_img + 1;
    
      if (num_img < 10) 
        IO = imread(strcat('IEntrenamiento_0',num2str(num_img),'.bmp')); 
      else
        IO = imread(strcat('IEntrenamiento_',num2str(num_img),'.bmp')); 
      end

      %leer los datos de la imagen original
      step = 1;
      [M,N,s] = size(IO);
      I = IO(1:step:M,1:step:N,:);
      figure; imshow(I);
  
      prompt = {'s/n'};
      dlg_title = '¿seleccionar muestras con el ratón?';
      num_lines = 1;
      def = {'s'};

      answer = inputdlg(prompt,dlg_title,num_lines,def,options);
      n = answer{1};
      if (strcmp(n,'s') ~= 0)
        [filas, columnas, XO] = impixel;
        [a,b] = size(filas);
        hold on
        for i=1:1:a
          plot(filas(i,:),columnas(i,:),'.b');
        end
      else
        % en esta opción se seleccionan las muestras pro filas y columnas
        % de forma automática eleigiendo filas y columnas según el paso de
        % muestreo elegido
        
        barraprogreso % Crea la figura y determina el progreso
        
        prompt = {'Paso de muestreo'};
        dlg_title = 'Paso de muestreo';
        num_lines = 1;
        def = {'4'};
        answer = inputdlg(prompt,dlg_title,num_lines,def,options);
        paso = str2num(answer{1});
   
        [m,n,s] = size(I);
        h = 0;
        for i=1:paso:m
            for j=1:paso:n
               h = h + 1;
               XO(h,:) = double(I(i,j,:)); 
            end
            barraprogreso(i/m) % actualiza figura
        end
      end
  
      X = [X;XO]; %concatenación
      Xo = X;     % datos originales
      clear XO;
      prompt = {'s/n'};
      dlg_title = '¿Más imágenes de entrenamiento?';
      num_lines = 1;
      def = {'s'};
      answer = inputdlg(prompt,dlg_title,num_lines,def,options);
      n = answer{1};
      
      if num_img == 5
        % se ha alcanzado el número máximo de imágenes disponibles
        disp('Se ha alcanzado el número máximo de imágenes de entrenamiento disponibles');
        n = 'n';
      end

  end %while(strcmp(n,'s') ~= 0) % pregunta de más imágenes

  % quitamos la primera fila de X, que se añadió al crearla
  if d1 == 1
    [m,n] = size(X);
    X = X(2:1:m,:);
  end
  save MuestrasAprendizaje X
end %nuevos datos

% en X1 se mantienen los datos originales sin normalizar
X1 = X;
% %normalizacion de datos en el rango [0,1]
MI=min(X); %minimo
MA=max(X); %maximo
X=(X-repmat(min(X),size(X,1),1))./(repmat(max(X),...
     size(X,1),1)-repmat(min(X),size(X,1),1));
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Paso 3: proceso de partición de los datos   %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% una vez obtenidos los datos, se realiza una partición de los mismos existen dos opciones:
% 1) opción ---> s: mediante el algoritmo pseudoaleatorio de Balasko    
% 2) opción ---> n: mediante el procedimeinto de cuantización vectorial
% En el primer caso hay que proporcionar el número de clases que se desean
% obtener y en el segundo hay que proporcionar un umbral T


% la variable nueva partición se establece para decidir si volver a
% realizar una nueva particón o no en función de los coeficientes de
% validación obtenido posterioremente.

NuevaParticion = true;

% se pide un número de clases, para distribuir las muestras en las
% respectivas clases, según el procedimeinto pseudoaleatorio incluido en
% la función fuzzy_clustering
%
while NuevaParticion
  n = 's';
  prompt = {'s/n'};
  dlg_title = 's -> Pseudo-aleatorio, n -> Cuantización Vectorial'; % s: fuzzy Clustering; n: Cuantización Vectorial 
  num_lines = 1;
  def = {'s'};
  answer = inputdlg(prompt,dlg_title,num_lines,def,options);
  n = answer{1};
  if (strcmp(n,'s') ~= 0)  % punto A
     % se ha elegido la opción de fijar el número de clústeres para dárselo
     % al algoritmo de Fuzzy clustering para realizar la partición de los
     % datos
     metodo_particion = 1; % método pseudoaleatorio de Balasko

     prompt = {'Número de Clases','Peso exponencial','Criterio de terminación','Número Máximo de iteraciones'};
     dlg_title = 'Párametros Agrupamiento Borroso (pseudoaleatorio)';
     num_lines = 1;
     def = {'3','2','1e-2','20'};
     answer = inputdlg(prompt,dlg_title,num_lines,def,options);
     c = str2num(answer{1});
     m = str2num(answer{2});
     e = str2num(answer{3});
     iterMax = str2num(answer{4});
     
     % los datos están todavía sin clasificar, por tanto los centros de las clases se
     % inicializarán de forma aleatoria
     
     vEstimados = zeros (1,3);  % no se usará realmente es sólo para concordancia de parámetros con
                                % la función posterior
     DatosClasificados = false; 
     [v,U, distancia, iteraciones, Jcost] = fuzzy_clustering (X,c,m,e,vEstimados,DatosClasificados,iterMax);

     %desnormalización de v
     v = (repmat(MA,size(v,1),1) - repmat(MI,size(v,1),1)).*v + ...
        repmat(MI,size(v,1),1);
    
     idx_clas = zeros(1,c);   
     [mx,nx] = size(X1);
     for hh=1:1:mx
       A = X1(hh,:);
       % calculo de los grados de pertenencia a cada una de las clases
       for i=1:1:c
          dd(i) = (1/sum((A - v(i,:)).^2))^(2/(m-1));
       end

       SUMA = sum(dd);
       for i=1:1:c
         mu(i) = dd(i)/SUMA;
       end

       [fila, columna, valor] = find(mu == max(mu));
        
       idx_clas (columna) = idx_clas(columna) + 1;
       j = idx_clas (columna);   
       DATOS_CLASES(columna).muestras(j,:) = A;   
       DATOS_CLASES(columna).numero_muestras = j; 
     end
  else %relativo al punto A
     % se ha elegido la opción de determinar el número de clústeres de forma
     % automática a través del algoritmo de cuantización vectorial
     metodo_particion = 2; % cuantización vectorial
     
     prompt = {'Umbral'};
     dlg_title = 'Umbral diferenciación de Clases';
     num_lines = 1;
     def = {'80'};
     answer = inputdlg(prompt,dlg_title,num_lines,def,options);
     T = str2num(answer{1});
     
     [N,n] = size(X);

     % inicialización
     DATOS_CLASES(1).centroCVectorial = X1(1,:);  
     DATOS_CLASES(1).muestras(1,:)    = X1(1,:);   
     DATOS_CLASES(1).numero_muestras = 1; %numero de patrones
     c = 1; %numero de clases

     for j=2:1:N % punto B: número de muestras
       % se computa la distancia del nuevo patron a los centros
       nueva_clase = false;
       for h=1:1:c
         d(h) = norm(X1(j,:)-DATOS_CLASES(h).centroCVectorial);
       end

       %calcular la distancia mínima a los centros de cada clase, con indicación de la clase (h) 
       dmin = min(d);
       h = find(d == dmin);

       % si esa distancia mínima es inferior a un umbral determinado, añadir el
       % patrón a la clase, de otro modo crear una nueva clase
       if dmin < T
          DATOS_CLASES(h).numero_muestras = DATOS_CLASES(h).numero_muestras + 1;
          DATOS_CLASES(h).muestras(DATOS_CLASES(h).numero_muestras,:) = X1(j,:); 
          % en este caso se recalculan los centros de la clase como el
          % valor medio
          DATOS_CLASES(h).centroCVectorial = mean(DATOS_CLASES(h).muestras);

       else % hay que crear una nueva clase 
          nueva_clase = true;
       end
  
       % creación de la nueva clase
       if nueva_clase
         c = c + 1;
         DATOS_CLASES(c).numero_muestras  = 1;
         DATOS_CLASES(c).muestras(1,:)    = X1(j,:);
         DATOS_CLASES(c).centroCVectorial = X1(j,:);
         
       end

     end % punto B
 end  %relativo al punto A if (strcmp(n,'s') ~= 0)
   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%  Paso 4: validación de la partición   %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
  % Generacion del vector X, de nuevo a partir de las clases
  % Se obtiene además una estima (vEstimados, media para cada clase)
  hh = 0;
  for ii=1:1:c
     muestras_clase = DATOS_CLASES(ii).numero_muestras;
     for jj=1:muestras_clase
        hh = hh + 1;
        X(hh,:) = DATOS_CLASES(ii).muestras(jj,:);
     end
     vEstimados(ii,:) = mean(DATOS_CLASES(ii).muestras); 
  end

  % normalizacion al rango [0,1]
  MI=min(X); %minimo
  MA=max(X); %maximo
  X=(X-repmat(MI,size(X,1),1))./(repmat(MA,...
     size(X,1),1)-repmat(MI,size(X,1),1));

  vEstimados = (vEstimados-repmat(MI,size(vEstimados,1),1))./(repmat(MA,... %para usarlos posteriormente
                size(vEstimados,1),1)-repmat(MI,size(vEstimados,1),1));
  
  % los datos han sido ya clasificados, por tanto los centros de las clases se
  % inicializarán a las medias de la clasificación obtenida
  if metodo_particion == 2
    % se hizo la partición mediante el método de cuantización vectorial y
    % es necesario pasar los parámetros del clasificador fuzzy, que no los
    % tenemos
     prompt = {'Peso exponencial','Criterio de terminación','Número Maximo iteraciones'};
     dlg_title = 'Párametros para el algoritmo borroso';
     num_lines = 1;
     def = {'2','1e-2','20'};
     answer = inputdlg(prompt,dlg_title,num_lines,def,options);
     m = str2num(answer{1});
     e = str2num(answer{2});
     iterMax = str2num(answer{3});
  end
  
  DatosClasificados = true; 
  [v,U, distancia, iteraciones, Jcost] = fuzzy_clustering (X,c,m,e,vEstimados,DatosClasificados,iterMax);
    
  CvalidezFuzzy = validezFuzzy(X,v,U,distancia,m,c);
  
  disp(' ');
  disp('El número de clases obtenido ha sido: '); disp(c);
  disp(' ');
  
  % se muestran por pantalla los diferentes coeficientes de validez Fuzzy
  disp('Coeficientes Fuzzy ....');
  disp('Coeficiente de partición   (CP):'); disp(CvalidezFuzzy.CP);
  disp('Coeficiente de entropía    (CE):'); disp(CvalidezFuzzy.CE);
  disp('Índice de partición        (SC):'); disp(CvalidezFuzzy.SC);
  disp('Índice de partición        (S) :'); disp(CvalidezFuzzy.S);
  disp('Índice de Xie-Beni         (XB):'); disp(CvalidezFuzzy.XB);
  disp('Índice de Dunn             (ID):'); disp(CvalidezFuzzy.ID);
  disp('Índice de Dunn Alternativo (IDA):'); disp(CvalidezFuzzy.IDA);
 
    %desnormalización de v
  v = (repmat(MA,size(v,1),1) - repmat(MI,size(v,1),1)).*v + ...
        repmat(MI,size(v,1),1);

  for i=1:1:c
    DATOS_CLASES(i).vFuzzy = v(i,:);
    DATOS_CLASES(i).UFuzzy = U(i,:);
  end

  % Clasificador Bayesiano 
  DATOS_CLASES = Bayesian(DATOS_CLASES, c);
  [Divergencia, Coseno_alpha, Jeffries_Matusita] = validezBayesian(DATOS_CLASES,c);
  
  disp('Coeficientes de Bayes ....');
  disp('Divergencia           :'); disp(Divergencia);
  disp('Coseno Alpha          :'); disp(Coseno_alpha);
  disp('Jeffries Matussita    :'); disp(Jeffries_Matusita);

  % En función de los coeficentes obtenidos previamente se pregunta sobre
  % una nueva aprtición. En este punto, dependiendo de si se ha optado por
  % el procedimiento de Balasko o por la cuantización vectorial se vuelve a
  % preguntar por los nuevos parámetros: 
  % 1) Balasko: m (peso exponencial, epsilon y número máximo de iteraciones)
  % 2) Cuantización vectorial: T (umbral)
  
  
  prompt = {'¿Nueva partición (s/n)?'};
  dlg_title = 'Intentar nueva partición';
  num_lines = 1;
  def = {'n'};
  answer = inputdlg(prompt,dlg_title,num_lines,def,options);
  n = answer{1};
 
  if (strcmp(n,'n') ~= 0)
    NuevaParticion = false; % no se intentan más particiones
  else
    % borramos lo anterior
    clear DATOS_CLASES;
  end
end % bucle while NuevaParticion  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Paso 5: alamacenar los datos en la Base de Conocimiento %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
save ParamAprendizaje DATOS_CLASES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Paso 6: representación de los datos y centros Fuzzy      %%%%%%%
%%%%%%%%          en 3 y 2 dimensiones                             %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot3(X1(:,1),X1(:,2),X1(:,3),'k.','LineWidth',3);
hold on

if c == 1
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
end

if c == 2
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
end

if c == 3
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
end

if c == 4
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
  plot3(v(4,1),v(4,2),v(4,3),'yd','LineWidth',5);
end

if c == 5
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
  plot3(v(4,1),v(4,2),v(4,3),'yd','LineWidth',5);
  plot3(v(5,1),v(5,2),v(5,3),'cd','LineWidth',5);
end
if c == 6
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
  plot3(v(4,1),v(4,2),v(4,3),'yd','LineWidth',5);
  plot3(v(5,1),v(5,2),v(5,3),'cd','LineWidth',5);
  plot3(v(6,1),v(6,2),v(6,3),'md','LineWidth',5);
end
if c == 7
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
  plot3(v(4,1),v(4,2),v(4,3),'yd','LineWidth',5);
  plot3(v(5,1),v(5,2),v(5,3),'cd','LineWidth',5);
  plot3(v(6,1),v(6,2),v(6,3),'md','LineWidth',5);
  plot3(v(7,1),v(7,2),v(7,3),'kd','LineWidth',5);
end
if c == 8
  plot3(v(1,1),v(1,2),v(1,3),'gd','LineWidth',5);
  plot3(v(2,1),v(2,2),v(2,3),'rd','LineWidth',5);
  plot3(v(3,1),v(3,2),v(3,3),'bd','LineWidth',5);
  plot3(v(4,1),v(4,2),v(4,3),'yd','LineWidth',5);
  plot3(v(5,1),v(5,2),v(5,3),'cd','LineWidth',5);
  plot3(v(6,1),v(6,2),v(6,3),'md','LineWidth',5);
  plot3(v(7,1),v(7,2),v(7,3),'kd','LineWidth',5);
  plot3(v(8,1),v(8,2),v(8,3),'wd','LineWidth',5);
end
grid on

  
q = 2; % proyección a la dimensión q (en este caso bidimensional)
[y,vp,V,D] = Componentes_principales(X1,v,q,c);
figure; plot(y(:,1),y(:,2),'k.','LineWidth',3)

if c == 1
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
end

if c == 2
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
end

if c == 3
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
end

if c == 4
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
  hold on; plot(vp(4,1),vp(4,2),'yd','LineWidth',5)   
end
if c == 5
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
  hold on; plot(vp(4,1),vp(4,2),'yd','LineWidth',5)   
  hold on; plot(vp(5,1),vp(5,2),'cd','LineWidth',5)   
end
if c == 6
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
  hold on; plot(vp(4,1),vp(4,2),'yd','LineWidth',5)   
  hold on; plot(vp(5,1),vp(5,2),'cd','LineWidth',5)   
  hold on; plot(vp(6,1),vp(6,2),'md','LineWidth',5)   
end
if c == 7
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
  hold on; plot(vp(4,1),vp(4,2),'yd','LineWidth',5)   
  hold on; plot(vp(5,1),vp(5,2),'cd','LineWidth',5)   
  hold on; plot(vp(6,1),vp(6,2),'md','LineWidth',5)   
  hold on; plot(vp(7,1),vp(7,2),'kd','LineWidth',5)   
end
if c == 8
  hold on; plot(vp(1,1),vp(1,2),'gd','LineWidth',5)
  hold on; plot(vp(2,1),vp(2,2),'rd','LineWidth',5)
  hold on; plot(vp(3,1),vp(3,2),'bd','LineWidth',5)
  hold on; plot(vp(4,1),vp(4,2),'yd','LineWidth',5)   
  hold on; plot(vp(5,1),vp(5,2),'cd','LineWidth',5)   
  hold on; plot(vp(6,1),vp(6,2),'md','LineWidth',5)   
  hold on; plot(vp(7,1),vp(7,2),'kd','LineWidth',5)   
  hold on; plot(vp(8,1),vp(8,2),'wd','LineWidth',5)   
end
grid on
