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-09-29	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	3429	TaylorGenerator::TaylorGenerator()
29	3429	: d_nm(NodeManager::currentNM()),
30	3429	d_taylor_real_fv(d_nm->mkBoundVar("x", d_nm->realType()))
31		{
32	3429	}
33
34	2574	TNode TaylorGenerator::getTaylorVariable() { return d_taylor_real_fv; }
35
36	3109	std::pair<Node, Node> TaylorGenerator::getTaylor(Kind k, std::uint64_t n)
37		{
38	3109	Assert(n > 0);
39		// check if we have already computed this Taylor series
40	3109	auto itt = d_taylor_terms[k].find(n);
41	3109	if (itt != d_taylor_terms[k].end())
42		{
43	3038	return itt->second;
44		}
45
46	71	NodeManager* nm = NodeManager::currentNM();
47
48		// the current factorial `counter!`
49	142	Integer factorial = 1;
50		// the current variable power `x^counter`
51	142	Node varpow = nm->mkConst(Rational(1));
52	142	std::vector<Node> sum;
53	573	for (std::uint64_t counter = 1; counter <= n; ++counter)
54		{
55	502	if (k == Kind::EXPONENTIAL)
56		{
57		// Maclaurin series for exponential:
58		// \sum_{n=0}^\infty x^n / n!
59	328	sum.push_back(
60	656	nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
61		}
62	174	else if (k == Kind::SINE)
63		{
64		// Maclaurin series for exponential:
65		// \sum_{n=0}^\infty (-1)^n / (2n+1)! * x^(2n+1)
66	174	if (counter % 2 == 0)
67		{
68	87	int sign = (counter % 4 == 0 ? -1 : 1);
69	261	sum.push_back(nm->mkNode(Kind::MULT,
70	348	nm->mkNode(Kind::DIVISION,
71	174	nm->mkConst<Rational>(sign),
72	174	nm->mkConst<Rational>(factorial)),
73		varpow));
74		}
75		}
76	502	factorial *= counter;
77	502	varpow =
78	1004	Rewriter::rewrite(nm->mkNode(Kind::MULT, d_taylor_real_fv, varpow));
79		}
80		Node taylor_sum =
81	142	Rewriter::rewrite(sum.size() == 1 ? sum[0] : nm->mkNode(Kind::PLUS, sum));
82		Node taylor_rem = Rewriter::rewrite(
83	142	nm->mkNode(Kind::DIVISION, varpow, nm->mkConst<Rational>(factorial)));
84
85	142	auto res = std::make_pair(taylor_sum, taylor_rem);
86
87		// put result in cache
88	71	d_taylor_terms[k][n] = res;
89
90	71	return res;
91		}
92
93	2526	void TaylorGenerator::getPolynomialApproximationBounds(
94		Kind k, std::uint64_t d, ApproximationBounds& pbounds)
95		{
96	2526	auto it = d_poly_bounds[k].find(d);
97	2526	if (it == d_poly_bounds[k].end())
98		{
99	71	NodeManager* nm = NodeManager::currentNM();
100		// n is the Taylor degree we are currently considering
101	71	std::uint64_t n = 2 * d;
102		// n must be even
103	142	std::pair<Node, Node> taylor = getTaylor(k, n);
104	142	Node taylor_sum = taylor.first;
105		// ru is x^{n+1}/(n+1)!
106	142	Node ru = taylor.second;
107	142	Trace("nl-trans") << "Taylor for " << k << " is : " << taylor.first
108	71	<< std::endl;
109	142	Trace("nl-trans") << "Taylor remainder for " << k << " is " << taylor.second
110	71	<< std::endl;
111	71	if (k == Kind::EXPONENTIAL)
112		{
113	41	pbounds.d_lower = taylor_sum;
114	41	pbounds.d_upperNeg =
115	82	Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
116	82	pbounds.d_upperPos = Rewriter::rewrite(
117	82	nm->mkNode(Kind::MULT,
118		taylor_sum,
119	123	nm->mkNode(Kind::PLUS, nm->mkConst(Rational(1)), ru)));
120		}
121		else
122		{
123	30	Assert(k == Kind::SINE);
124	60	Node l = Rewriter::rewrite(nm->mkNode(Kind::MINUS, taylor_sum, ru));
125	60	Node u = Rewriter::rewrite(nm->mkNode(Kind::PLUS, taylor_sum, ru));
126	30	pbounds.d_lower = l;
127	30	pbounds.d_upperNeg = u;
128	30	pbounds.d_upperPos = u;
129		}
130	142	Trace("nl-trans") << "Polynomial approximation for " << k
131	71	<< " is: " << std::endl;
132	71	Trace("nl-trans") << " Lower: " << pbounds.d_lower << std::endl;
133	71	Trace("nl-trans") << " Upper (neg): " << pbounds.d_upperNeg << std::endl;
134	71	Trace("nl-trans") << " Upper (pos): " << pbounds.d_upperPos << std::endl;
135	71	d_poly_bounds[k].emplace(d, pbounds);
136		}
137		else
138		{
139	2455	pbounds = it->second;
140		}
141	2526	}
142
143	1907	std::uint64_t TaylorGenerator::getPolynomialApproximationBoundForArg(
144		Kind k, Node c, std::uint64_t d, ApproximationBounds& pbounds)
145		{
146	1907	getPolynomialApproximationBounds(k, d, pbounds);
147	1907	Trace("nl-trans") << "c = " << c << std::endl;
148	1907	Assert(c.isConst());
149	1907	if (k == Kind::EXPONENTIAL && c.getConst<Rational>().sgn() == 1)
150		{
151	1433	bool success = false;
152	1433	std::uint64_t ds = d;
153	2866	TNode ttrf = getTaylorVariable();
154	2866	TNode tc = c;
155	1605	do
156		{
157	3038	success = true;
158	3038	std::uint64_t n = 2 * ds;
159	6076	std::pair<Node, Node> taylor = getTaylor(k, n);
160		// check that 1-c^{n+1}/(n+1)! > 0
161	6076	Node ru = taylor.second;
162	6076	Node rus = ru.substitute(ttrf, tc);
163	3038	rus = Rewriter::rewrite(rus);
164	3038	Assert(rus.isConst());
165	3038	if (rus.getConst<Rational>() > 1)
166		{
167	1605	success = false;
168	1605	ds = ds + 1;
169		}
170	3038	} while (!success);
171	1433	if (ds > d)
172		{
173	1238	Trace("nl-ext-exp-taylor")
174	619	<< "*** Increase Taylor bound to " << ds << " > " << d << " for ("
175	619	<< k << " " << c << ")" << std::endl;
176		// must use sound upper bound
177	1238	ApproximationBounds pboundss;
178	619	getPolynomialApproximationBounds(k, ds, pboundss);
179	619	pbounds.d_upperPos = pboundss.d_upperPos;
180		}
181	1433	return ds;
182		}
183	474	return d;
184		}
185
186	1037	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	2074	Node c = model.computeAbstractModelValue(tf[0]);
192	1037	Assert(c.isConst());
193	1037	int csign = c.getConst<Rational>().sgn();
194	1037	Kind k = tf.getKind();
195	1037	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	1037	bool isNeg = csign == -1;
209
210	2074	ApproximationBounds pbounds;
211	1037	getPolynomialApproximationBoundForArg(k, c, d, pbounds);
212
213	2074	std::vector<Node> bounds;
214	2074	TNode tfv = getTaylorVariable();
215	2074	TNode tfs = tf[0];
216	3111	for (unsigned d2 = 0; d2 < 2; d2++)
217		{
218		Node pab = (d2 == 0 ? pbounds.d_lower
219	4148	: (isNeg ? pbounds.d_upperNeg : pbounds.d_upperPos));
220	2074	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	4148	Node mtfs = model.computeAbstractModelValue(tfs);
230	2074	pab = pab.substitute(tfv, mtfs);
231	2074	pab = Rewriter::rewrite(pab);
232	2074	Assert(pab.isConst());
233	2074	bounds.push_back(pab);
234		}
235		else
236		{
237		bounds.push_back(Node::null());
238		}
239		}
240	1037	return std::pair<Node, Node>(bounds[0], bounds[1]);
241		}
242
243		} // namespace transcendental
244		} // namespace nl
245		} // namespace arith
246		} // namespace theory
247	22746	} // namespace cvc5