vn.m

function [k, X, Y, pasos, distancias, sal1, sal2] = vn(F, signoX, signoXF, algoritmo, par1, par2, par3)
%vn
%   [k, X, Y] = vn(F)
%
%   Resuelve el modelo Von Neumann
%      Max k
%      X F(k) ~ 0      <->     -Y ~ 0               X F(k) .* Y = 0
%      X ~ 0           <->     F(k) Y ~ 0           X .* F(k) Y = 0
%      k > 0                   1 + X F'(k) Y = 0
%
%   [k, X, Y, pasos, distancias] = vn(F, signoX, signoXF, algoritmo)
%      F          recetas en tiempo discreto o continuo
%      signoX     signo de las intensidades; por defecto X >= 0
%      signoXF    signo de los balances materiales; por defecto X F(k) = 0
%      algoritmo  especifica el algoritmo que se usar�; por defecto 0
%      k          factor de expansi�n
%      X          intensidades-VN
%      Y          valores-VN o precios
%      pasos      n�mero de iteraciones efectuadas
%      distancias distancia de la soluci�n a las condiciones
%
%   Los c�digos para los signos son   0 =    1 >=    2 <=    3 >=<
%   Si se escribe un escalar para signoX o signoXF se aplica a todos; as�
%     vn(F, 1, 0) resuelve la forma est�ndar, X >= 0, X F(k) = 0.
%     vn(F, 1, 1) resuelve la forma can�nica, X >= 0, X F(k) >= 0.
%
%   Seg�n definamos algoritmo se usar�:
%      0   selecci�n autom�tica (por defecto)
%      1   vnautovalor
%      2   vnautovalor mayor orden
%      3   vnbiseccion
%      4   vnnewton
%      5   vnpls
%      6   vnsimplex
%      7   vnnolineal
%      8   vnbrody
%      9   cpbusqueda -> vnbiseccion
%     10   cpbusqueda -> vnpls
%     11   vnpls; cpbusqueda -> vnpls; cpbusqueda -> vnbiseccion
%
%   Cuando F son recetas en tiempo discreto la selecci�n autom�tica escoge
%   el algoritmo 1 para los problemas peque�os, el 2 para los problemas
%   casi-cuadrados medianos y para los problemas cuadrados grandes en la
%   forma est�ndar, el 11 para los problemas rectangulares grandes, y el 6
%   para los muy grandes. 
%   Cuando F son unas recetas en tiempo continuo la selecci�n autom�tica 
%   escoge el 11 en cualquier caso. Si con unas recetas en tiempo continuo
%   pedimos la soluci�n con los algoritmos 1,2,6,7,8 primero se usa 
%   formapolinomioinverso. 
%   Los algoritmos 9 a 11 ejecutan sucesivamente los algoritmos indicados
%   hasta encontrar una soluci�n.

%necesita forma0.m, formaproblema.m, formasolucion.m, condicionesvn.m,
%forma1.m, cpbusqueda.m, vnautovalor.m, vnbiseccion.m, vnnewton.m, vnpls.m,
%vnnolineal.m, vnbrody.m, vnsimplex

if nargin < 2, signoX = []; end
if nargin < 3, signoXF = []; end
if nargin < 4, algoritmo = []; end
if nargin < 5, par1 = []; end
if nargin < 6, par2 = []; end
if nargin < 7, par3 = []; end

k = []; X = []; Y = []; pasos = []; distancias = []; sal1 = []; sal2 = [];         %define las variables para la sailda
if isempty(F), disp('Los procesos de producci�n no est�n definidos'); return, end  %avisa de que F = []
[filas, columnas, ancho] = size(F);                           %establece el tama�o del problema
if ancho == 1, [filas, columnas] = size(forma0(F, 1)); end    %si las recetas est�n en tiempo continuo, establece el tama�o
if isempty(signoX), signoX = 1; end                           %por defecto X >= 0
if size(signoX, 2) == 1, signoX = signoX.*ones(1,filas); end  %si signoX es un escalar lo aplica a todas las intensidades
if isempty(signoXF), signoXF = 0; end                         %por defecto X F(k) = 0
if size(signoXF, 2) == 1, signoXF = signoXF.*ones(1,columnas); end                 %si signoXF es un escalar lo aplica a todos los balances materiales
if isempty(algoritmo), algoritmo = 0; end                     %si algoritmo no est� definido selecci�n autom�tica

