1 /********************* */
2 /*! \file nonlinear_extension.cpp
4 ** Top contributors (to current version):
5 ** Tim King, Andrew Reynolds, Andres Noetzli
6 ** This file is part of the CVC4 project.
7 ** Copyright (c) 2009-2017 by the authors listed in the file AUTHORS
8 ** in the top-level source directory) and their institutional affiliations.
9 ** All rights reserved. See the file COPYING in the top-level source
10 ** directory for licensing information.\endverbatim
12 ** \brief [[ Add one-line brief description here ]]
14 ** [[ Add lengthier description here ]]
15 ** \todo document this file
18 #include "theory/arith/nonlinear_extension.h"
23 #include "expr/node_builder.h"
24 #include "options/arith_options.h"
25 #include "theory/arith/arith_msum.h"
26 #include "theory/arith/arith_utilities.h"
27 #include "theory/arith/theory_arith.h"
28 #include "theory/quantifiers/quant_util.h"
29 #include "theory/theory_model.h"
31 using namespace CVC4::kind
;
39 // Return true if a collection c contains an elem k. Compatible with map and set
41 template <class Container
, class Key
>
42 bool Contains(const Container
& c
, const Key
& k
) {
43 return c
.find(k
) != c
.end();
46 // Inserts value into the set/map c if the value was not present there. Returns
47 // true if the value was inserted.
48 template <class Container
, class Value
>
49 bool InsertIfNotPresent(Container
* c
, const Value
& value
) {
50 return (c
->insert(value
)).second
;
53 // Returns true if a vector c contains an elem t.
55 bool IsInVector(const std::vector
<T
>& c
, const T
& t
) {
56 return std::find(c
.begin(), c
.end(), t
) != c
.end();
59 // Returns the a[key] and assertion fails in debug mode.
60 inline unsigned getCount(const NodeMultiset
& a
, Node key
) {
61 NodeMultiset::const_iterator it
= a
.find(key
);
62 Assert(it
!= a
.end());
66 // Returns a[key] if key is in a or value otherwise.
67 unsigned getCountWithDefault(const NodeMultiset
& a
, Node key
, unsigned value
) {
68 NodeMultiset::const_iterator it
= a
.find(key
);
69 return (it
== a
.end()) ? value
: it
->second
;
72 // Returns true if for any key then a[key] <= b[key] where the value for any key
73 // not present is interpreted as 0.
74 bool isSubset(const NodeMultiset
& a
, const NodeMultiset
& b
) {
75 for (NodeMultiset::const_iterator it_a
= a
.begin(); it_a
!= a
.end(); ++it_a
) {
76 Node key
= it_a
->first
;
77 const unsigned a_value
= it_a
->second
;
78 const unsigned b_value
= getCountWithDefault(b
, key
, 0);
79 if (a_value
> b_value
) {
86 // Given two multisets return the multiset difference a \ b.
87 NodeMultiset
diffMultiset(const NodeMultiset
& a
, const NodeMultiset
& b
) {
88 NodeMultiset difference
;
89 for (NodeMultiset::const_iterator it_a
= a
.begin(); it_a
!= a
.end(); ++it_a
) {
90 Node key
= it_a
->first
;
91 const unsigned a_value
= it_a
->second
;
92 const unsigned b_value
= getCountWithDefault(b
, key
, 0);
93 if (a_value
> b_value
) {
94 difference
[key
] = a_value
- b_value
;
100 // Return a vector containing a[key] repetitions of key in a multiset a.
101 std::vector
<Node
> ExpandMultiset(const NodeMultiset
& a
) {
102 std::vector
<Node
> expansion
;
103 for (NodeMultiset::const_iterator it_a
= a
.begin(); it_a
!= a
.end(); ++it_a
) {
104 expansion
.insert(expansion
.end(), it_a
->second
, it_a
->first
);
109 void debugPrintBound(const char* c
, Node coeff
, Node x
, Kind type
, Node rhs
) {
110 Node t
= ArithMSum::mkCoeffTerm(coeff
, x
);
111 Trace(c
) << t
<< " " << type
<< " " << rhs
;
114 struct SortNonlinearExtension
{
115 SortNonlinearExtension()
116 : d_nla(NULL
), d_order_type(0), d_reverse_order(false) {}
117 NonlinearExtension
* d_nla
;
118 unsigned d_order_type
;
119 bool d_reverse_order
;
120 bool operator()(Node i
, Node j
) {
121 int cv
= d_nla
->compare(i
, j
, d_order_type
);
125 return d_reverse_order
? cv
< 0 : cv
> 0;
130 bool hasNewMonomials(Node n
, const std::vector
<Node
>& existing
) {
131 std::set
<Node
> visited
;
133 std::vector
<Node
> worklist
;
134 worklist
.push_back(n
);
135 while (!worklist
.empty()) {
136 Node current
= worklist
.back();
138 if (!Contains(visited
, current
)) {
139 visited
.insert(current
);
140 if (current
.getKind() == NONLINEAR_MULT
)
142 if (!IsInVector(existing
, current
)) {
146 worklist
.insert(worklist
.end(), current
.begin(), current
.end());
155 NonlinearExtension::NonlinearExtension(TheoryArith
& containing
,
156 eq::EqualityEngine
* ee
)
157 : d_builtModel(containing
.getSatContext(), false),
158 d_lemmas(containing
.getUserContext()),
159 d_zero_split(containing
.getUserContext()),
160 d_skolem_atoms(containing
.getUserContext()),
161 d_containing(containing
),
163 d_needsLastCall(false)
165 d_true
= NodeManager::currentNM()->mkConst(true);
166 d_false
= NodeManager::currentNM()->mkConst(false);
167 d_zero
= NodeManager::currentNM()->mkConst(Rational(0));
168 d_one
= NodeManager::currentNM()->mkConst(Rational(1));
169 d_neg_one
= NodeManager::currentNM()->mkConst(Rational(-1));
170 d_two
= NodeManager::currentNM()->mkConst(Rational(2));
171 d_order_points
.push_back(d_neg_one
);
172 d_order_points
.push_back(d_zero
);
173 d_order_points
.push_back(d_one
);
174 d_taylor_real_fv
= NodeManager::currentNM()->mkBoundVar(
175 "x", NodeManager::currentNM()->realType());
176 d_taylor_real_fv_base
= NodeManager::currentNM()->mkBoundVar(
177 "a", NodeManager::currentNM()->realType());
178 d_taylor_real_fv_base_rem
= NodeManager::currentNM()->mkBoundVar(
179 "b", NodeManager::currentNM()->realType());
180 d_taylor_degree
= options::nlExtTfTaylorDegree();
181 d_used_approx
= false;
184 NonlinearExtension::~NonlinearExtension() {}
186 // Returns a reference to either map[key] if it exists in the map
187 // or to a default value otherwise.
189 // Warning: sped_cial care must be taken if value is a temporary object.
190 template <class MapType
, class Key
, class Value
>
191 const Value
& FindWithDefault(const MapType
& map
, const Key
& key
,
192 const Value
& value
) {
193 typename
MapType::const_iterator it
= map
.find(key
);
194 if (it
== map
.end()) {
200 const NodeMultiset
& NonlinearExtension::getMonomialExponentMap(
201 Node monomial
) const {
202 MonomialExponentMap::const_iterator it
= d_m_exp
.find(monomial
);
203 Assert(it
!= d_m_exp
.end());
207 bool NonlinearExtension::isMonomialSubset(Node a
, Node b
) const {
208 const NodeMultiset
& a_exponent_map
= getMonomialExponentMap(a
);
209 const NodeMultiset
& b_exponent_map
= getMonomialExponentMap(b
);
211 return isSubset(a_exponent_map
, b_exponent_map
);
214 void NonlinearExtension::registerMonomialSubset(Node a
, Node b
) {
215 Assert(isMonomialSubset(a
, b
));
217 const NodeMultiset
& a_exponent_map
= getMonomialExponentMap(a
);
218 const NodeMultiset
& b_exponent_map
= getMonomialExponentMap(b
);
220 std::vector
<Node
> diff_children
=
221 ExpandMultiset(diffMultiset(b_exponent_map
, a_exponent_map
));
222 Assert(!diff_children
.empty());
224 d_m_contain_parent
[a
].push_back(b
);
225 d_m_contain_children
[b
].push_back(a
);
227 Node mult_term
= safeConstructNary(MULT
, diff_children
);
228 Node nlmult_term
= safeConstructNary(NONLINEAR_MULT
, diff_children
);
229 d_m_contain_mult
[a
][b
] = mult_term
;
230 d_m_contain_umult
[a
][b
] = nlmult_term
;
231 Trace("nl-ext-mindex") << "..." << a
<< " is a subset of " << b
232 << ", difference is " << mult_term
<< std::endl
;
235 bool NonlinearExtension::getCurrentSubstitution(
236 int effort
, const std::vector
<Node
>& vars
, std::vector
<Node
>& subs
,
237 std::map
<Node
, std::vector
<Node
> >& exp
) {
238 // get the constant equivalence classes
239 std::map
<Node
, std::vector
<int> > rep_to_subs_index
;
242 for (unsigned i
= 0; i
< vars
.size(); i
++) {
244 if (d_ee
->hasTerm(n
)) {
245 Node nr
= d_ee
->getRepresentative(n
);
248 Trace("nl-subs") << "Basic substitution : " << n
<< " -> " << nr
250 exp
[n
].push_back(n
.eqNode(nr
));
253 rep_to_subs_index
[nr
].push_back(i
);
261 // return true if the substitution is non-trivial
264 // d_containing.getValuation().getModel()->getRepresentative( n );
267 std::pair
<bool, Node
> NonlinearExtension::isExtfReduced(
268 int effort
, Node n
, Node on
, const std::vector
<Node
>& exp
) const {
270 Kind k
= n
.getKind();
271 return std::make_pair(k
!= NONLINEAR_MULT
&& !isTranscendentalKind(k
),
275 if (on
.getKind() == NONLINEAR_MULT
)
277 Trace("nl-ext-zero-exp") << "Infer zero : " << on
<< " == " << n
279 // minimize explanation if a substitution+rewrite results in zero
280 const std::set
<Node
> vars(on
.begin(), on
.end());
282 for (unsigned i
= 0, size
= exp
.size(); i
< size
; i
++)
284 Trace("nl-ext-zero-exp") << " exp[" << i
<< "] = " << exp
[i
]
286 std::vector
<Node
> eqs
;
287 if (exp
[i
].getKind() == EQUAL
)
289 eqs
.push_back(exp
[i
]);
291 else if (exp
[i
].getKind() == AND
)
293 for (const Node
& ec
: exp
[i
])
295 if (ec
.getKind() == EQUAL
)
302 for (unsigned j
= 0; j
< eqs
.size(); j
++)
304 for (unsigned r
= 0; r
< 2; r
++)
306 if (eqs
[j
][r
] == d_zero
&& vars
.find(eqs
[j
][1 - r
]) != vars
.end())
308 Trace("nl-ext-zero-exp") << "...single exp : " << eqs
[j
]
310 return std::make_pair(true, eqs
[j
]);
316 return std::make_pair(true, Node::null());
319 Node
NonlinearExtension::computeModelValue(Node n
, unsigned index
) {
320 std::map
<Node
, Node
>::iterator it
= d_mv
[index
].find(n
);
321 if (it
!= d_mv
[index
].end()) {
324 Trace("nl-ext-mv-debug") << "computeModelValue " << n
<< ", index=" << index
<< std::endl
;
329 else if (index
== 1 && (n
.getKind() == NONLINEAR_MULT
330 || isTranscendentalKind(n
.getKind())))
332 if (d_containing
.getValuation().getModel()->hasTerm(n
)) {
333 // use model value for abstraction
334 ret
= d_containing
.getValuation().getModel()->getRepresentative(n
);
336 // abstraction does not exist, use model value
337 //ret = computeModelValue(n, 0);
338 ret
= d_containing
.getValuation().getModel()->getValue(n
);
340 //Assert( ret.isConst() );
341 } else if (n
.getNumChildren() == 0) {
342 if (n
.getKind() == PI
)
344 // we are interested in the exact value of PI, which cannot be computed.
345 // hence, we return PI itself when asked for the concrete value.
348 ret
= d_containing
.getValuation().getModel()->getValue(n
);
351 // otherwise, compute true value
352 std::vector
<Node
> children
;
353 if (n
.getMetaKind() == metakind::PARAMETERIZED
)
355 children
.push_back(n
.getOperator());
357 for (unsigned i
= 0; i
< n
.getNumChildren(); i
++) {
358 Node mc
= computeModelValue(n
[i
], index
);
359 children
.push_back(mc
);
361 ret
= NodeManager::currentNM()->mkNode(n
.getKind(), children
);
362 if (n
.getKind() == APPLY_UF
)
364 ret
= d_containing
.getValuation().getModel()->getValue(ret
);
366 ret
= Rewriter::rewrite(ret
);
369 if (!ret.isConst()) {
370 Trace("nl-ext-debug") << "...got non-constant : " << ret << " for "
371 << n << ", ask model directly." << std::endl;
372 ret = d_containing.getValuation().getModel()->getValue(ret);
376 //if (ret.getType().isReal() && !isArithKind(n.getKind())) {
377 // Trace("nl-ext-mv-debug") << ( index==0 ? "M" : "M_A" ) << "[ " << n
378 // << " ] -> " << ret << std::endl;
379 //may involve transcendental functions
380 //Assert(ret.isConst());
382 Trace("nl-ext-mv-debug") << "computed " << (index
== 0 ? "M" : "M_A") << "["
383 << n
<< "] = " << ret
<< std::endl
;
384 d_mv
[index
][n
] = ret
;
389 void NonlinearExtension::registerMonomial(Node n
) {
390 if (!IsInVector(d_monomials
, n
)) {
391 d_monomials
.push_back(n
);
392 Trace("nl-ext-debug") << "Register monomial : " << n
<< std::endl
;
393 if (n
.getKind() == NONLINEAR_MULT
)
395 // get exponent count
396 for (unsigned k
= 0; k
< n
.getNumChildren(); k
++) {
398 if (k
== 0 || n
[k
] != n
[k
- 1]) {
399 d_m_vlist
[n
].push_back(n
[k
]);
402 d_m_degree
[n
] = n
.getNumChildren();
403 } else if (n
== d_one
) {
405 d_m_vlist
[n
].clear();
408 Assert(!isArithKind(n
.getKind()));
410 d_m_vlist
[n
].push_back(n
);
413 // TODO: sort necessary here?
414 std::sort(d_m_vlist
[n
].begin(), d_m_vlist
[n
].end());
415 Trace("nl-ext-mindex") << "Add monomial to index : " << n
<< std::endl
;
416 d_m_index
.addTerm(n
, d_m_vlist
[n
], this);
420 void NonlinearExtension::setMonomialFactor(Node a
, Node b
,
421 const NodeMultiset
& common
) {
422 // Could not tell if this was being inserted intentionally or not.
423 std::map
<Node
, Node
>& mono_diff_a
= d_mono_diff
[a
];
424 if (!Contains(mono_diff_a
, b
)) {
425 Trace("nl-ext-mono-factor")
426 << "Set monomial factor for " << a
<< "/" << b
<< std::endl
;
427 mono_diff_a
[b
] = mkMonomialRemFactor(a
, common
);
431 void NonlinearExtension::registerConstraint(Node atom
) {
432 if (!IsInVector(d_constraints
, atom
)) {
433 d_constraints
.push_back(atom
);
434 Trace("nl-ext-debug") << "Register constraint : " << atom
<< std::endl
;
435 std::map
<Node
, Node
> msum
;
436 if (ArithMSum::getMonomialSumLit(atom
, msum
))
438 Trace("nl-ext-debug") << "got monomial sum: " << std::endl
;
439 if (Trace
.isOn("nl-ext-debug")) {
440 ArithMSum::debugPrintMonomialSum(msum
, "nl-ext-debug");
442 unsigned max_degree
= 0;
443 std::vector
<Node
> all_m
;
444 std::vector
<Node
> max_deg_m
;
445 for (std::map
<Node
, Node
>::iterator itm
= msum
.begin(); itm
!= msum
.end();
447 if (!itm
->first
.isNull()) {
448 all_m
.push_back(itm
->first
);
449 registerMonomial(itm
->first
);
450 Trace("nl-ext-debug2")
451 << "...process monomial " << itm
->first
<< std::endl
;
452 Assert(d_m_degree
.find(itm
->first
) != d_m_degree
.end());
453 unsigned d
= d_m_degree
[itm
->first
];
454 if (d
> max_degree
) {
458 if (d
>= max_degree
) {
459 max_deg_m
.push_back(itm
->first
);
463 // isolate for each maximal degree monomial
464 for (unsigned i
= 0; i
< all_m
.size(); i
++) {
467 int res
= ArithMSum::isolate(m
, msum
, coeff
, rhs
, atom
.getKind());
469 Kind type
= atom
.getKind();
471 type
= reverseRelationKind(type
);
473 Trace("nl-ext-constraint") << "Constraint : " << atom
<< " <=> ";
474 if (!coeff
.isNull()) {
475 Trace("nl-ext-constraint") << coeff
<< " * ";
477 Trace("nl-ext-constraint")
478 << m
<< " " << type
<< " " << rhs
<< std::endl
;
479 d_c_info
[atom
][m
].d_rhs
= rhs
;
480 d_c_info
[atom
][m
].d_coeff
= coeff
;
481 d_c_info
[atom
][m
].d_type
= type
;
484 for (unsigned i
= 0; i
< max_deg_m
.size(); i
++) {
485 Node m
= max_deg_m
[i
];
486 d_c_info_maxm
[atom
][m
] = true;
489 Trace("nl-ext-debug") << "...failed to get monomial sum." << std::endl
;
494 bool NonlinearExtension::isArithKind(Kind k
) {
495 return k
== PLUS
|| k
== MULT
|| k
== NONLINEAR_MULT
;
498 Node
NonlinearExtension::mkLit(Node a
, Node b
, int status
, int orderType
) {
500 Node a_eq_b
= a
.eqNode(b
);
501 if (orderType
== 0) {
504 // return mkAbs( a ).eqNode( mkAbs( b ) );
505 Node negate_b
= NodeManager::currentNM()->mkNode(UMINUS
, b
);
506 return a_eq_b
.orNode(a
.eqNode(negate_b
));
508 } else if (status
< 0) {
509 return mkLit(b
, a
, -status
);
511 Assert(status
== 1 || status
== 2);
512 NodeManager
* nm
= NodeManager::currentNM();
513 Kind greater_op
= status
== 1 ? GEQ
: GT
;
514 if (orderType
== 0) {
515 return nm
->mkNode(greater_op
, a
, b
);
517 // return nm->mkNode( greater_op, mkAbs( a ), mkAbs( b ) );
518 Node zero
= mkRationalNode(0);
519 Node a_is_nonnegative
= nm
->mkNode(GEQ
, a
, zero
);
520 Node b_is_nonnegative
= nm
->mkNode(GEQ
, b
, zero
);
521 Node negate_a
= nm
->mkNode(UMINUS
, a
);
522 Node negate_b
= nm
->mkNode(UMINUS
, b
);
523 return a_is_nonnegative
.iteNode(
524 b_is_nonnegative
.iteNode(nm
->mkNode(greater_op
, a
, b
),
525 nm
->mkNode(greater_op
, a
, negate_b
)),
526 b_is_nonnegative
.iteNode(nm
->mkNode(greater_op
, negate_a
, b
),
527 nm
->mkNode(greater_op
, negate_a
, negate_b
)));
532 Node
NonlinearExtension::mkAbs(Node a
) {
534 return mkRationalNode(a
.getConst
<Rational
>().abs());
536 NodeManager
* nm
= NodeManager::currentNM();
537 Node a_is_nonnegative
= nm
->mkNode(GEQ
, a
, mkRationalNode(0));
538 return a_is_nonnegative
.iteNode(a
, nm
->mkNode(UMINUS
, a
));
542 Node
NonlinearExtension::mkValidPhase(Node a
, Node pi
) {
544 NodeManager::currentNM()->mkNode(MULT
, mkRationalNode(-1), pi
), a
, pi
);
547 Node
NonlinearExtension::mkBounded( Node l
, Node a
, Node u
) {
548 return NodeManager::currentNM()->mkNode(
550 NodeManager::currentNM()->mkNode(GEQ
, a
, l
),
551 NodeManager::currentNM()->mkNode(LEQ
, a
, u
));
554 // by a <k1> b, a <k2> b, we know a <ret> b
555 Kind
NonlinearExtension::joinKinds(Kind k1
, Kind k2
) {
557 return joinKinds(k2
, k1
);
558 } else if (k1
== k2
) {
561 Assert(k1
== EQUAL
|| k1
== LT
|| k1
== LEQ
|| k1
== GT
|| k1
== GEQ
);
562 Assert(k2
== EQUAL
|| k2
== LT
|| k2
== LEQ
|| k2
== GT
|| k2
== GEQ
);
565 if (k2
== LEQ
|| k2
== GEQ
)
591 return UNDEFINED_KIND
;
595 // by a <k1> b, b <k2> c, we know a <ret> c
596 Kind
NonlinearExtension::transKinds(Kind k1
, Kind k2
) {
598 return transKinds(k2
, k1
);
599 } else if (k1
== k2
) {
602 Assert(k1
== EQUAL
|| k1
== LT
|| k1
== LEQ
|| k1
== GT
|| k1
== GEQ
);
603 Assert(k2
== EQUAL
|| k2
== LT
|| k2
== LEQ
|| k2
== GT
|| k2
== GEQ
);
622 return UNDEFINED_KIND
;
626 bool NonlinearExtension::isTranscendentalKind(Kind k
) {
627 // many operators are eliminated during rewriting
628 Assert(k
!= TANGENT
&& k
!= COSINE
&& k
!= COSECANT
&& k
!= SECANT
630 return k
== EXPONENTIAL
|| k
== SINE
|| k
== PI
;
633 Node
NonlinearExtension::mkMonomialRemFactor(
634 Node n
, const NodeMultiset
& n_exp_rem
) const {
635 std::vector
<Node
> children
;
636 const NodeMultiset
& exponent_map
= getMonomialExponentMap(n
);
637 for (NodeMultiset::const_iterator itme2
= exponent_map
.begin();
638 itme2
!= exponent_map
.end(); ++itme2
) {
639 Node v
= itme2
->first
;
640 unsigned inc
= itme2
->second
;
641 Trace("nl-ext-mono-factor")
642 << "..." << inc
<< " factors of " << v
<< std::endl
;
643 unsigned count_in_n_exp_rem
= getCountWithDefault(n_exp_rem
, v
, 0);
644 Assert(count_in_n_exp_rem
<= inc
);
645 inc
-= count_in_n_exp_rem
;
646 Trace("nl-ext-mono-factor")
647 << "......rem, now " << inc
<< " factors of " << v
<< std::endl
;
648 children
.insert(children
.end(), inc
, v
);
650 Node ret
= safeConstructNary(MULT
, children
);
651 ret
= Rewriter::rewrite(ret
);
652 Trace("nl-ext-mono-factor") << "...return : " << ret
<< std::endl
;
656 int NonlinearExtension::flushLemma(Node lem
) {
657 Trace("nl-ext-lemma-debug")
658 << "NonlinearExtension::Lemma pre-rewrite : " << lem
<< std::endl
;
659 lem
= Rewriter::rewrite(lem
);
660 if (Contains(d_lemmas
, lem
)) {
661 Trace("nl-ext-lemma-debug")
662 << "NonlinearExtension::Lemma duplicate : " << lem
<< std::endl
;
663 // should not generate duplicates
667 d_lemmas
.insert(lem
);
668 Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : " << lem
<< std::endl
;
669 d_containing
.getOutputChannel().lemma(lem
);
673 int NonlinearExtension::flushLemmas(std::vector
<Node
>& lemmas
) {
674 if (options::nlExtEntailConflicts()) {
675 // check if any are entailed to be false
676 for (unsigned i
= 0; i
< lemmas
.size(); i
++) {
677 Node ch_lemma
= lemmas
[i
].negate();
678 ch_lemma
= Rewriter::rewrite(ch_lemma
);
679 Trace("nl-ext-et-debug")
680 << "Check entailment of " << ch_lemma
<< "..." << std::endl
;
681 std::pair
<bool, Node
> et
= d_containing
.getValuation().entailmentCheck(
682 THEORY_OF_TYPE_BASED
, ch_lemma
);
683 Trace("nl-ext-et-debug") << "entailment test result : " << et
.first
<< " "
684 << et
.second
<< std::endl
;
686 Trace("nl-ext-et") << "*** Lemma entailed to be in conflict : "
687 << lemmas
[i
] << std::endl
;
688 // return just this lemma
689 if (flushLemma(lemmas
[i
])) {
698 for (unsigned i
= 0; i
< lemmas
.size(); i
++) {
699 sum
+= flushLemma(lemmas
[i
]);
705 void NonlinearExtension::getAssertions(std::vector
<Node
>& assertions
)
707 Trace("nl-ext") << "Getting assertions..." << std::endl
;
708 NodeManager
* nm
= NodeManager::currentNM();
709 // get the assertions
710 std::map
<Node
, Rational
> init_bounds
[2];
711 std::map
<Node
, Node
> init_bounds_lit
[2];
712 unsigned nassertions
= 0;
713 std::unordered_set
<Node
, NodeHashFunction
> init_assertions
;
714 for (Theory::assertions_iterator it
= d_containing
.facts_begin();
715 it
!= d_containing
.facts_end();
719 const Assertion
& assertion
= *it
;
720 Node lit
= assertion
.assertion
;
721 init_assertions
.insert(lit
);
722 // check for concrete bounds
723 bool pol
= lit
.getKind() != NOT
;
724 Node atom_orig
= lit
.getKind() == NOT
? lit
[0] : lit
;
726 std::vector
<Node
> atoms
;
727 if (atom_orig
.getKind() == EQUAL
)
731 // t = s is ( t >= s ^ t <= s )
732 for (unsigned i
= 0; i
< 2; i
++)
734 Node atom_new
= nm
->mkNode(GEQ
, atom_orig
[i
], atom_orig
[1 - i
]);
735 atom_new
= Rewriter::rewrite(atom_new
);
736 atoms
.push_back(atom_new
);
742 atoms
.push_back(atom_orig
);
745 for (const Node
& atom
: atoms
)
747 // non-strict bounds only
748 if (atom
.getKind() == GEQ
|| (!pol
&& atom
.getKind() == GT
))
751 Assert(atom
[1].isConst());
752 Rational bound
= atom
[1].getConst
<Rational
>();
755 if (atom
[0].getType().isInteger())
757 // ~( p >= c ) ---> ( p <= c-1 )
758 bound
= bound
- Rational(1);
761 unsigned bindex
= pol
? 0 : 1;
762 bool setBound
= true;
763 std::map
<Node
, Rational
>::iterator itb
= init_bounds
[bindex
].find(p
);
764 if (itb
!= init_bounds
[bindex
].end())
766 if (itb
->second
== bound
)
768 setBound
= atom_orig
.getKind() == EQUAL
;
772 setBound
= pol
? itb
->second
< bound
: itb
->second
> bound
;
776 // the bound is subsumed
777 init_assertions
.erase(init_bounds_lit
[bindex
][p
]);
782 Trace("nl-ext-init") << (pol
? "Lower" : "Upper") << " bound for "
783 << p
<< " : " << bound
<< std::endl
;
784 init_bounds
[bindex
][p
] = bound
;
785 init_bounds_lit
[bindex
][p
] = lit
;
790 // for each bound that is the same, ensure we've inferred the equality
791 for (std::pair
<const Node
, Rational
>& ib
: init_bounds
[0])
794 Node lit1
= init_bounds_lit
[0][p
];
795 if (lit1
.getKind() != EQUAL
)
797 std::map
<Node
, Rational
>::iterator itb
= init_bounds
[1].find(p
);
798 if (itb
!= init_bounds
[1].end())
800 if (ib
.second
== itb
->second
)
802 Node eq
= p
.eqNode(nm
->mkConst(ib
.second
));
803 eq
= Rewriter::rewrite(eq
);
804 Node lit2
= init_bounds_lit
[1][p
];
805 Assert(lit2
.getKind() != EQUAL
);
806 // use the equality instead, thus these are redundant
807 init_assertions
.erase(lit1
);
808 init_assertions
.erase(lit2
);
809 init_assertions
.insert(eq
);
815 for (const Node
& a
: init_assertions
)
817 assertions
.push_back(a
);
819 Trace("nl-ext") << "...keep " << assertions
.size() << " / " << nassertions
820 << " assertions." << std::endl
;
823 std::vector
<Node
> NonlinearExtension::checkModelEval(
824 const std::vector
<Node
>& assertions
)
826 std::vector
<Node
> false_asserts
;
827 for (size_t i
= 0; i
< assertions
.size(); ++i
) {
828 Node lit
= assertions
[i
];
829 Node atom
= lit
.getKind()==NOT
? lit
[0] : lit
;
830 if( d_skolem_atoms
.find( atom
)==d_skolem_atoms
.end() ){
831 Node litv
= computeModelValue(lit
);
832 Trace("nl-ext-mv-assert") << "M[[ " << lit
<< " ]] -> " << litv
;
833 if (litv
!= d_true
) {
834 Trace("nl-ext-mv-assert") << " [model-false]" << std::endl
;
835 //Assert(litv == d_false);
836 false_asserts
.push_back(lit
);
838 Trace("nl-ext-mv-assert") << std::endl
;
842 return false_asserts
;
845 bool NonlinearExtension::checkModel(const std::vector
<Node
>& assertions
,
846 const std::vector
<Node
>& false_asserts
)
848 Trace("nl-ext-cm") << "--- check-model ---" << std::endl
;
849 d_check_model_solved
.clear();
850 d_check_model_bounds
.clear();
851 d_check_model_vars
.clear();
852 d_check_model_subs
.clear();
854 // get the presubstitution
855 Trace("nl-ext-cm-debug") << " apply pre-substitution..." << std::endl
;
856 std::vector
<Node
> pvars
;
857 std::vector
<Node
> psubs
;
858 for (std::pair
<const Node
, Node
>& tb
: d_trig_base
)
860 pvars
.push_back(tb
.first
);
861 psubs
.push_back(tb
.second
);
863 // initialize representation of assertions
864 std::vector
<Node
> passertions
;
865 for (const Node
& a
: assertions
)
871 pa
.substitute(pvars
.begin(), pvars
.end(), psubs
.begin(), psubs
.end());
872 pa
= Rewriter::rewrite(pa
);
874 if (!pa
.isConst() || !pa
.getConst
<bool>())
876 Trace("nl-ext-cm-assert") << "- assert : " << pa
<< std::endl
;
877 passertions
.push_back(pa
);
881 // get model bounds for all transcendental functions
882 Trace("nl-ext-cm-debug") << " get bounds for transcendental functions..."
884 for (std::pair
<const Kind
, std::map
<Node
, Node
> >& tfs
: d_tf_rep_map
)
887 for (std::pair
<const Node
, Node
>& tfr
: tfs
.second
)
889 // Figure 3 : tf( x )
890 Node tf
= tfr
.second
;
891 Node atf
= computeModelValue(tf
, 0);
894 addCheckModelBound(atf
, d_pi_bound
[0], d_pi_bound
[1]);
896 else if (isRefineableTfFun(tf
))
898 d_used_approx
= true;
899 std::pair
<Node
, Node
> bounds
= getTfModelBounds(tf
, d_taylor_degree
);
900 addCheckModelBound(atf
, bounds
.first
, bounds
.second
);
902 if (Trace
.isOn("nl-ext-cm"))
904 std::map
<Node
, std::pair
<Node
, Node
> >::iterator it
=
905 d_check_model_bounds
.find(tf
);
906 if (it
!= d_check_model_bounds
.end())
908 Trace("nl-ext-cm") << "check-model-bound : approximate (taylor) : ";
909 printRationalApprox("nl-ext-cm", it
->second
.first
);
910 Trace("nl-ext-cm") << " <= " << tf
<< " <= ";
911 printRationalApprox("nl-ext-cm", it
->second
.second
);
912 Trace("nl-ext-cm") << std::endl
;
918 Trace("nl-ext-cm-debug") << " solve for equalities..." << std::endl
;
919 for (const Node
& atom
: false_asserts
)
921 // see if it corresponds to a univariate polynomial equation of degree two
922 if (atom
.getKind() == EQUAL
)
924 if (!solveEqualitySimple(atom
))
926 // no chance we will satisfy this equality
927 Trace("nl-ext-cm") << "...check-model : failed to solve equality : "
928 << atom
<< std::endl
;
933 // all remaining variables are constrained to their exact model values
934 Trace("nl-ext-cm-debug") << " set exact bounds for remaining variables..."
936 std::unordered_set
<TNode
, TNodeHashFunction
> visited
;
937 std::vector
<TNode
> visit
;
939 for (const Node
& a
: passertions
)
946 if (visited
.find(cur
) == visited
.end())
949 if (cur
.getType().isReal() && !cur
.isConst())
951 Kind k
= cur
.getKind();
952 if (k
!= MULT
&& k
!= PLUS
&& k
!= NONLINEAR_MULT
953 && !isTranscendentalKind(k
))
955 // if we have not set an approximate bound for it
956 if (!hasCheckModelAssignment(cur
))
958 // set its exact model value in the substitution
959 Node curv
= computeModelValue(cur
);
961 << "check-model-bound : exact : " << cur
<< " = ";
962 printRationalApprox("nl-ext-cm", curv
);
963 Trace("nl-ext-cm") << std::endl
;
964 addCheckModelSubstitution(cur
, curv
);
968 for (const Node
& cn
: cur
)
973 } while (!visit
.empty());
976 Trace("nl-ext-cm-debug") << " check assertions..." << std::endl
;
977 std::vector
<Node
> check_assertions
;
978 for (const Node
& a
: passertions
)
980 if (d_check_model_solved
.find(a
) == d_check_model_solved
.end())
983 // apply the substitution to a
984 if (!d_check_model_vars
.empty())
986 av
= av
.substitute(d_check_model_vars
.begin(),
987 d_check_model_vars
.end(),
988 d_check_model_subs
.begin(),
989 d_check_model_subs
.end());
990 av
= Rewriter::rewrite(av
);
992 // simple check literal
993 if (!simpleCheckModelLit(av
))
995 Trace("nl-ext-cm") << "...check-model : assertion failed : " << a
997 check_assertions
.push_back(av
);
998 Trace("nl-ext-cm-debug")
999 << "...check-model : failed assertion, value : " << av
<< std::endl
;
1004 if (!check_assertions
.empty())
1006 Trace("nl-ext-cm") << "...simple check failed." << std::endl
;
1007 // TODO (#1450) check model for general case
1010 Trace("nl-ext-cm") << "...simple check succeeded!" << std::endl
;
1012 // must assert and re-check if produce models is true
1013 if (options::produceModels())
1015 NodeManager
* nm
= NodeManager::currentNM();
1016 // model guard whose semantics is "the model we constructed holds"
1017 Node mg
= nm
->mkSkolem("model", nm
->booleanType());
1018 mg
= Rewriter::rewrite(mg
);
1019 mg
= d_containing
.getValuation().ensureLiteral(mg
);
1020 d_containing
.getOutputChannel().requirePhase(mg
, true);
1021 // assert the constructed model as assertions
1022 for (const std::pair
<const Node
, std::pair
<Node
, Node
> > cb
:
1023 d_check_model_bounds
)
1025 Node l
= cb
.second
.first
;
1026 Node u
= cb
.second
.second
;
1028 Node pred
= nm
->mkNode(AND
, nm
->mkNode(GEQ
, v
, l
), nm
->mkNode(GEQ
, u
, v
));
1029 pred
= nm
->mkNode(OR
, mg
.negate(), pred
);
1030 Trace("nl-ext-lemma-model") << "Assert : " << pred
<< std::endl
;
1031 d_containing
.getOutputChannel().lemma(pred
);
1033 d_builtModel
= true;
1038 void NonlinearExtension::addCheckModelSubstitution(TNode v
, TNode s
)
1040 // should not substitute the same variable twice
1041 Assert(!hasCheckModelAssignment(v
));
1042 for (unsigned i
= 0, size
= d_check_model_subs
.size(); i
< size
; i
++)
1044 Node ms
= d_check_model_subs
[i
];
1045 Node mss
= ms
.substitute(v
, s
);
1048 mss
= Rewriter::rewrite(mss
);
1050 d_check_model_subs
[i
] = mss
;
1052 d_check_model_vars
.push_back(v
);
1053 d_check_model_subs
.push_back(s
);
1056 void NonlinearExtension::addCheckModelBound(TNode v
, TNode l
, TNode u
)
1058 Assert(!hasCheckModelAssignment(v
));
1059 Assert(l
.isConst());
1060 Assert(u
.isConst());
1061 Assert(l
.getConst
<Rational
>() <= u
.getConst
<Rational
>());
1062 d_check_model_bounds
[v
] = std::pair
<Node
, Node
>(l
, u
);
1065 bool NonlinearExtension::hasCheckModelAssignment(Node v
) const
1067 if (d_check_model_bounds
.find(v
) != d_check_model_bounds
.end())
1071 return std::find(d_check_model_vars
.begin(), d_check_model_vars
.end(), v
)
1072 != d_check_model_vars
.end();
1075 bool NonlinearExtension::solveEqualitySimple(Node eq
)
1078 if (!d_check_model_vars
.empty())
1080 seq
= eq
.substitute(d_check_model_vars
.begin(),
1081 d_check_model_vars
.end(),
1082 d_check_model_subs
.begin(),
1083 d_check_model_subs
.end());
1084 seq
= Rewriter::rewrite(seq
);
1087 if (seq
.getConst
<bool>())
1089 d_check_model_solved
[eq
] = Node::null();
1095 Trace("nl-ext-cms") << "simple solve equality " << seq
<< "..." << std::endl
;
1096 Assert(seq
.getKind() == EQUAL
);
1097 std::map
<Node
, Node
> msum
;
1098 if (!ArithMSum::getMonomialSumLit(seq
, msum
))
1100 Trace("nl-ext-cms") << "...fail, could not determine monomial sum."
1104 bool is_valid
= true;
1105 // the variable we will solve a quadratic equation for
1110 NodeManager
* nm
= NodeManager::currentNM();
1111 // the list of variables that occur as a monomial in msum, and whose value
1112 // is so far unconstrained in the model.
1113 std::unordered_set
<Node
, NodeHashFunction
> unc_vars
;
1114 // the list of variables that occur as a factor in a monomial, and whose
1115 // value is so far unconstrained in the model.
1116 std::unordered_set
<Node
, NodeHashFunction
> unc_vars_factor
;
1117 for (std::pair
<const Node
, Node
>& m
: msum
)
1120 Node coeff
= m
.second
.isNull() ? d_one
: m
.second
;
1125 else if (v
.getKind() == NONLINEAR_MULT
)
1127 if (v
.getNumChildren() == 2 && v
[0].isVar() && v
[0] == v
[1]
1128 && (var
.isNull() || var
== v
[0]))
1130 // may solve quadratic
1137 Trace("nl-ext-cms-debug")
1138 << "...invalid due to non-linear monomial " << v
<< std::endl
;
1139 // may wish to set an exact bound for a factor and repeat
1140 for (const Node
& vc
: v
)
1142 unc_vars_factor
.insert(vc
);
1146 else if (!v
.isVar() || (!var
.isNull() && var
!= v
))
1148 Trace("nl-ext-cms-debug")
1149 << "...invalid due to factor " << v
<< std::endl
;
1150 // cannot solve multivariate
1154 // if b is non-zero, then var is also an unconstrained variable
1157 unc_vars
.insert(var
);
1158 unc_vars_factor
.insert(var
);
1161 // if v is unconstrained, we may turn this equality into a substitution
1163 unc_vars_factor
.insert(v
);
1167 // set the variable to solve for
1174 // see if we can solve for a variable?
1175 for (const Node
& uv
: unc_vars
)
1177 Trace("nl-ext-cm-debug") << "check subs var : " << uv
<< std::endl
;
1178 // cannot already have a bound
1179 if (uv
.isVar() && !hasCheckModelAssignment(uv
))
1183 if (ArithMSum::isolate(uv
, msum
, veqc
, slv
, EQUAL
) != 0)
1185 Assert(!slv
.isNull());
1186 // currently do not support substitution-with-coefficients
1187 if (veqc
.isNull() && !slv
.hasSubterm(uv
))
1190 << "check-model-subs : " << uv
<< " -> " << slv
<< std::endl
;
1191 addCheckModelSubstitution(uv
, slv
);
1192 Trace("nl-ext-cms") << "...success, model substitution " << uv
1193 << " -> " << slv
<< std::endl
;
1194 d_check_model_solved
[eq
] = uv
;
1200 // see if we can assign a variable to a constant
1201 for (const Node
& uvf
: unc_vars_factor
)
1203 Trace("nl-ext-cm-debug") << "check set var : " << uvf
<< std::endl
;
1204 // cannot already have a bound
1205 if (uvf
.isVar() && !hasCheckModelAssignment(uvf
))
1207 Node uvfv
= computeModelValue(uvf
);
1208 Trace("nl-ext-cm") << "check-model-bound : exact : " << uvf
<< " = ";
1209 printRationalApprox("nl-ext-cm", uvfv
);
1210 Trace("nl-ext-cm") << std::endl
;
1211 addCheckModelSubstitution(uvf
, uvfv
);
1213 return solveEqualitySimple(eq
);
1216 Trace("nl-ext-cms") << "...fail due to constrained invalid terms."
1220 else if (var
.isNull() || var
.getType().isInteger())
1222 // cannot solve quadratic equations for integer variables
1223 Trace("nl-ext-cms") << "...fail due to variable to solve for." << std::endl
;
1227 // we are linear, it is simple
1232 Trace("nl-ext-cms") << "...fail due to zero a/b." << std::endl
;
1236 Node val
= nm
->mkConst(-c
.getConst
<Rational
>() / b
.getConst
<Rational
>());
1237 Trace("nl-ext-cm") << "check-model-bound : exact : " << var
<< " = ";
1238 printRationalApprox("nl-ext-cm", val
);
1239 Trace("nl-ext-cm") << std::endl
;
1240 addCheckModelSubstitution(var
, val
);
1241 Trace("nl-ext-cms") << "...success, solved linear." << std::endl
;
1242 d_check_model_solved
[eq
] = var
;
1245 Trace("nl-ext-quad") << "Solve quadratic : " << seq
<< std::endl
;
1246 Trace("nl-ext-quad") << " a : " << a
<< std::endl
;
1247 Trace("nl-ext-quad") << " b : " << b
<< std::endl
;
1248 Trace("nl-ext-quad") << " c : " << c
<< std::endl
;
1249 Node two_a
= nm
->mkNode(MULT
, d_two
, a
);
1250 two_a
= Rewriter::rewrite(two_a
);
1251 Node sqrt_val
= nm
->mkNode(
1252 MINUS
, nm
->mkNode(MULT
, b
, b
), nm
->mkNode(MULT
, d_two
, two_a
, c
));
1253 sqrt_val
= Rewriter::rewrite(sqrt_val
);
1254 Trace("nl-ext-quad") << "Will approximate sqrt " << sqrt_val
<< std::endl
;
1255 Assert(sqrt_val
.isConst());
1256 // if it is negative, then we are in conflict
1257 if (sqrt_val
.getConst
<Rational
>().sgn() == -1)
1259 Node conf
= seq
.negate();
1260 Trace("nl-ext-lemma") << "NonlinearExtension::Lemma : quadratic no root : "
1261 << conf
<< std::endl
;
1262 d_containing
.getOutputChannel().lemma(conf
);
1263 Trace("nl-ext-cms") << "...fail due to negative discriminant." << std::endl
;
1266 if (hasCheckModelAssignment(var
))
1268 Trace("nl-ext-cms") << "...fail due to bounds on variable to solve for."
1270 // two quadratic equations for same variable, give up
1273 // approximate the square root of sqrt_val
1275 if (!getApproximateSqrt(sqrt_val
, l
, u
, 15 + d_taylor_degree
))
1277 Trace("nl-ext-cms") << "...fail, could not approximate sqrt." << std::endl
;
1280 d_used_approx
= true;
1281 Trace("nl-ext-quad") << "...got " << l
<< " <= sqrt(" << sqrt_val
1282 << ") <= " << u
<< std::endl
;
1283 Node negb
= nm
->mkConst(-b
.getConst
<Rational
>());
1284 Node coeffa
= nm
->mkConst(Rational(1) / two_a
.getConst
<Rational
>());
1285 // two possible bound regions
1288 Node m_var
= computeModelValue(var
, 0);
1289 Assert(m_var
.isConst());
1290 for (unsigned r
= 0; r
< 2; r
++)
1292 for (unsigned b
= 0; b
< 2; b
++)
1294 Node val
= b
== 0 ? l
: u
;
1295 // (-b +- approx_sqrt( b^2 - 4ac ))/2a
1296 Node approx
= nm
->mkNode(
1297 MULT
, coeffa
, nm
->mkNode(r
== 0 ? MINUS
: PLUS
, negb
, val
));
1298 approx
= Rewriter::rewrite(approx
);
1299 bounds
[r
][b
] = approx
;
1300 Assert(approx
.isConst());
1302 if (bounds
[r
][0].getConst
<Rational
>() > bounds
[r
][1].getConst
<Rational
>())
1304 // ensure bound is (lower, upper)
1305 Node tmp
= bounds
[r
][0];
1306 bounds
[r
][0] = bounds
[r
][1];
1313 nm
->mkConst(Rational(1) / Rational(2)),
1314 nm
->mkNode(PLUS
, bounds
[r
][0], bounds
[r
][1])));
1315 Trace("nl-ext-cm-debug") << "Bound option #" << r
<< " : ";
1316 printRationalApprox("nl-ext-cm-debug", bounds
[r
][0]);
1317 Trace("nl-ext-cm-debug") << "...";
1318 printRationalApprox("nl-ext-cm-debug", bounds
[r
][1]);
1319 Trace("nl-ext-cm-debug") << std::endl
;
1320 diff
= Rewriter::rewrite(diff
);
1321 Assert(diff
.isConst());
1322 diff
= nm
->mkConst(diff
.getConst
<Rational
>().abs());
1323 diff_bound
[r
] = diff
;
1324 Trace("nl-ext-cm-debug") << "...diff from model value (";
1325 printRationalApprox("nl-ext-cm-debug", m_var
);
1326 Trace("nl-ext-cm-debug") << ") is ";
1327 printRationalApprox("nl-ext-cm-debug", diff_bound
[r
]);
1328 Trace("nl-ext-cm-debug") << std::endl
;
1330 // take the one that var is closer to in the model
1331 Node cmp
= nm
->mkNode(GEQ
, diff_bound
[0], diff_bound
[1]);
1332 cmp
= Rewriter::rewrite(cmp
);
1333 Assert(cmp
.isConst());
1334 unsigned r_use_index
= cmp
== d_true
? 1 : 0;
1335 Trace("nl-ext-cm") << "check-model-bound : approximate (sqrt) : ";
1336 printRationalApprox("nl-ext-cm", bounds
[r_use_index
][0]);
1337 Trace("nl-ext-cm") << " <= " << var
<< " <= ";
1338 printRationalApprox("nl-ext-cm", bounds
[r_use_index
][1]);
1339 Trace("nl-ext-cm") << std::endl
;
1340 addCheckModelBound(var
, bounds
[r_use_index
][0], bounds
[r_use_index
][1]);
1341 d_check_model_solved
[eq
] = var
;
1342 Trace("nl-ext-cms") << "...success, solved quadratic." << std::endl
;
1346 bool NonlinearExtension::simpleCheckModelLit(Node lit
)
1348 Trace("nl-ext-cms") << "*** Simple check-model lit for " << lit
<< "..."
1352 Trace("nl-ext-cms") << " return constant." << std::endl
;
1353 return lit
.getConst
<bool>();
1355 NodeManager
* nm
= NodeManager::currentNM();
1356 bool pol
= lit
.getKind() != kind::NOT
;
1357 Node atom
= lit
.getKind() == kind::NOT
? lit
[0] : lit
;
1359 if (atom
.getKind() == EQUAL
)
1361 // x = a is ( x >= a ^ x <= a )
1362 for (unsigned i
= 0; i
< 2; i
++)
1364 Node lit
= nm
->mkNode(GEQ
, atom
[i
], atom
[1 - i
]);
1369 lit
= Rewriter::rewrite(lit
);
1370 bool success
= simpleCheckModelLit(lit
);
1373 // false != true -> one conjunct of equality is false, we fail
1374 // true != false -> one disjunct of disequality is true, we succeed
1378 // both checks passed and polarity is true, or both checks failed and
1379 // polarity is false
1382 else if (atom
.getKind() != GEQ
)
1384 Trace("nl-ext-cms") << " failed due to unknown literal." << std::endl
;
1387 // get the monomial sum
1388 std::map
<Node
, Node
> msum
;
1389 if (!ArithMSum::getMonomialSumLit(atom
, msum
))
1391 Trace("nl-ext-cms") << " failed due to get msum." << std::endl
;
1394 // simple interval analysis
1395 if (simpleCheckModelMsum(msum
, pol
))
1399 // can also try reasoning about univariate quadratic equations
1400 Trace("nl-ext-cms-debug")
1401 << "* Try univariate quadratic analysis..." << std::endl
;
1402 std::vector
<Node
> vs_invalid
;
1403 std::unordered_set
<Node
, NodeHashFunction
> vs
;
1404 std::map
<Node
, Node
> v_a
;
1405 std::map
<Node
, Node
> v_b
;
1406 // get coefficients...
1407 for (std::pair
<const Node
, Node
>& m
: msum
)
1414 v_b
[v
] = m
.second
.isNull() ? d_one
: m
.second
;
1417 else if (v
.getKind() == NONLINEAR_MULT
&& v
.getNumChildren() == 2
1418 && v
[0] == v
[1] && v
[0].isVar())
1420 v_a
[v
[0]] = m
.second
.isNull() ? d_one
: m
.second
;
1425 vs_invalid
.push_back(v
);
1429 // solve the valid variables...
1430 Node invalid_vsum
= vs_invalid
.empty() ? d_zero
1431 : (vs_invalid
.size() == 1
1433 : nm
->mkNode(PLUS
, vs_invalid
));
1434 // substitution to try
1435 std::vector
<Node
> qvars
;
1436 std::vector
<Node
> qsubs
;
1437 for (const Node
& v
: vs
)
1439 // is it a valid variable?
1440 std::map
<Node
, std::pair
<Node
, Node
> >::iterator bit
=
1441 d_check_model_bounds
.find(v
);
1442 if (!invalid_vsum
.hasSubterm(v
) && bit
!= d_check_model_bounds
.end())
1444 std::map
<Node
, Node
>::iterator it
= v_a
.find(v
);
1445 if (it
!= v_a
.end())
1447 Node a
= it
->second
;
1448 Assert(a
.isConst());
1449 int asgn
= a
.getConst
<Rational
>().sgn();
1451 Node t
= nm
->mkNode(MULT
, a
, v
, v
);
1454 if (it
!= v_b
.end())
1457 t
= nm
->mkNode(PLUS
, t
, nm
->mkNode(MULT
, b
, v
));
1459 t
= Rewriter::rewrite(t
);
1460 Trace("nl-ext-cms-debug") << "Trying to find min/max for quadratic "
1461 << t
<< "..." << std::endl
;
1462 Trace("nl-ext-cms-debug") << " a = " << a
<< std::endl
;
1463 Trace("nl-ext-cms-debug") << " b = " << b
<< std::endl
;
1464 // find maximal/minimal value on the interval
1465 Node apex
= nm
->mkNode(
1466 DIVISION
, nm
->mkNode(UMINUS
, b
), nm
->mkNode(MULT
, d_two
, a
));
1467 apex
= Rewriter::rewrite(apex
);
1468 Assert(apex
.isConst());
1469 // for lower, upper, whether we are greater than the apex
1472 for (unsigned r
= 0; r
< 2; r
++)
1474 boundn
[r
] = r
== 0 ? bit
->second
.first
: bit
->second
.second
;
1475 Node cmpn
= nm
->mkNode(GT
, boundn
[r
], apex
);
1476 cmpn
= Rewriter::rewrite(cmpn
);
1477 Assert(cmpn
.isConst());
1478 cmp
[r
] = cmpn
.getConst
<bool>();
1480 Trace("nl-ext-cms-debug") << " apex " << apex
<< std::endl
;
1481 Trace("nl-ext-cms-debug")
1482 << " lower " << boundn
[0] << ", cmp: " << cmp
[0] << std::endl
;
1483 Trace("nl-ext-cms-debug")
1484 << " upper " << boundn
[1] << ", cmp: " << cmp
[1] << std::endl
;
1485 Assert(boundn
[0].getConst
<Rational
>()
1486 <= boundn
[1].getConst
<Rational
>());
1489 if (cmp
[0] != cmp
[1])
1491 Assert(cmp
[0] && !cmp
[1]);
1492 // does the sign match the bound?
1493 if ((asgn
== 1) == pol
)
1495 // the apex is the max/min value
1497 Trace("nl-ext-cms-debug") << " ...set to apex." << std::endl
;
1501 // it is one of the endpoints, plug in and compare
1503 for (unsigned r
= 0; r
< 2; r
++)
1505 qsubs
.push_back(boundn
[r
]);
1506 Node ts
= t
.substitute(
1507 qvars
.begin(), qvars
.end(), qsubs
.begin(), qsubs
.end());
1508 tcmpn
[r
] = Rewriter::rewrite(ts
);
1511 Node tcmp
= nm
->mkNode(LT
, tcmpn
[0], tcmpn
[1]);
1512 Trace("nl-ext-cms-debug")
1513 << " ...both sides of apex, compare " << tcmp
<< std::endl
;
1514 tcmp
= Rewriter::rewrite(tcmp
);
1515 Assert(tcmp
.isConst());
1516 unsigned bindex_use
= (tcmp
.getConst
<bool>() == pol
) ? 1 : 0;
1517 Trace("nl-ext-cms-debug")
1518 << " ...set to " << (bindex_use
== 1 ? "upper" : "lower")
1520 s
= boundn
[bindex_use
];
1525 // both to one side of the apex
1526 // we figure out which bound to use (lower or upper) based on
1528 // (1) whether a's sign is positive,
1529 // (2) whether we are greater than the apex of the parabola,
1530 // (3) the polarity of the constraint, i.e. >= or <=.
1531 // there are 8 cases of these factors, which we test here.
1532 unsigned bindex_use
= (((asgn
== 1) == cmp
[0]) == pol
) ? 0 : 1;
1533 Trace("nl-ext-cms-debug")
1534 << " ...set to " << (bindex_use
== 1 ? "upper" : "lower")
1536 s
= boundn
[bindex_use
];
1538 Assert(!s
.isNull());
1540 Trace("nl-ext-cms") << "* set bound based on quadratic : " << v
1541 << " -> " << s
<< std::endl
;
1547 Assert(qvars
.size() == qsubs
.size());
1549 lit
.substitute(qvars
.begin(), qvars
.end(), qsubs
.begin(), qsubs
.end());
1550 slit
= Rewriter::rewrite(slit
);
1551 return simpleCheckModelLit(slit
);
1556 bool NonlinearExtension::simpleCheckModelMsum(const std::map
<Node
, Node
>& msum
,
1559 Trace("nl-ext-cms-debug") << "* Try simple interval analysis..." << std::endl
;
1560 NodeManager
* nm
= NodeManager::currentNM();
1561 // map from transcendental functions to whether they were set to lower
1563 bool simpleSuccess
= true;
1564 std::map
<Node
, bool> set_bound
;
1565 std::vector
<Node
> sum_bound
;
1566 for (const std::pair
<const Node
, Node
>& m
: msum
)
1571 sum_bound
.push_back(m
.second
.isNull() ? d_one
: m
.second
);
1575 Trace("nl-ext-cms-debug") << "- monomial : " << v
<< std::endl
;
1576 // --- whether we should set a lower bound for this monomial
1578 (m
.second
.isNull() || m
.second
.getConst
<Rational
>().sgn() == 1)
1580 Trace("nl-ext-cms-debug")
1581 << "set bound to " << (set_lower
? "lower" : "upper") << std::endl
;
1583 // --- Collect variables and factors in v
1584 std::vector
<Node
> vars
;
1585 std::vector
<unsigned> factors
;
1586 if (v
.getKind() == NONLINEAR_MULT
)
1588 unsigned last_start
= 0;
1589 for (unsigned i
= 0, nchildren
= v
.getNumChildren(); i
< nchildren
; i
++)
1591 // are we at the end?
1592 if (i
+ 1 == nchildren
|| v
[i
+ 1] != v
[i
])
1594 unsigned vfact
= 1 + (i
- last_start
);
1595 last_start
= (i
+ 1);
1596 vars
.push_back(v
[i
]);
1597 factors
.push_back(vfact
);
1604 factors
.push_back(1);
1607 // --- Get the lower and upper bounds and sign information.
1608 // Whether we have an (odd) number of negative factors in vars, apart
1609 // from the variable at choose_index.
1610 bool has_neg_factor
= false;
1611 int choose_index
= -1;
1612 std::vector
<Node
> ls
;
1613 std::vector
<Node
> us
;
1614 // the relevant sign information for variables with odd exponents:
1615 // 1: both signs of the interval of this variable are positive,
1616 // -1: both signs of the interval of this variable are negative.
1617 std::vector
<int> signs
;
1618 Trace("nl-ext-cms-debug") << "get sign information..." << std::endl
;
1619 for (unsigned i
= 0, size
= vars
.size(); i
< size
; i
++)
1622 unsigned vcfact
= factors
[i
];
1623 if (Trace
.isOn("nl-ext-cms-debug"))
1625 Trace("nl-ext-cms-debug") << "-- " << vc
;
1628 Trace("nl-ext-cms-debug") << "^" << vcfact
;
1630 Trace("nl-ext-cms-debug") << " ";
1632 std::map
<Node
, std::pair
<Node
, Node
> >::iterator bit
=
1633 d_check_model_bounds
.find(vc
);
1634 // if there is a model bound for this term
1635 if (bit
!= d_check_model_bounds
.end())
1637 Node l
= bit
->second
.first
;
1638 Node u
= bit
->second
.second
;
1642 if (vcfact
% 2 == 1)
1645 int lsgn
= l
.getConst
<Rational
>().sgn();
1646 int usgn
= u
.getConst
<Rational
>().sgn();
1647 Trace("nl-ext-cms-debug")
1648 << "bound_sign(" << lsgn
<< "," << usgn
<< ") ";
1653 // must have a negative factor
1654 has_neg_factor
= !has_neg_factor
;
1657 else if (choose_index
== -1)
1659 // set the choose index to this
1665 // ambiguous, can't determine the bound
1667 << " failed due to ambiguious monomial." << std::endl
;
1672 Trace("nl-ext-cms-debug") << " -> " << vsign
<< std::endl
;
1673 signs
.push_back(vsign
);
1677 Trace("nl-ext-cms-debug") << std::endl
;
1679 << " failed due to unknown bound for " << vc
<< std::endl
;
1680 // should either assign a model bound or eliminate the variable
1686 // whether we will try to minimize/maximize (-1/1) the absolute value
1687 int setAbs
= (set_lower
== has_neg_factor
) ? 1 : -1;
1688 Trace("nl-ext-cms-debug")
1689 << "set absolute value to " << (setAbs
== 1 ? "maximal" : "minimal")
1692 std::vector
<Node
> vbs
;
1693 Trace("nl-ext-cms-debug") << "set bounds..." << std::endl
;
1694 for (unsigned i
= 0, size
= vars
.size(); i
< size
; i
++)
1697 unsigned vcfact
= factors
[i
];
1701 int vcsign
= signs
[i
];
1702 Trace("nl-ext-cms-debug")
1703 << "Bounds for " << vc
<< " : " << l
<< ", " << u
1704 << ", sign : " << vcsign
<< ", factor : " << vcfact
<< std::endl
;
1707 // by convention, always say it is lower if they are the same
1708 vc_set_lower
= true;
1709 Trace("nl-ext-cms-debug")
1710 << "..." << vc
<< " equal bound, set to lower" << std::endl
;
1714 if (vcfact
% 2 == 0)
1716 // minimize or maximize its absolute value
1717 Rational la
= l
.getConst
<Rational
>().abs();
1718 Rational ua
= u
.getConst
<Rational
>().abs();
1721 // by convention, always say it is lower if abs are the same
1722 vc_set_lower
= true;
1723 Trace("nl-ext-cms-debug")
1724 << "..." << vc
<< " equal abs, set to lower" << std::endl
;
1728 vc_set_lower
= (la
> ua
) == (setAbs
== 1);
1731 else if (signs
[i
] == 0)
1733 // we choose this index to match the overall set_lower
1734 vc_set_lower
= set_lower
;
1738 vc_set_lower
= (signs
[i
] != setAbs
);
1740 Trace("nl-ext-cms-debug")
1741 << "..." << vc
<< " set to " << (vc_set_lower
? "lower" : "upper")
1744 // check whether this is a conflicting bound
1745 std::map
<Node
, bool>::iterator itsb
= set_bound
.find(vc
);
1746 if (itsb
== set_bound
.end())
1748 set_bound
[vc
] = vc_set_lower
;
1750 else if (itsb
->second
!= vc_set_lower
)
1753 << " failed due to conflicting bound for " << vc
<< std::endl
;
1756 // must over/under approximate based on vc_set_lower, computed above
1757 Node vb
= vc_set_lower
? l
: u
;
1758 for (unsigned i
= 0; i
< vcfact
; i
++)
1767 Node vbound
= vbs
.size() == 1 ? vbs
[0] : nm
->mkNode(MULT
, vbs
);
1768 sum_bound
.push_back(ArithMSum::mkCoeffTerm(m
.second
, vbound
));
1771 // if the exact bound was computed via simple analysis above
1774 if (sum_bound
.size() > 1)
1776 bound
= nm
->mkNode(kind::PLUS
, sum_bound
);
1778 else if (sum_bound
.size() == 1)
1780 bound
= sum_bound
[0];
1786 // make the comparison
1787 Node comp
= nm
->mkNode(kind::GEQ
, bound
, d_zero
);
1790 comp
= comp
.negate();
1792 Trace("nl-ext-cms") << " comparison is : " << comp
<< std::endl
;
1793 comp
= Rewriter::rewrite(comp
);
1794 Assert(comp
.isConst());
1795 Trace("nl-ext-cms") << " returned : " << comp
<< std::endl
;
1796 return comp
== d_true
;
1799 std::vector
<Node
> NonlinearExtension::checkSplitZero() {
1800 std::vector
<Node
> lemmas
;
1801 for (unsigned i
= 0; i
< d_ms_vars
.size(); i
++) {
1802 Node v
= d_ms_vars
[i
];
1803 if (d_zero_split
.insert(v
)) {
1804 Node eq
= v
.eqNode(d_zero
);
1805 eq
= Rewriter::rewrite(eq
);
1806 Node literal
= d_containing
.getValuation().ensureLiteral(eq
);
1807 d_containing
.getOutputChannel().requirePhase(literal
, true);
1808 lemmas
.push_back(literal
.orNode(literal
.negate()));
1814 int NonlinearExtension::checkLastCall(const std::vector
<Node
>& assertions
,
1815 const std::vector
<Node
>& false_asserts
,
1816 const std::vector
<Node
>& xts
)
1822 d_m_nconst_factor
.clear();
1823 d_tplane_refine
.clear();
1827 d_tf_rep_map
.clear();
1828 d_tf_region
.clear();
1829 d_waiting_lemmas
.clear();
1831 int lemmas_proc
= 0;
1832 std::vector
<Node
> lemmas
;
1833 NodeManager
* nm
= NodeManager::currentNM();
1835 Trace("nl-ext-mv") << "Extended terms : " << std::endl
;
1836 // register the extended function terms
1837 std::map
< Node
, Node
> mvarg_to_term
;
1838 std::vector
<Node
> trig_no_base
;
1839 bool needPi
= false;
1840 for( unsigned i
=0; i
<xts
.size(); i
++ ){
1842 computeModelValue(a
, 0);
1843 computeModelValue(a
, 1);
1844 printModelValue("nl-ext-mv", a
);
1845 //Assert(d_mv[1][a].isConst());
1846 //Assert(d_mv[0][a].isConst());
1847 if (a
.getKind() == NONLINEAR_MULT
)
1849 d_ms
.push_back( a
);
1851 //context-independent registration
1852 registerMonomial(a
);
1854 std::map
<Node
, std::vector
<Node
> >::iterator itvl
= d_m_vlist
.find(a
);
1855 Assert(itvl
!= d_m_vlist
.end());
1856 for (unsigned k
= 0; k
< itvl
->second
.size(); k
++) {
1857 if (!IsInVector(d_ms_vars
, itvl
->second
[k
])) {
1858 d_ms_vars
.push_back(itvl
->second
[k
]);
1860 Node mvk
= computeModelValue( itvl
->second
[k
], 1 );
1861 if( !mvk
.isConst() ){
1862 d_m_nconst_factor
[a
] = true;
1866 //mark processed if has a "one" factor (will look at reduced monomial)
1867 std::map< Node, std::map< Node, unsigned > >::iterator itme =
1868 d_m_exp.find( a ); Assert( itme!=d_m_exp.end() ); for( std::map< Node,
1869 unsigned >::iterator itme2 = itme->second.begin(); itme2 !=
1870 itme->second.end(); ++itme2 ){ Node v = itme->first; Assert(
1871 d_mv[0].find( v )!=d_mv[0].end() ); Node mvv = d_mv[0][ v ]; if(
1872 mvv==d_one || mvv==d_neg_one ){ ms_proc[ a ] = true;
1874 << "...mark " << a << " reduced since has 1 factor." << std::endl;
1879 }else if( a
.getNumChildren()==1 ){
1880 bool consider
= true;
1881 // get shifted version
1882 if (a
.getKind() == SINE
)
1884 if( d_trig_is_base
.find( a
)==d_trig_is_base
.end() ){
1886 trig_no_base
.push_back(a
);
1891 Node r
= d_containing
.getValuation().getModel()->getRepresentative(a
[0]);
1892 std::map
< Node
, Node
>::iterator itrm
= d_tf_rep_map
[a
.getKind()].find( r
);
1893 if( itrm
!=d_tf_rep_map
[a
.getKind()].end() ){
1894 //verify they have the same model value
1895 if( d_mv
[1][a
]!=d_mv
[1][itrm
->second
] ){
1896 // if not, add congruence lemma
1897 Node cong_lemma
= nm
->mkNode(
1898 IMPLIES
, a
[0].eqNode(itrm
->second
[0]), a
.eqNode(itrm
->second
));
1899 lemmas
.push_back( cong_lemma
);
1903 d_tf_rep_map
[a
.getKind()][r
] = a
;
1907 else if (a
.getKind() == PI
)
1910 d_tf_rep_map
[a
.getKind()][a
] = a
;
1917 // initialize pi if necessary
1918 if (needPi
&& d_pi
.isNull())
1921 getCurrentPiBounds(lemmas
);
1924 lemmas_proc
= flushLemmas(lemmas
);
1925 if (lemmas_proc
> 0) {
1926 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas during registration." << std::endl
;
1930 // process SINE phase shifting
1931 for (const Node
& a
: trig_no_base
)
1933 if (d_trig_base
.find(a
) == d_trig_base
.end())
1936 nm
->mkSkolem("y", nm
->realType(), "phase shifted trigonometric arg");
1937 Node new_a
= nm
->mkNode(a
.getKind(), y
);
1938 d_trig_is_base
[new_a
] = true;
1939 d_trig_base
[a
] = new_a
;
1940 Trace("nl-ext-tf") << "Basis sine : " << new_a
<< " for " << a
1942 Assert(!d_pi
.isNull());
1943 Node shift
= nm
->mkSkolem("s", nm
->integerType(), "number of shifts");
1944 // FIXME : do not introduce shift here, instead needs model-based
1945 // refinement for constant shifts (#1284)
1946 Node shift_lem
= nm
->mkNode(
1948 mkValidPhase(y
, d_pi
),
1951 mkValidPhase(a
[0], d_pi
),
1953 a
[0].eqNode(nm
->mkNode(
1956 nm
->mkNode(MULT
, nm
->mkConst(Rational(2)), shift
, d_pi
)))),
1958 // must do preprocess on this one
1959 Trace("nl-ext-lemma")
1960 << "NonlinearExtension::Lemma : shift : " << shift_lem
<< std::endl
;
1961 d_containing
.getOutputChannel().lemma(shift_lem
, false, true);
1965 if (lemmas_proc
> 0)
1967 Trace("nl-ext") << " ...finished with " << lemmas_proc
1968 << " new lemmas SINE phase shifting." << std::endl
;
1971 Trace("nl-ext") << "We have " << d_ms
.size() << " monomials." << std::endl
;
1973 // register constants
1974 registerMonomial(d_one
);
1975 for (unsigned j
= 0; j
< d_order_points
.size(); j
++) {
1976 Node c
= d_order_points
[j
];
1977 computeModelValue(c
, 0);
1978 computeModelValue(c
, 1);
1981 // register variables
1982 Trace("nl-ext-mv") << "Variables in monomials : " << std::endl
;
1983 for (unsigned i
= 0; i
< d_ms_vars
.size(); i
++) {
1984 Node v
= d_ms_vars
[i
];
1985 registerMonomial(v
);
1986 computeModelValue(v
, 0);
1987 computeModelValue(v
, 1);
1988 printModelValue("nl-ext-mv", v
);
1990 if (Trace
.isOn("nl-ext-mv"))
1992 Trace("nl-ext-mv") << "Arguments of trancendental functions : "
1994 for (std::map
<Kind
, std::map
<Node
, Node
> >::iterator it
=
1995 d_tf_rep_map
.begin();
1996 it
!= d_tf_rep_map
.end();
2000 if (k
== SINE
|| k
== EXPONENTIAL
)
2002 for (std::map
<Node
, Node
>::iterator itt
= it
->second
.begin();
2003 itt
!= it
->second
.end();
2006 Node v
= itt
->second
[0];
2007 computeModelValue(v
, 0);
2008 computeModelValue(v
, 1);
2009 printModelValue("nl-ext-mv", v
);
2015 //----------------------------------- possibly split on zero
2016 if (options::nlExtSplitZero()) {
2017 Trace("nl-ext") << "Get zero split lemmas..." << std::endl
;
2018 lemmas
= checkSplitZero();
2019 lemmas_proc
= flushLemmas(lemmas
);
2020 if (lemmas_proc
> 0) {
2021 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2026 //-----------------------------------initial lemmas for transcendental functions
2027 lemmas
= checkTranscendentalInitialRefine();
2028 lemmas_proc
= flushLemmas(lemmas
);
2029 if (lemmas_proc
> 0) {
2030 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2034 //-----------------------------------lemmas based on sign (comparison to zero)
2035 lemmas
= checkMonomialSign();
2036 lemmas_proc
= flushLemmas(lemmas
);
2037 if (lemmas_proc
> 0) {
2038 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2042 //-----------------------------------monotonicity of transdental functions
2043 lemmas
= checkTranscendentalMonotonic();
2044 lemmas_proc
= flushLemmas(lemmas
);
2045 if (lemmas_proc
> 0) {
2046 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2050 //-----------------------------------lemmas based on magnitude of non-zero monomials
2051 Trace("nl-ext-proc") << "Assign order ids..." << std::endl
;
2053 assignOrderIds(d_ms_vars
, d_order_vars
, r
);
2055 // sort individual variable lists
2056 Trace("nl-ext-proc") << "Assign order var lists..." << std::endl
;
2057 SortNonlinearExtension smv
;
2059 smv
.d_order_type
= r
;
2060 smv
.d_reverse_order
= true;
2061 for (unsigned j
= 0; j
< d_ms
.size(); j
++) {
2062 std::sort(d_m_vlist
[d_ms
[j
]].begin(), d_m_vlist
[d_ms
[j
]].end(), smv
);
2064 for (unsigned c
= 0; c
< 3; c
++) {
2065 // c is effort level
2066 lemmas
= checkMonomialMagnitude( c
);
2067 unsigned nlem
= lemmas
.size();
2068 lemmas_proc
= flushLemmas(lemmas
);
2069 if (lemmas_proc
> 0) {
2070 Trace("nl-ext") << " ...finished with " << lemmas_proc
2071 << " new lemmas (out of possible " << nlem
<< ")."
2077 // sort monomials by degree
2078 Trace("nl-ext-proc") << "Sort monomials by degree..." << std::endl
;
2079 SortNonlinearExtension snlad
;
2081 snlad
.d_order_type
= 4;
2082 snlad
.d_reverse_order
= false;
2083 std::sort(d_ms
.begin(), d_ms
.end(), snlad
);
2085 d_mterms
.insert(d_mterms
.end(), d_ms_vars
.begin(), d_ms_vars
.end());
2086 d_mterms
.insert(d_mterms
.end(), d_ms
.begin(), d_ms
.end());
2088 //-----------------------------------inferred bounds lemmas
2089 // e.g. x >= t => y*x >= y*t
2090 std::vector
< Node
> nt_lemmas
;
2091 lemmas
= checkMonomialInferBounds(nt_lemmas
, assertions
, false_asserts
);
2092 // Trace("nl-ext") << "Bound lemmas : " << lemmas.size() << ", " <<
2093 // nt_lemmas.size() << std::endl; prioritize lemmas that do not
2094 // introduce new monomials
2095 lemmas_proc
= flushLemmas(lemmas
);
2097 if (options::nlExtTangentPlanes() && options::nlExtTangentPlanesInterleave())
2099 lemmas
= checkTangentPlanes();
2100 lemmas_proc
+= flushLemmas(lemmas
);
2103 if (lemmas_proc
> 0) {
2104 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas."
2109 // from inferred bound inferences : now do ones that introduce new terms
2110 lemmas_proc
= flushLemmas(nt_lemmas
);
2111 if (lemmas_proc
> 0) {
2112 Trace("nl-ext") << " ...finished with " << lemmas_proc
2113 << " new (monomial-introducing) lemmas." << std::endl
;
2117 //------------------------------------factoring lemmas
2118 // x*y + x*z >= t => exists k. k = y + z ^ x*k >= t
2119 if( options::nlExtFactor() ){
2120 lemmas
= checkFactoring(assertions
, false_asserts
);
2121 lemmas_proc
= flushLemmas(lemmas
);
2122 if (lemmas_proc
> 0) {
2123 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2128 //------------------------------------resolution bound inferences
2129 // e.g. ( y>=0 ^ s <= x*z ^ x*y <= t ) => y*s <= z*t
2130 if (options::nlExtResBound()) {
2131 lemmas
= checkMonomialInferResBounds();
2132 lemmas_proc
= flushLemmas(lemmas
);
2133 if (lemmas_proc
> 0) {
2134 Trace("nl-ext") << " ...finished with " << lemmas_proc
<< " new lemmas." << std::endl
;
2139 //------------------------------------tangent planes
2140 if (options::nlExtTangentPlanes() && !options::nlExtTangentPlanesInterleave())
2142 lemmas
= checkTangentPlanes();
2143 d_waiting_lemmas
.insert(
2144 d_waiting_lemmas
.end(), lemmas
.begin(), lemmas
.end());
2147 if (options::nlExtTfTangentPlanes())
2149 lemmas
= checkTranscendentalTangentPlanes();
2150 d_waiting_lemmas
.insert(
2151 d_waiting_lemmas
.end(), lemmas
.begin(), lemmas
.end());
2154 Trace("nl-ext") << " ...finished with " << d_waiting_lemmas
.size()
2155 << " waiting lemmas." << std::endl
;
2160 void NonlinearExtension::check(Theory::Effort e
) {
2161 Trace("nl-ext") << std::endl
;
2162 Trace("nl-ext") << "NonlinearExtension::check, effort = " << e
<< std::endl
;
2163 if (d_builtModel
.get())
2165 if (e
== Theory::EFFORT_FULL
)
2169 // now, record the approximations we used
2170 NodeManager
* nm
= NodeManager::currentNM();
2171 for (const std::pair
<const Node
, std::pair
<Node
, Node
> >& cb
:
2172 d_check_model_bounds
)
2174 Node l
= cb
.second
.first
;
2175 Node u
= cb
.second
.second
;
2180 nm
->mkNode(AND
, nm
->mkNode(GEQ
, v
, l
), nm
->mkNode(GEQ
, u
, v
));
2181 pred
= Rewriter::rewrite(pred
);
2182 d_containing
.getValuation().getModel()->recordApproximation(v
, pred
);
2187 if (e
== Theory::EFFORT_FULL
) {
2188 d_containing
.getExtTheory()->clearCache();
2189 d_needsLastCall
= true;
2190 if (options::nlExtRewrites()) {
2191 std::vector
<Node
> nred
;
2192 if (!d_containing
.getExtTheory()->doInferences(0, nred
)) {
2193 Trace("nl-ext") << "...sent no lemmas, # extf to reduce = "
2194 << nred
.size() << std::endl
;
2196 d_needsLastCall
= false;
2199 Trace("nl-ext") << "...sent lemmas." << std::endl
;
2203 // get the assertions
2204 std::vector
<Node
> assertions
;
2205 getAssertions(assertions
);
2207 // reset cached information
2211 Trace("nl-ext-mv-assert")
2212 << "Getting model values... check for [model-false]" << std::endl
;
2213 // get the assertions that are false in the model
2214 const std::vector
<Node
> false_asserts
= checkModelEval(assertions
);
2216 // get the extended terms belonging to this theory
2217 std::vector
<Node
> xts
;
2218 d_containing
.getExtTheory()->getTerms(xts
);
2220 if (Trace
.isOn("nl-ext-debug"))
2222 Trace("nl-ext-debug")
2223 << " processing NonlinearExtension::check : " << std::endl
;
2224 Trace("nl-ext-debug") << " " << false_asserts
.size()
2225 << " false assertions" << std::endl
;
2226 Trace("nl-ext-debug")
2227 << " " << xts
.size() << " extended terms: " << std::endl
;
2228 Trace("nl-ext-debug") << " ";
2229 for (unsigned j
= 0; j
< xts
.size(); j
++)
2231 Trace("nl-ext-debug") << xts
[j
] << " ";
2233 Trace("nl-ext-debug") << std::endl
;
2236 // compute whether shared terms have correct values
2237 unsigned num_shared_wrong_value
= 0;
2238 std::vector
<Node
> shared_term_value_splits
;
2239 // must ensure that shared terms are equal to their concrete value
2240 Trace("nl-ext-mv") << "Shared terms : " << std::endl
;
2241 for (context::CDList
<TNode
>::const_iterator its
=
2242 d_containing
.shared_terms_begin();
2243 its
!= d_containing
.shared_terms_end();
2246 TNode shared_term
= *its
;
2247 // compute its value in the model, and its evaluation in the model
2248 Node stv0
= computeModelValue(shared_term
, 0);
2249 Node stv1
= computeModelValue(shared_term
, 1);
2250 printModelValue("nl-ext-mv", shared_term
);
2253 num_shared_wrong_value
++;
2254 Trace("nl-ext-mv") << "Bad shared term value : " << shared_term
2256 if (shared_term
!= stv0
)
2258 // split on the value, this is non-terminating in general, TODO :
2260 Node eq
= shared_term
.eqNode(stv0
);
2261 shared_term_value_splits
.push_back(eq
);
2265 // this can happen for transcendental functions
2266 // the problem is that we cannot evaluate transcendental functions
2267 // (they don't have a rewriter that returns constants)
2268 // thus, the actual value in their model can be themselves, hence we
2269 // have no reference point to rule out the current model. In this
2270 // case, we may set incomplete below.
2274 Trace("nl-ext-debug") << " " << num_shared_wrong_value
2275 << " shared terms with wrong model value."
2280 d_used_approx
= false;
2281 needsRecheck
= false;
2282 Assert(e
== Theory::EFFORT_LAST_CALL
);
2284 // 1 : we may answer SAT, -1 : we may not answer SAT, 0 : unknown
2285 int complete_status
= 1;
2286 int num_added_lemmas
= 0;
2287 // we require a check either if an assertion is false or a shared term has
2289 if (!false_asserts
.empty() || num_shared_wrong_value
> 0)
2291 complete_status
= num_shared_wrong_value
> 0 ? -1 : 0;
2292 num_added_lemmas
= checkLastCall(assertions
, false_asserts
, xts
);
2293 if (num_added_lemmas
> 0)
2298 Trace("nl-ext") << "Finished check with status : " << complete_status
2301 // if we did not add a lemma during check and there is a chance for SAT
2302 if (complete_status
== 0)
2305 << "Checking model based on bounds for transcendental functions..."
2307 // check the model based on simple solving of equalities and using
2308 // error bounds on the Taylor approximation of transcendental functions.
2309 if (checkModel(assertions
, false_asserts
))
2311 complete_status
= 1;
2315 // if we have not concluded SAT
2316 if (complete_status
!= 1)
2318 // flush the waiting lemmas
2319 num_added_lemmas
= flushLemmas(d_waiting_lemmas
);
2320 if (num_added_lemmas
> 0)
2322 Trace("nl-ext") << "...added " << num_added_lemmas
2323 << " waiting lemmas." << std::endl
;
2326 // resort to splitting on shared terms with their model value
2327 // if we did not add any lemmas
2328 if (num_shared_wrong_value
> 0)
2330 complete_status
= -1;
2331 if (!shared_term_value_splits
.empty())
2333 std::vector
<Node
> shared_term_value_lemmas
;
2334 for (const Node
& eq
: shared_term_value_splits
)
2336 Node req
= Rewriter::rewrite(eq
);
2337 Node literal
= d_containing
.getValuation().ensureLiteral(req
);
2338 d_containing
.getOutputChannel().requirePhase(literal
, true);
2339 Trace("nl-ext-debug") << "Split on : " << literal
<< std::endl
;
2340 shared_term_value_lemmas
.push_back(
2341 literal
.orNode(literal
.negate()));
2343 num_added_lemmas
= flushLemmas(shared_term_value_lemmas
);
2344 if (num_added_lemmas
> 0)
2347 << "...added " << num_added_lemmas
2348 << " shared term value split lemmas." << std::endl
;
2354 // this can happen if we are trying to do theory combination with
2355 // trancendental functions
2356 // since their model value cannot even be computed exactly
2360 // we are incomplete
2361 if (options::nlExtIncPrecision() && d_used_approx
)
2364 d_used_approx
= false;
2365 needsRecheck
= true;
2366 // increase precision for PI?
2367 // Difficult since Taylor series is very slow to converge
2368 Trace("nl-ext") << "...increment Taylor degree to " << d_taylor_degree
2373 Trace("nl-ext") << "...failed to send lemma in "
2374 "NonLinearExtension, set incomplete"
2376 d_containing
.getOutputChannel().setIncomplete();
2380 } while (needsRecheck
);
2384 void NonlinearExtension::presolve()
2386 Trace("nl-ext") << "NonlinearExtension::presolve" << std::endl
;
2389 void NonlinearExtension::assignOrderIds(std::vector
<Node
>& vars
,
2390 NodeMultiset
& order
,
2391 unsigned orderType
) {
2392 SortNonlinearExtension smv
;
2394 smv
.d_order_type
= orderType
;
2395 smv
.d_reverse_order
= false;
2396 std::sort(vars
.begin(), vars
.end(), smv
);
2399 // assign ordering id's
2400 unsigned counter
= 0;
2401 unsigned order_index
= (orderType
== 0 || orderType
== 2) ? 0 : 1;
2403 for (unsigned j
= 0; j
< vars
.size(); j
++) {
2405 Node v
= get_compare_value(x
, orderType
);
2407 Trace("nl-ext-mvo") << "..do not assign order to " << x
<< " : " << v
<< std::endl
;
2408 //don't assign for non-constant values (transcendental function apps)
2411 Trace("nl-ext-mvo") << " order " << x
<< " : " << v
<< std::endl
;
2417 if (order_index
< d_order_points
.size()) {
2418 Node vv
= get_compare_value(d_order_points
[order_index
], orderType
);
2419 if (compare_value(v
, vv
, orderType
) <= 0) {
2422 << "O_" << orderType
<< "[" << d_order_points
[order_index
]
2423 << "] = " << counter
<< std::endl
;
2424 order
[d_order_points
[order_index
]] = counter
;
2432 if (prev
.isNull() || compare_value(v
, prev
, orderType
) != 0) {
2435 Trace("nl-ext-mvo") << "O_" << orderType
<< "[" << x
<< "] = " << counter
2440 while (order_index
< d_order_points
.size()) {
2442 Trace("nl-ext-mvo") << "O_" << orderType
<< "["
2443 << d_order_points
[order_index
] << "] = " << counter
2445 order
[d_order_points
[order_index
]] = counter
;
2450 int NonlinearExtension::compare(Node i
, Node j
, unsigned orderType
) const {
2451 Assert(orderType
>= 0);
2452 if (orderType
<= 3) {
2453 Node ci
= get_compare_value(i
, orderType
);
2454 Node cj
= get_compare_value(j
, orderType
);
2457 return compare_value(ci
, cj
, orderType
);
2469 } else if (orderType
== 4) {
2470 unsigned i_count
= getCount(d_m_degree
, i
);
2471 unsigned j_count
= getCount(d_m_degree
, j
);
2472 return i_count
== j_count
? 0 : (i_count
< j_count
? 1 : -1);
2480 void NonlinearExtension::mkPi(){
2481 if( d_pi
.isNull() ){
2482 d_pi
= NodeManager::currentNM()->mkNullaryOperator(
2483 NodeManager::currentNM()->realType(), PI
);
2484 d_pi_2
= Rewriter::rewrite(NodeManager::currentNM()->mkNode(
2487 NodeManager::currentNM()->mkConst(Rational(1) / Rational(2))));
2488 d_pi_neg_2
= Rewriter::rewrite(NodeManager::currentNM()->mkNode(
2491 NodeManager::currentNM()->mkConst(Rational(-1) / Rational(2))));
2492 d_pi_neg
= Rewriter::rewrite(NodeManager::currentNM()->mkNode(
2493 MULT
, d_pi
, NodeManager::currentNM()->mkConst(Rational(-1))));
2495 d_pi_bound
[0] = NodeManager::currentNM()->mkConst( Rational(333)/Rational(106) );
2496 d_pi_bound
[1] = NodeManager::currentNM()->mkConst( Rational(355)/Rational(113) );
2500 void NonlinearExtension::getCurrentPiBounds( std::vector
< Node
>& lemmas
) {
2501 Node pi_lem
= NodeManager::currentNM()->mkNode(
2503 NodeManager::currentNM()->mkNode(GEQ
, d_pi
, d_pi_bound
[0]),
2504 NodeManager::currentNM()->mkNode(LEQ
, d_pi
, d_pi_bound
[1]));
2505 lemmas
.push_back( pi_lem
);
2508 Node
NonlinearExtension::getApproximateConstant(Node c
,
2510 unsigned prec
) const
2512 Assert(c
.isConst());
2513 Rational cr
= c
.getConst
<Rational
>();
2516 unsigned upper
= pow(10, prec
);
2518 Rational den
= Rational(upper
);
2519 if (cr
.getDenominator() < den
.getNumerator())
2521 // denominator is not more than precision, we return it
2525 int csign
= cr
.sgn();
2531 Rational one
= Rational(1);
2532 Rational ten
= Rational(10);
2533 Rational pow_ten
= Rational(1);
2534 // inefficient for large numbers
2538 pow_ten
= pow_ten
* ten
;
2540 Rational allow_err
= one
/ den
;
2542 Trace("nl-ext-approx") << "Compute approximation for " << c
<< ", precision "
2543 << prec
<< "..." << std::endl
;
2544 // now do binary search
2545 Rational two
= Rational(2);
2546 NodeManager
* nm
= NodeManager::currentNM();
2550 unsigned curr
= (lower
+ upper
) / 2;
2551 Rational curr_r
= Rational(curr
) / den
;
2552 Rational err
= cr
- curr_r
;
2553 int esign
= err
.sgn();
2554 if (err
.abs() <= allow_err
)
2556 if (esign
== 1 && !isLower
)
2558 curr_r
= Rational(curr
+ 1) / den
;
2560 else if (esign
== -1 && isLower
)
2562 curr_r
= Rational(curr
- 1) / den
;
2564 curr_r
= curr_r
* pow_ten
;
2565 cret
= nm
->mkConst(csign
== 1 ? curr_r
: -curr_r
);
2570 // update lower/upper
2575 else if (esign
== 1)
2580 } while (cret
.isNull());
2581 Trace("nl-ext-approx") << "Approximation for " << c
<< " for precision "
2582 << prec
<< " is " << cret
<< std::endl
;
2586 bool NonlinearExtension::getApproximateSqrt(Node c
,
2589 unsigned iter
) const
2591 Assert(c
.isConst());
2592 if (c
== d_one
|| c
== d_zero
)
2598 Rational rc
= c
.getConst
<Rational
>();
2600 Rational rl
= rc
< Rational(1) ? rc
: Rational(1);
2601 Rational ru
= rc
< Rational(1) ? Rational(1) : rc
;
2603 Rational half
= Rational(1) / Rational(2);
2604 while (count
< iter
)
2606 Rational curr
= half
* (rl
+ ru
);
2607 Rational curr_sq
= curr
* curr
;
2614 else if (curr_sq
< rc
)
2625 NodeManager
* nm
= NodeManager::currentNM();
2626 l
= nm
->mkConst(rl
);
2627 u
= nm
->mkConst(ru
);
2631 void NonlinearExtension::printRationalApprox(const char* c
,
2633 unsigned prec
) const
2635 Assert(cr
.isConst());
2636 Node ca
= getApproximateConstant(cr
, true, prec
);
2644 Trace(c
) << " [0,10^" << prec
<< "])";
2648 void NonlinearExtension::printModelValue(const char* c
,
2650 unsigned prec
) const
2654 Trace(c
) << " " << n
<< " -> ";
2655 for (unsigned i
= 0; i
< 2; i
++)
2657 std::map
<Node
, Node
>::const_iterator it
= d_mv
[1 - i
].find(n
);
2658 Assert(it
!= d_mv
[1 - i
].end());
2659 if (it
->second
.isConst())
2661 printRationalApprox(c
, it
->second
, prec
);
2665 Trace(c
) << "?"; // it->second;
2667 Trace(c
) << (i
== 0 ? " [actual: " : " ]");
2669 Trace(c
) << std::endl
;
2673 int NonlinearExtension::compare_value(Node i
, Node j
,
2674 unsigned orderType
) const {
2675 Assert(orderType
>= 0 && orderType
<= 3);
2676 Assert( i
.isConst() && j
.isConst() );
2677 Trace("nl-ext-debug") << "compare value " << i
<< " " << j
2678 << ", o = " << orderType
<< std::endl
;
2682 } else if (orderType
== 0 || orderType
== 2) {
2683 ret
= i
.getConst
<Rational
>() < j
.getConst
<Rational
>() ? 1 : -1;
2685 Trace("nl-ext-debug") << "vals : " << i
.getConst
<Rational
>() << " "
2686 << j
.getConst
<Rational
>() << std::endl
;
2687 Trace("nl-ext-debug") << i
.getConst
<Rational
>().abs() << " "
2688 << j
.getConst
<Rational
>().abs() << std::endl
;
2689 ret
= (i
.getConst
<Rational
>().abs() == j
.getConst
<Rational
>().abs()
2691 : (i
.getConst
<Rational
>().abs() < j
.getConst
<Rational
>().abs()
2695 Trace("nl-ext-debug") << "...return " << ret
<< std::endl
;
2699 Node
NonlinearExtension::get_compare_value(Node i
, unsigned orderType
) const {
2704 Trace("nl-ext-debug") << "Compare variable " << i
<< " " << orderType
2706 Assert(orderType
>= 0 && orderType
<= 3);
2707 unsigned mindex
= orderType
<= 1 ? 0 : 1;
2708 std::map
<Node
, Node
>::const_iterator iti
= d_mv
[mindex
].find(i
);
2709 Assert(iti
!= d_mv
[mindex
].end());
2713 // show a <> 0 by inequalities between variables in monomial a w.r.t 0
2714 int NonlinearExtension::compareSign(Node oa
, Node a
, unsigned a_index
,
2715 int status
, std::vector
<Node
>& exp
,
2716 std::vector
<Node
>& lem
) {
2717 Trace("nl-ext-debug") << "Process " << a
<< " at index " << a_index
2718 << ", status is " << status
<< std::endl
;
2719 Assert(d_mv
[1].find(oa
) != d_mv
[1].end());
2720 if (a_index
== d_m_vlist
[a
].size()) {
2721 if (d_mv
[1][oa
].getConst
<Rational
>().sgn() != status
) {
2723 safeConstructNary(AND
, exp
).impNode(mkLit(oa
, d_zero
, status
* 2));
2724 lem
.push_back(lemma
);
2728 Assert(a_index
< d_m_vlist
[a
].size());
2729 Node av
= d_m_vlist
[a
][a_index
];
2730 unsigned aexp
= d_m_exp
[a
][av
];
2731 // take current sign in model
2732 Assert( d_mv
[1].find(av
) != d_mv
[0].end() );
2733 Assert( d_mv
[1][av
].isConst() );
2734 int sgn
= d_mv
[1][av
].getConst
<Rational
>().sgn();
2735 Trace("nl-ext-debug") << "Process var " << av
<< "^" << aexp
2736 << ", model sign = " << sgn
<< std::endl
;
2738 if (d_mv
[1][oa
].getConst
<Rational
>().sgn() != 0) {
2739 Node lemma
= av
.eqNode(d_zero
).impNode(oa
.eqNode(d_zero
));
2740 lem
.push_back(lemma
);
2744 if (aexp
% 2 == 0) {
2745 exp
.push_back(av
.eqNode(d_zero
).negate());
2746 return compareSign(oa
, a
, a_index
+ 1, status
, exp
, lem
);
2749 NodeManager::currentNM()->mkNode(sgn
== 1 ? GT
: LT
, av
, d_zero
));
2750 return compareSign(oa
, a
, a_index
+ 1, status
* sgn
, exp
, lem
);
2756 bool NonlinearExtension::compareMonomial(
2757 Node oa
, Node a
, NodeMultiset
& a_exp_proc
, Node ob
, Node b
,
2758 NodeMultiset
& b_exp_proc
, std::vector
<Node
>& exp
, std::vector
<Node
>& lem
,
2759 std::map
<int, std::map
<Node
, std::map
<Node
, Node
> > >& cmp_infers
) {
2760 Trace("nl-ext-comp-debug")
2761 << "Check |" << a
<< "| >= |" << b
<< "|" << std::endl
;
2762 unsigned pexp_size
= exp
.size();
2763 if (compareMonomial(oa
, a
, 0, a_exp_proc
, ob
, b
, 0, b_exp_proc
, 0, exp
, lem
,
2767 exp
.resize(pexp_size
);
2768 Trace("nl-ext-comp-debug")
2769 << "Check |" << b
<< "| >= |" << a
<< "|" << std::endl
;
2770 if (compareMonomial(ob
, b
, 0, b_exp_proc
, oa
, a
, 0, a_exp_proc
, 0, exp
, lem
,
2779 bool NonlinearExtension::cmp_holds(
2780 Node x
, Node y
, std::map
<Node
, std::map
<Node
, Node
> >& cmp_infers
,
2781 std::vector
<Node
>& exp
, std::map
<Node
, bool>& visited
) {
2784 } else if (visited
.find(x
) == visited
.end()) {
2786 std::map
<Node
, std::map
<Node
, Node
> >::iterator it
= cmp_infers
.find(x
);
2787 if (it
!= cmp_infers
.end()) {
2788 for (std::map
<Node
, Node
>::iterator itc
= it
->second
.begin();
2789 itc
!= it
->second
.end(); ++itc
) {
2790 exp
.push_back(itc
->second
);
2791 if (cmp_holds(itc
->first
, y
, cmp_infers
, exp
, visited
)) {
2801 // trying to show a ( >, = ) b by inequalities between variables in
2803 bool NonlinearExtension::compareMonomial(
2804 Node oa
, Node a
, unsigned a_index
, NodeMultiset
& a_exp_proc
, Node ob
,
2805 Node b
, unsigned b_index
, NodeMultiset
& b_exp_proc
, int status
,
2806 std::vector
<Node
>& exp
, std::vector
<Node
>& lem
,
2807 std::map
<int, std::map
<Node
, std::map
<Node
, Node
> > >& cmp_infers
) {
2808 Trace("nl-ext-comp-debug")
2809 << "compareMonomial " << oa
<< " and " << ob
<< ", indices = " << a_index
2810 << " " << b_index
<< std::endl
;
2811 Assert(status
== 0 || status
== 2);
2812 if (a_index
== d_m_vlist
[a
].size() && b_index
== d_m_vlist
[b
].size()) {
2813 // finished, compare abstract values
2814 int modelStatus
= compare(oa
, ob
, 3) * -2;
2815 Trace("nl-ext-comp") << "...finished comparison with " << oa
<< " <"
2816 << status
<< "> " << ob
2817 << ", model status = " << modelStatus
<< std::endl
;
2818 if (status
!= modelStatus
) {
2819 Trace("nl-ext-comp-infer")
2820 << "infer : " << oa
<< " <" << status
<< "> " << ob
<< std::endl
;
2822 // must state that all variables are non-zero
2823 for (unsigned j
= 0; j
< d_m_vlist
[a
].size(); j
++) {
2824 exp
.push_back(d_m_vlist
[a
][j
].eqNode(d_zero
).negate());
2827 Node clem
= NodeManager::currentNM()->mkNode(
2828 IMPLIES
, safeConstructNary(AND
, exp
), mkLit(oa
, ob
, status
, 1));
2829 Trace("nl-ext-comp-lemma") << "comparison lemma : " << clem
<< std::endl
;
2830 lem
.push_back(clem
);
2831 cmp_infers
[status
][oa
][ob
] = clem
;
2835 // get a/b variable information
2839 if (a_index
< d_m_vlist
[a
].size()) {
2840 av
= d_m_vlist
[a
][a_index
];
2841 Assert(a_exp_proc
[av
] <= d_m_exp
[a
][av
]);
2842 aexp
= d_m_exp
[a
][av
] - a_exp_proc
[av
];
2844 return compareMonomial(oa
, a
, a_index
+ 1, a_exp_proc
, ob
, b
, b_index
,
2845 b_exp_proc
, status
, exp
, lem
, cmp_infers
);
2847 Assert(d_order_vars
.find(av
) != d_order_vars
.end());
2848 avo
= d_order_vars
[av
];
2853 if (b_index
< d_m_vlist
[b
].size()) {
2854 bv
= d_m_vlist
[b
][b_index
];
2855 Assert(b_exp_proc
[bv
] <= d_m_exp
[b
][bv
]);
2856 bexp
= d_m_exp
[b
][bv
] - b_exp_proc
[bv
];
2858 return compareMonomial(oa
, a
, a_index
, a_exp_proc
, ob
, b
, b_index
+ 1,
2859 b_exp_proc
, status
, exp
, lem
, cmp_infers
);
2861 Assert(d_order_vars
.find(bv
) != d_order_vars
.end());
2862 bvo
= d_order_vars
[bv
];
2864 // get "one" information
2865 Assert(d_order_vars
.find(d_one
) != d_order_vars
.end());
2866 unsigned ovo
= d_order_vars
[d_one
];
2867 Trace("nl-ext-comp-debug") << "....vars : " << av
<< "^" << aexp
<< " "
2868 << bv
<< "^" << bexp
<< std::endl
;
2873 Trace("nl-ext-comp-debug") << "...take leading " << bv
<< std::endl
;
2874 // can multiply b by <=1
2875 exp
.push_back(mkLit(d_one
, bv
, bvo
== ovo
? 0 : 2, 1));
2876 return compareMonomial(oa
, a
, a_index
, a_exp_proc
, ob
, b
, b_index
+ 1,
2877 b_exp_proc
, bvo
== ovo
? status
: 2, exp
, lem
,
2880 Trace("nl-ext-comp-debug")
2881 << "...failure, unmatched |b|>1 component." << std::endl
;
2884 } else if (bv
.isNull()) {
2886 Trace("nl-ext-comp-debug") << "...take leading " << av
<< std::endl
;
2887 // can multiply a by >=1
2888 exp
.push_back(mkLit(av
, d_one
, avo
== ovo
? 0 : 2, 1));
2889 return compareMonomial(oa
, a
, a_index
+ 1, a_exp_proc
, ob
, b
, b_index
,
2890 b_exp_proc
, avo
== ovo
? status
: 2, exp
, lem
,
2893 Trace("nl-ext-comp-debug")
2894 << "...failure, unmatched |a|<1 component." << std::endl
;
2898 Assert(!av
.isNull() && !bv
.isNull());
2900 if (bvo
< ovo
&& avo
>= ovo
) {
2901 Trace("nl-ext-comp-debug") << "...take leading " << av
<< std::endl
;
2902 // do avo>=1 instead
2903 exp
.push_back(mkLit(av
, d_one
, avo
== ovo
? 0 : 2, 1));
2904 return compareMonomial(oa
, a
, a_index
+ 1, a_exp_proc
, ob
, b
, b_index
,
2905 b_exp_proc
, avo
== ovo
? status
: 2, exp
, lem
,
2908 unsigned min_exp
= aexp
> bexp
? bexp
: aexp
;
2909 a_exp_proc
[av
] += min_exp
;
2910 b_exp_proc
[bv
] += min_exp
;
2911 Trace("nl-ext-comp-debug")
2912 << "...take leading " << min_exp
<< " from " << av
<< " and "
2914 exp
.push_back(mkLit(av
, bv
, avo
== bvo
? 0 : 2, 1));
2915 bool ret
= compareMonomial(oa
, a
, a_index
, a_exp_proc
, ob
, b
, b_index
,
2916 b_exp_proc
, avo
== bvo
? status
: 2, exp
,
2918 a_exp_proc
[av
] -= min_exp
;
2919 b_exp_proc
[bv
] -= min_exp
;
2924 Trace("nl-ext-comp-debug") << "...take leading " << bv
<< std::endl
;
2925 // try multiply b <= 1
2926 exp
.push_back(mkLit(d_one
, bv
, bvo
== ovo
? 0 : 2, 1));
2927 return compareMonomial(oa
, a
, a_index
, a_exp_proc
, ob
, b
, b_index
+ 1,
2928 b_exp_proc
, bvo
== ovo
? status
: 2, exp
, lem
,
2931 Trace("nl-ext-comp-debug")
2932 << "...failure, leading |b|>|a|>1 component." << std::endl
;
2940 // status 0 : n equal, -1 : n superset, 1 : n subset
2941 void NonlinearExtension::MonomialIndex::addTerm(Node n
,
2942 const std::vector
<Node
>& reps
,
2943 NonlinearExtension
* nla
,
2944 int status
, unsigned argIndex
) {
2946 if (argIndex
== reps
.size()) {
2947 d_monos
.push_back(n
);
2949 d_data
[reps
[argIndex
]].addTerm(n
, reps
, nla
, status
, argIndex
+ 1);
2952 for (std::map
<Node
, MonomialIndex
>::iterator it
= d_data
.begin();
2953 it
!= d_data
.end(); ++it
) {
2954 if (status
!= 0 || argIndex
== reps
.size() || it
->first
!= reps
[argIndex
]) {
2955 // if we do not contain this variable, then if we were a superset,
2956 // fail (-2), otherwise we are subset. if we do contain this
2957 // variable, then if we were equal, we are superset since variables
2958 // are ordered, otherwise we remain the same.
2960 std::find(reps
.begin(), reps
.end(), it
->first
) == reps
.end()
2961 ? (status
>= 0 ? 1 : -2)
2962 : (status
== 0 ? -1 : status
);
2963 if (new_status
!= -2) {
2964 it
->second
.addTerm(n
, reps
, nla
, new_status
, argIndex
);
2968 // compare for subsets
2969 for (unsigned i
= 0; i
< d_monos
.size(); i
++) {
2970 Node m
= d_monos
[i
];
2972 // we are superset if we are equal and haven't traversed all variables
2973 int cstatus
= status
== 0 ? (argIndex
== reps
.size() ? 0 : -1) : status
;
2974 Trace("nl-ext-mindex-debug") << " compare " << n
<< " and " << m
2975 << ", status = " << cstatus
<< std::endl
;
2976 if (cstatus
<= 0 && nla
->isMonomialSubset(m
, n
)) {
2977 nla
->registerMonomialSubset(m
, n
);
2978 Trace("nl-ext-mindex-debug") << "...success" << std::endl
;
2979 } else if (cstatus
>= 0 && nla
->isMonomialSubset(n
, m
)) {
2980 nla
->registerMonomialSubset(n
, m
);
2981 Trace("nl-ext-mindex-debug") << "...success (rev)" << std::endl
;
2987 std::vector
<Node
> NonlinearExtension::checkMonomialSign() {
2988 std::vector
<Node
> lemmas
;
2989 std::map
<Node
, int> signs
;
2990 Trace("nl-ext") << "Get monomial sign lemmas..." << std::endl
;
2991 for (unsigned j
= 0; j
< d_ms
.size(); j
++) {
2993 if (d_ms_proc
.find(a
) == d_ms_proc
.end()) {
2994 std::vector
<Node
> exp
;
2995 Trace("nl-ext-debug") << " process " << a
<< ", mv=" << d_mv
[0][a
] << "..." << std::endl
;
2996 if( d_m_nconst_factor
.find( a
)==d_m_nconst_factor
.end() ){
2997 signs
[a
] = compareSign(a
, a
, 0, 1, exp
, lemmas
);
2998 if (signs
[a
] == 0) {
2999 d_ms_proc
[a
] = true;
3000 Trace("nl-ext-debug") << "...mark " << a
3001 << " reduced since its value is 0." << std::endl
;
3004 Trace("nl-ext-debug") << "...can't conclude sign lemma for " << a
<< " since model value of a factor is non-constant." << std::endl
;
3012 std::vector
<Node
> NonlinearExtension::checkMonomialMagnitude( unsigned c
) {
3014 std::vector
<Node
> lemmas
;
3015 // if (x,y,L) in cmp_infers, then x > y inferred as conclusion of L
3017 std::map
<int, std::map
<Node
, std::map
<Node
, Node
> > > cmp_infers
;
3018 Trace("nl-ext") << "Get monomial comparison lemmas (order=" << r
3019 << ", compare=" << c
<< ")..." << std::endl
;
3020 for (unsigned j
= 0; j
< d_ms
.size(); j
++) {
3022 if (d_ms_proc
.find(a
) == d_ms_proc
.end() &&
3023 d_m_nconst_factor
.find( a
)==d_m_nconst_factor
.end()) {
3025 // compare magnitude against 1
3026 std::vector
<Node
> exp
;
3027 NodeMultiset a_exp_proc
;
3028 NodeMultiset b_exp_proc
;
3029 compareMonomial(a
, a
, a_exp_proc
, d_one
, d_one
, b_exp_proc
, exp
,
3030 lemmas
, cmp_infers
);
3032 std::map
<Node
, NodeMultiset
>::iterator itmea
= d_m_exp
.find(a
);
3033 Assert(itmea
!= d_m_exp
.end());
3035 // TODO : not just against containing variables?
3036 // compare magnitude against variables
3037 for (unsigned k
= 0; k
< d_ms_vars
.size(); k
++) {
3038 Node v
= d_ms_vars
[k
];
3039 Assert( d_mv
[0].find( v
)!=d_mv
[0].end() );
3040 if( d_mv
[0][v
].isConst() ){
3041 std::vector
<Node
> exp
;
3042 NodeMultiset a_exp_proc
;
3043 NodeMultiset b_exp_proc
;
3044 if (itmea
->second
.find(v
) != itmea
->second
.end()) {
3047 setMonomialFactor(a
, v
, a_exp_proc
);
3048 setMonomialFactor(v
, a
, b_exp_proc
);
3049 compareMonomial(a
, a
, a_exp_proc
, v
, v
, b_exp_proc
, exp
,
3050 lemmas
, cmp_infers
);
3055 // compare magnitude against other non-linear monomials
3056 for (unsigned k
= (j
+ 1); k
< d_ms
.size(); k
++) {
3058 //(signs[a]==signs[b])==(r==0)
3059 if (d_ms_proc
.find(b
) == d_ms_proc
.end() &&
3060 d_m_nconst_factor
.find( b
)==d_m_nconst_factor
.end()) {
3061 std::map
<Node
, NodeMultiset
>::iterator itmeb
=
3063 Assert(itmeb
!= d_m_exp
.end());
3065 std::vector
<Node
> exp
;
3066 // take common factors of monomials, set minimum of
3067 // common exponents as processed
3068 NodeMultiset a_exp_proc
;
3069 NodeMultiset b_exp_proc
;
3070 for (NodeMultiset::iterator itmea2
= itmea
->second
.begin();
3071 itmea2
!= itmea
->second
.end(); ++itmea2
) {
3072 NodeMultiset::iterator itmeb2
=
3073 itmeb
->second
.find(itmea2
->first
);
3074 if (itmeb2
!= itmeb
->second
.end()) {
3075 unsigned min_exp
= itmea2
->second
> itmeb2
->second
3078 a_exp_proc
[itmea2
->first
] = min_exp
;
3079 b_exp_proc
[itmea2
->first
] = min_exp
;
3080 Trace("nl-ext-comp")
3081 << "Common exponent : " << itmea2
->first
<< " : "
3082 << min_exp
<< std::endl
;
3085 if (!a_exp_proc
.empty()) {
3086 setMonomialFactor(a
, b
, a_exp_proc
);
3087 setMonomialFactor(b
, a
, b_exp_proc
);
3090 if( !a_exp_proc.empty() ){
3091 //reduction based on common exponents a > 0 => ( a * b
3092 <> a * c <=> b <> c ), a < 0 => ( a * b <> a * c <=> b
3093 !<> c ) ? }else{ compareMonomial( a, a, a_exp_proc, b,
3094 b, b_exp_proc, exp, lemmas );
3097 compareMonomial(a
, a
, a_exp_proc
, b
, b
, b_exp_proc
, exp
,
3098 lemmas
, cmp_infers
);
3105 // remove redundant lemmas, e.g. if a > b, b > c, a > c were
3106 // inferred, discard lemma with conclusion a > c
3107 Trace("nl-ext-comp") << "Compute redundancies for " << lemmas
.size()
3108 << " lemmas." << std::endl
;
3110 std::vector
<Node
> r_lemmas
;
3111 for (std::map
<int, std::map
<Node
, std::map
<Node
, Node
> > >::iterator itb
=
3113 itb
!= cmp_infers
.end(); ++itb
) {
3114 for (std::map
<Node
, std::map
<Node
, Node
> >::iterator itc
=
3115 itb
->second
.begin();
3116 itc
!= itb
->second
.end(); ++itc
) {
3117 for (std::map
<Node
, Node
>::iterator itc2
= itc
->second
.begin();
3118 itc2
!= itc
->second
.end(); ++itc2
) {
3119 std::map
<Node
, bool> visited
;
3120 for (std::map
<Node
, Node
>::iterator itc3
= itc
->second
.begin();
3121 itc3
!= itc
->second
.end(); ++itc3
) {
3122 if (itc3
->first
!= itc2
->first
) {
3123 std::vector
<Node
> exp
;
3124 if (cmp_holds(itc3
->first
, itc2
->first
, itb
->second
, exp
,
3126 r_lemmas
.push_back(itc2
->second
);
3127 Trace("nl-ext-comp")
3128 << "...inference of " << itc
->first
<< " > "
3129 << itc2
->first
<< " was redundant." << std::endl
;
3137 std::vector
<Node
> nr_lemmas
;
3138 for (unsigned i
= 0; i
< lemmas
.size(); i
++) {
3139 if (std::find(r_lemmas
.begin(), r_lemmas
.end(), lemmas
[i
]) ==
3141 nr_lemmas
.push_back(lemmas
[i
]);
3144 // TODO: only take maximal lower/minimial lower bounds?
3146 Trace("nl-ext-comp") << nr_lemmas
.size() << " / " << lemmas
.size()
3147 << " were non-redundant." << std::endl
;
3151 std::vector
<Node
> NonlinearExtension::checkTangentPlanes() {
3152 std::vector
< Node
> lemmas
;
3153 Trace("nl-ext") << "Get monomial tangent plane lemmas..." << std::endl
;
3154 unsigned kstart
= d_ms_vars
.size();
3155 for (unsigned k
= kstart
; k
< d_mterms
.size(); k
++) {
3156 Node t
= d_mterms
[k
];
3157 // if this term requires a refinement
3158 if (d_tplane_refine
.find(t
) != d_tplane_refine
.end())
3160 Trace("nl-ext-tplanes")
3161 << "Look at monomial requiring refinement : " << t
<< std::endl
;
3162 // get a decomposition
3163 std::map
<Node
, std::vector
<Node
> >::iterator it
=
3164 d_m_contain_children
.find(t
);
3165 if (it
!= d_m_contain_children
.end()) {
3166 std::map
<Node
, std::map
<Node
, bool> > dproc
;
3167 for (unsigned j
= 0; j
< it
->second
.size(); j
++) {
3168 Node tc
= it
->second
[j
];
3170 Node tc_diff
= d_m_contain_umult
[tc
][t
];
3171 Assert(!tc_diff
.isNull());
3172 Node a
= tc
< tc_diff
? tc
: tc_diff
;
3173 Node b
= tc
< tc_diff
? tc_diff
: tc
;
3174 if (dproc
[a
].find(b
) == dproc
[a
].end()) {
3176 Trace("nl-ext-tplanes")
3177 << " decomposable into : " << a
<< " * " << b
<< std::endl
;
3178 Node a_v_c
= computeModelValue(a
, 1);
3179 Node b_v_c
= computeModelValue(b
, 1);
3180 // points we will add tangent planes for
3181 std::vector
< Node
> pts
[2];
3182 pts
[0].push_back( a_v_c
);
3183 pts
[1].push_back( b_v_c
);
3184 // if previously refined
3185 bool prevRefine
= d_tangent_val_bound
[0][a
].find( b
)!=d_tangent_val_bound
[0][a
].end();
3186 // a_min, a_max, b_min, b_max
3187 for( unsigned p
=0; p
<4; p
++ ){
3188 Node curr_v
= p
<=1 ? a_v_c
: b_v_c
;
3190 Node pt_v
= d_tangent_val_bound
[p
][a
][b
];
3191 Assert( !pt_v
.isNull() );
3193 Node do_extend
= NodeManager::currentNM()->mkNode( ( p
==1 || p
==3 ) ? GT
: LT
, curr_v
, pt_v
);
3194 do_extend
= Rewriter::rewrite( do_extend
);
3195 if( do_extend
==d_true
){
3196 for( unsigned q
=0; q
<2; q
++ ){
3197 pts
[ p
<=1 ? 0 : 1 ].push_back( curr_v
);
3198 pts
[ p
<=1 ? 1 : 0 ].push_back( d_tangent_val_bound
[ p
<=1 ? 2+q
: q
][a
][b
] );
3203 d_tangent_val_bound
[p
][a
][b
] = curr_v
;
3207 for( unsigned p
=0; p
<pts
[0].size(); p
++ ){
3208 Node a_v
= pts
[0][p
];
3209 Node b_v
= pts
[1][p
];
3212 Node tplane
= NodeManager::currentNM()->mkNode(
3214 NodeManager::currentNM()->mkNode(
3216 NodeManager::currentNM()->mkNode(MULT
, b_v
, a
),
3217 NodeManager::currentNM()->mkNode(MULT
, a_v
, b
)),
3218 NodeManager::currentNM()->mkNode(MULT
, a_v
, b_v
));
3219 for (unsigned d
= 0; d
< 4; d
++) {
3220 Node aa
= NodeManager::currentNM()->mkNode(
3221 d
== 0 || d
== 3 ? GEQ
: LEQ
, a
, a_v
);
3222 Node ab
= NodeManager::currentNM()->mkNode(
3223 d
== 1 || d
== 3 ? GEQ
: LEQ
, b
, b_v
);
3224 Node conc
= NodeManager::currentNM()->mkNode(
3225 d
<= 1 ? LEQ
: GEQ
, t
, tplane
);
3226 Node tlem
= NodeManager::currentNM()->mkNode(
3227 OR
, aa
.negate(), ab
.negate(), conc
);
3228 Trace("nl-ext-tplanes")
3229 << "Tangent plane lemma : " << tlem
<< std::endl
;
3230 lemmas
.push_back(tlem
);
3239 Trace("nl-ext") << "...trying " << lemmas
.size()
3240 << " tangent plane lemmas..." << std::endl
;
3244 std::vector
<Node
> NonlinearExtension::checkMonomialInferBounds(
3245 std::vector
<Node
>& nt_lemmas
,
3246 const std::vector
<Node
>& asserts
,
3247 const std::vector
<Node
>& false_asserts
)
3249 std::vector
< Node
> lemmas
;
3250 // register constraints
3251 Trace("nl-ext-debug") << "Register bound constraints..." << std::endl
;
3252 for (const Node
& lit
: asserts
)
3254 bool polarity
= lit
.getKind() != NOT
;
3255 Node atom
= lit
.getKind() == NOT
? lit
[0] : lit
;
3256 registerConstraint(atom
);
3258 std::find(false_asserts
.begin(), false_asserts
.end(), lit
)
3259 != false_asserts
.end();
3260 // add information about bounds to variables
3261 std::map
<Node
, std::map
<Node
, ConstraintInfo
> >::iterator itc
=
3262 d_c_info
.find(atom
);
3263 std::map
<Node
, std::map
<Node
, bool> >::iterator itcm
=
3264 d_c_info_maxm
.find(atom
);
3265 if (itc
!= d_c_info
.end()) {
3266 Assert(itcm
!= d_c_info_maxm
.end());
3267 for (std::map
<Node
, ConstraintInfo
>::iterator itcc
= itc
->second
.begin();
3268 itcc
!= itc
->second
.end(); ++itcc
) {
3269 Node x
= itcc
->first
;
3270 Node coeff
= itcc
->second
.d_coeff
;
3271 Node rhs
= itcc
->second
.d_rhs
;
3272 Kind type
= itcc
->second
.d_type
;
3278 // we will take the strict inequality in the direction of the
3280 Node lhs
= ArithMSum::mkCoeffTerm(coeff
, x
);
3281 Node query
= NodeManager::currentNM()->mkNode(GT
, lhs
, rhs
);
3282 Node query_mv
= computeModelValue(query
, 1);
3283 if (query_mv
== d_true
) {
3287 Assert(query_mv
== d_false
);
3288 exp
= NodeManager::currentNM()->mkNode(LT
, lhs
, rhs
);
3292 type
= negateKind(type
);
3295 // add to status if maximal degree
3296 d_ci_max
[x
][coeff
][rhs
] = itcm
->second
.find(x
) != itcm
->second
.end();
3297 if (Trace
.isOn("nl-ext-bound-debug2")) {
3298 Node t
= ArithMSum::mkCoeffTerm(coeff
, x
);
3299 Trace("nl-ext-bound-debug2")
3300 << "Add Bound: " << t
<< " " << type
<< " " << rhs
<< " by "
3301 << exp
<< std::endl
;
3303 bool updated
= true;
3304 std::map
<Node
, Kind
>::iterator its
= d_ci
[x
][coeff
].find(rhs
);
3305 if (its
== d_ci
[x
][coeff
].end()) {
3306 d_ci
[x
][coeff
][rhs
] = type
;
3307 d_ci_exp
[x
][coeff
][rhs
] = exp
;
3308 } else if (type
!= its
->second
) {
3309 Trace("nl-ext-bound-debug2")
3310 << "Joining kinds : " << type
<< " " << its
->second
<< std::endl
;
3311 Kind jk
= joinKinds(type
, its
->second
);
3312 if (jk
== UNDEFINED_KIND
)
3315 } else if (jk
!= its
->second
) {
3317 d_ci
[x
][coeff
][rhs
] = type
;
3318 d_ci_exp
[x
][coeff
][rhs
] = exp
;
3320 d_ci
[x
][coeff
][rhs
] = jk
;
3321 d_ci_exp
[x
][coeff
][rhs
] = NodeManager::currentNM()->mkNode(
3322 AND
, d_ci_exp
[x
][coeff
][rhs
], exp
);
3328 if (Trace
.isOn("nl-ext-bound")) {
3330 Trace("nl-ext-bound") << "Bound: ";
3331 debugPrintBound("nl-ext-bound", coeff
, x
, d_ci
[x
][coeff
][rhs
], rhs
);
3332 Trace("nl-ext-bound") << " by " << d_ci_exp
[x
][coeff
][rhs
];
3333 if (d_ci_max
[x
][coeff
][rhs
]) {
3334 Trace("nl-ext-bound") << ", is max degree";
3336 Trace("nl-ext-bound") << std::endl
;
3339 // compute if bound is not satisfied, and store what is required
3340 // for a possible refinement
3341 if (options::nlExtTangentPlanes()) {
3343 d_tplane_refine
.insert(x
);
3349 // reflexive constraints
3351 for (unsigned j
= 0; j
< d_mterms
.size(); j
++) {
3352 Node n
= d_mterms
[j
];
3353 d_ci
[n
][null_coeff
][n
] = EQUAL
;
3354 d_ci_exp
[n
][null_coeff
][n
] = d_true
;
3355 d_ci_max
[n
][null_coeff
][n
] = false;
3358 Trace("nl-ext") << "Get inferred bound lemmas..." << std::endl
;
3360 for (unsigned k
= 0; k
< d_mterms
.size(); k
++) {
3361 Node x
= d_mterms
[k
];
3362 Trace("nl-ext-bound-debug")
3363 << "Process bounds for " << x
<< " : " << std::endl
;
3364 std::map
<Node
, std::vector
<Node
> >::iterator itm
=
3365 d_m_contain_parent
.find(x
);
3366 if (itm
!= d_m_contain_parent
.end()) {
3367 Trace("nl-ext-bound-debug") << "...has " << itm
->second
.size()
3368 << " parent monomials." << std::endl
;
3369 // check derived bounds
3370 std::map
<Node
, std::map
<Node
, std::map
<Node
, Kind
> > >::iterator itc
=
3372 if (itc
!= d_ci
.end()) {
3373 for (std::map
<Node
, std::map
<Node
, Kind
> >::iterator itcc
=
3374 itc
->second
.begin();
3375 itcc
!= itc
->second
.end(); ++itcc
) {
3376 Node coeff
= itcc
->first
;
3377 Node t
= ArithMSum::mkCoeffTerm(coeff
, x
);
3378 for (std::map
<Node
, Kind
>::iterator itcr
= itcc
->second
.begin();
3379 itcr
!= itcc
->second
.end(); ++itcr
) {
3380 Node rhs
= itcr
->first
;
3381 // only consider this bound if maximal degree
3382 if (d_ci_max
[x
][coeff
][rhs
]) {
3383 Kind type
= itcr
->second
;
3384 for (unsigned j
= 0; j
< itm
->second
.size(); j
++) {
3385 Node y
= itm
->second
[j
];
3386 Assert(d_m_contain_mult
[x
].find(y
) !=
3387 d_m_contain_mult
[x
].end());
3388 Node mult
= d_m_contain_mult
[x
][y
];
3389 // x <k> t => m*x <k'> t where y = m*x
3390 // get the sign of mult
3391 Node mmv
= computeModelValue(mult
);
3392 Trace("nl-ext-bound-debug2")
3393 << "Model value of " << mult
<< " is " << mmv
<< std::endl
;
3395 int mmv_sign
= mmv
.getConst
<Rational
>().sgn();
3396 Trace("nl-ext-bound-debug2")
3397 << " sign of " << mmv
<< " is " << mmv_sign
<< std::endl
;
3398 if (mmv_sign
!= 0) {
3399 Trace("nl-ext-bound-debug")
3400 << " from " << x
<< " * " << mult
<< " = " << y
3401 << " and " << t
<< " " << type
<< " " << rhs
3402 << ", infer : " << std::endl
;
3404 mmv_sign
== -1 ? reverseRelationKind(type
) : type
;
3406 NodeManager::currentNM()->mkNode(MULT
, mult
, t
);
3408 NodeManager::currentNM()->mkNode(MULT
, mult
, rhs
);
3409 Node infer
= NodeManager::currentNM()->mkNode(
3410 infer_type
, infer_lhs
, infer_rhs
);
3411 Trace("nl-ext-bound-debug") << " " << infer
<< std::endl
;
3412 infer
= Rewriter::rewrite(infer
);
3413 Trace("nl-ext-bound-debug2")
3414 << " ...rewritten : " << infer
<< std::endl
;
3415 // check whether it is false in model for abstraction
3416 Node infer_mv
= computeModelValue(infer
, 1);
3417 Trace("nl-ext-bound-debug")
3418 << " ...infer model value is " << infer_mv
3420 if (infer_mv
== d_false
) {
3421 Node exp
= NodeManager::currentNM()->mkNode(
3423 NodeManager::currentNM()->mkNode(
3424 mmv_sign
== 1 ? GT
: LT
, mult
, d_zero
),
3425 d_ci_exp
[x
][coeff
][rhs
]);
3427 NodeManager::currentNM()->mkNode(IMPLIES
, exp
, infer
);
3428 Node pr_iblem
= iblem
;
3429 iblem
= Rewriter::rewrite(iblem
);
3430 bool introNewTerms
= hasNewMonomials(iblem
, d_ms
);
3431 Trace("nl-ext-bound-lemma")
3432 << "*** Bound inference lemma : " << iblem
3433 << " (pre-rewrite : " << pr_iblem
<< ")" << std::endl
;
3434 // Trace("nl-ext-bound-lemma") << " intro new
3435 // monomials = " << introNewTerms << std::endl;
3436 if (!introNewTerms
) {
3437 lemmas
.push_back(iblem
);
3439 nt_lemmas
.push_back(iblem
);
3443 Trace("nl-ext-bound-debug") << " ...coefficient " << mult
3444 << " is zero." << std::endl
;
3447 Trace("nl-ext-bound-debug") << " ...coefficient " << mult
3448 << " is non-constant (probably transcendental)." << std::endl
;
3456 Trace("nl-ext-bound-debug") << "...has no parent monomials." << std::endl
;
3462 std::vector
<Node
> NonlinearExtension::checkFactoring(
3463 const std::vector
<Node
>& asserts
, const std::vector
<Node
>& false_asserts
)
3465 std::vector
< Node
> lemmas
;
3466 Trace("nl-ext") << "Get factoring lemmas..." << std::endl
;
3467 for (const Node
& lit
: asserts
)
3469 bool polarity
= lit
.getKind() != NOT
;
3470 Node atom
= lit
.getKind() == NOT
? lit
[0] : lit
;
3471 Node litv
= computeModelValue(lit
);
3472 bool considerLit
= false;
3473 if( d_skolem_atoms
.find(atom
) != d_skolem_atoms
.end() )
3475 //always consider skolem literals
3480 // Only consider literals that are in false_asserts.
3481 considerLit
= std::find(false_asserts
.begin(), false_asserts
.end(), lit
)
3482 != false_asserts
.end();
3487 std::map
<Node
, Node
> msum
;
3488 if (ArithMSum::getMonomialSumLit(atom
, msum
))
3490 Trace("nl-ext-factor") << "Factoring for literal " << lit
<< ", monomial sum is : " << std::endl
;
3491 if (Trace
.isOn("nl-ext-factor")) {
3492 ArithMSum::debugPrintMonomialSum(msum
, "nl-ext-factor");
3494 std::map
< Node
, std::vector
< Node
> > factor_to_mono
;
3495 std::map
< Node
, std::vector
< Node
> > factor_to_mono_orig
;
3496 for( std::map
<Node
, Node
>::iterator itm
= msum
.begin(); itm
!= msum
.end(); ++itm
){
3497 if( !itm
->first
.isNull() ){
3498 if( itm
->first
.getKind()==NONLINEAR_MULT
){
3499 std::vector
< Node
> children
;
3500 for( unsigned i
=0; i
<itm
->first
.getNumChildren(); i
++ ){
3501 children
.push_back( itm
->first
[i
] );
3503 std::map
< Node
, bool > processed
;
3504 for( unsigned i
=0; i
<itm
->first
.getNumChildren(); i
++ ){
3505 if( processed
.find( itm
->first
[i
] )==processed
.end() ){
3506 processed
[itm
->first
[i
]] = true;
3507 children
[i
] = d_one
;
3508 if( !itm
->second
.isNull() ){
3509 children
.push_back( itm
->second
);
3511 Node val
= NodeManager::currentNM()->mkNode(MULT
, children
);
3512 if( !itm
->second
.isNull() ){
3513 children
.pop_back();
3515 children
[i
] = itm
->first
[i
];
3516 val
= Rewriter::rewrite( val
);
3517 factor_to_mono
[itm
->first
[i
]].push_back( val
);
3518 factor_to_mono_orig
[itm
->first
[i
]].push_back( itm
->first
);
3524 for( std::map
< Node
, std::vector
< Node
> >::iterator itf
= factor_to_mono
.begin(); itf
!= factor_to_mono
.end(); ++itf
){
3525 Node x
= itf
->first
;
3526 if (itf
->second
.size() == 1)
3528 std::map
<Node
, Node
>::iterator itm
= msum
.find(x
);
3529 if (itm
!= msum
.end())
3531 itf
->second
.push_back(itm
->second
.isNull() ? d_one
: itm
->second
);
3532 factor_to_mono_orig
[x
].push_back(x
);
3535 if( itf
->second
.size()>1 ){
3536 Node sum
= NodeManager::currentNM()->mkNode(PLUS
, itf
->second
);
3537 sum
= Rewriter::rewrite( sum
);
3538 Trace("nl-ext-factor")
3539 << "* Factored sum for " << x
<< " : " << sum
<< std::endl
;
3540 Node kf
= getFactorSkolem( sum
, lemmas
);
3541 std::vector
< Node
> poly
;
3542 poly
.push_back(NodeManager::currentNM()->mkNode(MULT
, x
, kf
));
3543 std::map
<Node
, std::vector
<Node
> >::iterator itfo
=
3544 factor_to_mono_orig
.find(x
);
3545 Assert( itfo
!=factor_to_mono_orig
.end() );
3546 for( std::map
<Node
, Node
>::iterator itm
= msum
.begin(); itm
!= msum
.end(); ++itm
){
3547 if( std::find( itfo
->second
.begin(), itfo
->second
.end(), itm
->first
)==itfo
->second
.end() ){
3548 poly
.push_back(ArithMSum::mkCoeffTerm(
3549 itm
->second
, itm
->first
.isNull() ? d_one
: itm
->first
));
3552 Node polyn
= poly
.size() == 1
3554 : NodeManager::currentNM()->mkNode(PLUS
, poly
);
3555 Trace("nl-ext-factor") << "...factored polynomial : " << polyn
<< std::endl
;
3556 Node conc_lit
= NodeManager::currentNM()->mkNode( atom
.getKind(), polyn
, d_zero
);
3557 conc_lit
= Rewriter::rewrite( conc_lit
);
3558 d_skolem_atoms
.insert( conc_lit
);
3560 conc_lit
= conc_lit
.negate();
3563 std::vector
< Node
> lemma_disj
;
3564 lemma_disj
.push_back( lit
.negate() );
3565 lemma_disj
.push_back( conc_lit
);
3566 Node flem
= NodeManager::currentNM()->mkNode(OR
, lemma_disj
);
3567 Trace("nl-ext-factor") << "...lemma is " << flem
<< std::endl
;
3568 lemmas
.push_back( flem
);
3577 Node
NonlinearExtension::getFactorSkolem( Node n
, std::vector
< Node
>& lemmas
) {
3578 std::map
< Node
, Node
>::iterator itf
= d_factor_skolem
.find( n
);
3579 if( itf
==d_factor_skolem
.end() ){
3580 Node k
= NodeManager::currentNM()->mkSkolem( "kf", n
.getType() );
3581 Node k_eq
= Rewriter::rewrite( k
.eqNode( n
) );
3582 d_skolem_atoms
.insert( k_eq
);
3583 lemmas
.push_back( k_eq
);
3584 d_factor_skolem
[n
] = k
;
3591 std::vector
<Node
> NonlinearExtension::checkMonomialInferResBounds() {
3592 std::vector
< Node
> lemmas
;
3593 Trace("nl-ext") << "Get monomial resolution inferred bound lemmas..." << std::endl
;
3594 for (unsigned j
= 0; j
< d_mterms
.size(); j
++) {
3595 Node a
= d_mterms
[j
];
3596 std::map
<Node
, std::map
<Node
, std::map
<Node
, Kind
> > >::iterator itca
=
3598 if (itca
!= d_ci
.end()) {
3599 for (unsigned k
= (j
+ 1); k
< d_mterms
.size(); k
++) {
3600 Node b
= d_mterms
[k
];
3601 std::map
<Node
, std::map
<Node
, std::map
<Node
, Kind
> > >::iterator
3602 itcb
= d_ci
.find(b
);
3603 if (itcb
!= d_ci
.end()) {
3604 Trace("nl-ext-rbound-debug") << "resolution inferences : compare "
3605 << a
<< " and " << b
<< std::endl
;
3606 // if they have common factors
3607 std::map
<Node
, Node
>::iterator ita
= d_mono_diff
[a
].find(b
);
3608 if (ita
!= d_mono_diff
[a
].end()) {
3609 std::map
<Node
, Node
>::iterator itb
= d_mono_diff
[b
].find(a
);
3610 Assert(itb
!= d_mono_diff
[b
].end());
3611 Node mv_a
= computeModelValue(ita
->second
, 1);
3612 Assert(mv_a
.isConst());
3613 int mv_a_sgn
= mv_a
.getConst
<Rational
>().sgn();
3614 Assert(mv_a_sgn
!= 0);
3615 Node mv_b
= computeModelValue(itb
->second
, 1);
3616 Assert(mv_b
.isConst());
3617 int mv_b_sgn
= mv_b
.getConst
<Rational
>().sgn();
3618 Assert(mv_b_sgn
!= 0);
3619 Trace("nl-ext-rbound") << "Get resolution inferences for [a] "
3620 << a
<< " vs [b] " << b
<< std::endl
;
3621 Trace("nl-ext-rbound")
3622 << " [a] factor is " << ita
->second
3623 << ", sign in model = " << mv_a_sgn
<< std::endl
;
3624 Trace("nl-ext-rbound")
3625 << " [b] factor is " << itb
->second
3626 << ", sign in model = " << mv_b_sgn
<< std::endl
;
3628 std::vector
<Node
> exp
;
3630 for (std::map
<Node
, std::map
<Node
, Kind
> >::iterator itcac
=
3631 itca
->second
.begin();
3632 itcac
!= itca
->second
.end(); ++itcac
) {
3633 Node coeff_a
= itcac
->first
;
3634 for (std::map
<Node
, Kind
>::iterator itcar
=
3635 itcac
->second
.begin();
3636 itcar
!= itcac
->second
.end(); ++itcar
) {
3637 Node rhs_a
= itcar
->first
;
3638 Node rhs_a_res_base
=
3639 NodeManager::currentNM()->mkNode(MULT
, itb
->second
, rhs_a
);
3640 rhs_a_res_base
= Rewriter::rewrite(rhs_a_res_base
);
3641 if (!hasNewMonomials(rhs_a_res_base
, d_ms
)) {
3642 Kind type_a
= itcar
->second
;
3643 exp
.push_back(d_ci_exp
[a
][coeff_a
][rhs_a
]);
3646 for (std::map
<Node
, std::map
<Node
, Kind
> >::iterator itcbc
=
3647 itcb
->second
.begin();
3648 itcbc
!= itcb
->second
.end(); ++itcbc
) {
3649 Node coeff_b
= itcbc
->first
;
3651 ArithMSum::mkCoeffTerm(coeff_b
, rhs_a_res_base
);
3652 for (std::map
<Node
, Kind
>::iterator itcbr
=
3653 itcbc
->second
.begin();
3654 itcbr
!= itcbc
->second
.end(); ++itcbr
) {
3655 Node rhs_b
= itcbr
->first
;
3656 Node rhs_b_res
= NodeManager::currentNM()->mkNode(
3657 MULT
, ita
->second
, rhs_b
);
3658 rhs_b_res
= ArithMSum::mkCoeffTerm(coeff_a
, rhs_b_res
);
3659 rhs_b_res
= Rewriter::rewrite(rhs_b_res
);
3660 if (!hasNewMonomials(rhs_b_res
, d_ms
)) {
3661 Kind type_b
= itcbr
->second
;
3662 exp
.push_back(d_ci_exp
[b
][coeff_b
][rhs_b
]);
3663 if (Trace
.isOn("nl-ext-rbound")) {
3664 Trace("nl-ext-rbound") << "* try bounds : ";
3665 debugPrintBound("nl-ext-rbound", coeff_a
, a
, type_a
,
3667 Trace("nl-ext-rbound") << std::endl
;
3668 Trace("nl-ext-rbound") << " ";
3669 debugPrintBound("nl-ext-rbound", coeff_b
, b
, type_b
,
3671 Trace("nl-ext-rbound") << std::endl
;
3674 for (unsigned r
= 0; r
< 2; r
++) {
3676 r
== 0 ? itb
->second
: ita
->second
;
3677 int pivot_factor_sign
=
3678 r
== 0 ? mv_b_sgn
: mv_a_sgn
;
3679 types
[r
] = r
== 0 ? type_a
: type_b
;
3680 if (pivot_factor_sign
== (r
== 0 ? 1 : -1)) {
3681 types
[r
] = reverseRelationKind(types
[r
]);
3683 if (pivot_factor_sign
== 1) {
3684 exp
.push_back(NodeManager::currentNM()->mkNode(
3685 GT
, pivot_factor
, d_zero
));
3687 exp
.push_back(NodeManager::currentNM()->mkNode(
3688 LT
, pivot_factor
, d_zero
));
3691 Kind jk
= transKinds(types
[0], types
[1]);
3692 Trace("nl-ext-rbound-debug")
3693 << "trans kind : " << types
[0] << " + "
3694 << types
[1] << " = " << jk
<< std::endl
;
3695 if (jk
!= UNDEFINED_KIND
)
3697 Node conc
= NodeManager::currentNM()->mkNode(
3698 jk
, rhs_a_res
, rhs_b_res
);
3699 Node conc_mv
= computeModelValue(conc
, 1);
3700 if (conc_mv
== d_false
) {
3701 Node rblem
= NodeManager::currentNM()->mkNode(
3703 NodeManager::currentNM()->mkNode(AND
, exp
),
3705 Trace("nl-ext-rbound-lemma-debug")
3706 << "Resolution bound lemma "
3709 << rblem
<< std::endl
;
3710 rblem
= Rewriter::rewrite(rblem
);
3711 Trace("nl-ext-rbound-lemma")
3712 << "Resolution bound lemma : " << rblem
3714 lemmas
.push_back(rblem
);
3735 std::vector
<Node
> NonlinearExtension::checkTranscendentalInitialRefine() {
3736 std::vector
< Node
> lemmas
;
3737 Trace("nl-ext") << "Get initial refinement lemmas for transcendental functions..." << std::endl
;
3738 for( std::map
< Kind
, std::map
< Node
, Node
> >::iterator it
= d_tf_rep_map
.begin(); it
!= d_tf_rep_map
.end(); ++it
){
3739 for( std::map
< Node
, Node
>::iterator itt
= it
->second
.begin(); itt
!= it
->second
.end(); ++itt
){
3740 Node t
= itt
->second
;
3741 Assert( d_mv
[1].find( t
)!=d_mv
[1].end() );
3742 //initial refinements
3743 if( d_tf_initial_refine
.find( t
)==d_tf_initial_refine
.end() ){
3744 d_tf_initial_refine
[t
] = true;
3746 if (it
->first
== SINE
)
3748 Node symn
= NodeManager::currentNM()->mkNode(
3749 SINE
, NodeManager::currentNM()->mkNode(MULT
, d_neg_one
, t
[0]));
3750 symn
= Rewriter::rewrite( symn
);
3751 //can assume its basis since phase is split over 0
3752 d_trig_is_base
[symn
] = true;
3753 Assert( d_trig_is_base
.find( t
) != d_trig_is_base
.end() );
3754 std::vector
< Node
> children
;
3756 lem
= NodeManager::currentNM()->mkNode(
3759 NodeManager::currentNM()->mkNode(
3761 NodeManager::currentNM()->mkNode(LEQ
, t
, d_one
),
3762 NodeManager::currentNM()->mkNode(GEQ
, t
, d_neg_one
)),
3764 NodeManager::currentNM()->mkNode(PLUS
, t
, symn
).eqNode(d_zero
),
3766 NodeManager::currentNM()->mkNode(
3768 NodeManager::currentNM()->mkNode(LT
, t
[0], d_zero
),
3769 NodeManager::currentNM()->mkNode(LT
, t
, d_zero
)),
3771 NodeManager::currentNM()->mkNode(
3773 NodeManager::currentNM()->mkNode(GT
, t
[0], d_zero
),
3774 NodeManager::currentNM()->mkNode(GT
, t
, d_zero
)));
3775 lem
= NodeManager::currentNM()->mkNode(
3779 NodeManager::currentNM()->mkNode(
3781 NodeManager::currentNM()->mkNode(
3783 NodeManager::currentNM()->mkNode(GT
, t
[0], d_zero
),
3784 NodeManager::currentNM()->mkNode(LT
, t
, t
[0])),
3785 NodeManager::currentNM()->mkNode(
3787 NodeManager::currentNM()->mkNode(LT
, t
[0], d_zero
),
3788 NodeManager::currentNM()->mkNode(GT
, t
, t
[0]))),
3790 NodeManager::currentNM()->mkNode(
3792 NodeManager::currentNM()->mkNode(
3794 NodeManager::currentNM()->mkNode(LT
, t
[0], d_pi
),
3795 NodeManager::currentNM()->mkNode(
3798 NodeManager::currentNM()->mkNode(MINUS
, d_pi
, t
[0]))),
3799 NodeManager::currentNM()->mkNode(
3801 NodeManager::currentNM()->mkNode(GT
, t
[0], d_pi_neg
),
3802 NodeManager::currentNM()->mkNode(
3805 NodeManager::currentNM()->mkNode(
3806 MINUS
, d_pi_neg
, t
[0])))));
3808 else if (it
->first
== EXPONENTIAL
)
3810 // ( exp(x) > 0 ) ^ ( x=0 <=> exp( x ) = 1 ) ^ ( x < 0 <=> exp( x ) <
3811 // 1 ) ^ ( x <= 0 V exp( x ) > x + 1 )
3812 lem
= NodeManager::currentNM()->mkNode(
3814 NodeManager::currentNM()->mkNode(GT
, t
, d_zero
),
3815 NodeManager::currentNM()->mkNode(
3816 EQUAL
, t
[0].eqNode(d_zero
), t
.eqNode(d_one
)),
3817 NodeManager::currentNM()->mkNode(
3819 NodeManager::currentNM()->mkNode(LT
, t
[0], d_zero
),
3820 NodeManager::currentNM()->mkNode(LT
, t
, d_one
)),
3821 NodeManager::currentNM()->mkNode(
3823 NodeManager::currentNM()->mkNode(LEQ
, t
[0], d_zero
),
3824 NodeManager::currentNM()->mkNode(
3827 NodeManager::currentNM()->mkNode(PLUS
, t
[0], d_one
))));
3829 if( !lem
.isNull() ){
3830 lemmas
.push_back( lem
);
3839 std::vector
<Node
> NonlinearExtension::checkTranscendentalMonotonic() {
3840 std::vector
< Node
> lemmas
;
3841 Trace("nl-ext") << "Get monotonicity lemmas for transcendental functions..." << std::endl
;
3843 //sort arguments of all transcendentals
3844 std::map
< Kind
, std::vector
< Node
> > sorted_tf_args
;
3845 std::map
< Kind
, std::map
< Node
, Node
> > tf_arg_to_term
;
3847 for( std::map
< Kind
, std::map
< Node
, Node
> >::iterator it
= d_tf_rep_map
.begin(); it
!= d_tf_rep_map
.end(); ++it
){
3849 if (k
== EXPONENTIAL
|| k
== SINE
)
3851 for (std::map
<Node
, Node
>::iterator itt
= it
->second
.begin();
3852 itt
!= it
->second
.end();
3855 Node a
= itt
->second
[0];
3856 computeModelValue(a
, 1);
3857 Assert(d_mv
[1].find(a
) != d_mv
[1].end());
3858 if (d_mv
[1][a
].isConst())
3860 Trace("nl-ext-tf-mono-debug") << "...tf term : " << a
<< std::endl
;
3861 sorted_tf_args
[k
].push_back(a
);
3862 tf_arg_to_term
[k
][a
] = itt
->second
;
3868 SortNonlinearExtension smv
;
3870 //sort by concrete values
3871 smv
.d_order_type
= 0;
3872 smv
.d_reverse_order
= true;
3873 for( std::map
< Kind
, std::map
< Node
, Node
> >::iterator it
= d_tf_rep_map
.begin(); it
!= d_tf_rep_map
.end(); ++it
){
3875 if( !sorted_tf_args
[k
].empty() ){
3876 std::sort( sorted_tf_args
[k
].begin(), sorted_tf_args
[k
].end(), smv
);
3877 Trace("nl-ext-tf-mono") << "Sorted transcendental function list for " << k
<< " : " << std::endl
;
3878 for( unsigned i
=0; i
<sorted_tf_args
[it
->first
].size(); i
++ ){
3879 Node targ
= sorted_tf_args
[k
][i
];
3880 Assert( d_mv
[1].find( targ
)!=d_mv
[1].end() );
3881 Trace("nl-ext-tf-mono") << " " << targ
<< " -> " << d_mv
[1][targ
] << std::endl
;
3882 Node t
= tf_arg_to_term
[k
][targ
];
3883 Assert( d_mv
[1].find( t
)!=d_mv
[1].end() );
3884 Trace("nl-ext-tf-mono") << " f-val : " << d_mv
[1][t
] << std::endl
;
3886 std::vector
< Node
> mpoints
;
3887 std::vector
< Node
> mpoints_vals
;
3888 if (it
->first
== SINE
)
3890 mpoints
.push_back( d_pi
);
3891 mpoints
.push_back( d_pi_2
);
3892 mpoints
.push_back(d_zero
);
3893 mpoints
.push_back( d_pi_neg_2
);
3894 mpoints
.push_back( d_pi_neg
);
3896 else if (it
->first
== EXPONENTIAL
)
3898 mpoints
.push_back( Node::null() );
3900 if( !mpoints
.empty() ){
3901 //get model values for points
3902 for( unsigned i
=0; i
<mpoints
.size(); i
++ ){
3904 if( !mpoints
[i
].isNull() ){
3905 mpv
= computeModelValue( mpoints
[i
], 1 );
3906 Assert( mpv
.isConst() );
3908 mpoints_vals
.push_back( mpv
);
3911 unsigned mdir_index
= 0;
3912 int monotonic_dir
= -1;
3913 Node mono_bounds
[2];
3914 Node targ
, targval
, t
, tval
;
3915 for( unsigned i
=0; i
<sorted_tf_args
[it
->first
].size(); i
++ ){
3916 Node sarg
= sorted_tf_args
[k
][i
];
3917 Assert( d_mv
[1].find( sarg
)!=d_mv
[1].end() );
3918 Node sargval
= d_mv
[1][sarg
];
3919 Assert( sargval
.isConst() );
3920 Node s
= tf_arg_to_term
[k
][ sarg
];
3921 Assert( d_mv
[1].find( s
)!=d_mv
[1].end() );
3922 Node sval
= d_mv
[1][s
];
3923 Assert( sval
.isConst() );
3925 //increment to the proper monotonicity region
3926 bool increment
= true;
3927 while (increment
&& mdir_index
< mpoints
.size())
3930 if( mpoints
[mdir_index
].isNull() ){
3933 Node pval
= mpoints_vals
[mdir_index
];
3934 Assert( pval
.isConst() );
3935 if( sargval
.getConst
<Rational
>() < pval
.getConst
<Rational
>() ){
3937 Trace("nl-ext-tf-mono") << "...increment at " << sarg
<< " since model value is less than " << mpoints
[mdir_index
] << std::endl
;
3941 tval
= Node::null();
3942 mono_bounds
[1] = mpoints
[mdir_index
];
3944 monotonic_dir
= regionToMonotonicityDir(it
->first
, mdir_index
);
3945 if (mdir_index
< mpoints
.size())
3947 mono_bounds
[0] = mpoints
[mdir_index
];
3949 mono_bounds
[0] = Node::null();
3953 // store the concavity region
3954 d_tf_region
[s
] = mdir_index
;
3955 Trace("nl-ext-concavity") << "Transcendental function " << s
3956 << " is in region #" << mdir_index
;
3957 Trace("nl-ext-concavity") << ", arg model value = " << sargval
3960 if( !tval
.isNull() ){
3962 if( monotonic_dir
==1 && sval
.getConst
<Rational
>() > tval
.getConst
<Rational
>() ){
3963 mono_lem
= NodeManager::currentNM()->mkNode(
3965 NodeManager::currentNM()->mkNode(GEQ
, targ
, sarg
),
3966 NodeManager::currentNM()->mkNode(GEQ
, t
, s
));
3967 }else if( monotonic_dir
==-1 && sval
.getConst
<Rational
>() < tval
.getConst
<Rational
>() ){
3968 mono_lem
= NodeManager::currentNM()->mkNode(
3970 NodeManager::currentNM()->mkNode(LEQ
, targ
, sarg
),
3971 NodeManager::currentNM()->mkNode(LEQ
, t
, s
));
3973 if( !mono_lem
.isNull() ){
3974 if( !mono_bounds
[0].isNull() ){
3975 Assert( !mono_bounds
[1].isNull() );
3976 mono_lem
= NodeManager::currentNM()->mkNode(
3978 NodeManager::currentNM()->mkNode(
3980 mkBounded(mono_bounds
[0], targ
, mono_bounds
[1]),
3981 mkBounded(mono_bounds
[0], sarg
, mono_bounds
[1])),
3984 Trace("nl-ext-tf-mono") << "Monotonicity lemma : " << mono_lem
<< std::endl
;
3985 lemmas
.push_back( mono_lem
);
3988 // store the previous values
4000 std::vector
<Node
> NonlinearExtension::checkTranscendentalTangentPlanes()
4002 std::vector
<Node
> lemmas
;
4003 Trace("nl-ext") << "Get tangent plane lemmas for transcendental functions..."
4005 // this implements Figure 3 of "Satisfiaility Modulo Transcendental Functions
4006 // via Incremental Linearization" by Cimatti et al
4007 for (std::pair
<const Kind
, std::map
<Node
, Node
> >& tfs
: d_tf_rep_map
)
4012 // We do not use Taylor approximation for PI currently.
4013 // This is because the convergence is extremely slow, and hence an
4014 // initial approximation is superior.
4017 Trace("nl-ext-tftp-debug2") << "Taylor variables: " << std::endl
;
4018 Trace("nl-ext-tftp-debug2")
4019 << " taylor_real_fv : " << d_taylor_real_fv
<< std::endl
;
4020 Trace("nl-ext-tftp-debug2")
4021 << " taylor_real_fv_base : " << d_taylor_real_fv_base
<< std::endl
;
4022 Trace("nl-ext-tftp-debug2")
4023 << " taylor_real_fv_base_rem : " << d_taylor_real_fv_base_rem
4025 Trace("nl-ext-tftp-debug2") << std::endl
;
4027 // we substitute into the Taylor sum P_{n,f(0)}( x )
4029 for (std::pair
<const Node
, Node
>& tfr
: tfs
.second
)
4031 // Figure 3 : tf( x )
4032 Node tf
= tfr
.second
;
4033 if (isRefineableTfFun(tf
))
4035 Trace("nl-ext-tftp") << "Compute tangent planes " << tf
<< std::endl
;
4036 // go until max degree is reached, or we don't meet bound criteria
4037 for (unsigned d
= 1; d
<= d_taylor_degree
; d
++)
4039 Trace("nl-ext-tftp") << "- run at degree " << d
<< "..." << std::endl
;
4040 unsigned prev
= lemmas
.size();
4041 if (!checkTfTangentPlanesFun(tf
, d
, lemmas
))
4043 Trace("nl-ext-tftp")
4044 << "...fail, #lemmas = " << (lemmas
.size() - prev
) << std::endl
;
4049 Trace("nl-ext-tftp") << "...success" << std::endl
;
4059 bool NonlinearExtension::isRefineableTfFun(Node tf
)
4061 Assert(tf
.getKind() == SINE
|| tf
.getKind() == EXPONENTIAL
);
4062 if (tf
.getKind() == SINE
)
4064 // we do not consider e.g. sin( -1*x ), since considering sin( x ) will
4065 // have the same effect
4072 Node c
= computeModelValue(tf
[0], 1);
4073 Assert(c
.isConst());
4074 int csign
= c
.getConst
<Rational
>().sgn();
4082 bool NonlinearExtension::checkTfTangentPlanesFun(Node tf
,
4084 std::vector
<Node
>& lemmas
)
4086 Assert(isRefineableTfFun(tf
));
4088 NodeManager
* nm
= NodeManager::currentNM();
4089 Kind k
= tf
.getKind();
4090 // Figure 3: P_l, P_u
4091 // mapped to for signs of c
4092 std::map
<int, Node
> poly_approx_bounds
[2];
4093 std::vector
<Node
> pbounds
;
4094 getPolynomialApproximationBounds(k
, d
, pbounds
);
4095 poly_approx_bounds
[0][1] = pbounds
[0];
4096 poly_approx_bounds
[0][-1] = pbounds
[1];
4097 poly_approx_bounds
[1][1] = pbounds
[2];
4098 poly_approx_bounds
[1][-1] = pbounds
[3];
4101 Node c
= computeModelValue(tf
[0], 1);
4102 int csign
= c
.getConst
<Rational
>().sgn();
4103 Assert(csign
== 1 || csign
== -1);
4106 Node v
= computeModelValue(tf
, 1);
4108 // check value of tf
4109 Trace("nl-ext-tftp-debug") << "Process tangent plane refinement for " << tf
4110 << ", degree " << d
<< "..." << std::endl
;
4111 Trace("nl-ext-tftp-debug") << " value in model : " << v
<< std::endl
;
4112 Trace("nl-ext-tftp-debug") << " arg value in model : " << c
<< std::endl
;
4114 std::vector
<Node
> taylor_vars
;
4115 taylor_vars
.push_back(d_taylor_real_fv
);
4117 // compute the concavity
4119 std::unordered_map
<Node
, int, NodeHashFunction
>::iterator itr
=
4120 d_tf_region
.find(tf
);
4121 if (itr
!= d_tf_region
.end())
4123 region
= itr
->second
;
4124 Trace("nl-ext-tftp-debug") << " region is : " << region
<< std::endl
;
4127 int concavity
= regionToConcavity(k
, itr
->second
);
4128 Trace("nl-ext-tftp-debug") << " concavity is : " << concavity
<< std::endl
;
4133 // bounds for which we are this concavity
4134 // Figure 3: < l, u >
4138 bounds
[0] = regionToLowerBound(k
, region
);
4139 Assert(!bounds
[0].isNull());
4140 bounds
[1] = regionToUpperBound(k
, region
);
4141 Assert(!bounds
[1].isNull());
4147 // compute whether this is a tangent refinement or a secant refinement
4148 bool is_tangent
= false;
4149 bool is_secant
= false;
4150 std::pair
<Node
, Node
> mvb
= getTfModelBounds(tf
, d
);
4151 for (unsigned r
= 0; r
< 2; r
++)
4153 Node pab
= poly_approx_bounds
[r
][csign
];
4154 Node v_pab
= r
== 0 ? mvb
.first
: mvb
.second
;
4155 if (!v_pab
.isNull())
4157 Assert(v_pab
.isConst());
4158 Trace("nl-ext-tftp-debug2") << "...model value of " << pab
<< " is "
4159 << v_pab
<< std::endl
;
4160 Node comp
= nm
->mkNode(r
== 0 ? LT
: GT
, v
, v_pab
);
4161 Trace("nl-ext-tftp-debug2") << "...compare : " << comp
<< std::endl
;
4162 Node compr
= Rewriter::rewrite(comp
);
4163 Trace("nl-ext-tftp-debug2") << "...got : " << compr
<< std::endl
;
4164 if (compr
== d_true
)
4166 // beyond the bounds
4169 poly_approx
= poly_approx_bounds
[r
][csign
];
4170 is_tangent
= concavity
== 1;
4171 is_secant
= concavity
== -1;
4175 poly_approx
= poly_approx_bounds
[r
][csign
];
4176 is_tangent
= concavity
== -1;
4177 is_secant
= concavity
== 1;
4179 if (Trace
.isOn("nl-ext-tftp"))
4181 Trace("nl-ext-tftp") << "*** Outside boundary point (";
4182 Trace("nl-ext-tftp") << (r
== 0 ? "low" : "high") << ") ";
4183 printRationalApprox("nl-ext-tftp", v_pab
);
4184 Trace("nl-ext-tftp") << ", will refine..." << std::endl
;
4185 Trace("nl-ext-tftp") << " poly_approx = " << poly_approx
4187 Trace("nl-ext-tftp") << " is_tangent = " << is_tangent
4189 Trace("nl-ext-tftp") << " is_secant = " << is_secant
<< std::endl
;
4195 Trace("nl-ext-tftp") << " ...within " << (r
== 0 ? "low" : "high")
4197 printRationalApprox("nl-ext-tftp", v_pab
);
4198 Trace("nl-ext-tftp") << std::endl
;
4205 if (is_tangent
|| is_secant
)
4207 Assert(!poly_approx
.isNull());
4208 std::vector
<Node
> taylor_subs
;
4209 taylor_subs
.push_back(c
);
4210 Assert(taylor_vars
.size() == taylor_subs
.size());
4211 poly_approx_c
= poly_approx
.substitute(taylor_vars
.begin(),
4213 taylor_subs
.begin(),
4215 Trace("nl-ext-tftp-debug2") << "...poly approximation at c is "
4216 << poly_approx_c
<< std::endl
;
4220 // we may want to continue getting better bounds
4226 // compute tangent plane
4229 Node poly_approx_deriv
= getDerivative(poly_approx
, d_taylor_real_fv
);
4230 Assert(!poly_approx_deriv
.isNull());
4231 poly_approx_deriv
= Rewriter::rewrite(poly_approx_deriv
);
4232 Trace("nl-ext-tftp-debug2") << "...derivative of " << poly_approx
<< " is "
4233 << poly_approx_deriv
<< std::endl
;
4234 std::vector
<Node
> taylor_subs
;
4235 taylor_subs
.push_back(c
);
4236 Assert(taylor_vars
.size() == taylor_subs
.size());
4237 Node poly_approx_c_deriv
= poly_approx_deriv
.substitute(taylor_vars
.begin(),
4239 taylor_subs
.begin(),
4241 tplane
= nm
->mkNode(
4244 nm
->mkNode(MULT
, poly_approx_c_deriv
, nm
->mkNode(MINUS
, tf
[0], c
)));
4246 Node lem
= nm
->mkNode(concavity
== 1 ? GEQ
: LEQ
, tf
, tplane
);
4247 std::vector
<Node
> antec
;
4248 for (unsigned i
= 0; i
< 2; i
++)
4250 if (!bounds
[i
].isNull())
4252 Node ant
= nm
->mkNode(i
== 0 ? GEQ
: LEQ
, tf
[0], bounds
[i
]);
4253 antec
.push_back(ant
);
4258 Node antec_n
= antec
.size() == 1 ? antec
[0] : nm
->mkNode(AND
, antec
);
4259 lem
= nm
->mkNode(IMPLIES
, antec_n
, lem
);
4261 Trace("nl-ext-tftp-debug2")
4262 << "*** Tangent plane lemma (pre-rewrite): " << lem
<< std::endl
;
4263 lem
= Rewriter::rewrite(lem
);
4264 Trace("nl-ext-tftp-lemma") << "*** Tangent plane lemma : " << lem
4266 Assert(computeModelValue(lem
, 1) == d_false
);
4267 // Figure 3 : line 9
4268 lemmas
.push_back(lem
);
4272 // bounds are the minimum and maximum previous secant points
4273 // should not repeat secant points: secant lemmas should suffice to
4274 // rule out previous assignment
4276 d_secant_points
[tf
][d
].begin(), d_secant_points
[tf
][d
].end(), c
)
4277 == d_secant_points
[tf
][d
].end());
4278 // insert into the vector
4279 d_secant_points
[tf
][d
].push_back(c
);
4281 SortNonlinearExtension smv
;
4283 smv
.d_order_type
= 0;
4285 d_secant_points
[tf
][d
].begin(), d_secant_points
[tf
][d
].end(), smv
);
4286 // get the resulting index of c
4289 d_secant_points
[tf
][d
].begin(), d_secant_points
[tf
][d
].end(), c
)
4290 - d_secant_points
[tf
][d
].begin();
4291 // bounds are the next closest upper/lower bound values
4294 bounds
[0] = d_secant_points
[tf
][d
][index
- 1];
4298 // otherwise, we use the lower boundary point for this concavity
4302 Assert(!bounds
[0].isNull());
4304 else if (k
== EXPONENTIAL
)
4307 bounds
[0] = Rewriter::rewrite(nm
->mkNode(MINUS
, c
, d_one
));
4310 if (index
< d_secant_points
[tf
][d
].size() - 1)
4312 bounds
[1] = d_secant_points
[tf
][d
][index
+ 1];
4316 // otherwise, we use the upper boundary point for this concavity
4320 Assert(!bounds
[1].isNull());
4322 else if (k
== EXPONENTIAL
)
4325 bounds
[1] = Rewriter::rewrite(nm
->mkNode(PLUS
, c
, d_one
));
4328 Trace("nl-ext-tftp-debug2") << "...secant bounds are : " << bounds
[0]
4329 << " ... " << bounds
[1] << std::endl
;
4331 for (unsigned s
= 0; s
< 2; s
++)
4333 // compute secant plane
4334 Assert(!poly_approx
.isNull());
4335 Assert(!bounds
[s
].isNull());
4336 // take the model value of l or u (since may contain PI)
4337 Node b
= computeModelValue(bounds
[s
], 1);
4338 Trace("nl-ext-tftp-debug2") << "...model value of bound " << bounds
[s
]
4339 << " is " << b
<< std::endl
;
4340 Assert(b
.isConst());
4343 // Figure 3 : P(l), P(u), for s = 0,1
4345 std::vector
<Node
> taylor_subs
;
4346 taylor_subs
.push_back(b
);
4347 Assert(taylor_vars
.size() == taylor_subs
.size());
4348 poly_approx_b
= poly_approx
.substitute(taylor_vars
.begin(),
4350 taylor_subs
.begin(),
4352 // Figure 3: S_l( x ), S_u( x ) for s = 0,1
4354 Node rcoeff_n
= Rewriter::rewrite(nm
->mkNode(MINUS
, b
, c
));
4355 Assert(rcoeff_n
.isConst());
4356 Rational rcoeff
= rcoeff_n
.getConst
<Rational
>();
4357 Assert(rcoeff
.sgn() != 0);
4358 splane
= nm
->mkNode(
4362 nm
->mkNode(MINUS
, poly_approx_b
, poly_approx_c
),
4363 nm
->mkConst(Rational(1) / rcoeff
),
4364 nm
->mkNode(MINUS
, tf
[0], b
)));
4366 Node lem
= nm
->mkNode(concavity
== 1 ? LEQ
: GEQ
, tf
, splane
);
4367 // With respect to Figure 3, this is slightly different.
4368 // In particular, we chose b to be the model value of bounds[s],
4369 // which is a constant although bounds[s] may not be (e.g. if it
4371 // To ensure that c...b does not cross an inflection point,
4372 // we guard with the symbolic version of bounds[s].
4373 // This leads to lemmas e.g. of this form:
4374 // ( c <= x <= PI/2 ) => ( sin(x) < ( P( b ) - P( c ) )*( x -
4376 // where b = (PI/2)^M, the current value of PI/2 in the model.
4377 // This is sound since we are guarded by the symbolic
4378 // representation of PI/2.
4381 nm
->mkNode(GEQ
, tf
[0], s
== 0 ? bounds
[s
] : c
),
4382 nm
->mkNode(LEQ
, tf
[0], s
== 0 ? c
: bounds
[s
]));
4383 lem
= nm
->mkNode(IMPLIES
, antec_n
, lem
);
4384 Trace("nl-ext-tftp-debug2")
4385 << "*** Secant plane lemma (pre-rewrite) : " << lem
<< std::endl
;
4386 lem
= Rewriter::rewrite(lem
);
4387 Trace("nl-ext-tftp-lemma") << "*** Secant plane lemma : " << lem
4389 // Figure 3 : line 22
4390 lemmas
.push_back(lem
);
4391 Assert(computeModelValue(lem
, 1) == d_false
);
4398 int NonlinearExtension::regionToMonotonicityDir(Kind k
, int region
)
4400 if (k
== EXPONENTIAL
)
4409 if (region
== 1 || region
== 4)
4413 else if (region
== 2 || region
== 3)
4421 int NonlinearExtension::regionToConcavity(Kind k
, int region
)
4423 if (k
== EXPONENTIAL
)
4432 if (region
== 1 || region
== 2)
4436 else if (region
== 3 || region
== 4)
4444 Node
NonlinearExtension::regionToLowerBound(Kind k
, int region
)
4452 else if (region
== 2)
4456 else if (region
== 3)
4460 else if (region
== 4)
4465 return Node::null();
4468 Node
NonlinearExtension::regionToUpperBound(Kind k
, int region
)
4476 else if (region
== 2)
4480 else if (region
== 3)
4484 else if (region
== 4)
4489 return Node::null();
4492 Node
NonlinearExtension::getDerivative(Node n
, Node x
)
4495 // only handle the cases of the taylor expansion of d
4496 if (n
.getKind() == EXPONENTIAL
)
4503 else if (n
.getKind() == SINE
)
4507 Node na
= NodeManager::currentNM()->mkNode(MINUS
, d_pi_2
, n
[0]);
4508 Node ret
= NodeManager::currentNM()->mkNode(SINE
, na
);
4509 ret
= Rewriter::rewrite(ret
);
4513 else if (n
.getKind() == PLUS
)
4515 std::vector
<Node
> dchildren
;
4516 for (unsigned i
= 0; i
< n
.getNumChildren(); i
++)
4518 // PLUS is flattened in rewriter, recursion depth is bounded by 1
4519 Node dc
= getDerivative(n
[i
], x
);
4524 dchildren
.push_back(dc
);
4527 return NodeManager::currentNM()->mkNode(PLUS
, dchildren
);
4529 else if (n
.getKind() == MULT
)
4531 Assert(n
[0].isConst());
4532 Node dc
= getDerivative(n
[1], x
);
4535 return NodeManager::currentNM()->mkNode(MULT
, n
[0], dc
);
4538 else if (n
.getKind() == NONLINEAR_MULT
)
4540 unsigned xcount
= 0;
4541 std::vector
<Node
> children
;
4542 unsigned xindex
= 0;
4543 for (unsigned i
= 0, size
= n
.getNumChildren(); i
< size
; i
++)
4550 children
.push_back(n
[i
]);
4558 children
[xindex
] = NodeManager::currentNM()->mkConst(Rational(xcount
));
4560 return NodeManager::currentNM()->mkNode(MULT
, children
);
4564 return n
== x
? d_one
: d_zero
;
4566 else if (n
.isConst())
4570 Trace("nl-ext-debug") << "No derivative computed for " << n
;
4571 Trace("nl-ext-debug") << " for d/d{" << x
<< "}" << std::endl
;
4572 return Node::null();
4575 std::pair
<Node
, Node
> NonlinearExtension::getTaylor(Node fa
, unsigned n
)
4578 Node fac
; // what term we cache for fa
4579 if (fa
[0] == d_zero
)
4581 // optimization : simpler to compute (x-fa[0])^n if we are centered around 0
4586 // otherwise we use a standard factor a in (x-a)^n
4587 fac
= NodeManager::currentNM()->mkNode(fa
.getKind(), d_taylor_real_fv_base
);
4591 // check if we have already computed this Taylor series
4592 std::unordered_map
<unsigned, Node
>::iterator itt
= d_taylor_sum
[fac
].find(n
);
4593 if (itt
== d_taylor_sum
[fac
].end())
4596 if (fa
[0] == d_zero
)
4598 i_exp_base
= d_taylor_real_fv
;
4602 i_exp_base
= Rewriter::rewrite(NodeManager::currentNM()->mkNode(
4603 MINUS
, d_taylor_real_fv
, d_taylor_real_fv_base
));
4606 Node i_fact
= d_one
;
4608 int i_derv_status
= 0;
4609 unsigned counter
= 0;
4610 std::vector
<Node
> sum
;
4614 if (fa
.getKind() == EXPONENTIAL
)
4618 else if (fa
.getKind() == SINE
)
4620 if (i_derv_status
% 2 == 1)
4622 Node arg
= NodeManager::currentNM()->mkNode(
4623 PLUS
, d_pi_2
, d_taylor_real_fv_base
);
4624 i_derv
= NodeManager::currentNM()->mkNode(SINE
, arg
);
4630 if (i_derv_status
>= 2)
4632 i_derv
= NodeManager::currentNM()->mkNode(MINUS
, d_zero
, i_derv
);
4634 i_derv
= Rewriter::rewrite(i_derv
);
4635 i_derv_status
= i_derv_status
== 3 ? 0 : i_derv_status
+ 1;
4637 if (counter
== (n
+ 1))
4639 TNode x
= d_taylor_real_fv_base
;
4640 i_derv
= i_derv
.substitute(x
, d_taylor_real_fv_base_rem
);
4642 Node curr
= NodeManager::currentNM()->mkNode(
4644 NodeManager::currentNM()->mkNode(DIVISION
, i_derv
, i_fact
),
4646 if (counter
== (n
+ 1))
4652 sum
.push_back(curr
);
4653 i_fact
= Rewriter::rewrite(NodeManager::currentNM()->mkNode(
4655 NodeManager::currentNM()->mkConst(Rational(counter
)),
4657 i_exp
= Rewriter::rewrite(
4658 NodeManager::currentNM()->mkNode(MULT
, i_exp_base
, i_exp
));
4660 } while (counter
<= n
);
4662 sum
.size() == 1 ? sum
[0] : NodeManager::currentNM()->mkNode(PLUS
, sum
);
4664 if (fac
[0] != d_taylor_real_fv_base
)
4666 TNode x
= d_taylor_real_fv_base
;
4667 taylor_sum
= taylor_sum
.substitute(x
, fac
[0]);
4671 d_taylor_sum
[fac
][n
] = taylor_sum
;
4672 d_taylor_rem
[fac
][n
] = taylor_rem
;
4676 taylor_sum
= itt
->second
;
4677 Assert(d_taylor_rem
[fac
].find(n
) != d_taylor_rem
[fac
].end());
4678 taylor_rem
= d_taylor_rem
[fac
][n
];
4681 // must substitute for the argument if we were using a different lookup
4682 if (fa
[0] != fac
[0])
4684 TNode x
= d_taylor_real_fv_base
;
4685 taylor_sum
= taylor_sum
.substitute(x
, fa
[0]);
4687 return std::pair
<Node
, Node
>(taylor_sum
, taylor_rem
);
4690 void NonlinearExtension::getPolynomialApproximationBounds(
4691 Kind k
, unsigned d
, std::vector
<Node
>& pbounds
)
4693 if (d_poly_bounds
[k
][d
].empty())
4695 NodeManager
* nm
= NodeManager::currentNM();
4696 Node tft
= nm
->mkNode(k
, d_zero
);
4697 // n is the Taylor degree we are currently considering
4700 std::pair
<Node
, Node
> taylor
= getTaylor(tft
, n
);
4701 Trace("nl-ext-tftp-debug2") << "Taylor for " << k
4702 << " is : " << taylor
.first
<< std::endl
;
4703 Node taylor_sum
= Rewriter::rewrite(taylor
.first
);
4704 Trace("nl-ext-tftp-debug2") << "Taylor for " << k
4705 << " is (post-rewrite) : " << taylor_sum
4707 Assert(taylor
.second
.getKind() == MULT
);
4708 Assert(taylor
.second
.getNumChildren() == 2);
4709 Assert(taylor
.second
[0].getKind() == DIVISION
);
4710 Trace("nl-ext-tftp-debug2") << "Taylor remainder for " << k
<< " is "
4711 << taylor
.second
<< std::endl
;
4712 // ru is x^{n+1}/(n+1)!
4713 Node ru
= nm
->mkNode(DIVISION
, taylor
.second
[1], taylor
.second
[0][1]);
4714 ru
= Rewriter::rewrite(ru
);
4715 Trace("nl-ext-tftp-debug2")
4716 << "Taylor remainder factor is (post-rewrite) : " << ru
<< std::endl
;
4717 if (k
== EXPONENTIAL
)
4719 pbounds
.push_back(taylor_sum
);
4720 pbounds
.push_back(taylor_sum
);
4721 pbounds
.push_back(Rewriter::rewrite(
4722 nm
->mkNode(MULT
, taylor_sum
, nm
->mkNode(PLUS
, d_one
, ru
))));
4723 pbounds
.push_back(Rewriter::rewrite(nm
->mkNode(PLUS
, taylor_sum
, ru
)));
4728 Node l
= Rewriter::rewrite(nm
->mkNode(MINUS
, taylor_sum
, ru
));
4729 Node u
= Rewriter::rewrite(nm
->mkNode(PLUS
, taylor_sum
, ru
));
4730 pbounds
.push_back(l
);
4731 pbounds
.push_back(l
);
4732 pbounds
.push_back(u
);
4733 pbounds
.push_back(u
);
4735 Trace("nl-ext-tf-tplanes") << "Polynomial approximation for " << k
4736 << " is: " << std::endl
;
4737 Trace("nl-ext-tf-tplanes") << " Lower (pos): " << pbounds
[0] << std::endl
;
4738 Trace("nl-ext-tf-tplanes") << " Upper (pos): " << pbounds
[2] << std::endl
;
4739 Trace("nl-ext-tf-tplanes") << " Lower (neg): " << pbounds
[1] << std::endl
;
4740 Trace("nl-ext-tf-tplanes") << " Upper (neg): " << pbounds
[3] << std::endl
;
4741 d_poly_bounds
[k
][d
].insert(
4742 d_poly_bounds
[k
][d
].end(), pbounds
.begin(), pbounds
.end());
4747 pbounds
.end(), d_poly_bounds
[k
][d
].begin(), d_poly_bounds
[k
][d
].end());
4751 std::pair
<Node
, Node
> NonlinearExtension::getTfModelBounds(Node tf
, unsigned d
)
4753 // compute the model value of the argument
4754 Node c
= computeModelValue(tf
[0], 1);
4755 Assert(c
.isConst());
4756 int csign
= c
.getConst
<Rational
>().sgn();
4758 bool isNeg
= csign
== -1;
4760 std::vector
<Node
> pbounds
;
4761 getPolynomialApproximationBounds(tf
.getKind(), d
, pbounds
);
4763 std::vector
<Node
> bounds
;
4764 TNode tfv
= d_taylor_real_fv
;
4766 for (unsigned d
= 0; d
< 2; d
++)
4768 int index
= d
== 0 ? (isNeg
? 1 : 0) : (isNeg
? 3 : 2);
4769 Node pab
= pbounds
[index
];
4773 pab
= pab
.substitute(tfv
, tfs
);
4774 pab
= Rewriter::rewrite(pab
);
4775 Node v_pab
= computeModelValue(pab
, 1);
4776 bounds
.push_back(v_pab
);
4780 bounds
.push_back(Node::null());
4783 return std::pair
<Node
, Node
>(bounds
[0], bounds
[1]);
4786 } // namespace arith
4787 } // namespace theory