function [amax,umax,envs,ccod] = lintregs(A,b,varargin) 
//  Построение линейной интервальной регрессии. 
//  LINTREGS(A,b) выдаёт  значения  коэффициентов  линейной  регрессионной
//  зависимости по матрице  входных значений  A  и  интервальному  вектору 
//  выходных данных b путём вычисления максимума распознающего функционала
//  множества решений интервальной линейной системы Ax = b. 
//  
//  Синтаксис вызова:
//      [amax,umax,envs,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 = 1.0;               //  q1, q2 - параметры адаптивной регулировки
q2 = 1.1;               //  шагового множителя
  
iprn = 0;               //  печать о ходе процесса через каждые iprn-итераций
                        //  (если iprn < 0, то печать подавляется)
    
////////////////////////////////////////////////////////////////////////////////
// 
//  переназначение параметров алгоритма, заданных пользователем 
// 
nargin = argn(2); 
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,tt] = calcfg(x)
// 
//  функция, вычисляющая значения f максимизируемого распознающего 
//  функционала объединённого множества решений и его суперградиента g;
//  кроме того, выводится вектор tt значений образующих функционала
//
////////////////////////////////////////////////////////////////////////////////
//
//  предварительное размещение рабочих массивов 
  
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  
  
    //  выбираем минимальную образующую функционала
    //  и присваиваем общий суперградиент g
    [f,mc] = min(tt);
    g = dd(:,mc);
  
endfunction 
  
////////////////////////////////////////////////////////////////////////////////
// 
//  формируем начальное приближение 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 = %inf*ones(nsims,1);

////////////////////////////////////////////////////////////////////////////////
//
//  Присваивание строк оформления протокола работы программы 

TitLine = 'Протокол максимизации распознающего функционала Uni';
HorLine = '-----------------------------------------------------------';
TabLine = ' Шаг          Uni(x)          Uni(xx)  ВычФун/шаг  ВычФун';
  
////////////////////////////////////////////////////////////////////////////////
//  
//  установка начальных параметров, вывод заголовка протокола
//
w = 1./alpha - 1.;
lp = iprn;
  
[f,g0,tt] = calcfg(x);
ff = f;   xx = x;
cal = 1;  ncals = 1;
  
if iprn > 0 
    printf('\n%100s',TitLine);
    printf('\n%60s',HorLine);
    printf('\n%77s',TabLine);
    printf('\n%60s',HorLine);
    printf('\n%5u%17g%17g%9u%9u',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,tt] = 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
        printf('\n%5u%17g%17g%9u%9u',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] = gsort(tt(:,2),'g','i');
envs = tt(ind,:);
   
////////////////////////////////////////////////////////////////////////////////
//
//  вывод результатов работы 
  
if iprn > 0
    if modulo(itn,iprn)~=0
        printf('\n%5u%17g%17g%9u%9u',itn,f,ff,cal,ncals);
    end
    printf('\n%60s\n',HorLine);
end
  
////////////////////////////////////////////////////////////////////////////////
  
if umax >= 0
    printf('\n%48s\n\n','Информационное множество непусто')
else 
    printf('\n%54s\n\n','Информационное множество пусто')
end 
  
////////////////////////////////////////////////////////////////////////////////
  
if ( umax < 0. & abs(umax/epsf) < 10 ) 
    printf('%55s\n','Найденное значение максимума сравнимо с заданной точностью')
    printf('%50s\n','- можно перезапустить программу с меньшими epsf и epsx,')
    printf('%55s\n\n','чтобы получить больше информации о совместности системы')
end 
  
////////////////////////////////////////////////////////////////////////////////
  
endfunction
  
////////////////////////////////////////////////////////////////////////////////