function [amax,umax,envs,ccod] = lintregs(A,b,varargin) 
%   Построение линейной интервальной регрессии. 
%   LINTREGS(A,b) выдаёт  значения  коэффициентов  линейной  регрессионной
%   зависимости  по матрице  входных значений A  и  интервальному  вектору 
%   выходных данных b путём вычисления максимума распознающего функционала
%   множества решений интервальной линейной системы Ax = b. 
%  
%   Синтаксис вызова:
%       [amax,umax,ccod] = lintregs(A,b,iprn,epsf,epsx,epsg,maxitn)
%
%   Обязательные входные аргументы функции:
%               A - точечная  mxn-матрица  интервальной  системы  линейных 
%                   алгебраических уравнений, в общем случае прямоугольна;
%               b - mx2-вектор  нижних  и  верхних  границ  правой  части 
%                   интервальной системы линейных алгебраических уравнений.
% 
%   Необязательные входные аргументы функции:
%            iprn - выдача протокола поиска;  если iprn > 0  - информация 
%                   о ходе процесса печатается через каждые iprn-итераций;
%                   если iprn <= 0 (значение по умолчанию), печати нет;
%            epsf - допуск на точность по значению целевого функционала,
%                   по умолчанию устанавливается 1.e-6;
%            epsx - допуск на точность по аргументу целевого функционала,
%                   по умолчанию устанавливается 1.e-6;
%            epsg - допуск на малость нормы суперградиента функционала,
%                   по умолчанию устанавливается 1.e-6;
%          maxitn - ограничение на количество шагов алгоритма, 
%                   по умолчанию устанавливается 2000.
%
%   Выходные аргументы функции: 
%            umax - значение  максимума  распознающего  функционала;
%            amax - доставляющий его вектор значений аргумента, который 
%                   лежит в информационном множестве задачи при umax>=0;
%            envs - вектор значений образующих распознающего функционала
%                   в точке максимума, упорядоченный по возрастанию;
%            ccod - код завершения алгоритма (1 - по допуску epsf на 
%                   изменения значений функционала, 2 - по допуску epsg 
%                   на суперградиент, 3 - по допуску epsx на вариацию 
%                   аргумента, 4 - по числу итераций, 5 - не найден 
%                   максимум по направлению).
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Эта программа выполняет исследование разрешимости интервальной линейной  
%   системы уравнений Ax = b с точечной матрицей A и интервальным вектором 
%   правой части b с помощью максимизации распознающего функционала Uni 
%   объединённого множества решений этой системы. См. подробности в 
%       Шарый С.П. Конечномерный интервальный анализ. - Новосибирск: XYZ, 
%       2010. - Электронная книга, доступная на http://www.nsc.ru/interval,
%       параграф 5.5;
% 
%   Для  максимизации  вогнутого  распознающего  функционала  используется 
%   вариант алгоритма суперградиентного подъёма с растяжением пространства 
%   в направлении разности  последовательных суперградиентов, предложенный
%   (для случая минимизации) в работе 
%       Шор Н.З., Журбенко Н.Г.  Метод минимизации, использующий  операцию
%       растяжения пространства в направлении разности двух последовательных 
%       градинетов // Кибернетика. - 1971. - №3. - С. 51-59.  
%   В качестве основы этой части программы использована процедура негладкой 
%   оптимизации ralgb5, разработанная и реализованная П.И.Стецюком (Институт 
%   кибернетики НАН Украины, Киев). 
%  
%   С.П. Шарый, 2011 год. 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%  
%   проверка корректности входных данных  
%
m = size(A,1);  n = size(A,2);  %  n - размерность пространства переменных
k = size(b,1);  l = size(b,2);
  
if ( l ~= 2 )
    error('Неправильны размеры вектора правой части')
end
  
if ( k ~= m )
    error('Размеры матрицы не соответствуют размерам правой части')
end
  
for i=1:m;
    if ( b(i,1) > b(i,2) ) 
        ErrStr = msprintf('%s %s %s','В компоненте', string(i), ... 
                                     'интервал правой части неправилен');
        error(ErrStr)
    end
end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%   задание параметров алгоритма суперградиентного подъёма 
%
maxitn = 2000;          %  ограничение на количество шагов алгоритма
nsims  = 30;            %  допустимое количество одинаковых шагов
epsf = 1.e-6;           %  допуск на изменение значения функционала 
epsx = 1.e-6;           %  допуск на изменение аргумента функционала
epsg = 1.e-6;           %  допуск на норму суперградиента функционала
  
alpha = 2.3;            %  коэффициент растяжения пространства в алгоритме
hs = 1.;                %  начальная величина шага одномерного поиска
nh = 3;                 %  число одинаковых шагов одномерного поиска 
q1 = 0.9;               %  q1, q2 - параметры адаптивной регулировки
q2 = 1.1;               %  шагового множителя
  
iprn = 0;               %  печать о ходе процесса через каждые iprn-итераций
                        %  (если iprn < 0, то печать подавляется)
format short g;         %   формат вывода данных                          
                          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   строки для оформления вывода результатов работы 
%
HorLine = '------------------------------------------------------------';
TitLine = 'Протокол максимизации распознающего функционала Uni';
TabLine = 'Шаг       Uni(x)      Uni(xx)      ФунВыч/шаг    ФунВыч';
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  переназначение параметров алгоритма, заданных пользователем 
% 
if nargin >= 3
    iprn = ceil(varargin{1});
    if nargin >= 4
        epsf = varargin{2};
        if nargin >= 5
            epsx = varargin{3};
            if nargin >= 6
                epsg = varargin{4};
                if nargin >= 7 
                    maxitn = varargin{5};
                end
            end
        end
    end 
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function [f,g] = calcfg(x)
%
%  функция, вычисляющая значения f максимизируемого распознающего 
%  функционала объединённого множества решений и его суперградиента g
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  предварительное размещение рабочих массивов 
  
