Head

GCC Code Coverage Report

Directory:	.		Exec	Total	Coverage
File:	src/theory/arith/nl/transcendental/taylor_generator.cpp	Lines:	118	126	93.7 %
Date:	2021-08-16	Branches:	238	550	43.3 %


/******************************************************************************
 * Top contributors (to current version):
 *   Gereon Kremer, Andrew Reynolds
 *
 * This file is part of the cvc5 project.
 *
 * Copyright (c) 2009-2021 by the authors listed in the file AUTHORS
 * in the top-level source directory and their institutional affiliations.
 * All rights reserved.  See the file COPYING in the top-level source
 * directory for licensing information.
 * ****************************************************************************
 *
 * Generate taylor approximations transcendental lemmas.
 */

#include "theory/arith/nl/transcendental/taylor_generator.h"

#include "theory/arith/arith_utilities.h"
#include "theory/arith/nl/nl_model.h"
#include "theory/rewriter.h"

namespace cvc5 {
namespace theory {
namespace arith {
namespace nl {
namespace transcendental {

TaylorGenerator::TaylorGenerator()
    : d_nm(NodeManager::currentNM()),
      d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType()))
{
}

TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }

std::pair<Node, Node> TaylorGenerator::getTaylor(Kind k, std::uint64_t n)
{
  Assert(n > 0);
  // check if we have already computed this Taylor series
  auto itt = d_taylor_terms[k].find(n);
  if (itt != d_taylor_terms[k].end())
  {
    return itt->second;
  }

  NodeManager* nm = NodeManager::currentNM();

  // the current factorial `counter!`
  Integer factorial = 1;
  // the current variable power `x^counter`
  Node varpow = nm->mkConst(Rational(1));
  std::vector<Node> sum;
  for (std::uint64_t counter = 1; counter <= n; ++counter)
  {
    if (k == Kind::EXPONENTIAL)
    {
      // Maclaurin series for exponential:
      //   \sum_{n=0}^\infty x^n / n!
      sum.push_back(
          nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
    }
    else if (k == Kind::SINE)
    {
      // Maclaurin series for exponential:
      //   \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1)
      if (counter % 2 == 0)
      {
        int sign = (counter % 4 == 0 ? -1 : 1);
        sum.push_back(nm->mkNode(Kind::MULT,
                                 nm->mkNode(Kind::DIVISION,
                                            nm->mkConst<Rational>(sign),
                                            nm->mkConst<Rational>(factorial)),
                                 varpow));
      }
    }
    factorial *= counter;
    varpow =
        Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow));
  }
  Node taylor_sum =
      Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum));
  Node taylor_rem = Rewriter::rewrite(
      nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));

  auto res = std::make_pair(taylor_sum, taylor_rem);

  // put result in cache
  d_taylor_terms[k][n] = res;

  return res;
}

void TaylorGenerator::getPolynomialApproximationBounds(
    Kind k, std::uint64_t d, ApproximationBounds& pbounds)
{
  auto it = d_poly_bounds[k].find(d);
  if (it == d_poly_bounds[k].end())
  {
    NodeManager* nm = NodeManager::currentNM();
    // n is the Taylor degree we are currently considering
    std::uint64_t n = 2 * d;
    // n must be even
    std::pair<Node, Node> taylor = getTaylor(k, n);
    Node taylor_sum = taylor.first;
    // ru is x^{n+1}/(n+1)!
    Node ru = taylor.second;
    Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first
                      << std::endl;
    Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second
                      << std::endl;
    if (k == Kind::EXPONENTIAL)
    {
      pbounds.d_lower = taylor_sum;
      pbounds.d_upperNeg =
          Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
      pbounds.d_upperPos = Rewriter::rewrite(
          nm->mkNode(Kind::MULT,
                     taylor_sum,
                     nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)));
    }
    else
    {
      Assert(k == Kind::SINE);
      Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
      Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
      pbounds.d_lower = l;
      pbounds.d_upperNeg = u;
      pbounds.d_upperPos = u;
    }
    Trace("nl-trans") << "Polynomial approximation for " << k
                      << " is: " << std::endl;
    Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl;
    Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl;
    Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl;
    d_poly_bounds[k].emplace(d, pbounds);
  }
  else
  {
    pbounds = it->second;
  }
}

