Parma_Polyhedra_Library::MIP_Problem Class Reference
[C++ Language Interface]

A Mixed Integer (linear) Programming problem. More...

#include <MIP_Problem.defs.hh>

Collaboration diagram for Parma_Polyhedra_Library::MIP_Problem:

Collaboration graph
[legend]

List of all members.

Public Types

enum  Control_Parameter_Name { PRICING }
 Names of MIP problems' control parameters. More...
enum  Control_Parameter_Value { PRICING_STEEPEST_EDGE_FLOAT, PRICING_STEEPEST_EDGE_EXACT, PRICING_TEXTBOOK }
 Possible values for MIP problem's control parameters. More...
typedef
Constraint_Sequence::const_iterator 
const_iterator
 A type alias for the read-only iterator on the constraints defining the feasible region.

Public Member Functions

 MIP_Problem (dimension_type dim=0)
 Builds a trivial MIP problem.
template<typename In>
 MIP_Problem (dimension_type dim, In first, In last, const Variables_Set &int_vars, const Linear_Expression &obj=Linear_Expression::zero(), Optimization_Mode mode=MAXIMIZATION)
 Builds an MIP problem having space dimension dim from the sequence of constraints in the range $[\mathrm{first}, \mathrm{last})$, the objective function obj and optimization mode mode; those dimensions whose indices occur in int_vars are constrained to take an integer value.
template<typename In>
 MIP_Problem (dimension_type dim, In first, In last, const Linear_Expression &obj=Linear_Expression::zero(), Optimization_Mode mode=MAXIMIZATION)
 Builds an MIP problem having space dimension dim from the sequence of constraints in the range $[\mathrm{first}, \mathrm{last})$, the objective function obj and optimization mode mode.
 MIP_Problem (dimension_type dim, const Constraint_System &cs, const Linear_Expression &obj=Linear_Expression::zero(), Optimization_Mode mode=MAXIMIZATION)
 Builds an MIP problem having space dimension dim from the constraint system cs, the objective function obj and optimization mode mode.
 MIP_Problem (const MIP_Problem &y)
 Ordinary copy-constructor.
 ~MIP_Problem ()
 Destructor.
MIP_Problemoperator= (const MIP_Problem &y)
 Assignment operator.
dimension_type space_dimension () const
 Returns the space dimension of the MIP problem.
const Variables_Set & integer_space_dimensions () const
 Returns a set containing all the variables' indexes constrained to be integral.
const_iterator constraints_begin () const
 Returns a read-only iterator to the first constraint defining the feasible region.
const_iterator constraints_end () const
 Returns a past-the-end read-only iterator to the sequence of constraints defining the feasible region.
const Linear_Expressionobjective_function () const
 Returns the objective function.
Optimization_Mode optimization_mode () const
 Returns the optimization mode.
void clear ()
 Resets *this to be equal to the trivial MIP problem.
void add_space_dimensions_and_embed (dimension_type m)
 Adds m new space dimensions and embeds the old MIP problem in the new vector space.
void add_to_integer_space_dimensions (const Variables_Set &i_vars)
 Sets the variables whose indexes are in set i_vars to be integer space dimensions.
void add_constraint (const Constraint &c)
 Adds a copy of constraint c to the MIP problem.
void add_constraints (const Constraint_System &cs)
 Adds a copy of the constraints in cs to the MIP problem.
void set_objective_function (const Linear_Expression &obj)
 Sets the objective function to obj.
void set_optimization_mode (Optimization_Mode mode)
 Sets the optimization mode to mode.
bool is_satisfiable () const
 Checks satisfiability of *this.
MIP_Problem_Status solve () const
 Optimizes the MIP problem.
void evaluate_objective_function (const Generator &evaluating_point, Coefficient &num, Coefficient &den) const
 Sets num and den so that $\frac{num}{den}$ is the result of evaluating the objective function on evaluating_point.
const Generatorfeasible_point () const
 Returns a feasible point for *this, if it exists.
const Generatoroptimizing_point () const
 Returns an optimal point for *this, if it exists.
void optimal_value (Coefficient &num, Coefficient &den) const
 Sets num and den so that $\frac{num}{den}$ is the solution of the optimization problem.
bool OK () const
 Checks if all the invariants are satisfied.
void ascii_dump () const
 Writes to std::cerr an ASCII representation of *this.
void ascii_dump (std::ostream &s) const
 Writes to s an ASCII representation of *this.
void print () const
 Prints *this to std::cerr using operator<<.
bool ascii_load (std::istream &s)
 Loads from s an ASCII representation (as produced by ascii_dump(std::ostream&) const) and sets *this accordingly. Returns true if successful, false otherwise.
memory_size_type total_memory_in_bytes () const
 Returns the total size in bytes of the memory occupied by *this.
memory_size_type external_memory_in_bytes () const
 Returns the size in bytes of the memory managed by *this.
void swap (MIP_Problem &y)
 Swaps *this with y.
Control_Parameter_Value get_control_parameter (Control_Parameter_Name name) const
 Returns the value of the control parameter name.
void set_control_parameter (Control_Parameter_Value value)
 Sets control parameter value.

Static Public Member Functions

static dimension_type max_space_dimension ()
 Returns the maximum space dimension an MIP_Problem can handle.

Private Types

enum  Status {
  UNSATISFIABLE, SATISFIABLE, UNBOUNDED, OPTIMIZED,
  PARTIALLY_SATISFIABLE
}
 An enumerated type describing the internal status of the MIP problem. More...
typedef std::vector< ConstraintConstraint_Sequence
 A type alias for a sequence of constraints.

Private Member Functions

bool process_pending_constraints ()
 Processes the pending constraints of *this.
void second_phase ()
 Optimizes the MIP problem using the second phase of the primal simplex algorithm.
MIP_Problem_Status compute_tableau (std::vector< dimension_type > &worked_out_row)
 Assigns to this->tableau a simplex tableau representing the MIP problem, inserting into this->mapping the information that is required to recover the original MIP problem.
bool parse_constraints (dimension_type &new_num_rows, dimension_type &num_slack_variables, std::deque< bool > &is_tableau_constraint, std::deque< bool > &nonnegative_variable, std::vector< dimension_type > &unfeasible_tableau_rows, std::deque< bool > &satisfied_ineqs)
 Parses the pending constraints to gather information on how to resize the tableau.
dimension_type get_exiting_base_index (dimension_type entering_var_index) const
 Computes the row index of the variable exiting the base of the MIP problem. Implemented with anti-cycling rule.
void pivot (dimension_type entering_var_index, dimension_type exiting_base_index)
 Performs the pivoting operation on the tableau.
dimension_type textbook_entering_index () const
 Computes the column index of the variable entering the base, using the textbook algorithm with anti-cycling rule.
dimension_type steepest_edge_exact_entering_index () const
 Computes the column index of the variable entering the base, using an exact steepest-edge algorithm with anti-cycling rule.
dimension_type steepest_edge_float_entering_index () const
 Same as steepest_edge_exact_entering_index, but using floating points.
bool compute_simplex_using_exact_pricing ()
 Returns true if and if only the algorithm successfully computed a feasible solution.
bool compute_simplex_using_steepest_edge_float ()
 Returns true if and if only the algorithm successfully computed a feasible solution.
void erase_artificials (dimension_type begin_artificials, dimension_type end_artificials)
 Drop unnecessary artificial variables from the tableau and get ready for the second phase of the simplex algorithm.
bool is_in_base (dimension_type var_index, dimension_type &row_index) const
void compute_generator () const
 Computes a valid generator that satisfies all the constraints of the Linear Programming problem associated to *this.
void merge_split_variables (dimension_type var_index, std::vector< dimension_type > &nonfeasible_cs)
 Merges previously split variables in the tableau if a nonnegativity constraint is detected.
bool is_lp_satisfiable () const

Static Private Member Functions

static void linear_combine (Row &x, const Row &y, const dimension_type k)
 Linearly combines x with y so that *this[k] is 0.
static bool is_satisfied (const Constraint &c, const Generator &g)
 Returns true if and only if c is satisfied by g.
static bool is_saturated (const Constraint &c, const Generator &g)
 Returns true if and only if c is saturated by g.
static MIP_Problem_Status solve_mip (bool &have_incumbent_solution, mpq_class &incumbent_solution_value, Generator &incumbent_solution_point, MIP_Problem &mip, const Variables_Set &i_vars)
 Returns a status that encodes the solution of the MIP problem.
static bool is_mip_satisfiable (MIP_Problem &mip, Generator &p, const Variables_Set &i_vars)
 Used with MIP_Problems with a non empty `i_vars', returns true if and if only a MIP problem is satisfiable, returns false otherwise.
static bool choose_branching_variable (const MIP_Problem &mip, const Variables_Set &i_vars, dimension_type &branching_index)
 Returns true if and if only `last_generator' satisfies all the integrality coditions.

Private Attributes

dimension_type external_space_dim
 The dimension of the vector space.
dimension_type internal_space_dim
 The space dimension of the current (partial) solution of the MIP problem; it may be smaller than external_space_dim.
Matrix tableau
 The matrix encoding the current feasible region in tableau form.
Row working_cost
 The working cost function.
std::vector< std::pair
< dimension_type,
dimension_type > > 
mapping
 A map between the variables of `input_cs' and `tableau'.
std::vector< dimension_typebase
 The current basic solution.
Status status
 The internal state of the MIP problem.
Control_Parameter_Value pricing
 The pricing method in use.
bool initialized
 A Boolean encoding whether or not internal data structures have already been properly sized and populated: useful to allow for deeper checks in method OK().
Constraint_Sequence input_cs
 The sequence of constraints describing the feasible region.
