function [x,h,f]=Simplexe(c,A,b,op)
%
% [x,h,f] = Simplexe(c,A,b,op)
%
%        Usa el método Simplex de las dos fases, para resolver cualquiera de los
%        siguientes problemas estandarizados (y sus duales):
%
%        op = 0 (Valor por omisión)                 op = 1
%        Problema de maximización:                  Problema de minimización:
%
%        Max    f = c*x                             Min    f = c*x
%
%        s.a.   A*x <= b                            s.a.   A*x >= b  
%                 x >= 0                                     x >= 0
%
%        Los valores de salida son:
%
%        x  = Vector solución al problema de maximización o minimización
%        h  = Vector de holguras del problema de maximización o minimización
%        f  = Valor óptimo de la función
%
% Simplexe(c,A,b) por si mismo, produce únicamente el vector solución al problema de
%                 maximización
%
%        En la entrada, deben darse los siguientes datos:
%
%        c  = Vector renglón con los coeficientes de la función a maximizar o minimizar
%        A  = Matriz de coeficientes de las restricciones, escritas en forma estándar
%        b  = Vector columna con los lados derechos de las restricciones
%        op = Opción que opera del siguiente modo:
%             op = 0, x y h corresponden al problema de maximización (valor por omisión)
%             op = 1, x y h corresponden al problema de minimización
%             op = 2, se muestran en la pantalla las etapas del algoritmo
%             (op puede considerar la suma de las opciones anteriores)
%
% Cuando no es posible calcular la solución, Simplexe detecta problemas NO FACTIBLES
% y problemas NO ACOTADOS, emitiendo el diagnóstico que corresponde.
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% Ver funciones relacionadas:
%
% Simplex, Simplexd, Simplexg
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%    Revisión 1.0
%    Sergio A. Pérez Ruiz
%    Instituto Tecnológico Autónomo de México
%    noviembre de 1996
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

%
% Se declara global el arreglo de trabajo (Tabla Símplex), para permitir
% interacciones posteriores con la Tabla Final
global W
%
% Valores por omisión
if nargin < 4, op=0; end
%
% VALIDACION Y DIAGNOSTICO DE LAS OPCIONES
err=0;
if op < 0 | op > 3
   disp('>>> Opción invalida (0 <= op <= 3)')
   disp('    op = 0, x y h corresponden al problema de maximización')
   disp('    op = 1, x y h corresponden al problema de minimización')
   disp('    op = 2, se muestran en la pantalla las etapas del algoritmo')
   disp('    (op puede considerar la suma de las opciones anteriores)')
   err=1;
end
%
% ANALISIS DE LAS OPCIONES
maxi=1; mini=0; show=0;
lop = op;
if lop-2 >= 0
   lop = lop - 2;
   show = 1;
end;
if lop-1 >= 0
   lop = lop - 1;
   maxi = 0; mini = 1;
end;
%
% PREPARACION DE LOS DATOS DE ACUERDO CON LA OPCION
% (SE TRANSPONEN LOS DATOS DE ENTRADA CUANDO SE REQUIERE LA SOLUCION
%  DEL PROBLEMA DUAL, DE MINIMIZACION)
if maxi
   cc=c;  AA=A;  bb=b;
end
if mini
   cc=b'; AA=A'; bb=c';
end
%
% DIMENSIONAMIENTO DEL PROBLEMA
[mc,nc]=size(cc); [mb,nb]=size(bb); [m,n]=size(AA);
%
% VALIDACION Y DIAGNOSTICO DE LOS PARAMETROS
if mc > 1
   disp('>>> Los coeficientes de la función deben ser vector renglón')
   err=1;
end
if nb > 1
   disp('>>> El lado derecho de las restricciones debe ser vector columna')
   err=1;
end
if mb ~= m | nc ~= n
   disp('>>> Matriz de restricciones mal dimensionada')
   err=1;
end
if err == 1
   disp('>>> POR FAVOR VERIFIQUE SUS PARAMETROS')
   return
