formapoligonal.m

function [F, T, gruposseparables, materias] = formapoligonal(nombrefuncion, msist, malet, maletodo, LB, UB, K, grado, potencia)
%formapoligonal
%   F = formapoligonal(funcion)
%
%   Obtiene unas recetas en tiempo discreto con la aproximaci�n poligonal a
%   unas funciones de producci�n (en tiempo discreto o continuo) 
%
%   [F, T] = formapoligonal(funcion, msist, malet, maletodo, LB, UB, K, grado)
%      funcion   archivo con la definici�n de las funciones de producci�n
%      msist     muestras sistem�ticas por cada insumo
%      malet     muestras aleatorias por grupo separable
%      maletodo  muestras aleatorias para todos los insumos
%      LB        l�mite inferior para cada insumo (por defecto 0)
%      UB        limite superior para cada insumo (por defecto 1)
%      K         factores con formapolinomioinverso (s�lo con tiempo continuo)
%      grado     grado con formapolinomioinverso (s�lo con tiempo continuo)
%      F         aproximaci�n poligonal como recetas en tiempo discreto
%      T         matriz de conversi�n de intensidades a insumos
%   msist y malet tambi�n pueden definirse como vectores para cada grupo
%   separable, y LB y UB para cada insumo. En Octave con funciones en
%   tiempo continuo hay que definir alg�n K.  
%
%   A partir de la X soluci�n de VN para F los insumos se obtienen con X*T. 
%   En el caso en tiempo discreto los flujos netos en el tiempo se obtienen 
%   con funcion(X*T); las materias se corresponden con las filas y el
%   tiempo con las columnas. En el caso en tiempo continuo los flujos
%   actualizados se obtienen con funcion(k, X*T). 

%necesita formapolinomioinverso.m

if nargin < 2, msist = []; end
if nargin < 3, malet = []; end
if nargin < 4, maletodo = []; end
if nargin < 5, LB = []; end
if nargin < 6, UB = []; end
if nargin < 7, K = []; end
if nargin < 8, grado = []; end
if nargin < 9, potencia = []; end

v = version;                                                 %determina la versi�n del software
if and(exist('OCTAVE_VERSION'), str2num(v([1,3])) < 36)      %con Octave en versiones inferiores a 3.6 (nargin necesita que la funci�n est� cargada)
   if isempty(K), numarg = 1; else numarg = 2; end               %define el caso en tiempo discreto si K no est� definido
else                                                         %si se usa Matlab u Octave 3.6 o superior
   numarg = nargin(nombrefuncion);                               %define el caso en tiempo discreto por el n�mero de argumentos
end                                                          %
                                                             %
switch numarg                                                %dependiendo de si el problema est� en tiempo discreto o continuo 
case 1                                                       %si el problema est� en tiempo discreto
   [F, gruposseparables] = feval(nombrefuncion, zeros(1,10000)); %obtiene la variable gruposseparables
   [materias, grado] = size(F);                                  %determina el n�mero de materias y de pasos temporales
case 2                                                       %si el problema est� en tiempo continuo
   [F, gruposseparables] = feval(nombrefuncion,1,zeros(1,10000));%obtiene la variable gruposseparables
   materias = size(F, 2);                                        %determina el n�mero de materias   
   if isempty(K), K = 1; end                                     %si K no est� definido es 1
   if size(K(:),1) == 1, K = (K/2:K/5:K*3/2); end                %si K es un escalar la aproximaci�n se efectua en [K/2, K*3/2] con 6 puntos
   puntos = size(K, 2);                                          %determina el n�mero de factores para los que se hace la aproximaci�n
   if isempty(grado), grado = fix((puntos+3)/2); end             %define el grado del polinomio inverso por defecto
end                                                          %
insumos = sum(gruposseparables);                             %establece el n�mero de insumos
numgrupos = size(gruposseparables, 2);                       %establece el n�mero de grupos de insumos separables
                                                             %
if isempty(msist), msist = 5; end                            %muestras sistem�ticas para cada insumo a variar en cada grupo separable por defecto
if isempty(malet), malet = 0; end                            %muestras aleatorias para cada grupo separable por defecto
if isempty(maletodo), maletodo = 0; end                      %muestras aleatorias para todos los insumos por defecto
if isempty(LB), LB = 0; end                                  %establece el l�mite inferior para cada insumo por defecto
if isempty(UB), UB = 1; end                                  %establece el l�mite superior para cada insumo por defecto
if isempty(potencia), potencia = 10; end                     %potencia a la que se elevan los insumos aleatorios por defecto
if size(msist,2) == 1, msist = msist.*ones(1,numgrupos); end %si las muestras sistem�ticas es un escalar lo aplica a todos los grupos
if size(malet,2) == 1, malet = malet.*ones(1,numgrupos); end %si las muestras aleatorias es un escalar lo aplica a todos los grupos
if size(LB, 2) == 1, LB = LB .* ones(1, insumos); end        %si el l�mite inferior es un escalar lo aplica a todos los insumos
if size(UB, 2) == 1, UB = UB .* ones(1, insumos); end        %si el l�mite superior es un escalar lo aplica a todos los insumos
UL = UB - LB;                                                %determina la amplitud de los l�mites
T = [];                                                      %define la variable para almacenar los resultados
                                                             %========= aproximaci�n sistem�tica por cada grupo de insumos separables
