Line data Source code
1 : // Author(s): Jeroen Keiren and Jan Friso Groote
2 : // Copyright: see the accompanying file COPYING or copy at
3 : // https://github.com/mCRL2org/mCRL2/blob/master/COPYING
4 : //
5 : // Distributed under the Boost Software License, Version 1.0.
6 : // (See accompanying file LICENSE_1_0.txt or copy at
7 : // http://www.boost.org/LICENSE_1_0.txt)
8 : //
9 : /// \file linear_inequalities.h
10 : /// \brief Contains a class linear_inequality to represent mcrl2 data
11 : /// expressions that are linear equalities, or inequalities.
12 : /// Furthermore, it contains some operations on these linear
13 : /// inequalities, such as Fourier-Motzkin elimination.
14 :
15 :
16 : #ifndef MCRL2_DATA_LINEAR_INEQUALITY_H
17 : #define MCRL2_DATA_LINEAR_INEQUALITY_H
18 :
19 : #include "mcrl2/data/rewriter.h"
20 : #include "mcrl2/data/substitutions/map_substitution.h"
21 :
22 : namespace mcrl2
23 : {
24 :
25 : namespace data
26 : {
27 :
28 : // Functions below should be made available in the data library.
29 : data_expression& real_zero();
30 : data_expression& real_one();
31 : data_expression& real_minus_one();
32 : data_expression min(const data_expression& e1,const data_expression& e2,const rewriter&);
33 : data_expression max(const data_expression& e1,const data_expression& e2,const rewriter&);
34 : bool is_closed_real_number(const data_expression& e);
35 : bool is_negative(const data_expression& e,const rewriter& r);
36 : bool is_positive(const data_expression& e,const rewriter& r);
37 : bool is_zero(const data_expression& e);
38 :
39 :
40 : // Efficient construction of times on reals.
41 : // The functions times, divide, plus, minus, negate and abs in the data library are not efficient since they determine the sort at runtime.
42 100 : inline application real_times(const data_expression& arg0, const data_expression& arg1)
43 : {
44 100 : static function_symbol times_f(sort_real::times(sort_real::real_(), sort_real::real_()));
45 200 : assert(arg0.sort()==sort_real::real_() && arg1.sort()==sort_real::real_());
46 100 : return application(times_f,arg0,arg1);
47 : }
48 :
49 101 : inline application real_plus(const data_expression& arg0, const data_expression& arg1)
50 : {
51 101 : static function_symbol plus_f(sort_real::plus(sort_real::real_(), sort_real::real_()));
52 202 : assert(arg0.sort()==sort_real::real_() && arg1.sort()==sort_real::real_());
53 101 : return application(plus_f,arg0,arg1);
54 : }
55 :
56 50 : inline application real_divides(const data_expression& arg0, const data_expression& arg1)
57 : {
58 50 : static function_symbol divides_f(sort_real::divides(sort_real::real_(), sort_real::real_()));
59 100 : assert(arg0.sort()==sort_real::real_() && arg1.sort()==sort_real::real_());
60 50 : return application(divides_f,arg0,arg1);
61 : }
62 :
63 19 : inline application real_minus(const data_expression& arg0, const data_expression& arg1)
64 : {
65 19 : static function_symbol minus_f(sort_real::minus(sort_real::real_(), sort_real::real_()));
66 38 : assert(arg0.sort()==sort_real::real_() && arg1.sort()==sort_real::real_());
67 19 : return application(minus_f,arg0,arg1);
68 : }
69 :
70 46 : inline application real_negate(const data_expression& arg)
71 : {
72 46 : static function_symbol negate_f(sort_real::negate(sort_real::real_()));
73 46 : assert(arg.sort()==sort_real::real_());
74 46 : return application(negate_f,arg);
75 : }
76 :
77 1 : inline application real_abs(const data_expression& arg)
78 : {
79 1 : static function_symbol abs_f(sort_real::abs(sort_real::real_()));
80 1 : assert(arg.sort()==sort_real::real_());
81 1 : return application(abs_f,arg);
82 : }
83 :
84 : // End of functions that ought to be defined elsewhere.
85 :
86 :
87 : inline data_expression rewrite_with_memory(const data_expression& t,const rewriter& r);
88 :
89 : // prototype
90 : class linear_inequality;
91 : namespace detail
92 : {
93 : class lhs_t;
94 : }
95 :
96 : inline std::string pp(const linear_inequality& l);
97 : template <class TYPE>
98 : inline std::string pp_vector(const TYPE& inequalities);
99 :
100 : namespace detail
101 : {
102 : enum comparison_t { less, less_eq, equal };
103 :
104 : inline std::string pp(const detail::lhs_t& lhs);
105 :
106 : inline detail::comparison_t negate(const detail::comparison_t t)
107 : {
108 : switch (t)
109 : {
110 : case detail::less: return detail::less_eq;
111 : case detail::less_eq: return detail::less;
112 : case detail::equal: return detail::equal;
113 : };
114 : return detail::equal; // This return statement should be unreachable. It is added to suppress a compiler warning.
115 : }
116 :
117 0 : inline std::string pp(const detail::comparison_t t)
118 : {
119 0 : switch (t)
120 : {
121 0 : case detail::less: return "<";
122 0 : case detail::less_eq: return "<=";
123 0 : case detail::equal: return "==";
124 : };
125 0 : return "##"; // This return statement should be unreachable. It is added to suppress a compiler warning.
126 : }
127 :
128 817 : inline atermpp::function_symbol f_variable_with_a_rational_factor()
129 : {
130 817 : static atermpp::function_symbol f("variable_with_a_rational_factor",2);
131 817 : return f;
132 : }
133 :
134 : class variable_with_a_rational_factor: public atermpp::aterm_appl
135 : {
136 : public:
137 : // \brief default constructor
138 94 : variable_with_a_rational_factor(const variable& v, const data_expression& f)
139 94 : : atermpp::aterm_appl(f_variable_with_a_rational_factor(),v,f)
140 : {
141 94 : assert(f!=real_zero());
142 94 : }
143 :
144 : // \brief Get the variable in a variable/rational factor pair.
145 463 : const variable& variable_name() const
146 : {
147 463 : assert(is_variable_with_a_rational_factor());
148 463 : return atermpp::down_cast<variable>((*this)[0]);
149 : }
150 :
151 : // \brief Get the rational factor in a variable/rational factor pair.
152 260 : const data_expression& factor() const
153 : {
154 260 : assert(is_variable_with_a_rational_factor());
155 260 : return atermpp::down_cast<data_expression>((*this)[1]);
156 : }
157 :
158 : // \brief Check that a term is indeed a variable with a rational factor pair.
159 723 : bool is_variable_with_a_rational_factor() const
160 : {
161 723 : return function()==f_variable_with_a_rational_factor();
162 : }
163 :
164 : // \brief Transform the variable with factor to a data_expression.
165 3 : data_expression transform_to_data_expression() const
166 : {
167 3 : return real_times(factor(),variable_name());
168 : }
169 : };
170 :
171 : // typedef atermpp::term_list<variable_with_a_rational_factor> lhs_t;
172 : typedef std::map < variable, data_expression > map_based_lhs_t;
173 :
174 : class lhs_t: public atermpp::term_list<variable_with_a_rational_factor>
175 : {
176 : public:
177 : /// \brief Constructor.
178 74 : lhs_t()
179 74 : : atermpp::term_list<variable_with_a_rational_factor>()
180 74 : {}
181 :
182 : /// \brief Constructor from an aterm.
183 639 : explicit lhs_t(const aterm& t)
184 639 : : atermpp::term_list<variable_with_a_rational_factor>(t)
185 639 : {}
186 :
187 : /// \brief Constructor
188 : template <class ITERATOR>
189 38 : lhs_t(const ITERATOR begin, const ITERATOR end)
190 38 : : atermpp::term_list<variable_with_a_rational_factor>(begin, end)
191 38 : {}
192 :
193 : /// \brief Constructor
194 : template <class ITERATOR, class TRANSFORMER>
195 1 : lhs_t(const ITERATOR begin, const ITERATOR end, TRANSFORMER f)
196 1 : : atermpp::term_list<variable_with_a_rational_factor>(begin, end, f)
197 1 : {}
198 :
199 : /// \brief Give an iterator of the factor/variable pair for v, or end() if v does not occur.
200 40 : lhs_t::const_iterator find(const variable& v) const
201 : {
202 40 : return std::find_if(begin(),
203 40 : end(),
204 70 : [&](const detail::variable_with_a_rational_factor& p)
205 190 : {return p.variable_name()==v;});
206 : }
207 :
208 : /// \brief Erase a variable and its factor.
209 8 : lhs_t erase(const variable& v) const
210 : {
211 8 : std::vector <variable_with_a_rational_factor> result;
212 25 : for(const variable_with_a_rational_factor& p: *this)
213 : {
214 17 : if (p.variable_name()!=v)
215 : {
216 9 : result.push_back(p);
217 : }
218 : }
219 16 : return lhs_t(result.begin(),result.end());
220 8 : }
221 :
222 : /// \brief Give the factor of variable v.
223 14 : const data_expression& operator[](const variable& v) const
224 : {
225 14 : lhs_t::const_iterator i=find(v);
226 14 : if (i==end()) // Not found.
227 : {
228 0 : return real_zero();
229 : }
230 14 : return i->factor();
231 : }
232 :
233 : /// \brief Give the factor of variable v.
234 12 : std::size_t count(const variable& v) const
235 : {
236 12 : lhs_t::const_iterator i=find(v);
237 12 : if (i==end()) // Not found.
238 : {
239 4 : return 0;
240 : }
241 8 : return 1;
242 : }
243 :
244 : /// \brief Evaluate the variables in this lhs_t according to the subsitution function.
245 : template <typename SubstitutionFunction>
246 34 : data_expression evaluate(const SubstitutionFunction& beta, const rewriter& r) const
247 : {
248 34 : data_expression result=real_zero();
249 34 : bool result_defined=false; // save adding the initial zero.
250 136 : for (const variable_with_a_rational_factor& p: *this)
251 : {
252 68 : if (result_defined)
253 : {
254 34 : result=rewrite_with_memory(real_plus(result,real_times(beta(p.variable_name()),p.factor())), r);
255 : }
256 : else
257 : {
258 34 : result=rewrite_with_memory(real_times(beta(p.variable_name()),p.factor()), r);
259 34 : result_defined=true;
260 : }
261 :
262 : }
263 34 : return result;
264 0 : }
265 :
266 : // \brief Transform this lhs to a data_expression.
267 3 : data_expression transform_to_data_expression() const
268 : {
269 3 : if (size()==0)
270 : {
271 0 : return real_zero();
272 : }
273 3 : data_expression result=real_zero();
274 :
275 6 : for(const variable_with_a_rational_factor& p: *this)
276 : {
277 3 : if (result==real_zero())
278 : {
279 3 : result=p.transform_to_data_expression();
280 : }
281 : else
282 : {
283 0 : result=real_plus(result,p.transform_to_data_expression());
284 : }
285 : }
286 3 : return result;
287 3 : }
288 : };
289 :
290 51 : inline void set_factor_for_a_variable(detail::map_based_lhs_t& new_lhs, const variable& x, const data_expression& e)
291 : {
292 51 : assert(is_closed_real_number(e));
293 51 : if (e==real_zero())
294 : {
295 3 : detail::map_based_lhs_t::iterator i=new_lhs.find(x);
296 3 : if (i!=new_lhs.end())
297 : {
298 3 : new_lhs.erase(i);
299 : }
300 : }
301 : else
302 : {
303 48 : new_lhs[x]=e;
304 : }
305 51 : }
306 :
307 4 : inline lhs_t set_factor_for_a_variable(const lhs_t& lhs, const variable& x, const data_expression& e)
308 : {
309 4 : assert(is_closed_real_number(e));
310 4 : bool inserted=false;
311 4 : std::vector<variable_with_a_rational_factor> result;
312 9 : for(const variable_with_a_rational_factor& p: lhs)
313 : {
314 5 : if (!inserted && x<=p.variable_name())
315 : {
316 1 : result.emplace_back(x,e);
317 1 : inserted=true;
318 1 : if (x!=p.variable_name())
319 : {
320 1 : result.emplace_back(p.variable_name(),p.factor());
321 : }
322 : }
323 4 : else result.emplace_back(p.variable_name(),p.factor());
324 : }
325 4 : if (!inserted)
326 : {
327 3 : result.emplace_back(x,e);
328 : }
329 4 : assert(std::find(result.begin(),result.end(),variable_with_a_rational_factor(x,e))!=result.end());
330 8 : return lhs_t(result.begin(),result.end());
331 4 : }
332 :
333 33 : inline const lhs_t map_to_lhs_type(const map_based_lhs_t& lhs)
334 : {
335 33 : lhs_t result;
336 77 : for(map_based_lhs_t::const_reverse_iterator i=lhs.rbegin(); i!=lhs.rend(); ++i)
337 : {
338 44 : result.push_front(variable_with_a_rational_factor(i->first,i->second));
339 : }
340 33 : return result;
341 0 : }
342 :
343 0 : inline const lhs_t map_to_lhs_type(const map_based_lhs_t& lhs, const data_expression& factor, const rewriter& r)
344 : {
345 0 : assert(factor!=real_one() && factor!=real_minus_one());
346 0 : lhs_t result;
347 0 : for(map_based_lhs_t::const_reverse_iterator i=lhs.rbegin(); i!=lhs.rend(); ++i)
348 : {
349 0 : result.push_front(variable_with_a_rational_factor(i->first, rewrite_with_memory(real_divides(i->second,factor), r)));
350 : }
351 0 : return result;
352 0 : }
353 :
354 68 : inline bool is_well_formed(const lhs_t& lhs)
355 : {
356 68 : if (lhs.empty())
357 : {
358 30 : return true;
359 : }
360 38 : if (lhs.front().factor()!=real_one() && lhs.front().factor()!=real_minus_one())
361 : {
362 0 : return false;
363 : }
364 38 : variable last_variable_seen=lhs.front().variable_name();
365 38 : bool first=true;
366 88 : for(const variable_with_a_rational_factor& p: lhs)
367 : {
368 50 : if (!first)
369 : {
370 0 : first=false;
371 0 : if (p.variable_name()<=last_variable_seen)
372 : {
373 0 : return false;
374 : }
375 0 : last_variable_seen=p.variable_name();
376 : }
377 : }
378 38 : return true;
379 38 : }
380 :
381 8 : inline const lhs_t remove_variable_and_divide(const lhs_t& lhs, const variable& v, const data_expression& f, const rewriter& r)
382 : {
383 8 : std::vector <variable_with_a_rational_factor> result;
384 20 : for(const variable_with_a_rational_factor& p: lhs)
385 : {
386 12 : if (p.variable_name()!=v)
387 : {
388 4 : result.emplace_back(p.variable_name(),rewrite_with_memory(real_divides(p.factor(),f), r));
389 : }
390 : }
391 16 : return lhs_t(result.begin(),result.end());
392 8 : }
393 :
394 36 : inline void emplace_back_if_not_zero(std::vector<detail::variable_with_a_rational_factor>& r, const variable& v, const data_expression& f)
395 : {
396 36 : if (f!=real_zero())
397 : {
398 32 : r.emplace_back(v,f);
399 : }
400 36 : }
401 :
402 :
403 : template < application Operation(const data_expression&, const data_expression&) >
404 9 : inline const lhs_t meta_operation_constant(const lhs_t& argument, const data_expression& v, const rewriter& r)
405 : {
406 9 : std::vector<detail::variable_with_a_rational_factor> result;
407 9 : result.reserve(argument.size());
408 38 : for (const detail::variable_with_a_rational_factor& p: argument)
409 : {
410 20 : emplace_back_if_not_zero(result,p.variable_name(), rewrite_with_memory(Operation(v,p.factor()),r));
411 : }
412 18 : return lhs_t(result.begin(),result.end());
413 9 : }
414 :
415 : // Template method to add or subtract lhs_t's
416 : // template < application Operation(const data_expression&, const data_expression&) >
417 : template <class OPERATION>
418 9 : inline lhs_t meta_operation_lhs(const lhs_t& argument1,
419 : const lhs_t& argument2,
420 : OPERATION operation,
421 : const rewriter& r)
422 : {
423 9 : std::vector<detail::variable_with_a_rational_factor> result;
424 9 : result.reserve(argument1.size()+argument2.size());
425 :
426 9 : lhs_t::const_iterator i1=argument1.begin();
427 9 : lhs_t::const_iterator i2=argument2.begin();
428 :
429 17 : while (i1!=argument1.end() && i2!=argument2.end())
430 : {
431 8 : if (i1->variable_name()<i2->variable_name())
432 : {
433 0 : emplace_back_if_not_zero(result,i1->variable_name(), rewrite_with_memory(operation(i1->factor(),real_zero()),r));
434 0 : ++i1;
435 : }
436 8 : else if (i1->variable_name()>i2->variable_name())
437 : {
438 4 : emplace_back_if_not_zero(result,i2->variable_name(), rewrite_with_memory(operation(real_zero(),i2->factor()),r));
439 4 : ++i2;
440 : }
441 : else
442 : {
443 4 : assert(i1->variable_name()==i2->variable_name());
444 4 : emplace_back_if_not_zero(result,i1->variable_name(), rewrite_with_memory(operation(i1->factor(),i2->factor()),r));
445 4 : ++i1;
446 4 : ++i2;
447 : }
448 : }
449 :
450 9 : if (i1==argument1.end())
451 : {
452 15 : while (i2!=argument2.end())
453 : {
454 7 : emplace_back_if_not_zero(result,i2->variable_name(), rewrite_with_memory(operation(real_zero(),i2->factor()),r));
455 7 : ++i2;
456 : }
457 : }
458 : else
459 : {
460 2 : while (i1!=argument1.end())
461 : {
462 1 : emplace_back_if_not_zero(result,i1->variable_name(), rewrite_with_memory(operation(i1->factor(),real_zero()),r));
463 1 : ++i1;
464 : }
465 : }
466 18 : return lhs_t(result.begin(),result.end());
467 9 : }
468 :
469 : inline const lhs_t add(const data_expression& v, const lhs_t& argument, const rewriter& r)
470 : {
471 : return meta_operation_constant<real_plus>(argument,v,r);
472 : }
473 :
474 : inline const lhs_t subtract(const lhs_t& argument, const data_expression& v, const rewriter& r)
475 : {
476 : return meta_operation_constant<real_minus>(argument, v, r);
477 : }
478 :
479 8 : inline const lhs_t multiply(const lhs_t& argument, const data_expression& v, const rewriter& r)
480 : {
481 8 : return meta_operation_constant<real_times>(argument, v, r);
482 : }
483 :
484 1 : inline const lhs_t divide(const lhs_t& argument, const data_expression& v, const rewriter& r)
485 : {
486 1 : return meta_operation_constant<real_divides>(argument, v, r);
487 : }
488 :
489 4 : inline const lhs_t add(const lhs_t& argument, const lhs_t& e, const rewriter& r)
490 : {
491 : return meta_operation_lhs(argument,
492 : e,
493 10 : [](const data_expression& d1, const data_expression& d2)->data_expression
494 10 : { return real_plus(d1,d2); },
495 4 : r);
496 : }
497 :
498 4 : inline const lhs_t subtract(const lhs_t& argument, const lhs_t& e, const rewriter& r)
499 : {
500 : return meta_operation_lhs(argument,
501 : e,
502 4 : [](const data_expression& d1, const data_expression& d2)->data_expression
503 4 : { return real_minus(d1,d2); },
504 4 : r);
505 : }
506 :
507 1 : inline std::string pp(const lhs_t& lhs)
508 : {
509 1 : std::string s;
510 1 : if (lhs.begin()==lhs.end())
511 : {
512 1 : s="0";
513 : }
514 1 : for (lhs_t::const_iterator i=lhs.begin(); i!=lhs.end(); ++i)
515 : {
516 0 : s=s + (i==lhs.begin()?"":" + ") ;
517 :
518 0 : if (i->factor()==real_one())
519 : {
520 0 : s=s + pp(i->variable_name());
521 : }
522 0 : else if (i->factor()==real_minus_one())
523 : {
524 0 : s=s + "-" + pp(i->variable_name());
525 : }
526 : else
527 : {
528 0 : s=s + data::pp(i->factor()) + "*" + data::pp(i->variable_name());
529 : }
530 : }
531 1 : return s;
532 0 : }
533 :
534 160 : inline atermpp::function_symbol linear_inequality_less()
535 : {
536 160 : static atermpp::function_symbol f("linear_inequality_less",2);
537 160 : return f;
538 : }
539 :
540 :
541 90 : inline atermpp::function_symbol linear_inequality_less_equal()
542 : {
543 90 : static atermpp::function_symbol f("linear_inequality_less_equal",2);
544 90 : return f;
545 : }
546 :
547 :
548 2 : inline atermpp::function_symbol linear_inequality_equal()
549 : {
550 2 : static atermpp::function_symbol f("linear_inequality_equal",2);
551 2 : return f;
552 : }
553 :
554 : } // end namespace detail
555 :
556 : class linear_inequality: public atermpp::aterm_appl
557 : {
558 : // The structure of a linear equality is a function application
559 : // to two arguments. The function application is either linear_inequality_less,
560 : // linear_inequality_less_equal, or linear_inequality_equal. There are two arguments,
561 : // namely an ordered list of pairs of a variable and a factor. This list is ordered
562 : // on the variables, and a closed expression which forms the right hand side.
563 : // The first variable in the list has a factor one or minus one.
564 :
565 : protected:
566 :
567 116 : static void parse_and_store_expression(
568 : const data_expression& e,
569 : detail::map_based_lhs_t& new_lhs,
570 : data_expression& new_rhs,
571 : const rewriter& r,
572 : const bool negate = false,
573 : const data_expression& factor = real_one())
574 : {
575 116 : if (sort_real::is_minus_application(e))
576 : {
577 5 : parse_and_store_expression(binary_left1(e),new_lhs,new_rhs,r,negate,factor);
578 5 : parse_and_store_expression(binary_right1(e),new_lhs,new_rhs,r,!negate,factor);
579 : }
580 111 : else if (sort_real::is_negate_application(e))
581 : {
582 13 : parse_and_store_expression(unary_operand1(e),new_lhs,new_rhs,r,!negate,factor);
583 : }
584 98 : else if (sort_real::is_plus_application(e))
585 : {
586 9 : parse_and_store_expression(binary_left1(e),new_lhs,new_rhs,r,negate,factor);
587 9 : parse_and_store_expression(binary_right1(e),new_lhs,new_rhs,r,negate,factor);
588 : }
589 89 : else if (sort_real::is_times_application(e))
590 : {
591 1 : data_expression lhs = rewrite_with_memory(binary_left1(e),r);
592 1 : data_expression rhs = rewrite_with_memory(binary_right1(e),r);
593 1 : if (is_closed_real_number(lhs))
594 : {
595 0 : parse_and_store_expression(rhs,new_lhs,new_rhs,r,negate,real_times(lhs,factor));
596 : }
597 1 : else if (is_closed_real_number(rhs))
598 : {
599 0 : parse_and_store_expression(lhs,new_lhs,new_rhs,r,negate,real_times(rhs,factor));
600 : }
601 : else
602 : {
603 1 : throw mcrl2::runtime_error("Expect constant multiplies expression: " + pp(e) + "\n");
604 : }
605 2 : }
606 88 : else if (is_variable(e))
607 : {
608 51 : const variable& v = atermpp::down_cast<variable>(e);
609 51 : if(v.sort() != sort_real::real_())
610 : {
611 0 : throw mcrl2::runtime_error("Encountered a variable in a real expression which is not of sort real: " + pp(e) + "\n");
612 : }
613 51 : const detail::map_based_lhs_t::const_iterator var_factor = new_lhs.find(v);
614 :
615 51 : const data_expression neg_factor(negate ? real_negate(factor) : factor);
616 51 : const data_expression new_factor(var_factor == new_lhs.end() ? neg_factor : real_plus(var_factor->second, neg_factor));
617 51 : detail::set_factor_for_a_variable(new_lhs, v, rewrite_with_memory(new_factor, r));
618 51 : }
619 37 : else if (is_closed_real_number(rewrite_with_memory(e,r)))
620 : {
621 : // We are reasoning about the rhs, so apply double negation in case 'negate' is true
622 37 : data_expression add_to_rhs(negate ? e : real_negate(e));
623 37 : if(factor != real_one())
624 : {
625 0 : add_to_rhs = real_times(factor, add_to_rhs);
626 : }
627 37 : new_rhs = rewrite_with_memory(real_plus(new_rhs, add_to_rhs), r);
628 37 : }
629 : else
630 : {
631 0 : throw mcrl2::runtime_error("Expect linear expression over reals: " + pp(e) + "\n");
632 : }
633 115 : }
634 :
635 :
636 : public:
637 :
638 : /// \brief Constructor yielding an inconsistent inequality.
639 25 : linear_inequality()
640 25 : : linear_inequality(detail::lhs_t(),real_zero(),detail::less)
641 25 : {}
642 :
643 : /// Basic constructor.
644 68 : linear_inequality(const detail::lhs_t& lhs, const data_expression& r, detail::comparison_t t)
645 136 : : atermpp::aterm_appl((t==detail::less?
646 : detail::linear_inequality_less():
647 : (t==detail::less_eq?detail::linear_inequality_less_equal():detail::linear_inequality_equal())),
648 136 : lhs,r)
649 : {
650 68 : assert(detail::is_well_formed(lhs));
651 68 : assert(t==detail::less || t==detail::less_eq || t==detail::equal);
652 68 : }
653 :
654 : /// \brief constructor.
655 5 : linear_inequality(const detail::lhs_t& lhs, const data_expression& rhs, detail::comparison_t comparison, const rewriter& r)
656 5 : {
657 5 : if (lhs.empty())
658 : {
659 1 : *this=linear_inequality(detail::lhs_t(),rhs,comparison);
660 4 : return;
661 : }
662 : // Normalize the linear_inequality such that the first term has factor 1 or -1.
663 4 : data_expression factor=lhs.begin()->factor();
664 4 : if (factor==real_one() || factor==real_minus_one())
665 : {
666 3 : *this=linear_inequality(lhs,rhs,comparison);
667 3 : return;
668 : }
669 1 : factor=rewrite_with_memory(real_abs(factor), r);
670 :
671 1 : *this=linear_inequality(divide(lhs,factor,r),rewrite_with_memory(real_divides(rhs,factor), r),comparison);
672 4 : }
673 :
674 : /// \brief constructor.
675 38 : linear_inequality(const data_expression& lhs,
676 : const data_expression& rhs,
677 : const detail::comparison_t comparison,
678 : const rewriter& r,
679 : const bool negate=false)
680 38 : {
681 38 : detail::map_based_lhs_t new_lhs;
682 38 : data_expression new_rhs(real_zero());
683 38 : parse_and_store_expression(lhs,new_lhs,new_rhs,r,negate);
684 37 : parse_and_store_expression(rhs,new_lhs,new_rhs,r,!negate);
685 :
686 37 : if (new_lhs.empty())
687 : {
688 4 : if ((comparison==detail::equal && new_rhs==real_zero()) ||
689 9 : (comparison==detail::less_eq && (new_rhs == real_zero() || is_positive(new_rhs,r))) ||
690 1 : (comparison==detail::less && is_positive(new_rhs,r)))
691 : {
692 : // The linear inequality represents true.
693 3 : *this=linear_inequality(detail::lhs_t(),real_one(),detail::less);
694 3 : return;
695 : }
696 : // The linear inequality represents false.
697 1 : *this=linear_inequality(detail::lhs_t(),real_minus_one(),detail::less);
698 1 : return;
699 : }
700 :
701 : // Normalize the linear_inequality such that the first term has factor 1 or -1.
702 33 : data_expression factor=new_lhs.begin()->second;
703 33 : if (factor==real_one() || factor==real_minus_one())
704 : {
705 33 : *this=linear_inequality(detail::map_to_lhs_type(new_lhs),new_rhs,comparison);
706 33 : return;
707 : }
708 0 : factor=rewrite_with_memory(real_abs(factor), r);
709 :
710 0 : *this=linear_inequality(detail::map_to_lhs_type(new_lhs,factor,r),rewrite_with_memory(real_divides(new_rhs,factor), r),comparison);
711 110 : }
712 :
713 :
714 :
715 : /// \brief Constructor that constructs a linear inequality out of a data expression.
716 : /// \details The data expression e is expected to have the form
717 : /// lhs op rhs where op is one of <=,<,==,>,>= and lhs and rhs
718 : /// satisfy the syntax t ::= x | c*t | t*c | t+t | t-t | -t where x is
719 : /// a variable and c is a real constant.
720 : /// \param e Contains the expression to become a linear inequality.
721 : /// \param r A rewriter used to evaluate the expression.
722 :
723 39 : linear_inequality(const data_expression& e,
724 : const rewriter& r)
725 39 : {
726 : detail::comparison_t comparison;
727 39 : bool negate(false);
728 39 : if (is_equal_to_application(e))
729 : {
730 1 : comparison=detail::equal;
731 : }
732 38 : else if (is_less_application(e))
733 : {
734 14 : comparison=detail::less;
735 : }
736 24 : else if (is_less_equal_application(e))
737 : {
738 17 : comparison=detail::less_eq;
739 : }
740 7 : else if (is_greater_application(e))
741 : {
742 0 : comparison=detail::less;
743 0 : negate=true;
744 : }
745 7 : else if (is_greater_equal_application(e))
746 : {
747 6 : comparison=detail::less_eq;
748 6 : negate=true;
749 : }
750 : else
751 : {
752 1 : throw mcrl2::runtime_error("Unexpected equality or inequality: " + pp(e) + "\n") ;
753 : }
754 :
755 38 : data_expression lhs=data::binary_left(atermpp::down_cast<application>(e));
756 38 : data_expression rhs=data::binary_right(atermpp::down_cast<application>(e));
757 38 : *this=linear_inequality(lhs,rhs,comparison,r,negate);
758 :
759 41 : }
760 :
761 :
762 80 : detail::lhs_t::const_iterator lhs_begin() const
763 : {
764 80 : return lhs().begin();
765 : }
766 :
767 83 : detail::lhs_t::const_iterator lhs_end() const
768 : {
769 83 : return lhs().end();
770 : }
771 :
772 639 : const detail::lhs_t& lhs() const
773 : {
774 639 : return atermpp::down_cast<detail::lhs_t>((*this)[0]);
775 : }
776 :
777 49 : const data_expression& rhs() const
778 : {
779 49 : return atermpp::down_cast<data_expression>((*this)[1]);
780 : }
781 :
782 2 : data_expression get_factor_for_a_variable(const variable& x)
783 : {
784 2 : return lhs()[x];
785 : }
786 :
787 116 : detail::comparison_t comparison() const
788 : {
789 116 : if (this->function()==detail::linear_inequality_less())
790 : {
791 49 : return detail::less;
792 : }
793 67 : else if (this->function()==detail::linear_inequality_less_equal())
794 : {
795 66 : return detail::less_eq;
796 : }
797 1 : assert(this->function()==detail::linear_inequality_equal());
798 1 : return detail::equal;
799 : }
800 :
801 3 : data_expression transform_to_data_expression() const
802 : {
803 3 : const detail::comparison_t c=comparison();
804 3 : if (c==detail::less_eq)
805 : {
806 2 : return data::less_equal(lhs().transform_to_data_expression(),rhs());
807 : }
808 1 : if (c==detail::less)
809 : {
810 1 : return data::less(lhs().transform_to_data_expression(),rhs());
811 : }
812 0 : assert(c==detail::equal);
813 0 : return data::equal_to(lhs().transform_to_data_expression(),rhs());
814 : }
815 :
816 101 : bool is_false(const rewriter& r) const
817 : {
818 103 : return lhs().empty() &&
819 3 : ((comparison()==detail::less_eq)?is_negative(rhs(),r):
820 102 : ((comparison()==detail::equal)?!is_zero(rhs()):!is_positive(rhs(),r)));
821 : }
822 :
823 75 : bool is_true(const rewriter& r) const
824 : {
825 79 : return lhs().empty() &&
826 7 : ((comparison()==detail::less_eq)?!is_negative(rhs(),r):
827 78 : ((comparison()==detail::equal)?is_zero(rhs()):is_positive(rhs(),r)));
828 : }
829 :
830 : /// \brief Return this inequality as a typical pair of terms of the form <x1+c2 x2+...+cn xn, d> where c2,...,cn, d are real constants.
831 : /// \return The return value indicates whether the left and right hand side have been negated
832 : /// when yielding the typical pair.
833 : bool typical_pair(
834 : data_expression& lhs_expression,
835 : data_expression& rhs_expression,
836 : detail::comparison_t& comparison_operator,
837 : const rewriter& r) const
838 : {
839 : if (lhs_begin()==lhs_end())
840 : {
841 : lhs_expression=real_zero();
842 : rhs_expression=rhs();
843 : comparison_operator=comparison();
844 : return false;
845 : }
846 :
847 : data_expression factor=lhs_begin()->factor();
848 :
849 : for (detail::lhs_t::const_iterator i=lhs_begin(); i!=lhs_end(); ++i)
850 : {
851 : variable v=i->variable_name();
852 : data_expression e=real_times(rewrite_with_memory(real_divides(i->factor(),factor),r),
853 : data_expression(v));
854 : if (i==lhs_begin())
855 : {
856 : lhs_expression=e;
857 : }
858 : else
859 : {
860 : lhs_expression=real_plus(lhs_expression,e);
861 : }
862 : }
863 :
864 : rhs_expression=rewrite_with_memory(real_divides(rhs(),factor),r);
865 : if (is_negative(factor,r))
866 : {
867 : comparison_operator=negate(comparison());
868 : return true;
869 : }
870 : else
871 : {
872 : comparison_operator=comparison();
873 : return false;
874 : }
875 : }
876 :
877 1 : linear_inequality invert(const rewriter& r)
878 : {
879 1 : const detail::lhs_t new_lhs(lhs().begin(),
880 1 : lhs().end(),
881 1 : [&](const detail::variable_with_a_rational_factor& p)
882 : { return detail::variable_with_a_rational_factor(
883 : p.variable_name(),
884 4 : rewrite_with_memory(real_negate(p.factor()),r));});
885 :
886 1 : const data_expression new_rhs=rewrite_with_memory(real_negate(rhs()),r);
887 1 : if (comparison()==detail::less)
888 : {
889 : // set_comparison(detail::less_eq);
890 1 : return linear_inequality(new_lhs,new_rhs,detail::less_eq);
891 : }
892 0 : else if (comparison()==detail::less_eq)
893 : {
894 : // set_comparison(detail::less);
895 0 : return linear_inequality(new_lhs,new_rhs,detail::less);
896 : }
897 0 : return linear_inequality(new_lhs,new_rhs,detail::equal);
898 1 : }
899 :
900 59 : void add_variables(std::set < variable >& variable_set) const
901 : {
902 133 : for (detail::lhs_t::const_iterator i=lhs().begin(); i!=lhs().end(); ++i)
903 : {
904 74 : variable_set.insert(i->variable_name());
905 : }
906 59 : }
907 : };
908 :
909 0 : inline std::string pp(const linear_inequality& l)
910 : {
911 0 : return pp(l.lhs()) + " " + detail::pp(l.comparison()) + " " + pp(l.rhs());
912 : }
913 :
914 : /// \brief Subtract the given equality, multiplied by f1/f2. The result is e1-(f1/f2)e2,
915 : // where the comparison operator of e1 is kept.
916 1 : inline linear_inequality subtract(const linear_inequality& e1,
917 : const linear_inequality& e2,
918 : const data_expression& f1,
919 : const data_expression& f2,
920 : const rewriter& r)
921 : {
922 1 : const data_expression f=rewrite_with_memory(real_divides(f1,f2),r);
923 : return linear_inequality(
924 2 : detail::meta_operation_lhs(e1.lhs(),e2.lhs(),
925 2 : [&](const data_expression& d1, const data_expression& d2)->data_expression
926 2 : { return real_minus(d1,real_times(f,d2)); },r),
927 2 : rewrite_with_memory(real_minus(e1.rhs(),real_times(f,e2.rhs())),r),
928 : e1.comparison(),
929 3 : r);
930 1 : }
931 :
932 : // Real zero and real one are an ad hoc solution. They should be provided by
933 : // the data type library.
934 :
935 414 : inline data_expression& real_zero()
936 : {
937 414 : static data_expression real_zero=sort_real::real_("0");
938 414 : return real_zero;
939 : }
940 :
941 195 : inline data_expression& real_one()
942 : {
943 195 : static data_expression real_one=sort_real::real_("1");
944 195 : return real_one;
945 : }
946 :
947 59 : inline data_expression& real_minus_one()
948 : {
949 59 : static data_expression real_minus_one=sort_real::real_("-1");
950 59 : return real_minus_one;
951 : }
952 :
953 0 : inline data_expression min(const data_expression& e1,const data_expression& e2,const rewriter& r)
954 : {
955 0 : if (rewrite_with_memory(less_equal(e1,e2),r)==sort_bool::true_())
956 : {
957 0 : return e1;
958 : }
959 0 : if (rewrite_with_memory(less_equal(e2,e1),r)==sort_bool::true_())
960 : {
961 0 : return e2;
962 : }
963 0 : throw mcrl2::runtime_error("Fail to determine the minimum of: " + pp(e1) + " and " + pp(e2) + "\n");
964 : }
965 :
966 0 : inline data_expression max(const data_expression& e1,const data_expression& e2,const rewriter& r)
967 : {
968 0 : if (rewrite_with_memory(less_equal(e2,e1),r)==sort_bool::true_())
969 : {
970 0 : return e1;
971 : }
972 0 : if (rewrite_with_memory(less_equal(e1,e2),r)==sort_bool::true_())
973 : {
974 0 : return e2;
975 : }
976 0 : throw mcrl2::runtime_error("Fail to determine the maximum of: " + pp(e1) + " and " + pp(e2) + "\n");
977 : }
978 :
979 94 : inline bool is_closed_real_number(const data_expression& e)
980 : {
981 94 : if (e.sort()!=sort_real::real_())
982 : {
983 0 : return false;
984 : }
985 :
986 94 : std::set < variable > s=find_all_variables(e);
987 94 : return s.empty();
988 94 : }
989 :
990 15 : inline bool is_negative(const data_expression& e,const rewriter& r)
991 : {
992 15 : data_expression result=rewrite_with_memory(less(e,real_zero()),r);
993 15 : if (result==sort_bool::true_())
994 : {
995 9 : return true;
996 : }
997 6 : if (result==sort_bool::false_())
998 : {
999 6 : return false;
1000 : }
1001 0 : throw mcrl2::runtime_error("Cannot determine that " + pp(e) + " is smaller than 0");
1002 15 : }
1003 :
1004 51 : inline bool is_positive(const data_expression& e,const rewriter& r)
1005 : {
1006 51 : data_expression result=rewrite_with_memory(greater(e,real_zero()),r);
1007 51 : if (result==sort_bool::true_())
1008 : {
1009 22 : return true;
1010 : }
1011 29 : if (result==sort_bool::false_())
1012 : {
1013 29 : return false;
1014 : }
1015 0 : throw mcrl2::runtime_error("Cannot determine that " + pp(e) + " is larger than or equal to 0");
1016 51 : }
1017 :
1018 0 : inline bool is_zero(const data_expression& e)
1019 : {
1020 : // Assume data_expression is in normal form.
1021 0 : assert(is_closed_real_number(e));
1022 0 : return (e==real_zero());
1023 : }
1024 :
1025 : /// \brief Print the vector of inequalities to stderr in readable form.
1026 : template <class TYPE>
1027 0 : inline std::string pp_vector(const TYPE& inequalities)
1028 : {
1029 0 : std::string s="[";
1030 0 : bool first=true;
1031 0 : for (const linear_inequality& l: inequalities)
1032 : {
1033 0 : s=s+ (first?"":", ") + pp(l);
1034 0 : first=false;
1035 : }
1036 0 : s=s+ "]";
1037 0 : return s;
1038 0 : }
1039 :
1040 : bool is_inconsistent(
1041 : const std::vector < linear_inequality >& inequalities_in,
1042 : const rewriter& r,
1043 : const bool use_cache=true);
1044 :
1045 :
1046 : // Count the occurrences of variables that occur in inequalities.
1047 : inline
1048 3 : void count_occurrences(
1049 : const std::vector < linear_inequality >& inequalities,
1050 : std::map < variable, std::size_t>& nr_positive_occurrences,
1051 : std::map < variable, std::size_t>& nr_negative_occurrences,
1052 : const rewriter& r)
1053 : {
1054 3 : for (std::vector < linear_inequality >::const_iterator i=inequalities.begin();
1055 10 : i!=inequalities.end(); ++i)
1056 : {
1057 18 : for (detail::lhs_t::const_iterator j=i->lhs_begin(); j!=i->lhs_end(); ++j)
1058 : {
1059 11 : if (is_positive(j->factor(),r))
1060 : {
1061 5 : nr_positive_occurrences[j->variable_name()]=nr_positive_occurrences[j->variable_name()]+1;
1062 : }
1063 : else
1064 : {
1065 6 : nr_negative_occurrences[j->variable_name()]=nr_negative_occurrences[j->variable_name()]+1;
1066 : }
1067 : }
1068 : }
1069 3 : }
1070 :
1071 : template < class Variable_iterator >
1072 : std::set < variable > gauss_elimination(
1073 : const std::vector < linear_inequality >& inequalities,
1074 : std::vector < linear_inequality >& resulting_equalities,
1075 : std::vector < linear_inequality >& resulting_inequalities,
1076 : Variable_iterator variables_begin,
1077 : Variable_iterator variables_end,
1078 : const rewriter& r);
1079 :
1080 :
1081 :
1082 : /// \brief Indicate whether an inequality from a set of inequalities is redundant.
1083 : /// \details Return whether the inequality referred to by i is inconsistent.
1084 : /// It is expected that i refers to an equality in the vector inequalities.
1085 : /// The vector inequalities might be changed within the procedure, but
1086 : /// will be restored to its original value when this function terminates.
1087 : /// \param inequalities A list of inequalities
1088 : /// \param i An iterator pointing into \a inequalities.
1089 : /// \param r A rewriter
1090 : /// \return An indication whether the inequality referred to by i is inconsistent
1091 : /// in the context of inequalities.
1092 1 : inline bool is_a_redundant_inequality(
1093 : const std::vector < linear_inequality >& inequalities,
1094 : const std::vector < linear_inequality > :: iterator i,
1095 : const rewriter& r)
1096 : {
1097 : #ifndef NDEBUG
1098 : // Check that i points to some position in inequalities.
1099 1 : bool found=false;
1100 1 : for (std::vector < linear_inequality >:: const_iterator j=inequalities.begin() ;
1101 1 : j!=inequalities.end() ; ++j)
1102 : {
1103 1 : if (j==i)
1104 : {
1105 1 : found=true;
1106 1 : break;
1107 : }
1108 : }
1109 1 : assert(found);
1110 : #endif
1111 : // Check whether the inequalities, with the i-th equality with a reversed comparison operator is inconsistent.
1112 : // If yes, the i-th inequality is redundant.
1113 1 : if (i->comparison()==detail::equal)
1114 : {
1115 : // An inequality t==u is only redundant for equalities if
1116 : // t<u and t>u are both inconsistent
1117 0 : const linear_inequality old_inequality=*i;
1118 0 : *i=linear_inequality(i->lhs(),i->rhs(),detail::less);
1119 0 : if (is_inconsistent(inequalities,r))
1120 : {
1121 0 : *i=i->invert(r);
1122 0 : if (is_inconsistent(inequalities,r))
1123 : {
1124 0 : *i=old_inequality;
1125 0 : return true;
1126 : }
1127 : }
1128 0 : *i=old_inequality;
1129 0 : return false;
1130 0 : }
1131 : else
1132 : {
1133 : // an inequality t<u, t<=u, t>u and t>=u is redundant in equalities
1134 : // if its inversion is inconsistent.
1135 1 : const linear_inequality old_inequality=*i;
1136 1 : *i=i->invert(r);
1137 1 : if (is_inconsistent(inequalities,r))
1138 : {
1139 0 : *i=old_inequality;
1140 0 : return true;
1141 : }
1142 : else
1143 : {
1144 1 : *i=old_inequality;
1145 1 : return false;
1146 : }
1147 1 : }
1148 : }
1149 :
1150 : /// \brief Remove every redundant inequality from a vector of inequalities.
1151 : /// \details If inequalities is inconsistent, [false] is returned. Otherwise
1152 : /// a list of inequalities is returned, from which no inequality can
1153 : /// be removed without changing the vector of solutions of the inequalities.
1154 : /// Redundancy of equalities is not checked, because this is quite expensive.
1155 : /// \param inequalities A list of inequalities
1156 : /// \param resulting_inequalities A list of inequalities to which the result is stored.
1157 : // Initially this list must be empty.
1158 : /// \param r A rewriter
1159 :
1160 1 : inline void remove_redundant_inequalities(
1161 : const std::vector < linear_inequality >& inequalities,
1162 : std::vector < linear_inequality >& resulting_inequalities,
1163 : const rewriter& r)
1164 : {
1165 1 : assert(resulting_inequalities.empty());
1166 1 : if (inequalities.empty())
1167 : {
1168 0 : return;
1169 : }
1170 :
1171 : // If false is among the inequalities, [false] is the minimal result.
1172 1 : if (is_inconsistent(inequalities,r))
1173 : {
1174 0 : resulting_inequalities.push_back(linear_inequality());
1175 0 : return;
1176 : }
1177 :
1178 1 : resulting_inequalities=inequalities;
1179 2 : for (std::size_t i=0; i<resulting_inequalities.size();)
1180 : {
1181 : // Check whether the inequalities, with the i-th equality with a reversed comparison operator is inconsistent.
1182 : // If yes, the i-th inequality is redundant.
1183 1 : if (resulting_inequalities[i].comparison()==detail::equal)
1184 : {
1185 : // Do nothing, as removing redundant inequalities is expensive.
1186 0 : ++i;
1187 : }
1188 : else
1189 : {
1190 1 : if (is_a_redundant_inequality(resulting_inequalities,
1191 2 : resulting_inequalities.begin()+i,
1192 : r))
1193 : {
1194 : /* The code below does not preserve the ordering in inequalities.
1195 : if (i+1<resulting_inequalities.size())
1196 : {
1197 : // Copy the last element to the current position.
1198 : resulting_inequalities[i].swap(resulting_inequalities.back());
1199 : }
1200 : resulting_inequalities.pop_back(); */
1201 0 : resulting_inequalities.erase(resulting_inequalities.begin()+i);
1202 : }
1203 : else
1204 : {
1205 1 : ++i;
1206 : }
1207 : }
1208 : }
1209 : }
1210 :
1211 : //---------------------------------------------------------------------------------------------------
1212 :
1213 4 : static void pivot_and_update(
1214 : const variable& xi, // a basic variable
1215 : const variable& xj, // a non basic variable
1216 : const data_expression& v,
1217 : const data_expression& v_delta_correction,
1218 : std::map < variable,data_expression >& beta,
1219 : std::map < variable,data_expression >& beta_delta_correction,
1220 : std::set < variable >& basic_variables,
1221 : std::map < variable, detail::lhs_t >& working_equalities,
1222 : const rewriter& r)
1223 : {
1224 4 : mCRL2log(log::trace) << "Pivoting " << pp(xi) << " " << pp(xj) << "\n";
1225 4 : const data_expression aij=working_equalities[xi][xj];
1226 8 : const data_expression theta=rewrite_with_memory(real_divides(real_minus(v,beta[xi]),aij),r);
1227 8 : const data_expression theta_delta_correction=rewrite_with_memory(real_divides(real_minus(v,beta_delta_correction[xi]),aij),r);
1228 4 : beta[xi]=v;
1229 4 : beta_delta_correction[xi]=v_delta_correction;
1230 4 : beta[xj]=rewrite_with_memory(real_plus(beta[xj],theta),r);
1231 4 : beta_delta_correction[xj]=rewrite_with_memory(real_plus(beta_delta_correction[xj],theta_delta_correction),r);
1232 :
1233 4 : mCRL2log(log::trace) << "Pivoting phase 0\n";
1234 4 : for (std::set < variable >::const_iterator k=basic_variables.begin();
1235 14 : k!=basic_variables.end(); ++k)
1236 : {
1237 10 : mCRL2log(log::trace) << "Working equalities " << *k << ": " << pp(working_equalities[*k]) << "\n";
1238 10 : if ((*k!=xi) && (working_equalities[*k].count(xj)>0))
1239 : {
1240 4 : const data_expression akj=working_equalities[*k][xj];
1241 4 : beta[*k]=rewrite_with_memory(real_plus(beta[*k],real_times(akj ,theta)),r);
1242 4 : beta_delta_correction[*k]=rewrite_with_memory(real_plus(beta_delta_correction[*k],real_times(akj ,theta_delta_correction)),r);
1243 4 : }
1244 : }
1245 : // Apply pivoting on variables xi and xj;
1246 4 : mCRL2log(log::trace) << "Pivoting phase 1\n";
1247 4 : basic_variables.erase(xi);
1248 4 : basic_variables.insert(xj);
1249 :
1250 4 : detail::lhs_t expression_for_xj=working_equalities[xi];
1251 4 : expression_for_xj=expression_for_xj.erase(xj);
1252 4 : expression_for_xj=set_factor_for_a_variable(expression_for_xj,xi,real_minus_one());
1253 4 : expression_for_xj=multiply(expression_for_xj,real_divides(real_minus_one(),aij),r);
1254 4 : mCRL2log(log::trace) << "Expression for xj:" << pp(expression_for_xj) << "\n";
1255 4 : mCRL2log(log::trace) << "Pivoting phase 2\n";
1256 4 : working_equalities.erase(xi);
1257 :
1258 4 : for (std::map < variable, detail::lhs_t >::iterator j=working_equalities.begin();
1259 10 : j!=working_equalities.end(); ++j)
1260 : {
1261 6 : if (j->second.count(xj)>0)
1262 : {
1263 4 : const data_expression factor=j->second[xj];
1264 4 : mCRL2log(log::trace) << "VAR: " << pp(j->first) << " Factor " << pp(factor) << "\n";
1265 4 : j->second=j->second.erase(xj);
1266 : // We need a temporary copy of expression_for_xj as otherwise the multiply
1267 : // below will change this expression.
1268 : //detail::lhs_t temporary_expression_for_xj(expression_for_xj);
1269 4 : j->second=add(j->second,multiply(expression_for_xj,factor,r),r);
1270 4 : }
1271 : }
1272 :
1273 4 : working_equalities[xj]=expression_for_xj;
1274 :
1275 4 : mCRL2log(log::trace) << "Pivoting phase 3\n";
1276 : // Calculate the values for beta and beta_delta_correction for the basic variables
1277 4 : for (std::map < variable, detail::lhs_t >::const_iterator i=working_equalities.begin();
1278 14 : i!=working_equalities.end() ; ++i)
1279 : {
1280 10 : beta[i->first]=i->second.evaluate(make_map_substitution(beta),r);
1281 10 : beta_delta_correction[i->first]=i->second.evaluate(make_map_substitution(beta_delta_correction),r);
1282 : }
1283 :
1284 4 : mCRL2log(log::trace) << "End pivoting " << pp(xj) << "\n";
1285 4 : if (mCRL2logEnabled(log::trace))
1286 : {
1287 0 : for (std::map < variable,data_expression >::const_iterator i=beta.begin();
1288 0 : i!=beta.end(); ++i)
1289 : {
1290 0 : mCRL2log(log::trace) << "(1) beta[" << pp(i->first) << "]= " << pp(beta[i->first]) << "+ delta* " << pp(beta_delta_correction[i->first]) << "\n";
1291 : }
1292 0 : for (std::map < variable, detail::lhs_t >::const_iterator i=working_equalities.begin();
1293 0 : i!=working_equalities.end() ; ++i)
1294 : {
1295 0 : mCRL2log(log::trace) << "EQ: " << pp(i->first) << " := " << pp(i->second) << "\n";
1296 : }
1297 : }
1298 4 : }
1299 :
1300 : namespace detail
1301 : {
1302 : /* False end nodes could be associated with NULL */
1303 : typedef enum { true_end_node, false_end_node, intermediate_node } node_type;
1304 : class inequality_inconsistency_cache;
1305 : class inequality_consistency_cache;
1306 :
1307 : class inequality_inconsistency_cache_base
1308 : {
1309 : friend inequality_inconsistency_cache;
1310 : friend inequality_consistency_cache;
1311 :
1312 : protected:
1313 : node_type m_node;
1314 : linear_inequality m_inequality;
1315 : inequality_inconsistency_cache_base* m_present_branch;
1316 : inequality_inconsistency_cache_base* m_non_present_branch;
1317 :
1318 : public:
1319 :
1320 : inequality_inconsistency_cache_base(const inequality_inconsistency_cache_base& )=delete;
1321 : inequality_inconsistency_cache_base& operator=(const inequality_inconsistency_cache_base& )=delete;
1322 :
1323 24 : inequality_inconsistency_cache_base(const node_type node)
1324 24 : : m_node(node), m_present_branch(), m_non_present_branch()
1325 24 : {}
1326 :
1327 20 : inequality_inconsistency_cache_base(
1328 : const node_type node,
1329 : const linear_inequality& inequality,
1330 : inequality_inconsistency_cache_base* present_branch,
1331 : inequality_inconsistency_cache_base* non_present_branch)
1332 20 : : m_node(node),
1333 20 : m_inequality(inequality),
1334 20 : m_present_branch(present_branch),
1335 20 : m_non_present_branch(non_present_branch)
1336 20 : {}
1337 :
1338 44 : ~inequality_inconsistency_cache_base()
1339 : {
1340 44 : if (m_present_branch!=nullptr)
1341 : {
1342 20 : delete m_present_branch;
1343 : }
1344 44 : if (m_non_present_branch!=nullptr)
1345 : {
1346 20 : delete m_non_present_branch;
1347 : }
1348 44 : }
1349 : };
1350 :
1351 : class inequality_inconsistency_cache
1352 : {
1353 : protected:
1354 : inequality_inconsistency_cache_base* m_cache;
1355 :
1356 : inequality_inconsistency_cache(const inequality_inconsistency_cache& )=delete;
1357 : inequality_inconsistency_cache& operator=(const inequality_consistency_cache& )=delete;
1358 :
1359 : public:
1360 :
1361 1 : inequality_inconsistency_cache()
1362 1 : : m_cache(new inequality_inconsistency_cache_base(false_end_node))
1363 1 : {}
1364 :
1365 1 : ~inequality_inconsistency_cache()
1366 : {
1367 1 : if (m_cache!=nullptr)
1368 : {
1369 1 : delete m_cache;
1370 : }
1371 1 : }
1372 :
1373 11 : bool is_inconsistent(const std::vector < linear_inequality >& inequalities_in_) const
1374 : {
1375 11 : std::set < linear_inequality > inequalities_in(inequalities_in_.begin(),inequalities_in_.end());
1376 11 : const inequality_inconsistency_cache_base* current_root=m_cache;
1377 17 : for(const linear_inequality& l: inequalities_in)
1378 : {
1379 : /* First walk down the three until an endnode is found
1380 : that with l<=current_root.m_inequality. */
1381 24 : while (current_root->m_node==intermediate_node && l>current_root->m_inequality)
1382 : {
1383 9 : current_root=current_root->m_non_present_branch;
1384 : }
1385 15 : if (current_root->m_node==intermediate_node)
1386 : {
1387 6 : if (l==current_root->m_inequality)
1388 : {
1389 6 : current_root=current_root->m_present_branch;
1390 : }
1391 6 : assert(current_root->m_node!=intermediate_node || l<current_root->m_inequality);
1392 : }
1393 : else
1394 : {
1395 9 : return current_root->m_node==true_end_node;
1396 : }
1397 : }
1398 2 : return current_root->m_node==true_end_node;
1399 11 : }
1400 :
1401 2 : void add_inconsistent_inequality_set(const std::vector < linear_inequality >& inequalities_in_)
1402 : {
1403 2 : std::set < linear_inequality > inequalities_in(inequalities_in_.begin(),inequalities_in_.end());
1404 2 : inequality_inconsistency_cache_base** current_root=&m_cache;
1405 6 : for(const linear_inequality& l: inequalities_in)
1406 : {
1407 : /* First walk down the tree until an endnode is found
1408 : that with l<=current_root->m_inequality. */
1409 5 : while ((*current_root)->m_node==intermediate_node && l>(*current_root)->m_inequality)
1410 : {
1411 1 : current_root=&((*current_root)->m_non_present_branch);
1412 : }
1413 4 : if ((*current_root)->m_node==intermediate_node)
1414 : {
1415 1 : if (l==(*current_root)->m_inequality)
1416 : {
1417 1 : current_root=&((*current_root)->m_present_branch);
1418 1 : assert((*current_root)->m_node!=intermediate_node || l<(*current_root)->m_inequality);
1419 : }
1420 : else
1421 : {
1422 : // Add the node.
1423 0 : inequality_inconsistency_cache_base* new_false_node = new inequality_inconsistency_cache_base(false_end_node);
1424 0 : inequality_inconsistency_cache_base* new_node = new inequality_inconsistency_cache_base(intermediate_node,l,new_false_node,*current_root);
1425 0 : *current_root=new_node;
1426 0 : current_root = &(new_node->m_present_branch);
1427 : }
1428 : }
1429 : else
1430 : {
1431 3 : if ((*current_root)->m_node==true_end_node)
1432 : {
1433 : // A shorter sequence than the linear inequality set is already inconsistent.
1434 : // This should not occur, assuming that the linear inequality sequence is checked
1435 : // in in this cache, before being proven inconsistent.
1436 0 : assert(0);
1437 : }
1438 : else
1439 : {
1440 : // Add the remaining sequence.
1441 3 : inequality_inconsistency_cache_base* new_false_node= new inequality_inconsistency_cache_base(false_end_node);
1442 3 : inequality_inconsistency_cache_base* new_node = new inequality_inconsistency_cache_base(intermediate_node,l,new_false_node,*current_root);
1443 3 : *current_root=new_node;
1444 3 : current_root = &(new_node->m_present_branch);
1445 : }
1446 : }
1447 : }
1448 : // At this point the sequence of inequalities has been explored, but the tree has not ended.
1449 : // We expect the current node to be a true_end_node. If not, we replace it by one.
1450 2 : if ((*current_root)->m_node!=true_end_node)
1451 : {
1452 2 : assert(*current_root!=nullptr);
1453 2 : delete *current_root;
1454 2 : *current_root=new inequality_inconsistency_cache_base(true_end_node);
1455 : }
1456 2 : }
1457 : };
1458 :
1459 : class inequality_consistency_cache
1460 : {
1461 : protected:
1462 : inequality_inconsistency_cache_base* m_cache;
1463 :
1464 : inequality_consistency_cache(const inequality_consistency_cache& )=delete;
1465 : inequality_consistency_cache& operator=(const inequality_consistency_cache& )=delete;
1466 :
1467 : public:
1468 :
1469 1 : inequality_consistency_cache()
1470 1 : : m_cache(new inequality_inconsistency_cache_base(false_end_node))
1471 : {
1472 1 : }
1473 :
1474 1 : ~inequality_consistency_cache()
1475 : {
1476 1 : if (m_cache!=nullptr)
1477 : {
1478 1 : delete m_cache;
1479 : }
1480 1 : }
1481 :
1482 : // Sort the vector inequalities_in if not sorted.
1483 16 : bool is_consistent(const std::vector < linear_inequality >& inequalities_in_) const
1484 : {
1485 16 : std::set < linear_inequality > inequalities_in(inequalities_in_.begin(),inequalities_in_.end());
1486 16 : inequality_inconsistency_cache_base* current_root=m_cache;
1487 40 : for(std::set < linear_inequality >::const_iterator i=inequalities_in.begin(); i!=inequalities_in.end(); ++i)
1488 : {
1489 64 : while (current_root->m_node==intermediate_node && *i>current_root->m_inequality)
1490 : {
1491 31 : current_root=current_root->m_non_present_branch;
1492 : }
1493 33 : if (current_root->m_node==intermediate_node)
1494 : {
1495 25 : if (*i==current_root->m_inequality)
1496 : {
1497 24 : current_root=current_root->m_present_branch;
1498 : }
1499 : else
1500 : {
1501 9 : return false; // there are more inequalities than available in the tree. We know nothing about it being consistent.
1502 : }
1503 24 : assert(current_root->m_node!=intermediate_node || *i<current_root->m_inequality);
1504 : }
1505 : else
1506 : {
1507 8 : return current_root->m_node!=false_end_node && i==inequalities_in.end();
1508 : }
1509 : }
1510 7 : return true;
1511 16 : }
1512 :
1513 7 : void add_consistent_inequality_set(const std::vector < linear_inequality >& inequalities_in_)
1514 : {
1515 7 : std::set < linear_inequality > inequalities_in(inequalities_in_.begin(),inequalities_in_.end());
1516 7 : inequality_inconsistency_cache_base** current_root=&m_cache;
1517 27 : for(const linear_inequality& l: inequalities_in)
1518 : {
1519 : /* First walk down the three until an endnode is found
1520 : with l<=current_root->m_inequality. */
1521 35 : while ((*current_root)->m_node==intermediate_node && l>(*current_root)->m_inequality)
1522 : {
1523 15 : current_root=&((*current_root)->m_non_present_branch);
1524 : }
1525 20 : if ((*current_root)->m_node==intermediate_node)
1526 : {
1527 4 : if (l==(*current_root)->m_inequality)
1528 : {
1529 3 : current_root=&((*current_root)->m_present_branch);
1530 3 : assert((*current_root)->m_node!=intermediate_node || l<(*current_root)->m_inequality);
1531 : }
1532 : else
1533 : {
1534 : // Add the node.
1535 1 : inequality_inconsistency_cache_base* new_true_node = new inequality_inconsistency_cache_base(true_end_node);
1536 1 : inequality_inconsistency_cache_base* new_node = new inequality_inconsistency_cache_base(intermediate_node,l,new_true_node,*current_root);
1537 1 : *current_root=new_node;
1538 1 : current_root = &(new_node->m_present_branch);
1539 : }
1540 : }
1541 : else
1542 : {
1543 : // Add the remaining sequence.
1544 16 : inequality_inconsistency_cache_base* new_true_node=new inequality_inconsistency_cache_base(true_end_node);
1545 16 : inequality_inconsistency_cache_base* new_node = new inequality_inconsistency_cache_base(intermediate_node,l,new_true_node,*current_root);
1546 16 : *current_root=new_node;
1547 16 : current_root = &(new_node->m_present_branch);
1548 : }
1549 : }
1550 7 : }
1551 : };
1552 :
1553 : }
1554 :
1555 :
1556 : /// \brief Determine whether a list of data expressions is inconsistent.
1557 : /// \details First it is checked whether false is among the input. If
1558 : /// not, Fourier-Motzkin is applied to all variables in the
1559 : /// inequalities. If the empty vector of equalities is the result,
1560 : /// the input was consistent. Otherwise the resulting vector contains
1561 : /// an inconsistent inequality.
1562 : /// The implementation uses a feasible point detection algorithm as described by
1563 : /// Bruno Dutertre and Leonardo de Moura. Integrating Simplex with DPLL(T).
1564 : /// CSL Technical Report SRI-CSL-06-01, 2006.
1565 : /// \param inequalities_in A list of inequalities.
1566 : /// \param r A rewriter.
1567 : /// \param use_cache A boolean indicating whether results can be cahced.
1568 : /// \return true if the system of inequalities can be determined to be
1569 : /// inconsistent, false otherwise.
1570 11 : inline bool is_inconsistent(
1571 : const std::vector < linear_inequality >& inequalities_in,
1572 : const rewriter& r,
1573 : const bool use_cache/*=true*/)
1574 : {
1575 : // Transform the linear inequalities into a vector of equalities and a
1576 : // sequence of constraints on variables. All variables, including
1577 : // those that will be generated as slack variables will have values indicated
1578 : // by beta, which must lie between the lowerbounds and the upperbounds.
1579 :
1580 : // First remove all equalities by Gauss elimination and make a fresh variable
1581 : // generator.
1582 :
1583 11 : mCRL2log(log::trace) << "Starting an inconsistency check on " + pp_vector(inequalities_in) << "\n";
1584 :
1585 11 : static detail::inequality_inconsistency_cache inconsistency_cache;
1586 11 : static detail::inequality_consistency_cache consistency_cache;
1587 :
1588 11 : if (use_cache && consistency_cache.is_consistent(inequalities_in))
1589 : {
1590 0 : assert(!is_inconsistent(inequalities_in,r,false));
1591 0 : return false;
1592 : }
1593 11 : if (use_cache && inconsistency_cache.is_inconsistent(inequalities_in))
1594 : {
1595 0 : assert(is_inconsistent(inequalities_in,r,false));
1596 0 : return true;
1597 : }
1598 : // The required data structures
1599 11 : std::map < variable,data_expression > lowerbounds;
1600 11 : std::map < variable,data_expression > upperbounds;
1601 11 : std::map < variable,data_expression > beta;
1602 11 : std::map < variable,data_expression > lowerbounds_delta_correction;
1603 11 : std::map < variable,data_expression > upperbounds_delta_correction;
1604 11 : std::map < variable,data_expression > beta_delta_correction;
1605 11 : std::set < variable > non_basic_variables;
1606 11 : std::set < variable > basic_variables;
1607 11 : std::map < variable, detail::lhs_t > working_equalities;
1608 :
1609 11 : set_identifier_generator fresh_variable_name;
1610 :
1611 11 : for (std::vector < linear_inequality >::const_iterator i=inequalities_in.begin();
1612 40 : i!=inequalities_in.end(); ++i)
1613 : {
1614 29 : if (i->is_false(r))
1615 : {
1616 0 : mCRL2log(log::trace) << "Inconsistent, because linear inequalities contains an inconsistent inequality\n";
1617 0 : if (use_cache)
1618 : {
1619 0 : inconsistency_cache.add_inconsistent_inequality_set(inequalities_in); // Necessary? Only false should be added.
1620 0 : assert(inconsistency_cache.is_inconsistent(inequalities_in));
1621 : }
1622 0 : return true;
1623 : }
1624 29 : i->add_variables(non_basic_variables);
1625 29 : for (detail::lhs_t::const_iterator j=i->lhs_begin();
1626 65 : j!=i->lhs_end(); ++j)
1627 : {
1628 36 : fresh_variable_name.add_identifier(j->variable_name().name());
1629 : }
1630 : }
1631 :
1632 11 : std::vector < linear_inequality > inequalities;
1633 11 : std::vector < linear_inequality > equalities;
1634 : non_basic_variables=
1635 22 : gauss_elimination(inequalities_in,
1636 : equalities, // Store all resulting equalities here.
1637 : inequalities, // Store all resulting non equalities here.
1638 : non_basic_variables.begin(),
1639 : non_basic_variables.end(),
1640 11 : r);
1641 :
1642 11 : assert(equalities.size()==0); // There are no resulting equalities left.
1643 11 : non_basic_variables.clear(); // gauss_elimination has removed certain variables. So, we reconstruct the non
1644 : // basic variables again below.
1645 :
1646 11 : mCRL2log(log::trace) << "Resulting equalities " << pp_vector(equalities) << "\n";
1647 11 : mCRL2log(log::trace) << "Resulting inequalities " << pp_vector(inequalities) << "\n";
1648 :
1649 : // Now bring the linear equalities in the basic form described
1650 : // in the article by Bruno Dutertre and Leonardo de Moura.
1651 :
1652 : // First set lower and upperbounds, and introduce slack variables
1653 : // if required.
1654 11 : for (std::vector < linear_inequality >::iterator i=inequalities.begin();
1655 40 : i!=inequalities.end(); ++i)
1656 : {
1657 29 : mCRL2log(log::trace) << "Investigate inequality: " << ":=" << pp(*i) << " || " << pp(i->lhs()) << "\n";
1658 29 : if (i->is_false(r))
1659 : {
1660 0 : mCRL2log(log::trace) << "Inconsistent, because linear inequalities contains an inconsistent inequality after Gauss elimination\n";
1661 0 : if (use_cache)
1662 : {
1663 0 : inconsistency_cache.add_inconsistent_inequality_set(inequalities_in); // Necessary? Only false should be added.
1664 0 : assert(inconsistency_cache.is_inconsistent(inequalities_in));
1665 : }
1666 0 : return true;
1667 : }
1668 29 : if (!i->is_true(r)) // This inequality is redundant and can be skipped.
1669 : {
1670 29 : assert(i->comparison()!=detail::equal);
1671 29 : assert(i->lhs().size()>0); // this signals a redundant or an inconsistent inequality.
1672 29 : i->add_variables(non_basic_variables);
1673 :
1674 29 : if (i->lhs().size()==1) // the left hand side consists of a single variable.
1675 : {
1676 22 : variable v=i->lhs_begin()->variable_name();
1677 22 : data_expression factor=i->lhs_begin()->factor();
1678 22 : assert(factor!=real_zero());
1679 22 : data_expression bound=rewrite_with_memory(real_divides(i->rhs(),factor),r);
1680 22 : if (is_positive(factor,r))
1681 : {
1682 : // The inequality has the shape factor*v<=c or factor*v<c with factor positive
1683 6 : if ((upperbounds.count(v)==0) ||
1684 6 : (rewrite_with_memory(less(bound,upperbounds[v]),r)==sort_bool::true_()))
1685 : {
1686 6 : upperbounds[v]=bound;
1687 6 : upperbounds_delta_correction[v]=
1688 12 : ((i->comparison()==detail::less)?real_minus_one():real_zero());
1689 :
1690 : }
1691 : else
1692 : {
1693 0 : if (bound==upperbounds[v])
1694 : {
1695 0 : upperbounds_delta_correction[v]=
1696 0 : min(upperbounds_delta_correction[v],
1697 0 : ((i->comparison()==detail::less)?real_minus_one():real_zero()),r);
1698 : }
1699 : }
1700 : }
1701 : else
1702 : {
1703 : // The inequality has the shape factor*v<=c or factor*v<c with factor negative
1704 20 : if ((lowerbounds.count(v)==0) ||
1705 20 : (rewrite_with_memory(less(lowerbounds[v],bound),r)==sort_bool::true_()))
1706 : {
1707 13 : lowerbounds[v]=bound;
1708 13 : lowerbounds_delta_correction[v]=
1709 26 : ((i->comparison()==detail::less)?real_one():real_zero());
1710 : }
1711 : else
1712 : {
1713 3 : if (bound==lowerbounds[v])
1714 : {
1715 0 : lowerbounds_delta_correction[v]=
1716 0 : max(lowerbounds_delta_correction[v],
1717 0 : ((i->comparison()==detail::less)?real_one():real_zero()),r);
1718 : }
1719 : }
1720 : }
1721 22 : }
1722 : else
1723 : {
1724 : // The inequality has more than one variable at the left hand side.
1725 : // We transform it into an equation with a new slack variable.
1726 14 : variable new_basic_variable(fresh_variable_name("slack_var"), sort_real::real_());
1727 7 : basic_variables.insert(new_basic_variable);
1728 7 : upperbounds[new_basic_variable]=i->rhs();
1729 7 : upperbounds_delta_correction[new_basic_variable]=
1730 14 : ((i->comparison()==detail::less)?real_minus_one():real_zero());
1731 7 : working_equalities[new_basic_variable]=i->lhs();
1732 7 : mCRL2log(log::trace) << "New slack variable: " << pp(new_basic_variable) << ":=" << pp(*i) << " " << pp(i->lhs()) << "\n";
1733 7 : }
1734 : }
1735 : }
1736 : // Now set the values for beta:
1737 : // The beta values for the non basic variables must satisfy the lower and
1738 : // upperbounds.
1739 11 : for (std::set < variable >::const_iterator i=non_basic_variables.begin();
1740 23 : i!=non_basic_variables.end(); ++i)
1741 : {
1742 15 : if (lowerbounds.count(*i)>0)
1743 : {
1744 21 : if ((upperbounds.count(*i)>0) &&
1745 21 : ((rewrite_with_memory(less(upperbounds[*i],lowerbounds[*i]),r)==sort_bool::true_()) ||
1746 7 : ((upperbounds[*i]==lowerbounds[*i]) &&
1747 15 : (rewrite_with_memory(less(upperbounds_delta_correction[*i],lowerbounds_delta_correction[*i]),r)==sort_bool::true_()))))
1748 : {
1749 3 : mCRL2log(log::trace) << "Inconsistent, preprocessing " << pp(*i) << "\n";
1750 3 : if (use_cache)
1751 : {
1752 2 : inconsistency_cache.add_inconsistent_inequality_set(inequalities_in);
1753 2 : assert(inconsistency_cache.is_inconsistent(inequalities_in));
1754 : }
1755 3 : return true; // Inconsistent.
1756 : }
1757 9 : beta[*i]=lowerbounds[*i];
1758 9 : beta_delta_correction[*i]=lowerbounds_delta_correction[*i];
1759 : }
1760 3 : else if (upperbounds.count(*i)>0)
1761 : {
1762 1 : beta[*i]=upperbounds[*i];
1763 1 : beta_delta_correction[*i]=upperbounds_delta_correction[*i];
1764 : }
1765 : else // *i has neither a lower or an upperbound
1766 : {
1767 2 : beta[*i]=real_zero();
1768 2 : beta_delta_correction[*i]=real_zero();
1769 : }
1770 12 : mCRL2log(log::trace) << "(2) beta[" << pp(*i) << "]=" << pp(beta[*i])<< "+delta*" << pp(beta_delta_correction[*i]) <<"\n";
1771 : }
1772 :
1773 : // Subsequently set the values for the basic variables
1774 8 : for (std::set < variable >::const_iterator i=basic_variables.begin();
1775 15 : i!=basic_variables.end(); ++i)
1776 : {
1777 7 : beta[*i]=working_equalities[*i].evaluate(make_map_substitution(beta),r);
1778 14 : beta_delta_correction[*i]=working_equalities[*i].
1779 14 : evaluate(make_map_substitution(beta_delta_correction),r);
1780 7 : mCRL2log(log::trace) << "(3) beta[" << pp(*i) << "]=" << pp(beta[*i])<< "+delta*" << pp(beta_delta_correction[*i]) <<"\n";
1781 : }
1782 :
1783 : // Now the basic data structure has been set up.
1784 : // We must find the first basic variable that does not satisfy its
1785 : // upper and bounds. This is essentially the check algorithm in the
1786 : // article by Bruno Dutertre and Leonardo de Moura.
1787 :
1788 : for (; true ;)
1789 : {
1790 : // select the smallest basic variable that does not satisfy the bounds.
1791 12 : bool found=false;
1792 12 : bool lowerbound_violation = false;
1793 12 : variable xi;
1794 12 : for (std::set < variable > :: const_iterator i=basic_variables.begin() ;
1795 22 : i!=basic_variables.end() ; ++i)
1796 : {
1797 14 : mCRL2log(log::trace) << "Evaluate start\n";
1798 14 : assert(!found);
1799 14 : data_expression value=beta[*i]; // working_equalities[*i].evaluate(beta,r);
1800 14 : data_expression value_delta_correction=beta_delta_correction[*i]; // working_equalities[*i].evaluate(beta_delta_correction,r);
1801 14 : mCRL2log(log::trace) << "Evaluate end\n";
1802 28 : if ((upperbounds.count(*i)>0) &&
1803 28 : ((rewrite_with_memory(less(upperbounds[*i],value),r)==sort_bool::true_()) ||
1804 6 : ((upperbounds[*i]==value) &&
1805 15 : (rewrite_with_memory(less(upperbounds_delta_correction[*i],value_delta_correction),r)==sort_bool::true_()))))
1806 : {
1807 : // The value of variable *i does not satisfy its upperbound. This must
1808 : // be corrected using pivoting.
1809 4 : mCRL2log(log::trace) << "Upperbound violation " << pp(*i) << " bound: " << pp(upperbounds[*i]) << "\n";
1810 4 : found=true;
1811 4 : xi=*i;
1812 4 : lowerbound_violation=false;
1813 4 : break;
1814 : }
1815 16 : else if ((lowerbounds.count(*i)>0) &&
1816 16 : ((rewrite_with_memory(less(value,lowerbounds[*i]),r)==sort_bool::true_()) ||
1817 3 : ((lowerbounds[*i]==value) &&
1818 10 : (rewrite_with_memory(less(value_delta_correction,lowerbounds_delta_correction[*i]),r)==sort_bool::true_()))))
1819 : {
1820 : // The value of variable *i does not satisfy its upperbound. This must
1821 : // be corrected using pivoting.
1822 0 : mCRL2log(log::trace) << "Lowerbound violation " << pp(*i) << " bound: " << pp(lowerbounds[*i]) << "\n";
1823 0 : found=true;
1824 0 : xi=*i;
1825 0 : lowerbound_violation=true;
1826 0 : break;
1827 : }
1828 18 : }
1829 12 : if (!found)
1830 : {
1831 : // The inequalities are consistent. Return false.
1832 8 : mCRL2log(log::trace) << "Consistent while pivoting\n";
1833 8 : if (use_cache)
1834 : {
1835 7 : consistency_cache.add_consistent_inequality_set(inequalities_in);
1836 7 : assert(consistency_cache.is_consistent(inequalities_in));
1837 : }
1838 8 : return false;
1839 : }
1840 :
1841 4 : mCRL2log(log::trace) << "The smallest basic variable that does not satisfy the bounds is " << pp(xi) << "\n";
1842 4 : if (lowerbound_violation)
1843 : {
1844 0 : mCRL2log(log::trace) << "Lowerbound violation \n";
1845 : // select the smallest non-basic variable with which pivoting can take place.
1846 0 : bool found=false;
1847 0 : const detail::lhs_t& lhs=working_equalities[xi];
1848 0 : for (detail::lhs_t::const_iterator xj_it=lhs.begin(); xj_it!=lhs.end(); ++xj_it)
1849 : {
1850 0 : const variable xj=xj_it->variable_name();
1851 0 : mCRL2log(log::trace) << pp(xj) << " -- " << pp(xj_it->factor()) << "\n";
1852 0 : if ((is_positive(xj_it->factor(),r) &&
1853 0 : ((upperbounds.count(xj)==0) ||
1854 0 : ((rewrite_with_memory(less(beta[xj],upperbounds[xj]),r)==sort_bool::true_())||
1855 0 : ((beta[xj]==upperbounds[xj])&& (rewrite_with_memory(less(beta_delta_correction[xj],upperbounds_delta_correction[xj]),r)==sort_bool::true_()))))) ||
1856 0 : (is_negative(xj_it->factor(),r) &&
1857 0 : ((lowerbounds.count(xj)==0) ||
1858 0 : ((rewrite_with_memory(greater(beta[xj],lowerbounds[xj]),r)==sort_bool::true_())||
1859 0 : ((beta[xj]==lowerbounds[xj]) && (rewrite_with_memory(greater(beta_delta_correction[xj],lowerbounds_delta_correction[xj]),r)==sort_bool::true_()))))))
1860 : {
1861 0 : found=true;
1862 0 : pivot_and_update(xi,xj,lowerbounds[xi],lowerbounds_delta_correction[xi],
1863 : beta, beta_delta_correction,
1864 : basic_variables,working_equalities,r);
1865 0 : break;
1866 : }
1867 0 : }
1868 0 : if (!found)
1869 : {
1870 : // The inequalities are inconsistent.
1871 0 : mCRL2log(log::trace) << "Inconsistent while pivoting\n";
1872 0 : if (use_cache)
1873 : {
1874 0 : inconsistency_cache.add_inconsistent_inequality_set(inequalities_in);
1875 0 : assert(inconsistency_cache.is_inconsistent(inequalities_in));
1876 : }
1877 0 : return true;
1878 : }
1879 :
1880 : }
1881 : else // Upperbound violation.
1882 : {
1883 4 : mCRL2log(log::trace) << "Upperbound violation \n";
1884 : // select the smallest non-basic variable with which pivoting can take place.
1885 4 : bool found=false;
1886 4 : for (detail::lhs_t::const_iterator xj_it=working_equalities[xi].begin();
1887 9 : xj_it!=working_equalities[xi].end(); ++xj_it)
1888 : {
1889 9 : const variable xj=xj_it->variable_name();
1890 9 : mCRL2log(log::trace) << pp(xj) << " -- " << pp(xj_it->factor()) << " POS " << is_positive(xj_it->factor(),r) << "\n";
1891 9 : if ((is_negative(xj_it->factor(),r) &&
1892 6 : ((upperbounds.count(xj)==0) ||
1893 11 : ((rewrite_with_memory(less(beta[xj],upperbounds[xj]),r)==sort_bool::true_()) ||
1894 38 : ((beta[xj]==upperbounds[xj])&& (rewrite_with_memory(less(beta_delta_correction[xj],upperbounds_delta_correction[xj]),r)==sort_bool::true_()))))) ||
1895 5 : (is_positive(xj_it->factor(),r) &&
1896 8 : ((lowerbounds.count(xj)==0) ||
1897 17 : ((rewrite_with_memory(greater(beta[xj],lowerbounds[xj]),r)==sort_bool::true_()) ||
1898 17 : ((beta[xj]==lowerbounds[xj]) && (rewrite_with_memory(greater(beta_delta_correction[xj],lowerbounds_delta_correction[xj]),r)==sort_bool::true_()))))))
1899 : {
1900 4 : found=true;
1901 4 : pivot_and_update(xi,xj,upperbounds[xi],upperbounds_delta_correction[xi],
1902 : beta,beta_delta_correction,
1903 : basic_variables,working_equalities,r);
1904 4 : break;
1905 : }
1906 9 : }
1907 4 : if (!found)
1908 : {
1909 : // The inequalities are inconsistent.
1910 0 : mCRL2log(log::trace) << "Inconsistent while pivoting (1)\n";
1911 0 : if (use_cache)
1912 : {
1913 0 : inconsistency_cache.add_inconsistent_inequality_set(inequalities_in);
1914 0 : assert(inconsistency_cache.is_inconsistent(inequalities_in));
1915 : }
1916 0 : return true;
1917 : }
1918 : }
1919 16 : }
1920 11 : }
1921 :
1922 : //---------------------------------------------------------------------------------------------------
1923 :
1924 :
1925 : /// \brief Try to eliminate variables from a system of inequalities using Gauss elimination.
1926 : /// \details For all variables yi in y1,...,yn indicated by variables_begin to variables_end, it
1927 : /// attempted to find and equation among inequalities of the form yi==expression. All
1928 : /// occurrences of yi in equalities are subsequently replaced by yi. If no equation of
1929 : /// the form yi can be found, yi is added to the list of variables that is returned by
1930 : /// this function. If the input contains an inconsistent inequality, resulting_equalities
1931 : /// becomes empty, resulting_inequalities contains false and the returned list of variables
1932 : /// is also empty. The resulting equalities and inequalities do not contain linear inequalites
1933 : /// equivalent to true.
1934 : /// \param inequalities A list of inequalities over real numbers
1935 : /// \param resulting_equalities A list with the resulting equalities.
1936 : /// \param resulting_inequalities A list of the resulting inequalities
1937 : /// \param variables_begin An iterator indicating the beginning of the eliminatable variables.
1938 : /// \param variables_end An iterator indicating the end of the eliminatable variables.
1939 : /// \param r A rewriter.
1940 : /// \post variables contains the list of variables that have not been eliminated
1941 : /// \return The variables that could not be removed by gauss elimination.
1942 :
1943 : template < class Variable_iterator >
1944 14 : std::set < variable > gauss_elimination(
1945 : const std::vector < linear_inequality >& inequalities,
1946 : std::vector < linear_inequality >& resulting_equalities,
1947 : std::vector < linear_inequality >& resulting_inequalities,
1948 : Variable_iterator variables_begin,
1949 : Variable_iterator variables_end,
1950 : const rewriter& r)
1951 : {
1952 14 : std::set < variable > remaining_variables;
1953 :
1954 : // First copy equalities to the resulting_equalities and the inequalites to resulting_inequalities.
1955 51 : for (std::vector < linear_inequality > ::const_iterator j = inequalities.begin(); j != inequalities.end(); ++j)
1956 : {
1957 37 : if (j->is_false(r))
1958 : {
1959 : // The input contains false. Return false and stop.
1960 0 : resulting_equalities.clear();
1961 0 : resulting_inequalities.clear();
1962 0 : resulting_inequalities.push_back(linear_inequality());
1963 0 : return remaining_variables;
1964 : }
1965 37 : else if (!j->is_true(r)) // Do not consider redundant equations.
1966 : {
1967 37 : if (j->comparison()==detail::equal)
1968 : {
1969 1 : resulting_equalities.push_back(*j);
1970 : }
1971 : else
1972 : {
1973 36 : resulting_inequalities.push_back(*j);
1974 : }
1975 : }
1976 : }
1977 :
1978 : // Now find out whether there are variables that occur in an equality, so
1979 : // that we can perform gauss elimination.
1980 32 : for (Variable_iterator i = variables_begin; i != variables_end; ++i)
1981 : {
1982 : std::size_t j;
1983 20 : for (j=0; j<resulting_equalities.size(); ++j)
1984 : {
1985 1 : bool check_equalities_for_redundant_inequalities(false);
1986 1 : std::set < variable > vars;
1987 1 : resulting_equalities[j].add_variables(vars);
1988 1 : if (vars.count(*i)>0)
1989 : {
1990 : // Equality *j contains data variable *i.
1991 : // Perform gauss elimination, and break the loop.
1992 :
1993 2 : for (std::size_t k = 0; k < resulting_inequalities.size();)
1994 : {
1995 2 : resulting_inequalities[k]=subtract(resulting_inequalities[k],
1996 : resulting_equalities[j],
1997 1 : resulting_inequalities[k].get_factor_for_a_variable(*i),
1998 1 : resulting_equalities[j].get_factor_for_a_variable(*i),
1999 : r);
2000 1 : if (resulting_inequalities[k].is_false(r))
2001 : {
2002 : // The input is inconsistent. Return false.
2003 0 : resulting_equalities.clear();
2004 0 : resulting_inequalities.clear();
2005 0 : resulting_inequalities.push_back(linear_inequality());
2006 0 : remaining_variables.clear();
2007 0 : return remaining_variables;
2008 : }
2009 1 : else if (resulting_inequalities[k].is_true(r))
2010 : {
2011 : // Inequality k has become redundant, and can be removed.
2012 0 : if ((k+1)<resulting_inequalities.size())
2013 : {
2014 0 : resulting_inequalities[k].swap(resulting_inequalities.back());
2015 : }
2016 0 : resulting_inequalities.pop_back();
2017 : }
2018 : else
2019 : {
2020 1 : ++k;
2021 : }
2022 : }
2023 :
2024 2 : for (std::size_t k = 0; k<resulting_equalities.size();)
2025 : {
2026 1 : if (k==j)
2027 : {
2028 1 : ++k;
2029 : }
2030 : else
2031 : {
2032 0 : resulting_equalities[k]=subtract(
2033 : resulting_equalities[k],
2034 : resulting_equalities[j],
2035 0 : resulting_equalities[k].get_factor_for_a_variable(*i),
2036 0 : resulting_equalities[j].get_factor_for_a_variable(*i),
2037 : r);
2038 0 : if (resulting_equalities[k].is_false(r))
2039 : {
2040 : // The input is inconsistent. Return false.
2041 0 : resulting_equalities.clear();
2042 0 : resulting_inequalities.clear();
2043 0 : resulting_inequalities.push_back(linear_inequality());
2044 0 : remaining_variables.clear();
2045 0 : return remaining_variables;
2046 : }
2047 0 : else if (resulting_equalities[k].is_true(r))
2048 : {
2049 : // Equality k has become redundant, and can be removed.
2050 0 : if (j+1==resulting_equalities.size())
2051 : {
2052 : // It is not possible to move move the last element of resulting
2053 : // inequalities to position k, because j is at this last position.
2054 : // Hence, we must recall to check the resulting_equalities for inequalities
2055 : // that are true.
2056 0 : check_equalities_for_redundant_inequalities=true;
2057 : }
2058 : else
2059 : {
2060 0 : if ((k+1)<resulting_equalities.size())
2061 : {
2062 0 : resulting_equalities[k].swap(resulting_equalities.back());
2063 : }
2064 0 : resulting_equalities.pop_back();
2065 : }
2066 : }
2067 : else
2068 : {
2069 0 : ++k;
2070 : }
2071 : }
2072 : }
2073 :
2074 : // Remove equation j.
2075 :
2076 1 : if (j+1<resulting_equalities.size())
2077 : {
2078 0 : resulting_equalities[j].swap(resulting_equalities.back());
2079 : }
2080 1 : resulting_equalities.pop_back();
2081 :
2082 : // If there are unremoved resulting equalities, remove them now.
2083 1 : if (check_equalities_for_redundant_inequalities)
2084 : {
2085 0 : for (std::size_t k = 0; k<resulting_equalities.size();)
2086 : {
2087 0 : if (resulting_equalities[k].is_true(r))
2088 : {
2089 : // Equality k is redundant, and can be removed.
2090 0 : if ((k+1)<resulting_equalities.size())
2091 : {
2092 0 : resulting_equalities[k].swap(resulting_equalities.back());
2093 : }
2094 0 : resulting_equalities.pop_back();
2095 : }
2096 : else
2097 : {
2098 0 : ++k;
2099 : }
2100 : }
2101 : }
2102 : }
2103 : }
2104 18 : remaining_variables.insert(*i);
2105 : }
2106 :
2107 14 : return remaining_variables;
2108 0 : }
2109 :
2110 : // The introduction of the function rewrite_with_memory using a
2111 : // hash table here is a temporary trick, to boost
2112 : // performance, which is slow due to translations necessary from and to
2113 : // rewrite format.
2114 :
2115 392 : inline data_expression rewrite_with_memory(
2116 : const data_expression& t,const rewriter& r)
2117 : {
2118 392 : static std::map < data_expression, data_expression > rewrite_hash_table;
2119 392 : std::map < data_expression, data_expression > :: iterator i=rewrite_hash_table.find(t);
2120 392 : if (i==rewrite_hash_table.end())
2121 : {
2122 110 : data_expression t1=r(t);
2123 110 : rewrite_hash_table.insert(std::make_pair(t, t1));
2124 110 : return t1;
2125 110 : }
2126 282 : return i->second;
2127 : }
2128 :
2129 : } // namespace data
2130 :
2131 : } // namespace mcrl2
2132 :
2133 : #endif // MCRL2_DATA_LINEAR_INEQUALITY_H
|