end
%
% LOS DATOS DE ENTRADA FUERON ACEPTADOS E INTERPRETADOS
%
% VALORES INICIALES Y APUNTADORES
np1=n+1; mp1=m+1; mpn=m+n; mpnp1=mpn+1;
x=zeros(mpn,1); y=zeros(mpn,1); ix=[np1:mpn]'; fmax=0;
%
% CONSTRUYE LA TABLA SIMPLEX INICIAL
W=[
   AA  eye(m)      bb
   cc  zeros(1,m)   0
                     ];
W=(~(abs(W)<100*eps)).*W;
if show
   clc;
   if maxi
      disp('TABLA SIMPLEX INICIAL (MAX)')
   else
      disp('TABLA SIMPLEX INICIAL (MIN)')
   end      
   disp(W)
   disp('Para continuar omprima <enter>'); pause;
end
%
% CARACTERIZA EL ESTAR FUERA O DENTRO DEL CONJUNTO FACTIBLE
outcf=W(1:m,mpnp1)<0;
% ==================
% COMIENZA LA FASE I
% ==================
imq=0; colpiv=0; ishead=1; kc=1;
while max(outcf) ~= 0
   if show & ishead
      disp(' '); disp('EVOLUCION DE LA FASE I')
      ishead=0;
   end
%
%  IDENTIFICA UN RENGLON CON ELEMENTO NEGATIVO EN EL LADO DERECHO
   k=1;
   while k <= m
      if W(k,mpnp1)<0
         colpiv=0; imq=0;
%
%        BUSCA UN ELEMENTO NEGATIVO EN EL RENGLON
         j=1;
         while j <= mpn
            if W(k,j)<0
               imq=0; colpiv=j; mq=1/eps;
%
%              BUSCA EL COCIENTE MINIMO
               for i=1:m
                  if W(i,mpnp1)*W(i,colpiv) > 0
                     mq1=W(i,mpnp1)/W(i,colpiv);
                     if mq1 < mq
                        imq=i; mq=mq1;
                     end
                  end
               end
%              
%              SE LOCALIZO UN PIVOTE EN LA POSICION (imq,colpiv)
%              ELIMINACION GAUSSIANA EN COLUMNA colpiv
               W=eg(W,imq,colpiv);
%
%              PREVIENE LOOP INFINITO
               cycle(:,kc)=W(1:m,mpnp1);
               jc=1;
               while jc <= kc-1
                  if abs(cycle(:,jc)-W(1:m,mpnp1)) < 100*eps
                     imq=0; k=m;
                  end
                  jc=jc+1;
               end
               kc=kc+1;
%              
%              IDENTIFICA LA NUEVA VARIABLE BASICA x(colpiv)
               if imq ~= 0, ix(imq)=colpiv; end
               j=mpn; k=m;
            end   
            j=j+1;
         end
      end
      k=k+1;
   end
%              
%  PUNTOS DE SALIDA EN LA FASE I
   W=(~(abs(W)<100*eps)).*W;
   if show
      disp(W)
      disp('Para continuar omprima <enter>'); pause;
   end
%       
%  >PROBLEMA NO FACTIBLE
%   (NO HAY CANDIDATOS PARA PIVOTES)
   if imq == 0
      if maxi
         disp('>>> PROBLEMA NO FACTIBLE: La Fase I no puede continuar')
      else
         disp('>>> PROBLEMA DUAL NO FACTIBLE: La Fase I no puede continuar')
      end
      x=[]; y=[]; f='NF';
      return
   end
   outcf=W(1:m,mpnp1)<0;
end
% ===================
% COMIENZA LA FASE II
% ===================
% ITERACION PRINCIPAL
% PARA MOVERSE A UN PUNTO ESQUINA ADYACENTE E INCREMENTAR LA FUNCION
if show
   disp(' '); disp('EVOLUCION DE LA FASE II')
end
true=1;
while true == 1
%
%  ANALIZA LOS INDICADORES PARA BUSCAR UN PIVOTE
%  PUEDEN PRESENTARSE DOS SITUACIONES:
%  1) colpiv = 0, CUANDO YA NO HAY INCICADORES POSITIVOS
%  2) colpiv > 0, PERO imq = 0, CUANDO YA NO ES POSIBLE PIVOTEAR
   colpiv=0;
   j=1;
   while j <= mpn
      if W(mp1,j)>0
         colpiv=j;