std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg(
    Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds)
{
  getPolynomialApproximationBounds(k, d, pbounds);
  Trace("nl-trans") << "c = " << c << std::endl;
  Assert(c.isConst());
  if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
  {
    bool success = false;
    std::uint64_t ds = d;
    TNode ttrf = getTaylorVariable();
    TNode tc = c;
    do
    {
      success = true;
      std::uint64_t n = 2 * ds;
      std::pair<Node, Node> taylor = getTaylor(k, n);
      // check that 1-c^{n+1}/(n+1)! > 0
      Node ru = taylor.second;
      Node rus = ru.substitute(ttrf, tc);
      rus = Rewriter::rewrite(rus);
      Assert(rus.isConst());
      if (rus.getConst<Rational>() > 1)
      {
        success = false;
        ds = ds + 1;
      }
    } while (!success);
    if (ds > d)
    {
      Trace("nl-ext-exp-taylor")
          << "*** Increase Taylor bound to " << ds << " > " << d << " for ("
          << k << " " << c << ")" << std::endl;
      // must use sound upper bound
      ApproximationBounds pboundss;
      getPolynomialApproximationBounds(k, ds, pboundss);
      pbounds.d_upperPos = pboundss.d_upperPos;
    }
    return ds;
  }
  return d;
}

std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf,
                                                        std::uint64_t d,
                                                        NlModel& model)
{
  // compute the model value of the argument
  Node c = model.computeAbstractModelValue(tf[0]);
  Assert(c.isConst());
  int csign = c.getConst<Rational>().sgn();
  Kind k = tf.getKind();
  if (csign == 0)
  {
    NodeManager* nm = NodeManager::currentNM();
    // at zero, its trivial
    if (k == Kind::SINE)
    {
      Node zero = nm->mkConst(Rational(0));
      return std::pair<Node, Node>(zero, zero);
    }
    Assert(k == Kind::EXPONENTIAL);
    Node one = nm->mkConst(Rational(1));
    return std::pair<Node, Node>(one, one);
  }
  bool isNeg = csign == -1;

  ApproximationBounds pbounds;
  getPolynomialApproximationBoundForArg(k, c, d, pbounds);

  std::vector<Node> bounds;
  TNode tfv = getTaylorVariable();
  TNode tfs = tf[0];
  for (unsigned d2 = 0; d2 < 2; d2++)
  {
    Node pab = (d2 == 0 ? pbounds.d_lower
                        : (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos));
    if (!pab.isNull())
    {
      // { x -> M_A(tf[0]) }
      // Notice that we compute the model value of tfs first, so that
      // the call to rewrite below does not modify the term, where notice that
      // rewrite( x*x { x -> M_A(t) } ) = M_A(t)*M_A(t)
      // is not equal to
      // M_A( x*x { x -> t } ) = M_A( t*t )
      // where M_A denotes the abstract model.
      Node mtfs = model.computeAbstractModelValue(tfs);
      pab = pab.substitute(tfv, mtfs);
      pab = Rewriter::rewrite(pab);
      Assert(pab.isConst());
      bounds.push_back(pab);
    }
    else
    {
      bounds.push_back(Node::null());
    }
  }
  return std::pair<Node, Node>(bounds[0], bounds[1]);
}

}  // namespace transcendental
}  // namespace nl
}  // namespace arith
}  // namespace theory
}  // namespace cvc5


Generated by: GCOVR (Version 3.2)

