This file is indexed.

/usr/share/octave/packages/interval-3.1.0/verlinprog.m is in octave-interval 3.1.0-5.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
## Copyright 2007 Jiří Rohn
## Copyright 2016 Oliver Heimlich
##
## This program is derived from verlinprog in VERSOFT, published on
## 2016-07-26, which is distributed under the terms of the Expat license,
## a.k.a. the MIT license.  Original Author is Jiří Rohn.  Migration to Octave
## code has been performed by Oliver Heimlich.
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @documentencoding UTF-8
## @deftypefun {[@var{flag}, @var{x}, @var{y}, @var{h}] =} verlinprog (@var{A}, @var{b}, @var{c})
## Verified linear programming.
##
## For a real matrix @var{A} (full or sparse) and matching real vectors
## @var{b}, @var{c}, this function either computes verified optimal solution
## @var{x}, verified dual optimal solution @var{y} and verified optimal value
## @var{h} of the linear programming problem
## @display
## min @var{c}' * @var{x} subject to @var{A} * @var{x} = @var{b}, x @geq{} 0,
## @end display
## or verifies (in)feasibility, or verifies unboundedness, or yields no
## verified result. The respective outcome is always described verbally in the
## variable @var{flag}.
##
## Possible values of @var{flag}:
## @table @option
## @item verified optimum
## @var{x} is verified to enclose a primal optimal solution,
## @var{y} is verified to enclose a dual optimal solution,
## @var{h} is verified to enclose the optimal value,
##
## @item verified unbounded
## @var{x} is verified to enclose a primal feasible solution @var{xo}, and
## @var{y} is verified to enclose a vector @var{yo} such that the objective
## tends to -Inf along the feasible half-line
## @{@var{xo} + @var{t} * @var{yo} | @var{t} @geq{} 0@},
## @var{h} is empty,
##
## @item verified feasible
## @var{x} is verified to enclose a primal feasible solution
## (neither optimality nor unboundedness could be verified),
## @var{y}, @var{h} are empty,
##
## @item verified infeasible
## @var{y} is verified to enclose a Farkas vector @var{yo} satisfying
## @var{A}' * @var{yo} @geq{} 0, @var{b}' * @var{yo} < 0
## (whose existence proves primal infeasibility),
## @var{x}, @var{h} are empty,
##
## @item no verified result
## @var{x}, @var{y}, and @var{h} are empty (no verified result could be found).
## @end table
##
## Complexity: The algorithm solves at most four linear programming problems
## (independently of the size of the original problem) and uses a verification
## procedure which runs approximately in O(@var{m}³) time, where
## @var{m} = rows (@var{A}).
##
## This work was supported by the Czech Republic National Research
## Program “Information Society”, project 1ET400300415.
## @seealso{linprog}
## @end deftypefun

## Author: Jiří Rohn
## Keywords: interval
## Created: 2007