dimension_type first_pending_constraint
 The first index of `input_cs' containing a pending constraint.
Linear_Expression input_obj_function
 The objective function to be optimized.
Optimization_Mode opt_mode
 The optimization mode requested.
Generator last_generator
 The last successfully computed feasible or optimizing point.
Variables_Set i_variables
 A set containing all the indexes of variables that are constrained to have an integer value.

Related Functions

(Note that these are not member functions.)

std::ostream & operator<< (std::ostream &s, const MIP_Problem &lp)
 Output operator.
void swap (Parma_Polyhedra_Library::MIP_Problem &x, Parma_Polyhedra_Library::MIP_Problem &y)
 Specializes std::swap.


Detailed Description

A Mixed Integer (linear) Programming problem.

An object of this class encodes a mixed integer (linear) programming problem. The MIP problem is specified by providing:

The class provides support for the (incremental) solution of the MIP problem based on variations of the revised simplex method and on branch-and-bound techniques. The result of the resolution process is expressed in terms of an enumeration, encoding the feasibility and the unboundedness of the optimization problem. The class supports simple feasibility tests (i.e., no optimization), as well as the extraction of an optimal (resp., feasible) point, provided the MIP_Problem is optimizable (resp., feasible).

By exploiting the incremental nature of the solver, it is possible to reuse part of the computational work already done when solving variants of a given MIP_Problem: currently, incremental resolution supports the addition of space dimensions, the addition of constraints, the change of objective function and the change of optimization mode.

Definition at line 79 of file MIP_Problem.defs.hh.


Member Typedef Documentation

A type alias for a sequence of constraints.

Definition at line 229 of file MIP_Problem.defs.hh.

typedef Constraint_Sequence::const_iterator Parma_Polyhedra_Library::MIP_Problem::const_iterator

A type alias for the read-only iterator on the constraints defining the feasible region.

Definition at line 236 of file MIP_Problem.defs.hh.


Member Enumeration Documentation

Names of MIP problems' control parameters.

Enumerator:
PRICING  The pricing rule.

Definition at line 403 of file MIP_Problem.defs.hh.

00403                               {
00405     PRICING
00406   };

Possible values for MIP problem's control parameters.

Enumerator:
PRICING_STEEPEST_EDGE_FLOAT  Steepest edge pricing method, using floating points (default).
PRICING_STEEPEST_EDGE_EXACT  Steepest edge pricing method, using Coefficient.
PRICING_TEXTBOOK  Textbook pricing method.

Definition at line 409 of file MIP_Problem.defs.hh.

An enumerated type describing the internal status of the MIP problem.

Enumerator:
UNSATISFIABLE  The MIP problem is unsatisfiable.
SATISFIABLE  The MIP problem is satisfiable; a feasible solution has been computed.
UNBOUNDED  The MIP problem is unbounded; a feasible solution has been computed.
OPTIMIZED  The MIP problem is optimized; an optimal solution has been computed.
PARTIALLY_SATISFIABLE  The feasible region of the MIP problem has been changed by adding new space dimensions or new constraints; a feasible solution for the old feasible region is still available.

Definition at line 455 of file MIP_Problem.defs.hh.

00455               {
00457     UNSATISFIABLE,
00459     SATISFIABLE,
00461     UNBOUNDED,
00463     OPTIMIZED,
00469     PARTIALLY_SATISFIABLE
00470   };


Constructor & Destructor Documentation

Parma_Polyhedra_Library::MIP_Problem::MIP_Problem ( dimension_type  dim = 0  )  [explicit]

Builds a trivial MIP problem.

A trivial MIP problem requires to maximize the objective function $0$ on a vector space under no constraints at all: the origin of the vector space is an optimal solution.

Parameters:
dim The dimension of the vector space enclosing *this (optional argument with default value $0$).
Exceptions:
std::length_error Thrown if dim exceeds max_space_dimension().

Definition at line 62 of file MIP_Problem.cc.

References max_space_dimension(), and OK().

00063   : external_space_dim(dim),
00064     internal_space_dim(0),
00065     tableau(),
00066     working_cost(0, Row::Flags()),
00067     mapping(),
00068     base(),
00069     status(PARTIALLY_SATISFIABLE),
00070     pricing(PRICING_STEEPEST_EDGE_FLOAT),
00071     initialized(false),
00072     input_cs(),
00073     first_pending_constraint(0),
00074     input_obj_function(),
00075     opt_mode(MAXIMIZATION),
00076     last_generator(point()),
00077     i_variables() {
00078   // Check for space dimension overflow.
00079   if (dim > max_space_dimension())
00080     throw std::length_error("PPL::MIP_Problem::MIP_Problem(dim, cs, obj, "
00081                             "mode):\n"
00082                             "dim exceeds the maximum allowed "
00083                             "space dimension.");
00084   assert(OK());
00085 }

template<typename In>
Parma_Polyhedra_Library::MIP_Problem::MIP_Problem ( dimension_type  dim,
In  first,
In  last,
const Variables_Set &  int_vars,
const Linear_Expression obj = Linear_Expression::zero(),
Optimization_Mode  mode = MAXIMIZATION 
) [inline]

Builds an MIP problem having space dimension dim from the sequence of constraints in the range $[\mathrm{first}, \mathrm{last})$, the objective function obj and optimization mode mode; those dimensions whose indices occur in int_vars are constrained to take an integer value.

Parameters:
dim The dimension of the vector space enclosing *this.
first An input iterator to the start of the sequence of constraints.
last A past-the-end input iterator to the sequence of constraints.
int_vars The set of variables' indexes that are constrained to take integer values.
obj The objective function (optional argument with default value $0$).
mode The optimization mode (optional argument with default value MAXIMIZATION).
Exceptions:
std::length_error Thrown if dim exceeds max_space_dimension().
std::invalid_argument Thrown if a constraint in the sequence is a strict inequality, if the space dimension of a constraint (resp., of the objective function or of the integer variables) or the space dimension of the integer variable set is strictly greater than dim.

template<typename In>
Parma_Polyhedra_Library::MIP_Problem::MIP_Problem ( dimension_type  dim,
In  first,
In  last,
const Linear_Expression obj = Linear_Expression::zero(),
Optimization_Mode  mode = MAXIMIZATION 
) [inline]

Builds an MIP problem having space dimension dim from the sequence of constraints in the range $[\mathrm{first}, \mathrm{last})$, the objective function obj and optimization mode mode.

Parameters:
dim The dimension of the vector space enclosing *this.
first An input iterator to the start of the sequence of constraints.
last A past-the-end input iterator to the sequence of constraints.
obj The objective function (optional argument with default value $0$).
mode The optimization mode (optional argument with default value MAXIMIZATION).
Exceptions:
std::length_error Thrown if dim exceeds max_space_dimension().
std::invalid_argument Thrown if a constraint in the sequence is a strict inequality or if the space dimension of a constraint (resp., of the objective function or of the integer variables) is strictly greater than dim.

Parma_Polyhedra_Library::MIP_Problem::MIP_Problem ( dimension_type  dim,
const Constraint_System cs,
const Linear_Expression obj = Linear_Expression::zero(),
Optimization_Mode  mode = MAXIMIZATION 
)

Builds an MIP problem having space dimension dim from the constraint system cs, the objective function obj and optimization mode mode.

Parameters:
dim The dimension of the vector space enclosing *this.
cs The constraint system defining the feasible region.
obj The objective function (optional argument with default value $0$).
mode The optimization mode (optional argument with default value MAXIMIZATION).
Exceptions:
std::length_error Thrown if dim exceeds max_space_dimension().
std::invalid_argument Thrown if the constraint system contains any strict inequality or if the space dimension of the constraint system (resp., the objective function) is strictly greater than dim.

Definition at line 87 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Constraint_System::begin(), Parma_Polyhedra_Library::Constraint_System::end(), Parma_Polyhedra_Library::Constraint_System::has_strict_inequalities(), input_cs, max_space_dimension(), OK(), Parma_Polyhedra_Library::Constraint_System::space_dimension(), and Parma_Polyhedra_Library::Linear_Expression::space_dimension().

00091   : external_space_dim(dim),
00092     internal_space_dim(0),
00093     tableau(),
00094     working_cost(0, Row::Flags()),
00095     mapping(),
00096     base(),
00097     status(PARTIALLY_SATISFIABLE),
00098     pricing(PRICING_STEEPEST_EDGE_FLOAT),
00099     initialized(false),
00100     input_cs(),
00101     first_pending_constraint(0),
00102     input_obj_function(obj),
00103     opt_mode(mode),
00104     last_generator(point()),
00105     i_variables() {
00106   // Check for space dimension overflow.
00107   if (dim > max_space_dimension())
00108     throw std::length_error("PPL::MIP_Problem::MIP_Problem(dim, cs, obj, "
00109                             "mode):\n"
00110                             "dim exceeds the maximum allowed"
00111                             "space dimension.");
00112   // Check the objective function.
00113   if (obj.space_dimension() > dim) {
00114     std::ostringstream s;
00115     s << "PPL::MIP_Problem::MIP_Problem(dim, cs, obj,"
00116       << " mode):\n"
00117       << "obj.space_dimension() == "<< obj.space_dimension()
00118       << " exceeds dim == "<< dim << ".";
00119     throw std::invalid_argument(s.str());
00120   }
00121   // Check the constraint system.
00122   if (cs.space_dimension() > dim) {
00123     std::ostringstream s;
00124     s << "PPL::MIP_Problem::MIP_Problem(dim, cs, obj, mode):\n"
00125       << "cs.space_dimension == " << cs.space_dimension() << " exceeds dim == "
00126       << dim << ".";
00127     throw std::invalid_argument(s.str());
00128   }
00129   if (cs.has_strict_inequalities())
00130     throw std::invalid_argument("PPL::MIP_Problem::"
00131                                 "MIP_Problem(d, cs, obj, m):\n"
00132                                 "cs contains strict inequalities.");
00133   // Actually copy the constraints.
00134   input_cs.insert(input_cs.end(), cs.begin(), cs.end());
00135   assert(OK());
00136 }

Parma_Polyhedra_Library::MIP_Problem::MIP_Problem ( const MIP_Problem y  )  [inline]

Ordinary copy-constructor.

Definition at line 43 of file MIP_Problem.inlines.hh.

References OK().

00044   : external_space_dim(y.external_space_dim),
00045     internal_space_dim(y.internal_space_dim),
00046     tableau(y.tableau),
00047     working_cost(y.working_cost),
00048     mapping(y.mapping),
00049     base(y.base),
00050     status(y.status),
00051     pricing(y.pricing),
00052     initialized(y.initialized),
00053     input_cs(y.input_cs),
00054     first_pending_constraint(y.first_pending_constraint),
00055     input_obj_function(y.input_obj_function),
00056     opt_mode(y.opt_mode),
00057     last_generator(y.last_generator),
00058     i_variables(y.i_variables) {
00059   assert(OK());
00060 }

Parma_Polyhedra_Library::MIP_Problem::~MIP_Problem (  )  [inline]

Destructor.

Definition at line 63 of file MIP_Problem.inlines.hh.

00063                           {
00064 }


Member Function Documentation

MIP_Problem & Parma_Polyhedra_Library::MIP_Problem::operator= ( const MIP_Problem y  )  [inline]

Assignment operator.

Definition at line 140 of file MIP_Problem.inlines.hh.

References swap().

00140                                            {
00141   MIP_Problem tmp(y);
00142   swap(tmp);
00143   return *this;
00144 }

dimension_type Parma_Polyhedra_Library::MIP_Problem::max_space_dimension (  )  [inline, static]

Returns the maximum space dimension an MIP_Problem can handle.

Definition at line 32 of file MIP_Problem.inlines.hh.

References Parma_Polyhedra_Library::Constraint::max_space_dimension().

Referenced by add_space_dimensions_and_embed(), and MIP_Problem().

00032                                  {
00033   return Constraint::max_space_dimension();
00034 }

dimension_type Parma_Polyhedra_Library::MIP_Problem::space_dimension (  )  const [inline]

const Variables_Set & Parma_Polyhedra_Library::MIP_Problem::integer_space_dimensions (  )  const [inline]

Returns a set containing all the variables' indexes constrained to be integral.

Definition at line 104 of file MIP_Problem.inlines.hh.

References i_variables.

Referenced by is_satisfiable(), operator<<(), and solve().

00104                                             {
00105   return i_variables;
00106 }

MIP_Problem::const_iterator Parma_Polyhedra_Library::MIP_Problem::constraints_begin (  )  const [inline]

Returns a read-only iterator to the first constraint defining the feasible region.

Definition at line 94 of file MIP_Problem.inlines.hh.

References input_cs.

Referenced by operator<<().

00094                                      {
00095   return input_cs.begin();
00096 }

MIP_Problem::const_iterator Parma_Polyhedra_Library::MIP_Problem::constraints_end (  )  const [inline]

Returns a past-the-end read-only iterator to the sequence of constraints defining the feasible region.

Definition at line 99 of file MIP_Problem.inlines.hh.

References input_cs.

Referenced by operator<<().

00099                                    {
00100   return input_cs.end();
00101 }

const Linear_Expression & Parma_Polyhedra_Library::MIP_Problem::objective_function (  )  const [inline]

Returns the objective function.

Definition at line 78 of file MIP_Problem.inlines.hh.

References input_obj_function.

Referenced by operator<<().

00078                                       {
00079   return input_obj_function;
00080 }

Optimization_Mode Parma_Polyhedra_Library::MIP_Problem::optimization_mode (  )  const [inline]

Returns the optimization mode.

Definition at line 83 of file MIP_Problem.inlines.hh.

References opt_mode.

Referenced by operator<<(), and solve_mip().

00083                                      {
00084   return opt_mode;
00085 }

void Parma_Polyhedra_Library::MIP_Problem::clear (  )  [inline]

Resets *this to be equal to the trivial MIP problem.

The space dimension is reset to $0$.

Definition at line 147 of file MIP_Problem.inlines.hh.

References swap().

00147                    {
00148   MIP_Problem tmp;
00149   swap(tmp);
00150 }

void Parma_Polyhedra_Library::MIP_Problem::add_space_dimensions_and_embed ( dimension_type  m  ) 

Adds m new space dimensions and embeds the old MIP problem in the new vector space.

Parameters:
m The number of dimensions to add.
Exceptions:
std::length_error Thrown if adding m new space dimensions would cause the vector space to exceed dimension max_space_dimension().
The new space dimensions will be those having the highest indexes in the new MIP problem; they are initially unconstrained.

Definition at line 336 of file MIP_Problem.cc.

References external_space_dim, max_space_dimension(), OK(), PARTIALLY_SATISFIABLE, space_dimension(), status, and UNSATISFIABLE.

Referenced by Parma_Polyhedra_Library::Polyhedron::simplify_using_context_assign(), and Parma_Polyhedra_Library::Polyhedron::strongly_minimize_constraints().

00336                                                                      {
00337   // The space dimension of the resulting MIP problem should not
00338   // overflow the maximum allowed space dimension.
00339   if (m > max_space_dimension() - space_dimension())
00340     throw std::length_error("PPL::MIP_Problem::"
00341                             "add_space_dimensions_and_embed(m):\n"
00342                             "adding m new space dimensions exceeds "
00343                             "the maximum allowed space dimension.");
00344   external_space_dim += m;
00345   if (status != UNSATISFIABLE)
00346     status = PARTIALLY_SATISFIABLE;
00347   assert(OK());
00348 }

void Parma_Polyhedra_Library::MIP_Problem::add_to_integer_space_dimensions ( const Variables_Set &  i_vars  ) 

Sets the variables whose indexes are in set i_vars to be integer space dimensions.

Exceptions:
std::invalid_argument Thrown if some index in i_vars does not correspond to a space dimension in *this.

Definition at line 352 of file MIP_Problem.cc.

References external_space_dim, i_variables, PARTIALLY_SATISFIABLE, status, and UNSATISFIABLE.

00352                                                              {
00353   if (i_vars.space_dimension() > external_space_dim)
00354     throw std::invalid_argument("PPL::MIP_Problem::"
00355                                 "add_to_integer_space_dimension(i_vars):\n"
00356                                 "*this and i_vars are dimension"
00357                                 "incompatible.");
00358   const dimension_type original_size = i_variables.size();
00359   i_variables.insert(i_vars.begin(), i_vars.end());
00360   // If a new integral variable was inserted, set the internal status to
00361   // PARTIALLY_SATISFIABLE.
00362   if (i_variables.size() != original_size && status != UNSATISFIABLE)
00363     status = PARTIALLY_SATISFIABLE;
00364 }

void Parma_Polyhedra_Library::MIP_Problem::add_constraint ( const Constraint c  ) 

Adds a copy of constraint c to the MIP problem.

Exceptions:
std::invalid_argument Thrown if the constraint c is a strict inequality or if its space dimension is strictly greater than the space dimension of *this.

Definition at line 139 of file MIP_Problem.cc.

References input_cs, Parma_Polyhedra_Library::Constraint::is_strict_inequality(), OK(), PARTIALLY_SATISFIABLE, Parma_Polyhedra_Library::Constraint::space_dimension(), space_dimension(), status, and UNSATISFIABLE.

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), is_mip_satisfiable(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), Parma_Polyhedra_Library::Polyhedron::simplify_using_context_assign(), and solve_mip().

00139                                                   {
00140   if (space_dimension() < c.space_dimension()) {
00141     std::ostringstream s;
00142     s << "PPL::MIP_Problem::add_constraint(c):\n"
00143       << "c.space_dimension() == "<< c.space_dimension() << " exceeds "
00144       "this->space_dimension == " << space_dimension() << ".";
00145     throw std::invalid_argument(s.str());
00146   }
00147   if (c.is_strict_inequality())
00148     throw std::invalid_argument("PPL::MIP_Problem::add_constraint(c):\n"
00149                                 "c is a strict inequality.");
00150   input_cs.push_back(c);
00151   if (status != UNSATISFIABLE)
00152     status = PARTIALLY_SATISFIABLE;
00153   assert(OK());
00154 }

void Parma_Polyhedra_Library::MIP_Problem::add_constraints ( const Constraint_System cs  ) 

Adds a copy of the constraints in cs to the MIP problem.

Exceptions:
std::invalid_argument Thrown if the constraint system cs contains any strict inequality or if its space dimension is strictly greater than the space dimension of *this.

Definition at line 157 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Constraint_System::begin(), Parma_Polyhedra_Library::Constraint_System::end(), Parma_Polyhedra_Library::Constraint_System::has_strict_inequalities(), input_cs, OK(), PARTIALLY_SATISFIABLE, Parma_Polyhedra_Library::Constraint_System::space_dimension(), space_dimension(), status, and UNSATISFIABLE.

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), Parma_Polyhedra_Library::Polyhedron::simplify_using_context_assign(), and Parma_Polyhedra_Library::Polyhedron::strongly_minimize_constraints().

00157                                                            {
00158   if (space_dimension() < cs.space_dimension()) {
00159     std::ostringstream s;
00160     s << "PPL::MIP_Problem::add_constraints(cs):\n"
00161       << "cs.space_dimension() == " << cs.space_dimension()
00162       << " exceeds this->space_dimension() == " << this->space_dimension()
00163       << ".";
00164     throw std::invalid_argument(s.str());
00165   }
00166   if (cs.has_strict_inequalities())
00167     throw std::invalid_argument("PPL::MIP_Problem::add_constraints(cs):\n"
00168                                 "cs contains strict inequalities.");
00169   input_cs.insert(input_cs.end(), cs.begin(), cs.end());
00170   if (status != UNSATISFIABLE)
00171     status = PARTIALLY_SATISFIABLE;
00172   assert(OK());
00173 }

void Parma_Polyhedra_Library::MIP_Problem::set_objective_function ( const Linear_Expression obj  ) 

Sets the objective function to obj.

Exceptions:
std::invalid_argument Thrown if the space dimension of obj is strictly greater than the space dimension of *this.

Definition at line 176 of file MIP_Problem.cc.

References input_obj_function, OK(), OPTIMIZED, SATISFIABLE, Parma_Polyhedra_Library::Linear_Expression::space_dimension(), space_dimension(), status, and UNBOUNDED.

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), and Parma_Polyhedra_Library::Polyhedron::strongly_minimize_constraints().

00176                                                                    {
00177   if (space_dimension() < obj.space_dimension()) {
00178     std::ostringstream s;
00179     s << "PPL::MIP_Problem::set_objective_function(obj):\n"
00180       << "obj.space_dimension() == " << obj.space_dimension()
00181       << " exceeds this->space_dimension == " << space_dimension()
00182       << ".";
00183     throw std::invalid_argument(s.str());
00184   }
00185   input_obj_function = obj;
00186   if (status == UNBOUNDED || status == OPTIMIZED)
00187     status = SATISFIABLE;
00188   assert(OK());
00189 }

void Parma_Polyhedra_Library::MIP_Problem::set_optimization_mode ( Optimization_Mode  mode  )  [inline]

bool Parma_Polyhedra_Library::MIP_Problem::is_satisfiable (  )  const

Checks satisfiability of *this.

Returns:
true if and only if the MIP problem is satisfiable.

Definition at line 210 of file MIP_Problem.cc.

References i_variables, integer_space_dimensions(), is_lp_satisfiable(), is_mip_satisfiable(), last_generator, OK(), OPTIMIZED, PARTIALLY_SATISFIABLE, SATISFIABLE, status, UNBOUNDED, and UNSATISFIABLE.

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), feasible_point(), and Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape().

00210                                      {
00211   // Check `status' to filter out trivial cases.
00212   switch (status) {
00213   case UNSATISFIABLE:
00214     assert(OK());
00215     return false;
00216   case SATISFIABLE:
00217     // Intentionally fall through
00218   case UNBOUNDED:
00219     // Intentionally fall through.
00220   case OPTIMIZED:
00221     assert(OK());
00222     return true;
00223   case PARTIALLY_SATISFIABLE:
00224     { // LP case.
00225       if (i_variables.empty())
00226         return is_lp_satisfiable();
00227       // MIP Case.
00228       const Variables_Set this_variables_set = integer_space_dimensions();
00229       MIP_Problem& x = const_cast<MIP_Problem&>(*this);
00230       Generator p = point();
00231       // This disable the Variable integrality check in OK() until we will
00232       // find a feasible point.
00233       x.i_variables.clear();
00234       x.is_lp_satisfiable();
00235          if (is_mip_satisfiable(x, p, this_variables_set)) {
00236         x.last_generator = p;
00237         x.status = SATISFIABLE;
00238         // Restore i_variables;
00239         x.i_variables = this_variables_set;
00240         return true;
00241          }
00242       else {
00243         x.status = UNSATISFIABLE;
00244         // Restore i_variables;
00245         x.i_variables = this_variables_set;
00246         return false;
00247       }
00248     }
00249   }
00250   // We should not be here!
00251   throw std::runtime_error("PPL internal error");
00252 }

PPL::MIP_Problem_Status Parma_Polyhedra_Library::MIP_Problem::solve (  )  const

Optimizes the MIP problem.

Returns:
An MIP_Problem_Status flag indicating the outcome of the optimization attempt (unfeasible, unbounded or optimized problem).

Definition at line 255 of file MIP_Problem.cc.