tt = zeros(m,1);
dd = zeros(n,m);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   вычисляем значение распознающего функционала и матрицу dd, 
    %   составленную из суперградиентов его образующих 
    
    s = 0.5*(b(:,1)+b(:,2)) - A*x; 
    
    for i = 1:m;
        %  вычисление значения i-ой образующей и её суперградиента
        tt(i) = 0.5*(b(i,2)-b(i,1)) - abs(s(i));
        di = A(i,:);
        if s(i) < 0
            di = -di;
        end
        dd(:,i) = di';
    end  
  
    %  выбираем минимальную по значению образующую
    %  и конструируем общий суперградиент
    [f,mc] = min(tt);
    g = dd(:,mc);
  
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  формируем начальное приближение x как решение либо псевдорешение 
%  'средней' точечной системы, если она не слишком плохо обусловлена,
%  иначе берём начальным приближением нулевой вектор 
% 
sv = svd(A);
minsv = min(sv);
maxsv = max(sv);
  
if ( minsv~=0 && maxsv/minsv < 1.e+23 ) 
    x = A\(0.5*(b(:,1)+b(:,2))); 
else
    x = zeros(n,1);
end  
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Рабочие массивы:
%      B - матрица обратного преобразования пространства аргументов
%      g - вектор суперградиента максимизируемого функционала 
%      g1,g2 - используются для хранения вспомогательных векторов
%      vf - вектор приращений функционала на последних nsims шагах алгоритма
  
B  = eye(n,n);
vf = realmax*ones(nsims,1);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%  установка начальных параметров, вывод заголовка протокола
%
w = 1./alpha - 1.;
lp = iprn;
  
[f,g0] = calcfg(x);
ff = f;   xx = x;
cal = 1;  ncals = 1;
  
if iprn > 0 
    disp(' '); 
    disp([ blanks(11), TitLine ]);
    disp([ blanks(7), HorLine ]);
    disp([ blanks(10), TabLine ]);
    disp([ blanks(7), HorLine ]);
    disp([ 0, f, ff, cal, ncals ]);
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  основной цикл программы:
%      itn - счётчик числа итераций
%      xx  - найденная точка максимума функционала
%      ff  - значение функционала в точке максимума
%      cal - количество вычислений функционала на текущем шаге 
%    ncals - общее количество вычислений целевого функционала
% 
for itn = 1:maxitn;
    pf = ff;
    %  критерий останова по норме суперградиента
    if  norm(g0) < epsg
        ccod = 2;  
        break
    end
    %  вычисляем суперградиент в преобразованном пространстве
    g1 = B' * g0;
    g = B * g1/norm(g1); 
    normg = norm(g);
    %  одномерный подъём по направлению dx:
    %      cal - счётчик шагов одномерного поиска,
    %      deltax - вариация аргумента в процессе поиска
    r = 1;  
    cal = 0;  
    deltax = 0;
    while ( r > 0. && cal <= 500 )
        cal = cal + 1; 
        x = x + hs*g;  
        deltax = deltax + hs*normg;
        [f, g1] = calcfg(x);
        if f > ff  
            ff = f;
            xx = x; 
        end 
        if cal > nh 
            hs = hs*q2;
        end  
        r = g'*g1; 
    end
    %  если превышено число шагов одномерного подъёма, то выход
    if cal == 500 
        ccod = 5;  
        break; 
    end 
    %  если одномерный подъём занял один шаг, уменьшаем величину шага
    if cal == 1
        hs = hs*q1;
    end 
    %  уточняем статистику и при необходимости выводим её
    ncals = ncals + cal;
    if itn==lp
        disp([ itn, f, ff, cal, ncals ]);        
        lp = lp + iprn;
    end
    %  если вариация аргумента в одномерном поиске мала, то выход
    if deltax < epsx 
        ccod = 3;  
        break; 
    end
    %   пересчитываем матрицу преобразования пространства 
    dg = B' * (g1 - g0);
    xi = dg / norm(dg);
    B = B + w*B*xi*xi';
    g0 = g1;
    %   проверка изменения значения функционала на последних nsims шагах,
    %   относительного либо абсолютного 
    vf(2:nsims) = vf(1:nsims-1);
    vf(1) = abs(ff - pf);
    if abs(ff) > 1 
        deltaf = sum(vf)/abs(ff);
    else
        deltaf = sum(vf);  
    end
    if deltaf < epsf 
        ccod = 1;  break
    end
    ccod = 4;
end
   
umax = ff;
amax = xx;
%   переупорядочение образующих по возрастанию 
tt = [(1:m)', tt];
[z,ind] = sort(tt(:,2));
envs = tt(ind,:);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  вывод результатов работы 
  
if iprn > 0
    if rem(itn,iprn)~=0
        disp([ itn, f, ff, cal, ncals ]);
    end
    disp([ blanks(7), HorLine ]);
end
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
disp(' ');
if umax >= 0
    disp(' Информационное множество задачи непусто');
else 
    disp(' Информационное множество задачи пусто');
end 
disp(' ');
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
if ( umax < 0. && abs(umax/epsf) < 10 ) 
    disp(' Вычисленный максимум находится в пределах заданной точности');
    disp(' - можно перезапустить программу с меньшими epsf и epsx, чтобы');
    disp(' получить больше информации о пустоте/непустоте информационного');
    disp(' множества'); 
    disp(' ');
end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