function [flag, x, y, h] = verlinprog (A, b, c)

  if (nargin ~= 3)
    print_usage ();
    return
  endif

  b = b(:); c = c(:);
  [m, n] = size(A);
  p = length (b); q = length (c);

  flag = "no verified result";
  x = repmat (infsup, n, 1);
  y = repmat (infsup, m, 1);
  h = infsup;

  if (~(m == p && n == q) || (m > n))
    error ("verlinprog: sizes do not match");
  endif
  if (~isreal (A) || ~isreal (b) || ~isreal (c))
    error("verlinprog: data not real");
  endif
  if issparse (b)
    b = full (b);
  end
  if issparse (c)
    c = full (c);
  end

                                # verifying infeasibility
  yi = verinfeas (A, b);
  if (~isempty(yi(1))) # verified Farkas vector found
    y = yi;
    flag = "verified infeasible";
    return
  endif

                                # verifying feasibility
  xf = veropt (A, b, ones (n, 1));
  if (isempty (xf(1))) # verified feasible solution not found
    flag = "no verified result";
    return
  endif

                                # verifying unboundedness
  yu = verunbound (A, c);
  if (~isempty (yu(1))) # verified descent direction found
    x = xf;
    y = yu;
    flag = "verified unbounded";
    return
  endif

                                # verifying optimality
  [xo, B, N] = veropt (A, b, c);
  if (isempty (xo(1))) % verified feasible primal solution with basis B not found
    x = xf; # previous feasible solution outputted
    flag = "verified feasible";
    return
  endif

  AB = A(:, B);
  if (issparse (AB))
    AB = full (AB); # only the square submatrix taken full
  endif
  yB = mldivide (infsup (AB'), infsup (c(B)));
  if (isempty (yB(1))) # verified feasible dual solution not found
    x = xo; # candidate for optimum outputted as feasible solution
    flag = "verified feasible";
    return
  endif

  c = infsup (c);
  A = infsup (A);
  crit = c' - yB' * A; # criterial row (dual feasibility)
  crit = crit(N);      # nonbasic part of it
  if (~all (crit.inf >= 0)) # verified feasible dual solution not found
    x = xo; % candidate for optimum outputted as feasible solution
    flag = "verified feasible";
    return
  endif

# verified quantities     # verified primal and dual feasible solutions found
  x = xo;                   # x is a verified primal optimal solution
  y = yB;                   # y is a verified dual optimal solution
  if (nargout >= 3)
    h1 = c' * x;
    h2 = b' * y;
    h = intersect (h1, h2);   # h is a verified optimal value (duality theorem)
    if (isempty (h))
      h = h1;
    end
  endif
  flag = "verified optimum";
endfunction


function [x, B, N] = veropt (A, b, c)
  ## B is the "basis index set" of an optimal solution of the LP problem
  ## min c'*x  subject to  A*x=b, x>=0,
  ## x is a verified basic feasible solution with this basis
  ## N is the set of nonbasic indices

  persistent GLP_MSG_OFF = 0;
  [m, n] = size (A);
  x = repmat(infsup, n, 1);
  B = nan (m, 1);
  N = nan (n, 1);

  [xopt, ~, exitflag] = glpk (c, A, b, ...
                              [], [], ... # 0 <= x <= inf
                              repmat ("S", 1, m), ... # equality constraints Ax = b
                              repmat ("C", 1, n), ... # continuous variable x
                              1, ... # minimization
                              struct ("msglev", GLP_MSG_OFF));

  if (exitflag ~= 0)
    return
  endif

  [xx, J] = sort (xopt);
  B = J(n - m + 1 : n); # B is set of "basic" indices,
  N = J(1 : n - m);     # N of "nonbasic" ones

  AB = A(:, B);
  if (issparse (AB))
    AB = full (AB); # only the square submatrix taken full (because of mldivide)
  endif
  xB = mldivide (infsup (AB), infsup (b));
  if (isempty (xB(1)) || ~all (xB.inf >= 0))
                               # verified "optimal" solution not found
    return
  endif

                                # verified "optimal" solution found
  x = infsup (zeros (n, 1));
  x(B) = xB;
endfunction


function y = verinfeas (A, b)
  # y verified to enclose a Farkas vector yo (i.e., satisfying A'*yo>=0, b'*yo<0)
  # its existence implies infeasibility of A*x=b

  [m, n] = size (A);
  y = repmat(infsup, m, 1);
  ep = max (1e-08, max ([m n 100]) * max ([norm(A, inf) norm(b, inf)]) * eps);
  Afv = [A' -A'    -eye(n) zeros(n,1);    # Afv is (n+1)x(2*m+n+1)
         b' -b' zeros(1,n)         1];
  bfv = [zeros(n, 1)' -1]';               # bfv is (n+1)x1

        # perturbation to compensate roundoff errors (so that A'*y>=0)
  bfv = bfv + ep * [ones(1, n) -1]';
  yf = veropt (Afv, bfv, ones (2 * m + n + 1, 1)); % system: A'*y>=0, b'*y<=-1, y written as y=y1-y2
  if (~isempty (yf(1)))
    yf = mid (yf);
    y1 = yf(1 : m);
    y2 = yf(m + 1 : 2 * m);
    yf = y1 - y2;                      # would-be Farkas vector
    A = infsup (A);
    b = infsup (b);
    yf = infsup (yf);                  # (i.e., should satisfy A'*y>=0, b'*y<0)
    alpha = A' * yf;
    beta = b' * yf;
    if (all (alpha.inf >= 0)) && (beta.sup < 0)
                                # infeasibility verified
      y=yf; # Farkas vector outputted
    endif
  endif
endfunction


function y = verunbound (A, c)
  # y verified to enclose a vector yo satisfying A*yo=0, yo>=0, c'*yo<=-1
  # under feasibility its existence implies unboundedness

  [m, n] = size (A);
  y = repmat (infsup, n, 1);
  Aunb = [A zeros(m, 1);                       # Aunb is (m+1)x(n+1)
          c'         1];
  bunb = [zeros(1, m) -1]';                    # bunb is (m+1)x1
  yunb = veropt (Aunb, bunb, ones (n + 1, 1)); # yunb is (n+1)x1
  if (~isempty (yunb(1)))
                                # y satisfies A*y=0, y>=0, c'*y=-1
    y = yunb(1:n);
    return
  endif
endfunction

%!test
%! A = [-2, -3; -2, -1];
%! b = [-1500, -1000];
%! c = [1; 1];
%! [flag, x, y, h] = verlinprog (A, b, c);
%! assert (flag, "verified optimum");
%! assert (ismember ([375; 250], x));
%! assert (wid (x) < 1e-12);
%! assert (ismember ([-0.25; -0.25], y));
%! assert (wid (y) < 1e-16);
%! assert (ismember (625, h));
%! assert (wid (h) < 1e-12);