simplex = (algoritmo == -1);                                  %si algoritmo == -1 lanza para simplex
if algoritmo <= 0                                             %si seleccion automatica o lanzado desde simplex
    if ancho == 1                                                 %si las recetas est�n definidas inline
        algoritmo = 11;                                               %algoritmo 11: vnpls; cpbusqueda -> vnpls; cpbusqueda -> vnbiseccion
    else                                                          %si las recetas est�n definidas como matrices
        filasestandar = filas - sum(signoX == 0) + sum(signoXF == 1) + sum(signoXF == 2) + 2*sum(signoXF == 3); %calcula el n�mero de filas en la forma est�ndar
        if (filas*columnas) < 50                                      %si problema peque�o
            if simplex, par1 = []; par2 = 1; end                          %si desde simplex, busca en todos los menores y una soluci�n
            algoritmo = 1;                                                    %algoritmo 1: vnautovalor
        elseif and(abs(filasestandar-columnas) < 2, columnas < 100)   %si problema casi-cuadrado mediano (en forma est�ndar)
            if simplex, par1 = 1; end                                     %si desde simplex, busca una solucion
            algoritmo = 2;                                                    %algoritmo 2: vnautovalor mayor orden
        elseif and(filasestandar==columnas, columnas < 500)           %si problema cuadrado grande (en forma est�ndar)
            if simplex, par1 = 1; end                                     %si desde simplex, busca una solucion
            algoritmo = 2;                                                    %algoritmo 2: vnautovalor mayor orden
        elseif and(filas < 3000, columnas < 200)                      %si problema rectangular grande
            algoritmo = 11;                                                   %algoritmo 11: vnpls; cpbusqueda -> vnpls; cpbusqueda -> vnbiseccion
        else                                                          %si problema muy grande
            if simplex                                                    %si desde simplex
                algoritmo = 11;                                               %algoritmo 11: vnpls; cpbusqueda -> vnpls; cpbusqueda -> vnbiseccion
            else                                                          %en caso contrario
                algoritmo = 6;                                                %algoritmo 6: vnsimplex
            end
        end
    end
end

if and(ancho == 1, sum(algoritmo==[1,2,6,7,8]))
    F = formapolinomioinverso(F);
    disp('Matrices de flujos netos aproximadas con formapolinomioinverso');
end

switch algoritmo
    case 1
        disp('Algoritmo 1: vnautovalor');
        [k, X, Y, pasos] = vnautovalor(F, signoX, signoXF, par1, par2);
        [X, Y, k] = normalizacion(k, X, Y, F, 1);
        [k, X, Y] = quitarsolucionesrepetidas(k, X, Y);
    case 2
        disp('Algoritmo 2: vnautovalor mayor orden');
        [Fn, signoX, signoXF] = formaproblema(F, signoX, signoXF);  %convierte a la forma estandar
        [k, X, Y, pasos] = vnautovalor(Fn, 1, 0, 0, par1);          %busca en los menores de mayor orden
        [X, Y] = formasolucion(X, Y, signoX, signoXF);              %convierte desde la forma estandar
        [X, Y, k] = normalizacion(k, X, Y, F, 1);
        [k, X, Y] = quitarsolucionesrepetidas(k, X, Y);
    case 3
        disp('Algoritmo 3: vnbiseccion');
        [k, X, Y, pasos] = vnbiseccion(F, signoX, signoXF, par1, par2);
        [X, Y] = normalizacion(k, X, Y, F);
    case 4
        disp('Algoritmo 4: vnnewton');
        [k, X, Y, pasos] = vnnewton(F, signoX, signoXF, par1);
        [X, Y] = normalizacion(k, X, Y, F);
    case 5
        disp('Algoritmo 5: vnpls');
        [k, X, Y, pasos] = vnpls(F, signoX, signoXF, par1, par2);
    case 6
        disp('Algoritmo 6: vnsimplex');
        [k, X, Y, pasos, sal2] = vnsimplex(F, signoX, signoXF, par1);
    case 7
        disp('Algoritmo 7: vnnolineal');
        [k, X, Y, pasos, sal1] = vnnolineal(F, signoX, signoXF, par1, par2);
        [X, Y] = normalizacion(k, X, Y, F);
    case 8
        disp('Algoritmo 8: vnbrody');
        [k, X, Y] = vnbrody(F, signoX, signoXF);
        [X, Y] = normalizacion(k, X, Y, F);
    case 9
        [k, Xsi, Ysi, u, pasos, kno, Xno, Yno] = cpbusqueda(F, signoX, signoXF, par1);
        [X, Y] = normalizacion(k, Xsi, Yno, F);
        if condicionesvn(k, X, Y, F, signoX, signoXF)
            disp('Algoritmo 9: cpbusqueda');
        else
            [k, X, Y, pasos2] = vnbiseccion(F, signoX, signoXF, k, kno);
            pasos = [pasos, pasos2];
            [X, Y] = normalizacion(k, X, Y, F);
            disp('Algoritmo 9: cpbusqueda -> vnbiseccion');
        end
    case 10
        [k, X, Y, u, pasos] = cpbusqueda(F, signoX, signoXF, par1);
        [X, Y] = normalizacion(k, X, Y, F);
        if condicionesvn(k, X, Y, F, signoX, signoXF)
            disp('Algoritmo 10: cpbusqueda');
        else
            [k, X, Y, pasos2] = vnpls(F, signoX, signoXF, k, X);
            pasos = [pasos, pasos2];
            disp('Algoritmo 10: cpbusqueda -> vnpls');
        end
    case 11
        k0 = par1;
        X = par2;
        [k1, X, Y, pasos] = vnpls(F, signoX, signoXF, k0, X);
        if condicionesvn(k1, X, Y, F, signoX, signoXF)
            disp('Algoritmo 11: vnpls');
            k = k1;
        else
            [ksi, Xsi, Ysi, usi, pasos2] = cpbusqueda(F, signoX, signoXF, k0);
            pasos = [pasos, pasos2];
            [X, Y] = normalizacion(ksi, Xsi, Ysi, F);
            if condicionesvn(ksi, X, Y, F, signoX, signoXF)
                disp('Algoritmo 11: vnpls; cpbusqueda');
                k = ksi;
            else
                [k3, X, Y, pasos3] = vnpls(F, signoX, signoXF, ksi, Xsi);
                pasos = [pasos, pasos3];
                if condicionesvn(k3, X, Y, F, signoX, signoXF)
                    disp('Algoritmo 11: vnpls; cpbusqueda -> vnpls');
                    k = k3;
                else
                    [kw, Xw, Yw, uw, pasos4, kno, Xno, Yno] = cpbusqueda(F, signoX, signoXF, ksi, 1); %busca un kno
                    pasos = [pasos, pasos4];
                    [X, Y] = normalizacion(ksi, Xsi, Yno, F);
                    if condicionesvn(ksi, X, Y, F, signoX, signoXF)
                        disp('Algoritmo 11: vnpls; cpbusqueda -> vnpls; cpbusqueda');
                        k = ksi;
                    else
                        [k, X, Y, pasos5] = vnbiseccion(F, signoX, signoXF, ksi, kno);
                        pasos = [pasos, pasos5];
                        [X, Y] = normalizacion(k, X, Y, F);
                        disp('Algoritmo 11: vnpls; cpbusqueda -> vnpls; cpbusqueda -> vnbiseccion');
                    end
                end
            end
        end
    otherwise
        disp('Atenci�n, no hay ning�n algoritmo con ese n�mero');