References i_variables, integer_space_dimensions(), is_lp_satisfiable(), last_generator, OK(), OPTIMIZED, Parma_Polyhedra_Library::OPTIMIZED_MIP_PROBLEM, PARTIALLY_SATISFIABLE, SATISFIABLE, second_phase(), solve_mip(), status, UNBOUNDED, Parma_Polyhedra_Library::UNBOUNDED_MIP_PROBLEM, Parma_Polyhedra_Library::UNFEASIBLE_MIP_PROBLEM, and UNSATISFIABLE.

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Octagonal_Shape< T >::bounds(), Parma_Polyhedra_Library::BD_Shape< T >::bounds(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Octagonal_Shape< T >::max_min(), Parma_Polyhedra_Library::BD_Shape< T >::max_min(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), optimizing_point(), Parma_Polyhedra_Library::Polyhedron::simplify_using_context_assign(), and Parma_Polyhedra_Library::Polyhedron::strongly_minimize_constraints().

00255                            {
00256   switch (status) {
00257   case UNSATISFIABLE:
00258     assert(OK());
00259     return UNFEASIBLE_MIP_PROBLEM;
00260   case UNBOUNDED:
00261     assert(OK());
00262     return UNBOUNDED_MIP_PROBLEM;
00263   case OPTIMIZED:
00264     assert(OK());
00265     return OPTIMIZED_MIP_PROBLEM;
00266   case SATISFIABLE:
00267     // Intentionally fall through
00268   case PARTIALLY_SATISFIABLE:
00269     {
00270       MIP_Problem& x = const_cast<MIP_Problem&>(*this);
00271       if (i_variables.empty()) {
00272         // LP Problem case.
00273         if (is_lp_satisfiable()) {
00274           x.second_phase();
00275           if (x.status == UNBOUNDED)
00276             return UNBOUNDED_MIP_PROBLEM;
00277           else {
00278             assert(x.status == OPTIMIZED);
00279             return OPTIMIZED_MIP_PROBLEM;
00280           }
00281         }
00282         return UNFEASIBLE_MIP_PROBLEM;
00283       }
00284       // MIP Problem case.
00285       // This disable the Variable integrality check in OK() until we will find
00286       // an optimizing point.
00287       const Variables_Set this_variables_set = integer_space_dimensions();
00288       x.i_variables.clear();
00289       if (x.is_lp_satisfiable())
00290         x.second_phase();
00291       else {
00292         x.status = UNSATISFIABLE;
00293         // Restore i_variables;
00294         x.i_variables = this_variables_set;
00295         return UNFEASIBLE_MIP_PROBLEM;
00296       }
00297       DIRTY_TEMP0(mpq_class, incumbent_solution);
00298       Generator g = point();
00299       bool have_incumbent_solution = false;
00300 
00301       MIP_Problem mip_copy(*this);
00302       // Treat this MIP_Problem as an LP one: we have to deal with
00303       // the relaxation in solve_mip().
00304       mip_copy.i_variables.clear();
00305       MIP_Problem_Status mip_status = solve_mip(have_incumbent_solution,
00306                                                 incumbent_solution, g,
00307                                                 mip_copy,
00308                                                 this_variables_set);
00309       // Restore i_variables;
00310       x.i_variables = this_variables_set;
00311       switch (mip_status) {
00312       case UNFEASIBLE_MIP_PROBLEM:
00313         x.status = UNSATISFIABLE;
00314         break;
00315       case UNBOUNDED_MIP_PROBLEM:
00316         x.status = UNBOUNDED;
00317         // A feasible point has been set in `solve_mip()', so that
00318         // a call to `feasible_point' will be successful.
00319         x.last_generator = g;
00320         break;
00321       case OPTIMIZED_MIP_PROBLEM:
00322         x.status = OPTIMIZED;
00323         // Set the internal generator.
00324         x.last_generator = g;
00325         break;
00326       }
00327       assert(OK());
00328       return mip_status;
00329     }
00330   }
00331   // We should not be here!
00332   throw std::runtime_error("PPL internal error");
00333 }

void Parma_Polyhedra_Library::MIP_Problem::evaluate_objective_function ( const Generator evaluating_point,
Coefficient num,
Coefficient den 
) const

Sets num and den so that $\frac{num}{den}$ is the result of evaluating the objective function on evaluating_point.

Parameters:
evaluating_point The point on which the objective function will be evaluated.
num On exit will contain the numerator of the evaluated value.
den On exit will contain the denominator of the evaluated value.
Exceptions:
std::invalid_argument Thrown if *this and evaluating_point are dimension-incompatible or if the generator evaluating_point is not a point.

Definition at line 1453 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::Generator::coefficient(), Parma_Polyhedra_Library::Generator::divisor(), Parma_Polyhedra_Library::Linear_Expression::inhomogeneous_term(), input_obj_function, Parma_Polyhedra_Library::Generator::is_point(), Parma_Polyhedra_Library::normalize2(), Parma_Polyhedra_Library::Linear_Expression::space_dimension(), space_dimension(), and Parma_Polyhedra_Library::Generator::space_dimension().

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Octagonal_Shape< T >::max_min(), Parma_Polyhedra_Library::BD_Shape< T >::max_min(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), optimal_value(), and solve_mip().

01455                                                         {
01456   const dimension_type ep_space_dim = evaluating_point.space_dimension();
01457   if (space_dimension() < ep_space_dim)
01458     throw std::invalid_argument("PPL::MIP_Problem::"
01459                                 "evaluate_objective_function(p, n, d):\n"
01460                                 "*this and p are dimension incompatible.");
01461   if (!evaluating_point.is_point())
01462     throw std::invalid_argument("PPL::MIP_Problem::"
01463                                 "evaluate_objective_function(p, n, d):\n"
01464                                 "p is not a point.");
01465 
01466   // Compute the smallest space dimension  between `input_obj_function'
01467   // and `evaluating_point'.
01468   const dimension_type working_space_dim
01469     = std::min(ep_space_dim, input_obj_function.space_dimension());
01470   // Compute the optimal value of the cost function.
01471   const Coefficient& divisor = evaluating_point.divisor();
01472   ext_n = input_obj_function.inhomogeneous_term() * divisor;
01473   for (dimension_type i = working_space_dim; i-- > 0; )
01474     ext_n += evaluating_point.coefficient(Variable(i))
01475       * input_obj_function.coefficient(Variable(i));
01476   // Numerator and denominator should be coprime.
01477   normalize2(ext_n, divisor, ext_n, ext_d);
01478 }

const PPL::Generator & Parma_Polyhedra_Library::MIP_Problem::feasible_point (  )  const

Returns a feasible point for *this, if it exists.

Exceptions:
std::domain_error Thrown if the MIP problem is not satisfiable.

Definition at line 192 of file MIP_Problem.cc.

References is_satisfiable(), and last_generator.

00192                                      {
00193   if (is_satisfiable())
00194     return last_generator;
00195   else
00196     throw std::domain_error("PPL::MIP_Problem::feasible_point():\n"
00197                             "*this is not satisfiable.");
00198 }

const PPL::Generator & Parma_Polyhedra_Library::MIP_Problem::optimizing_point (  )  const

Returns an optimal point for *this, if it exists.

Exceptions:
std::domain_error Thrown if *this doesn't not have an optimizing point, i.e., if the MIP problem is unbounded or not satisfiable.

Definition at line 201 of file MIP_Problem.cc.

References last_generator, Parma_Polyhedra_Library::OPTIMIZED_MIP_PROBLEM, and solve().

Referenced by Parma_Polyhedra_Library::BD_Shape< T >::BD_Shape(), Parma_Polyhedra_Library::Box< ITV >::Box(), Parma_Polyhedra_Library::Octagonal_Shape< T >::max_min(), Parma_Polyhedra_Library::BD_Shape< T >::max_min(), Parma_Polyhedra_Library::Octagonal_Shape< T >::Octagonal_Shape(), and optimal_value().

00201                                        {
00202   if (solve() == OPTIMIZED_MIP_PROBLEM)
00203     return last_generator;
00204   else
00205     throw std::domain_error("PPL::MIP_Problem::optimizing_point():\n"
00206                             "*this doesn't have an optimizing point.");
00207 }

void Parma_Polyhedra_Library::MIP_Problem::optimal_value ( Coefficient num,
Coefficient den 
) const [inline]

Sets num and den so that $\frac{num}{den}$ is the solution of the optimization problem.

Exceptions:
std::domain_error Thrown if *this doesn't not have an optimizing point, i.e., if the MIP problem is unbounded or not satisfiable.

Definition at line 88 of file MIP_Problem.inlines.hh.

References evaluate_objective_function(), and optimizing_point().

Referenced by Parma_Polyhedra_Library::Octagonal_Shape< T >::max_min(), and Parma_Polyhedra_Library::BD_Shape< T >::max_min().

00088                                                                    {
00089   const Generator& g = optimizing_point();
00090   evaluate_objective_function(g, num, den);
00091 }

bool Parma_Polyhedra_Library::MIP_Problem::OK (  )  const

Checks if all the invariants are satisfied.

Definition at line 1735 of file MIP_Problem.cc.

References ascii_dump(), base, Parma_Polyhedra_Library::Generator::coefficient(), Parma_Polyhedra_Library::Generator::divisor(), external_space_dim, first_pending_constraint, Parma_Polyhedra_Library::gcd_assign(), i_variables, initialized, input_cs, input_obj_function, is_satisfied(), last_generator, mapping, Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), Parma_Polyhedra_Library::Generator::OK(), Parma_Polyhedra_Library::Linear_Expression::OK(), Parma_Polyhedra_Library::Matrix::OK(), Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::Generator::space_dimension(), Parma_Polyhedra_Library::Linear_Expression::space_dimension(), status, tableau, TEMP_INTEGER, UNSATISFIABLE, and working_cost.

Referenced by add_constraint(), add_constraints(), add_space_dimensions_and_embed(), ascii_load(), is_lp_satisfiable(), is_satisfiable(), MIP_Problem(), process_pending_constraints(), second_phase(), set_objective_function(), set_optimization_mode(), and solve().

01735                          {
01736 #ifndef NDEBUG
01737   using std::endl;
01738   using std::cerr;
01739 #endif
01740   const dimension_type input_cs_num_rows = input_cs.size();
01741   // Check that every member used is OK.
01742 
01743   for (dimension_type i = input_cs_num_rows; i-- > 0; )
01744     if (!input_cs[i].OK())
01745       return false;
01746 
01747   if (!tableau.OK() || !input_obj_function.OK() || !last_generator.OK())
01748     return false;
01749 
01750   // Constraint system should contain no strict inequalities.
01751   for (dimension_type i = input_cs_num_rows; i-- > 0; )
01752     if (input_cs[i].is_strict_inequality()) {
01753 #ifndef NDEBUG
01754       cerr << "The feasible region of the MIP_Problem is defined by "
01755            << "a constraint system containing strict inequalities."
01756            << endl;
01757       ascii_dump(cerr);
01758 #endif
01759       return false;
01760     }
01761 
01762   // Constraint system and objective function should be dimension compatible.
01763   if (external_space_dim < input_obj_function.space_dimension()) {
01764 #ifndef NDEBUG
01765     cerr << "The MIP_Problem and the objective function have "
01766          << "incompatible space dimensions ("
01767          << external_space_dim << " < "
01768          << input_obj_function.space_dimension() << ")."
01769          << endl;
01770     ascii_dump(cerr);
01771 #endif
01772     return false;
01773   }
01774 
01775   if (status != UNSATISFIABLE && initialized) {
01776     // Here `last_generator' has to be meaningful.
01777     // Check for dimension compatibility and actual feasibility.
01778     if (external_space_dim != last_generator.space_dimension()) {
01779 #ifndef NDEBUG
01780       cerr << "The MIP_Problem and the cached feasible point have "
01781            << "incompatible space dimensions ("
01782            << external_space_dim << " != "
01783            << last_generator.space_dimension() << ")."
01784            << endl;
01785       ascii_dump(cerr);
01786 #endif
01787       return false;
01788     }
01789 
01790     for (dimension_type i = 0; i < first_pending_constraint; ++i)
01791       if (!is_satisfied(input_cs[i], last_generator)) {
01792 #ifndef NDEBUG
01793         cerr << "The cached feasible point does not belong to "
01794              << "the feasible region of the MIP_Problem."
01795              << endl;
01796         ascii_dump(cerr);
01797 #endif
01798         return false;
01799       }
01800 
01801     // Check that every integer declared variable is really integer.
01802     // in the solution found.
01803     if (!i_variables.empty()) {
01804       TEMP_INTEGER(gcd);
01805       for (Variables_Set::const_iterator v_it = i_variables.begin(),
01806              v_end = i_variables.end(); v_it != v_end; ++v_it) {
01807         gcd_assign(gcd, last_generator.coefficient(Variable(*v_it)),
01808                    last_generator.divisor());
01809         if (gcd != last_generator.divisor())
01810           return false;
01811       }
01812     }
01813 
01814     const dimension_type tableau_nrows = tableau.num_rows();
01815     const dimension_type tableau_ncols = tableau.num_columns();
01816 
01817     // The number of rows in the tableau and base should be equal.
01818     if (tableau_nrows != base.size()) {
01819 #ifndef NDEBUG
01820       cerr << "tableau and base have incompatible sizes" << endl;
01821       ascii_dump(cerr);
01822 #endif
01823       return false;
01824     }
01825     // The size of `mapping' should be equal to the space dimension
01826     // of `input_cs' plus one.
01827     if (mapping.size() != external_space_dim + 1) {
01828 #ifndef NDEBUG
01829       cerr << "`input_cs' and `mapping' have incompatible sizes" << endl;
01830       ascii_dump(cerr);
01831 #endif
01832       return false;
01833     }
01834 
01835     // The number of columns in the tableau and working_cost should be equal.
01836     if (tableau_ncols != working_cost.size()) {
01837 #ifndef NDEBUG
01838       cerr << "tableau and working_cost have incompatible sizes" << endl;
01839       ascii_dump(cerr);
01840 #endif
01841       return false;
01842     }
01843 
01844     // The vector base should contain indices of tableau's columns.
01845     for (dimension_type i = base.size(); i-- > 0; ) {
01846       if (base[i] > tableau_ncols) {
01847 #ifndef NDEBUG
01848         cerr << "base contains an invalid column index" << endl;
01849         ascii_dump(cerr);
01850 #endif
01851         return false;
01852       }
01853       // tableau[i][base[i] must be different from zero.
01854       // tableau[i][base[j], with i different from j, must not be a zero.
01855       for (dimension_type j = tableau_nrows; j-- > 0; )
01856         if (i != j && tableau[j][base[i]] != 0) {
01857 #ifndef NDEBUG
01858           cerr << "tableau[i][base[i] must be different from zero" << endl;
01859           ascii_dump(cerr);
01860 #endif
01861           return false;
01862         }
01863       if (tableau[i][base[i]] == 0) {
01864 #ifndef NDEBUG
01865         cerr << "tableau[i][base[j], with i different from j, must not be "
01866              << "a zero" << endl;
01867         ascii_dump(cerr);
01868 #endif
01869         return false;
01870       }
01871     }
01872 
01873     // The last column of the tableau must contain only zeroes.
01874     for (dimension_type i = tableau_nrows; i-- > 0; )
01875       if (tableau[i][tableau_ncols-1] != 0) {
01876 #ifndef NDEBUG
01877         cerr << "the last column of the tableau must contain only"
01878           "zeroes"<< endl;
01879         ascii_dump(cerr);
01880 #endif
01881         return false;
01882       }
01883    }
01884 
01885    // All checks passed.
01886    return true;
01887 }

void Parma_Polyhedra_Library::MIP_Problem::ascii_dump (  )  const

Writes to std::cerr an ASCII representation of *this.

Referenced by ascii_dump(), and OK().

void Parma_Polyhedra_Library::MIP_Problem::ascii_dump ( std::ostream &  s  )  const

Writes to s an ASCII representation of *this.

Definition at line 1890 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Generator::ascii_dump(), Parma_Polyhedra_Library::Row::ascii_dump(), Parma_Polyhedra_Library::Matrix::ascii_dump(), Parma_Polyhedra_Library::Linear_Expression::ascii_dump(), ascii_dump(), base, external_space_dim, first_pending_constraint, i_variables, initialized, input_cs, input_obj_function, internal_space_dim, last_generator, mapping, Parma_Polyhedra_Library::MAXIMIZATION, opt_mode, OPTIMIZED, PARTIALLY_SATISFIABLE, pricing, PRICING_STEEPEST_EDGE_EXACT, PRICING_STEEPEST_EDGE_FLOAT, PRICING_TEXTBOOK, SATISFIABLE, Parma_Polyhedra_Library::Row::size(), status, tableau, UNBOUNDED, UNSATISFIABLE, and working_cost.

01890                                               {
01891   using namespace IO_Operators;
01892   s << "\nexternal_space_dim: " << external_space_dim << " \n";
01893   s << "\ninternal_space_dim: " << internal_space_dim << " \n";
01894 
01895   const dimension_type input_cs_size = input_cs.size();
01896 
01897   s << "\ninput_cs( " << input_cs_size << " )\n";
01898   for (dimension_type i = 0; i < input_cs_size; ++i)
01899     input_cs[i].ascii_dump(s);
01900 
01901   s << "\nfirst_pending_constraint: " <<  first_pending_constraint
01902     << std::endl;
01903 
01904   s << "\ninput_obj_function\n";
01905   input_obj_function.ascii_dump(s);
01906   s << "\nopt_mode "
01907     << (opt_mode == MAXIMIZATION ? "MAXIMIZATION" : "MINIMIZATION") << "\n";
01908 
01909   s << "\ninitialized: " << (initialized ? "YES" : "NO") << "\n";
01910   s << "\npricing: ";
01911   switch (pricing) {
01912   case PRICING_STEEPEST_EDGE_FLOAT:
01913     s << "PRICING_STEEPEST_EDGE_FLOAT";
01914     break;
01915   case PRICING_STEEPEST_EDGE_EXACT:
01916     s << "PRICING_STEEPEST_EDGE_EXACT";
01917     break;
01918   case PRICING_TEXTBOOK:
01919     s << "PRICING_TEXTBOOK";
01920     break;
01921   }
01922   s << "\n";
01923 
01924   s << "\nstatus: ";
01925   switch (status) {
01926   case UNSATISFIABLE:
01927     s << "UNSATISFIABLE";
01928     break;
01929   case SATISFIABLE:
01930     s << "SATISFIABLE";
01931     break;
01932   case UNBOUNDED:
01933     s << "UNBOUNDED";
01934     break;
01935   case OPTIMIZED:
01936     s << "OPTIMIZED";
01937     break;
01938   case PARTIALLY_SATISFIABLE:
01939     s << "PARTIALLY_SATISFIABLE";
01940     break;
01941   }
01942   s << "\n";
01943 
01944   s << "\ntableau\n";
01945   tableau.ascii_dump(s);
01946   s << "\nworking_cost( " << working_cost.size()<< " )\n";
01947   working_cost.ascii_dump(s);
01948 
01949   const dimension_type base_size = base.size();
01950   s << "\nbase( " << base_size << " )\n";
01951   for (dimension_type i = 0; i != base_size; ++i)
01952     s << base[i] << ' ';
01953 
01954   s << "\nlast_generator\n";
01955   last_generator.ascii_dump(s);
01956 
01957   const dimension_type mapping_size = mapping.size();
01958   s << "\nmapping( " << mapping_size << " )\n";
01959   for (dimension_type i = 1; i < mapping_size; ++i)
01960     s << "\n"<< i << " -> " << mapping[i].first << " -> " << mapping[i].second
01961       << ' ';
01962 
01963   s << "\n\ninteger_variables";
01964   i_variables.ascii_dump(s);
01965 }

