function [tolmax,argmax] = tolinprog(infA,supA,infb,supb,varargin)
%  
%   Нахождение максимума распознающего функционала допускового множества 
%   решений для интервальной системы линейных алгебраических уравнений (ИСЛАУ)  
%   через решение задач линейного программирования (ЛП). 
%  
%   В зависимости от значения параметра lpmode процедура реализует либо способ 
%   нахождения максимума через решение одной задачи ЛП (которая формируется по 
%   рецепту С.И.Жилина  на основе теоремы И.Рона  о представлении допускового 
%   множества решений), либо через вариации уровня - частный случай метода 
%   дробления графика.  См. подробности этого подхода в книге 
%               С.П. Шарый, Конечномерный интервальный анализ.
%               – Новосибирск: XYZ, 2022, параграф 3.6.
%  
%   В реализованом ниже  алгоритме  "вариация уровня"  выполняется  неявно,
%   с помощью вариации правой части - либо её равномерного раздувания, либо
%   равномерного сужения. Контроль сильной совместности, т.е. непустоты или
%   пустоты допускового множества решений, обеспечивается решением системы
%   линейных неравеств, которая указывается теоремой И.Рона о представлении
%   допускового множества решений  ИСЛАУ  (J.Rohn, 1985).  Варьируя  правую
%   часть ИСЛАУ, мы нащупываем тот "уровень" на котором происходит  переход
%   от совместности ИСЛАУ  к её несовместности, т.е. максимум распознающего
%   функционала.
%  
%   Обязательные входные аргументы функции:
%        infA, supA - матрицы левых и правых концов интервальных коэффициентов 
%                     при  неизвестных   в   интервальной  системе  линейных 
%                     алгебраических уравнений; они могут быть прямоугольными, 
%                     но должны иметь одинаковые размеры; 
%        infb, supb - векторы левых и правых концов интервалов  правой части 
%                     интервальной системы линейных алгебраических уравнений. 
%  
%   Необязательные входные аргументы функции:
%            lpmode - задание способа использования задачи ЛП при нахождении 
%                     максимума распознающего функционала; 
%                     если lpmode = 0 (значение по умолчанию), то максимум  
%                     находится  посредством решения  одной  задачи  ЛП; 
%                     если  lpmode ~= 0,  то  максимум находится  с помощью 
%                     бисекции величины раздутия правой части, где задача ЛП 
%                     решается на каждом шаге процесса. 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   С.П. Шарый, ФИЦ ИВТ, НГУ, 2022 г. 
%   С.И. Жилин, "СиСорт", АлтГУ, 2022 г. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
lpmode = 0;         %   задание способа решения задачи максимизации 
IterLim = 100;      %   ограничение на количество бисекций
delta = 1.e-8;      %   порог точности нахождения максимума
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   проверка корректности входных данных
%
mi = size(infA,1);  ni = size(infA,2);
ms = size(supA,1);  ns = size(supA,2);
if mi==ms   %   m - количество уравнений в системе
    m = ms;
else
    error('Количество строк в матрицах левых и правых концов неодинаково')
end
if ni==ns
    n = ns; %   n - количество неизвестных переменных в системе
else
    error('Количество столбцов в матрицах левых и правых концов неодинаково')
end
  
ki = size(infb,1);
ks = size(supb,1);
if ki==ks
    k = ks;
else
    error('Количество компонент у векторов левых и правых концов неодинаково')
end
if k~=m
    error('Размеры матрицы системы не соответствуют размерам правой части')
end
  
if ~all(all(infA <= supA))
    error('В матрице системы задан неправильный интервальный элемент')
end

if ~all(infb <= supb)
    error('В векторе правой части задана неправильная интервальная компонента')
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%   переназначение параметров алгоритма, заданных пользователем 
% 
if nargin >= 5 
    lpmode = ceil(varargin{1}); 
end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~lpmode 
    %   находим максимум с помощью единственной задачи ЛП
    [argmax, tolmax] = SingleLPP; 