end

numsoluciones = size(k, 2);
distancias = zeros(filas+columnas+1, 3, numsoluciones);
if numsoluciones == 0
    disp('Atenci�n, no se ha encontrado ninguna soluci�n');
else
    for c1 = 1:numsoluciones
        [solucion, distancia, distanciatotal] = condicionesvn(k(c1), X(c1,:), Y(:,c1), F, signoX, signoXF);
        if nargout > 4, distancias(:,:,c1) = distancia; end
        if ~solucion
            if numsoluciones > 1, q = [' ', num2str(c1), ' ']; else q = ' '; end
            disp(['Atenci�n, la distancia de la soluci�n', q, 'a las condiciones de m�ximo, ', num2str(distanciatotal) , ', supera la prefijada.']);
        end
    end
end










function [X, Y, k] = normalizacion(k, X, Y, F, tipo)
if nargin < 5, tipo = 0; end
v = version;
for c1 = 1:size(k, 2)
    if and(tipo == 1, v([1,2]) == '5.')  %si vnautovalor y Matlab 5.0 ajustar a n�meros reales      
        w = max(abs(imag([k(c1), X(c1,:), Y(:,c1)'])));
        if w > 1e-10, disp(['Atenci�n, la soluci�n ', num2str(c1) , ' conten�a t�rm�nos complejos de magnitud ', num2str(w), '.']); end
        k(c1) = real(k(c1)); X(c1,:) = real(X(c1,:)); Y(:,c1) = real(Y(:,c1)); 
    end 
    s = sum(abs(X(c1,:)));
    if s > 0
        X(c1,:) = X(c1,:) ./ s;
        Y(:,c1) = Y(:,c1) .* s;
    end
    if tipo == 0        %si no vnautovalor
       q = -X(c1,:) * forma1(F, k(c1)) * Y(:,c1);
       if abs(q) > 0, Y(:,c1) = Y(:,c1) ./ q; end
    end
end




function [k, X, Y] = quitarsolucionesrepetidas(kn, Xn, Yn)
if ~isempty(kn)
    k = kn(1); X = Xn(1, :); Y = Yn(:, 1);
    for c1 = 2:size(kn, 2)
        poner = 1;
        for c2 = 1:size(k, 2)
            if sum(abs([kn(c1), Xn(c1,:), Yn(:,c1)'] - [k(c2), X(c2,:), Y(:,c2)'])) < 1e-12
                poner = 0; break;
            end
        end
        if poner, k = [k, kn(c1)]; X = [X; Xn(c1,:)]; Y = [Y, Yn(:,c1)]; end
    end
else
    k = []; X = []; Y = [];
end