void Parma_Polyhedra_Library::MIP_Problem::print (  )  const

Prints *this to std::cerr using operator<<.

bool Parma_Polyhedra_Library::MIP_Problem::ascii_load ( std::istream &  s  ) 

Loads from s an ASCII representation (as produced by ascii_dump(std::ostream&) const) and sets *this accordingly. Returns true if successful, false otherwise.

Definition at line 1970 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Generator::ascii_load(), Parma_Polyhedra_Library::Row::ascii_load(), Parma_Polyhedra_Library::Matrix::ascii_load(), Parma_Polyhedra_Library::Linear_Expression::ascii_load(), Parma_Polyhedra_Library::Constraint::ascii_load(), base, external_space_dim, first_pending_constraint, i_variables, initialized, input_cs, input_obj_function, internal_space_dim, last_generator, mapping, Parma_Polyhedra_Library::MAXIMIZATION, Parma_Polyhedra_Library::MINIMIZATION, Parma_Polyhedra_Library::Matrix::num_columns(), OK(), OPTIMIZED, PARTIALLY_SATISFIABLE, pricing, PRICING_STEEPEST_EDGE_EXACT, PRICING_STEEPEST_EDGE_FLOAT, PRICING_TEXTBOOK, SATISFIABLE, set_optimization_mode(), status, tableau, UNBOUNDED, UNSATISFIABLE, working_cost, and Parma_Polyhedra_Library::Constraint::zero_dim_positivity().

01970                                         {
01971   std::string str;
01972 if (!(s >> str) || str != "external_space_dim:")
01973     return false;
01974 
01975 if (!(s >> external_space_dim))
01976     return false;
01977 
01978 if (!(s >> str) || str != "internal_space_dim:")
01979     return false;
01980 
01981 if (!(s >> internal_space_dim))
01982     return false;
01983 
01984  if (!(s >> str) || str != "input_cs(")
01985     return false;
01986 
01987   dimension_type input_cs_size;
01988 
01989   if (!(s >> input_cs_size))
01990     return false;
01991 
01992   if (!(s >> str) || str != ")")
01993     return false;
01994 
01995   Constraint c(Constraint::zero_dim_positivity());
01996   for (dimension_type i = 0; i < input_cs_size; ++i) {
01997     if (!c.ascii_load(s))
01998       return false;
01999     input_cs.push_back(c);
02000   }
02001 
02002   if (!(s >> str) || str != "first_pending_constraint:")
02003     return false;
02004 
02005   if (!(s >> first_pending_constraint))
02006     return false;
02007 
02008   if (!(s >> str) || str != "input_obj_function")
02009     return false;
02010 
02011   if (!input_obj_function.ascii_load(s))
02012     return false;
02013 
02014   if (!(s >> str) || str != "opt_mode")
02015     return false;
02016 
02017   if (!(s >> str))
02018     return false;
02019 
02020   if (str == "MAXIMIZATION")
02021     set_optimization_mode(MAXIMIZATION);
02022   else {
02023     if (str != "MINIMIZATION")
02024       return false;
02025     set_optimization_mode(MINIMIZATION);
02026   }
02027 
02028   if (!(s >> str) || str != "initialized:")
02029     return false;
02030   if (!(s >> str))
02031     return false;
02032   if (str == "YES")
02033     initialized = true;
02034   else if (str == "NO")
02035     initialized = false;
02036   else
02037     return false;
02038 
02039   if (!(s >> str) || str != "pricing:")
02040     return false;
02041   if (!(s >> str))
02042     return false;
02043   if (str == "PRICING_STEEPEST_EDGE_FLOAT")
02044     pricing = PRICING_STEEPEST_EDGE_FLOAT;
02045   else if (str == "PRICING_STEEPEST_EDGE_EXACT")
02046     pricing = PRICING_STEEPEST_EDGE_EXACT;
02047   else if (str == "PRICING_TEXTBOOK")
02048     pricing = PRICING_TEXTBOOK;
02049   else
02050     return false;
02051 
02052   if (!(s >> str) || str != "status:")
02053     return false;
02054 
02055   if (!(s >> str))
02056     return false;
02057 
02058   if (str == "UNSATISFIABLE")
02059     status = UNSATISFIABLE;
02060   else if (str == "SATISFIABLE")
02061     status = SATISFIABLE;
02062   else if (str == "UNBOUNDED")
02063     status = UNBOUNDED;
02064   else if (str == "OPTIMIZED")
02065     status = OPTIMIZED;
02066   else if (str == "PARTIALLY_SATISFIABLE")
02067     status = PARTIALLY_SATISFIABLE;
02068   else
02069     return false;
02070 
02071   if (!(s >> str) || str != "tableau")
02072     return false;
02073 
02074   if (!tableau.ascii_load(s))
02075     return false;
02076 
02077   if (!(s >> str) || str != "working_cost(")
02078     return false;
02079 
02080   dimension_type working_cost_dim;
02081 
02082   if (!(s >> working_cost_dim))
02083     return false;
02084 
02085   if (!(s >> str) || str != ")")
02086     return false;
02087 
02088   if (!working_cost.ascii_load(s))
02089     return false;
02090 
02091   if (!(s >> str) || str != "base(")
02092     return false;
02093 
02094   dimension_type base_size;
02095   if (!(s >> base_size))
02096     return false;
02097 
02098   if (!(s >> str) || str != ")")
02099     return false;
02100 
02101   dimension_type base_value;
02102   for (dimension_type i = 0; i != base_size; ++i) {
02103     if (!(s >> base_value))
02104       return false;
02105     base.push_back(base_value);
02106   }
02107 
02108   if (!(s >> str) || str != "last_generator")
02109     return false;
02110 
02111   if (!last_generator.ascii_load(s))
02112     return false;
02113 
02114   if (!(s >> str) || str != "mapping(")
02115     return false;
02116 
02117   dimension_type mapping_size;
02118   if (!(s >> mapping_size))
02119     return false;
02120 
02121   if (!(s >> str) || str != ")")
02122     return false;
02123 
02124   dimension_type first_value;
02125   dimension_type second_value;
02126   dimension_type index;
02127 
02128   // The first `mapping' index is never used, so we initialize
02129   // it pushing back a dummy value.
02130   if (tableau.num_columns() != 0)
02131     mapping.push_back(std::make_pair(0, 0));
02132 
02133   for (dimension_type i = 1; i < mapping_size; ++i) {
02134     if (!(s >> index))
02135       return false;
02136     if (!(s >> str) || str != "->")
02137       return false;
02138     if (!(s >> first_value))
02139       return false;
02140     if (!(s >> str) || str != "->")
02141       return false;
02142     if (!(s >> second_value))
02143       return false;
02144     mapping.push_back(std::make_pair(first_value, second_value));
02145   }
02146 
02147   if (!(s >> str) || str != "integer_variables")
02148     return false;
02149 
02150   if (!i_variables.ascii_load(s))
02151     return false;
02152 
02153   assert(OK());
02154   return true;
02155 }

memory_size_type Parma_Polyhedra_Library::MIP_Problem::total_memory_in_bytes (  )  const [inline]

Returns the total size in bytes of the memory occupied by *this.

Definition at line 172 of file MIP_Problem.inlines.hh.

References external_memory_in_bytes().

00172                                          {
00173   return sizeof(*this) + external_memory_in_bytes();
00174 }

memory_size_type Parma_Polyhedra_Library::MIP_Problem::external_memory_in_bytes (  )  const [inline]

Returns the size in bytes of the memory managed by *this.

Definition at line 153 of file MIP_Problem.inlines.hh.

References base, Parma_Polyhedra_Library::Generator::external_memory_in_bytes(), Parma_Polyhedra_Library::Linear_Expression::external_memory_in_bytes(), Parma_Polyhedra_Library::Row::external_memory_in_bytes(), Parma_Polyhedra_Library::Matrix::external_memory_in_bytes(), input_cs, input_obj_function, last_generator, mapping, tableau, and working_cost.

Referenced by total_memory_in_bytes().

00153                                             {
00154   memory_size_type n
00155     = tableau.external_memory_in_bytes()
00156     + working_cost.external_memory_in_bytes()
00157     + input_obj_function.external_memory_in_bytes()
00158     + last_generator.external_memory_in_bytes();
00159   // Adding the external memory for `input_cs'.
00160   n += input_cs.capacity() * sizeof(Constraint);
00161   for (const_iterator i = input_cs.begin(),
00162          i_end = input_cs.end(); i != i_end; ++i)
00163     n += (i->external_memory_in_bytes());
00164   // Adding the external memory for `base'.
00165   n += base.capacity() * sizeof(dimension_type);
00166   // Adding the external memory for `mapping'.
00167   n += mapping.capacity() * sizeof(std::pair<dimension_type, dimension_type>);
00168   return n;
00169 }

void Parma_Polyhedra_Library::MIP_Problem::swap ( MIP_Problem y  )  [inline]

Swaps *this with y.

Definition at line 121 of file MIP_Problem.inlines.hh.

References base, external_space_dim, first_pending_constraint, i_variables, initialized, input_cs, input_obj_function, internal_space_dim, last_generator, mapping, opt_mode, pricing, status, Parma_Polyhedra_Library::swap(), tableau, and working_cost.

Referenced by clear(), operator=(), and swap().

00121                                 {
00122   std::swap(external_space_dim, y.external_space_dim);
00123   std::swap(internal_space_dim, y.internal_space_dim);
00124   std::swap(tableau, y.tableau);
00125   std::swap(working_cost, y.working_cost);
00126   std::swap(mapping, y.mapping);
00127   std::swap(initialized, y.initialized);
00128   std::swap(base, y.base);
00129   std::swap(status, y.status);
00130   std::swap(pricing, y.pricing);
00131   std::swap(input_cs, y.input_cs);
00132   std::swap(first_pending_constraint, y.first_pending_constraint);
00133   std::swap(input_obj_function, y.input_obj_function);
00134   std::swap(opt_mode, y.opt_mode);
00135   std::swap(last_generator, y.last_generator);
00136   std::swap(i_variables, y.i_variables);
00137 }

MIP_Problem::Control_Parameter_Value Parma_Polyhedra_Library::MIP_Problem::get_control_parameter ( Control_Parameter_Name  name  )  const [inline]

Returns the value of the control parameter name.

Definition at line 109 of file MIP_Problem.inlines.hh.

References pricing, and PRICING.

Referenced by compute_simplex_using_exact_pricing(), process_pending_constraints(), and second_phase().

00109                                                                     {
00110   used(name);
00111   assert(name == PRICING);
00112   return pricing;
00113 }

void Parma_Polyhedra_Library::MIP_Problem::set_control_parameter ( Control_Parameter_Value  value  )  [inline]

Sets control parameter value.

Definition at line 116 of file MIP_Problem.inlines.hh.

References pricing.

00116                                                                 {
00117   pricing = value;
00118 }

bool Parma_Polyhedra_Library::MIP_Problem::process_pending_constraints (  )  [private]

Processes the pending constraints of *this.

Returns:
true if and only if the MIP problem is satisfiable after processing the pending constraints, false otherwise.

Definition at line 629 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Matrix::add_zero_columns(), Parma_Polyhedra_Library::Matrix::add_zero_rows(), base, Parma_Polyhedra_Library::Linear_Expression::coefficient(), Parma_Polyhedra_Library::Constraint::coefficient(), compute_generator(), compute_simplex_using_exact_pricing(), compute_simplex_using_steepest_edge_float(), erase_artificials(), external_space_dim, first_pending_constraint, get_control_parameter(), Parma_Polyhedra_Library::Constraint::inhomogeneous_term(), input_cs, input_obj_function, internal_space_dim, Parma_Polyhedra_Library::Constraint::is_inequality(), last_generator, linear_combine(), mapping, Parma_Polyhedra_Library::MAXIMIZATION, Parma_Polyhedra_Library::MINIMIZATION, Parma_Polyhedra_Library::neg_assign(), Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), OK(), opt_mode, OPTIMIZED, parse_constraints(), PRICING, PRICING_STEEPEST_EDGE_FLOAT, SATISFIABLE, Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::Linear_Expression::space_dimension(), space_dimension(), Parma_Polyhedra_Library::Constraint::space_dimension(), status, tableau, UNBOUNDED, UNSATISFIABLE, and working_cost.

Referenced by is_lp_satisfiable().