Xinicio = 1;                                                 %Xinicio es la variable que describe el primer insumo del grupo separable
for c1 = 1:numgrupos                                         %para cada grupo de insumos separables
    Xfin = Xinicio + gruposseparables(c1) - 1;                   %Xfin es la variable que describe el �ltimo insumo del grupo separable
    if msist(c1) > 0                                             %si queremos obtener alguna muestra del grupo separable
        V = varia(gruposseparables(c1),msist(c1)) ./ msist(c1);      %V es la variaci�n de msist muestras para los insumos en el grupo
        sV = size(V, 1);                                             %establece el n�mero de muestras
        Xn = zeros(sV, insumos);                                     %define Xn
        for c2 = 1:sV                                                %para cada muestra
            Xn(c2,Xinicio:Xfin)=LB(Xinicio:Xfin)+V(c2,:).*UL(Xinicio:Xfin);%ajusta a LB y UB
        end                                                          %
        T = [T; Xn];                                                 %almacena los insumos
    end                                                          %
    Xinicio = Xfin + 1;                                          %establece Xinicio para el pr�ximo grupo separable
end                                                          %
                                                             %========= aproximaci�n aleatoria por cada grupo de insumos separables
Xinicio = 1;                                                 %Xinicio es la variable que describe el primer insumo del grupo separable
for c1 = 1:numgrupos                                         %para cada grupo de insumos separables
    Xfin = Xinicio + gruposseparables(c1) - 1;                   %Xfin es la variable que describe el �ltimo insumo del grupo separable
    if malet(c1) > 0                                             %si queremos obtener alguna muestra del grupo separable
        V = rand(malet(c1), Xfin-Xinicio+1) .^ potencia;              %V es la variaci�n aleatoria de malet muestras para los insumos en el grupo
        Xn = zeros(malet(c1), insumos);                               %define Xn
        for c2 = 1:malet(c1)                                          %para cada muestra
            V(c2, :) = V(c2, :) ./ sum(V(c2, :));                          %normaliza para que sume 1
            Xn(c2,Xinicio:Xfin)=LB(Xinicio:Xfin)+V(c2,:).*UL(Xinicio:Xfin);%ajusta a LB y UB
        end                                                           %
        T = [T; Xn];                                                  %almacena los insumos
    end                                                          %
    Xinicio = Xfin + 1;                                          %establece Xinicio para el pr�ximo grupo separable
end                                                          %
                                                             %========= aproximaci�n aleatoria de todos los insumos
Xn = rand(maletodo, insumos) .^ potencia;                    %obtiene las muestras aleatorias para todos los inusmos
for c1 = 1:maletodo                                          %para cada muestra
   Xn(c1, :) = Xn(c1, :) ./ sum(Xn(c1, :));                      %normaliza cada muestra para que sume 1
   Xn(c1, :) = LB + Xn(c1, :) .* UL;                             %ajusta a LB y UB
end                                                          %
T = [T; Xn];                                                 %almacena los insumos
                                                             %========= calcula las recetas con los flujos netos
procesos = size(T, 1);                                       %establece el n�mero de procesos
F = zeros(procesos, materias, grado);                        %define las matrices de flujos netos
for c1 = 1:procesos                                          %para cada proceso
   switch numarg                                                 %dependiendo de si el problema est� en tiempo discreto o continuo 
   case 1, F(c1,:,:) = feval(nombrefuncion, T(c1, :));                           %establece los flujos netos para el caso en tiempo discreto
   case 2, F(c1,:,:) = formapolinomioinverso(nombrefuncion, K, grado, T(c1, :)); %aproxima los flujos netos para el caso en tiempo continuo
   end                                                           %
end                                                          %




function P = varia(columnas, n)
if columnas == 1
    P = n;
elseif columnas == 2     %no es estrictamente necesario
    P = zeros(n+1, 2);
    R = (0:n)';
    P(:, 1) = n - R;
    P(:, 2) = R;
else
    P = [];
    for c1 = 0:n
        P2 = varia(columnas-1, n-c1);
        q = size(P2, 1);
        P = [c1*ones(q,1), P2; P];
    end
end