# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.1 `mip` Chapter.
``mip`` - Operations Research
This module provides functions to solve certain integer programming, transportation and shortest path problems.
Additionally 'best subset' functions are included.
See Also
--------
``naginterfaces.library.examples.mip`` :
This subpackage contains examples for the ``mip`` module.
See also the :ref:`library_mip_ex` subsection.
Functionality Index
-------------------
**Mixed integer linear programming (MILP)**
dense
branch and bound method: :meth:`ilp_dense`
large-scale
branch and bound method: :meth:`handle_solve_milp`
**Mixed integer quadratic programming (MIQP)**
dense
branch and bound method: :meth:`iqp_dense`
sparse
branch and bound method: :meth:`iqp_sparse`
**Mixed integer nonlinear programming (MINLP)**
dense
mixed integer sequential quadratic programming (MISQP): :meth:`sqp`
**Operations Research (OR)**
feature selection
best subset of given size
direct communication: :meth:`best_subset_given_size`
reverse communication: :meth:`best_subset_given_size_revcomm`
shortest path through directed or undirected network: :meth:`shortestpath`
transportation problem: :meth:`transportation`
travelling salesman problem, simulated annealing: :meth:`tsp_simann`
**Service functions**
input and output (I/O)
print solution of a dense MILP problem: :meth:`ilp_print`
read MILP problem from MPS file and solve it by branch and bound: :meth:`ilp_mpsx`
read MPS file defining dense MILP problem: :meth:`ilp_mpsx_convert`
option setting functions
:meth:`iqp_dense`
supply option values from a character string: :meth:`iqp_dense_optstr`
supply option values from external file: :meth:`iqp_dense_optfile`
:meth:`iqp_sparse`
supply option values from a character string: :meth:`iqp_sparse_optstr`
supply option values from external file: :meth:`iqp_sparse_optfile`
:meth:`sqp`
supply option values from a character string: :meth:`optset`
retrieve option values: :meth:`optget`
miscellaneous
extract further information on the solution obtained from :meth:`ilp_dense`: :meth:`ilp_info`
For full information please refer to the NAG Library document
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/hintro.html
"""
# NAG Copyright 2017-2024.
[docs]def ilp_dense(itmax, msglvl, a, bl, bu, intvar, cvec, maxnod, intfst, toliv, tolfes, bigbnd, x, comm, maxdpt=None, io_manager=None):
r"""
``ilp_dense`` solves 'zero-one', 'general', 'mixed' or 'all' integer programming problems using a branch and bound method.
The function may also be used to find either the first integer solution or the optimum integer solution.
It is not intended for large sparse problems.
.. deprecated:: 29.3.0.0
``ilp_dense`` is deprecated.
Please use :meth:`handle_solve_milp` instead.
See also the :ref:`Replacement Calls <replace>` document.
.. _h02bb-py2-py-doc:
For full information please refer to the NAG Library document for h02bb
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bbf.html
.. _h02bb-py2-py-parameters:
**Parameters**
**itmax** : int
An upper bound on the number of iterations for each LP problem.
**msglvl** : int
The amount of printout produced by ``ilp_dense``, as indicated below (see `Description of Printed Output <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bbf.html#parameters1>`__ for a description of the printed output). All output is written to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`).
.. rst-class:: nag-rules-none nag-align-left
+-----+---------------------------------------------------------------------------------------------------------------+
|Value|Definition |
+=====+===============================================================================================================+
|0 |No output. |
+-----+---------------------------------------------------------------------------------------------------------------+
|1 |The final IP solution only. |
+-----+---------------------------------------------------------------------------------------------------------------+
|5 |One line of output for each node investigated and the final IP solution. |
+-----+---------------------------------------------------------------------------------------------------------------+
|10 |The original LP solution (first node), one line of output for each node investigated and the final IP solution.|
+-----+---------------------------------------------------------------------------------------------------------------+
**a** : float, array-like, shape :math:`\left(m, :\right)`
Note: the required extent for this argument in dimension 2 is determined as follows: if :math:`m > 0`: :math:`n`; if :math:`m=0`: :math:`1`; otherwise: :math:`0`.
The :math:`\textit{i}`\ th row of :math:`\mathrm{a}` must contain the coefficients of the :math:`\textit{i}`\ th general constraint, for :math:`\textit{i} = 1,2,\ldots,m`.
If :math:`m = 0` then the array :math:`\mathrm{a}` is not referenced.
**bl** : float, array-like, shape :math:`\left(n+m\right)`
:math:`\mathrm{bl}` must contain the lower bounds and :math:`\mathrm{bu}` the upper bounds, for all the constraints in the following order. The first :math:`n` elements of each array must contain the bounds on the variables, and the next :math:`m` elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., :math:`l_j = {-\infty }`), set :math:`\mathrm{bl}[j]\leq -\mathrm{bigbnd}` and to specify a nonexistent upper bound (i.e., :math:`u_j = {+\infty }`), set :math:`\mathrm{bu}[j]\geq \mathrm{bigbnd}`. To specify the :math:`j`\ th constraint as an equality, set :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = \beta`, say, where :math:`\left\lvert \beta \right\rvert < \mathrm{bigbnd}`.
**bu** : float, array-like, shape :math:`\left(n+m\right)`
:math:`\mathrm{bl}` must contain the lower bounds and :math:`\mathrm{bu}` the upper bounds, for all the constraints in the following order. The first :math:`n` elements of each array must contain the bounds on the variables, and the next :math:`m` elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., :math:`l_j = {-\infty }`), set :math:`\mathrm{bl}[j]\leq -\mathrm{bigbnd}` and to specify a nonexistent upper bound (i.e., :math:`u_j = {+\infty }`), set :math:`\mathrm{bu}[j]\geq \mathrm{bigbnd}`. To specify the :math:`j`\ th constraint as an equality, set :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = \beta`, say, where :math:`\left\lvert \beta \right\rvert < \mathrm{bigbnd}`.
**intvar** : int, array-like, shape :math:`\left(n\right)`
Indicates which are the integer variables in the problem. For example, if :math:`x_{\textit{j}}` is an integer variable then :math:`\mathrm{intvar}[\textit{j}]` must be set to :math:`1`, and :math:`0` otherwise.
**cvec** : float, array-like, shape :math:`\left(n\right)`
The coefficients :math:`c_j` of the objective function :math:`F\left(x\right) = c_1x_1+c_2x_2 + \cdots +c_nx_n`. The function attempts to find a minimum of :math:`F\left(x\right)`. If a maximum of :math:`F\left(x\right)` is desired, :math:`\mathrm{cvec}[\textit{j}]` should be set to :math:`{-c_{\textit{j}}}`, for :math:`\textit{j} = 1,2,\ldots,n`, so that the function will find a minimum of :math:`{-F}\left(x\right)`.
**maxnod** : int
The maximum number of nodes that are to be searched in order to find a solution (optimum integer solution). If :math:`\mathrm{maxnod}\leq 0` and :math:`\mathrm{intfst}\leq 0`, then the B&B tree search is continued until all the nodes have been investigated.
**intfst** : int
Specifies whether to terminate the B&B tree search after the first integer solution (if any) is obtained. If :math:`\mathrm{intfst} > 0` then the B&B tree search is terminated at node :math:`k` say, which contains the first integer solution. For :math:`\mathrm{maxnod} > 0` this applies only if :math:`k\leq \mathrm{maxnod}`.
**toliv** : float
The integer feasibility tolerance; i.e., an integer variable is considered to take an integer value if its violation does not exceed :math:`\mathrm{toliv}`. For example, if the integer variable :math:`x_j` is near unity then :math:`x_j` is considered to be integer only if :math:`\left(1-\mathrm{toliv}\right)\leq x_j\leq \left(1+\mathrm{toliv}\right)`.
**tolfes** : float
The maximum acceptable absolute violation in each constraint at a 'feasible' point (feasibility tolerance); i.e., a constraint is considered satisfied if its violation does not exceed :math:`\mathrm{tolfes}`.
**bigbnd** : float
The 'infinite' bound size in the definition of the problem constraints. More precisely, any upper bound greater than or equal to :math:`\mathrm{bigbnd}` will be regarded as :math:`{+\infty }` and any lower bound less than or equal to :math:`{-\mathrm{bigbnd}}` will be regarded as :math:`{-\infty }`.
**x** : float, array-like, shape :math:`\left(n\right)`
An initial estimate of the original LP solution.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.lp_solve'}`.
**maxdpt** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:`3\times n/2`.
The maximum depth of the B&B tree used for branch and bound.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**itmax** : int
Unchanged if on entry :math:`\mathrm{itmax} > 0`, else :math:`\mathrm{itmax} = \mathrm{max}\left(50, {5\times \left(n+m\right)}\right)`.
**toliv** : float
Unchanged if on entry :math:`\mathrm{toliv} > 0.0`, else :math:`\mathrm{toliv} = 10^{-5}`.
**tolfes** : float
Unchanged if on entry :math:`\mathrm{tolfes} > 0.0`, else :math:`\mathrm{tolfes} = \sqrt{\epsilon }` (where :math:`\epsilon` is the machine precision).
**bigbnd** : float
Unchanged if on entry :math:`\mathrm{bigbnd} > 0.0`, else :math:`\mathrm{bigbnd} = 10^{20}`.
**x** : float, ndarray, shape :math:`\left(n\right)`
With no exception or warning is raised, :math:`\mathrm{x}` contains a solution which will be an estimate of either the optimum integer solution or the first integer solution, depending on the value of :math:`\mathrm{intfst}`. If :math:`\mathrm{errno}` = 9, then :math:`\mathrm{x}` contains a solution which will be an estimate of the best integer solution that was obtained by searching :math:`\mathrm{maxnod}` nodes.
**objmip** : float
With the function exits successfully or :math:`\mathrm{errno}` = 9, :math:`\mathrm{objmip}` contains the value of the objective function for the IP solution.
.. _h02bb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
The problem does not have a feasible integer solution.
(`errno` :math:`2`)
The :math:`\mathrm{LP}` solution is unbounded.
(`errno` :math:`3`)
The :math:`\mathrm{LP}` does not have a feasible solution.
(`errno` :math:`4`)
Iteration limit reached without finding a solution.
(`errno` :math:`6`)
On entry, there were :math:`\langle\mathit{\boldsymbol{value}}\rangle` infinite or inconsistent bounds given in arrays :math:`\mathrm{bl}` and :math:`\mathrm{bu}`. For further advisory details set :math:`\mathrm{msglvl} > 2`.
(`errno` :math:`6`)
On entry, the dimension of :math:`\mathrm{comm}`\ ['iwork'] is too small: :math:`\textit{liwork} = \langle\mathit{\boldsymbol{value}}\rangle`. :math:`\mathrm{comm}`\ ['iwork'] must be of dimension (at least) :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`6`)
On entry, the dimension of :math:`\mathrm{comm}`\ ['rwork'] is too small: :math:`\textit{lrwork} = \langle\mathit{\boldsymbol{value}}\rangle`. :math:`\mathrm{comm}`\ ['rwork'] must be of dimension (at least) :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{maxdpt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxdpt}\geq 2`.
(`errno` :math:`6`)
On entry, :math:`i = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{intvar}[i-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{intvar} = 0` or :math:`1`.
(`errno` :math:`6`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 0`.
(`errno` :math:`6`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
(`errno` :math:`8`)
Increase :math:`\mathrm{maxdpt}` and rerun ``ilp_dense``.
(`errno` :math:`8`)
:math:`\mathrm{maxdpt}` is too small to solve the problem: :math:`\mathrm{maxdpt} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`10`)
No feasible solution was found for the number of nodes investigated in the B&B tree.
(`errno` :math:`11`)
Not enough workspace to solve problem.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`7`)
Search of a branch was terminated due to iteration limit.
(`errno` :math:`9`)
The :math:`\mathrm{IP}` solution returned is the best solution for the number of nodes investigated in the B&B tree.
.. _h02bb-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``ilp_dense`` is capable of solving certain types of integer programming (IP) problems using a branch and bound (B&B) method, see Taha (1987).
In order to describe these types of integer programs and to briefly state the B&B method, we define the following Linear Programming (LP) problem:
Minimize
.. math::
F\left(x\right) = c_1x_1+c_2x_2 + \cdots +c_nx_n
subject to
.. math::
\sum_{{j = 1}}^na_{{ij}}x_j\left\{\begin{array}{c} = \\\leq \\\geq \end{array}\right\}b_i\text{, }\quad i = 1,2,\ldots,m
.. math::
l_j\leq x_j\leq u_j\text{, }\quad j = 1,2,\ldots,n
If, in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bbf.html#eqn1>`__, it is required that (some or) all the variables take integer values, then the integer program is of type (`mixed` or) `all` general IP problem.
If additionally, the integer variables are restricted to take only :math:`0`--:math:`1` values (i.e., :math:`l_j = 0` and :math:`u_j = 1`) then the integer program is of type (mixed or all) `zero-one` IP problem.
The B&B method applies directly to these integer programs.
The general idea of B&B (for a full description see Dakin (1965) and Mitra (1973)) is to solve the problem without the integral restrictions as an LP problem (first `node`).
If in the optimal solution an integer variable :math:`x_k` takes a noninteger value :math:`x_k^*`, two LP sub-problems are created by `branching`, imposing :math:`x_k\leq \left[x_k^*\right]` and :math:`x_k\geq \left[x_k^*\right]+1` respectively, where :math:`\left[x_k^*\right]` denotes the integer part of :math:`x_k^*`.
This method of branching continues until the first integer solution (`bound`) is obtained.
The hanging nodes are then solved and investigated in order to prove the optimality of the solution.
At each node, an LP problem is solved using :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>`.
.. _h02bb-py2-py-references:
**References**
Dakin, R J, 1965, `A tree search algorithm for mixed integer programming problems`, Comput. J. (8), 250--255
Mitra, G, 1973, `Investigation of some branch and bound strategies for the solution of mixed integer linear programs`, Math. Programming (4), 155--170
Taha, H A, 1987, `Operations Research: An Introduction`, Macmillan, New York
See Also
--------
:meth:`naginterfaces.library.examples.mip.ilp_dense_ex.main`
"""
raise NotImplementedError
[docs]def ilp_mpsx(infile, maxn, maxm, optim, xbldef, xbudef, maxdpt, msglvl, liwork, lrwork, comm, io_manager=None):
r"""
``ilp_mpsx`` solves linear or integer programming problems specified in MPSX input format.
It is not intended for large sparse problems.
.. _h02bf-py2-py-doc:
For full information please refer to the NAG Library document for h02bf
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bff.html
.. _h02bf-py2-py-parameters:
**Parameters**
**infile** : int
The unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) associated with the MPSX data file.
**maxn** : int
An upper limit for the number of variables in the problem.
**maxm** : int
An upper limit for the number of constraints (including the objective) in the problem.
**optim** : str, length 3
Specifies the direction of the optimization. :math:`\mathrm{optim}` must be set to 'MIN' for minimization and to 'MAX' for maximization.
**xbldef** : float
The default lower bound to be used for the variables in the problem when none is specified in the BOUNDS section of the MPSX data file. For a standard LP or IP problem :math:`\mathrm{xbldef}` would normally be set to zero.
**xbudef** : float
The default upper bound to be used for the variables in the problem when none is specified in the BOUNDS section of the MPSX data file. For a standard LP or IP problem :math:`\mathrm{xbudef}` would normally be set to 'infinity' (i.e., :math:`\mathrm{xbudef}\geq 10^{20}`).
**maxdpt** : int
For an IP problem, :math:`\mathrm{maxdpt}` must specify the maximum depth of the branch and bound tree.
**msglvl** : int
The amount of printout produced by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_dense`, as indicated below. For a description of the printed output see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/e04/e04mff.html#fc-printedoutput>`__ or `Description of Printed Output for ilp_dense <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bbf.html#parameters1>`__ (as appropriate). All output is written to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`).
For an LP problem (:meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>`):
.. rst-class:: nag-rules-none nag-align-left
+----------+--------------------------------------------------------------------------+
|Value |Definition |
+==========+==========================================================================+
|:math:`0` |No output. |
+----------+--------------------------------------------------------------------------+
|:math:`1` |The final solution only. |
+----------+--------------------------------------------------------------------------+
|:math:`5` |One line of output for each iteration (no printout of the final solution).|
+----------+--------------------------------------------------------------------------+
|:math:`10`|The final solution and one line of output for each iteration. |
+----------+--------------------------------------------------------------------------+
For an IP problem (:meth:`ilp_dense`):
.. rst-class:: nag-rules-none nag-align-left
+----------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|Value |Definition |
+==========+==================================================================================================================================================================================================+
|:math:`0` |No output. |
+----------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |The final IP solution only. |
+----------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`5` |One line of output for each node investigated and the final IP solution. |
+----------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`10`|The original LP solution (first node) with dummy names for the rows and columns, one line of output for each node investigated and the final IP solution with MPSX names for the rows and columns.|
+----------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
**liwork** : int
The dimension of the array :math:`\mathrm{iwork}`.
**lrwork** : int
The dimension of the array :math:`\mathrm{rwork}`.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.lp_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**n** : int
:math:`n`, the actual number of variables in the problem.
**m** : int
:math:`m`, the actual number of general linear constraints in the problem.
**x** : float, ndarray, shape :math:`\left(\mathrm{maxn}\right)`
The solution to the problem, stored in :math:`\mathrm{x}[0],\mathrm{x}[1],\ldots,\mathrm{x}[\mathrm{n}-1]`. :math:`\mathrm{x}[\textit{i}-1]` is the value of the variable whose MPSX name is stored in :math:`\mathrm{crname}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}`.
**crname** : str, length 8, ndarray, shape :math:`\left(\mathrm{maxn}+\mathrm{maxm}\right)`
The first :math:`\mathrm{n}` elements contain the MPSX names for the variables in the problem.
**iwork** : int, ndarray, shape :math:`\left(\mathrm{liwork}\right)`
The first (:math:`\mathrm{n}+\mathrm{m}`) elements contain ISTATE (the status of the constraints in the working set at the solution). Further details can be found in :ref:`Parameters for opt.lp_solve <e04mf-py2-py-parameters>` and :ref:`Parameters for ilp_info <h02bz-py2-py-parameters>` (as appropriate).
**rwork** : float, ndarray, shape :math:`\left(\mathrm{lrwork}\right)`
The first (:math:`\mathrm{n}+\mathrm{m}`) elements contain BL (the lower bounds), the next (:math:`\mathrm{n}+\mathrm{m}`) elements contain BU (the upper bounds) and the next (:math:`\mathrm{n}+\mathrm{m}`) elements contain CLAMDA (the Lagrange-multipliers). Further details can be found in :ref:`Parameters for opt.lp_solve <e04mf-py2-py-parameters>` and :ref:`Parameters for ilp_info <h02bz-py2-py-parameters>` (as appropriate). Note that for an IP problem the contents of BL and BU may not be the same as those originally specified in the MPSX data file and/or via the arguments :math:`\mathrm{xbldef}` and :math:`\mathrm{xbudef}`.
.. _h02bf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`i < 0`)
Either :math:`\mathrm{maxm}` and/or :math:`\mathrm{maxn}` are too small or the MPSX data file is nonstandard and/or corrupt.
(`errno` :math:`2`)
The LP solution is unbounded.
(`errno` :math:`3`)
The LP does not have a feasible solution.
(`errno` :math:`4`)
Iteration limit reached without finding a solution.
(`errno` :math:`5`)
On entry, not enough real workspace to solve problem: :math:`\mathrm{lrwork} = \langle\mathit{\boldsymbol{value}}\rangle` :math:`\mathrm{lrwork}` must be at least :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, not enough integer workspace to solve problem: :math:`\mathrm{liwork} = \langle\mathit{\boldsymbol{value}}\rangle` :math:`\mathrm{liwork}` must be at least :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{maxdpt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxdpt}\geq 2`.
(`errno` :math:`5`)
On entry, not enough real workspace to read data file: :math:`\mathrm{lrwork} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, not enough integer workspace to read data file: :math:`\mathrm{liwork} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{optim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{optim} = \texttt{'MIN'}` or :math:`\texttt{'MAX'}`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{xbldef} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xbudef} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xbldef}\leq \mathrm{xbudef}`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{maxm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxm}\geq 1`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{maxn} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxn}\geq 1`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{infile} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{infile}\leq 2147483647`.
(`errno` :math:`6`)
A serious error has occurred. Check all function calls and array dimensions.
(`errno` :math:`8`)
Increase :math:`\mathrm{maxdpt}` and rerun ``ilp_mpsx``.
(`errno` :math:`8`)
:math:`\mathrm{maxdpt}` is too small to solve the problem: :math:`\mathrm{maxdpt} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`9`)
The IP solution returned is the best solution for the number of nodes investigated in the branch and bound tree.
(`errno` :math:`10`)
No feasible solution was found for the number of nodes investigated in the branch and bound tree.
(`errno` :math:`11`)
Not enough workspace to solve problem.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1`)
Weak LP solution.
(`errno` :math:`1`)
The problem does not have a feasible integer solution.
(`errno` :math:`7`)
Search of a branch was terminated due to iteration limit.
.. _h02bf-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``ilp_mpsx`` solves Linear Programming (LP) or integer programming (IP) problems specified in MPSX (see IBM (1971)) input format.
It calls either :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` (to solve an LP problem) or :meth:`ilp_dense` and :meth:`ilp_info` (to solve an IP problem); these functions are designed to solve problems of the form
.. math::
\mathrm{minimize}_{{x \in R^n}}c^\mathrm{T}x\quad \text{ subject to }\quad l\leq \begin{pmatrix}x\\Ax\end{pmatrix}\leq u
where :math:`c` is an :math:`n`-element vector and :math:`A` is an :math:`m\times n` matrix (i.e., there are :math:`n` variables and :math:`m` general linear constraints). :meth:`ilp_dense` is used if at least one of the variables is restricted to take an integer value at the optimum solution.
The document for :meth:`ilp_mpsx_convert` should be consulted for a detailed description of the MPSX format.
In the MPSX data file the first free row, that is a row defined with the row type :math:`\mathrm{n}`, is taken as the objective row.
Similarly, if there are more than one RHS, RANGES or BOUNDS sets, then the first set is used for the optimization. ``ilp_mpsx`` also prints the solution to the problem using the row and column names specified in the MPSX data file (by calling :meth:`ilp_print`).
.. _h02bf-py2-py-references:
**References**
IBM, 1971, `MPSX -- Mathematical programming system`, Program Number 5734 XM4, IBM Trade Corporation, New York
"""
raise NotImplementedError
[docs]def handle_solve_milp(handle, x=None, monit=None, data=None, io_manager=None):
r"""
``handle_solve_milp`` is a solver from the NAG optimization modelling suite for large-scale Mixed Integer Linear Programming (MILP) problems.
It is a branch-and-cut method optimization solver based on the HiGHS software package.
Note: this function uses optional algorithmic parameters, see also: :meth:`opt.handle_opt_set <naginterfaces.library.opt.handle_opt_set>`, :meth:`opt.handle_opt_get <naginterfaces.library.opt.handle_opt_get>`.
.. _h02bk-py2-py-doc:
For full information please refer to the NAG Library document for h02bk
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bkf.html
.. _h02bk-py2-py-parameters:
**Parameters**
**handle** : Handle
The handle to the problem. It needs to be initialized (e.g., by :meth:`opt.handle_init <naginterfaces.library.opt.handle_init>`) and to hold a problem formulation compatible with ``handle_solve_milp``. It **must not** be changed between calls to the NAG optimization modelling suite.
**x** : None or float, array-like, shape :math:`\left(\textit{nvar}\right)`, optional
The input of :math:`\mathrm{x}` is reserved for future releases of the NAG Library and it is ignored at the moment.
**monit** : None or callable monit(handle, rinfo, stats, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{monit}` is reserved for future releases of the NAG Library which will allow you to monitor the progress of the optimization.
It will never be called in the current implementation and :math:`\mathrm{monit}` may be **None**.
**Parameters**
**handle** : Handle
The handle to the problem as provided on entry to ``handle_solve_milp``. It may be used to query the model during the solve, and extract the current approximation of the solution by :meth:`opt.handle_set_get_real <naginterfaces.library.opt.handle_set_get_real>`.
**rinfo** : float, ndarray, shape :math:`\left(100\right)`
Error measures and various indicators at the end of the current iteration as described in :math:`\mathrm{rinfo}`.
**stats** : float, ndarray, shape :math:`\left(100\right)`
Solver statistics at the end of the current iteration as described in :math:`\mathrm{stats}`, however, elements :math:`2`, :math:`3`, :math:`5`, :math:`9`, :math:`10`, :math:`11` and :math:`15` refer to the quantities in the last iteration rather than accumulated over all iterations through the whole algorithm run.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**x** : float, ndarray, shape :math:`\left(\textit{nvar}\right)`
The final values of the variables :math:`x`.
**rinfo** : float, ndarray, shape :math:`\left(100\right)`
Error measures and various indicators of the algorithm as given in the table below:
.. rst-class:: nag-rules-none nag-align-left
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`0` |Value of the primal objective. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |Value of the dual objective bound. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`2` |The absolute value of the gap between the primal and dual bounds, relative to the primal bound. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`3` |The maximum deviation from an integer value over all the discrete variables. This is an indication of the integrality violation.|
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`4` |The maximum violation of a bound on a variable. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`5` |The sum of violations of bounds by variables. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
|:math:`6`--:math:`99`|Reserved for future use. |
+---------------------+--------------------------------------------------------------------------------------------------------------------------------+
**stats** : float, ndarray, shape :math:`\left(100\right)`
Solver statistics as given in the table below.
.. rst-class:: nag-rules-none nag-align-left
+---------------------+---------------------------------------------+
|:math:`0` |Total number of nodes generated. |
+---------------------+---------------------------------------------+
|:math:`1` |Total number of simplex iterations performed.|
+---------------------+---------------------------------------------+
|:math:`2` |Total time spent in the solver. |
+---------------------+---------------------------------------------+
|:math:`3`--:math:`99`|Reserved for future use. |
+---------------------+---------------------------------------------+
.. _h02bk-py2-py-other_params:
**Other Parameters**
**'Defaults'** : valueless
This special keyword may be used to reset all options to their default values.
Any value given with this keyword will be ignored.
**'Infinite Bound Size'** : float
Default :math:`\text{} = 10^{20}`
This defines the 'infinite' bound :math:`\textit{bigbnd}` in the definition of the problem constraints.
Any upper bound greater than or equal to :math:`\textit{bigbnd}` will be regarded as :math:`{+\infty }` (and similarly any lower bound less than or equal to :math:`{-\textit{bigbnd}}` will be regarded as :math:`{-\infty }`).
Note that a modification of this option does not influence constraints which have already been defined; only the constraints formulated after the change will be affected.
Constraint: :math:`\text{‘Infinite Bound Size'}\geq 1000`.
**'Milp Abs Gap'** : float
Default :math:`\text{} = 10^{-6}`
Tolerance on absolute gap, the absolute value of the difference between the primal objective value and the dual objective bound, to determine whether optimality has been reached.
Constraint: :math:`\text{‘Milp Abs Gap'}\geq \epsilon`.
**'Milp Detect Symmetry'** : str
Default :math:`= \texttt{'YES'}`
This argument controls whether symmetry should be detected.
Symmetry arises when there are multiple equivalent solutions that have the same objective function value but different variable assignments.
Identifying and handling symmetry can significantly improve the efficiency and speed of the solver, and reduce memory usage.
However, it's worth noting for smaller problems, the time spent on symmetry detection and breaking may outweigh the benefits.
Constraint: :math:`\text{‘Milp Detect Symmetry'} = \texttt{'YES'}` or :math:`\texttt{'NO'}`.
**'Milp Feasibility Tol'** : float
Default :math:`= 10^{-6}`
The maximum acceptable absolute violation in each constraint (bound, linear or integrality constraint) at a 'feasible' point; i.e., a constraint is considered satisfied if its violation does not exceed 'Milp Feasibility Tol' :math:`r`.
For example, if the integer variable :math:`x_j` is near unity then :math:`x_j` is considered to be integer only if :math:`\left(1-r\right)\leq x_j\leq \left(1+r\right)`.
Constraint: :math:`\text{‘Milp Feasibility Tol'} > 10^{-10}`.
**'Milp Max Nodes'** : int
Default :math:`= \textit{Imax}`
The maximum number of nodes that are to be searched in the B&B tree in order to find a solution.
Constraint: :math:`\text{‘Milp Max Nodes'}\geq 0`.
**'Milp Presolve'** : str
Default :math:`= \texttt{'YES'}`
This argument allows you to turn the presolving of the problem off completely.
If the presolver is turned off, the solver will try to handle the original problem you have given.
In such a case, the presence of linear dependencies in the constraint matrix can cause numerical instabilities to occur.
In normal circumstances, it is recommended to use the presolve which is the default.
Constraint: :math:`\text{‘Milp Presolve'} = \texttt{'YES'}` or :math:`\texttt{'NO'}`.
**'Milp Random Seed'** : int
Default :math:`= 0`
Initial seed used for random selection in tasks such as primal heuristics and cut generation.
Constraint: :math:`\text{‘Milp Random Seed'}\geq 0`.
**'Milp Rel Gap'** : float
Default :math:`\text{} = 10^{-4}`
Tolerance on relative gap, to determine whether optimality has been reached.
Constraint: :math:`\text{‘Milp Rel Gap'}\geq \epsilon`.
**'Milp Small Matrix Value'** : float
Default :math:`\text{} = 10^{-9}`
Lower limit on the absolute value of the linear constraint coefficients in the matrix :math:`A` defined in `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bkf.html#milp>`__.
Values smaller than this will be treated as zero.
Constraint: :math:`\text{‘Milp Small Matrix Value'}\geq 10^{-12}`.
**'Monitoring File'** : int
Default :math:`= -1`
If :math:`i\geq 0`, the unit number for the secondary (monitoring) output.
If set to :math:`-1`, no secondary output is provided.
The following information is output to the unit:
- a listing of the options if set by 'Print Options';
- problem statistics, the iteration log and the final status as set by 'Monitoring Level';
- the solution if set by 'Print Solution'.
Constraint: :math:`\text{‘Monitoring File'}\geq -1`.
**'Monitoring Level'** : int
Default :math:`= 4`
This argument sets the amount of information detail that will be printed by the solver to the secondary output.
The meaning of the levels is the same as with 'Print Level'.
Constraint: :math:`0\leq \text{‘Monitoring Level'}\leq 5`.
**'Print File'** : int
Default :math:`= \text{advisory message unit number}`
If :math:`i\geq 0`, the unit number for the primary output of the solver.
If :math:`\text{‘Print File'} = -1`, the primary output is completely turned off independently of other settings. The default value is the advisory message unit number at the time of the options initialization, e.g., at the initialization of the handle. The following information is output to the unit:
- a listing of options if set by 'Print Options';
- problem statistics, the iteration log and the final status from the solver as set by 'Print Level';
- the solution if set by 'Print Solution'.
Constraint: :math:`\text{‘Print File'} \geq -1`.
**'Print Level'** : int
Default :math:`= 2`
This argument defines how detailed information should be printed by the solver to the primary output.
.. rst-class:: nag-rules-none nag-align-left
+------------------------------------------+-------------------------------+
|:math:`i` |Output |
+==========================================+===============================+
|:math:`0` |No output from the solver |
+------------------------------------------+-------------------------------+
|:math:`1` |The Header and Summary |
+------------------------------------------+-------------------------------+
|:math:`2`, :math:`3`, :math:`4`, :math:`5`|Additionally, the Iteration log|
+------------------------------------------+-------------------------------+
Constraint: :math:`0\leq \text{‘Print Level'}\leq 5`.
**'Print Options'** : str
Default :math:`= \texttt{'YES'}`
If :math:`\text{‘Print Options'} = \texttt{'YES'}`, a listing of options will be printed to the primary and secondary output.
Constraint: :math:`\text{‘Print Options'} = \texttt{'YES'}` or :math:`\texttt{'NO'}`.
**'Print Solution'** : str
Default :math:`= \texttt{'NO'}`
If :math:`\text{‘Print Solution'} = \texttt{'YES'}`, the final values of the primal variables are printed on the primary and secondary outputs.
Constraint: :math:`\text{‘Print Solution'} = \texttt{'YES'}` or :math:`\texttt{'NO'}`.
**'Task'** : str
Default :math:`= \texttt{'MINIMIZE'}`
This argument specifies the required direction of the optimization.
If :math:`\text{‘Task'} = \texttt{'FEASIBLE POINT'}`, the objective function (if set) is ignored and the algorithm stops as soon as a feasible point is found with respect to the given tolerance.
If no objective function is set, 'Task' reverts to 'FEASIBLE POINT' automatically.
Constraint: :math:`\text{‘Task'} = \texttt{'MINIMIZE'}`, :math:`\texttt{'MAXIMIZE'}` or :math:`\texttt{'FEASIBLE POINT'}`.
**'Time Limit'** : float
Default :math:`\text{} = 10^6`
This argument specifies a limit in seconds that the solver can use to solve one problem.
If during the convergence check this limit is exceeded, the solver will terminate with :math:`\mathrm{errno}` = 23 error message.
Constraint: :math:`\text{‘Time Limit'} > 0`.
.. _h02bk-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
:math:`\mathrm{handle}` has not been initialized.
(`errno` :math:`1`)
:math:`\mathrm{handle}` does not belong to the NAG optimization modelling suite, has not been initialized properly or is corrupted.
(`errno` :math:`1`)
:math:`\mathrm{handle}` has not been initialized properly or is corrupted.
(`errno` :math:`2`)
This solver does not support the model defined in the handle.
(`errno` :math:`2`)
The problem is already being solved.
(`errno` :math:`4`)
On entry, :math:`\textit{nvar} = \langle\mathit{\boldsymbol{value}}\rangle`, expected :math:`\mathrm{value} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nvar}` must match the current number of variables of the model in the :math:`\mathrm{handle}`.
(`errno` :math:`51`)
The problem was found to be primal infeasible.
(`errno` :math:`52`)
The problem was found to be dual infeasible.
(`errno` :math:`53`)
The problem seems to be primal or dual infeasible, the algorithm was stopped.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`22`)
Maximum number of nodes exceeded.
**NagAlgorithmicMajorWarning**
(`errno` :math:`23`)
The solver terminated after the maximum time allowed was exhausted.
(`errno` :math:`27`)
No feasible solution was found for the number of nodes investigated in the B&B tree.
.. _h02bk-py2-py-notes:
**Notes**
``handle_solve_milp`` solves a large-scale MILP problem in the following form:
.. math::
\begin{array}{ll} \mathrm{minimize}_{{x \in ℝ^n}} & c^\mathrm{T}x \\ \text{subject to} & l_A \leq Ax \leq u_A \text{,} \\& l_x \leq x \leq u_x \text{,} \\& x \in \mathcal{P} \text{,} \end{array}
where :math:`\mathcal{P} = ℤ^{n_{\textit{int}}}\times ℝ^{n_l}` is a Cartesian product of :math:`n_{\textit{int}}`-dimensional integer space and :math:`n_l`-dimensional real space, and :math:`n = n_{\textit{int}}+n_l` is the number of decision variables.
Here :math:`c`, :math:`x`, :math:`l_x`, :math:`u_x` are :math:`n`-dimensional vectors, :math:`A` is an :math:`m\times n` sparse matrix and :math:`l_A`, :math:`u_A` are :math:`m`-dimensional vectors.
``handle_solve_milp`` solves mixed integer linear programming problems stored as a handle.
The handle points to an internal data structure which defines the problem and serves as a means of communication for functions in the NAG optimization modelling suite.
First, the problem handle is initialized by calling :meth:`opt.handle_init <naginterfaces.library.opt.handle_init>`.
Then some of the functions :meth:`opt.handle_set_linobj <naginterfaces.library.opt.handle_set_linobj>`, :meth:`opt.handle_set_quadobj <naginterfaces.library.opt.handle_set_quadobj>`, :meth:`opt.handle_set_simplebounds <naginterfaces.library.opt.handle_set_simplebounds>` or :meth:`opt.handle_set_linconstr <naginterfaces.library.opt.handle_set_linconstr>` may be called to formulate the objective function, bounds of the variables, and the block of linear constraints, respectively.
Marking a set of variables as being integer-valued is achieved by calling :meth:`opt.handle_set_property <naginterfaces.library.opt.handle_set_property>`.
Once the problem is fully set, the handle may be passed to the solver.
When the handle is not needed anymore, :meth:`opt.handle_free <naginterfaces.library.opt.handle_free>` should be called to destroy it and deallocate the memory held within.
See `the E04 Introduction <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/e04/e04intro.html#optsuite>`__ for more details about the NAG optimization modelling suite.
The solver method can be modified by various options (see :ref:`Other Parameters <h02bk-py2-py-other_params>`) which can be set by :meth:`opt.handle_opt_set <naginterfaces.library.opt.handle_opt_set>` and :meth:`opt.handle_opt_set_file <naginterfaces.library.opt.handle_opt_set_file>` anytime between the initialization of the handle by :meth:`opt.handle_init <naginterfaces.library.opt.handle_init>` and a call to the solver.
Once the solver has finished, options may be modified for the next solve.
The solver may be called repeatedly with various options.
The option 'Task' may be used to switch the problem to maximization or to ignore the objective function and find only a feasible point.
Several options may have significant impact on the performance of the solver.
Even if the defaults were chosen to suit the majority of problems, it is recommended to experiment in order to find the most suitable set of options for a particular problem, see :ref:`Other Parameters <h02bk-py2-py-other_params>` for further details.
.. _h02bk-py2-py-references:
**References**
Dakin, R J, 1965, `A tree search algorithm for mixed integer programming problems`, Comput. J. (8), 250--255
Huangfu, Q, and Hall, J.A. J., 2018, `Parallelizing the dual revised simplex method`, Mathematical Programming Computation (10(1)), 119--142
Mitra, G, 1973, `Investigation of some branch and bound strategies for the solution of mixed integer linear programs`, Math. Programming (4), 155--170
Taha, H A, 1987, `Operations Research: An Introduction`, Macmillan, New York
See Also
--------
:meth:`naginterfaces.library.examples.mip.handle_solve_milp_ex.main`
"""
raise NotImplementedError
[docs]def ilp_mpsx_convert(infile, maxn, maxm, optim, xbldef, xbudef, nmobj, nmrhs, nmrng, nmbnd, mpslst, io_manager=None):
r"""
``ilp_mpsx_convert`` reads data for a linear or integer programming problem from an external file which is in standard or compatible MPSX input format.
.. _h02bu-py2-py-doc:
For full information please refer to the NAG Library document for h02bu
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02buf.html
.. _h02bu-py2-py-parameters:
**Parameters**
**infile** : int
The unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) associated with the MPSX data file.
**maxn** : int
An upper limit for the number of variables in the problem.
**maxm** : int
An upper limit for the number of constraints (including the objective) in the problem.
**optim** : str, length 3
Specifies the direction of the optimization. :math:`\mathrm{optim}` must be set to 'MIN' for minimization and to 'MAX' for maximization.
**xbldef** : float
The default lower bound to be used for the variables in the problem when none is specified in the BOUNDS section of the MPSX data file. For a standard LP or IP problem :math:`\mathrm{xbldef}` would normally be set to zero.
**xbudef** : float
The default upper bound to be used for the variables in the problem when none is specified in the BOUNDS section of the MPSX data file. For a standard LP or IP problem :math:`\mathrm{xbudef}` would normally be set to 'infinity' (i.e., :math:`\mathrm{xbudef}\geq 10^{20}`).
**nmobj** : str, length 8
Either the name of the objective function to be used for the optimization, or blank (in which case the first objective (free) row in the file is used).
**nmrhs** : str, length 8
Either the name of the RHS set to be used for the optimization, or blank (in which case the first RHS set is used).
**nmrng** : str, length 8
Either the name of the RANGE set to be used for the optimization, or blank (in which case the first RANGE set (if any) is used).
**nmbnd** : str, length 8
Either the name of the BOUNDS set to be used for the optimization, or blank (in which case the first BOUNDS set (if any) is used).
**mpslst** : bool
If :math:`\mathrm{mpslst} = \mathbf{True}`, then a listing of the input data is sent to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`). This can be useful for debugging the MPSX data file.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**nmobj** : str, length 8
The name of the objective row as defined in the MPSX data file.
**nmrhs** : str, length 8
The name of the RHS set read in the MPSX data file.
**nmrng** : str, length 8
The name of the RANGE set read in the MPSX data file. This is blank if the MPSX data file does not have a RANGE set.
**nmbnd** : str, length 8
The name of the BOUNDS set read in the MPSX data file. This is blank if the MPSX data file does not have a BOUNDS set.
**n** : int
:math:`n`, the actual number of variables in the problem.
**m** : int
:math:`m`, the actual number of general linear constraints in the problem.
**a** : float, ndarray, shape :math:`\left(\mathrm{maxm}, \mathrm{maxn}\right)`
:math:`A`, the matrix of general linear constraints.
**bl** : float, ndarray, shape :math:`\left(\mathrm{maxn}+\mathrm{maxm}\right)`
:math:`l`, the lower bounds for all the variables and constraints in the following order. The first :math:`\mathrm{n}` elements of :math:`\mathrm{bl}` contain the bounds on the variables and the next :math:`\mathrm{m}` elements contain the bounds for the general linear constraints (if any). Note that an 'infinite' lower bound is indicated by :math:`\mathrm{bl}[j-1] = {-1.0e+20}` and an equality constraint by :math:`\mathrm{bl}[j-1] = \mathrm{bu}[j-1]`.
**bu** : float, ndarray, shape :math:`\left(\mathrm{maxn}+\mathrm{maxm}\right)`
:math:`u`, the upper bounds for all the variables and constraints in the following order. The first :math:`\mathrm{n}` elements of :math:`\mathrm{bu}` contain the bounds on the variables and the next :math:`\mathrm{m}` elements contain the bounds for the general linear constraints (if any). Note that an 'infinite' upper bound is indicated by :math:`\mathrm{bu}[j-1] = 1.0e+20` and an equality constraint by :math:`\mathrm{bu}[j-1] = \mathrm{bl}[j-1]`.
**cvec** : float, ndarray, shape :math:`\left(\mathrm{maxn}\right)`
:math:`c`, the coefficients of the objective function. The signs of these coefficients are determined by the problem (either LP or IP) and the direction of the optimization (see :math:`\mathrm{optim}` above).
**x** : float, ndarray, shape :math:`\left(\mathrm{maxn}\right)`
An initial estimate of the solution to the problem. More precisely, :math:`\mathrm{x}[\textit{j}-1] = 1.0` if :math:`\textit{j}` is odd and :math:`0.0` otherwise, for :math:`\textit{j} = 1,2,\ldots,\mathrm{n}`.
**intvar** : int, ndarray, shape :math:`\left(\mathrm{maxn}\right)`
Indicates which are the integer variables in the problem. More precisely, :math:`\mathrm{intvar}[\textit{k}-1] = 1` if :math:`x_{\textit{k}}` is an integer variable, and :math:`0` otherwise, for :math:`\textit{k} = 1,2,\ldots,\mathrm{n}`.
**crname** : str, length 8, ndarray, shape :math:`\left(\mathrm{maxn}+\mathrm{maxm}\right)`
The MPSX names of all the variables and constraints in the problem in the following order. The first :math:`\mathrm{n}` elements contain the MPSX names for the variables and the next :math:`\mathrm{m}` elements contain the MPSX names for the general linear constraints (if any).
**nmprob** : str, length 8
The name of the problem as defined in the MPSX data file.
.. _h02bu-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
Too many rows. Limit is :math:`\langle\mathit{\boldsymbol{value}}\rangle`, but the actual number required is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`2`)
Too many columns. Limit is :math:`\langle\mathit{\boldsymbol{value}}\rangle`, but the actual number required is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`3`)
No objective function row found.
(`errno` :math:`4`)
No rows specified in the ROWS section.
(`errno` :math:`5`)
Illegal constraint type :math:`\langle\mathit{\boldsymbol{value}}\rangle` was detected in the ROWS section.
(`errno` :math:`6`)
Row name with leading blank or non-alphanumeric character :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`7`)
Column name with leading blank or non-alphanumeric character :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`8`)
Illegal bound type :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`9`)
Column name :math:`\langle\mathit{\boldsymbol{value}}\rangle` is not defined in the COLUMNS section.
(`errno` :math:`10`)
The last line must be the ENDATA indicator line.
(`errno` :math:`11`)
Line :math:`\langle\mathit{\boldsymbol{value}}\rangle` is not a comment nor a valid line.
(`errno` :math:`12`)
Row name :math:`\langle\mathit{\boldsymbol{value}}\rangle` is not defined in the ROWS section.
(`errno` :math:`13`)
No columns specified in the COLUMNS section.
(`errno` :math:`14`)
On entry, :math:`\mathrm{optim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{optim} = \texttt{'MIN'}` or :math:`\texttt{'MAX'}`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{xbldef} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{xbudef} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{xbldef}\leq \mathrm{xbudef}`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{maxm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxm}\geq 1`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{maxn} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxn}\geq 1`.
(`errno` :math:`14`)
On entry, :math:`\mathrm{infile} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{infile}\leq 2147483647`.
(`errno` :math:`15`)
Integer marker is not in a correct position.
.. _h02bu-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``ilp_mpsx_convert`` reads Linear Programming (LP) or integer programming (IP) problem data from an external file which is prepared in standard or compatible MPSX (see IBM (1971)) input format and then initializes :math:`n` (the number of variables), :math:`m` (the number of general linear constraints), the vectors :math:`c`, :math:`l` and :math:`u` and the :math:`m\times n` matrix :math:`A` for use with :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_dense`, which are designed to solve problems of the form
.. math::
\mathrm{minimize}_{{x \in R^n}}c^{\mathrm{T}}x\quad \text{ subject to }\quad l\leq \begin{pmatrix}x\\Ax\end{pmatrix}\leq u\text{.}
``ilp_mpsx_convert`` may be followed by calls to either :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` (to solve an LP problem) or :meth:`ilp_dense` and :meth:`ilp_info` (to solve an IP problem), possibly followed by a call to :meth:`ilp_print` (to print the solution using MPSX names).
Note that ``ilp_mpsx_convert`` uses an 'infinite' bound size of :math:`10^{20}` in the definition of :math:`l` and :math:`u`.
In other words, any element of :math:`u` greater than or equal to :math:`10^{20}` will be regarded as :math:`{+\infty }` (and similarly any element of :math:`l` less than or equal to :math:`{-10^{20}}` will be regarded as :math:`{-\infty }`).
If this value is deemed to be 'inappropriate', you are recommended to reset the value of either the option 'Infinite Bound Size' (if an LP problem is being solved) or the argument BIGBND (if an IP problem is being solved) and make any necessary changes to :math:`\mathrm{bl}` and/or :math:`\mathrm{bu}` prior to calling :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_dense` (as appropriate).
The documents for :meth:`ilp_print` , :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` and/or :meth:`ilp_dense` and :meth:`ilp_info` should be consulted for further details.
**MPSX input format**
The input file of data may only contain two types of lines.
(1) Indicator lines (specifying the type of data which is to follow).
(#) Data lines (specifying the actual data).
The input file must not contain any blank lines.
Any characters beyond column 80 are ignored.
Indicator lines must not contain leading blank characters (in other words they must begin in column 1).
The following displays the order in which the indicator lines must appear in the file:
.. rst-class:: nag-rules-none nag-align-left
+-----------------+---------------+
|NAME |user-given name|
+-----------------+---------------+
|ROWS |data line(s) |
+-----------------+---------------+
|COLUMNS |data line(s) |
+-----------------+---------------+
|RHS |data line(s) |
+-----------------+---------------+
|RANGES (optional)|data line(s) |
+-----------------+---------------+
|BOUNDS (optional)|data line(s) |
+-----------------+---------------+
|ENDATA | |
+-----------------+---------------+
The 'user-given name' specifies a name for the problem and must occupy columns 15--22.
The name can either be blank or up to a maximum of :math:`8` characters.
A data line follows the same fixed format made up of fields defined below.
The contents of the fields may have different significance depending upon the section of data in which they appear.
+--------+-----------+------------+-------------+-------------+-------------+-------------+
| |Field 1 |Field 2 |Field 3 |Field 4 |Field 5 |Field 6 |
+========+===========+============+=============+=============+=============+=============+
|Columns |:math:`2-3`|:math:`5-12`|:math:`15-22`|:math:`25-36`|:math:`40-47`|:math:`50-61`|
+--------+-----------+------------+-------------+-------------+-------------+-------------+
|Contents|Code |Name |Name |Value |Name |Value |
+--------+-----------+------------+-------------+-------------+-------------+-------------+
The names and codes consist of 'alphanumeric' characters (i.e., a--z, A--Z, :math:`0`--:math:`9`, :math:`+`, :math:`-`, asterisk (*), blank (), colon (:), dollar sign ($) or fullstop (.) only) and the names must not contain leading blank characters.
Values are read using Fortran format E12.0.
This allows values to be entered in several equivalent forms.
For example, :math:`1.2345678`, :math:`1.2345678E+0`, :math:`123.45678e-2` and :math:`12345678e-07` all represent the same number.
It is safest to include an explicit decimal point.
Note that in order to ensure numeric values are interpreted as intended, they should be `right-justified` in the :math:`12`-character field, with no trailing blanks.
This is because in some situations trailing blanks may be interpreted as zeros and this can dramatically affect the interpretation of the value.
This is relevant if the value contains an exponent, or if it contains neither an exponent nor an explicit decimal point.
For example, the fields
::
%%%%1.23e-2%
%%%%%%%123%%
may be interpreted as :math:`1.23e-20` and :math:`12300` respectively (where % is used to denote a blank).
The actual behaviour is system-dependent.
Comment lines are allowed in the data file.
These must have an asterisk (``*``) in column 1 and any characters in columns 2--80.
In any data line, a dollar sign (``$``) as the first character in Field 3 or 5 indicates that the information from that point through column 80 consists of comments.
Columns outside the six fields must be blank, except for columns 72--80, whose contents are ignored by the function.
These columns may be used to enter a sequence number.
A non-blank character outside the predefined six fields and columns 72--80 is considered to be a major error (:math:`\mathrm{errno}` = 11; see :ref:`Exceptions <h02bu-py2-py-errors>`), unless it is part of a comment.
**ROWS data line(s)**
These lines specify row (constraint) names and their inequality types (i.e., :math:`=`, :math:`\geq` or :math:`\leq`).
.. rst-class:: nag-rules-none nag-align-left
+--------+----------------------------------------------------------------------------+
|Field 1:|defines the constraint type. It may be in column 2 or column 3. |
+--------+----------------------------------------------------------------------------+
|``N`` |free row, that is no constraint. It may be used to define the objective row.|
+--------+----------------------------------------------------------------------------+
|``G`` |greater than or equal to (i.e., :math:`\geq`). |
+--------+----------------------------------------------------------------------------+
|``L`` |less than or equal to (i.e., :math:`\leq`). |
+--------+----------------------------------------------------------------------------+
|``E`` |exactly equal to (i.e., :math:`=`). |
+--------+----------------------------------------------------------------------------+
|Field 2:|defines the row name. |
+--------+----------------------------------------------------------------------------+
Row type ``N`` stands for 'Not binding', also known as 'Free'.
It can be used to define the objective row.
The objective row is a free row that specifies the vector :math:`c` in the objective function.
It is taken to be the first free row, unless some other free row name is specified by the argument :math:`\mathrm{nmobj}` (see :ref:`Parameters <h02bu-py2-py-parameters>`).
Note that the objective function must be included in the MPSX data file.
Thus the maximum number of constraints (:math:`\mathrm{maxm}`; see :ref:`Parameters <h02bu-py2-py-parameters>`) in the problem must be :math:`m+1`.
**COLUMNS data line(s)**
These lines specify the names to be assigned to the variables (columns) in the constraint matrix :math:`A`, and define, in terms of column vectors, the actual values of the corresponding matrix elements.
.. rst-class:: nag-rules-none nag-align-left
+--------+--------------------------------------------------------------------------------------------+
|Field 1:|blank (ignored) |
+--------+--------------------------------------------------------------------------------------------+
|Field 2:|gives the name of the column associated with the elements specified in the following fields.|
+--------+--------------------------------------------------------------------------------------------+
|Field 3:|contains the name of a row. |
+--------+--------------------------------------------------------------------------------------------+
|Field 4:|used in conjunction with Field 3 contains the value of the matrix element. |
+--------+--------------------------------------------------------------------------------------------+
|Field 5:|is optional (may be used like Field 3). |
+--------+--------------------------------------------------------------------------------------------+
|Field 6:|is optional (may be used like Field 4). |
+--------+--------------------------------------------------------------------------------------------+
Note that only nonzero elements of :math:`A` need to be specified in the COLUMNS section, as any unspecified elements are assumed to be zero.
**RHS data line(s)**
This section specifies the right-hand side values of the constraint matrix :math:`A`.
The lines specify the name of the RHS (right-hand side) vector given to the problem, the numerical values of the elements of the vector are also defined by the data lines and may appear in any order.
The data lines have exactly the same format as the COLUMNS data lines, except that the column name is replaced by the RHS name.
Note that any unspecified elements are assumed to be zero.
**RANGES data line(s) (optional)**
Ranges are used for constraints of the form :math:`l\leq Ax\leq u`, where :math:`l` and :math:`u` are finite.
The range of the constraint is :math:`r = u-l`.
Either :math:`l` or :math:`u` must be specified in the RHS section and :math:`r` must be defined in this section.
The data lines have exactly the same format as the COLUMNS data lines, except that the column name is replaced by the RANGES name.
**BOUNDS data line(s) (optional)**
These lines specify limits on the values of the variables (:math:`l` and :math:`u` in :math:`l\leq x\leq u`).
If the variable is not specified in the bound set then it is automatically assumed to lie between default lower and upper bounds (usually :math:`0` and :math:`{+\infty }`).
Like an RHS column which is given a name, the set of variables in one bound set is also given a name.
.. rst-class:: nag-rules-none nag-align-left
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 1:|specifies the type of bound or defines the variable type. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|LO |lower bound |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|UP |upper bound |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|FX |fixed variable |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|FR |free variable (:math:`{-\infty }` to :math:`{+\infty }`) |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|MI |lower bound is :math:`{-\infty }` |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|PL |upper bound is :math:`{+\infty }`. This is the default variable type. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 2:|identifies a name for the bound set. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 3:|identifies the column name of the variable belonging to this set. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 4:|identifies the value of the bound; this has a numerical value only in association with LO, UP, FX in Field 1, otherwise it is blank.|
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 5:|is blank and ignored. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
|Field 6:|is blank and ignored. |
+--------+------------------------------------------------------------------------------------------------------------------------------------+
Note that if RANGES and BOUNDS sections are both present, the RANGES section must appear first.
**Integer Problems**
In IP problems there are two common integer variable types.
(1) :math:`0`--1 integer variables which represent 'on' or 'off' situations and
(#) General integer variables which are forced to take an integer value, in a specified range, at the optimal integer solution.
Integer variables can be defined in the following compatible and standard MPSX forms.
In the compatible MPSX format, the type of integer variables is defined in Field 1 of the BOUNDS section, that is:
.. rst-class:: nag-rules-none nag-align-left
+--------+----------------------------------------------------------+
|Field 1:|specifies the type of the integer variable. |
+--------+----------------------------------------------------------+
|BV |:math:`0-1` integer variable (bound value is :math:`1.0`).|
+--------+----------------------------------------------------------+
|UI |general integer variable (bound value is in Field 4). |
+--------+----------------------------------------------------------+
In the standard MPSX format, the integer variables are treated the same as the 'ordinary' bounded variables, in the BOUNDS section.
Integer markers are, however, introduced in the COLUMNS section to specify the integer variables.
The indicator lines for these markers are:
+--------+-----------+------------+-------------+-------------+-------------+-------------+
| |Field 1 |Field 2 |Field 3 |Field 4 |Field 5 |Field 6 |
+========+===========+============+=============+=============+=============+=============+
|Columns |:math:`2-3`|:math:`5-12`|:math:`15-22`|:math:`25-36`|:math:`40-47`|:math:`50-61`|
+--------+-----------+------------+-------------+-------------+-------------+-------------+
|Contents| |INTEGER |'MARKER' | |'INTORG' | |
+--------+-----------+------------+-------------+-------------+-------------+-------------+
to mark the beginning of the integer variables and
+--------+-----------+------------+-------------+-------------+-------------+-------------+
| |Field 1 |Field 2 |Field 3 |Field 4 |Field 5 |Field 6 |
+========+===========+============+=============+=============+=============+=============+
|Columns |:math:`2-3`|:math:`5-12`|:math:`15-22`|:math:`25-36`|:math:`40-47`|:math:`50-61`|
+--------+-----------+------------+-------------+-------------+-------------+-------------+
|Contents| |INTEGER |'MARKER' | |'INTEND' | |
+--------+-----------+------------+-------------+-------------+-------------+-------------+
to mark the end.
That is, any variables between these markers are treated as integer variables.
Note that if the (INTEND) indicator line is not specified in the file then all the variables between the (INTORG) indicator line and the end of the COLUMNS section are assumed to be integer variables.
The function accepts both standard and/or compatible MPSX format as a means of specifying integer variables.
.. _h02bu-py2-py-references:
**References**
IBM, 1971, `MPSX -- Mathematical programming system`, Program Number 5734 XM4, IBM Trade Corporation, New York
"""
raise NotImplementedError
[docs]def ilp_print(a, bl, bu, x, clamda, istate, crname, comm, io_manager=None):
r"""
``ilp_print`` prints the solution to a linear or integer programming problem computed by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_dense` and :meth:`ilp_info`, with user-supplied names for the rows and columns.
.. _h02bv-py2-py-doc:
For full information please refer to the NAG Library document for h02bv
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bvf.html
.. _h02bv-py2-py-parameters:
**Parameters**
**a** : float, array-like, shape :math:`\left(m, :\right)`
Note: the required extent for this argument in dimension 2 is determined as follows: if :math:`m > 0`: :math:`n`; if :math:`m=0`: :math:`1`; otherwise: :math:`0`.
The matrix of general linear constraints, as returned by :meth:`ilp_mpsx_convert`.
**bl** : float, array-like, shape :math:`\left(n+m\right)`
The lower bounds for all the constraints, as returned by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_info`.
**bu** : float, array-like, shape :math:`\left(n+m\right)`
The upper bounds for all the constraints, as returned by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_info`.
**x** : float, array-like, shape :math:`\left(n\right)`
The solution to the problem, as returned by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_dense`.
**clamda** : float, array-like, shape :math:`\left(n+m\right)`
The Lagrange-multipliers (reduced costs) for each constraint with respect to the working set, as returned by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_info`.
**istate** : int, array-like, shape :math:`\left(n+m\right)`
The status of every constraint in the working set at the solution, as returned by :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` or :meth:`ilp_info`.
**crname** : str, length 8, array-like, shape :math:`\left(n+m\right)`
The user-defined names for all the variables and constraints, as returned by :meth:`ilp_mpsx_convert`.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.lp_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
.. _h02bv-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m\geq 0`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
.. _h02bv-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``ilp_print`` prints the solution to a linear or integer programming problem with user-supplied names for the rows and columns.
All output is written to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`).
The function must be preceded in the same program by calls to :meth:`ilp_mpsx_convert` and either :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` (if an LP problem has been solved) or :meth:`ilp_dense` and :meth:`ilp_info` (if an IP problem has been solved).
The documents for :meth:`ilp_mpsx_convert` , :meth:`opt.lp_solve <naginterfaces.library.opt.lp_solve>` and/or :meth:`ilp_dense` and :meth:`ilp_info` should be consulted for further details.
.. _h02bv-py2-py-references:
**References**
IBM, 1971, `MPSX -- Mathematical programming system`, Program Number 5734 XM4, IBM Trade Corporation, New York
"""
raise NotImplementedError
[docs]def ilp_info(n, m, comm):
r"""
``ilp_info`` extracts more information associated with the solution of an integer programming problem computed by :meth:`ilp_dense`.
.. deprecated:: 29.3.0.0
``ilp_info`` is deprecated.
Please use :meth:`handle_solve_milp` instead.
See also the :ref:`Replacement Calls <replace>` document.
.. _h02bz-py2-py-doc:
For full information please refer to the NAG Library document for h02bz
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02bzf.html
.. _h02bz-py2-py-parameters:
**Parameters**
**n** : int
This **must** be the same argument :math:`\mathrm{n}` as supplied to :meth:`ilp_dense`.
**m** : int
This **must** be the same argument :math:`\mathrm{m}` as supplied to :meth:`ilp_dense`.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by a prior call to :meth:`ilp_dense`.
**Returns**
**bl** : float, ndarray, shape :math:`\left(\mathrm{n}+\mathrm{m}\right)`
If :meth:`ilp_dense` exits with the function exits successfully or :math:`\mathrm{errno}` = 7 or 9, the values in the array :math:`\mathrm{bl}` contain the lower bounds imposed on the integer solution for all the constraints. The first :math:`\mathrm{n}` elements contain the lower bounds on the variables, and the next :math:`\mathrm{m}` elements contain the lower bounds for the general linear constraints (if any).
**bu** : float, ndarray, shape :math:`\left(\mathrm{n}+\mathrm{m}\right)`
If :meth:`ilp_dense` exits with the function exits successfully or :math:`\mathrm{errno}` = 7 or 9, the values in the array :math:`\mathrm{bu}` contain the upper bounds imposed on the integer solution for all the constraints. The first :math:`\mathrm{n}` elements contain the upper bounds on the variables, and the next :math:`\mathrm{m}` elements contain the upper bounds for the general linear constraints (if any).
**clamda** : float, ndarray, shape :math:`\left(\mathrm{n}+\mathrm{m}\right)`
If :meth:`ilp_dense` exits with the function exits successfully or :math:`\mathrm{errno}` = 7 or 9, the values in the array :math:`\mathrm{clamda}` contain the values of the Lagrange-multipliers for each constraint with respect to the current working set. The first :math:`\mathrm{n}` elements contain the multipliers (reduced costs) for the bound constraints on the variables, and the next :math:`\mathrm{m}` elements contain the multipliers (shadow costs) for the general linear constraints (if any).
**istate** : int, ndarray, shape :math:`\left(\mathrm{n}+\mathrm{m}\right)`
If :meth:`ilp_dense` exits with the function exits successfully or :math:`\mathrm{errno}` = 7 or 9, the values in the array :math:`\mathrm{istate}` indicate the status of the constraints in the working set at an integer solution. Otherwise, :math:`\mathrm{istate}` indicates the composition of the working set at the final iterate. The significance of each possible value of :math:`\mathrm{istate}[j]` is as follows.
.. rst-class:: nag-rules-none nag-align-left
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\mathrm{istate}[j]`|Meaning |
+==========================+==================================================================================================================================================================================================================================================================================+
|:math:`-2` |The constraint violates its lower bound by more than :math:`\textit{tolfes}` (the feasibility tolerance, see :meth:`ilp_dense`). |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`-1` |The constraint violates its upper bound by more than :math:`\textit{tolfes}`. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`0` |The constraint is satisfied to within :math:`\textit{tolfes}`, but is not in the working set. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |This inequality constraint is included in the working set at its lower bound. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`2` |This inequality constraint is included in the working set at its upper bound. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`3` |This constraint is included in the working set as an equality. This value of :math:`\mathrm{istate}` can occur only when :math:`\mathrm{bl}[j] = \mathrm{bu}[j]`. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`4` |This corresponds to an integer solution being declared with :math:`x_j` being temporarily fixed at its current value. This value of :math:`\mathrm{istate}` can occur only when the function exits successfully or :math:`\mathrm{errno}` = 7 or 9 on exit from :meth:`ilp_dense`.|
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
.. _h02bz-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n} > 0`.
.. _h02bz-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``ilp_info`` extracts the following information associated with the solution of an integer programming problem computed by :meth:`ilp_dense`.
The upper and lower bounds used for the solution, the Lagrange-multipliers (costs), and the status of the variables at the solution.
In the branch and bound method employed by :meth:`ilp_dense`, the arrays :math:`\mathrm{bl}` and :math:`\mathrm{bu}` are used to impose restrictions on the values of the integer variables in each sub-problem.
That is, if the variable :math:`x_j` is restricted to take value :math:`v_j` in a particular sub-problem, then :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = v_j` is set in the sub-problem.
Thus, on exit from this function, some of the elements of :math:`\mathrm{bl}` and :math:`\mathrm{bu}` which correspond to integer variables may contain these imposed values, rather than those originally supplied to :meth:`ilp_dense`.
"""
raise NotImplementedError
[docs]def iqp_dense(bl, bu, h, intvar, xs, strtgy, comm, a=None, cvec=None, qphess=None, mdepth=None, istate=None, monit=None, data=None, io_manager=None):
r"""
``iqp_dense`` solves general quadratic programming problems with integer constraints on the variables.
It is not intended for large sparse problems.
Note: this function uses optional algorithmic parameters, see also: :meth:`iqp_dense_optfile`, :meth:`iqp_dense_optstr`.
.. _h02cb-py2-py-doc:
For full information please refer to the NAG Library document for h02cb
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html
.. _h02cb-py2-py-parameters:
**Parameters**
**bl** : float, array-like, shape :math:`\left(n+\textit{nclin}\right)`
:math:`\mathrm{bl}` must contain the lower bounds and :math:`\mathrm{bu}` the upper bounds, for all the constraints in the following order. The first :math:`n` elements of each array must contain the bounds on the variables, and the next :math:`m_L` elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., :math:`l_j = {-\infty }`), set :math:`\mathrm{bl}[j]\leq -\textit{bigbnd}`, and to specify a nonexistent upper bound (i.e., :math:`u_j = {+\infty }`), set :math:`\mathrm{bu}[j]\geq \textit{bigbnd}`; the default value of :math:`\textit{bigbnd}` is :math:`10^{20}`, but this may be changed by the 'Infinite Bound Size'. To specify the :math:`j`\ th constraint as an `equality`, set :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = \beta`, say, where :math:`\left\lvert \beta \right\rvert < \textit{bigbnd}`.
**bu** : float, array-like, shape :math:`\left(n+\textit{nclin}\right)`
:math:`\mathrm{bl}` must contain the lower bounds and :math:`\mathrm{bu}` the upper bounds, for all the constraints in the following order. The first :math:`n` elements of each array must contain the bounds on the variables, and the next :math:`m_L` elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., :math:`l_j = {-\infty }`), set :math:`\mathrm{bl}[j]\leq -\textit{bigbnd}`, and to specify a nonexistent upper bound (i.e., :math:`u_j = {+\infty }`), set :math:`\mathrm{bu}[j]\geq \textit{bigbnd}`; the default value of :math:`\textit{bigbnd}` is :math:`10^{20}`, but this may be changed by the 'Infinite Bound Size'. To specify the :math:`j`\ th constraint as an `equality`, set :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = \beta`, say, where :math:`\left\lvert \beta \right\rvert < \textit{bigbnd}`.
**h** : float, array-like, shape :math:`\left(:, :\right)`
May be used to store the quadratic term :math:`H` of the QP objective function if desired. In some cases, you need not use :math:`\mathrm{h}` to store :math:`H` explicitly (see the specification of :math:`\mathrm{qphess}`). The elements of :math:`\mathrm{h}` are referenced only by :math:`\mathrm{qphess}`. The number of rows of :math:`H` is denoted by :math:`m`, whose default value is :math:`n`. (The 'Hessian Rows' may be used to specify a value of :math:`m < n`.)
If the default version of :math:`\mathrm{qphess}` is used and the problem is of type QP1 or QP2 (the default), the first :math:`m` rows and columns of :math:`\mathrm{h}` must contain the leading :math:`m\times m` rows and columns of the symmetric Hessian matrix :math:`H`.
Only the diagonal and upper triangular elements of the leading :math:`m` rows and columns of :math:`\mathrm{h}` are referenced.
The remaining elements need not be assigned.
If the default version of :math:`\mathrm{qphess}` is used and the problem is of type QP3 or QP4, the first :math:`m` rows of :math:`\mathrm{h}` must contain an :math:`m\times n` upper trapezoidal factor of the symmetric Hessian matrix :math:`H^\mathrm{T}H`.
The factor need not be of full rank, i.e., some of the diagonal elements may be zero.
However, as a general rule, the larger the dimension of the leading nonsingular sub-matrix of :math:`\mathrm{h}`, the fewer iterations will be required.
Elements outside the upper trapezoidal part of the first :math:`m` rows of :math:`\mathrm{h}` need not be assigned.
In other situations, it may be desirable to compute :math:`Hx` or :math:`H^\mathrm{T}Hx` without accessing :math:`\mathrm{h}` -- for example, if :math:`H` or :math:`H^\mathrm{T}H` is sparse or has special structure.
The arguments :math:`\mathrm{h}` and :math:`\textit{ldh}` may then refer to any convenient array.
If the problem is of type FP or LP, :math:`\mathrm{h}` is not referenced.
**intvar** : int, array-like, shape :math:`\left(\textit{lintvr}\right)`
:math:`\mathrm{intvar}[i]` must contain the index of the solution vector :math:`x` which is required to be integer. For example, if :math:`x_1` and :math:`x_3` are constrained to take integer values then :math:`\mathrm{intvar}[0]` might be set to :math:`1` and :math:`\mathrm{intvar}[1]` to :math:`3`. The order in which the indices are specified is important, since this determines the order in which the sub-problems are generated. As a rule-of-thumb, the important variables should always be specified first. Thus, in the above example, if :math:`x_3` relates to a more important quantity than :math:`x_1`, then it might be advantageous to set :math:`\mathrm{intvar}[0] = 3` and :math:`\mathrm{intvar}[1] = 1`. If :math:`k` is the smallest integer such that :math:`\mathrm{intvar}[k]` is less than or equal to zero then ``iqp_dense`` assumes that :math:`k-1` variables are constrained to be integer; components :math:`\mathrm{intvar}[k+1]`, :math:`\ldots`, :math:`\mathrm{intvar}[\textit{lintvr}-1]` are `not` referenced.
**xs** : float, array-like, shape :math:`\left(n\right)`
An initial estimate of the solution.
**strtgy** : int
Determines a branching strategy to be used throughout the computation, as follows:
.. rst-class:: nag-rules-none nag-align-left
+-----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\mathrm{strtgy}`|Meaning |
+=======================+================================================================================================================================================================================================================+
|:math:`0` |Always left branch first, i.e., impose an upper bound constraint on the variable first. |
+-----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |Always right branch first, i.e., impose a lower bound constraint on the variable first. |
+-----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`2` |Branch towards the nearest integer, i.e., if :math:`x_k = 2.4` then impose an upper bound constraint :math:`x_k\leq 2`, whereas if :math:`x_k = 2.6` then impose the lower bound constraint :math:`x_k\geq 3.0`.|
+-----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`3` |A random choice is made between a left-hand and a right-hand branch. |
+-----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qp_dense_solve'}`.
**a** : None or float, array-like, shape :math:`\left(\textit{nclin}, :\right)`, optional
Note: the required extent for this argument in dimension 2 is determined as follows: if :math:`\textit{nclin} > 0`: :math:`n`; if :math:`\textit{nclin}=0`: :math:`1`; otherwise: :math:`0`.
The :math:`\textit{i}`\ th row of :math:`\mathrm{a}` must contain the coefficients of the :math:`\textit{i}`\ th general linear constraint, for :math:`\textit{i} = 1,2,\ldots,m_L`.
**cvec** : None or float, array-like, shape :math:`\left(n\right)`, optional
The coefficients of the explicit linear term of the objective function when the problem is of type LP, QP2 (the default) and QP4.
If the problem is of type FP, QP1, or QP3, :math:`\mathrm{cvec}` is not referenced.
**qphess** : None or callable hx = qphess(jthcol, h, x, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
In general, you need not provide a version of :math:`\mathrm{qphess}`, because supplying **None** will cause a 'default' function to be used.
However, the algorithm of ``iqp_dense`` requires only the product of :math:`H` or :math:`H^\mathrm{T}H` and a vector :math:`x`; and in some cases you may obtain increased efficiency by providing a version of :math:`\mathrm{qphess}` that avoids the need to define the elements of the matrices :math:`H` or :math:`H^\mathrm{T}H` explicitly. :math:`\mathrm{qphess}` is not referenced if the problem is of type FP or LP, in which case :math:`\mathrm{qphess}` may be **None**.
**Parameters**
**jthcol** : int
Specifies whether or not the vector :math:`x` is a column of the identity matrix.
:math:`\mathrm{jthcol} = j > 0`
The vector :math:`x` is the :math:`j`\ th column of the identity matrix, and hence :math:`Hx` or :math:`\mathrm{h}^\mathrm{T}Hx` is the :math:`j`\ th column of :math:`\mathrm{h}` or :math:`\mathrm{h}^\mathrm{T}H`, respectively, which may in some cases require very little computation and :math:`\mathrm{qphess}` may be coded to take advantage of this. However special code is not necessary because :math:`x` is always stored explicitly in the array :math:`\mathrm{x}`.
:math:`\mathrm{jthcol} = 0`
:math:`x` has no special form.
**h** : float, ndarray, shape :math:`\left(:, :\right)`
This is the same argument :math:`\textit{h}` as supplied to ``iqp_dense``.
**x** : float, ndarray, shape :math:`\left(n\right)`
The vector :math:`x`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**hx** : float, array-like, shape :math:`\left(n\right)`
The product :math:`Hx` if the problem is of type QP1 or QP2 (the default), or the product :math:`\mathrm{h}^\mathrm{T}Hx` if the problem is of type QP3 or QP4.
**mdepth** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:`3\times n/2`.
The maximum depth (i.e., number of extra constraints) that ``iqp_dense`` may insert before admitting failure.
**istate** : None or int, array-like, shape :math:`\left(n+\textit{nclin}\right)`, optional
Need not be set if the (default) option 'Cold Start' is used.
If the option 'Warm Start' has been chosen, :math:`\mathrm{istate}` specifies the desired status of the constraints at the start of the feasibility phase.
More precisely, the first :math:`n` elements of :math:`\mathrm{istate}` refer to the upper and lower bounds on the variables, and the next :math:`m_L` elements refer to the general linear constraints (if any).
Possible values for :math:`\mathrm{istate}[j]` are as follows:
.. rst-class:: nag-rules-none nag-align-left
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\mathrm{istate}[j]`|Meaning |
+==========================+====================================================================================================================================================+
|0 |The corresponding constraint should `not` be in the initial working set. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+
|1 |The constraint should be in the initial working set at its lower bound. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+
|2 |The constraint should be in the initial working set at its upper bound. |
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+
|3 |The constraint should be in the initial working set as an equality. This value must not be specified unless :math:`\mathrm{bl}[j] = \mathrm{bu}[j]`.|
+--------------------------+----------------------------------------------------------------------------------------------------------------------------------------------------+
The values :math:`-2`, :math:`-1` and :math:`4` are also acceptable but will be reset to zero by the function.
If ``iqp_dense`` has been called previously with the same values of :math:`\textit{n}` and :math:`\textit{nclin}`, :math:`\mathrm{istate}` already contains satisfactory information. (See also the description of the option 'Warm Start'.) The function also adjusts (if necessary) the values supplied in :math:`\mathrm{xs}` to be consistent with :math:`\mathrm{istate}`.
**monit** : None or callable (bstval, halt, count) = monit(intfnd, nodes, depth, obj, x, bstval, bstsol, bl, bu, halt, count, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{monit}` may be used to print out intermediate output and to affect the course of the computation.
Specifically, it allows you to specify a realistic value for the cut-off value (see :ref:`Notes <h02cb-py2-py-notes>`) and to terminate the algorithm.
If you do not require any intermediate output, have no estimate of the cut-off value and require an exhaustive tree search then :math:`\mathrm{monit}` may be **None**.
**Parameters**
**intfnd** : int
Specifies the number of integer solutions obtained so far.
**nodes** : int
Specifies the number of nodes (sub-problems) solved so far.
**depth** : int
Specifies the depth in the tree of sub-problems the algorithm has now reached.
**obj** : float
Specifies the value of the objective function of the end of the latest sub-problem.
**x** : float, ndarray, shape :math:`\left(n\right)`
Specifies the values of the independent variables at the end of the latest sub-problem.
**bstval** : float
Normally specifies the value of the best integer solution found so far.
**bstsol** : float, ndarray, shape :math:`\left(n\right)`
Specifies the solution vector which gives rise to the best integer solution value so far discovered.
**bl** : float, ndarray, shape :math:`\left(n\right)`
:math:`\mathrm{bl}[i]` specifies the current lower bounds on the variable :math:`x_i`.
**bu** : float, ndarray, shape :math:`\left(n\right)`
:math:`\mathrm{bu}[i]` specifies the current upper bounds on the variable :math:`x_i`.
**halt** : bool
Will have the value :math:`\mathbf{False}`.
**count** : int
Unchanged from previous call.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**bstval** : float
May be set to a cut-off value if you are an experienced user as follows. Before an integer solution has been found :math:`\mathrm{bstval}` will be set by ``iqp_dense`` to the largest machine representable number (see :meth:`machine.real_largest <naginterfaces.library.machine.real_largest>`). If you know that the solution being sought is a much smaller number, then :math:`\mathrm{bstval}` may be set to this number as a cut-off value (see :ref:`Notes <h02cb-py2-py-notes>`). Beware of setting :math:`\mathrm{bstval}` too small, since then no integer solutions will be discovered. Also, on entry to :math:`\mathrm{monit}`, :math:`\mathrm{bstval}` should be set using a statement of the form
::
if intfnd == 0: bstval = cut_off_value
Setting :math:`\mathrm{bstval}` in this way will not prevent the normal operation of the algorithm when subsequent integer solutions are found. It would be a grievous mistake to unconditionally set :math:`\mathrm{bstval}` and if you have any doubts whatsoever about the correct use of this argument then you are strongly recommended to leave it unchanged.
**halt** : bool
If :math:`\mathrm{halt}` is set to :math:`\mathbf{True}`, :meth:`opt.qp_dense_solve <naginterfaces.library.opt.qp_dense_solve>` will be brought to a halt with :math:`\mathrm{errno}` = -1. This facility may be useful if you are content with `any` integer solution, or with any integer solution that fits certain criteria. Under these circumstances setting :math:`\mathrm{halt} = \mathbf{True}` can save considerable unnecessary computation.
**count** : int
May be used by you to save the last value of :math:`\mathrm{intfnd}`. If a subsequent call of :math:`\mathrm{monit}` has a value of :math:`\mathrm{intfnd}` which is greater than :math:`\mathrm{count}`, then you know that a new integer solution has been found at this node.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**istate** : int, ndarray, shape :math:`\left(n+\textit{nclin}\right)`
The status of the constraints in the working set at the point returned in :math:`\mathrm{xs}`. The significance of each possible value of :math:`\mathrm{istate}[j]` is as follows:
.. rst-class:: nag-rules-none nag-align-left
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\mathrm{istate}[j]`|Meaning |
+==========================+=====================================================================================================================================================================================================================+
|:math:`-2` |The constraint violates its lower bound by more than the feasibility tolerance. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`-1` |The constraint violates its upper bound by more than the feasibility tolerance. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`0` |The constraint is satisfied to within the feasibility tolerance, but is not in the working set. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |This inequality constraint is included in the working set at its lower bound. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`2` |This inequality constraint is included in the working set at its upper bound. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`3` |This constraint is included in the working set as an equality. This value of :math:`\mathrm{istate}` can occur only when :math:`\mathrm{bl}[j] = \mathrm{bu}[j]`. |
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`4` |This corresponds to optimality being declared with :math:`\mathrm{xs}[j]` being temporarily fixed at its current value. This value of :math:`\mathrm{istate}` can occur only when :math:`\mathrm{errno}` = 1 on exit.|
+--------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
**xs** : float, ndarray, shape :math:`\left(n\right)`
The point at which ``iqp_dense`` terminated. If the function exits successfully or :math:`\mathrm{errno}` = 1 or 3, :math:`\mathrm{xs}` contains an estimate of the solution.
**obj** : float
The value of the objective function at :math:`x` if :math:`x` is feasible, or the sum of infeasibilities at :math:`x` otherwise. If the problem is of type FP and :math:`x` is feasible, :math:`\mathrm{obj}` is set to zero.
**ax** : None or float, ndarray, shape :math:`\left(\max\left(1,\textit{nclin}\right)\right)`
The final values of the linear constraints :math:`Ax`.
**clamda** : float, ndarray, shape :math:`\left(n+\textit{nclin}\right)`
The values of the Lagrange-multipliers for each constraint with respect to the current working set. The first :math:`n` elements contain the multipliers for the bound constraints on the variables, and the next :math:`m_L` elements contain the multipliers for the general linear constraints (if any). If :math:`\mathrm{istate}[j] = 0` (i.e., constraint :math:`j` is not in the working set), :math:`\mathrm{clamda}[j]` is zero. If :math:`x` is optimal, :math:`\mathrm{clamda}[j]` should be non-negative if :math:`\mathrm{istate}[j] = 1`, non-positive if :math:`\mathrm{istate}[j] = 2` and zero if :math:`\mathrm{istate}[j] = 4`.
.. _h02cb-py2-py-other_params:
**Other Parameters**
**'Check Frequency'** : int
Default :math:`\text{} = 50`
Every :math:`i`\ th iteration, a numerical test is made to see if the current solution :math:`x` satisfies the constraints in the working set.
If the largest residual of the constraints in the working set is judged to be too large, the current working set is refactorized and the variables are recomputed to satisfy the constraints more accurately.
If :math:`i\leq 0`, the default value is used.
**'Cold Start'** : valueless
Default
This option specifies how the initial working set is chosen.
With a 'Cold Start', ``iqp_dense`` chooses the initial working set based on the values of the variables and constraints at the initial point.
Broadly speaking, the initial working set will include equality constraints and bounds or inequality constraints that violate or 'nearly' satisfy their bounds (to within 'Crash Tolerance').
With a 'Warm Start', you must provide a valid definition of every element of the array :math:`\mathrm{istate}` (see :ref:`Parameters <h02cb-py2-py-parameters>` for the definition of this array). ``iqp_dense`` will override your specification of :math:`\mathrm{istate}` if necessary, so that a poor choice of the working set will not cause a fatal error.
For instance, any elements of :math:`\mathrm{istate}` which are set to :math:`-2`, :math:`-1` or :math:`4` will be reset to zero, as will any elements which are set to :math:`3` when the corresponding elements of :math:`\mathrm{bl}` and :math:`\mathrm{bu}` are not equal.
A warm start will be advantageous if a good estimate of the initial working set is available -- for example, when ``iqp_dense`` is called repeatedly to solve related problems.
**'Warm Start'** : valueless
This option specifies how the initial working set is chosen.
With a 'Cold Start', ``iqp_dense`` chooses the initial working set based on the values of the variables and constraints at the initial point.
Broadly speaking, the initial working set will include equality constraints and bounds or inequality constraints that violate or 'nearly' satisfy their bounds (to within 'Crash Tolerance').
With a 'Warm Start', you must provide a valid definition of every element of the array :math:`\mathrm{istate}` (see :ref:`Parameters <h02cb-py2-py-parameters>` for the definition of this array). ``iqp_dense`` will override your specification of :math:`\mathrm{istate}` if necessary, so that a poor choice of the working set will not cause a fatal error.
For instance, any elements of :math:`\mathrm{istate}` which are set to :math:`-2`, :math:`-1` or :math:`4` will be reset to zero, as will any elements which are set to :math:`3` when the corresponding elements of :math:`\mathrm{bl}` and :math:`\mathrm{bu}` are not equal.
A warm start will be advantageous if a good estimate of the initial working set is available -- for example, when ``iqp_dense`` is called repeatedly to solve related problems.
**'Crash Tolerance'** : float
Default :math:`\text{} = 0.01`
This value is used in conjunction with the option 'Cold Start' (the default value) when ``iqp_dense`` selects an initial working set.
If :math:`0\leq r\leq 1`, the initial working set will include (if possible) bounds or general inequality constraints that lie within :math:`r` of their bounds.
In particular, a constraint of the form :math:`a_j^\mathrm{T}x\geq l` will be included in the initial working set if :math:`\left\lvert a_j^\mathrm{T}x-l\right\rvert \leq r\left(1+\left\lvert l\right\rvert \right)`.
If :math:`r < 0` or :math:`r > 1`, the default value is used.
**'Defaults'** : valueless
This special keyword may be used to reset all options to their default values.
**'Expand Frequency'** : int
Default :math:`\text{} = 5`
This option is part of an anti-cycling procedure designed to guarantee progress even on highly degenerate problems.
The strategy is to force a positive step at every iteration, at the expense of violating the constraints by a small amount.
Suppose that the value of the option 'Feasibility Tolerance' is :math:`\delta`.
Over a period of :math:`i` iterations, the feasibility tolerance actually used by ``iqp_dense`` (i.e., the `working` feasibility tolerance) increases from :math:`0.5\delta` to :math:`\delta` (in steps of :math:`0.5\delta /i`).
At certain stages the following 'resetting procedure' is used to remove constraint infeasibilities.
First, all variables whose upper or lower bounds are in the working set are moved exactly onto their bounds.
A count is kept of the number of nontrivial adjustments made.
If the count is positive, iterative refinement is used to give variables that satisfy the working set to (essentially) machine precision.
Finally, the working feasibility tolerance is reinitialized to :math:`0.5\delta`.
If a problem requires more than :math:`i` iterations, the resetting procedure is invoked and a new cycle of :math:`i` iterations is started with :math:`i` incremented by :math:`10`. (The decision to resume the feasibility phase or optimality phase is based on comparing any constraint infeasibilities with :math:`\delta`.)
The resetting procedure is also invoked when ``iqp_dense`` reaches an apparently optimal, infeasible or unbounded solution, unless this situation has already occurred twice.
If any nontrivial adjustments are made, iterations are continued.
If :math:`i\leq 0`, the default value is used.
If :math:`i\geq 9999999`, no anti-cycling procedure is invoked.
**'Feasibility Phase Iteration Limit'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m_L\right)}\right)`
The scalars :math:`i_1` and :math:`i_2` specify the maximum number of iterations allowed in the feasibility and optimality phases. 'Optimality Phase Iteration Limit' is equivalent to 'Iteration Limit'.
Setting :math:`i_1 = 0` and :math:`\text{‘Print Level'} > 0` means that the workspace needed will be computed and printed, but no iterations will be performed.
If :math:`i_1 < 0` or :math:`i_2 < 0`, the default value is used.
**'Optimality Phase Iteration Limit'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m_L\right)}\right)`
The scalars :math:`i_1` and :math:`i_2` specify the maximum number of iterations allowed in the feasibility and optimality phases. 'Optimality Phase Iteration Limit' is equivalent to 'Iteration Limit'.
Setting :math:`i_1 = 0` and :math:`\text{‘Print Level'} > 0` means that the workspace needed will be computed and printed, but no iterations will be performed.
If :math:`i_1 < 0` or :math:`i_2 < 0`, the default value is used.
**'Feasibility Tolerance'** : float
Default :math:`\text{} = \sqrt{\epsilon }`
If :math:`r\geq \epsilon`, :math:`r` defines the maximum acceptable `absolute` violation in each constraint at a 'feasible' point.
For example, if the variables and the coefficients in the general constraints are of order unity, and the latter are correct to about :math:`6` decimal digits, it would be appropriate to specify :math:`r` as :math:`10^{-6}`.
If :math:`0\leq r < \epsilon`, the default value is used.
``iqp_dense`` attempts to find a feasible solution before optimizing the objective function.
If the sum of infeasibilities cannot be reduced to zero, the 'Minimum Sum of Infeasibilities' can be used to find the minimum value of the sum.
Let Sinf be the corresponding sum of infeasibilities.
If Sinf is quite small, it may be appropriate to raise :math:`r` by a factor of :math:`10` or :math:`100`.
Otherwise, some error in the data should be suspected.
Note that a 'feasible solution' is a solution that satisfies the current constraints to within the tolerance :math:`r`.
**'Hessian Rows'** : int
Default :math:`\text{} = n`
Note that this option does not apply to problems of type FP or LP.
This specifies :math:`m`, the number of rows of the Hessian matrix :math:`H`.
The default value of :math:`m` is :math:`n`, the number of variables of the problem.
If the problem is of type QP, :math:`m` will usually be :math:`n`, the number of variables.
However, a value of :math:`m` less than :math:`n` is appropriate for QP3 or QP4 if :math:`\mathrm{h}` is an upper trapezoidal matrix with :math:`m` rows.
Similarly, :math:`m` may be used to define the dimension of a leading block of nonzeros in the Hessian matrices of QP1 or QP2, in which case the last :math:`n-m` rows and columns of :math:`\mathrm{h}` are assumed to be zero.
In the QP case, :math:`m` should not be greater than :math:`n`; if it is, the last :math:`m-n` rows of :math:`\mathrm{h}` are ignored.
If :math:`i < 0` or :math:`i > n`, the default value is used.
**'Infinite Bound Size'** : float
Default :math:`\text{} = 10^{20}`
If :math:`r > 0`, :math:`r` defines the 'infinite' bound :math:`\textit{bigbnd}` in the definition of the problem constraints.
Any upper bound greater than or equal to :math:`\textit{bigbnd}` will be regarded as :math:`+\infty` (and similarly any lower bound less than or equal to :math:`{-\textit{bigbnd}}` will be regarded as :math:`-\infty`).
If :math:`r\leq 0`, the default value is used.
**'Infinite Step Size'** : float
Default :math:`\text{} = \mathrm{max}\left(\textit{bigbnd}, 10^{20}\right)`
If :math:`r > 0`, :math:`r` specifies the magnitude of the change in variables that will be considered a step to an unbounded solution. (Note that an unbounded solution can occur only when the Hessian is not positive definite.) If the change in :math:`x` during an iteration would exceed the value of :math:`r`, the objective function is considered to be unbounded below in the feasible region.
If :math:`r\leq 0`, the default value is used.
**'Iteration Limit'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m_L\right)}\right)`
See option 'Feasibility Phase Iteration Limit'.
**'Iters'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m_L\right)}\right)`
See option 'Feasibility Phase Iteration Limit'.
**'Itns'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m_L\right)}\right)`
See option 'Feasibility Phase Iteration Limit'.
**'List'** : valueless
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'Nolist'** : valueless
Default
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'Maximum Degrees of Freedom'** : int
Default :math:`\text{} = n`
Note that this option does not apply to problems of type FP or LP.
This places a limit on the storage allocated for the triangular factor :math:`R` of the reduced Hessian :math:`H_R`.
Ideally, :math:`i` should be set slightly larger than the value of :math:`n_R` expected at the solution.
It need not be larger than :math:`m_n+1`, where :math:`m_n` is the number of variables that appear nonlinearly in the quadratic objective function.
For many problems it can be much smaller than :math:`m_n`.
For quadratic problems, a minimizer may lie on any number of constraints, so that :math:`n_R` may vary between :math:`1` and :math:`n`.
The default value of :math:`i` is, therefore, the number of variables :math:`n`.
If 'Hessian Rows' :math:`m` is specified, the default value of :math:`i` is the same number, :math:`m`.
**'Minimum Sum of Infeasibilities'** : str
Default :math:`= \texttt{'NO'}`
If no feasible point exists for the constraints, this option is used to control whether or not ``iqp_dense`` will calculate a point that minimizes the constraint violations.
If :math:`\text{‘Minimum Sum of Infeasibilities'} = \texttt{'[N]O'}`, ``iqp_dense`` will terminate as soon as it is evident that no feasible point exists for the constraints.
The final point will generally not be the point at which the sum of infeasibilities is minimized.
If :math:`\text{‘Minimum Sum of Infeasibilities'} = \texttt{'[Y]ES'}`, ``iqp_dense`` will continue until the sum of infeasibilities is minimized.
**'Monitoring File'** : int
Default :math:`\text{} = {-1}`
If :math:`i\geq 0` and :math:`\text{‘Print Level'} \geq 5`, monitoring information produced by ``iqp_dense`` at every iteration is sent to a file with logical unit number :math:`i`.
If :math:`i < 0` and/or :math:`\text{‘Print Level'} < 5`, no monitoring information is produced.
**'Optimality Tolerance'** : float
Default :math:`\text{} = \epsilon^{0.8}`
If :math:`r\geq \epsilon`, :math:`r` defines the tolerance used to determine if the bounds and general constraints have the right 'sign' for the solution to be judged to be optimal.
If :math:`0\leq r < \epsilon`, the default value is used.
**'Print Level'** : int
Default :math:`\text{} = 10`
The value of :math:`i` controls the amount of printout produced by ``iqp_dense``, as indicated below.
A detailed description of the printed output is given in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#fcomments2>`__ (summary output at each iteration and the final solution) and `Monitoring Information <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#monitoring>`__ (monitoring information at each iteration).
If :math:`i < 0`, the default value is used.
The following printout is sent to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`):
.. rst-class:: nag-rules-none nag-align-left
+----------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`i` |**Output** |
+======================+===============================================================================================================================================================================================================================+
|:math:`0` |No output. |
+----------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |The final solution only. |
+----------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`5` |One line of summary output (:math:`\text{} < 80` characters; see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#fcomments2>`__) for each iteration (no printout of the final solution).|
+----------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 10`|The final solution and one line of summary output for each iteration. |
+----------------------+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
The following printout is sent to the logical unit number defined by the 'Monitoring File':
.. rst-class:: nag-rules-none nag-align-left
+----------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`i` |**Output** |
+======================+====================================================================================================================================================================================================================================================================================================================================================================================================================================+
|:math:`\text{} < 5` |No output. |
+----------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 5` |One long line of output (:math:`\text{} > 80` characters; see `Monitoring Information <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#monitoring>`__) for each iteration (no printout of the final solution). |
+----------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 20`|At each iteration, the Lagrange-multipliers, the variables :math:`x`, the constraint values :math:`Ax` and the constraint status. |
+----------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 30`|At each iteration, the diagonal elements of the upper triangular matrix :math:`T` associated with the :math:`TQ` factorization `(3) <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#eqn3>`__ (see `Definition of the Search Direction <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#method2>`__) of the working set, and the diagonal elements of the upper triangular matrix :math:`R`.|
+----------------------+------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
If :math:`\text{‘Print Level'} \geq 5` and the unit number defined by 'Monitoring File' is the advisory unit number, then the summary output is suppressed.
**'Problem Type'** : str
Default :math:`=` QP2
This option specifies the type of objective function to be minimized during the optimality phase.
The following are the five optional keywords and the dimensions of the arrays that must be specified in order to define the objective function:
.. rst-class:: nag-rules-none nag-align-left
+---+--------------------------------------------------------------------------------------------------------------+
|LP |:math:`\mathrm{h}` not referenced, length-:math:`\textit{n}` :math:`\mathrm{cvec}` required; |
+---+--------------------------------------------------------------------------------------------------------------+
|QP1|rank-:math:`2` :math:`\mathrm{h}` symmetric, :math:`\mathrm{cvec}` not referenced; |
+---+--------------------------------------------------------------------------------------------------------------+
|QP2|rank-:math:`2` :math:`\mathrm{h}` symmetric, length-:math:`\textit{n}` :math:`\mathrm{cvec}` required; |
+---+--------------------------------------------------------------------------------------------------------------+
|QP3|rank-:math:`2` :math:`\mathrm{h}` upper trapezoidal, :math:`\mathrm{cvec}` not referenced; |
+---+--------------------------------------------------------------------------------------------------------------+
|QP4|rank-:math:`2` :math:`\mathrm{h}` upper trapezoidal, length-:math:`\textit{n}` :math:`\mathrm{cvec}` required.|
+---+--------------------------------------------------------------------------------------------------------------+
For problems of type FP, the objective function is omitted and neither :math:`\mathrm{h}` nor :math:`\mathrm{cvec}` are referenced.
The following keywords are also acceptable.
.. rst-class:: nag-rules-none nag-align-left
+---------+------+
|:math:`a`|Option|
+=========+======+
|Quadratic|QP2 |
+---------+------+
|Linear |LP |
+---------+------+
|Feasible |FP |
+---------+------+
In addition, the keyword QP is equivalent to the default option QP2.
If :math:`\mathrm{h} = 0`, i.e., the objective function is purely linear, the efficiency of ``iqp_dense`` may be increased by specifying :math:`a` as LP.
**'Rank Tolerance'** : float
Default :math:`\text{} = 100\epsilon`
Note that this option does not apply to problems of type FP or LP.
This argument enables you to control the condition number of the triangular factor :math:`R` (see `Algorithmic Details <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cbf.html#algdetails>`__).
If :math:`\rho_i` denotes the function :math:`\rho_i = \mathrm{max}\left\{\left\lvert R_{11}\right\rvert,\left\lvert R_{22}\right\rvert,\ldots,\left\lvert R_{{ii}}\right\rvert \right\}`, the dimension of :math:`R` is defined to be smallest index :math:`i` such that :math:`\left\lvert R_{{i+1,i+1}}\right\rvert \leq \sqrt{r}\left\lvert \rho_{{i+1}}\right\rvert`.
If :math:`r\leq 0`, the default value is used.
.. _h02cb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
Invalid input argument for basic problem.
(`errno` :math:`1`)
On entry, :math:`\mathrm{strtgy} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{strtgy} = 0`, :math:`1`, :math:`2` or :math:`3`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{mdepth} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mdepth}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{lintvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{lintvr}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{nclin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nclin}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
(`errno` :math:`2`)
No integer solutions found.
(`errno` :math:`3`)
Need to increase depth, :math:`\mathrm{mdepth} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
Basic problem (without integer constraints) unbounded.
(`errno` :math:`5`)
Basic problem is infeasible.
(`errno` :math:`6`)
Basic problem requires too many iterations.
(`errno` :math:`7`)
Basic problem has a reduced Hessian exceeding assigned dimensions.
(`errno` :math:`9`)
Error in designated problem type.
(`errno` :math:`12`)
Unexpected error in NAG function -- contact `NAG <https://www.nag.com>`__.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`-1`)
Halted at your request.
.. _h02cb-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_dense`` uses a 'Branch and Bound' algorithm in conjunction with :meth:`opt.qp_dense_solve <naginterfaces.library.opt.qp_dense_solve>` to try and determine integer solutions to a general quadratic programming problem.
The problem is assumed to be stated in the following general form:
.. math::
\mathrm{minimize}_{{x \in R^n}}f\left(x\right)\quad \text{ subject to }\quad l\leq \begin{pmatrix}x\\Ax\end{pmatrix}\leq u\text{,}
where :math:`A` is an :math:`m_L\times n` matrix and :math:`f\left(x\right)` may be specified in a variety of ways depending upon the particular problem to be solved.
The available forms for :math:`f\left(x\right)` are listed in Table [label omitted], in which the prefixes FP, LP and QP stand for 'feasible point', 'linear programming' and 'quadratic programming' respectively and :math:`c` is an :math:`n`-element vector.
+------------+-----------------------------------------------------------+-----------------------------------+
|Problem type|:math:`f\left(x\right)` |Matrix :math:`H` |
+============+===========================================================+===================================+
|FP |Not applicable |Not applicable |
+------------+-----------------------------------------------------------+-----------------------------------+
|LP |:math:`c^\mathrm{T}x` |Not applicable |
+------------+-----------------------------------------------------------+-----------------------------------+
|QP1 |:math:`\frac{1}{2}x^\mathrm{T}Hx` |symmetric |
+------------+-----------------------------------------------------------+-----------------------------------+
|QP2 |:math:`c^\mathrm{T}x+\frac{1}{2}x^\mathrm{T}Hx` |symmetric |
+------------+-----------------------------------------------------------+-----------------------------------+
|QP3 |:math:`\frac{1}{2}x^\mathrm{T}H^\mathrm{T}Hx` |:math:`m\times n` upper trapezoidal|
+------------+-----------------------------------------------------------+-----------------------------------+
|QP4 |:math:`c^\mathrm{T}x+\frac{1}{2}x^\mathrm{T}H^\mathrm{T}Hx`|:math:`m\times n` upper trapezoidal|
+------------+-----------------------------------------------------------+-----------------------------------+
Only when the problem is linear or the matrix :math:`H` is positive definite can the technique be guaranteed to work; but often useful results can be obtained for a wider class of problems.
The default problem type is QP2 and other objective functions are selected by using the option 'Problem Type'.
For problems of type FP, the objective function is omitted and ``iqp_dense`` attempts to find a feasible point for the set of constraints.
Branch and bound consists firstly of obtaining a solution without any of the variables :math:`x = \left(x_1,x_2,\ldots,x_n\right)^\mathrm{T}` constrained to be integer.
Suppose :math:`x_1` ought to be integer, but at the optimal value just computed :math:`x_1 = 2.4`.
A constraint :math:`x_1\leq 2` is added to the system and the second problem solved.
A constraint :math:`x_1\geq 3` gives rise to a third sub-problem.
In a similar manner a whole series of sub-problems may be generated, corresponding to integer constraints on the variables.
The sub-problems are all solved using :meth:`opt.qp_dense_solve <naginterfaces.library.opt.qp_dense_solve>`.
In practice the function tries to compute an integer solution as quickly as possible using a depth-first approach, since this helps determine a realistic cut-off value.
If we have a cut-off value, say the value of the function at this first integer solution, and any sub-problem, :math:`W` say, has a solution value greater than this cut-off value, then subsequent sub-problems of :math:`W` must have solutions greater than the value of the solution at :math:`W` and, therefore, need not be computed.
Thus a knowledge of a good cut-off value can result in fewer sub-problems being solved and thus speed up the operation of the function. (See the description of :math:`\mathrm{monit}` in :ref:`Parameters <h02cb-py2-py-parameters>` for details of how you can supply your own cut-off value.)
.. _h02cb-py2-py-references:
**References**
Gill, P E, Hammarling, S, Murray, W, Saunders, M A and Wright, M H, 1986, `Users' guide for LSSOL (Version 1.0)`, Report SOL 86-1, Department of Operations Research, Stanford University
Gill, P E and Murray, W, 1978, `Numerically stable methods for quadratic programming`, Math. Programming (14), 349--372
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1984, `Procedures for optimization problems with a mixture of bounds and general linear constraints`, ACM Trans. Math. Software (10), 282--298
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1989, `A practical anti-cycling procedure for linearly constrained optimization`, Math. Programming (45), 437--474
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1991, `Inertia-controlling methods for general quadratic programming`, SIAM Rev. (33), 1--36
Gill, P E, Murray, W and Wright, M H, 1981, `Practical Optimization`, Academic Press
Pardalos, P M and Schnitger, G, 1988, `Checking local optimality in constrained quadratic programming is NP-hard`, Operations Research Letters (7), 33--35
"""
raise NotImplementedError
[docs]def iqp_dense_optfile(ioptns, comm, io_manager=None):
r"""
``iqp_dense_optfile`` may be used to supply options to :meth:`iqp_dense` from an external file.
.. _h02cc-py2-py-doc:
For full information please refer to the NAG Library document for h02cc
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02ccf.html
.. _h02cc-py2-py-parameters:
**Parameters**
**ioptns** : int
The unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) of the options file to be read.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qp_dense_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
.. _h02cc-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{ioptns} = \langle\mathit{\boldsymbol{value}}\rangle`. Constraint: :math:`0\leq \mathrm{ioptns}\leq 2147483647`.
(`errno` :math:`2`)
``Begin`` was found, but end-of-file was found before ``End`` was found.
(`errno` :math:`3`)
End-of-file was found before ``Begin`` was found.
(`errno` :math:`5`)
One or more lines of the options file is invalid.
.. _h02cc-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_dense_optfile`` may be used to supply values for options to :meth:`iqp_dense`. ``iqp_dense_optfile`` reads an external file and each line of the file defines a single option.
It is only necessary to supply values for those arguments whose values are to be different from their default values.
Each option is defined by a single character string of up to :math:`72` characters, consisting of one or more items.
The items associated with a given option must be separated by spaces, or equal signs :math:`\left[ = \right]`.
Alphabetic characters may be upper or lower case.
The string
::
Print level = 1
is an example of a string used to set an option.
For each option the string contains one or more of the following items:
- a mandatory keyword;
- a phrase that qualifies the keyword;
- a number that specifies an `int` or `float` value. Such numbers may be up to :math:`16` contiguous characters in Fortran 77's I, F, E or D formats, terminated by a space if this is not the last item on the line.
Blank strings and comments are ignored.
A comment begins with an asterisk (*) and all subsequent characters in the string are regarded as part of the comment.
The file containing the options must start with ``Begin`` and must finish with ``End``.
An example of a valid options file is:
::
Begin * Example options file
Print level = 10
End
Printing of user-supplied options is turned off by default, but may be turned on at any time using the keyword 'List'.
Option settings are preserved following a call to :meth:`iqp_dense`, and so the keyword 'Defaults' is provided to allow you to reset all the options to their default values prior to a subsequent call to :meth:`iqp_dense`.
A complete list of options, their abbreviations, synonyms and default values is given in :ref:`Other Parameters for iqp_dense <h02cb-py2-py-other_params>`.
"""
raise NotImplementedError
[docs]def iqp_dense_optstr(optstr, comm, io_manager=None):
r"""
``iqp_dense_optstr`` may be used to supply individual options to :meth:`iqp_dense`.
.. _h02cd-py2-py-doc:
For full information please refer to the NAG Library document for h02cd
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cdf.html
.. _h02cd-py2-py-parameters:
**Parameters**
**optstr** : str
A single valid option string (as described in :ref:`Notes <h02cd-py2-py-notes>` above and in :ref:`Other Parameters for iqp_dense <h02cb-py2-py-other_params>`).
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qp_dense_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
.. _h02cd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`5`)
The supplied option string is invalid. Supplied value was: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _h02cd-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_dense_optstr`` may be used to supply values for options to :meth:`iqp_dense`.
It is only necessary to call ``iqp_dense_optstr`` for those arguments whose values are to be different from their default values.
One call to ``iqp_dense_optstr`` sets one argument value.
Each option is defined by a single character string of up to :math:`72` characters, consisting of one or more items.
The items associated with a given option must be separated by spaces, or equal signs :math:`\left[ = \right]`.
Alphabetic characters may be upper or lower case.
The string
::
Print level = 1
is an example of a string used to set an option.
For each option the string contains one or more of the following items:
- a mandatory keyword;
- a phrase that qualifies the keyword;
- a number that specifies an `int` or `float` value. Such numbers may be up to :math:`16` contiguous characters in Fortran 77's I, F, E or D formats, terminated by a space if this is not the last item on the line.
Blank strings and comments are ignored.
A comment begins with an asterisk (*) and all subsequent characters in the string are regarded as part of the comment.
Printing of user-specified options is turned off by default.
It may be turned on at any time using the keyword 'List'.
Option settings are preserved following a call to :meth:`iqp_dense`, and so the keyword 'Defaults' is provided to allow you to reset all the options to their default values.
A complete list of options, their abbreviations, synonyms and default values is given in :ref:`Other Parameters for iqp_dense <h02cb-py2-py-other_params>`.
"""
raise NotImplementedError
[docs]def iqp_sparse(m, iobj, ncolh, a, ha, ka, bl, bu, start, names, crname, ns, xs, intvar, istate, strtgy, leniz, lenz, comm, qphx=None, mdepth=None, monit=None, data=None, io_manager=None):
r"""
``iqp_sparse`` obtains integer solutions to sparse linear programming and quadratic programming problems.
Note: this function uses optional algorithmic parameters, see also: :meth:`iqp_sparse_optfile`, :meth:`iqp_sparse_optstr`.
.. _h02ce-py2-py-doc:
For full information please refer to the NAG Library document for h02ce
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html
.. _h02ce-py2-py-parameters:
**Parameters**
**m** : int
:math:`m`, the number of general linear constraints (or slacks). This is the number of rows in :math:`A`, including the free row (if any; see :math:`\mathrm{iobj}`).
**iobj** : int
If :math:`\mathrm{iobj} > 0`, row :math:`\mathrm{iobj}` of :math:`A` is a free row containing the nonzero elements of the vector :math:`c` appearing in the linear objective term :math:`c^{\mathrm{T}}x`.
If :math:`\mathrm{iobj} = 0`, there is no free row, i.e., the problem is either an FP problem (in which case :math:`\mathrm{iobj}` must be set to zero), or a QP problem with :math:`c = 0`.
**ncolh** : int
:math:`n_H`, the number of leading nonzero columns of the Hessian matrix :math:`H`. For FP and LP problems, :math:`\mathrm{ncolh}` must be set to zero.
**a** : float, array-like, shape :math:`\left(\textit{nnz}\right)`
`On entry`: the nonzero elements of :math:`A`, ordered by increasing column index. Note that multiple elements with the same row and column indices are not allowed.
`On exit`: used as internal workspace prior to being restored and hence is unchanged.
**ha** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
:math:`\mathrm{ha}[\textit{i}]` must contain the row index of the nonzero element stored in :math:`\mathrm{a}[\textit{i}]`, for :math:`\textit{i} = 1,2,\ldots,\textit{nnz}`. Note that the row indices for a column may be supplied in any order.
**ka** : int, array-like, shape :math:`\left(n+1\right)`
:math:`\mathrm{ka}[\textit{j}]` must contain the index in :math:`\mathrm{a}` of the start of the :math:`\textit{j}`\ th column, for :math:`\textit{j} = 1,2,\ldots,n`. To specify the :math:`j`\ th column as empty, set :math:`\mathrm{ka}[j] = \mathrm{ka}[j+1]`. Note that the first and last elements of :math:`\mathrm{ka}` must be such that :math:`\mathrm{ka}[0] = 1` and :math:`\mathrm{ka}[n] = \textit{nnz}+1`.
**bl** : float, array-like, shape :math:`\left(n+\mathrm{m}\right)`
:math:`l`, the lower bounds for all the variables and general constraints, in the following order. The first :math:`\textit{n}` elements of :math:`\mathrm{bl}` must contain the bounds on the variables :math:`x`, and the next :math:`\mathrm{m}` elements the bounds for the general linear constraints :math:`Ax` (or slacks :math:`s`) and the free row (if any). To specify a nonexistent lower bound (i.e., :math:`l_j = {-\infty }`), set :math:`\mathrm{bl}[j]\leq -\textit{bigbnd}`, where :math:`\textit{bigbnd}` is the value of the option 'Infinite Bound Size' (:math:`\text{default value} = 10^{20}`). To specify the :math:`j`\ th constraint as an `equality`, set :math:`\mathrm{bl}[j] = \mathrm{bu}[j] = \beta`, say, where :math:`\left\lvert \beta \right\rvert < \textit{bigbnd}`. Note that the lower bound corresponding to the free row must be set to :math:`{-\infty }` and stored in :math:`\mathrm{bl}[n+\mathrm{iobj}-1]`.
**bu** : float, array-like, shape :math:`\left(n+\mathrm{m}\right)`
`On entry`: :math:`u`, the upper bounds for all the variables and general constraints, in the following order. The first :math:`\textit{n}` elements of :math:`\mathrm{bl}` must contain the bounds on the variables :math:`x`, and the next :math:`\mathrm{m}` elements the bounds for the general linear constraints :math:`Ax` (or slacks :math:`s`) and the free row (if any). To specify a nonexistent upper bound (i.e., :math:`u_j = {+\infty }`), set :math:`\mathrm{bu}[j]\geq \textit{bigbnd}`. Note that the upper bound corresponding to the free row must be set to :math:`{+\infty }` and stored in :math:`\mathrm{bu}[n+\mathrm{iobj}-1]`.
`On exit`: used as internal workspace prior to being restored and hence is unchanged.
**start** : str, length 1
Indicates how a starting basis is to be obtained.
:math:`\mathrm{start} = \texttt{'C'}`
An internal crash procedure will be used to choose an initial basis matrix :math:`B`.
:math:`\mathrm{start} = \texttt{'W'}`
A basis is already defined in :math:`\mathrm{istate}` (probably from a previous call).
**names** : str, length 8, array-like, shape :math:`\left(5\right)`
A set of names associated with the so-called MPSX form of the problem.
:math:`\mathrm{names}[0]`
Must contain the name for the problem (or be blank).
:math:`\mathrm{names}[1]`
Must contain the name for the free row (or be blank).
:math:`\mathrm{names}[2]`
Must contain the name for the constraint right-hand side (or be blank).
:math:`\mathrm{names}[3]`
Must contain the name for the ranges (or be blank).
:math:`\mathrm{names}[4]`
Must contain the name for the bounds (or be blank).
(These names are used in the monitoring file output; see `Monitoring Information <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#monitoring>`__.)
**crname** : str, length 8, array-like, shape :math:`\left(\textit{nname}\right)`
The optional column and row names.
If :math:`\textit{nname} = 1`, :math:`\mathrm{crname}` is not referenced and the printed output will use default names for the columns and rows.
If :math:`\textit{nname} = n+\mathrm{m}`, the first :math:`\textit{n}` elements must contain the names for the columns and the next :math:`\mathrm{m}` elements must contain the names for the rows.
Note that the name for the free row (if any) must be stored in :math:`\mathrm{crname}[n+\mathrm{iobj}]`.
**ns** : int
:math:`n_S`, the number of superbasics. For QP problems, :math:`\mathrm{ns}` need not be specified if :math:`\mathrm{start} = \texttt{'C'}`, but must retain its value from a previous call when :math:`\mathrm{start} = \texttt{'W'}`. For FP and LP problems, :math:`\mathrm{ns}` need not be initialized.
**xs** : float, array-like, shape :math:`\left(n+\mathrm{m}\right)`
The initial values of the variables and slacks :math:`\left(x, s\right)`. (See the description for :math:`\mathrm{istate}`.)
**intvar** : int, array-like, shape :math:`\left(\textit{lintvr}\right)`
Specifies which components of the solution vector :math:`x` are constrained to be integer. Specifically, if :math:`k` elements of :math:`x` are required to take integer values then :math:`\mathrm{intvar}[\textit{i}] = l_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,k`, where :math:`l_i` is the integer index such that :math:`x_{l_i}` is integer. If :math:`k < \textit{lintvr}` then :math:`\mathrm{intvar}[k+1]` must be set to :math:`-1` to signal the end of the integer variable indices.
The order in which the indices of those components of :math:`x` required to be integer is presented determines the order in which the sub-problems are treated and solved.
As such it can be a powerful tool to assist the function in achieving a solution efficiently.
The general advice is to enter the important integer variables in the model early in :math:`\mathrm{intvar}`; secondary or less important variables should be entered near the end of the list.
However some experimentation might be required to find the optimal order.
**istate** : int, array-like, shape :math:`\left(n+\mathrm{m}\right)`
If :math:`\mathrm{start} = \texttt{'C'}`, the first :math:`\textit{n}` elements of :math:`\mathrm{istate}` and :math:`\mathrm{xs}` must specify the initial states and values, respectively, of the variables :math:`x`. (The slacks :math:`s` need not be initialized.) An internal crash procedure is then used to select an initial basis matrix :math:`B`. The initial basis matrix will be triangular (neglecting certain small elements in each column). It is chosen from various rows and columns of columns of :math:`\left(A-I\right)`. Possible values for :math:`\mathrm{istate}[j]` are as follows:
.. rst-class:: nag-rules-none nag-align-left
+--------------------------+---------------------------------------------------------------------+
|:math:`\mathrm{istate}[j]`|State of :math:`\mathrm{xs}[j]` during crash procedure |
+==========================+=====================================================================+
|0 or :math:`1` |Eligible for the basis |
+--------------------------+---------------------------------------------------------------------+
|2 |Ignored |
+--------------------------+---------------------------------------------------------------------+
|3 |Eligible for the basis (given preference over :math:`0` or :math:`1`)|
+--------------------------+---------------------------------------------------------------------+
|4 or :math:`5` |Ignored |
+--------------------------+---------------------------------------------------------------------+
If nothing special is known about the problem, or there is no wish to provide special information, you may set :math:`\mathrm{istate}[\textit{j}] = 0` and :math:`\mathrm{xs}[\textit{j}] = 0.0`, for :math:`\textit{j} = 1,2,\ldots,n`. All variables will then be eligible for the initial basis. Less trivially, to say that the :math:`j`\ th variable will probably be equal to one of its bounds, set :math:`\mathrm{istate}[j] = 4` and :math:`\mathrm{xs}[j] = \mathrm{bl}[j]` or :math:`\mathrm{istate}[j] = 5` and :math:`\mathrm{xs}[j] = \mathrm{bu}[j]` as appropriate.
Following the crash procedure, variables for which :math:`\mathrm{istate}[j] = 2` are made superbasic.
Other variables not selected for the basis are then made nonbasic at the value :math:`\mathrm{xs}[j]` if :math:`\mathrm{bl}[j]\leq \mathrm{xs}[j]\leq \mathrm{bu}[j]`, or at the value :math:`\mathrm{bl}[j]` or :math:`\mathrm{bu}[j]` closest to :math:`\mathrm{xs}[j]`.
If :math:`\mathrm{start} = \texttt{'W'}`, :math:`\mathrm{istate}` and :math:`\mathrm{xs}` must specify the initial states and values, respectively, of the variables and slacks :math:`\left(x, s\right)`.
If ``iqp_sparse`` has been called previously with the same values of :math:`\textit{n}` and :math:`\mathrm{m}`, :math:`\mathrm{istate}` already contains satisfactory information.
**strtgy** : int
Defines the branching strategy adopted by the function.
:math:`\mathrm{strtgy} = 0`
Each sub-problem first explored imposes a tighter upper bound on the component of :math:`x`.
:math:`\mathrm{strtgy} = 1`
Each sub-problem first explored imposes a tighter lower bound on the component of :math:`x`.
:math:`\mathrm{strtgy} = 2`
Each branch explored imposes a tighter upper bound on the component of :math:`x` if its fractional part is less than :math:`0.5`, otherwise it imposes a tighter lower bound.
:math:`\mathrm{strtgy} = 3`
Random choice is made between first exploring a tighter lower bound or a tighter upper bound sub-problem.
**leniz** : int
The dimension of the internal workspace array :math:`\textit{iz}`.
**lenz** : int
The dimension of the internal workspace array :math:`\textit{z}`.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qpconvex1_sparse_solve'}`.
**qphx** : None or callable hx = qphx(nstate, x, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
For QP problems, you must supply a version of :math:`\mathrm{qphx}` to compute the matrix product :math:`Hx`.
If :math:`H` has rows and columns consisting entirely of zeros, it is most efficient to order the variables :math:`x = \left(y\quad \text{ }\quad z\right)^\mathrm{T}` so that
.. math::
Hx = \begin{pmatrix}H_1&0\\0&0\end{pmatrix}\begin{pmatrix}y\\z\end{pmatrix} = \begin{pmatrix}H_1y\\0\end{pmatrix}\text{,}
where the nonlinear variables :math:`y` appear first as shown.
For LP problems, :math:`\mathrm{qphx}` will never be called by ``iqp_sparse``.
**Parameters**
**nstate** : int
If :math:`\mathrm{nstate} = 1`, then ``iqp_sparse`` is calling :math:`\mathrm{qphx}` for the first time on a sub-problem. This argument setting allows you to save computation time if certain data must be read or calculated only once.
If :math:`\mathrm{nstate}\geq 2`, then ``iqp_sparse`` is calling :math:`\mathrm{qphx}` for the last time.
This argument setting allows you to perform some additional computation on the final sub-problem solution.
In general, the last call to :math:`\mathrm{qphx}` is made with :math:`\mathrm{nstate} = 2+\textit{errno}` (see :ref:`Exceptions <h02ce-py2-py-errors>`).
Otherwise, :math:`\mathrm{nstate} = 0`.
**x** : float, ndarray, shape :math:`\left(\textit{ncolh}\right)`
The first :math:`\mathrm{ncolh}` elements of the vector :math:`x`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**hx** : float, array-like, shape :math:`\left(\textit{ncolh}\right)`
The product :math:`Hx`.
**mdepth** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:`2\times n+20`.
Specifies the maximum depth the tree of sub-problems may be developed.
**monit** : None or callable (bstval, halt) = monit(intfnd, nodes, depth, obj, x, bstval, bstsol, bl, bu, halt, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
To provide feed-back on the progress of the branch and bound process.
Additionally :math:`\mathrm{monit}` provides, via its argument :math:`\mathrm{halt}`, the ability to terminate the process. (You might choose to do this when an integer solution is found, rather than search for a better solution.) If you do not require any intermediate output then :math:`\mathrm{monit}` may be **None**.
**Parameters**
**intfnd** : int
Contains the number of integer solutions obtained so far.
**nodes** : int
Contains the number of nodes (sub-problems) solved so far.
**depth** : int
Contains the depth reached in the tree of problems.
**obj** : float
Contains the solution value to the sub-problem at this node.
**x** : float, ndarray, shape :math:`\left(n\right)`
Contains the solution vector to the sub-problem at this node.
**bstval** : float
Contains the value of the objective function corresponding to the best integer solution obtained so far. If no integer solution has been found :math:`\mathrm{bstval}` contains the largest machine representable number (see :meth:`machine.real_largest <naginterfaces.library.machine.real_largest>`).
**bstsol** : float, ndarray, shape :math:`\left(n\right)`
Contains the value of the best integer solution obtained so far.
**bl** : float, ndarray, shape :math:`\left(n\right)`
Contains the current lower bounds on the variables at this point.
**bu** : float, ndarray, shape :math:`\left(n\right)`
Contains the current upper bounds on the variables at this point.
**halt** : bool
Will have the value :math:`\mathbf{False}`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**bstval** : float
May be set to a cut-off value, if you are a sophisticated user, as follows. Before an integer solution has been found :math:`\mathrm{bstval}` will be set by ``iqp_sparse`` to the largest machine representable number (see :meth:`machine.real_largest <naginterfaces.library.machine.real_largest>`). If you know that the solution being sought is a much smaller number, then :math:`\mathrm{bstval}` may be set to this number as a cut-off value (see :ref:`Notes <h02ce-py2-py-notes>`). Beware of setting :math:`\mathrm{bstval}` too small, since then no integer solutions will be discovered. Also, on entry to :math:`\mathrm{monit}`, :math:`\mathrm{bstval}` should be set using a statement of the form
::
if intfnd == 0: bstval = cut_off_value
Setting :math:`\mathrm{bstval}` in this way will not prevent the normal operation of the algorithm when subsequent integer solutions are found. It would be a grievous mistake to unconditionally set :math:`\mathrm{bstval}` and if you have any doubts whatsoever about the correct use of this argument then you are strongly recommended to leave it unchanged.
**halt** : bool
If :math:`\mathrm{halt}` is set to :math:`\mathbf{True}`, :meth:`opt.qpconvex1_sparse_solve <naginterfaces.library.opt.qpconvex1_sparse_solve>` will be brought to a halt with :math:`\textit{errno}` exit :math:`-1`. This facility may be useful if you are content with `any` integer solution, or with any integer solution that fits certain criteria. Under these circumstances setting :math:`\mathrm{halt} = \mathbf{True}` can save considerable unnecessary computation.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**ns** : int
The final number of superbasics. This will be zero for FP and LP problems.
**xs** : float, ndarray, shape :math:`\left(n+\mathrm{m}\right)`
:math:`\mathrm{xs}[\textit{i}-1]` contains the final value of :math:`x_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`.
**istate** : int, ndarray, shape :math:`\left(n+\mathrm{m}\right)`
The final states of the variables and slacks :math:`\left(x, s\right)` from the solution of the last sub-problem tackled. The significance of each possible value of :math:`\mathrm{istate}[j]` is as follows:
.. rst-class:: nag-rules-none nag-align-left
+--------------------------+---------------------------+---------------------------------------------------------+
|:math:`\mathrm{istate}[j]`|State of variable :math:`j`|Normal value of :math:`\mathrm{xs}[j]` |
+==========================+===========================+=========================================================+
|:math:`0` |Nonbasic |:math:`\mathrm{bl}[j]` |
+--------------------------+---------------------------+---------------------------------------------------------+
|:math:`1` |Nonbasic |:math:`\mathrm{bu}[j]` |
+--------------------------+---------------------------+---------------------------------------------------------+
|:math:`2` |Superbasic |Between :math:`\mathrm{bl}[j]` and :math:`\mathrm{bu}[j]`|
+--------------------------+---------------------------+---------------------------------------------------------+
|:math:`3` |Basic |Between :math:`\mathrm{bl}[j]` and :math:`\mathrm{bu}[j]`|
+--------------------------+---------------------------+---------------------------------------------------------+
If :math:`\texttt{Ninf} = 0` (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#fc-printedoutput>`__), basic and superbasic variables may be outside their bounds by as much as the value of the option 'Feasibility Tolerance' (:math:`\text{default value} = \mathrm{max}\left(10^{-6}, {\sqrt{\epsilon }}\right)`, where :math:`\epsilon` is the machine precision). Note that unless the option :math:`\text{‘Scale Option'} = 0` (:math:`\text{default value} = 2`) is specified, the 'Feasibility Tolerance' applies to the variables of the scaled problem. In this case, the variables of the original problem may be as much as :math:`0.1` outside their bounds, but this is unlikely unless the problem is very badly scaled.
Very occasionally some nonbasic variables may be outside their bounds by as much as the 'Feasibility Tolerance', and there may be some nonbasic variables for which :math:`\mathrm{xs}[j]` lies strictly between its bounds.
If :math:`\texttt{Ninf} > 0`, some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by Sinf (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#fc-printedoutput>`__) if :math:`\text{‘Scale Option'} = 0`).
**miniz** : int
The minimum value of :math:`\mathrm{leniz}` required to start solving the problem. If :math:`\mathrm{errno}` = 14, ``iqp_sparse`` may be called again with :math:`\mathrm{leniz}` suitably larger than :math:`\mathrm{miniz}`. (The bigger the better, since it is not certain how much workspace the basis factors need.)
**minz** : int
The minimum value of :math:`\mathrm{lenz}` required to start solving the problem. If :math:`\mathrm{errno}` = 15, ``iqp_sparse`` may be called again with :math:`\mathrm{lenz}` suitably larger than :math:`\mathrm{minz}`. (The bigger the better, since it is not certain how much workspace the basis factors need.)
**obj** : float
The value of the objective function.
If :math:`\texttt{Ninf} = 0`, :math:`\mathrm{obj}` includes the quadratic objective term :math:`\frac{1}{2}x^{\mathrm{T}}Hx` (if any).
If :math:`\texttt{Ninf} > 0`, :math:`\mathrm{obj}` is just the linear objective term :math:`c^{\mathrm{T}}x` (if any).
For FP problems, :math:`\mathrm{obj}` is set to zero.
**clamda** : float, ndarray, shape :math:`\left(n+\mathrm{m}\right)`
A set of Lagrange-multipliers for the bounds on the variables and the general constraints. More precisely, the first :math:`\textit{n}` elements contain the multipliers (`reduced costs`) for the bounds on the variables, and the next :math:`\mathrm{m}` elements contain the multipliers (`shadow prices`) for the general linear constraints.
.. _h02ce-py2-py-other_params:
**Other Parameters**
**'Check Frequency'** : int
Default :math:`\text{} = 60`
Every :math:`i`\ th iteration after the most recent basis factorization, a numerical test is made to see if the current solution :math:`\left(x, s\right)` satisfies the linear constraints :math:`Ax-s = 0`.
If the largest element of the residual vector :math:`r = Ax-s` is judged to be too large, the current basis is refactorized and the basic variables recomputed to satisfy the constraints more accurately.
If :math:`i < 0`, the default value is used.
If :math:`i = 0`, the value :math:`i = 99999999` is used and effectively no checks are made.
**'Crash Option'** : int
Default :math:`\text{} = 2`
Note that this option does not apply when :math:`\mathrm{start} = \texttt{'W'}` (see :ref:`Parameters <h02ce-py2-py-parameters>`).
If :math:`\mathrm{start} = \texttt{'C'}`, an internal crash procedure is used to select an initial basis from various rows and columns of the constraint matrix :math:`\left(A-I\right)`.
The value of :math:`i` determines which rows and columns are initially eligible for the basis, and how many times the crash procedure is called.
If :math:`i = 0`, the all-slack basis :math:`B = -I` is chosen.
If :math:`i = 1`, the crash procedure is called once (looking for a triangular basis in all rows and columns of the linear constraint matrix :math:`A`).
If :math:`i = 2`, the crash procedure is called twice (looking at any `equality` constraints first followed by any `inequality` constraints).
If :math:`i < 0` or :math:`i > 2`, the default value is used.
If :math:`i = 1` or :math:`2`, certain slacks on inequality rows are selected for the basis first. (If :math:`i = 2`, numerical values are used to exclude slacks that are close to a bound.) The crash procedure then makes several passes through the columns of :math:`A`, searching for a basis matrix that is essentially triangular.
A column is assigned to 'pivot' on a particular row if the column contains a suitably large element in a row that has not yet been assigned. (The pivot elements ultimately form the diagonals of the triangular basis.) For remaining unassigned rows, slack variables are inserted to complete the basis.
**'Crash Tolerance'** : float
Default :math:`\text{} = 0.1`
This value allows the crash procedure to ignore certain 'small' nonzero elements in the constraint matrix :math:`A` while searching for a triangular basis.
For each column of :math:`A`, if :math:`a_{\mathrm{max}}` is the largest element in the column, other nonzeros in that column are ignored if they are less than (or equal to) :math:`a_{\mathrm{max}}\times r`.
When :math:`r > 0`, the basis obtained by the crash procedure may not be strictly triangular, but it is likely to be nonsingular and almost triangular.
The intention is to obtain a starting basis with more column variables and fewer (arbitrary) slacks.
A feasible solution may be reached earlier for some problems.
If :math:`r < 0` or :math:`r\geq 1`, the default value is used.
**'Defaults'** : valueless
This special keyword may be used to reset all options to their default values.
**'Expand Frequency'** : int
Default :math:`\text{} = 10000`
This option is part of an anti-cycling procedure (see `Miscellaneous <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#ad-misc>`__) designed to allow progress even on highly degenerate problems.
For LP problems, the strategy is to force a positive step at every iteration, at the expense of violating the constraints by a small amount.
Suppose that the value of the option 'Feasibility Tolerance' is :math:`\delta`.
Over a period of :math:`i` iterations, the feasibility tolerance actually used by ``iqp_sparse`` (i.e., the `working` feasibility tolerance) increases from :math:`0.5\delta` to :math:`\delta` (in steps of :math:`0.5\delta /i`).
For QP problems, the same procedure is used for iterations in which there is only one superbasic variable. (Cycling can only occur when the current solution is at a vertex of the feasible region.) Thus, zero steps are allowed if there is more than one superbasic variable, but otherwise positive steps are enforced.
Increasing the value of :math:`i` helps reduce the number of slightly infeasible nonbasic basic variables (most of which are eliminated during the resetting procedure).
However, it also diminishes the freedom to choose a large pivot element (see the description of the option 'Pivot Tolerance').
If :math:`i < 0`, the default value is used.
If :math:`i = 0`, the value :math:`i = 99999999` is used and effectively no anti-cycling procedure is invoked.
**'Factorization Frequency'** : int
Default :math:`\text{} = 100`
If :math:`i > 0`, at most :math:`i` basis changes will occur between factorizations of the basis matrix.
For LP problems, the basis factors are usually updated at every iteration.
For QP problems, fewer basis updates will occur as the solution is approached.
The number of iterations between basis factorizations will, therefore, increase.
During these iterations a test is made regularly according to the value of option 'Check Frequency' to ensure that the linear constraints :math:`Ax-s = 0` are satisfied.
If necessary, the basis will be refactorized before the limit of :math:`i` updates is reached.
If :math:`i\leq 0`, the default value is used.
**'Feasibility Tolerance'** : float
Default :math:`\text{} = \mathrm{max}\left(10^{{-6}}, {\sqrt{\epsilon }}\right)`
If :math:`r\geq \epsilon`, :math:`r` defines the maximum acceptable `absolute` violation in each constraint at a 'feasible' point (including slack variables).
For example, if the variables and the coefficients in the linear constraints are of order unity, and the latter are correct to about five decimal digits, it would be appropriate to specify :math:`r` as :math:`10^{-5}`.
If :math:`r < \epsilon`, the default value is used.
``iqp_sparse`` attempts to find a feasible solution before optimizing the objective function.
If the sum of infeasibilities cannot be reduced to zero, the problem is assumed to be `infeasible`.
Let Sinf be the corresponding sum of infeasibilities.
If Sinf is quite small, it may be appropriate to raise :math:`r` by a factor of :math:`10` or :math:`100`.
Otherwise, some error in the data should be suspected.
Note that the function does `not` attempt to find the minimum value of Sinf.
If the constraints and variables have been scaled (see the description of the option 'Scale Option'), then feasibility is defined in terms of the scaled problem (since it is more likely to be meaningful).
**'Infinite Bound Size'** : float
Default :math:`\text{} = 10^{20}`
If :math:`r > 0`, :math:`r` defines the 'infinite' bound :math:`\textit{bigbnd}` in the definition of the problem constraints.
Any upper bound greater than or equal to :math:`\textit{bigbnd}` will be regarded as :math:`+\infty` (and similarly any lower bound less than or equal to :math:`{-\textit{bigbnd}}` will be regarded as :math:`-\infty`).
If :math:`r\leq 0`, the default value is used.
**'Infinite Step Size'** : float
Default :math:`\text{} = \mathrm{max}\left(\textit{bigbnd}, {10^{20}}\right)`
If :math:`r > 0`, :math:`r` specifies the magnitude of the change in variables that will be considered a step to an unbounded solution. (Note that an unbounded solution can occur only when the Hessian is not positive definite.) If the change in :math:`x` during an iteration would exceed the value of :math:`r`, the objective function is considered to be unbounded below in the feasible region.
If :math:`r\leq 0`, the default value is used.
**'Iteration Limit'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m\right)}\right)`
The value of :math:`i` specifies the maximum number of iterations allowed before termination.
Setting :math:`i = 0` and :math:`\text{‘Print Level'} > 0` means that the workspace needed to start solving the problem will be computed and printed, but no iterations will be performed.
If :math:`i < 0`, the default value is used.
**'Iters'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m\right)}\right)`
The value of :math:`i` specifies the maximum number of iterations allowed before termination.
Setting :math:`i = 0` and :math:`\text{‘Print Level'} > 0` means that the workspace needed to start solving the problem will be computed and printed, but no iterations will be performed.
If :math:`i < 0`, the default value is used.
**'Itns'** : int
Default :math:`\text{} = \mathrm{max}\left(50, {5\left(n+m\right)}\right)`
The value of :math:`i` specifies the maximum number of iterations allowed before termination.
Setting :math:`i = 0` and :math:`\text{‘Print Level'} > 0` means that the workspace needed to start solving the problem will be computed and printed, but no iterations will be performed.
If :math:`i < 0`, the default value is used.
**'List'** : valueless
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'Nolist'** : valueless
Default
Option 'List' enables printing of each option specification as it is supplied. 'Nolist' suppresses this printing.
**'LU Factor Tolerance'** : float
Default :math:`\text{} = 100.0`
The values of :math:`r_1` and :math:`r_2` affect the stability and sparsity of the basis factorization :math:`B = LU`, during refactorization and updates respectively.
The lower triangular matrix :math:`L` is a product of matrices of the form
.. math::
\begin{pmatrix}1& \\\mu &1\end{pmatrix}
where the multipliers :math:`\mu` will satisfy :math:`\left\lvert \mu \right\rvert \leq r_i`.
The default values of :math:`r_1` and :math:`r_2` usually strike a good compromise between stability and sparsity.
For large and relatively dense problems, setting :math:`r_1` and :math:`r_2` to :math:`25` (say) may give a marked improvement in sparsity without impairing stability to a serious degree.
Note that for band matrices it may be necessary to set :math:`r_1` in the range :math:`1\leq r_1 < 2` in order to achieve stability.
If :math:`r_1 < 1` or :math:`r_2 < 1`, the default value is used.
**'LU Update Tolerance'** : float
Default :math:`\text{} = 10.0`
The values of :math:`r_1` and :math:`r_2` affect the stability and sparsity of the basis factorization :math:`B = LU`, during refactorization and updates respectively.
The lower triangular matrix :math:`L` is a product of matrices of the form
.. math::
\begin{pmatrix}1& \\\mu &1\end{pmatrix}
where the multipliers :math:`\mu` will satisfy :math:`\left\lvert \mu \right\rvert \leq r_i`.
The default values of :math:`r_1` and :math:`r_2` usually strike a good compromise between stability and sparsity.
For large and relatively dense problems, setting :math:`r_1` and :math:`r_2` to :math:`25` (say) may give a marked improvement in sparsity without impairing stability to a serious degree.
Note that for band matrices it may be necessary to set :math:`r_1` in the range :math:`1\leq r_1 < 2` in order to achieve stability.
If :math:`r_1 < 1` or :math:`r_2 < 1`, the default value is used.
**'LU Singularity Tolerance'** : float
Default :math:`\text{} = \epsilon^{0.67}`
If :math:`r > 0`, :math:`r` defines the singularity tolerance used to guard against ill-conditioned basis matrices.
Whenever the basis is refactorized, the diagonal elements of :math:`U` are tested as follows.
If :math:`\left\lvert u_{{jj}}\right\rvert \leq r` or :math:`\left\lvert u_{{jj}}\right\rvert < r\times \mathrm{max}_i\left\lvert u_{{ij}}\right\rvert`, the :math:`j`\ th column of the basis is replaced by the corresponding slack variable.
If :math:`r\leq 0`, the default value is used.
**'Minimize'** : valueless
Default
This option specifies the required direction of the optimization.
It applies to both linear and nonlinear terms (if any) in the objective function.
Note that if two problems are the same except that one minimizes :math:`f\left(x\right)` and the other maximizes :math:`{-f}\left(x\right)`, their solutions will be the same but the signs of the dual variables :math:`\pi_i` and the reduced gradients :math:`d_j` (see `The Main Iteration <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#ad-mainiteration>`__) will be reversed.
**'Maximize'** : valueless
This option specifies the required direction of the optimization.
It applies to both linear and nonlinear terms (if any) in the objective function.
Note that if two problems are the same except that one minimizes :math:`f\left(x\right)` and the other maximizes :math:`{-f}\left(x\right)`, their solutions will be the same but the signs of the dual variables :math:`\pi_i` and the reduced gradients :math:`d_j` (see `The Main Iteration <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#ad-mainiteration>`__) will be reversed.
**'Monitoring File'** : int
Default :math:`\text{} = {-1}`
If :math:`i\geq 0` and :math:`\text{‘Print Level'} > 0`, monitoring information produced by ``iqp_sparse`` is sent to a file with logical unit number :math:`i`.
If :math:`i < 0` and/or :math:`\text{‘Print Level'} = 0`, the default value is used and hence no monitoring information is produced.
**'Optimality Tolerance'** : float
Default :math:`\text{} = \mathrm{max}\left(10^{{-6}}, {\sqrt{\epsilon }}\right)`
If :math:`r\geq \epsilon`, :math:`r` is used to judge the size of the reduced gradients :math:`d_j = g_j-\pi^{\mathrm{T}}a_j`.
By definition, the reduced gradients for basic variables are always zero.
Optimality is declared if the reduced gradients for any nonbasic variables at their lower or upper bounds satisfy :math:`{-r}\times \mathrm{max}\left(1, {\left\lVert \pi \right\rVert }\right)\leq d_j\leq r\times \mathrm{max}\left(1, {\left\lVert \pi \right\rVert }\right)`, and if :math:`\left\lvert d_j\right\rvert \leq r\times \mathrm{max}\left(1, {\left\lVert \pi \right\rVert }\right)` for any superbasic variables.
If :math:`r < \epsilon`, the default value is used.
**'Partial Price'** : int
Default :math:`\text{} = 10`
Note that this option does not apply to QP problems.
This option is recommended for large FP or LP problems that have significantly more variables than constraints (i.e., :math:`n≫m`).
It reduces the work required for each pricing operation (i.e., when a nonbasic variable is selected to enter the basis).
If :math:`i = 1`, all columns of the constraint matrix :math:`\left(A-I\right)` are searched.
If :math:`i > 1`, :math:`A` and :math:`{-I}` are partitioned to give :math:`i` roughly equal segments :math:`A_{\textit{j}},K_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,p` (modulo :math:`p`).
If the previous pricing search was successful on :math:`A_{{j-1}},K_{{j-1}}`, the next search begins on the segments :math:`A_j,K_j`.
If a reduced gradient is found that is larger than some dynamic tolerance, the variable with the largest such reduced gradient (of appropriate sign) is selected to enter the basis.
If nothing is found, the search continues on the next segments :math:`A_{{j+1}},K_{{j+1}}`, and so on.
If :math:`i\leq 0`, the default value is used.
**'Pivot Tolerance'** : float
Default :math:`\text{} = \epsilon^{0.67}`
If :math:`r > 0`, :math:`r` is used to prevent columns entering the basis if they would cause the basis to become almost singular.
If :math:`r\leq 0`, the default value is used.
**'Print Level'** : int
Default :math:`\text{} = 10`
The value of :math:`i` controls the amount of printout produced by ``iqp_sparse``, as indicated below.
A detailed description of the printed output is given in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#fc-printedoutput>`__ (summary output at each iteration and the final solution) and `Monitoring Information <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cef.html#monitoring>`__ (monitoring information at each iteration).
Note that the summary output will not exceed :math:`80` characters per line and that the monitoring information will not exceed :math:`120` characters per line.
If :math:`i < 0`, the default value is used.
The following printout is sent to the file object associated with the advisory I/O unit (see :class:`~naginterfaces.base.utils.FileObjManager`):
.. rst-class:: nag-rules-none nag-align-left
+----------------------+----------------------------------------------------------------------------------+
|:math:`i` |**Output** |
+======================+==================================================================================+
|:math:`0` |No output. |
+----------------------+----------------------------------------------------------------------------------+
|:math:`1` |The final solution only. |
+----------------------+----------------------------------------------------------------------------------+
|:math:`5` |One line of summary output for each iteration (no printout of the final solution).|
+----------------------+----------------------------------------------------------------------------------+
|:math:`\text{}\geq 10`|The final solution and one line of summary output for each iteration. |
+----------------------+----------------------------------------------------------------------------------+
The following printout is sent to the logical unit number defined by the 'Monitoring File':
.. rst-class:: nag-rules-none nag-align-left
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`i` |**Output** |
+======================+================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================+
|:math:`0` |No output. |
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`1` |The final solution only. |
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`5` |One long line of output for each iteration (no printout of the final solution). |
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 10`|The final solution and one long line of output for each iteration. |
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`\text{}\geq 20`|The final solution, one long line of output for each iteration, matrix statistics (initial status of rows and columns, number of elements, density, biggest and smallest elements, etc.), details of the scale factors resulting from the scaling procedure (if :math:`\text{‘Scale Option'} = 1` or :math:`2`), basis factorization statistics and details of the initial basis resulting from the crash procedure (if :math:`\mathrm{start} = \texttt{'C'}`; see :ref:`Parameters <h02ce-py2-py-parameters>`).|
+----------------------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
If :math:`\text{‘Print Level'} > 0` and the unit number defined by 'Monitoring File' is the advisory unit number, then the summary output is suppressed.
**'Rank Tolerance'** : float
Default :math:`\text{} = 100\epsilon`
See above.
**'Scale Option'** : int
Default :math:`\text{} = 2`
This option enables you to scale the variables and constraints using an iterative procedure due to Fourer (see Hock and Schittkowski (1981)), which attempts to compute row scales :math:`r_i` and column scales :math:`c_j` such that the scaled matrix coefficients :math:`\bar{a}_{{ij}} = a_{{ij}}\times \left(c_j/r_i\right)` are as close as possible to unity.
This may improve the overall efficiency of the function on some problems. (The lower and upper bounds on the variables and slacks for the scaled problem are redefined as :math:`\bar{l}_j = l_j/c_j` and :math:`\bar{u}_j = u_j/c_j` respectively, where :math:`c_j\equiv r_{{j-n}}` if :math:`j > n`.)
If :math:`i = 0`, no scaling is performed.
If :math:`i = 1`, all rows and columns of the constraint matrix :math:`A` are scaled.
If :math:`i = 2`, an additional scaling is performed that may be helpful when the solution :math:`x` is large; it takes into account columns of :math:`\left(A-I\right)` that are fixed or have positive lower bounds or negative upper bounds.
If :math:`i < 0` or :math:`i > 2`, the default value is used.
**'Scale Tolerance'** : float
Default :math:`\text{} = 0.9`
Note that this option does not apply when :math:`\text{‘Scale Option'} = 0`.
If :math:`0 < r < 1`, :math:`r` is used to control the number of scaling passes to be made through the constraint matrix :math:`A`.
At least :math:`3` (and at most :math:`10`) passes will be made.
More precisely, let :math:`a_p` denote the largest column ratio (i.e., :math:`\frac{{\text{‘biggest'}\text{element}}}{{\text{‘smallest'}\text{element}}}` in some sense) after the :math:`p`\ th scaling pass through :math:`A`.
The scaling procedure is terminated if :math:`a_p\geq a_{{p-1}}\times r` for some :math:`p\geq 3`.
Thus, increasing the value of :math:`r` from :math:`0.9` to :math:`0.99` (say) will probably increase the number of passes through :math:`A`.
If :math:`r\leq 0` or :math:`r\geq 1`, the default value is used.
**'Superbasics Limit'** : int
Default :math:`\text{} = \mathrm{min}\left({n_H+1}, {n}\right)`
Note that this option does not apply to FP or LP problems.
The value of :math:`i` specifies 'how nonlinear' you expect the QP problem to be.
If :math:`i\leq 0`, the default value is used.
.. _h02ce-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{strtgy} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{strtgy} = 0`, :math:`1`, :math:`2` or :math:`3`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{mdepth} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mdepth}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{lintvr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{lintvr}\geq 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{istate}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{istate}[\textit{i}]\leq 3` with :math:`\mathrm{start} = \texttt{'W'}`.
(`errno` :math:`1`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{istate}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{istate}[\textit{i}]\leq 5` with :math:`\mathrm{start} = \texttt{'C'}`.
(`errno` :math:`1`)
On entry, :math:`\textit{nname} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nname} = 1` or :math:`\textit{nname} = n+\mathrm{m}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{start} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{start} = \texttt{'C'}` or :math:`\texttt{'W'}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ka}[\textit{i}] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{ka}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle` :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq {\mathrm{ka}[\textit{i}]-\mathrm{ka}[\textit{i}-1]}\leq \mathrm{m}`.
(`errno` :math:`1`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ka}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ka}[\textit{i}-1]\geq 1`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{ka}[n] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ka}[n] = \textit{nnz}+1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ka}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ka}[0] = 1`.
(`errno` :math:`1`)
On entry, :math:`\textit{i} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{ha}[\textit{i}-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{ha}[\textit{i}]\leq \mathrm{m}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ncolh} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{ncolh}\leq n`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iobj} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \mathrm{iobj}\leq \mathrm{m}`.
(`errno` :math:`1`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`0\leq \textit{nnz}\leq n\times \mathrm{m}`
(`errno` :math:`1`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m} > 0`.
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
(`errno` :math:`2`)
No integer solutions found.
(`errno` :math:`3`)
Need to increase depth, :math:`\mathrm{mdepth} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
The problem is unbounded (or badly scaled).
(`errno` :math:`5`)
The problem is infeasible.
(`errno` :math:`6`)
Too many iterations. Try increasing Iteration Limit.
(`errno` :math:`7`)
Reduced Hessian exceeds assigned dimension.
(`errno` :math:`8`)
Hessian appears to be indefinite.
(`errno` :math:`9`)
Invalid input argument for internal call to local optimizer.
(`errno` :math:`10`)
Numerical error trying to satisfy constraints. The basis is very ill-conditioned.
(`errno` :math:`13`)
Singular basis after :math:`15` attempts to factorize.
(`errno` :math:`16`)
Unexpected error in NAG function -- contact `NAG <https://www.nag.com>`__.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`-1`)
Halted at your request.
.. _h02ce-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_sparse`` is designed to obtain integer solutions to a class of quadratic programming problems addressed by :meth:`opt.qpconvex1_sparse_solve <naginterfaces.library.opt.qpconvex1_sparse_solve>`.
Specifically it solves the following problem:
.. math::
\begin{array}{l}\mathrm{minimize}_{{x \in R^n}}f\left(x\right)\quad \text{ subject to }\quad l\leq \left\{\begin{array}{c}x\\Ax\end{array}\right\} \leq u\text{,}\end{array}
where :math:`x = \left(x_1,x_2,\ldots,x_n\right)^\mathrm{T}` is a set of variables (some of which may be required to be integer), :math:`A` is an :math:`m\times n` matrix and the objective function :math:`f\left(x\right)` may be specified in a variety of ways depending upon the particular problem to be solved.
The option 'Maximize' may be used to specify an alternative problem in which :math:`f\left(x\right)` is maximized.
The possible forms for :math:`f\left(x\right)` are listed in Table [label omitted], in which the prefixes LP and QP stand for 'linear programming' and 'quadratic programming' respectively, :math:`c` is an :math:`n`-element vector and :math:`H` is the :math:`n\times n` second-derivative matrix :math:`\nabla^2f\left(x\right)` (the `Hessian matrix`).
+------------+---------------------------------------------------+-------------------------------+
|Problem type|Objective function :math:`f\left(x\right)` |Hessian matrix :math:`H` |
+============+===================================================+===============================+
|LP |:math:`c^{\mathrm{T}}x` |Not applicable |
+------------+---------------------------------------------------+-------------------------------+
|QP |:math:`c^{\mathrm{T}}x+\frac{1}{2}x^{\mathrm{T}}Hx`|Symmetric positive semidefinite|
+------------+---------------------------------------------------+-------------------------------+
For LP and QP problems, the unique global minimum value of :math:`f\left(x\right)` is found.
For QP problems, you must also provide a function that computes :math:`Hx` for any given vector :math:`x`. (:math:`H` need not be stored explicitly.)
(It is not expected that the feasibility problem of :meth:`opt.qpconvex1_sparse_solve <naginterfaces.library.opt.qpconvex1_sparse_solve>` would be relevant here.)
The function employs a 'Branch and Bound' technique to enforce the integer constraints.
In this technique the problem is first solved without the integer constraints.
If a variable is found to be non-integral when it is required to have an integer value then two additional problems are constructed.
One bounds the variable above by the nearest integer value below the optimal value previously obtained.
The second problem is formed by bounding the variable below by the nearest integer value above the optimal value.
This process is continued until an integer solution is found.
At this point you may elect to stop, or may continue to search for better integer solutions by examining any other sub-problems that remain to be explained.
In practice the function tries to compute an integer solution as quickly as possible using a depth-first approach, since this helps determine a realistic cut-off value.
If we have a cut-off value, say the value of the function at this first integer solution, and any sub-problem, :math:`W` say, has a solution value greater than this cut-off value, then subsequent sub-problems of :math:`W` must have solutions greater than the value of the solution at :math:`W` and, therefore, need not be computed.
Thus a knowledge of a good cut-off value can result in fewer sub-problems being solved and thus speed up the operation of the function. (See the description of :math:`\mathrm{monit}` in :ref:`Parameters <h02ce-py2-py-parameters>` for details of how you can supply your own cut-off value.)
Each sub-problem is solved using :meth:`opt.qpconvex1_sparse_solve <naginterfaces.library.opt.qpconvex1_sparse_solve>`.
You are referred to the function document for :meth:`opt.qpconvex1_sparse_solve <naginterfaces.library.opt.qpconvex1_sparse_solve>` for details of the algorithm used.
.. _h02ce-py2-py-references:
**References**
Gill, P E, Hammarling, S, Murray, W, Saunders, M A and Wright, M H, 1986, `Users' guide for LSSOL (Version 1.0)`, Report SOL 86-1, Department of Operations Research, Stanford University
Gill, P E and Murray, W, 1978, `Numerically stable methods for quadratic programming`, Math. Programming (14), 349--372
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1986, `Some theoretical properties of an augmented Lagrangian merit function`, Report SOL, 86--6R, Department of Operations Research, Stanford University
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1989, `A practical anti-cycling procedure for linearly constrained optimization`, Math. Programming (45), 437--474
Gill, P E, Murray, W, Saunders, M A and Wright, M H, 1991, `Inertia-controlling methods for general quadratic programming`, SIAM Rev. (33), 1--36
Hock, W and Schittkowski, K, 1981, `Test Examples for Nonlinear Programming Codes. Lecture Notes in Economics and Mathematical Systems` (187), Springer--Verlag
Lawson, C L, Hanson, R J, Kincaid, D R and Krogh, F T, 1979, `Basic linear algebra supbrograms for Fortran usage`, ACM Trans. Math. Software (5), 308--325
Murtagh, B A and Saunders, M A, 1983, `MINOS 5.0 user's guide`, Report SOL 83-20, Department of Operations Research, Stanford University
"""
raise NotImplementedError
[docs]def iqp_sparse_optfile(ioptns, comm, io_manager=None):
r"""
``iqp_sparse_optfile`` may be used to supply options to :meth:`iqp_sparse` from an external file.
.. _h02cf-py2-py-doc:
For full information please refer to the NAG Library document for h02cf
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cff.html
.. _h02cf-py2-py-parameters:
**Parameters**
**ioptns** : int
The unit number (see :meth:`~naginterfaces.base.utils.FileObjManager.unit_from_fileobj`) of the options file to be read.
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qpconvex1_sparse_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
.. _h02cf-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{ioptns} = \langle\mathit{\boldsymbol{value}}\rangle`. Constraint: :math:`0\leq \mathrm{ioptns}\leq 2147483647`.
(`errno` :math:`2`)
``Begin`` was found, but end-of-file was found before ``End`` was found.
(`errno` :math:`3`)
End-of-file was found before ``Begin`` was found.
(`errno` :math:`5`)
The supplied option string is invalid. Check that the keywords are neither ambiguous nor misspelt. The option string is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _h02cf-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_sparse_optfile`` may be used to supply values for options to :meth:`iqp_sparse`. ``iqp_sparse_optfile`` reads an external file and each line of the file defines a single option.
It is only necessary to supply values for those arguments whose values are to be different from their default values.
Each option is defined by a single character string of up to :math:`72` characters, consisting of one or more items.
The items associated with a given option must be separated by spaces, or equal signs :math:`\left[ = \right]`.
Alphabetic characters may be upper or lower case.
The string
::
Print level = 1
is an example of a string used to set an option.
For each option the string contains one or more of the following items:
- a mandatory keyword;
- a phrase that qualifies the keyword;
- a number that specifies an `int` or `float` value. Such numbers may be up to :math:`16` contiguous characters in Fortran 77's I, F, E or D formats, terminated by a space if this is not the last item on the line.
Blank strings and comments are ignored.
A comment begins with an asterisk (*) and all subsequent characters in the string are regarded as part of the comment.
The file containing the options must start with ``Begin`` and must finish with ``End``.
An example of a valid options file is:
::
Begin * Example options file
Print Level = 1
End
Printing of user-supplied options is turned off by default, but may be turned on at any time using the keyword 'List'.
Option settings are preserved following a call to :meth:`iqp_sparse`, and so the keyword 'Defaults' is provided to allow you to reset all the options to their default values prior to a subsequent call to :meth:`iqp_sparse`.
A complete list of options, their abbreviations, synonyms and default values is given in :ref:`Other Parameters for iqp_sparse <h02ce-py2-py-other_params>`.
"""
raise NotImplementedError
[docs]def iqp_sparse_optstr(optstr, comm, io_manager=None):
r"""
``iqp_sparse_optstr`` may be used to supply individual options to :meth:`iqp_sparse`.
.. _h02cg-py2-py-doc:
For full information please refer to the NAG Library document for h02cg
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02cgf.html
.. _h02cg-py2-py-parameters:
**Parameters**
**optstr** : str
A single valid option string (as described in :ref:`Notes <h02cg-py2-py-notes>` above and in :ref:`Other Parameters for iqp_sparse <h02ce-py2-py-other_params>`).
**comm** : dict, communication object, modified in place
Communication structure.
This argument must have been initialized by a prior call to :meth:`opt.nlp1_init <naginterfaces.library.opt.nlp1_init>` with :math:`{\textit{rname}} = \texttt{'opt.qpconvex1_sparse_solve'}`.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
.. _h02cg-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`5`)
The supplied option string is invalid. Supplied value was: :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _h02cg-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``iqp_sparse_optstr`` may be used to supply values for options to :meth:`iqp_sparse`.
It is only necessary to call ``iqp_sparse_optstr`` for those arguments whose values are to be different from their default values.
One call to ``iqp_sparse_optstr`` sets one argument value.
Each option is defined by a single character string of up to :math:`72` characters, consisting of one or more items.
The items associated with a given option must be separated by spaces, or equal signs :math:`\left[ = \right]`.
Alphabetic characters may be upper or lower case.
The string
::
Print level = 1
is an example of a string used to set an option.
For each option the string contains one or more of the following items:
- a mandatory keyword;
- a phrase that qualifies the keyword;
- a number that specifies an `int` or `float` value. Such numbers may be up to :math:`16` contiguous characters in Fortran 77's I, F, E or D formats, terminated by a space if this is not the last item on the line.
Blank strings and comments are ignored.
A comment begins with an asterisk (*) and all subsequent characters in the string are regarded as part of the comment.
Printing of user-specified options is turned off by default.
It may be turned on at any time using the keyword 'List'.
Option settings are preserved following a call to :meth:`iqp_sparse`, and so the keyword 'Defaults' is provided to allow you to reset all the options to their default values.
A complete list of options, their abbreviations, synonyms and default values is given in :ref:`Other Parameters for iqp_sparse <h02ce-py2-py-other_params>`.
"""
raise NotImplementedError
[docs]def sqp(ncnln, bl, bu, varcon, x, objfun, comm, a=None, d=None, confun=None, maxit=500, acc=1.0e-4, data=None, io_manager=None, spiked_sorder='C'):
r"""
``sqp`` solves general nonlinear programming problems with integer constraints on some of the variables.
Note: this function uses optional algorithmic parameters, see also: :meth:`optset`.
.. _h02da-py2-py-doc:
For full information please refer to the NAG Library document for h02da
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02daf.html
.. _h02da-py2-py-parameters:
**Parameters**
**ncnln** : int
:math:`n_N`, the number of constraints supplied by :math:`c\left(x\right)`.
**bl** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{bl}` must contain the lower bounds, :math:`l_i`, and :math:`\mathrm{bu}` the upper bounds, :math:`u_i`, for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.
**bu** : float, array-like, shape :math:`\left(n\right)`
:math:`\mathrm{bl}` must contain the lower bounds, :math:`l_i`, and :math:`\mathrm{bu}` the upper bounds, :math:`u_i`, for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.
**varcon** : int, array-like, shape :math:`\left(n+\textit{nclin}+\mathrm{ncnln}\right)`
:math:`\mathrm{varcon}` indicates the nature of each variable and constraint in the problem. The first :math:`n` elements of the array must describe the nature of the variables, the next :math:`n_L` elements the nature of the general linear constraints (if any) and the next :math:`n_N` elements the general constraints (if any).
:math:`\mathrm{varcon}[j-1] = 0`
A continuous variable.
:math:`\mathrm{varcon}[j-1] = 1`
A binary variable.
:math:`\mathrm{varcon}[j-1] = 2`
An integer variable.
:math:`\mathrm{varcon}[j-1] = 3`
An equality constraint.
:math:`\mathrm{varcon}[j-1] = 4`
An inequality constraint.
**x** : float, array-like, shape :math:`\left(n\right)`
An initial estimate of the solution, which need not be feasible. Values corresponding to integer variables are rounded; if an initial value less than half is supplied for a binary variable the value zero is used, otherwise the value one is used.
**objfun** : callable (objmip, objgrd) = objfun(mode, varcon, x, objgrd, nstate, data=None)
:math:`\mathrm{objfun}` must calculate the objective function :math:`f\left(x\right)` and its gradient for a specified :math:`n`-element vector :math:`x`.
**Parameters**
**mode** : int
Indicates which values must be assigned during each call of :math:`\mathrm{objfun}`. Only the following values need be assigned:
:math:`\mathrm{mode} = 0`
The objective function value, :math:`\mathrm{objmip}`.
:math:`\mathrm{mode} = 1`
The continuous elements of :math:`\mathrm{objgrd}`.
**varcon** : int, ndarray, shape :math:`\left(\textit{lvarcon}\right)`
The array :math:`\mathrm{varcon}` as supplied to ``sqp``.
**x** : float, ndarray, shape :math:`\left(n\right)`
The vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.
**objgrd** : float, ndarray, shape :math:`\left(n\right)`
Continuous elements of :math:`\mathrm{objgrd}` are set to the value of NaN.
**nstate** : int
If :math:`\mathrm{nstate} = 1`, ``sqp`` is calling :math:`\mathrm{objfun}` for the first time. This argument setting allows you to save computation time if certain data must be read or calculated only once.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**objmip** : float
Must be set to the objective function value, :math:`f`.
**objgrd** : float, array-like, shape :math:`\left(n\right)`
Must contain the gradient vector of the objective function if :math:`\mathrm{mode} = 1`, with :math:`\mathrm{objgrd}[j-1]` containing the partial derivative of :math:`f` with respect to continuous variable :math:`x_j`; otherwise :math:`\mathrm{objgrd}` is not referenced.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by a prior call to :meth:`optset`.
**a** : None or float, array-like, shape :math:`\left(\textit{nclin}, :\right)`, optional
Note: the required extent for this argument in dimension 2 is determined as follows: if :math:`\textit{nclin} > 0`: :math:`n`; otherwise: :math:`0`.
The :math:`\textit{i}`\ th row of :math:`\mathrm{a}` must contain the coefficients of the :math:`\textit{i}`\ th general linear constraint, for :math:`\textit{i} = 1,2,\ldots,n_l`. Any equality constraints must be specified first.
If :math:`\textit{nclin} = 0`, the array :math:`\mathrm{a}` is not referenced.
**d** : None or float, array-like, shape :math:`\left(\textit{nclin}\right)`, optional
:math:`d_i`, the constant for the :math:`i`\ th linear constraint.
If :math:`\textit{nclin} = 0`, the array :math:`\mathrm{d}` is not referenced.
**confun** : None or callable (c, cjac) = confun(mode, varcon, x, cjac, nstate, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.
:math:`\mathrm{confun}` must calculate the constraint functions supplied by :math:`c\left(x\right)` and their Jacobian at :math:`x`. If all constraints are supplied by :math:`A` and :math:`d` (i.e., :math:`\mathrm{ncnln} = 0`), :math:`\mathrm{confun}` will never be called by ``sqp`` and :math:`\mathrm{confun}` may be **None**. If :math:`\mathrm{ncnln} > 0`, the first call to :math:`\mathrm{confun}` will occur after the first call to :math:`\mathrm{objfun}`.
**Parameters**
**mode** : int
Indicates which values must be assigned during each call of :math:`\mathrm{objfun}`. Only the following values need be assigned:
:math:`\mathrm{mode} = 0`
Elements of :math:`\mathrm{c}` containing continuous variables.
:math:`\mathrm{mode} = 1`
Elements of :math:`\mathrm{cjac}` containing continuous variables.
**varcon** : int, ndarray, shape :math:`\left(\textit{lvarcon}\right)`
The array :math:`\mathrm{varcon}` as supplied to ``sqp``.
**x** : float, ndarray, shape :math:`\left(n\right)`
The vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.
**cjac** : float, ndarray, shape :math:`\left(\textit{ncnln}, n\right)`
Note: the derivative of the :math:`i`\ th constraint with respect to the :math:`j`\ th variable, :math:`\frac{{\partial c_i}}{{\partial x_j}}`, is stored in :math:`\mathrm{cjac}[i-1,j-1]`.
Continuous elements of :math:`\mathrm{cjac}` are set to the value of NaN.
**nstate** : int
If :math:`\mathrm{nstate} = 1`, ``sqp`` is calling :math:`\mathrm{confun}` for the first time. This argument setting allows you to save computation time if certain data must be read or calculated only once.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**c** : float, array-like, shape :math:`\left(\textit{ncnln}\right)`
Must contain :math:`\mathrm{ncnln}` constraint values, with the value of the :math:`j`\ th constraint :math:`c_j\left(x\right)` in :math:`\mathrm{c}[j-1]`.
**cjac** : float, array-like, shape :math:`\left(\textit{ncnln}, n\right)`
The :math:`i`\ th row of :math:`\mathrm{cjac}` must contain elements of :math:`\frac{{\partial c_i}}{{\partial x_j}}` for each continuous variable :math:`x_j`.
**maxit** : int, optional
The maximum number of iterations within which to find a solution. If :math:`\mathrm{maxit}` is less than or equal to zero, the suggested value below is used.
**acc** : float, optional
The requested accuracy for determining feasible points during iterations and for halting the method when the predicted improvement in objective function is less than :math:`\mathrm{acc}`. If :math:`\mathrm{acc}` is less than or equal to :math:`\epsilon` (:math:`\epsilon` being the machine precision as given by :meth:`machine.precision <naginterfaces.library.machine.precision>`), the below suggested value is used.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**spiked_sorder** : str, optional
If :math:`\mathrm{a}` is spiked (i.e., has unit extent in all but one dimension, or has size :math:`1`), :math:`\mathrm{spiked\_sorder}` selects the storage order to associate with it in the NAG Engine:
spiked_sorder = :math:`\texttt{'C'}`
row-major storage will be used;
spiked_sorder = :math:`\texttt{'F'}`
column-major storage will be used.
Two-dimensional arrays returned from callback functions in this routine must then use the same storage order.
**Returns**
**ax** : float, ndarray, shape :math:`\left(\textit{nclin}\right)`
The final residuals of the linear constraints :math:`Ax-d`.
If :math:`\textit{nclin} = 0`, :math:`\mathrm{ax}` is not referenced.
**x** : float, ndarray, shape :math:`\left(n\right)`
The final estimate of the solution.
**c** : float, ndarray, shape :math:`\left(\mathrm{ncnln}\right)`
If :math:`\mathrm{ncnln} > 0`, :math:`\mathrm{c}[\textit{j}-1]` contains the value of the :math:`\textit{j}`\ th constraint function :math:`c_{\textit{j}}\left(x\right)` at the final iterate, for :math:`\textit{j} = 1,2,\ldots,\mathrm{ncnln}`.
If :math:`\mathrm{ncnln} = 0`, the array :math:`\mathrm{c}` is not referenced.
**cjac** : float, ndarray, shape :math:`\left(\mathrm{ncnln}, n\right)`
If :math:`\mathrm{ncnln} > 0`, :math:`\mathrm{cjac}` contains the Jacobian matrix of the constraint functions at the final iterate, i.e., :math:`\mathrm{cjac}[\textit{i}-1,\textit{j}-1]` contains the partial derivative of the :math:`\textit{i}`\ th constraint function with respect to the :math:`\textit{j}`\ th variable, for :math:`\textit{j} = 1,2,\ldots,n`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{ncnln}`. (See the discussion of argument :math:`\mathrm{cjac}` under :math:`\mathrm{confun}`.)
If :math:`\mathrm{ncnln} = 0`, the array :math:`\mathrm{cjac}` is not referenced.
**objgrd** : float, ndarray, shape :math:`\left(n\right)`
The objective function gradient at the solution.
**objmip** : float
With no exception or warning is raised, :math:`\mathrm{objmip}` contains the value of the objective function for the MINLP solution.
.. _h02da-py2-py-other_params:
**Other Parameters**
**'Branch Bound Steps'** : int
Default :math:`\text{} = 500`
Maximum number of branch-and-bound steps for solving the mixed integer quadratic problems.
**'Branching Rule'** : str
Default :math:`\text{} = \texttt{'Maximum'}`
Branching rule for branch and bound search.
:math:`\text{‘Branching Rule'} = \texttt{'Maximum'}`
Maximum fractional branching.
:math:`\text{‘Branching Rule'} = \texttt{'Minimum'}`
Minimum fractional branching.
**'Check Gradients'** : str
Default :math:`\text{} = \texttt{'No'}`
Perform an internal check of supplied objective and constraint gradients.
It is advisable to set :math:`\text{‘Check Gradients'} = \texttt{'Yes'}` during code development to avoid difficulties associated with incorrect user-supplied gradients.
**'Continuous Trust Radius'** : float
Default :math:`\text{} = 10.0`
Initial continuous trust region radius, :math:`\Delta_0^c`; the initial trial step :math:`d \in R^{n_c}` for the SQP approximation must satisfy :math:`\left\lVert d\right\rVert_\infty\leq \Delta_0^c`.
**'Descent'** : float
Default :math:`\text{} = 0.05`
Initial descent parameter, :math:`\delta`, larger values of :math:`\delta` allow penalty option :math:`\sigma` to increase faster which can lead to faster convergence.
**'Descent Factor'** : float
Default :math:`\text{} = 0.1`
Factor for decreasing the internal descent parameter, :math:`\delta`, between iterations.
**'Feasible Steps'** : int
Default :math:`\text{} = 10`
Maximum number of feasible steps without improvements, where feasibility is measured by :math:`\left\lVert g\left(x\right)\right\rVert_\infty\leq \sqrt{\mathrm{acc}}`.
**'Improved Bounds'** : str
Default :math:`\text{} = \texttt{'No'}`
Calculate improved bounds in case of 'Best of all' node selection strategy.
**'Integer Trust Radius'** : float
Default :math:`\text{} = 10.0`
Initial integer trust region radius, :math:`\Delta_0^i`; the initial trial step :math:`e \in R^{n_i}` for the SQP approximation must satisfy :math:`\left\lVert e\right\rVert_\infty\leq \Delta_0^i`.
**'Maximum Restarts'** : int
Default :math:`\text{} = 2`
Maximum number of restarts that allow the mixed integer SQP algorithm to return to a better solution.
Setting a value higher than the default might lead to better results at the expense of function evaluations.
**'Minor Print Level'** : int
Default :math:`\text{} = 0`
Print level of the subproblem solver.
Active only if :math:`\text{‘Print Level'}\neq 0`.
**'Modify Hessian'** : str
Default :math:`\text{} = \texttt{'Yes'}`
Modify the Hessian approximation in an attempt to get more accurate search directions.
Calculation time is increased when the number of integer variables is large.
**'Node Selection'** : str
Default :math:`\text{} = \texttt{'Depth First'}`
Node selection strategy for branch and bound.
:math:`\text{‘Node Selection'} = \texttt{'Best of all'}`
Large tree search; this method is the slowest as it solves all subproblem QPs independently.
:math:`\text{‘Node Selection'} = \texttt{'Best of two'}`
Uses warm starts and less memory.
:math:`\text{‘Node Selection'} = \texttt{'Depth first'}`
Uses more warm starts. If warm starts are applied, they can speed up the solution of mixed integer subproblems significantly when solving almost identical QPs.
**'Non Monotone'** : int
Default :math:`\text{} = 10`
Maximum number of successive iterations considered for the non-monotone trust region algorithm, allowing the penalty function to increase.
**'Objective Scale Bound'** : float
Default :math:`\text{} = 1.0`
When :math:`\text{‘Scale Objective Function'} > 0` internally scale absolute function values greater than :math:`1.0` or 'Objective Scale Bound'.
**'Penalty'** : float
Default :math:`\text{} = 1000.0`
Initial penalty option, :math:`\sigma`.
**'Penalty Factor'** : float
Default :math:`\text{} = 10.0`
Factor for increasing penalty option :math:`\sigma` when the trust regions cannot be enlarged at a trial step.
**'Print Level'** : int
Default :math:`\text{} = 0`
Specifies the desired output level of printing.
:math:`\text{‘Print Level'} = \texttt{'0'}`
No output.
:math:`\text{‘Print Level'} = \texttt{'1'}`
Final convergence analysis.
:math:`\text{‘Print Level'} = \texttt{'2'}`
One line of intermediate results per iteration.
:math:`\text{‘Print Level'} = \texttt{'3'}`
Detailed information printed per iteration.
**'QP Accuracy'** : float
Default :math:`\text{} = 1.0e-10`
Termination tolerance of the relaxed quadratic program subproblems.
**'Scale Continuous Variables'** : str
Default :math:`\text{} = \texttt{'Yes'}`
Internally scale continuous variables values.
**'Scale Objective Function'** : int
Default :math:`\text{} = 1`
Internally scale objective function values.
:math:`\text{‘Scale Objective Function'} = \texttt{'0'}`
No scaling.
:math:`\text{‘Scale Objective Function'} = \texttt{'1'}`
Scale absolute values greater than 'Objective Scale Bound'.
**'Warm Starts'** : int
Default :math:`\text{} = 100`
Maximum number of warm starts within the mixed integer QP solver, see 'Node Selection'.
.. _h02da-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`n = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`n > 0`.
(`errno` :math:`2`)
On entry, :math:`\textit{nclin} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nclin}\geq 0`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{ncnln} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ncnln}\geq 0`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{bl}[\langle\mathit{\boldsymbol{value}}\rangle] > \mathrm{bu}[\langle\mathit{\boldsymbol{value}}\rangle]`.
Constraint: :math:`\mathrm{bl}[\textit{i}-1]\leq \mathrm{bu}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,n`.
(`errno` :math:`6`)
On entry, :math:`\mathrm{varcon}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{varcon}[\textit{i}-1] = 0`, :math:`1` or :math:`2`, for :math:`\textit{i} = 1,2,\ldots,n`.
(`errno` :math:`7`)
On entry, :math:`\mathrm{varcon}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{varcon}[\textit{i}-1] = 3` or :math:`4`, for :math:`\textit{i} = n+1,\ldots,n+\textit{nclin}+\mathrm{ncnln}`.
(`errno` :math:`8`)
The supplied :math:`\mathrm{objfun}` returned a NaN value.
(`errno` :math:`9`)
The supplied :math:`\mathrm{confun}` returned a NaN value.
(`errno` :math:`10`)
On entry, the option arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] have either not been initialized or been corrupted.
(`errno` :math:`11`)
On entry, there are no binary or integer variables.
(`errno` :math:`12`)
On entry, linear equality constraints do not precede linear inequality constraints.
(`errno` :math:`13`)
On entry, nonlinear equality constraints do not precede nonlinear inequality constraints.
(`errno` :math:`81`)
One or more objective gradients appear to be incorrect.
(`errno` :math:`91`)
One or more constraint gradients appear to be incorrect.
(`errno` :math:`1002`)
More than the maximum number of feasible steps without improvement in the objective function. If the maximum number of feasible steps is small, say less than :math:`5`, try increasing it. Option :math:`\text{‘Feasible Steps'} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1003`)
Penalty parameter tends to infinity in an underlying mixed-integer quadratic program; the problem may be infeasible. If :math:`\sigma` is relatively low value, try a higher one, for example :math:`10^{20}`.
Option :math:`\text{‘Penalty'} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1004`)
Termination at an infeasible iterate; if the problem is feasible, try a different starting value.
(`errno` :math:`1005`)
Termination with zero integer trust region for integer variables; try a different starting value.
Option :math:`\text{‘Integer Trust Radius'} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1008`)
The optimization failed due to numerical difficulties. Set option :math:`\text{‘Print Level'} = 3` for more information.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`1001`)
On entry, :math:`\mathrm{maxit} = \langle\mathit{\boldsymbol{value}}\rangle`. Exceeded the maximum number of iterations.
**NagCallbackTerminateWarning**
(`errno` :math:`i < 0`)
The optimization halted because you set :math:`\mathrm{mode}` negative in :math:`\mathrm{objfun}` or :math:`\mathrm{mode}` negative in :math:`\mathrm{confun}`, to :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
.. _h02da-py2-py-notes:
**Notes**
``sqp`` solves mixed integer nonlinear programming problems using a modified sequential quadratic programming method.
The problem is assumed to be stated in the following general form:
.. math::
\begin{array}{lll} \mathrm{minimize}_{{x \in \left\{R^{n_c}, Z^{n_i}\right\}}} & f\left(x\right) &\\\text{subject to}&c_j\left(x\right) = 0\text{,}&j = 1,2,\ldots,m_e\\&c_j\left(x\right)\geq 0\text{,}&j = m_e+1,m_e+2,\ldots,m\\&l\leq x_i\leq u\text{,}&i = 1,2,\ldots,n\end{array}
with :math:`n_c` continuous variables and :math:`n_i` binary and integer variables in a total of :math:`n` variables; :math:`m_e` equality constraints in a total of :math:`m` constraint functions.
Partial derivatives of :math:`f\left(x\right)` and :math:`c\left(x\right)` are not required for the :math:`n_i` integer variables.
Gradients with respect to integer variables are approximated by difference formulae.
No assumptions are made regarding :math:`f\left(x\right)` except that it is twice continuously differentiable with respect to continuous elements of :math:`x`.
It is not assumed that integer variables are relaxable.
In other words, problem functions are evaluated only at integer points.
The method seeks to minimize the exact penalty function:
.. math::
P_{\sigma }\left(x\right) = f\left(x\right)+\sigma \left\lVert g\left(x\right)\right\rVert_{\infty }
where :math:`\sigma` is adapted by the algorithm and :math:`g\left(x\right)` is given by:
.. math::
\begin{array}{lcll}g\left(x\right)& = &c_j\left(x\right)\text{,}&j = 1,2,\ldots,m_e\\& = &\mathrm{min}\left({c_j\left(x\right)}, 0\right)\text{,}&j = m_e+1,m_e+2,\ldots,m\text{.}\end{array}
Successive quadratic approximations are applied under the assumption that integer variables have a smooth influence on the model functions, that is function values do not change drastically when incrementing or decrementing an integer value.
In practice this requires integer variables to be ordinal not categorical.
The algorithm is stabilised by a trust region method including Yuan's second order corrections, see Yuan and Sun (2006).
The Hessian of the Lagrangian function is approximated by BFGS (see `The Quasi-Newton Update for opt.nlp1_solve <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/e04/e04ucf.html#ad-quasinewton>`__) updates subject to the continuous and integer variables.
The mixed-integer quadratic programming subproblems of the SQP-trust region method are solved by a branch and cut method with continuous subproblem solutions obtained by the primal-dual method of Goldfarb and Idnani, see Powell (1983).
Different strategies are available for selecting a branching variable:
`Maximal fractional branching.` Select an integer variable from the relaxed subproblem solution with largest distance from next integer value
`Minimal fractional branching.` Select an integer variable from the relaxed subproblem solution with smallest distance from next integer value
and a node from where branching, that is the generation of two new subproblems, begins:
`Best of two.` The optimal objective function values of the two child nodes are compared and the node with a lower value is chosen
`Best of all.` Select an integer variable from the relaxed subproblem solution with the smallest distance from the next integer value
`Depth first.` Select a child node whenever possible.
This implementation is based on the algorithm MISQP as described in Exler `et al.` (2013).
Linear constraints may optionally be supplied by a matrix :math:`A` and vector :math:`d` rather than the constraint functions :math:`c\left(x\right)` such that
.. math::
Ax = d\quad \text{ or }Ax\geq d\text{.}
Partial derivatives with respect to :math:`x` of these constraint functions are not requested by ``sqp``.
.. _h02da-py2-py-references:
**References**
Exler, O, Lehmann, T and Schittkowski, K, 2013, `A comparative study of SQP-type algorithms for nonlinear and nonconvex mixed-integer optimization`, Mathematical Programming Computation (4), 383--412
Mann, A, 1986, `GAMS/MINOS: Three examples`, Department of Operations Research Technical Report, Stanford University
Powell, M J D, 1983, `On the quadratic programming algorithm of Goldfarb and Idnani`, Report DAMTP 1983/Na 19, University of Cambridge, Cambridge
Yuan, Y-x and Sun, W, 2006, `Optimization Theory and Methods`, Springer--Verlag
See Also
--------
:meth:`naginterfaces.library.examples.mip.sqp_ex.main`
"""
raise NotImplementedError
[docs]def optset(optstr, comm):
r"""
``optset`` either initializes or resets the option arrays or sets a single option for supported problem solving functions in submodule ``mip``.
Currently, only :meth:`sqp` is supported.
.. _h02zk-py2-py-doc:
For full information please refer to the NAG Library document for h02zk
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02zkf.html
.. _h02zk-py2-py-parameters:
**Parameters**
**optstr** : str
A string identifying the option to be set.
Initialize = :math:`\textit{function name}`
Initialize the option arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] for use with function :math:`\textit{function name}`, where :math:`\textit{function name}` is the name associated with the function of interest.
Defaults
Resets all options to their default values.
:math:`\textit{option} = \textit{optval}`
See :ref:`Other Parameters for sqp <h02da-py2-py-other_params>` for details of valid values for :math:`\textit{option}` and :math:`\textit{optval}`. The equals sign (:math:`=`) delimiter must be used to separate the :math:`\textit{option}` from its :math:`\textit{optval}` value.
:math:`\mathrm{optstr}` is case insensitive.
Each token in the :math:`\textit{option}` and :math:`\textit{optval}` component must be separated by at least one space.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
.. _h02zk-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, the option in :math:`\mathrm{optstr}` was not recognized: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`12`)
On entry, the expected delimiter ':math:`=`' was not found in :math:`\mathrm{optstr}`: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`13`)
On entry, could not convert the specified :math:`\textit{optval}` to an integer: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`13`)
On entry, could not convert the specified :math:`\textit{optval}` to a real: :math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`14`)
On entry, attempting to initialize the option arrays but specified function name was not valid: :math:`\text{name} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`15`)
On entry, the :math:`\textit{optval}` supplied for the integer option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`16`)
On entry, the :math:`\textit{optval}` supplied for the real option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`17`)
On entry, the :math:`\textit{optval}` supplied for the character option is not valid.
:math:`\mathrm{optstr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`21`)
On entry, either the option arrays have not been initialized or they have been corrupted.
.. _h02zk-py2-py-notes:
**Notes**
``optset`` has three purposes: to initialize option arrays; to reset all options to their default values; or to set a single option to a user-supplied value.
Options and their values are, in general, presented as a character string, :math:`\mathrm{optstr}`, of the form ':math:`\textit{option} = \textit{optval}`'; alphabetic characters can be supplied in either upper or lower case.
Both :math:`\textit{option}` and :math:`\textit{optval}` may consist of one or more tokens separated by white space.
The tokens that comprise :math:`\textit{optval}` will normally be either an integer, real or character value as defined in the description of the specific optional argument.
In addition all options can take an :math:`\textit{optval}` DEFAULT which resets the option to its default value.
It is imperative that option arrays are initialized before any options are set, before the relevant problem solving function is called and before any options are queried using :meth:`optget`.
Information relating to available option names and their corresponding valid values is given in :ref:`Other Parameters for sqp <h02da-py2-py-other_params>`.
See Also
--------
:meth:`naginterfaces.library.examples.mip.sqp_ex.main`
"""
raise NotImplementedError
[docs]def optget(optstr, comm):
r"""
``optget`` is used to query the value of options available to supported problem solving functions in submodule ``mip``.
Currently, only :meth:`sqp` is supported.
.. _h02zl-py2-py-doc:
For full information please refer to the NAG Library document for h02zl
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h02zlf.html
.. _h02zl-py2-py-parameters:
**Parameters**
**optstr** : str
A string identifying the option whose current value is required. See :ref:`Other Parameters for sqp <h02da-py2-py-other_params>` for information on valid options. In addition, the following is a valid option:
Identify
``optget`` returns the function name supplied to :meth:`optset` when the option arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] were initialized.
**comm** : dict, communication object
Communication structure.
This argument must have been initialized by a prior call to :meth:`optset`.
**Returns**
**optvalue** : dict
The option-value ``dict``, with the following keys:
``'value'`` : float, int or str
The value of the requested option.
``'annotation'`` : None or str
Possible additional information about the option value.
.. _h02zl-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, the :math:`\textit{option}` in :math:`\mathrm{optstr}` has not been recognized.
(`errno` :math:`61`)
The arrays :math:`\mathrm{comm}`\ ['iopts'] and :math:`\mathrm{comm}`\ ['opts'] have either not been initialized, have become corrupted, or are not compatible with this option setting function.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`41`)
On entry, :math:`\mathrm{optstr}` indicates a character option, but :math:`\textit{cvalue}` is too short to hold the stored value. The returned value will be truncated.
.. _h02zl-py2-py-notes:
**Notes**
``optget`` is used to query the current values of options.
It is necessary to initialize option arrays using :meth:`optset` before any options are queried.
Information on option names and whether these options are real, integer or character can be found in :ref:`Other Parameters for sqp <h02da-py2-py-other_params>`.
"""
raise NotImplementedError
[docs]def transportation(kost, k15, maxit):
r"""
``transportation`` solves the classical Transportation ('Hitchcock') problem.
.. _h03ab-py2-py-doc:
For full information please refer to the NAG Library document for h03ab
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h03abf.html
.. _h03ab-py2-py-parameters:
**Parameters**
**kost** : int, array-like, shape :math:`\left(\textit{ma}, \textit{mb}\right)`
The coefficients :math:`c_{{\textit{i}\textit{j}}}`, for :math:`\textit{j} = 1,2,\ldots,m_b`, for :math:`\textit{i} = 1,2,\ldots,m_a`.
**k15** : int, array-like, shape :math:`\left(m\right)`
:math:`\mathrm{k15}[\textit{i}]` must be set to the availabilities :math:`A_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,m_a`; and :math:`\mathrm{k15}[m_a+j]` must be set to the requirements :math:`B_{\textit{j}}`, for :math:`\textit{j} = 1,2,\ldots,m_b`.
**maxit** : int
The maximum number of iterations allowed.
**Returns**
**k15** : int, ndarray, shape :math:`\left(m\right)`
The contents of :math:`\mathrm{k15}` are undefined.
**numit** : int
The number of iterations performed.
**k6** : int, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{k6}[\textit{k}-1]`, for :math:`\textit{k} = 1,2,\ldots,m_a+m_b-1`, contains the source indices of the optimal solution (see :math:`\mathrm{k11}`).
**k8** : int, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{k8}[\textit{k}-1]`, for :math:`\textit{k} = 1,2,\ldots,m_a+m_b-1`, contains the destination indices of the optimal solution (see :math:`\mathrm{k11}`).
**k11** : int, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{k11}[\textit{k}-1]`, for :math:`\textit{k} = 1,2,\ldots,m_a+m_b-1`, contains the optimal quantities :math:`x_{{ij}}` which, sent from source :math:`i = \mathrm{k6}[k-1]` to destination :math:`j = \mathrm{k8}[k-1]`, minimize :math:`z`.
**k12** : int, ndarray, shape :math:`\left(m\right)`
:math:`\mathrm{k12}[\textit{k}-1]`, for :math:`\textit{k} = 1,2,\ldots,m_a+m_b-1`, contains the unit cost :math:`c_{{ij}}` associated with the route from source :math:`i = \mathrm{k6}[k-1]` to destination :math:`j = \mathrm{k8}[k-1]`.
**z** : float
The value of the minimized total cost.
.. _h03ab-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, the sum of the availabilities does not equal the sum of the requirements.
(`errno` :math:`2`)
During computations :math:`\mathrm{maxit}` has been exceeded.
(`errno` :math:`3`)
On entry, :math:`\mathrm{maxit} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{maxit} > 0`.
(`errno` :math:`4`)
On entry, :math:`m = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ma} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{mb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`m = \textit{ma}+\textit{mb}`.
(`errno` :math:`4`)
On entry, :math:`\textit{ma} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\textit{ldkost} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ma}\leq \textit{ldkost}`.
(`errno` :math:`4`)
On entry, :math:`\textit{mb} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{mb} > 0`.
(`errno` :math:`4`)
On entry, :math:`\textit{ma} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ma} > 0`.
.. _h03ab-py2-py-notes:
**Notes**
`In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.`
``transportation`` solves the Transportation problem by minimizing
.. math::
z = \sum_i^{{m_a}}\sum_j^{{m_b}}c_{{ij}}x_{{ij}}\text{.}
subject to the constraints
.. math::
\begin{array}{l} \sum_j^{{m_b}} x_{{ij}} = A_i \quad \text{ (Availabilities)} \\ \sum_i^{{m_a}} \sum_i x_{{ij}} = B_j \quad \text{ (Requirements)} \end{array}
where the :math:`x_{{\textit{i}\textit{j}}}` can be interpreted as quantities of goods sent from source :math:`\textit{i}` to destination :math:`\textit{j}`, for :math:`\textit{j} = 1,2,\ldots,m_b`, for :math:`\textit{i} = 1,2,\ldots,m_a`, at a cost of :math:`c_{{\textit{i}\textit{j}}}` per unit, and it is assumed that :math:`\sum_i^{{m_a}}A_i = \sum_j^{{m_b}}\sum_jB_j` and :math:`x_{{ij}}\geq 0`.
``transportation`` uses the 'stepping stone' method, modified to accept degenerate cases.
.. _h03ab-py2-py-references:
**References**
Hadley, G, 1962, `Linear Programming`, Addison--Wesley
"""
raise NotImplementedError
[docs]def shortestpath(n, ns, ne, direct, d, irow, icol):
r"""
``shortestpath`` finds the shortest path through a directed or undirected acyclic network using Dijkstra's algorithm.
.. _h03ad-py2-py-doc:
For full information please refer to the NAG Library document for h03ad
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h03adf.html
.. _h03ad-py2-py-parameters:
**Parameters**
**n** : int
:math:`n`, the number of vertices.
**ns** : int
:math:`n_s` and :math:`n_e`, the labels of the first and last vertices, respectively, between which the shortest path is sought.
**ne** : int
:math:`n_s` and :math:`n_e`, the labels of the first and last vertices, respectively, between which the shortest path is sought.
**direct** : bool
Indicates whether the network is directed or undirected.
:math:`\mathrm{direct} = \mathbf{True}`
The network is directed.
:math:`\mathrm{direct} = \mathbf{False}`
The network is undirected.
**d** : float, array-like, shape :math:`\left(\textit{nnz}\right)`
The nonzero elements of the distance matrix :math:`D`, ordered by increasing row index and increasing column index within each row. More precisely, :math:`\mathrm{d}[k]` must contain the value of the nonzero element with indices (:math:`\mathrm{irow}[k],\mathrm{icol}[k]`); this is the length of the arc from the vertex with label :math:`\mathrm{irow}[k]` to the vertex with label :math:`\mathrm{icol}[k]`. Elements with the same row and column indices are not allowed. If :math:`\mathrm{direct} = \mathbf{False}`, then only those nonzero elements in the strict upper triangle of :math:`\mathrm{d}` need be supplied since :math:`d_{{ij}} = d_{{ji}}`. (:meth:`sparse.real_gen_sort <naginterfaces.library.sparse.real_gen_sort>` may be used to sort the elements of an arbitrarily ordered matrix into the required form.)
**irow** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
:math:`\mathrm{irow}[k]` and :math:`\mathrm{icol}[k]` must contain the row and column indices, respectively, for the nonzero element stored in :math:`\mathrm{d}[k]`.
**icol** : int, array-like, shape :math:`\left(\textit{nnz}\right)`
:math:`\mathrm{irow}[k]` and :math:`\mathrm{icol}[k]` must contain the row and column indices, respectively, for the nonzero element stored in :math:`\mathrm{d}[k]`.
**Returns**
**splen** : float
The length of the shortest path between the specified vertices :math:`n_s` and :math:`n_e`.
**path** : int, ndarray, shape :math:`\left(\mathrm{n}\right)`
Contains details of the shortest path between the specified vertices :math:`n_s` and :math:`n_e`. More precisely, :math:`\mathrm{ns} = \mathrm{path}[0]→\mathrm{path}[1]→ \cdots →\mathrm{path}[p] = \mathrm{ne}` for some :math:`p\leq n`. The remaining :math:`\left(n-p\right)` elements are set to zero.
.. _h03ad-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{ns} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ne} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{ns}\neq \mathrm{ne}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ne} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{ne}\leq \mathrm{n}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{ns} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{ns}\leq \mathrm{n}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 2`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{direct} = \mathbf{False}`, :math:`1\leq \textit{nnz}\leq \mathrm{n}\times \left(\mathrm{n}-1\right)/2`.
(`errno` :math:`2`)
On entry, :math:`\textit{nnz} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{direct} = \mathbf{True}`, :math:`1\leq \textit{nnz}\leq \mathrm{n}\times \left(\mathrm{n}-1\right)`.
(`errno` :math:`3`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{irow}[k-1]\leq \mathrm{n}`, :math:`1\leq \mathrm{icol}[k-1]\leq \mathrm{n}`; :math:`\mathrm{icol}[k-1]\neq \mathrm{irow}[k-1]` when :math:`\mathrm{direct} = \mathbf{True}`.
(`errno` :math:`4`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{irow}[k-1] < \mathrm{icol}[k-1]\leq \mathrm{n}` when :math:`\mathrm{direct} = \mathbf{False}`.
(`errno` :math:`5`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{d}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{d}[k-1] > 0.0`.
(`errno` :math:`6`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[k-2] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[k-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraints: :math:`\mathrm{irow}[k-2] < \mathrm{irow}[k-1]` or :math:`\mathrm{irow}[k-2] = \mathrm{irow}[k-1]` and :math:`\mathrm{icol}[k-2] < \mathrm{icol}[k-1]`.
(`errno` :math:`7`)
On entry, :math:`k = \langle\mathit{\boldsymbol{value}}\rangle` :math:`\mathrm{irow}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{irow}[k-2] = \langle\mathit{\boldsymbol{value}}\rangle` :math:`\mathrm{icol}[k-1] = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{icol}[k-2] = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irow}[k-1]\neq \mathrm{irow}[k-2]` or :math:`\mathrm{icol}[k-1]\neq \mathrm{icol}[k-2]`.
(`errno` :math:`8`)
On entry, :math:`\mathrm{ns} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{ne} = \langle\mathit{\boldsymbol{value}}\rangle`.
No connected network exists between vertices :math:`\mathrm{ns}` and :math:`\mathrm{ne}`.
.. _h03ad-py2-py-notes:
**Notes**
`No equivalent traditional C interface for this routine exists in the NAG Library.`
``shortestpath`` attempts to find the shortest path through a **directed** or **undirected** **acyclic** network, which consists of a set of points called **vertices** and a set of curves called **arcs** that connect certain pairs of distinct vertices.
An acyclic network is one in which there are no paths connecting a vertex to itself.
An arc whose origin vertex is :math:`i` and whose destination vertex is :math:`j` can be written as :math:`i→j`.
In an undirected network the arcs :math:`i→j` and :math:`j→i` are equivalent (i.e., :math:`i\leftrightarrow j`), whereas in a directed network they are different.
Note that the shortest path may not be unique and in some cases may not even exist (e.g., if the network is disconnected).
The network is assumed to consist of :math:`n` vertices which are labelled by the integers :math:`1,2,\ldots,n`.
The lengths of the arcs between the vertices are defined by the :math:`n\times n` **distance matrix** :math:`\mathrm{d}`, in which the element :math:`d_{{ij}}` gives the length of the arc :math:`i→j`; :math:`d_{{ij}} = 0` if there is no arc connecting vertices :math:`i` and :math:`j` (as is the case for an acyclic network when :math:`i = j`).
Thus the matrix :math:`D` is usually **sparse**.
For example, if :math:`n = 4` and the network is directed, then
.. math::
\mathrm{d} = \begin{pmatrix}0&d_{12}&d_{13}&d_{14}\\d_{21}&0&d_{23}&d_{24}\\d_{31}&d_{32}&0&d_{34}\\d_{41}&d_{42}&d_{43}&0\end{pmatrix}\text{.}
If the network is undirected, :math:`\mathrm{d}` is **symmetric** since :math:`d_{{ij}} = d_{{ji}}` (i.e., the length of the arc :math:`i→j\equiv \text{}` the length of the arc :math:`j→i`).
The method used by ``shortestpath`` is described in detail in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h03adf.html#fcomments>`__.
.. _h03ad-py2-py-references:
**References**
Dijkstra, E W, 1959, `A note on two problems in connection with graphs`, Numer. Math. (1), 269--271
"""
raise NotImplementedError
[docs]def tsp_simann(dm, bound, targc, statecomm):
r"""
``tsp_simann`` calculates an approximate solution to a symmetric travelling salesman problem using simulated annealing via a configuration free interface.
.. _h03bb-py2-py-doc:
For full information please refer to the NAG Library document for h03bb
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h03bbf.html
.. _h03bb-py2-py-parameters:
**Parameters**
**dm** : float, array-like, shape :math:`\left(\textit{nc}, \textit{nc}\right)`
The distance matrix; each :math:`\mathrm{dm}[\textit{i}-1,\textit{j}-1]` is the effective cost or weight between nodes :math:`\textit{i}` and :math:`\textit{j}`. Only the strictly upper half of the matrix is referenced.
**bound** : float
A lower bound on the solution. If the optimum is unknown set :math:`\mathrm{bound}` to zero or a negative value; the function will then calculate the minimum spanning tree for :math:`\mathrm{dm}` and use this as a lower bound (returned in :math:`\mathrm{alg\_stats}[5]`). If an optimal value for the cost is known then this should be used for the lower bound. A detailed discussion of relaxations for lower bounds, including the minimal spanning tree, can be found in Reinelt (1994).
**targc** : float
A measure of how close an approximation needs to be to the lower bound. The function terminates when a cost is found less than or equal to :math:`\mathrm{bound}+\mathrm{targc}`. This argument is useful when an optimal value for the cost is known and supplied in :math:`\mathrm{bound}`. It may be sufficient to obtain a path that is close enough (in terms of cost) to the optimal path; this allows the algorithm to terminate at that point and avoid further computation in attempting to find a better path.
If :math:`\mathrm{targc} < 0`, :math:`\mathrm{targc} = 0` is assumed.
**statecomm** : dict, RNG communication object, modified in place
RNG communication structure.
This argument must have been initialized by a prior call to :meth:`rand.init_repeat <naginterfaces.library.rand.init_repeat>` or :meth:`rand.init_nonrepeat <naginterfaces.library.rand.init_nonrepeat>`.
**Returns**
**path** : int, ndarray, shape :math:`\left(\textit{nc}\right)`
The best path discovered by the simulation. That is, :math:`\mathrm{path}` contains the city indices in path order.
**cost** : float
The cost or weight of :math:`\mathrm{path}`.
**tmode** : int
The termination mode of the function:
:math:`\mathrm{tmode} = 0`
Optimal solution found, :math:`\mathrm{cost} = \mathrm{bound}`.
:math:`\mathrm{tmode} = 1`
System temperature cooled. The algorithm returns a :math:`\mathrm{path}` and associated :math:`\mathrm{cost}` that does not attain, nor lie within :math:`\mathrm{targc}` of, the :math:`\mathrm{bound}`. This could be a sufficiently good approximation to the optimal :math:`\mathrm{path}`, particularly when :math:`\mathrm{bound}+\mathrm{targc}` lies below the optimal :math:`\mathrm{cost}`.
:math:`\mathrm{tmode} = 2`
Halted by :math:`\mathrm{cost}` falling within the desired :math:`\mathrm{targc}` range of the :math:`\mathrm{bound}`.
:math:`\mathrm{tmode} = 3`
System stalled following lack of improvement.
:math:`\mathrm{tmode} = 4`
Initial search failed to find a single improvement (the solution could be optimal).
**alg_stats** : float, ndarray, shape :math:`\left(6\right)`
An array of metrics collected during the initial search. These could be used as a basis for future optimization. If an exception is raised, the elements of :math:`\mathrm{alg\_stats}` are set to zero; the first five elements are also set to zero in the trival cases :math:`\textit{nc} = 1`, :math:`2` or :math:`3`.
:math:`\mathrm{alg\_stats}[0]`
Mean delta.
:math:`\mathrm{alg\_stats}[1]`
Standard deviation of deltas.
:math:`\mathrm{alg\_stats}[2]`
Cost at end of initial search phase.
:math:`\mathrm{alg\_stats}[3]`
Best cost encountered during search phase.
:math:`\mathrm{alg\_stats}[4]`
Initial system temperature. At the end of stage 1 of the algorithm, this is a function of the mean and variance of the deltas, and of the distance from best cost to the lower bound. It is a measure of the initial acceptance criteria for stage :math:`2`. The larger this value, the more iterations it will take to geometrically reduce it during stage 2 until the system is cooled (below a threshold value).
:math:`\mathrm{alg\_stats}[5]`
The lower bound used, which will be that computed internally when :math:`\mathrm{bound}\leq 0` on input. Subsequent calls with different random states can set :math:`\mathrm{bound}` to the value returned in :math:`\mathrm{alg\_stats}[5]` to avoid recomputation of the minimal spanning tree.
.. _h03bb-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\textit{nc} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nc} \geq 1`.
(`errno` :math:`2`)
On entry, the strictly upper triangle of :math:`\mathrm{dm}` had a negative element.
(`errno` :math:`9`)
On entry, :math:`\mathrm{statecomm}`\ ['state'] vector has been corrupted or not initialized.
.. _h03bb-py2-py-notes:
**Notes**
``tsp_simann`` provides a probabilistic strategy for the calculation of a near optimal path through a symmetric and fully connected distance matrix; that is, a matrix for which element :math:`\left(i, j\right)` is the pairwise distance (also called the cost, or weight) between nodes (cities) :math:`i` and :math:`j`.
This problem is better known as the Travelling Salesman Problem (TSP), and symmetric means that the distance to travel between two cities is independent of which is the destination city.
In the classical TSP, which this function addresses, a salesman wishes to visit a given set of cities once only by starting and finishing in a home city and travelling the minimum total distance possible.
It is one of the most intensively studied problems in computational mathematics and, as a result, has developed some fairly sophisticated techniques for getting near-optimal solutions for large numbers of cities. ``tsp_simann`` adopts a very simple approach to try to find a reasonable solution, for moderately large problems.
The function uses simulated annealing: a stochastic mechanical process in which the heating and controlled cooling of a material is used to optimally refine its molecular structure.
The material in the TSP is the distance matrix and a given state is represented by the order in which each city is visited---the path.
This system can move from one state to a neighbouring state by selecting two cities on the current path at random and switching their places; the order of the cities in the path between the switched cities is then reversed.
The cost of a state is the total cost of traversing its path; the resulting difference in cost between the current state and this new proposed state is called the delta; a negative delta indicates the proposal creates a more optimal path and a positive delta a less optimal path.
The random selection of cities to switch uses random number generators (RNGs) from submodule :mod:`~naginterfaces.library.rand`; it is thus necessary to initialize a state array for the RNG of choice (by a call to :meth:`rand.init_repeat <naginterfaces.library.rand.init_repeat>` or :meth:`rand.init_nonrepeat <naginterfaces.library.rand.init_nonrepeat>`) prior to calling ``tsp_simann``.
The simulation itself is executed in two stages.
In the first stage, a series of sample searches through the distance matrix is conducted where each proposed new state is accepted, regardless of the change in cost (delta) incurred by applying the switches, and statistics on the set of deltas are recorded.
These metrics are updated after each such sample search; the number of these searches and the number of switches applied in each search is dependent on the number of cities.
The final collated set of metrics for the deltas obtained by the first stage are used as control parameters for the second stage.
If no single improvement in cost is found during the first stage, the algorithm is terminated.
In the second stage, as before, neighbouring states are proposed.
If the resulting delta is negative or causes no change the proposal is accepted and the path updated; otherwise moves are accepted based on a probabilistic criterion, a modified version of the Metropolis--Hastings algorithm.
The acceptance of some positive deltas (increased cost) reduces the probability of a solution getting trapped at a non-optimal solution where any single switch causes an increase in cost.
Initially the acceptance criteria allow for relatively large positive deltas, but as the number of proposed changes increases, the criteria become more stringent, allowing fewer positive deltas of smaller size to be accepted; this process is, within the realm of the simulated annealing algorithm, referred to as 'cooling'.
Further exploration of the system is initially encouraged by accepting non-optimal routes, but is increasingly discouraged as the process continues.
The second stage will terminate when:
- a solution is obtained that is deemed acceptable (as defined by supplied values);
- the algorithm will accept no further positive deltas and a set of proposed changes have resulted in no improvements (has cooled);
- a number of consecutive sets of proposed changes has resulted in no improvement.
.. _h03bb-py2-py-references:
**References**
Applegate, D L, Bixby, R E, Chvátal, V and Cook, W J, 2006, `The Traveling Salesman Problem: A Computational Study`, Princeton University Press
Cook, W J, 2012, `In Pursuit of the Traveling Salesman`, Princeton University Press
Johnson, D S and McGeoch, L A, `The traveling salesman problem: A case study in local optimization`, Local search in combinatorial optimization, 1997, 215--310
Press, W H, Teukolsky, S A, Vetterling, W T and Flannery, B P, 2007, `Numerical Recipes`, The Art of Scientific Computing ((3rd Edition))
Rego, C, Gamboa, D, Glover, F and Osterman, C, 2011, `Traveling salesman problem heuristics: leading methods, implementations and latest advances`, European Journal of Operational Research (211 (3)), 427--441
Reinelt, G, 1994, `The Travelling Salesman. Computational Solutions for TSP Applications, Volume 840 of Lecture Notes in Computer Science`, Springer--Verlag, Berlin Heidelberg New York
"""
raise NotImplementedError
[docs]def best_subset_given_size_revcomm(irevcm, mincr, m, ip, nbest, drop, lz, z, la, a, bscore, bz, mincnt, gamma, acc, comm, io_manager=None):
r"""
Given a set of :math:`m` features and a scoring mechanism for any subset of those features, ``best_subset_given_size_revcomm`` selects the best :math:`n` subsets of size :math:`p` using a reverse communication branch and bound algorithm.
.. _h05aa-py2-py-doc:
For full information please refer to the NAG Library document for h05aa
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05aaf.html
.. _h05aa-py2-py-parameters:
**Parameters**
**irevcm** : int
`On initial entry`: must be set to :math:`0`.
`On intermediate entry`: :math:`\mathrm{irevcm}` must remain unchanged.
**mincr** : int
Flag indicating whether the scoring function :math:`f` is increasing or decreasing.
:math:`\mathrm{mincr} = 1`
:math:`f\left(S_i\right)\leq f\left(S_j\right)`, i.e., the subsets with the largest score will be selected.
:math:`\mathrm{mincr} = 0`
:math:`f\left(S_i\right)\geq f\left(S_j\right)`, i.e., the subsets with the smallest score will be selected.
For all :math:`S_j⊆\Omega` and :math:`S_i⊆S_j`.
**m** : int
:math:`m`, the number of features in the full feature set.
**ip** : int
:math:`p`, the number of features in the subset of interest.
**nbest** : int
:math:`n`, the maximum number of best subsets required. The actual number of subsets returned is given by :math:`\mathrm{la}` on final exit. If on final exit :math:`\mathrm{la}\neq \mathrm{nbest}` then :math:`\mathrm{errno}` = 53 is returned.
**drop** : int
`On initial entry`: :math:`\mathrm{drop}` need not be set.
`On intermediate entry`: :math:`\mathrm{drop}` must remain unchanged.
**lz** : int
`On initial entry`: :math:`\mathrm{lz}` need not be set.
`On intermediate entry`: :math:`\mathrm{lz}` must remain unchanged.
**z** : int, ndarray, shape :math:`\left(\mathrm{m}-\mathrm{ip}\right)`, modified in place
`On initial entry`: :math:`\mathrm{z}` need not be set.
`On intermediate exit`: :math:`\mathrm{z}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{lz}`, contains the list of features which, along with those specified in :math:`\mathrm{a}`, define the subsets whose score is required. See :math:`\mathrm{irevcm}` for additional details.
`On intermediate entry`: :math:`\mathrm{z}` must remain unchanged.
`On final exit`: :math:`\mathrm{z}` is undefined.
**la** : int
`On initial entry`: :math:`\mathrm{la}` need not be set.
`On intermediate entry`: :math:`\mathrm{la}` must remain unchanged.
**a** : int, ndarray, shape :math:`\left(\max\left(\mathrm{nbest},\mathrm{m}\right)\right)`, modified in place
`On initial entry`: :math:`\mathrm{a}` need not be set.
`On intermediate exit`: :math:`\mathrm{a}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{la}`, contains the list of features which, along with those specified in :math:`\mathrm{z}`, define the subsets whose score is required. See :math:`\mathrm{irevcm}` for additional details.
`On intermediate entry`: :math:`\mathrm{a}` must remain unchanged.
`On final exit`: :math:`\mathrm{a}` is undefined.
**bscore** : float, ndarray, shape :math:`\left(\max\left(\mathrm{nbest},\mathrm{m}\right)\right)`, modified in place
`On initial entry`: :math:`\mathrm{bscore}` need not be set.
`On intermediate exit`: :math:`\mathrm{bscore}` is undefined.
`On intermediate entry`: :math:`\mathrm{bscore}[j-1]` must hold the score for the :math:`j`\ th subset as described in :math:`\mathrm{irevcm}`.
`On final exit`: holds the score for the :math:`\mathrm{la}` best subsets returned in :math:`\mathrm{bz}`.
**bz** : int, ndarray, shape :math:`\left(\mathrm{m}-\mathrm{ip}, \mathrm{nbest}\right)`, modified in place
`On initial entry`: :math:`\mathrm{bz}` need not be set.
`On intermediate exit`: :math:`\mathrm{bz}` is used for storage between calls to ``best_subset_given_size_revcomm``.
`On intermediate entry`: :math:`\mathrm{bz}` must remain unchanged.
`On final exit`: the :math:`j`\ th best subset is constructed by dropping the features specified in :math:`\mathrm{bz}[\textit{i}-1,\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{la}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{m}-\mathrm{ip}`, from the set of all features, :math:`\Omega`. The score for the :math:`j`\ th best subset is given in :math:`\mathrm{bscore}[j-1]`.
**mincnt** : int
:math:`k`, the minimum number of times the effect of each feature, :math:`x_i`, must have been observed before :math:`f\left(S-\left\{x_i\right\}\right)` is estimated from :math:`f\left(S\right)` as opposed to being calculated directly.
If :math:`k = 0` then :math:`f\left(S-\left\{x_i\right\}\right)` is never estimated.
If :math:`\mathrm{mincnt} < 0` then :math:`k` is set to :math:`1`.
**gamma** : float
:math:`\gamma`, the scaling factor used when estimating scores. If :math:`\mathrm{gamma} < 0` then :math:`\gamma = 1` is used.
**acc** : float, array-like, shape :math:`\left(2\right)`
A measure of the accuracy of the scoring function, :math:`f`.
Letting :math:`a_i = \epsilon_1\left\lvert f\left(S_i\right)\right\rvert +\epsilon_2`, then when confirming whether the scoring function is strictly increasing or decreasing (as described in :math:`\mathrm{mincr}`), or when assessing whether a node defined by subset :math:`S_i` can be trimmed, then any values in the range :math:`f\left(S_i\right)\pm a_i` are treated as being numerically equivalent.
If :math:`0\leq \mathrm{acc}[0]\leq 1` then :math:`\epsilon_1 = \mathrm{acc}[0]`, otherwise :math:`\epsilon_1 = 0`.
If :math:`\mathrm{acc}[1]\geq 0` then :math:`\epsilon_2 = \mathrm{acc}[1]`, otherwise :math:`\epsilon_2 = 0`.
In most situations setting both :math:`\epsilon_1` and :math:`\epsilon_2` to zero should be sufficient.
Using a nonzero value, when one is not required, can significantly increase the number of subsets that need to be evaluated.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**irevcm** : int
`On intermediate exit`: :math:`\mathrm{irevcm} = 1` and before re-entry the scores associated with :math:`\mathrm{la}` subsets must be calculated and returned in :math:`\mathrm{bscore}`.
The :math:`\mathrm{la}` subsets are constructed as follows:
:math:`\mathrm{drop} = 1`
The :math:`j`\ th subset is constructed by `dropping` the features specified in the first :math:`\mathrm{lz}` elements of :math:`\mathrm{z}` and the single feature given in :math:`\mathrm{a}[j-1]` from the full set of features, :math:`\Omega`. The subset will, therefore, contain :math:`\mathrm{m}-\mathrm{lz}-1` features.
:math:`\mathrm{drop} = 0`
The :math:`j`\ th subset is constructed by `adding` the features specified in the first :math:`\mathrm{lz}` elements of :math:`\mathrm{z}` and the single feature specified in :math:`\mathrm{a}[j-1]` to the empty set, :math:`∅`. The subset will, therefore, contain :math:`\mathrm{lz}+1` features.
In both cases the individual features are referenced by the integers :math:`1` to :math:`\mathrm{m}` with :math:`1` indicating the first feature, :math:`2` the second, etc., for some arbitrary ordering of the features.
The same ordering must be used in all calls to ``best_subset_given_size_revcomm``.
If :math:`\mathrm{la} = 0`, the score for a single subset should be returned.
This subset is constructed by adding or removing only those features specified in the first :math:`\mathrm{lz}` elements of :math:`\mathrm{z}`.
If :math:`\mathrm{lz} = 0`, this subset will either be :math:`\Omega` or :math:`∅`.
The score associated with the :math:`j`\ th subset must be returned in :math:`\mathrm{bscore}[j-1]`.
`On final exit`: :math:`\mathrm{irevcm} = 0`, and the algorithm has terminated.
**drop** : int
`On intermediate exit`: flag indicating whether the intermediate subsets should be constructed by dropping features from the full set (:math:`\mathrm{drop} = 1`) or adding features to the empty set (:math:`\mathrm{drop} = 0`). See :math:`\mathrm{irevcm}` for details.
`On final exit`: :math:`\mathrm{drop}` is undefined.
**lz** : int
`On intermediate exit`: the number of features stored in :math:`\mathrm{z}`.
`On final exit`: :math:`\mathrm{lz}` is undefined.
**la** : int
`On intermediate exit`: if :math:`\mathrm{la} > 0`, the number of subsets for which a score must be returned.
If :math:`\mathrm{la} = 0`, the score for a single subset should be returned.
See :math:`\mathrm{irevcm}` for additional details.
`On final exit`: the number of best subsets returned.
.. _h05aa-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, :math:`\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{irevcm} = 0` or :math:`1`.
(`errno` :math:`21`)
On entry, :math:`\mathrm{mincr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mincr} = 0` or :math:`1`.
(`errno` :math:`22`)
:math:`\mathrm{mincr}` has changed between calls.
On intermediate entry, :math:`\mathrm{mincr} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{mincr} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`31`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 2`.
(`errno` :math:`32`)
:math:`\mathrm{m}` has changed between calls.
On intermediate entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`41`)
On entry, :math:`\mathrm{ip} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{ip}\leq \mathrm{m}`.
(`errno` :math:`42`)
:math:`\mathrm{ip}` has changed between calls.
On intermediate entry, :math:`\mathrm{ip} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{ip} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`51`)
On entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nbest}\geq 1`.
(`errno` :math:`52`)
:math:`\mathrm{nbest}` has changed between calls.
On intermediate entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`61`)
:math:`\mathrm{drop}` has changed between calls.
On intermediate entry, :math:`\mathrm{drop} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{drop} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`71`)
:math:`\mathrm{lz}` has changed between calls.
On entry, :math:`\mathrm{lz} = \langle\mathit{\boldsymbol{value}}\rangle`.
On previous exit, :math:`\mathrm{lz} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`91`)
:math:`\mathrm{la}` has changed between calls.
On entry, :math:`\mathrm{la} = \langle\mathit{\boldsymbol{value}}\rangle`.
On previous exit, :math:`\mathrm{la} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`111`)
:math:`\mathrm{bscore}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle`, which is inconsistent with the score for the parent node. Score for the parent node is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`131`)
:math:`\mathrm{mincnt}` has changed between calls.
On intermediate entry, :math:`\mathrm{mincnt} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{mincnt} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`141`)
:math:`\mathrm{gamma}` has changed between calls.
On intermediate entry, :math:`\mathrm{gamma} = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{gamma} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`151`)
:math:`\mathrm{acc}[0]` has changed between calls.
On intermediate entry, :math:`\mathrm{acc}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{acc}[0] = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`152`)
:math:`\mathrm{acc}[1]` has changed between calls.
On intermediate entry, :math:`\mathrm{acc}[1] = \langle\mathit{\boldsymbol{value}}\rangle`.
On initial entry, :math:`\mathrm{acc}[1] = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`161`)
:math:`\mathrm{comm}`\ ['icomm'] has been corrupted between calls.
(`errno` :math:`181`)
:math:`\mathrm{comm}`\ ['rcomm'] has been corrupted between calls.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`53`)
On entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
But only :math:`\langle\mathit{\boldsymbol{value}}\rangle` best subsets could be calculated.
.. _h05aa-py2-py-notes:
**Notes**
Given :math:`\Omega = \left\{x_i:i \in ℤ,1\leq i\leq m\right\}`, a set of :math:`m` unique features and a scoring mechanism :math:`f\left(S\right)` defined for all :math:`S⊆\Omega` then ``best_subset_given_size_revcomm`` is designed to find :math:`S_{{o1}}⊆\Omega,\left|S_{{o1}}\right| = p`, an optimal subset of size :math:`p`.
Here :math:`\left|S_{{o1}}\right|` denotes the cardinality of :math:`S_{{o1}}`, the number of elements in the set.
The definition of the optimal subset depends on the properties of the scoring mechanism, if
.. math::
\begin{array}{cc} f\left(S_i\right) \leq f\left(S_j\right) \text{,} & \text{for all } S_j ⊆ \Omega \text{ and } S_i ⊆ S_j \end{array}
then the optimal subset is defined as one of the solutions to
.. math::
\textit{maximize}_{{S⊆\Omega }}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
else if
.. math::
\begin{array}{cc} f\left(S_i\right) \geq f\left(S_j\right) \text{,} & \text{for all } S_j ⊆ \Omega \text{ and } S_i ⊆ S_j \end{array}
then the optimal subset is defined as one of the solutions to
.. math::
\textit{minimize}_{{S⊆\Omega }}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p\text{.}
If neither of these properties hold then ``best_subset_given_size_revcomm`` cannot be used.
As well as returning the optimal subset, :math:`S_{{o1}}`, ``best_subset_given_size_revcomm`` can return the best :math:`n` solutions of size :math:`p`.
If :math:`S_{{o\textit{i}}}` denotes the :math:`\textit{i}`\ th best subset, for :math:`\textit{i} = 1,2,\ldots,n-1`, then the :math:`\left(i+1\right)`\ th best subset is defined as the solution to either
.. math::
\textit{maximize}_{{S⊆\Omega -\left\{S_{{oj}}:j \in ℤ,1\leq j\leq i\right\}}}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
or
.. math::
\textit{minimize}_{{S⊆\Omega -\left\{S_{{oj}}:j \in ℤ,1\leq j\leq i\right\}}}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
depending on the properties of :math:`f`.
The solutions are found using a branch and bound method, where each node of the tree is a subset of :math:`\Omega`.
Assuming that `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05aaf.html#score_property_eqn1>`__ holds then a particular node, defined by subset :math:`S_i`, can be trimmed from the tree if :math:`f\left(S_i\right) < \hat{f}\left(S_{{on}}\right)` where :math:`\hat{f}\left(S_{{on}}\right)` is the :math:`n`\ th highest score we have observed so far for a subset of size :math:`p`, i.e., our current best guess of the score for the :math:`n`\ th best subset.
In addition, because of `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05aaf.html#score_property_eqn1>`__ we can also drop all nodes defined by any subset :math:`S_j` where :math:`S_j⊆S_i`, thus avoiding the need to enumerate the whole tree.
Similar short cuts can be taken if `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05aaf.html#score_property_eqn2>`__ holds.
A full description of this branch and bound algorithm can be found in Ridout (1988).
Rather than calculate the score at a given node of the tree ``best_subset_given_size_revcomm`` utilizes the fast branch and bound algorithm of Somol `et al.` (2004), and attempts to estimate the score where possible.
For each feature, :math:`x_i`, two values are stored, a count :math:`c_i` and :math:`\hat{\mu }_i`, an estimate of the contribution of that feature.
An initial value of zero is used for both :math:`c_i` and :math:`\hat{\mu }_i`.
At any stage of the algorithm where both :math:`f\left(S\right)` and :math:`f\left(S-\left\{x_i\right\}\right)` have been calculated (as opposed to estimated), the estimated contribution of the feature :math:`x_i` is updated to
.. math::
\frac{{c_i\hat{\mu }_i+\left[f\left(S\right)-f\left(S-\left\{x_j\right\}\right)\right]}}{{c_i+1}}
and :math:`c_i` is incremented by :math:`1`, therefore, at each stage :math:`\hat{\mu }_i` is the mean contribution of :math:`x_i` observed so far and :math:`c_i` is the number of observations used to calculate that mean.
As long as :math:`c_i\geq k`, for the user-supplied constant :math:`k`, then rather than calculating :math:`f\left(S-\left\{x_i\right\}\right)` this function estimates it using :math:`\hat{f}\left(S-\left\{x_i\right\}\right) = f\left(S\right)-\gamma \hat{\mu }_i` or :math:`\hat{f}\left(S\right)-\gamma \hat{\mu }_i` if :math:`f\left(S\right)` has been estimated, where :math:`\gamma` is a user-supplied scaling factor.
An estimated score is never used to trim a node or returned as the optimal score.
Setting :math:`k = 0` in this function will cause the algorithm to always calculate the scores, returning to the branch and bound algorithm of Ridout (1988).
In most cases it is preferable to use the fast branch and bound algorithm, by setting :math:`k > 0`, unless the score function is iterative in nature, i.e., :math:`f\left(S\right)` must have been calculated before :math:`f\left(S-\left\{x_i\right\}\right)` can be calculated.
.. _h05aa-py2-py-references:
**References**
Narendra, P M and Fukunaga, K, 1977, `A branch and bound algorithm for feature subset selection`, IEEE Transactions on Computers (9), 917--922
Ridout, M S, 1988, `Algorithm AS 233: An improved branch and bound algorithm for feature subset selection`, Journal of the Royal Statistics Society, Series C (Applied Statistics) (Volume 37) (1), 139--147
Somol, P, Pudil, P and Kittler, J, 2004, `Fast branch and bound algorithms for optimal feature selection`, IEEE Transactions on Pattern Analysis and Machine Intelligence (Volume 26) (7), 900--912
"""
raise NotImplementedError
[docs]def best_subset_given_size(mincr, m, ip, nbest, f, mincnt, gamma, acc, data=None, io_manager=None):
r"""
Given a set of :math:`m` features and a scoring mechanism for any subset of those features, ``best_subset_given_size`` selects the best :math:`n` subsets of size :math:`p` using a direct communication branch and bound algorithm.
.. _h05ab-py2-py-doc:
For full information please refer to the NAG Library document for h05ab
https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05abf.html
.. _h05ab-py2-py-parameters:
**Parameters**
**mincr** : int
Flag indicating whether the scoring function :math:`f` is increasing or decreasing.
:math:`\mathrm{mincr} = 1`
:math:`f\left(S_i\right)\leq f\left(S_j\right)`, i.e., the subsets with the largest score will be selected.
:math:`\mathrm{mincr} = 0`
:math:`f\left(S_i\right)\geq f\left(S_j\right)`, i.e., the subsets with the smallest score will be selected.
For all :math:`S_j⊆\Omega` and :math:`S_i⊆S_j`.
**m** : int
:math:`m`, the number of features in the full feature set.
**ip** : int
:math:`p`, the number of features in the subset of interest.
**nbest** : int
:math:`n`, the maximum number of best subsets required. The actual number of subsets returned is given by :math:`\mathrm{la}` on final exit. If on final exit :math:`\mathrm{la}\neq \mathrm{nbest}` then :math:`\mathrm{errno}` = 42 is returned.
**f** : callable score = f(m, drop, z, a, data=None)
:math:`\mathrm{f}` must evaluate the scoring function :math:`f`.
**Parameters**
**m** : int
:math:`m = \left|\Omega \right|`, the number of features in the full feature set.
**drop** : int
Flag indicating whether the intermediate subsets should be constructed by dropping features from the full set (:math:`\mathrm{drop} = 1`) or adding features to the empty set (:math:`\mathrm{drop} = 0`). See :math:`\mathrm{score}` for additional details.
**z** : int, ndarray, shape :math:`\left(\textit{lz}\right)`
:math:`\mathrm{z}[\textit{i}-1]`, for :math:`\textit{i} = 1,2,\ldots,\textit{lz}`, contains the list of features which, along with those specified in :math:`\mathrm{a}`, define the subsets whose score is required. See :math:`\mathrm{score}` for additional details.
**a** : int, ndarray, shape :math:`\left(\textit{la}\right)`
:math:`\mathrm{a}[\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{la}`, contains the list of features which, along with those specified in :math:`\mathrm{z}`, define the subsets whose score is required. See :math:`\mathrm{score}` for additional details.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**score** : float, array-like, shape :math:`\left(\max\left(\mathrm{la},1\right)\right)`
The value :math:`f\left(S_{\textit{j}}\right)`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{la}`, the score associated with the :math:`j`\ th subset. :math:`S_j` is constructed as follows:
:math:`\mathrm{drop} = 1`
:math:`S_j` is constructed by `dropping` the features specified in the first :math:`\textit{lz}` elements of :math:`\mathrm{z}` and the single feature given in :math:`\mathrm{a}[j-1]` from the full set of features, :math:`\Omega`. The subset will, therefore, contain :math:`\mathrm{m}-\textit{lz}-1` features.
:math:`\mathrm{drop} = 0`
:math:`S_j` is constructed by `adding` the features specified in the first :math:`\textit{lz}` elements of :math:`\mathrm{z}` and the single feature specified in :math:`\mathrm{a}[j-1]` to the empty set, :math:`∅`. The subset will, therefore, contain :math:`\textit{lz}+1` features.
In both cases the individual features are referenced by the integers :math:`1` to :math:`\mathrm{m}` with :math:`1` indicating the first feature, :math:`2` the second, etc., for some arbitrary ordering of the features, chosen by you prior to calling ``best_subset_given_size``.
For example, :math:`1` might refer to the first variable in a particular set of data, :math:`2` the second, etc..
If :math:`\mathrm{la} = 0`, the score for a single subset should be returned.
This subset is constructed by adding or removing only those features specified in the first :math:`\textit{lz}` elements of :math:`\mathrm{z}`.
If :math:`\textit{lz} = 0`, this subset will either be :math:`\Omega` or :math:`∅`.
**mincnt** : int
:math:`k`, the minimum number of times the effect of each feature, :math:`x_i`, must have been observed before :math:`f\left(S-\left\{x_i\right\}\right)` is estimated from :math:`f\left(S\right)` as opposed to being calculated directly.
If :math:`k = 0` then :math:`f\left(S-\left\{x_i\right\}\right)` is never estimated.
If :math:`\mathrm{mincnt} < 0` then :math:`k` is set to :math:`1`.
**gamma** : float
:math:`\gamma`, the scaling factor used when estimating scores. If :math:`\mathrm{gamma} < 0` then :math:`\gamma = 1` is used.
**acc** : float, array-like, shape :math:`\left(2\right)`
A measure of the accuracy of the scoring function, :math:`f`.
Letting :math:`a_i = \epsilon_1\left\lvert f\left(S_i\right)\right\rvert +\epsilon_2`, then when confirming whether the scoring function is strictly increasing or decreasing (as described in :math:`\mathrm{mincr}`), or when assessing whether a node defined by subset :math:`S_i` can be trimmed, then any values in the range :math:`f\left(S_i\right)\pm a_i` are treated as being numerically equivalent.
If :math:`0\leq \mathrm{acc}[0]\leq 1` then :math:`\epsilon_1 = \mathrm{acc}[0]`, otherwise :math:`\epsilon_1 = 0`.
If :math:`\mathrm{acc}[1]\geq 0` then :math:`\epsilon_2 = \mathrm{acc}[1]`, otherwise :math:`\epsilon_2 = 0`.
In most situations setting both :math:`\epsilon_1` and :math:`\epsilon_2` to zero should be sufficient.
Using a nonzero value, when one is not required, can significantly increase the number of subsets that need to be evaluated.
**data** : arbitrary, optional
User-communication data for callback functions.
**io_manager** : FileObjManager, optional
Manager for I/O in this routine.
**Returns**
**la** : int
The number of best subsets returned.
**bscore** : float, ndarray, shape :math:`\left(\mathrm{nbest}\right)`
Holds the score for the :math:`\mathrm{la}` best subsets returned in :math:`\mathrm{bz}`.
**bz** : int, ndarray, shape :math:`\left(\mathrm{m}-\mathrm{ip}, \mathrm{nbest}\right)`
The :math:`j`\ th best subset is constructed by dropping the features specified in :math:`\mathrm{bz}[\textit{i}-1,\textit{j}-1]`, for :math:`\textit{j} = 1,2,\ldots,\mathrm{la}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{m}-\mathrm{ip}`, from the set of all features, :math:`\Omega`. The score for the :math:`j`\ th best subset is given in :math:`\mathrm{bscore}[j-1]`.
.. _h05ab-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`11`)
On entry, :math:`\mathrm{mincr} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{mincr} = 0` or :math:`1`.
(`errno` :math:`21`)
On entry, :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{m}\geq 2`.
(`errno` :math:`31`)
On entry, :math:`\mathrm{ip} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`1\leq \mathrm{ip}\leq \mathrm{m}`.
(`errno` :math:`41`)
On entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nbest}\geq 1`.
(`errno` :math:`81`)
On exit from :math:`\mathrm{f}`, :math:`\mathrm{score}[\langle\mathit{\boldsymbol{value}}\rangle] = \langle\mathit{\boldsymbol{value}}\rangle`, which is inconsistent with the score for the parent node. Score for the parent node is :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`42`)
On entry, :math:`\mathrm{nbest} = \langle\mathit{\boldsymbol{value}}\rangle`.
But only :math:`\langle\mathit{\boldsymbol{value}}\rangle` best subsets could be calculated.
**NagCallbackTerminateWarning**
(`errno` :math:`82`)
A nonzero value for :math:`\mathrm{info}` has been returned: :math:`\mathrm{info} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _h05ab-py2-py-notes:
**Notes**
Given :math:`\Omega = \left\{x_i:i \in ℤ,1\leq i\leq m\right\}`, a set of :math:`m` unique features and a scoring mechanism :math:`f\left(S\right)` defined for all :math:`S⊆\Omega` then ``best_subset_given_size`` is designed to find :math:`S_{{o1}}⊆\Omega,\left|S_{{o1}}\right| = p`, an optimal subset of size :math:`p`.
Here :math:`\left|S_{{o1}}\right|` denotes the cardinality of :math:`S_{{o1}}`, the number of elements in the set.
The definition of the optimal subset depends on the properties of the scoring mechanism, if
.. math::
\begin{array}{cc} f\left(S_i\right) \leq f\left(S_j\right) \text{,} & \text{for all } S_j ⊆ \Omega \text{ and } S_i ⊆ S_j \end{array}
then the optimal subset is defined as one of the solutions to
.. math::
\textit{maximize}_{{S⊆\Omega }}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
else if
.. math::
\begin{array}{cc} f\left(S_i\right) \geq f\left(S_j\right) \text{,} & \text{for all } S_j ⊆ \Omega \text{ and } S_i ⊆ S_j \end{array}
then the optimal subset is defined as one of the solutions to
.. math::
\textit{minimize}_{{S⊆\Omega }}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p\text{.}
If neither of these properties hold then ``best_subset_given_size`` cannot be used.
As well as returning the optimal subset, :math:`S_{{o1}}`, ``best_subset_given_size`` can return the best :math:`n` solutions of size :math:`p`.
If :math:`S_{{o\textit{i}}}` denotes the :math:`\textit{i}`\ th best subset, for :math:`\textit{i} = 1,2,\ldots,n-1`, then the :math:`\left(i+1\right)`\ th best subset is defined as the solution to either
.. math::
\textit{maximize}_{{S⊆\Omega -\left\{S_{{oj}}:j \in ℤ,1\leq j\leq i\right\}}}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
or
.. math::
\textit{minimize}_{{S⊆\Omega -\left\{S_{{oj}}:j \in ℤ,1\leq j\leq i\right\}}}f\left(S\right)\quad \text{ subject to }\quad \left|S\right| = p
depending on the properties of :math:`f`.
The solutions are found using a branch and bound method, where each node of the tree is a subset of :math:`\Omega`.
Assuming that `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05abf.html#score_property_eqn1>`__ holds then a particular node, defined by subset :math:`S_i`, can be trimmed from the tree if :math:`f\left(S_i\right) < \hat{f}\left(S_{{on}}\right)` where :math:`\hat{f}\left(S_{{on}}\right)` is the :math:`n`\ th highest score we have observed so far for a subset of size :math:`p`, i.e., our current best guess of the score for the :math:`n`\ th best subset.
In addition, because of `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05abf.html#score_property_eqn1>`__ we can also drop all nodes defined by any subset :math:`S_j` where :math:`S_j⊆S_i`, thus avoiding the need to enumerate the whole tree.
Similar short cuts can be taken if `[equation] <https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/h/h05abf.html#score_property_eqn2>`__ holds.
A full description of this branch and bound algorithm can be found in Ridout (1988).
Rather than calculate the score at a given node of the tree ``best_subset_given_size`` utilizes the fast branch and bound algorithm of Somol `et al.` (2004), and attempts to estimate the score where possible.
For each feature, :math:`x_i`, two values are stored, a count :math:`c_i` and :math:`\hat{\mu }_i`, an estimate of the contribution of that feature.
An initial value of zero is used for both :math:`c_i` and :math:`\hat{\mu }_i`.
At any stage of the algorithm where both :math:`f\left(S\right)` and :math:`f\left(S-\left\{x_i\right\}\right)` have been calculated (as opposed to estimated), the estimated contribution of the feature :math:`x_i` is updated to
.. math::
\frac{{c_i\hat{\mu }_i+\left[f\left(S\right)-f\left(S-\left\{x_j\right\}\right)\right]}}{{c_i+1}}
and :math:`c_i` is incremented by :math:`1`, therefore, at each stage :math:`\hat{\mu }_i` is the mean contribution of :math:`x_i` observed so far and :math:`c_i` is the number of observations used to calculate that mean.
As long as :math:`c_i\geq k`, for the user-supplied constant :math:`k`, then rather than calculating :math:`f\left(S-\left\{x_i\right\}\right)` this function estimates it using :math:`\hat{f}\left(S-\left\{x_i\right\}\right) = f\left(S\right)-\gamma \hat{\mu }_i` or :math:`\hat{f}\left(S\right)-\gamma \hat{\mu }_i` if :math:`f\left(S\right)` has been estimated, where :math:`\gamma` is a user-supplied scaling factor.
An estimated score is never used to trim a node or returned as the optimal score.
Setting :math:`k = 0` in this function will cause the algorithm to always calculate the scores, returning to the branch and bound algorithm of Ridout (1988).
In most cases it is preferable to use the fast branch and bound algorithm, by setting :math:`k > 0`, unless the score function is iterative in nature, i.e., :math:`f\left(S\right)` must have been calculated before :math:`f\left(S-\left\{x_i\right\}\right)` can be calculated.
``best_subset_given_size`` is a direct communication version of :meth:`best_subset_given_size_revcomm`.
.. _h05ab-py2-py-references:
**References**
Narendra, P M and Fukunaga, K, 1977, `A branch and bound algorithm for feature subset selection`, IEEE Transactions on Computers (9), 917--922
Ridout, M S, 1988, `Algorithm AS 233: An improved branch and bound algorithm for feature subset selection`, Journal of the Royal Statistics Society, Series C (Applied Statistics) (Volume 37) (1), 139--147
Somol, P, Pudil, P and Kittler, J, 2004, `Fast branch and bound algorithms for optimal feature selection`, IEEE Transactions on Pattern Analysis and Machine Intelligence (Volume 26) (7), 900--912
"""
raise NotImplementedError