Line	Exec	Source
1		/******************************************************************************
2		* Top contributors (to current version):
3		* Gereon Kremer, Andrew Reynolds
4		*
5		* This file is part of the cvc5 project.
6		*
7		* Copyright (c) 2009-2021 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.
11		* ****************************************************************************
12		*
13		* Generate taylor approximations transcendental lemmas.
14		*/
15
16		#include "theory/arith/nl/transcendental/taylor_generator.h"
17
18		#include "theory/arith/arith_utilities.h"
19		#include "theory/arith/nl/nl_model.h"
20		#include "theory/rewriter.h"
21
22		namespace cvc5 {
23		namespace theory {
24		namespace arith {
25		namespace nl {
26		namespace transcendental {
27
28	5212	TaylorGenerator::TaylorGenerator()
29	5212	: d_nm(NodeManager::currentNM()),
30	5212	d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType()))
31		{
32	5212	}
33
34	3351	TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }
35
36	3496	std::pair<Node, Node> TaylorGenerator::getTaylor(Kind k, std::uint64_t n)
37		{
38	3496	Assert(n > 0);
39		// check if we have already computed this Taylor series
40	3496	auto itt = d_taylor_terms[k].find(n);
41	3496	if (itt != d_taylor_terms[k].end())
42		{
43	3254	return itt->second;
44		}
45
46	242	NodeManager* nm = NodeManager::currentNM();
47
48		// the current factorial `counter!`
49	484	Integer factorial = 1;
50		// the current variable power `x^counter`
51	484	Node varpow = nm->mkConst(Rational(1));
52	484	std::vector<Node> sum;
53	1642	for (std::uint64_t counter = 1; counter <= n; ++counter)
54		{
55	1400	if (k == Kind::EXPONENTIAL)
56		{
57		// Maclaurin series for exponential:
58		// \sum_{n=0}^\infty x^n / n!
59	788	sum.push_back(
60	1576	nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
61		}
62	612	else if (k == Kind::SINE)
63		{
64		// Maclaurin series for exponential:
65		// \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1)
66	612	if (counter % 2 == 0)
67		{
68	306	int sign = (counter % 4 == 0 ? -1 : 1);
69	918	sum.push_back(nm->mkNode(Kind::MULT,
70	1224	nm->mkNode(Kind::DIVISION,
71	612	nm->mkConst<Rational>(sign),
72	612	nm->mkConst<Rational>(factorial)),
73		varpow));
74		}
75		}
76	1400	factorial *= counter;
77	1400	varpow =
78	2800	Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow));
79		}
80		Node taylor_sum =
81	484	Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum));
82		Node taylor_rem = Rewriter::rewrite(
83	484	nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
84
85	484	auto res = std::make_pair(taylor_sum, taylor_rem);
86
87		// put result in cache
88	242	d_taylor_terms[k][n] = res;
89
90	242	return res;
91		}
92
93	3240	void TaylorGenerator::getPolynomialApproximationBounds(
94		Kind k, std::uint64_t d, ApproximationBounds& pbounds)
95		{
96	3240	auto it = d_poly_bounds[k].find(d);
97	3240	if (it == d_poly_bounds[k].end())
98		{
99	242	NodeManager* nm = NodeManager::currentNM();
100		// n is the Taylor degree we are currently considering
101	242	std::uint64_t n = 2 * d;
102		// n must be even
103	484	std::pair<Node, Node> taylor = getTaylor(k, n);
104	484	Node taylor_sum = taylor.first;
105		// ru is x^{n+1}/(n+1)!
106	484	Node ru = taylor.second;
107	484	Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first
108	242	<< std::endl;
109	484	Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second
110	242	<< std::endl;
111	242	if (k == Kind::EXPONENTIAL)
112		{
113	130	pbounds.d_lower = taylor_sum;
114	130	pbounds.d_upperNeg =
115	260	Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
116	260	pbounds.d_upperPos = Rewriter::rewrite(
117	260	nm->mkNode(Kind::MULT,
118		taylor_sum,
119	390	nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)));
120		}
121		else
122		{
123	112	Assert(k == Kind::SINE);
124	224	Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
125	224	Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
126	112	pbounds.d_lower = l;
127	112	pbounds.d_upperNeg = u;
128	112	pbounds.d_upperPos = u;
129		}
130	484	Trace("nl-trans") << "Polynomial approximation for " << k
131	242	<< " is: " << std::endl;
132	242	Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl;
133	242	Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl;
134	242	Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl;
135	242	d_poly_bounds[k].emplace(d, pbounds);
136		}
137		else
138		{
139	2998	pbounds = it->second;
140		}
141	3240	}
142
143	2533	std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg(
144		Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds)
145		{
146	2533	getPolynomialApproximationBounds(k, d, pbounds);
147	2533	Trace("nl-trans") << "c = " << c << std::endl;
148	2533	Assert(c.isConst());
149	2533	if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
150		{
151	1604	bool success = false;
152	1604	std::uint64_t ds = d;
153	3208	TNode ttrf = getTaylorVariable();
154	3208	TNode tc = c;
155	1650	do
156		{
157	3254	success = true;
158	3254	std::uint64_t n = 2 * ds;
159	6508	std::pair<Node, Node> taylor = getTaylor(k, n);
160		// check that 1-c^{n+1}/(n+1)! > 0
161	6508	Node ru = taylor.second;
162	6508	Node rus = ru.substitute(ttrf, tc);
163	3254	rus = Rewriter::rewrite(rus);
164	3254	Assert(rus.isConst());
165	3254	if (rus.getConst<Rational>() > 1)
166		{
167	1650	success = false;
168	1650	ds = ds + 1;
169		}
170	3254	} while (!success);
171	1604	if (ds > d)
172		{
173	1280	Trace("nl-ext-exp-taylor")
174	640	<< "*** Increase Taylor bound to " << ds << " > " << d << " for ("
175	640	<< k << " " << c << ")" << std::endl;
176		// must use sound upper bound
177	1280	ApproximationBounds pboundss;
178	640	getPolynomialApproximationBounds(k, ds, pboundss);
179	640	pbounds.d_upperPos = pboundss.d_upperPos;
180		}
181	1604	return ds;
182		}
183	929	return d;
184		}
185
186	1402	std::pair<Node, Node> TaylorGenerator::getTfModelBounds(Node tf,
187		std::uint64_t d,
188		NlModel& model)
189		{
190		// compute the model value of the argument
191	2804	Node c = model.computeAbstractModelValue(tf[0]);
192	1402	Assert(c.isConst());
193	1402	int csign = c.getConst<Rational>().sgn();
194	1402	Kind k = tf.getKind();
195	1402	if (csign == 0)
196		{
197		NodeManager* nm = NodeManager::currentNM();
198		// at zero, its trivial
199		if (k == Kind::SINE)
200		{
201		Node zero = nm->mkConst(Rational(0));
202		return std::pair<Node, Node>(zero, zero);
203		}
204		Assert(k == Kind::EXPONENTIAL);
205		Node one = nm->mkConst(Rational(1));
206		return std::pair<Node, Node>(one, one);
207		}
208	1402	bool isNeg = csign == -1;
209
210	2804	ApproximationBounds pbounds;
211	1402	getPolynomialApproximationBoundForArg(k, c, d, pbounds);
212
213	2804	std::vector<Node> bounds;
214	2804	TNode tfv = getTaylorVariable();
215	2804	TNode tfs = tf[0];
216	4206	for (unsigned d2 = 0; d2 < 2; d2++)
217		{
218		Node pab = (d2 == 0 ? pbounds.d_lower
219	5608	: (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos));
220	2804	if (!pab.isNull())
221		{
222		// { x -> M_A(tf[0]) }
223		// Notice that we compute the model value of tfs first, so that
224		// the call to rewrite below does not modify the term, where notice that
225		// rewrite( xx { x -> M_A(t) } ) = M_A(t)M_A(t)
226		// is not equal to
227		// M_A( xx { x -> t } ) = M_A( tt )
228		// where M_A denotes the abstract model.
229	5608	Node mtfs = model.computeAbstractModelValue(tfs);
230	2804	pab = pab.substitute(tfv, mtfs);
231	2804	pab = Rewriter::rewrite(pab);
232	2804	Assert(pab.isConst());
233	2804	bounds.push_back(pab);
234		}
235		else
236		{
237		bounds.push_back(Node::null());
238		}
239		}
240	1402	return std::pair<Node, Node>(bounds[0], bounds[1]);
241		}
242
243		} // namespace transcendental
244		} // namespace nl
245		} // namespace arith
246		} // namespace theory
247	29340	} // namespace cvc5