%        
%        LOCALIZA UN PIVOTE Y PREPARA EL RENGLON CORRESPONDIENTE
         imq=0; mq=1/eps;
         for i=1:m
            if W(i,colpiv) > 0
               mq1=W(i,mpnp1)/W(i,colpiv);
               if mq1 < mq
                  imq=i; mq=mq1;
               end
            end
         end
%        
%        AHORA SE TIENE UN PIVOTE EN LA POSICION (imq,colpiv)
         if imq > 0 
            j=mpn;
         end
      end
      j=j+1;
   end
%
%  PUNTOS DE SALIDA EN LA FASE II
   W=(~(abs(W)<100*eps)).*W;
   if show
      disp(W)
      disp('Para continuar omprima <enter>'); pause;
   end
%
%  >TERMINACION NORMAL DEL METODO SIMPLEX
%  (SE AGOTARON LOS INDICADORES POSITIVOS)
   if colpiv == 0
%
%     ACTUALIZA EL PUNTO ESQUINA EN EL PASO ACTUAL, PARA EL PROBLEMA DE
%     MAXIMOS Y PARA EL DE MINIMOS, Y EL VALOR DE LA FUNCION OBJETIVO
      x=zeros(mpn,1);
      for i=1:m
         if ix(i) ~= 0
             x(ix(i))=W(i,mpnp1);
         end
      end;
      y=-[W(mp1,np1:mpn) W(mp1,1:n)]';
      if maxi
         f=cc*x(1:n);
      end
      if mini
         f=bb'*y(1:m);
      end
%
%     AHORA, LOS VECTORES x E y CONTIENEN LAS SOLUCIONES A LOS PROBLEMAS
%     DE MAXIMIZACION Y MINIMIZACION, RESPECTIVAMENTE, E INCLUYEN LAS HOLGURAS
%     AL FINAL DE CADA UNO
      if maxi
         h=x(np1:mpn);
         x=x(1:n);
      end
      if mini
         h=y(mp1:mpn);
         x=y(1:m);
      end
      return
   end
%
%  >PROBLEMA NO ACOTADO
%  (AUN HAY INDICADORES POSITIVOS, PERO NO HAY CANDIDATOS PARA PIVOTES)
   if imq == 0
      if maxi
         disp('>>> PROBLEMA NO ACOTADO: La Fase II no puede continuar')
      else
         disp('>>> PROBLEMA DUAL NO ACOTADO: La Fase II no puede continuar')
      end
      x=[]; h=[]; f='NA';
      return
   end
%
%  NUEVA VARIABLE BASICA x(colpiv)
%  ELIMINACION GAUSSIANA EN COLUMNA colpiv
   W=eg(W,imq,colpiv);
%
%  IDENTIFICA LA NUEVA VARIABLE BASICA
   ix(imq)=colpiv;
%
%  ACTUALIZA EL PUNTO ESQUINA EN EL PASO ACTUAL, PARA EL PROBLEMA DE
%  MAXIMOS Y PARA EL DE MINIMOS, Y EL VALOR DE LA FUNCION OBJETIVO
   x=zeros(mpn,1);
   for i=1:m
      if ix(i) ~= 0
          x(ix(i))=W(i,mpnp1);
      end
   end;
   y=-[W(mp1,np1:mpn) W(mp1,1:n)]';
      if maxi
         f=cc*x(1:n);
      end
      if mini
         f=bb'*y(1:m);
      end
%
%  AHORA, LOS VECTORES x E y CONTIENEN LAS SOLUCIONES A LOS PROBLEMAS
%  DE MAXIMIZACION Y MINIMIZACION, RESPECTIVAMENTE, E INCLUYEN LAS HOLGURAS
%  AL FINAL DE CADA UNO
   if maxi
      h=x(np1:mpn);
      x=x(1:n);
   end
   if mini
      h=y(mp1:mpn);
      x=y(1:m);
   end
end
return