#include <MIP_Problem.defs.hh>
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 ![]() 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 ![]() 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_Problem & | operator= (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_Expression & | objective_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 ![]() evaluating_point . | |
const Generator & | feasible_point () const |
Returns a feasible point for *this , if it exists. | |
const Generator & | optimizing_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 ![]() | |
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< Constraint > | Constraint_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_type > | base |
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 . |
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.
typedef std::vector<Constraint> Parma_Polyhedra_Library::MIP_Problem::Constraint_Sequence [private] |
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.
Names of MIP problems' control parameters.
Definition at line 403 of file MIP_Problem.defs.hh.
00403 { 00405 PRICING 00406 };
Possible values for MIP problem's control parameters.
Definition at line 409 of file MIP_Problem.defs.hh.
00409 { 00411 PRICING_STEEPEST_EDGE_FLOAT, 00413 PRICING_STEEPEST_EDGE_EXACT, 00415 PRICING_TEXTBOOK 00416 };
enum Parma_Polyhedra_Library::MIP_Problem::Status [private] |
An enumerated type describing the internal status of the MIP problem.
Definition at line 455 of file MIP_Problem.defs.hh.
00455 { 00457 UNSATISFIABLE, 00459 SATISFIABLE, 00461 UNBOUNDED, 00463 OPTIMIZED, 00469 PARTIALLY_SATISFIABLE 00470 };
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 on a vector space under no constraints at all: the origin of the vector space is an optimal solution.
dim | The dimension of the vector space enclosing *this (optional argument with default value ![]() |
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 }
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 , the objective function
obj
and optimization mode mode
; those dimensions whose indices occur in int_vars
are constrained to take an integer value.
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 ![]() | |
mode | The optimization mode (optional argument with default value MAXIMIZATION ). |
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 . |
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 , the objective function
obj
and optimization mode mode
.
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 ![]() | |
mode | The optimization mode (optional argument with default value MAXIMIZATION ). |
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
.
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 ![]() | |
mode | The optimization mode (optional argument with default value MAXIMIZATION ). |
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] |
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] |
Returns the space dimension of the MIP problem.
Definition at line 37 of file MIP_Problem.inlines.hh.
References external_space_dim.
Referenced by add_constraint(), add_constraints(), add_space_dimensions_and_embed(), choose_branching_variable(), evaluate_objective_function(), is_mip_satisfiable(), process_pending_constraints(), set_objective_function(), and solve_mip().
00037 { 00038 return external_space_dim; 00039 }
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 .
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.
m | The number of dimensions to add. |
std::length_error | Thrown if adding m new space dimensions would cause the vector space to exceed dimension max_space_dimension() . |
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.
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.
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.
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
.
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] |
Sets the optimization mode to mode
.
Definition at line 68 of file MIP_Problem.inlines.hh.
References OK(), opt_mode, OPTIMIZED, SATISFIABLE, status, and UNBOUNDED.
Referenced by ascii_load(), 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().
00068 { 00069 if (opt_mode != mode) { 00070 opt_mode = mode; 00071 if (status == UNBOUNDED || status == OPTIMIZED) 00072 status = SATISFIABLE; 00073 assert(OK()); 00074 } 00075 }
bool Parma_Polyhedra_Library::MIP_Problem::is_satisfiable | ( | ) | const |
Checks satisfiability of *this
.
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.
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 is the result of evaluating the objective function on
evaluating_point
.
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. |
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.
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.
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 is the solution of the optimization problem.
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 |
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().
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
.
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.
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.
UNSATISFIABLE
if is detected a trivially false constraint, SATISFIABLE
otherwise.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.
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.
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 ![]() |
x
and y
having the element of index k
equal 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.
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.
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.
0
is returned.
Recall that, due to the exact integer implementation of the algorithm, our tableau doesn't contain the ``real'' 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.
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.
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.
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.
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.
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.
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.
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 }
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 }
void swap | ( | Parma_Polyhedra_Library::MIP_Problem & | x, | |
Parma_Polyhedra_Library::MIP_Problem & | y | |||
) | [related] |
Specializes std::swap
.
Definition at line 182 of file MIP_Problem.inlines.hh.
References swap().
00183 { 00184 x.swap(y); 00185 }
The dimension of the vector space.
Definition at line 427 of file MIP_Problem.defs.hh.
Referenced by add_space_dimensions_and_embed(), add_to_integer_space_dimensions(), ascii_dump(), ascii_load(), compute_generator(), is_lp_satisfiable(), OK(), parse_constraints(), process_pending_constraints(), space_dimension(), and swap().
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().
The matrix encoding the current feasible region in tableau form.
Definition at line 436 of file MIP_Problem.defs.hh.
Referenced by ascii_dump(), ascii_load(), compute_generator(), compute_simplex_using_exact_pricing(), compute_simplex_using_steepest_edge_float(), erase_artificials(), external_memory_in_bytes(), get_exiting_base_index(), is_lp_satisfiable(), merge_split_variables(), OK(), pivot(), process_pending_constraints(), second_phase(), steepest_edge_exact_entering_index(), steepest_edge_float_entering_index(), and swap().
The working cost function.
Definition at line 439 of file MIP_Problem.defs.hh.
Referenced by ascii_dump(), ascii_load(), compute_simplex_using_exact_pricing(), compute_simplex_using_steepest_edge_float(), erase_artificials(), external_memory_in_bytes(), OK(), pivot(), process_pending_constraints(), second_phase(), steepest_edge_exact_entering_index(), steepest_edge_float_entering_index(), swap(), and textbook_entering_index().
std::vector<std::pair<dimension_type, dimension_type> > Parma_Polyhedra_Library::MIP_Problem::mapping [private] |
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().
std::vector<dimension_type> Parma_Polyhedra_Library::MIP_Problem::base [private] |
The current basic solution.
Definition at line 452 of file MIP_Problem.defs.hh.
Referenced by ascii_dump(), ascii_load(), erase_artificials(), external_memory_in_bytes(), get_exiting_base_index(), is_in_base(), merge_split_variables(), OK(), pivot(), process_pending_constraints(), second_phase(), steepest_edge_exact_entering_index(), steepest_edge_float_entering_index(), and swap().
The internal state of the MIP problem.
Definition at line 473 of file MIP_Problem.defs.hh.
Referenced by add_constraint(), add_constraints(), add_space_dimensions_and_embed(), add_to_integer_space_dimensions(), ascii_dump(), ascii_load(), is_lp_satisfiable(), is_satisfiable(), OK(), process_pending_constraints(), second_phase(), set_objective_function(), set_optimization_mode(), solve(), solve_mip(), 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().
bool Parma_Polyhedra_Library::MIP_Problem::initialized [private] |
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 sequence of constraints describing the feasible region.
Definition at line 490 of file MIP_Problem.defs.hh.
Referenced by add_constraint(), add_constraints(), ascii_dump(), ascii_load(), choose_branching_variable(), constraints_begin(), constraints_end(), external_memory_in_bytes(), is_lp_satisfiable(), MIP_Problem(), OK(), parse_constraints(), process_pending_constraints(), 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 objective function to be optimized.
Definition at line 496 of file MIP_Problem.defs.hh.
Referenced by ascii_dump(), ascii_load(), evaluate_objective_function(), external_memory_in_bytes(), objective_function(), OK(), process_pending_constraints(), second_phase(), set_objective_function(), 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().
The last successfully computed feasible or optimizing point.
Definition at line 502 of file MIP_Problem.defs.hh.
Referenced by ascii_dump(), ascii_load(), choose_branching_variable(), compute_generator(), external_memory_in_bytes(), feasible_point(), is_mip_satisfiable(), is_satisfiable(), OK(), optimizing_point(), parse_constraints(), process_pending_constraints(), solve(), solve_mip(), and swap().
Variables_Set Parma_Polyhedra_Library::MIP_Problem::i_variables [private] |
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().