00629                                             {
00630   const dimension_type num_original_rows = tableau.num_rows();
00631   dimension_type new_rows = 0;
00632   dimension_type new_slacks = 0;
00633   dimension_type new_var_columns = 0;
00634   std::deque<bool> is_tableau_constraint;
00635   std::deque<bool> nonnegative_variable;
00636   std::vector<dimension_type> unfeasible_tableau_rows;
00637   std::deque<bool> satisfied_ineqs;
00638   // Check the new constraints to adjust the data structures.
00639   // If `false' is returned, the pending constraints are trivially
00640   // unfeasible.
00641   if (!parse_constraints(new_rows, new_slacks, is_tableau_constraint,
00642                          nonnegative_variable, unfeasible_tableau_rows,
00643                          satisfied_ineqs)) {
00644     status = UNSATISFIABLE;
00645     return false;
00646   };
00647 
00648   const dimension_type first_free_tableau_index = tableau.num_columns()-1;
00649 
00650   if (external_space_dim > internal_space_dim) {
00651     const dimension_type space_diff = external_space_dim - internal_space_dim;
00652     for (dimension_type i = 0, j = 0; i < space_diff; ++i, ++j) {
00653       // Set `mapping' properly to store that every variable is split.
00654       // In the following case the value of the original variable can be
00655       // negative.
00656       if (!nonnegative_variable[internal_space_dim+i]) {
00657         mapping.push_back(std::make_pair(first_free_tableau_index+j,
00658                                          first_free_tableau_index+1+j));
00659         ++j;
00660         new_var_columns += 2;
00661       }
00662       // The variable is nonnegative.
00663       else {
00664         mapping.push_back(std::make_pair(first_free_tableau_index+j, 0));
00665         ++new_var_columns;
00666       }
00667     }
00668   }
00669 
00670   // Resize the tableau and adding the necessary columns for artificial and
00671   // slack variables.
00672   dimension_type num_satisfied_ineqs = std::count(satisfied_ineqs.begin(),
00673                                                   satisfied_ineqs.end(),
00674                                                   true);
00675   const dimension_type unfeasible_tableau_rows_size
00676     = unfeasible_tableau_rows.size();
00677   const dimension_type artificial_cols
00678     = new_rows + unfeasible_tableau_rows_size - num_satisfied_ineqs;
00679   const dimension_type new_total_columns
00680     = new_var_columns + new_slacks + artificial_cols;
00681   if (new_rows > 0)
00682     tableau.add_zero_rows(new_rows, Row::Flags());
00683   if (new_total_columns > 0)
00684     tableau.add_zero_columns(new_total_columns);
00685   dimension_type tableau_num_rows = tableau.num_rows();
00686 
00687   // The following vector will be useful know if a constraint is feasible
00688   // and does not require an additional artificial variable.
00689   std::deque<bool> worked_out_row (tableau_num_rows, false);
00690   dimension_type tableau_num_columns = tableau.num_columns();
00691 
00692   // Sync the `base' vector size to the new tableau: fill with zeros to encode
00693   // that these rows are not OK and must be adjusted.
00694   base.insert(base.end(), new_rows, 0);
00695   const dimension_type base_size = base.size();
00696 
00697   // These indexes will be used to insert slack and artificial variables.
00698   dimension_type slack_index = tableau_num_columns - artificial_cols - 1;
00699   dimension_type artificial_index = slack_index;
00700 
00701   // The first column index of the tableau that contains an
00702   // artificial variable. Encode with 0 the fact the there are not
00703   // artificial variables.
00704   const dimension_type begin_artificials = artificial_cols > 0
00705     ? artificial_index : 0;
00706 
00707   // Proceed with the insertion of the constraints.
00708   for (dimension_type k = tableau_num_rows, i = input_cs.size();
00709        i-- > first_pending_constraint;  )
00710     if (is_tableau_constraint[i]) {
00711       // Copy the original constraint in the tableau.
00712       Row& tableau_k = tableau[--k];
00713       const Constraint& cs_i = input_cs[i];
00714       for (dimension_type sd = cs_i.space_dimension(); sd-- > 0; ) {
00715         tableau_k[mapping[sd+1].first] = cs_i.coefficient(Variable(sd));
00716         // Split if needed.
00717         if (mapping[sd+1].second != 0)
00718           neg_assign(tableau_k[mapping[sd+1].second],
00719                      tableau_k[mapping[sd+1].first]);
00720       }
00721       tableau_k[mapping[0].first] = cs_i.inhomogeneous_term();
00722       // Split if needed.
00723       if (mapping[0].second != 0)
00724         tableau_k[mapping[0].second] = -cs_i.inhomogeneous_term();
00725 
00726       // Add the slack variable, if needed.
00727       if (cs_i.is_inequality()) {
00728         tableau_k[--slack_index] = -1;
00729         // If the constraint is already satisfied, we will not use artificial
00730         // variables to compute a feasible base: this to speed up
00731         // the algorithm.
00732           if (satisfied_ineqs[i]) {
00733             base[k] = slack_index;
00734             worked_out_row[k] = true;
00735           }
00736       }
00737       for (dimension_type j = base_size; j-- > 0; )
00738         if (k != j && tableau_k[base[j]] != 0 && base[j] != 0)
00739           linear_combine(tableau_k, tableau[j], base[j]);
00740     }
00741 
00742   // We negate the row if tableau[i][0] <= 0 to get the inhomogeneous term > 0.
00743   // This simplifies the insertion of the artificial variables: the value of
00744   // each artificial variable will be 1.
00745   for (dimension_type i = tableau_num_rows; i-- > 0 ; ) {
00746     Row& tableau_i = tableau[i];
00747     if (tableau_i[0] > 0)
00748       for (dimension_type j = tableau_num_columns; j-- > 0; )
00749         neg_assign(tableau_i[j]);
00750   }
00751 
00752   // Set the working cost function with the right size.
00753   working_cost = Row(tableau_num_columns, Row::Flags());
00754 
00755   // Insert artificial variables for the nonfeasible constraints.
00756   for (dimension_type i = 0; i < unfeasible_tableau_rows_size; ++i) {
00757     tableau[unfeasible_tableau_rows[i]][artificial_index] = 1;
00758     working_cost[artificial_index] = -1;
00759     base[unfeasible_tableau_rows[i]] = artificial_index;
00760     ++artificial_index;
00761   }
00762 
00763   // Modify the tableau and the new cost function by adding
00764   // the artificial variables (which enter the base). Note that if an
00765   // inequality was satisfied by `last_generator', this will be not processed.
00766   // This information in encoded in `worked_out_row'.
00767   // As for the cost function, all the artificial variables should have
00768   // coefficient -1.
00769   for (dimension_type i = num_original_rows; i < tableau_num_rows; ++i) {
00770     if (worked_out_row[i])
00771       continue;
00772     tableau[i][artificial_index] = 1;
00773     working_cost[artificial_index] = -1;
00774     base[i] = artificial_index;
00775     ++artificial_index;
00776   }
00777   // The last column index of the tableau containing an artificial variable.
00778   const dimension_type end_artificials = artificial_index - 1;
00779 
00780   // Set the extra-coefficient of the cost functions to record its sign.
00781   // This is done to keep track of the possible sign's inversion.
00782   const dimension_type last_obj_index = working_cost.size() - 1;
00783   working_cost[last_obj_index] = 1;
00784 
00785   // Express the problem in terms of the variables in base.
00786   for (dimension_type i = tableau_num_rows; i-- > 0; )
00787     if (working_cost[base[i]] != 0)
00788       linear_combine(working_cost, tableau[i], base[i]);
00789 
00790   // Deal with zero dimensional problems.
00791   if (space_dimension() == 0) {
00792     status = OPTIMIZED;
00793     last_generator = point();
00794     return true;
00795   }
00796   // Deal with trivial cases.
00797   // If there is no constraint in the tableau, then the feasible region
00798   // is only delimited by non-negativity constraints. Therefore,
00799   // the problem is unbounded as soon as the cost function has
00800   // a variable with a positive coefficient.
00801   if (tableau_num_rows == 0) {
00802     const dimension_type input_obj_function_size
00803       = input_obj_function.space_dimension();
00804     for (dimension_type i = input_obj_function_size; i-- > 0; )
00805       // If a the value of a variable in the objective function is
00806       // different from zero, the final status is unbounded.
00807       // In the first part the variable is constrained to be greater or equal
00808       // than zero.
00809       if ((((input_obj_function.coefficient(Variable(i)) > 0
00810              && opt_mode == MAXIMIZATION)
00811             || (input_obj_function.coefficient(Variable(i)) < 0
00812                 && opt_mode == MINIMIZATION)) && mapping[i].second == 0)
00813           // In the following case the variable is unconstrained.
00814           || (input_obj_function.coefficient(Variable(i)) != 0
00815               && mapping[i].second != 0)) {
00816         // Ensure the right space dimension is obtained.
00817         last_generator = point(0 * Variable(space_dimension()-1));
00818         status = UNBOUNDED;
00819         return true;
00820       }
00821 
00822     // The problem is neither trivially unfeasible nor trivially unbounded.
00823     // The tableau was successful computed and the caller has to figure
00824     // out which case applies.
00825     status = OPTIMIZED;
00826     // Ensure the right space dimension is obtained.
00827     last_generator = point(0*Variable(space_dimension()-1));
00828     assert(OK());
00829     return true;
00830   }
00831 
00832   // Now we are ready to solve the first phase.
00833   bool first_phase_succesful
00834     = (get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_FLOAT)
00835     ? compute_simplex_using_steepest_edge_float()
00836     : compute_simplex_using_exact_pricing();
00837 
00838 #if PPL_NOISY_SIMPLEX
00839   std::cout << "MIP_Problem::solve: 1st phase ended at iteration "
00840             << num_iterations << "." << std::endl;
00841 #endif
00842 
00843   if (!first_phase_succesful || working_cost[0] != 0) {
00844     // The feasible region is empty.
00845     status = UNSATISFIABLE;
00846     return false;
00847   }
00848 
00849   // Prepare *this for a possible second phase.
00850   if (begin_artificials != 0)
00851     erase_artificials(begin_artificials, end_artificials);
00852   compute_generator();
00853   status = SATISFIABLE;
00854   assert(OK());
00855   return true;
00856 }

void Parma_Polyhedra_Library::MIP_Problem::second_phase (  )  [private]

Optimizes the MIP problem using the second phase of the primal simplex algorithm.

Definition at line 1396 of file MIP_Problem.cc.

References base, Parma_Polyhedra_Library::Linear_Expression::coefficient(), compute_generator(), compute_simplex_using_exact_pricing(), compute_simplex_using_steepest_edge_float(), get_control_parameter(), Parma_Polyhedra_Library::Linear_Expression::inhomogeneous_term(), input_obj_function, linear_combine(), mapping, Parma_Polyhedra_Library::MINIMIZATION, Parma_Polyhedra_Library::neg_assign(), Parma_Polyhedra_Library::Matrix::num_rows(), OK(), opt_mode, OPTIMIZED, PRICING, PRICING_STEEPEST_EDGE_FLOAT, SATISFIABLE, Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::Linear_Expression::space_dimension(), status, Parma_Polyhedra_Library::Row::swap(), tableau, UNBOUNDED, and working_cost.

Referenced by solve(), and solve_mip().

01396                              {
01397   // Second_phase requires that *this is satisfiable.
01398   assert(status == SATISFIABLE || status == UNBOUNDED || status == OPTIMIZED);
01399   // In the following cases the problem is already solved.
01400   if (status == UNBOUNDED || status == OPTIMIZED)
01401     return;
01402 
01403   // Build the objective function for the second phase.
01404   const dimension_type input_obj_function_sd
01405     = input_obj_function.space_dimension();
01406   Row new_cost(input_obj_function_sd + 1, Row::Flags());
01407   for (dimension_type i = input_obj_function_sd; i-- > 0; )
01408     new_cost[i+1] = input_obj_function.coefficient(Variable(i));
01409   new_cost[0] = input_obj_function.inhomogeneous_term();
01410 
01411   // Negate the cost function if we are minimizing.
01412   if (opt_mode == MINIMIZATION)
01413     for (dimension_type i = new_cost.size(); i-- > 0; )
01414       neg_assign(new_cost[i]);
01415 
01416   // Substitute properly the cost function in the `costs' matrix.
01417   const dimension_type cost_zero_size = working_cost.size();
01418   Row tmp_cost = Row(new_cost, cost_zero_size, cost_zero_size);
01419   tmp_cost.swap(working_cost);
01420   working_cost[cost_zero_size-1] = 1;
01421 
01422   // Split the variables the cost function.
01423   for (dimension_type i = new_cost.size(); i-- > 1; ) {
01424     const dimension_type original_var = mapping[i].first;
01425     const dimension_type split_var = mapping[i].second;
01426     working_cost[original_var] = new_cost[i];
01427     if (mapping[i].second != 0)
01428       working_cost[split_var] = - new_cost[i];
01429   }
01430   // Here the first phase problem succeeded with optimum value zero.
01431   // Express the old cost function in terms of the computed base.
01432   for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
01433     const dimension_type base_i = base[i];
01434     if (working_cost[base_i] != 0)
01435       linear_combine(working_cost, tableau[i], base_i);
01436   }
01437   // Solve the second phase problem.
01438   bool second_phase_successful
01439     = (get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_FLOAT)
01440     ? compute_simplex_using_steepest_edge_float()
01441     : compute_simplex_using_exact_pricing();
01442   compute_generator();
01443 #if PPL_NOISY_SIMPLEX
01444   std::cout << "MIP_Problem::solve: 2nd phase ended at iteration "
01445             << num_iterations << "." << std::endl;
01446 #endif
01447   status = second_phase_successful ? OPTIMIZED : UNBOUNDED;
01448   assert(OK());
01449 }

MIP_Problem_Status Parma_Polyhedra_Library::MIP_Problem::compute_tableau ( std::vector< dimension_type > &  worked_out_row  )  [private]

Assigns to this->tableau a simplex tableau representing the MIP problem, inserting into this->mapping the information that is required to recover the original MIP problem.

Returns:
UNFEASIBLE_MIP_PROBLEM if the constraint system contains any trivially unfeasible constraint (tableau was not computed); UNBOUNDED_MIP_PROBLEM if the problem is trivially unbounded (the computed tableau contains no constraints); OPTIMIZED_MIP_PROBLEM> if the problem is neither trivially unfeasible nor trivially unbounded (the tableau was computed successfully).

bool Parma_Polyhedra_Library::MIP_Problem::parse_constraints ( dimension_type new_num_rows,
dimension_type num_slack_variables,
std::deque< bool > &  is_tableau_constraint,
std::deque< bool > &  nonnegative_variable,
std::vector< dimension_type > &  unfeasible_tableau_rows,
std::deque< bool > &  satisfied_ineqs 
) [private]

Parses the pending constraints to gather information on how to resize the tableau.

Returns:
UNSATISFIABLE if is detected a trivially false constraint, SATISFIABLE otherwise.
Parameters:
new_num_rows This will store the number of rows that has to be added to the original tableau.
num_slack_variables This will store the number of slack variables that has to be added to the original tableau.
is_tableau_constraint Every element of this vector will be set to true if the associated pending constraint has to be inserted in the tableau, false otherwise.
nonnegative_variable This will encode for each variable if this one was split or not. Every element of this vector will be set to true if the associated variable is split, false otherwise.
unfeasible_tableau_rows This will contain all the row indexes of the tableau that are no more satisfied after adding more constraints to *this.
satisfied_ineqs This will contain all the row indexes of the tableau that are already satisfied by `last_generator' and do not require artificial variables to have a starting feasible base.

Definition at line 452 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Constraint::coefficient(), external_space_dim, first_pending_constraint, Parma_Polyhedra_Library::Constraint::inhomogeneous_term(), input_cs, Parma_Polyhedra_Library::Constraint::is_equality(), Parma_Polyhedra_Library::Constraint::is_inequality(), is_satisfied(), last_generator, mapping, merge_split_variables(), and Parma_Polyhedra_Library::Constraint::space_dimension().

Referenced by process_pending_constraints().

00458                                                                      {
00459   satisfied_ineqs.clear();
00460   satisfied_ineqs.insert(satisfied_ineqs.end(), input_cs.size(),
00461                          false);
00462 
00463   const dimension_type cs_num_rows = input_cs.size();
00464   const dimension_type cs_space_dim = external_space_dim;
00465 
00466   // Step 1:
00467   // determine variables that are constrained to be nonnegative,
00468   // detect (non-negativity or tautology) constraints that will not
00469   // be part of the tableau and count the number of slack variables.
00470 
00471   // Counters determining the dimensions of the tableau:
00472   // initialized here, they will be updated while examining `cs'.
00473   tableau_num_rows = cs_num_rows;
00474   dimension_type tableau_num_cols = 2*cs_space_dim;
00475   num_slack_variables = 0;
00476 
00477   // On exit, `is_tableau_constraint[i]' will be true if and only if
00478   // `cs[i]' is neither a tautology (e.g., 1 >= 0) nor a non-negativity
00479   // constraint (e.g., X >= 0).
00480   is_tableau_constraint = std::deque<bool> (cs_num_rows, true);
00481 
00482   // On exit, `nonnegative_variable[j]' will be true if and only if
00483   // Variable(j) is bound to be nonnegative in `cs'.
00484   nonnegative_variable = std::deque<bool> (cs_space_dim, false);
00485 
00486   // Check for already known information about space dimensions and
00487   // store them in `nonnegative_variable'.
00488   const dimension_type mapping_size = mapping.size();
00489   for (dimension_type i = std::min(mapping_size, cs_space_dim+1); i-- > 1; )
00490     if (mapping[i].second == 0) {
00491       nonnegative_variable[i-1] = true;
00492       --tableau_num_cols;
00493     }
00494 
00495   // Process each row of the `cs' matrix.
00496   for (dimension_type i = cs_num_rows; i-- > first_pending_constraint; ) {
00497     const Constraint& cs_i = input_cs[i];
00498     bool found_a_nonzero_coeff = false;
00499     bool found_many_nonzero_coeffs = false;
00500     dimension_type nonzero_coeff_column_index = 0;
00501     for (dimension_type sd = cs_i.space_dimension(); sd-- > 0; ) {
00502       if (cs_i.coefficient(Variable(sd)) != 0) {
00503         if (found_a_nonzero_coeff) {
00504           found_many_nonzero_coeffs = true;
00505           if (cs_i.is_inequality())
00506             ++num_slack_variables;
00507           break;
00508         }
00509         else {
00510           nonzero_coeff_column_index = sd + 1;
00511           found_a_nonzero_coeff = true;
00512         }
00513       }
00514     }
00515     // If more than one coefficient is nonzero,
00516     // continue with next constraint.
00517     if (found_many_nonzero_coeffs) {
00518       // CHECKME: Is it true that in the first phase we can apply
00519       // `is_satisfied()' with the generator `point()'?  If so, the following
00520       // code works even if we do not have a feasible point.
00521       // Check for satisfiability of the inequality. This can be done if we
00522       // have a feasible point of *this.
00523       if (cs_i.is_inequality() && is_satisfied(cs_i, last_generator))
00524         satisfied_ineqs[i] = true;
00525       continue;
00526     }
00527 
00528     if (!found_a_nonzero_coeff) {
00529       // All coefficients are 0.
00530       // The constraint is either trivially true or trivially false.
00531       if (cs_i.is_inequality()) {
00532         if (cs_i.inhomogeneous_term() < 0)
00533           // A constraint such as -1 >= 0 is trivially false.
00534           return false;
00535       }
00536       else
00537         // The constraint is an equality.
00538         if (cs_i.inhomogeneous_term() != 0)
00539           // A constraint such as 1 == 0 is trivially false.
00540           return false;
00541       // Here the constraint is trivially true.
00542       is_tableau_constraint[i] = false;
00543       --tableau_num_rows;
00544       continue;
00545     }
00546     else {
00547       // Here we have only one nonzero coefficient.
00548       /*
00549 
00550       We have the following methods:
00551       A) Do split the variable and do add the constraint
00552          in the tableau.
00553       B) Do not split the variable and do add the constraint
00554          in the tableau.
00555       C) Do not split the variable and do not add the constraint
00556          in the tableau.
00557 
00558       Let the constraint be (a*v + b relsym 0).
00559       These are the 12 possible combinations we can have:
00560                 a |  b | relsym | method
00561       ----------------------------------
00562       1)       >0 | >0 |   >=   |   A
00563       2)       >0 | >0 |   ==   |   A
00564       3)       <0 | <0 |   >=   |   A
00565       4)       >0 | =0 |   ==   |   B
00566       5)       >0 | <0 |   ==   |   B
00567       Note:    <0 | >0 |   ==   | impossible by strong normalization
00568       Note:    <0 | =0 |   ==   | impossible by strong normalization
00569       Note:    <0 | <0 |   ==   | impossible by strong normalization
00570       6)       >0 | <0 |   >=   |   B
00571       7)       >0 | =0 |   >=   |   C
00572       8)       <0 | >0 |   >=   |   A
00573       9)       <0 | =0 |   >=   |   A
00574 
00575       The next lines will apply the correct method to each case.
00576       */
00577 
00578       // The variable index is not equal to the column index.
00579       const dimension_type nonzero_var_index = nonzero_coeff_column_index - 1;
00580 
00581       const int sgn_a
00582         = sgn(cs_i.coefficient(Variable(nonzero_coeff_column_index-1)));
00583       const int sgn_b = sgn(cs_i.inhomogeneous_term());
00584       // Cases 1-3: apply method A.
00585       if (sgn_a == sgn_b) {
00586         if (cs_i.is_inequality())
00587           ++num_slack_variables;
00588       }
00589       // Cases 4-5: apply method B.
00590       else if (cs_i.is_equality()) {
00591         if (!nonnegative_variable[nonzero_var_index]) {
00592           nonnegative_variable[nonzero_var_index] = true;
00593           --tableau_num_cols;
00594         }
00595       }
00596       // Case 6: apply method B.
00597       else if (sgn_b < 0) {
00598         if (!nonnegative_variable[nonzero_var_index]) {
00599           nonnegative_variable[nonzero_var_index] = true;
00600           --tableau_num_cols;
00601         }
00602         ++num_slack_variables;
00603       }
00604       // Case 7: apply method C.
00605       else if (sgn_a > 0) {
00606         // This is the most important case in the incrementality solving:
00607         // merge two variables.
00608         if (!nonnegative_variable[nonzero_var_index]) {
00609           nonnegative_variable[nonzero_var_index] = true;
00610           --tableau_num_cols;
00611           if (nonzero_coeff_column_index < mapping_size)
00612             merge_split_variables(nonzero_coeff_column_index,
00613                                   unfeasible_tableau_rows);
00614           is_tableau_constraint[i] = false;
00615         }
00616         else
00617           is_tableau_constraint[i] = false;
00618         --tableau_num_rows;
00619       }
00620       // Cases 8-9: apply method A.
00621       else
00622         ++num_slack_variables;
00623     }
00624   }
00625   return true;
00626 }