else
    %   находим максимум с помощью бисекции по alpha 
  
    %   находим интервал значений параметра раздутия alpha правой части ИСЛАУ,
    %   на котором максимум  распознающего функционала  допускового множества
    %   решений меняет знак, иными словами, сильная совместность превращается
    %   в несовместность 
  
    %   сначала проверяем значение alpha = 0, т.е. исходную ИСЛАУ
    [LEs, UEs] = LevStat(0);

    %   в зависимости от результата
    %       либо раздуваем правую часть ИСЛАУ до появления разрешимости,
    %       либо сжимаем её до исчезновения разрешимости
    if isnan(LEs)
        %   допусковое множество решений пусто,
        %   раздуваем правую часть ИСЛАУ
        incr = 0.5*max(supb - infb);    %   шаг раздутия
        alpha = 0.;
        do
            alpha = alpha + incr;
            [LEs, UEs] = LevStat(alpha);
            incr = 2*incr;
        until ~isnan(LEs)
        Lower = 0; Upper = alpha;
    else
        %   допусковое множество решений непусто,
        %   суживаем правую часть ИСЛАУ
        incr = 0.3*min(supb - infb);    %   шаг сужения
        alpha = 0.;
        do
            alpha = alpha - incr;
            [LEs, UEs] = LevStat(alpha);
            incr = 0.5*incr;
        until isnan(LEs)
        Lower = alpha; Upper = 0;
    endif
    
printf("\n Интервал параметра раздутия правой части = [%d, %d]\n",Lower,Upper)
    
    %   после того, как найден интервал [Lower, Upper], который содержит точку
    %   перехода  параметра  раздутия  alpha  через границу совместности ИСЛАУ, 
    %   уточняем её положение  с помощью процесса дихотомии (бисекции) 
    k = 0;
    while (Upper-Lower) > delta && k < IterLim
        k = k + 1;
        mid = 0.5*(Lower + Upper);
        [LEs, UEs] = LevStat(mid);
        if isnan(LEs)
            Lower = mid;
        else
            Upper = mid;
            LoArg = LEs;
            UpArg = UEs;
        endif
    endwhile

    %   формируем решение задачи
    argmax = 0.5*(LoArg + UpArg);
    %   можно оценить максимум  как середину двух уровней,  между 
    %   которыми произошёл переход от соместности к несовместности 
    %tolmax = -0.5*(Lower + Upper);

    %   точнее максимум оценивается через вычисление значения функционала
    tolmax = calcf(argmax); 
endif 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   вывод диагности решения задачи 
  
if lpmode 
    printf(" Количество бисекций величины раздутия правой части = %u\n",k)
endif
  
if tolmax >= 0
    printf("\n Допусковое множество решений интервальной системы непусто\n\n") 
else
    printf("\n Допусковое множество решений интервальной системы пусто\n\n") 
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [LEs,UEs] = LevStat(t)
%
%   Функция, определяющая пустоту или непустоту допускового множества решений
%   (т.е. сильную совместность)  интервальной  линейной  системы  уравнений
%   с правой частью, равномерно "раздутой" на величину t, а также его оценки
%   по выделенной координате. Значение t может быть также отрицательным, что
%   соответствует равномерному сужению правой части.
%
%   Решение системы линейных неравенств, определяющее допусковое множество
%   решений в силу теоремы И.Рона, а также его оценивание выполняется с помощью
%   стандартной функции glpk из состава системы Octave, которая предназначена
%   для решения задачи линейного программирования
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %   формируем матрицу  и правую часть  расширенной системы  линейных
    %   неравенств (из теоремы И.Рона), определяющую допусковое множество
    %   решений ИСЛАУ с "раздутой" интервальной правой частью
    AA = [ supA, -infA; -infA, supA ];
    bb = [ supb + t; -infb + t ];

    %   формируем вектор градиента целевой функции задачи ЛП для оценки
    %   первой компоненты точек из допускового множества решений
    c = zeros(2*n,1);   c(1) = 1;   c(n+1) = -1;

    %   формируем остальные параметры функции glpk, описывающие
    %   тип задачи линейного программирования и опции её решения
    lb = zeros(2*n,1);              %   определяем нижнюю границу переменных
    ub = [];                        %   определяем верхнюю границу переменных
    ctype = repmat("U",[1,2*m]);    %   определяем тип ограничений
    vartype = repmat("C",[1,2*n]);  %   определяем тип переменных -
    sense = -1;                     %   определяем тип вычисляемого экстремума
    param.msglev = 0;               %   определяем тип диагностики
    %param.itlim = 100;              %   ограничиваем число итераций

    %   решаем задачу линейного прграммирования по оцениванию максимума
    %   первой компоненты  точек  из допускового  множества решений,
    %   если, конечно, система линейных неравенств из теоремы И.Рона
    %   разрешима, т.е. множество допустимых решений задачи ЛП непусто
    [xopt, fopt, errnum, extra] = ...
                glpk (c, AA, bb, lb, ub, ctype, vartype, sense, param);
    %   запоминаем точку оптимума
    UEs = xopt(1:n) - xopt((n+1):(2*n));
  
    %   если множество допустимых решений задачи ЛП пусто (errnum == 10),
    %   т.е. система линейных неравенств из теоремы И.Рона несовместна,
    %   то интервал оценки по первой компоненте полагается пустым
    if errnum ~= 0
        if errnum == 10
            LEs = NaN;
            UEs = NaN;
        else 
            warning(" Нетипичное завершение функции glpk ") 
        endif
        %   можно сразу выйти и больше не оценивать множество решений
        return
    endif
  
    %   если  множество  допустимых  решений  задачи ЛП  непусто, то
    %   решаем задачу оценивания минимума первой координаты его точек
    sense = 1;
    [xopt, fopt, errnum, extra] = ...
                glpk (c, AA, bb, lb, ub, ctype, vartype, sense);
    %   запоминаем точку оптимума
    LEs = xopt(1:n) - xopt((n+1):(2*n));