PPL::dimension_type Parma_Polyhedra_Library::MIP_Problem::get_exiting_base_index ( dimension_type  entering_var_index  )  const [private]

Computes the row index of the variable exiting the base of the MIP problem. Implemented with anti-cycling rule.

Returns:
The row index of the variable exiting the base.
Parameters:
entering_var_index The column index of the variable entering the base.

Definition at line 1060 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::abs_assign(), base, Parma_Polyhedra_Library::exact_div_assign(), Parma_Polyhedra_Library::lcm_assign(), Parma_Polyhedra_Library::Matrix::num_rows(), tableau, and TEMP_INTEGER.

Referenced by compute_simplex_using_exact_pricing(), and compute_simplex_using_steepest_edge_float().

01060                                                                        {
01061   // The variable exiting the base should be associated to a tableau
01062   // constraint such that the ratio
01063   // tableau[i][entering_var_index] / tableau[i][base[i]]
01064   // is strictly positive and minimal.
01065 
01066   // Find the first tableau constraint `c' having a positive value for
01067   // tableau[i][entering_var_index] / tableau[i][base[i]]
01068   const dimension_type tableau_num_rows = tableau.num_rows();
01069   dimension_type exiting_base_index = tableau_num_rows;
01070   for (dimension_type i = 0; i < tableau_num_rows; ++i) {
01071     const Row& t_i = tableau[i];
01072     const int num_sign = sgn(t_i[entering_var_index]);
01073     if (num_sign != 0 && num_sign == sgn(t_i[base[i]])) {
01074       exiting_base_index = i;
01075       break;
01076     }
01077   }
01078   // Check for unboundedness.
01079   if (exiting_base_index == tableau_num_rows)
01080     return tableau_num_rows;
01081 
01082   // Reaching this point means that a variable will definitely exit the base.
01083   TEMP_INTEGER(lcm);
01084   TEMP_INTEGER(current_min);
01085   TEMP_INTEGER(challenger);
01086   for (dimension_type i = exiting_base_index + 1; i < tableau_num_rows; ++i) {
01087     const Row& t_i = tableau[i];
01088     const Coefficient& t_ie = t_i[entering_var_index];
01089     const Coefficient& t_ib = t_i[base[i]];
01090     const int t_ie_sign = sgn(t_ie);
01091     if (t_ie_sign != 0 && t_ie_sign == sgn(t_ib)) {
01092       const Row& t_e = tableau[exiting_base_index];
01093       const Coefficient& t_ee = t_e[entering_var_index];
01094       lcm_assign(lcm, t_ee, t_ie);
01095       exact_div_assign(current_min, lcm, t_ee);
01096       current_min *= t_e[0];
01097       abs_assign(current_min);
01098       exact_div_assign(challenger, lcm, t_ie);
01099       challenger *= t_i[0];
01100       abs_assign(challenger);
01101       current_min -= challenger;
01102       const int sign = sgn(current_min);
01103       if (sign > 0
01104           || (sign == 0 && base[i] < base[exiting_base_index]))
01105         exiting_base_index = i;
01106     }
01107   }
01108   return exiting_base_index;
01109 }

void Parma_Polyhedra_Library::MIP_Problem::linear_combine ( Row x,
const Row y,
const dimension_type  k 
) [static, private]

Linearly combines x with y so that *this[k] is 0.

Parameters:
x The Row that will be combined with y object.
y The Row that will be combined with x object.
k The position of *this that have to be $0$.
Computes a linear combination of x and y having the element of index k equal to $0$. Then it assigns the resulting Linear_Row to x and normalizes it.

Definition at line 1012 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Row::normalize(), Parma_Polyhedra_Library::normalize2(), Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::sub_mul_assign(), and TEMP_INTEGER.

Referenced by pivot(), process_pending_constraints(), and second_phase().

01014                                                          {
01015   assert(x.size() == y.size());
01016   assert(y[k] != 0 && x[k] != 0);
01017   // Let g be the GCD between `x[k]' and `y[k]'.
01018   // For each i the following computes
01019   //   x[i] = x[i]*y[k]/g - y[i]*x[k]/g.
01020   TEMP_INTEGER(normalized_x_k);
01021   TEMP_INTEGER(normalized_y_k);
01022   normalize2(x[k], y[k], normalized_x_k, normalized_y_k);
01023   for (dimension_type i = x.size(); i-- > 0; )
01024     if (i != k) {
01025       Coefficient& x_i = x[i];
01026       x_i *= normalized_y_k;
01027 #if 1 // CHECKME: the test seems to speed up the GMP computation.
01028       const Coefficient& y_i = y[i];
01029       if (y_i != 0)
01030         sub_mul_assign(x_i, y_i, normalized_x_k);
01031 #else
01032       sub_mul_assign(x_i, y[i], normalized_x_k);
01033 #endif // 1
01034     }
01035   x[k] = 0;
01036   x.normalize();
01037 }

void Parma_Polyhedra_Library::MIP_Problem::pivot ( dimension_type  entering_var_index,
dimension_type  exiting_base_index 
) [private]

Performs the pivoting operation on the tableau.

Parameters:
entering_var_index The index of the variable entering the base.
exiting_base_index The index of the row exiting the base.

Definition at line 1041 of file MIP_Problem.cc.

References base, linear_combine(), Parma_Polyhedra_Library::Matrix::num_rows(), tableau, and working_cost.

Referenced by compute_simplex_using_exact_pricing(), compute_simplex_using_steepest_edge_float(), and erase_artificials().

01042                                                                  {
01043   const Row& tableau_out = tableau[exiting_base_index];
01044   // Linearly combine the constraints.
01045   for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
01046     Row& tableau_i = tableau[i];
01047     if (i != exiting_base_index && tableau_i[entering_var_index] != 0)
01048       linear_combine(tableau_i, tableau_out, entering_var_index);
01049   }
01050   // Linearly combine the cost function.
01051   if (working_cost[entering_var_index] != 0)
01052     linear_combine(working_cost, tableau_out, entering_var_index);
01053   // Adjust the base.
01054   base[exiting_base_index] = entering_var_index;
01055 }

PPL::dimension_type Parma_Polyhedra_Library::MIP_Problem::textbook_entering_index (  )  const [private]

Computes the column index of the variable entering the base, using the textbook algorithm with anti-cycling rule.

Returns:
The column index of the variable that enters the base. If no such variable exists, optimality was achieved and 0 is returned.

Definition at line 993 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Row::size(), and working_cost.

Referenced by compute_simplex_using_exact_pricing(), and compute_simplex_using_steepest_edge_float().

00993                                               {
00994   // The variable entering the base is the first one whose coefficient
00995   // in the cost function has the same sign the cost function itself.
00996   // If no such variable exists, then we met the optimality condition
00997   // (and return 0 to the caller).
00998 
00999   // Get the "sign" of the cost function.
01000   const dimension_type cost_sign_index = working_cost.size() - 1;
01001   const int cost_sign = sgn(working_cost[cost_sign_index]);
01002   assert(cost_sign != 0);
01003   for (dimension_type i = 1; i < cost_sign_index; ++i)
01004     if (sgn(working_cost[i]) == cost_sign)
01005       return i;
01006   // No variable has to enter the base:
01007   // the cost function was optimized.
01008   return 0;
01009 }

PPL::dimension_type Parma_Polyhedra_Library::MIP_Problem::steepest_edge_exact_entering_index (  )  const [private]

Computes the column index of the variable entering the base, using an exact steepest-edge algorithm with anti-cycling rule.

Returns:
The column index of the variable that enters the base. If no such variable exists, optimality was achieved and 0 is returned.
To compute the entering_index, the steepest edge algorithm chooses the index `j' such that $\frac{d_{j}}{\|\Delta x^{j} \|}$ is the largest in absolute value, where

\[ \|\Delta x^{j} \| = \left( 1+\sum_{i=1}^{m} \alpha_{ij}^2 \right)^{\frac{1}{2}}. \]

Recall that, due to the exact integer implementation of the algorithm, our tableau doesn't contain the ``real'' $\alpha$ values, but these can be computed dividing the value of the coefficient by the value of the variable in base. Obviously the result may not be an integer, so we will proceed in another way: we compute the lcm of all the variables in base to get the good ``weight'' of each Coefficient of the tableau.

Definition at line 920 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::add_mul_assign(), base, Parma_Polyhedra_Library::exact_div_assign(), Parma_Polyhedra_Library::lcm_assign(), Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::swap(), tableau, TEMP_INTEGER, and working_cost.

Referenced by compute_simplex_using_exact_pricing().

00920                                                          {
00921   const dimension_type tableau_num_rows = tableau.num_rows();
00922   assert(tableau_num_rows == base.size());
00923   // The square of the lcm of all the coefficients of variables in base.
00924   TEMP_INTEGER(squared_lcm_basis);
00925   // The normalization factor for each coefficient in the tableau.
00926   std::vector<Coefficient> norm_factor(tableau_num_rows);
00927   {
00928     // Compute the lcm of all the coefficients of variables in base.
00929     TEMP_INTEGER(lcm_basis);
00930     lcm_basis = 1;
00931     for (dimension_type i = tableau_num_rows; i-- > 0; )
00932       lcm_assign(lcm_basis, lcm_basis, tableau[i][base[i]]);
00933     // Compute normalization factors.
00934     for (dimension_type i = tableau_num_rows; i-- > 0; )
00935       exact_div_assign(norm_factor[i], lcm_basis, tableau[i][base[i]]);
00936     // Compute the square of `lcm_basis', exploiting the fact that
00937     // `lcm_basis' will no longer be needed.
00938     lcm_basis *= lcm_basis;
00939     std::swap(squared_lcm_basis, lcm_basis);
00940   }
00941 
00942   // Defined here to avoid repeated (de-)allocations.
00943   TEMP_INTEGER(challenger_num);
00944   TEMP_INTEGER(scalar_value);
00945   TEMP_INTEGER(challenger_den);
00946   TEMP_INTEGER(challenger_value);
00947   TEMP_INTEGER(current_value);
00948 
00949   TEMP_INTEGER(current_num);
00950   TEMP_INTEGER(current_den);
00951   dimension_type entering_index = 0;
00952   const int cost_sign = sgn(working_cost[working_cost.size() - 1]);
00953   for (dimension_type j = tableau.num_columns() - 1; j-- > 1; ) {
00954     const Coefficient& cost_j = working_cost[j];
00955     if (sgn(cost_j) == cost_sign) {
00956       // We cannot compute the (exact) square root of abs(\Delta x_j).
00957       // The workaround is to compute the square of `cost[j]'.
00958       challenger_num = cost_j * cost_j;
00959       // Due to our integer implementation, the `1' term in the denominator
00960       // of the original formula has to be replaced by `squared_lcm_basis'.
00961       challenger_den = squared_lcm_basis;
00962       for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00963         const Coefficient& tableau_ij = tableau[i][j];
00964         // FIXME: the test seems to speed up the GMP computation.
00965         if (tableau_ij != 0) {
00966           scalar_value = tableau_ij * norm_factor[i];
00967           add_mul_assign(challenger_den, scalar_value, scalar_value);
00968         }
00969       }
00970       // Initialization during the first loop.
00971       if (entering_index == 0) {
00972         std::swap(current_num, challenger_num);
00973         std::swap(current_den, challenger_den);
00974         entering_index = j;
00975         continue;
00976       }
00977       challenger_value = challenger_num * current_den;
00978       current_value = current_num * challenger_den;
00979       // Update the values, if the challenger wins.
00980       if (challenger_value > current_value) {
00981         std::swap(current_num, challenger_num);
00982         std::swap(current_den, challenger_den);
00983         entering_index = j;
00984       }
00985     }
00986   }
00987   return entering_index;
00988 }

PPL::dimension_type Parma_Polyhedra_Library::MIP_Problem::steepest_edge_float_entering_index (  )  const [private]

Same as steepest_edge_exact_entering_index, but using floating points.

Note:
Due to rounding errors, the index of the variable entering the base of the MIP problem is not predictable across different architectures. Hence, the overall simplex computation may differ in the path taken to reach the optimum. Anyway, the exact final result will be computed for the MIP_Problem.

Definition at line 875 of file MIP_Problem.cc.

References assign(), Parma_Polyhedra_Library::assign_r(), base, Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), Parma_Polyhedra_Library::Row::size(), tableau, and working_cost.

Referenced by compute_simplex_using_steepest_edge_float().

00875                                                          {
00876   DIRTY_TEMP0(mpq_class, real_coeff);
00877   const dimension_type tableau_num_rows = tableau.num_rows();
00878   assert(tableau_num_rows == base.size());
00879   double challenger_num = 0.0;
00880   double challenger_den = 0.0;
00881   double current_value = 0.0;
00882   dimension_type entering_index = 0;
00883   const int cost_sign = sgn(working_cost[working_cost.size() - 1]);
00884   for (dimension_type j = tableau.num_columns() - 1; j-- > 1; ) {
00885     const Coefficient& cost_j = working_cost[j];
00886     if (sgn(cost_j) == cost_sign) {
00887       // We cannot compute the (exact) square root of abs(\Delta x_j).
00888       // The workaround is to compute the square of `cost[j]'.
00889       assign(challenger_num, cost_j);
00890       challenger_num = fabs(challenger_num);
00891       // Due to our integer implementation, the `1' term in the denominator
00892       // of the original formula has to be replaced by `squared_lcm_basis'.
00893       challenger_den = 1.0;
00894       for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00895         const Row& tableau_i = tableau[i];
00896         const Coefficient& tableau_ij = tableau_i[j];
00897         if (tableau_ij != 0) {
00898           assert(tableau_i[base[i]] != 0);
00899           assign_r(real_coeff.get_num(), tableau_ij, ROUND_NOT_NEEDED);
00900           assign_r(real_coeff.get_den(), tableau_i[base[i]], ROUND_NOT_NEEDED);
00901           real_coeff.canonicalize();
00902           double float_tableau_value;
00903           assign(float_tableau_value, real_coeff);
00904           challenger_den += float_tableau_value * float_tableau_value;
00905         }
00906       }
00907       double challenger_value = sqrt(challenger_den);
00908       // Initialize `current_value' during the first iteration.
00909       // Otherwise update if the challenger wins.
00910       if (entering_index == 0 || challenger_value > current_value) {
00911         current_value = challenger_value;
00912         entering_index = j;
00913       }
00914     }
00915   }
00916   return entering_index;
00917 }

bool Parma_Polyhedra_Library::MIP_Problem::compute_simplex_using_exact_pricing (  )  [private]

Returns true if and if only the algorithm successfully computed a feasible solution.

Note:
Uses an exact pricing method (either textbook or exact steepest edge), so that the result is deterministic across different architectures.

Definition at line 1202 of file MIP_Problem.cc.

References get_control_parameter(), get_exiting_base_index(), Parma_Polyhedra_Library::maybe_abandon(), Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), pivot(), PRICING, PRICING_STEEPEST_EDGE_EXACT, PRICING_TEXTBOOK, Parma_Polyhedra_Library::Row::size(), steepest_edge_exact_entering_index(), tableau, textbook_entering_index(), and working_cost.

Referenced by process_pending_constraints(), and second_phase().

01202                                                     {
01203   assert(tableau.num_columns() == working_cost.size());
01204   assert(get_control_parameter(PRICING) == PRICING_STEEPEST_EDGE_EXACT
01205          || get_control_parameter(PRICING) == PRICING_TEXTBOOK);
01206 
01207   const dimension_type tableau_num_rows = tableau.num_rows();
01208   const bool textbook_pricing
01209     = (PRICING_TEXTBOOK == get_control_parameter(PRICING));
01210 
01211   while (true) {
01212     // Choose the index of the variable entering the base, if any.
01213     const dimension_type entering_var_index
01214       = textbook_pricing
01215       ? textbook_entering_index()
01216       : steepest_edge_exact_entering_index();
01217     // If no entering index was computed, the problem is solved.
01218     if (entering_var_index == 0)
01219       return true;
01220 
01221     // Choose the index of the row exiting the base.
01222     const dimension_type exiting_base_index
01223       = get_exiting_base_index(entering_var_index);
01224     // If no exiting index was computed, the problem is unbounded.
01225     if (exiting_base_index == tableau_num_rows)
01226       return false;
01227 
01228     // Check if the client has requested abandoning all expensive
01229     // computations. If so, the exception specified by the client
01230     // is thrown now.
01231     maybe_abandon();
01232 
01233     // We have not reached the optimality or unbounded condition:
01234     // compute the new base and the corresponding vertex of the
01235     // feasible region.
01236     pivot(entering_var_index, exiting_base_index);
01237 #if PPL_NOISY_SIMPLEX
01238     ++num_iterations;
01239     if (num_iterations % 200 == 0)
01240       std::cout << "Primal Simplex: iteration "
01241                 << num_iterations << "." << std::endl;
01242 #endif
01243   }
01244 }

bool Parma_Polyhedra_Library::MIP_Problem::compute_simplex_using_steepest_edge_float (  )  [private]

Returns true if and if only the algorithm successfully computed a feasible solution.

Note:
Uses a floating point implementation of the steepest edge pricing method, so that the result is correct, but not deterministic across different architectures.

Definition at line 1113 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::abs_assign(), get_exiting_base_index(), Parma_Polyhedra_Library::maybe_abandon(), Parma_Polyhedra_Library::neg_assign(), Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), pivot(), Parma_Polyhedra_Library::Row::size(), steepest_edge_float_entering_index(), tableau, TEMP_INTEGER, textbook_entering_index(), and working_cost.

Referenced by process_pending_constraints(), and second_phase().

01113                                                           {
01114   // We may need to temporarily switch to the textbook pricing.
01115   const unsigned long allowed_non_increasing_loops = 200;
01116   unsigned long non_increased_times = 0;
01117   bool textbook_pricing = false;
01118 
01119   TEMP_INTEGER(cost_sgn_coeff);
01120   TEMP_INTEGER(current_num);
01121   TEMP_INTEGER(current_den);
01122   TEMP_INTEGER(challenger);
01123   TEMP_INTEGER(current);
01124 
01125   cost_sgn_coeff = working_cost[working_cost.size()-1];
01126   current_num = working_cost[0];
01127   if (cost_sgn_coeff < 0)
01128     neg_assign(current_num);
01129   abs_assign(current_den, cost_sgn_coeff);
01130   assert(tableau.num_columns() == working_cost.size());
01131   const dimension_type tableau_num_rows = tableau.num_rows();
01132 
01133   while (true) {
01134     // Choose the index of the variable entering the base, if any.
01135     const dimension_type entering_var_index
01136       = textbook_pricing
01137       ? textbook_entering_index()
01138       : steepest_edge_float_entering_index();
01139 
01140     // If no entering index was computed, the problem is solved.
01141     if (entering_var_index == 0)
01142       return true;
01143 
01144     // Choose the index of the row exiting the base.
01145     const dimension_type exiting_base_index
01146       = get_exiting_base_index(entering_var_index);
01147     // If no exiting index was computed, the problem is unbounded.
01148     if (exiting_base_index == tableau_num_rows)
01149       return false;
01150 
01151     // Check if the client has requested abandoning all expensive
01152     // computations. If so, the exception specified by the client
01153     // is thrown now.
01154     maybe_abandon();
01155 
01156     // We have not reached the optimality or unbounded condition:
01157     // compute the new base and the corresponding vertex of the
01158     // feasible region.
01159     pivot(entering_var_index, exiting_base_index);
01160 
01161     // Now begins the objective function's value check to choose between
01162     // the `textbook' and the float `steepest-edge' technique.
01163     cost_sgn_coeff = working_cost[working_cost.size()-1];
01164 
01165     challenger = working_cost[0];
01166     if (cost_sgn_coeff < 0)
01167       neg_assign(challenger);
01168     challenger *= current_den;
01169     abs_assign(current, cost_sgn_coeff);
01170     current *= current_num;
01171 #if PPL_NOISY_SIMPLEX
01172     ++num_iterations;
01173     if (num_iterations % 200 == 0)
01174       std::cout << "Primal Simplex: iteration "
01175                 << num_iterations << "." << std::endl;
01176 #endif
01177     // If the following condition fails, probably there's a bug.
01178     assert(challenger >= current);
01179     // If the value of the objective function does not improve,
01180     // keep track of that.
01181     if (challenger == current) {
01182       ++non_increased_times;
01183       // In the following case we will proceed using the `textbook'
01184       // technique, until the objective function is not improved.
01185       if (non_increased_times > allowed_non_increasing_loops)
01186         textbook_pricing = true;
01187     }
01188     // The objective function has an improvement:
01189     // reset `non_increased_times' and `textbook_pricing'.
01190     else {
01191       non_increased_times = 0;
01192       textbook_pricing = false;
01193     }
01194     current_num = working_cost[0];
01195     if (cost_sgn_coeff < 0)
01196       neg_assign(current_num);
01197     abs_assign(current_den, cost_sgn_coeff);
01198   }
01199 }

void Parma_Polyhedra_Library::MIP_Problem::erase_artificials ( dimension_type  begin_artificials,
dimension_type  end_artificials 
) [private]

Drop unnecessary artificial variables from the tableau and get ready for the second phase of the simplex algorithm.

Definition at line 1249 of file MIP_Problem.cc.

References base, Parma_Polyhedra_Library::Matrix::erase_to_end(), Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), pivot(), Parma_Polyhedra_Library::Matrix::remove_trailing_columns(), Parma_Polyhedra_Library::Row::shrink(), Parma_Polyhedra_Library::Row::size(), Parma_Polyhedra_Library::Row::swap(), tableau, and working_cost.

Referenced by process_pending_constraints().

01250                                                                           {
01251   const dimension_type tableau_last_index = tableau.num_columns() - 1;
01252   dimension_type tableau_n_rows = tableau.num_rows();
01253   // Step 1: try to remove from the base all the remaining slack variables.
01254   for (dimension_type i = 0; i < tableau_n_rows; ++i)
01255     if (begin_artificials <= base[i] && base[i] <= end_artificials) {
01256       // Search for a non-zero element to enter the base.
01257       Row& tableau_i = tableau[i];
01258       bool redundant = true;
01259       for (dimension_type j = end_artificials+1; j-- > 1; )
01260         if (!(begin_artificials <= j && j <= end_artificials)
01261             && tableau_i[j] != 0) {
01262           pivot(j, i);
01263           redundant = false;
01264           break;
01265         }
01266       if (redundant) {
01267         // No original variable entered the base:
01268         // the constraint is redundant and should be deleted.
01269         --tableau_n_rows;
01270         if (i < tableau_n_rows) {
01271           // Replace the redundant row with the last one,
01272           // taking care of adjusting the iteration index.
01273           tableau_i.swap(tableau[tableau_n_rows]);
01274           base[i] = base[tableau_n_rows];
01275           --i;
01276         }
01277         tableau.erase_to_end(tableau_n_rows);
01278         base.pop_back();
01279       }
01280     }
01281 
01282 
01283   // Step 2: Adjust data structures so as to enter phase 2 of the simplex.
01284 
01285   // Compute the dimensions of the new tableau.
01286   dimension_type num_artificials = 0;
01287   for (dimension_type i = end_artificials + 1; i-- > 1; )
01288     if (begin_artificials <= i && i <= end_artificials) {
01289       ++num_artificials;
01290       tableau.remove_trailing_columns(1);
01291     }
01292 
01293   // Zero the last column of the tableau.
01294   for (dimension_type i = tableau_n_rows; i-- > 0; )
01295     tableau[i][tableau.num_columns()-1] = 0;
01296 
01297   // ... then properly set the element in the (new) last column,
01298   // encoding the kind of optimization; ...
01299   working_cost[tableau.num_columns()-1] = working_cost[tableau_last_index];
01300   // ... and finally remove redundant columns.
01301   const dimension_type working_cost_new_size = working_cost.size() -
01302     num_artificials;
01303   working_cost.shrink(working_cost_new_size);
01304 }

bool Parma_Polyhedra_Library::MIP_Problem::is_in_base ( dimension_type  var_index,
dimension_type row_index 
) const [private]

Definition at line 367 of file MIP_Problem.cc.

References base.

Referenced by compute_generator(), and merge_split_variables().

00368                                                               {
00369   for (row_index = base.size(); row_index-- > 0; )
00370     if (base[row_index] == var_index)
00371       return true;
00372   return false;
00373 }

void Parma_Polyhedra_Library::MIP_Problem::compute_generator (  )  const [private]

Computes a valid generator that satisfies all the constraints of the Linear Programming problem associated to *this.

Definition at line 1308 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::exact_div_assign(), external_space_dim, is_in_base(), last_generator, Parma_Polyhedra_Library::lcm_assign(), mapping, Parma_Polyhedra_Library::neg_assign(), Parma_Polyhedra_Library::sub_mul_assign(), tableau, and TEMP_INTEGER.

Referenced by process_pending_constraints(), and second_phase().

01308                                         {
01309   // We will store in num[] and in den[] the numerators and
01310   // the denominators of every variable of the original problem.
01311   std::vector<Coefficient> num(external_space_dim);
01312   std::vector<Coefficient> den(external_space_dim);
01313   dimension_type row = 0;
01314 
01315   TEMP_INTEGER(lcm);
01316   // Speculatively allocate temporaries out of loop.
01317   TEMP_INTEGER(split_num);
01318   TEMP_INTEGER(split_den);
01319 
01320   // We start to compute num[] and den[].
01321   for (dimension_type i = external_space_dim; i-- > 0; ) {
01322     Coefficient& num_i = num[i];
01323     Coefficient& den_i = den[i];
01324     // Get the value of the variable from the tableau
01325     // (if it is not a basic variable, the value is 0).
01326     const dimension_type original_var = mapping[i+1].first;
01327     if (is_in_base(original_var, row)) {
01328       const Row& t_row = tableau[row];
01329       if (t_row[original_var] > 0) {
01330         neg_assign(num_i, t_row[0]);
01331         den_i = t_row[original_var];
01332       }
01333       else {
01334         num_i = t_row[0];
01335         neg_assign(den_i, t_row[original_var]);
01336       }
01337     }
01338     else {
01339       num_i = 0;
01340       den_i = 1;
01341     }
01342     // Check whether the variable was split.
01343     const dimension_type split_var = mapping[i+1].second;
01344     if (split_var != 0) {
01345       // The variable was split: get the value for the negative component,
01346       // having index mapping[i+1].second .
01347       // Like before, we he have to check if the variable is in base.
01348       if (is_in_base(split_var, row)) {
01349         const Row& t_row = tableau[row];
01350         if (t_row[split_var] > 0) {
01351           split_num = -t_row[0];
01352           split_den = t_row[split_var];
01353         }
01354         else {
01355           split_num = t_row[0];
01356           split_den = -t_row[split_var];
01357         }
01358         // We compute the lcm to compute subsequently the difference
01359         // between the 2 variables.
01360         lcm_assign(lcm, den_i, split_den);
01361         exact_div_assign(den_i, lcm, den_i);
01362         exact_div_assign(split_den, lcm, split_den);
01363         num_i *= den_i;
01364         sub_mul_assign(num_i, split_num, split_den);
01365         if (num_i == 0)
01366           den_i = 1;
01367         else
01368           den_i = lcm;
01369       }
01370       // Note: if the negative component was not in base, then
01371       // it has value zero and there is nothing left to do.
01372     }
01373   }
01374 
01375   // Compute the lcm of all denominators.
01376   lcm = den[0];
01377   for (dimension_type i = 1; i < external_space_dim; ++i)
01378     lcm_assign(lcm, lcm, den[i]);
01379   // Use the denominators to store the numerators' multipliers
01380   // and then compute the normalized numerators.
01381   for (dimension_type i = external_space_dim; i-- > 0; ) {
01382     exact_div_assign(den[i], lcm, den[i]);
01383     num[i] *= den[i];
01384   }
01385 
01386   // Finally, build the generator.
01387   Linear_Expression expr;
01388   for (dimension_type i = external_space_dim; i-- > 0; )
01389     expr += num[i] * Variable(i);
01390 
01391   MIP_Problem& x = const_cast<MIP_Problem&>(*this);
01392   x.last_generator = point(expr, lcm);
01393 }

void Parma_Polyhedra_Library::MIP_Problem::merge_split_variables ( dimension_type  var_index,
std::vector< dimension_type > &  nonfeasible_cs 
) [private]

Merges previously split variables in the tableau if a nonnegativity constraint is detected.

Parameters:
var_index The index of the variable that has to be merged.
nonfeasible_cs This will contain all the row indexes that are no more satisfied by the computed generator after merging a variable.

Definition at line 376 of file MIP_Problem.cc.

References base, is_in_base(), mapping, Parma_Polyhedra_Library::Matrix::num_columns(), Parma_Polyhedra_Library::Matrix::num_rows(), Parma_Polyhedra_Library::Matrix::permute_columns(), Parma_Polyhedra_Library::Matrix::remove_trailing_columns(), and tableau.

Referenced by parse_constraints().

00378                                                                  {
00379   const dimension_type tableau_nrows = tableau.num_rows();
00380   const dimension_type column = mapping[var_index].second;
00381 
00382   for (dimension_type i = 0; i < tableau_nrows; ++i) {
00383     // In the following case the negative side of the split variable is
00384     // in base: this means that the constraint will be nonfeasible.
00385     if (base[i] == mapping[var_index].second) {
00386       // CHECKME: I do not know if is possible that the positive and
00387       // the negative part of a split variable can be together in
00388       // base: it seems that this case is not possible. The algorithm
00389       // requires that condition.
00390 #ifndef NDEBUG
00391       for (dimension_type j = 0; j < tableau_nrows; ++j) {
00392         dimension_type dummy = 0;
00393         assert(!is_in_base(mapping[var_index].first, dummy));
00394       }
00395 #endif
00396       // We set base[i] to zero to keep track that that the constraint is not
00397       // feasible by `last_generator'.
00398       base[i] = 0;
00399       unfeasible_tableau_rows.push_back(i);
00400     }
00401   }
00402 
00403   const dimension_type tableau_cols = tableau.num_columns();
00404   // Remove the column.
00405   if (column != tableau_cols - 1) {
00406     std::vector<dimension_type> cycle;
00407     for (dimension_type j = tableau_cols - 1; j >= column; --j)
00408       cycle.push_back(j);
00409     cycle.push_back(0);
00410     tableau.permute_columns(cycle);
00411   }
00412   tableau.remove_trailing_columns(1);
00413 
00414   // var_index is no longer split.
00415   mapping[var_index].second = 0;
00416 
00417   // Adjust data structured, `shifting' the proper columns to the left by 1.
00418   const dimension_type base_size = base.size();
00419   for (dimension_type i = base_size; i-- > 0; )
00420     if (base[i] > column)
00421       --base[i];
00422   const dimension_type mapping_size = mapping.size();
00423   for (dimension_type i = mapping_size; i-- > 0; ) {
00424     if (mapping[i].first > column)
00425       --mapping[i].first;
00426     if (mapping[i].second > column)
00427       --mapping[i].second;
00428   }
00429 }

bool Parma_Polyhedra_Library::MIP_Problem::is_satisfied ( const Constraint c,
const Generator g 
) [static, private]

Returns true if and only if c is satisfied by g.

Definition at line 432 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Constraint::is_inequality(), Parma_Polyhedra_Library::Scalar_Products::sign(), Parma_Polyhedra_Library::Constraint::space_dimension(), and Parma_Polyhedra_Library::Generator::space_dimension().

Referenced by OK(), and parse_constraints().

00432                                                                     {
00433   // Scalar_Products::sign() requires the second argument to be at least
00434   // as large as the first one.
00435   int sp_sign = g.space_dimension() <= c.space_dimension()
00436     ? Scalar_Products::sign(g, c)
00437     : Scalar_Products::sign(c, g);
00438   return c.is_inequality() ? sp_sign >= 0 : sp_sign == 0;
00439 }

bool Parma_Polyhedra_Library::MIP_Problem::is_saturated ( const Constraint c,
const Generator g 
) [static, private]

Returns true if and only if c is saturated by g.

Definition at line 442 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Scalar_Products::sign(), Parma_Polyhedra_Library::Constraint::space_dimension(), and Parma_Polyhedra_Library::Generator::space_dimension().

Referenced by choose_branching_variable().

00442                                                                     {
00443   // Scalar_Products::sign() requires the second argument to be at least
00444   // as large as the first one.
00445   int sp_sign = g.space_dimension() <= c.space_dimension()
00446     ? Scalar_Products::sign(g, c)
00447     : Scalar_Products::sign(c, g);
00448   return sp_sign == 0;
00449 }

PPL::MIP_Problem_Status Parma_Polyhedra_Library::MIP_Problem::solve_mip ( bool &  have_incumbent_solution,
mpq_class &  incumbent_solution_value,
Generator incumbent_solution_point,
MIP_Problem mip,
const Variables_Set &  i_vars 
) [static, private]

Returns a status that encodes the solution of the MIP problem.