endfunction
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function [f,tt] = calcf(x)
%
%   Функция, которая вычисляет значение f  максимизируемого распознающего
%   функционала в точке x и выдаёт также вектор tt из значений образующих
%   функционала в данной точке
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   предварительное размещение рабочих массивов
  
infs = zeros(m,1);
sups = zeros(m,1);
tt = zeros(m,1);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   вычисляем значение распознающего функционала и матрицу dd,
    %   составленную из суперградиентов его образующих
  
    for i = 1:m
        los = 0.5*(infb(i)+supb(i));
        sus = los;
        for j = 1:n
            if x(j) >= 0
                los = los - supA(i,j)*x(j);
                sus = sus - infA(i,j)*x(j);
            else
                los = los - infA(i,j)*x(j);
                sus = sus - supA(i,j)*x(j);
            end
        end
        %   нижний и верхний концы интервала, который получается
        %   под модулем в выражении для i-ой образующей
        infs(i) = los;
        sups(i) = sus;
    end
  
    %   вычисление значения i-ой образующей распознающего функционала
    for i = 1:m
        alos = abs(infs(i));  asus = abs(sups(i));
        if alos ~= asus
            if alos < asus
                mags = asus;
            else
                mags = alos;
            end
        else
            mags = alos;
        end
        %   нахождение и запоминание значения i-ой образующей
        tt(i) = 0.5*(supb(i) - infb(i)) - mags;
    end
  
    f = min(tt);
  
endfunction
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function [Xopt, topt] = SingleLPP
    %   формируем  специальную задачу  линейного программирования 
    %   с переменными  (x_1, ..., x_n, x_{n+1}, ..., x_{2n}, t) 
    AA = [ supA, -infA -ones(m,1);
          -infA,  supA -ones(m,1) ];
    bb = [ supb;
          -infb ];
  
    %   значение целевой функции совпадает со значением переменной t 
    c = [ zeros(1,2*n,1) 1 ];
  
    lb = [ zeros(1,2*n) -Inf ];     %   определяем нижнюю границу переменных
    ub = Inf * ones(1,2*n+1);       %   определяем верхнюю границу переменных
    ctype = repmat("U",[1,2*m]);    %   определяем тип ограничений
    vartype = repmat("C",[1,2*n+1]);%   определяем тип переменных -
    sense = 1;                      %   определяем тип вычисляемого экстремума
                                    %                        (max: -1, min: 1)
    param.msglev = 1;               %   определяем тип диагностики

    [xopt, fopt, errnum, extra] = ...
                glpk (c, AA, bb, lb, ub, ctype, vartype, sense, param);

    %   запоминаем точку оптимума
    Xopt = xopt(1:n) - xopt(n+1:2*n);
    topt = -xopt(2*n+1); % можно и как topt = -fopt;

    %   показываем дополнительные сведения о решении
    %   (двойственные переменные, время решения и т.д.)
    % extra
endfunction
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
endfunction
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  