Parameters:
have_incumbent_solution It is used to store if the solving process has found a provisional optimum point.
incumbent_solution_value Encodes the evaluated value of the provisional optimum point found.
incumbent_solution_point If the method returns `OPTIMIZED', this will contain the optimality point.
mip The problem that has to be solved.
i_vars The variables that are constrained to take an integer value.

Definition at line 1526 of file MIP_Problem.cc.

References add_constraint(), Parma_Polyhedra_Library::assign_r(), Parma_Polyhedra_Library::Generator::coefficient(), Parma_Polyhedra_Library::Generator::divisor(), evaluate_objective_function(), Parma_Polyhedra_Library::gcd_assign(), Parma_Polyhedra_Library::is_canonical(), is_lp_satisfiable(), last_generator, Parma_Polyhedra_Library::MAXIMIZATION, Parma_Polyhedra_Library::MINIMIZATION, optimization_mode(), OPTIMIZED, Parma_Polyhedra_Library::OPTIMIZED_MIP_PROBLEM, second_phase(), space_dimension(), status, TEMP_INTEGER, Parma_Polyhedra_Library::UNBOUNDED_MIP_PROBLEM, and Parma_Polyhedra_Library::UNFEASIBLE_MIP_PROBLEM.

Referenced by solve().

01530                                                          {
01531   // Solve the problem as a non MIP one, it must be done internally.
01532   PPL::MIP_Problem_Status lp_status;
01533   if (lp.is_lp_satisfiable()) {
01534     lp.second_phase();
01535     lp_status = (lp.status == OPTIMIZED) ? OPTIMIZED_MIP_PROBLEM
01536       : UNBOUNDED_MIP_PROBLEM;
01537   }
01538   else
01539     return UNFEASIBLE_MIP_PROBLEM;
01540 
01541   DIRTY_TEMP0(mpq_class, tmp_rational);
01542 
01543   Generator p = point();
01544   TEMP_INTEGER(tmp_coeff1);
01545   TEMP_INTEGER(tmp_coeff2);
01546 
01547   if (lp_status == UNBOUNDED_MIP_PROBLEM)
01548     p = lp.last_generator;
01549   else {
01550     assert(lp_status == OPTIMIZED_MIP_PROBLEM);
01551     // Do not call optimizing_point().
01552     p = lp.last_generator;
01553     lp.evaluate_objective_function(p, tmp_coeff1, tmp_coeff2);
01554     assign_r(tmp_rational.get_num(), tmp_coeff1, ROUND_NOT_NEEDED);
01555     assign_r(tmp_rational.get_den(), tmp_coeff2, ROUND_NOT_NEEDED);
01556     assert(is_canonical(tmp_rational));
01557     if (have_incumbent_solution
01558         && ((lp.optimization_mode() == MAXIMIZATION
01559              && tmp_rational <= incumbent_solution_value)
01560             || (lp.optimization_mode() == MINIMIZATION
01561                 && tmp_rational >= incumbent_solution_value)))
01562       // Abandon this path.
01563       return lp_status;
01564   }
01565 
01566   bool found_satisfiable_generator = true;
01567   TEMP_INTEGER(gcd);
01568   const Coefficient& p_divisor = p.divisor();
01569   dimension_type nonint_dim;
01570   for (Variables_Set::const_iterator v_begin = i_vars.begin(),
01571          v_end = i_vars.end(); v_begin != v_end; ++v_begin) {
01572     gcd_assign(gcd, p.coefficient(Variable(*v_begin)), p_divisor);
01573     if (gcd != p_divisor) {
01574       nonint_dim = *v_begin;
01575       found_satisfiable_generator = false;
01576       break;
01577     }
01578   }
01579   if (found_satisfiable_generator) {
01580     // All the coordinates of `point' are satisfiable.
01581     if (lp_status == UNBOUNDED_MIP_PROBLEM) {
01582       // This is a point that belongs to the MIP_Problem.
01583       // In this way we are sure that we will return every time
01584       // a feasible point if requested by the user.
01585       incumbent_solution_point = p;
01586       return lp_status;
01587     }
01588     if (!have_incumbent_solution
01589         || (lp.optimization_mode() == MAXIMIZATION
01590             && tmp_rational > incumbent_solution_value)
01591         || tmp_rational < incumbent_solution_value) {
01592       incumbent_solution_value = tmp_rational;
01593       incumbent_solution_point = p;
01594       have_incumbent_solution = true;
01595 #if PPL_NOISY_SIMPLEX
01596       TEMP_INTEGER(num);
01597       TEMP_INTEGER(den);
01598       lp.evaluate_objective_function(p, num, den);
01599       std::cerr << "new value found: " << num << "/" << den << std::endl;
01600 #endif
01601     }
01602     return lp_status;
01603   }
01604 
01605   assert(nonint_dim < lp.space_dimension());
01606 
01607   assign_r(tmp_rational.get_num(), p.coefficient(Variable(nonint_dim)),
01608            ROUND_NOT_NEEDED);
01609   assign_r(tmp_rational.get_den(), p_divisor, ROUND_NOT_NEEDED);
01610   tmp_rational.canonicalize();
01611   assign_r(tmp_coeff1, tmp_rational, ROUND_DOWN);
01612   assign_r(tmp_coeff2, tmp_rational, ROUND_UP);
01613   {
01614     MIP_Problem lp_aux = lp;
01615     lp_aux.add_constraint(Variable(nonint_dim) <= tmp_coeff1);
01616     solve_mip(have_incumbent_solution, incumbent_solution_value,
01617               incumbent_solution_point, lp_aux, i_vars);
01618   }
01619   // TODO: change this when we will be able to remove constraints.
01620   lp.add_constraint(Variable(nonint_dim) >= tmp_coeff2);
01621   solve_mip(have_incumbent_solution, incumbent_solution_value,
01622             incumbent_solution_point, lp, i_vars);
01623   return have_incumbent_solution ? lp_status : UNFEASIBLE_MIP_PROBLEM;
01624 }

bool Parma_Polyhedra_Library::MIP_Problem::is_lp_satisfiable (  )  const [private]

Definition at line 1481 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Matrix::add_zero_columns(), external_space_dim, first_pending_constraint, initialized, input_cs, internal_space_dim, mapping, Parma_Polyhedra_Library::Matrix::num_columns(), OK(), OPTIMIZED, PARTIALLY_SATISFIABLE, process_pending_constraints(), SATISFIABLE, status, tableau, UNBOUNDED, and UNSATISFIABLE.

Referenced by is_mip_satisfiable(), is_satisfiable(), solve(), and solve_mip().

01481                                         {
01482 #if PPL_NOISY_SIMPLEX
01483   num_iterations = 0;
01484 #endif
01485   switch (status) {
01486   case UNSATISFIABLE:
01487     return false;
01488   case SATISFIABLE:
01489    // Intentionally fall through.
01490   case UNBOUNDED:
01491     // Intentionally fall through.
01492   case OPTIMIZED:
01493     // Intentionally fall through.
01494     return true;
01495   case PARTIALLY_SATISFIABLE:
01496     {
01497     MIP_Problem& x = const_cast<MIP_Problem&>(*this);
01498   // This code tries to handle the case that happens if the tableau is
01499   // empty, so it must be initialized.
01500   if (tableau.num_columns() == 0) {
01501     // Add two columns, the first that handles the inhomogeneous term and
01502     // the second that represent the `sign'.
01503     x.tableau.add_zero_columns(2);
01504     // Sync `mapping' for the inhomogeneous term.
01505     x.mapping.push_back(std::make_pair(0, 0));
01506     // The internal data structures are ready, so prepare for more
01507     // assertion to be checked.
01508     x.initialized = true;
01509   }
01510 
01511   // Apply incrementality to the pending constraint system.
01512   x.process_pending_constraints();
01513   // Update `first_pending_constraint': no more pending.
01514   x.first_pending_constraint = input_cs.size();
01515   // Update also `internal_space_dim'.
01516   x.internal_space_dim = x.external_space_dim;
01517   assert(OK());
01518   return (status != UNSATISFIABLE);
01519     }
01520   }
01521   // We should not be here!
01522   throw std::runtime_error("PPL internal error");
01523 }

bool Parma_Polyhedra_Library::MIP_Problem::is_mip_satisfiable ( MIP_Problem mip,
Generator p,
const Variables_Set &  i_vars 
) [static, private]

Used with MIP_Problems with a non empty `i_vars', returns true if and if only a MIP problem is satisfiable, returns false otherwise.

Parameters:
mip The problem that has to be solved.
p This will encode the feasible point, only if true is returned.
i_vars The variables that are constrained to take an integer value.

Definition at line 1681 of file MIP_Problem.cc.

References add_constraint(), Parma_Polyhedra_Library::assign_r(), choose_branching_variable(), Parma_Polyhedra_Library::Generator::coefficient(), Parma_Polyhedra_Library::Generator::divisor(), Parma_Polyhedra_Library::gcd_assign(), is_lp_satisfiable(), last_generator, space_dimension(), and TEMP_INTEGER.

Referenced by is_satisfiable().

01682                                                                   {
01683   // Solve the problem as a non MIP one, it must be done internally.
01684   if (!lp.is_lp_satisfiable())
01685     return false;
01686   DIRTY_TEMP0(mpq_class, tmp_rational);
01687 
01688   TEMP_INTEGER(tmp_coeff1);
01689   TEMP_INTEGER(tmp_coeff2);
01690   p = lp.last_generator;
01691 
01692   bool found_satisfiable_generator = true;
01693   dimension_type nonint_dim;
01694   const Coefficient& p_divisor = p.divisor();
01695 
01696 #if PPL_SIMPLEX_USE_MIP_HEURISTIC
01697   found_satisfiable_generator
01698     = choose_branching_variable(lp, i_vars, nonint_dim);
01699 #else
01700   TEMP_INTEGER(gcd);
01701   for (Variables_Set::const_iterator v_begin = i_vars.begin(),
01702          v_end = i_vars.end(); v_begin != v_end; ++v_begin) {
01703     gcd_assign(gcd, p.coefficient(Variable(*v_begin)), p_divisor);
01704     if (gcd != p_divisor) {
01705       nonint_dim = *v_begin;
01706       found_satisfiable_generator = false;
01707       break;
01708     }
01709   }
01710 #endif
01711 
01712   if (found_satisfiable_generator)
01713     return true;
01714 
01715 
01716   assert(nonint_dim < lp.space_dimension());
01717 
01718   assign_r(tmp_rational.get_num(), p.coefficient(Variable(nonint_dim)),
01719            ROUND_NOT_NEEDED);
01720   assign_r(tmp_rational.get_den(), p_divisor, ROUND_NOT_NEEDED);
01721   tmp_rational.canonicalize();
01722   assign_r(tmp_coeff1, tmp_rational, ROUND_DOWN);
01723   assign_r(tmp_coeff2, tmp_rational, ROUND_UP);
01724   {
01725     MIP_Problem lp_aux = lp;
01726     lp_aux.add_constraint(Variable(nonint_dim) <= tmp_coeff1);
01727     if (is_mip_satisfiable(lp_aux, p, i_vars))
01728       return true;
01729   }
01730   lp.add_constraint(Variable(nonint_dim) >= tmp_coeff2);
01731   return is_mip_satisfiable(lp, p, i_vars);
01732 }

bool Parma_Polyhedra_Library::MIP_Problem::choose_branching_variable ( const MIP_Problem mip,
const Variables_Set &  i_vars,
dimension_type branching_index 
) [static, private]

Returns true if and if only `last_generator' satisfies all the integrality coditions.

Parameters:
mip The MIP problem.
i_vars The variables that are constrained to take an integer value.
branching_index If false is returned, this will encode the variable index on which must be applied the `branch and bound' algorithm.

Definition at line 1627 of file MIP_Problem.cc.

References Parma_Polyhedra_Library::Generator::coefficient(), Parma_Polyhedra_Library::Generator::divisor(), Parma_Polyhedra_Library::gcd_assign(), input_cs, is_saturated(), last_generator, space_dimension(), and TEMP_INTEGER.

Referenced by is_mip_satisfiable().

01629                                                                              {
01630   // Insert here the variables that don't satisfy the integrality condition.
01631   const Constraint_Sequence& input_cs = mip.input_cs;
01632   const Generator& last_generator = mip.last_generator;
01633   const Coefficient& last_generator_divisor = last_generator.divisor();
01634   Variables_Set candidate_variables;
01635 
01636   TEMP_INTEGER(gcd);
01637   for (Variables_Set::const_iterator v_it = i_vars.begin(),
01638          v_end = i_vars.end(); v_it != v_end; ++v_it) {
01639     gcd_assign(gcd,
01640                last_generator.coefficient(Variable(*v_it)),
01641                last_generator_divisor);
01642     if (gcd != last_generator_divisor)
01643       candidate_variables.insert(*v_it);
01644   }
01645   // If this set is empty, we have finished.
01646   if (candidate_variables.empty())
01647     return true;
01648 
01649   // Check how many `active constraints' we have and track them.
01650   const dimension_type input_cs_num_rows = input_cs.size();
01651   std::deque<bool> satisfiable_constraints (input_cs_num_rows, false);
01652   for (dimension_type i = input_cs_num_rows; i-- > 0; )
01653     // An equality is an `active constraint' by definition.
01654     // If we have an inequality, check if it is an `active constraint'.
01655     if (input_cs[i].is_equality()
01656         || is_saturated(input_cs[i], last_generator))
01657       satisfiable_constraints[i] = true;
01658 
01659   dimension_type current_num_appearances = 0;
01660   dimension_type winning_num_appearances = 0;
01661 
01662   // For every candidate variable, check how many times this appear in the
01663   // active constraints.
01664   for (Variables_Set::const_iterator v_it = candidate_variables.begin(),
01665          v_end = candidate_variables.end(); v_it != v_end; ++v_it) {
01666     current_num_appearances = 0;
01667     for (dimension_type i = input_cs_num_rows; i-- > 0; )
01668       if (satisfiable_constraints[i]
01669           && *v_it < input_cs[i].space_dimension()
01670           && input_cs[i].coefficient(Variable(*v_it)) != 0)
01671         ++current_num_appearances;
01672     if (current_num_appearances >= winning_num_appearances) {
01673       winning_num_appearances = current_num_appearances;
01674       branching_index = *v_it;
01675     }
01676   }
01677   return false;
01678 }


Friends And Related Function Documentation

std::ostream & operator<< ( std::ostream &  s,
const MIP_Problem lp 
) [related]

Output operator.

Definition at line 2159 of file MIP_Problem.cc.

References constraints_begin(), constraints_end(), integer_space_dimensions(), Parma_Polyhedra_Library::MAXIMIZATION, objective_function(), and optimization_mode().

02159                                                                 {
02160   s << "Constraints:";
02161   for (MIP_Problem::const_iterator i = lp.constraints_begin(),
02162          i_end = lp.constraints_end(); i != i_end; ++i)
02163     s << "\n" << *i;
02164   s << "\nObjective function: "
02165     << lp.objective_function()
02166     << "\nOptimization mode: "
02167     << (lp.optimization_mode() == MAXIMIZATION
02168         ? "MAXIMIZATION"
02169         : "MINIMIZATION");
02170   s << "\nInteger variables: " << lp.integer_space_dimensions();
02171  return s;
02172 }

Specializes std::swap.

Definition at line 182 of file MIP_Problem.inlines.hh.

References swap().

00183                                             {
00184   x.swap(y);
00185 }


Member Data Documentation

The space dimension of the current (partial) solution of the MIP problem; it may be smaller than external_space_dim.

Definition at line 433 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), ascii_load(), is_lp_satisfiable(), process_pending_constraints(), and swap().

A map between the variables of `input_cs' and `tableau'.

Contains all the pairs (i, j) such that mapping[i].first encodes the index of the column in the tableau where input_cs[i] is stored; mapping[i].second not a zero, encodes the split part of the tableau of input_cs[i]. The "positive" one is represented by mapping[i].first and the "negative" one is represented by mapping[i].second.

Definition at line 449 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), ascii_load(), compute_generator(), external_memory_in_bytes(), is_lp_satisfiable(), merge_split_variables(), OK(), parse_constraints(), process_pending_constraints(), second_phase(), and swap().

The pricing method in use.

Definition at line 480 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), ascii_load(), get_control_parameter(), set_control_parameter(), and swap().

A Boolean encoding whether or not internal data structures have already been properly sized and populated: useful to allow for deeper checks in method OK().

Definition at line 487 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), ascii_load(), is_lp_satisfiable(), OK(), and swap().

The first index of `input_cs' containing a pending constraint.

Definition at line 493 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), ascii_load(), is_lp_satisfiable(), OK(), parse_constraints(), process_pending_constraints(), and swap().

The optimization mode requested.

Definition at line 499 of file MIP_Problem.defs.hh.

Referenced by ascii_dump(), optimization_mode(), process_pending_constraints(), second_phase(), set_optimization_mode(), and swap().

A set containing all the indexes of variables that are constrained to have an integer value.

Definition at line 508 of file MIP_Problem.defs.hh.

Referenced by add_to_integer_space_dimensions(), ascii_dump(), ascii_load(), integer_space_dimensions(), is_satisfiable(), OK(), solve(), and swap().


The documentation for this class was generated from the following files:

Generated on Sat Oct 11 10:41:06 2008 for PPL by  doxygen 1.5.6