subs.h
1 #ifndef SYMENGINE_SUBS_H
2 #define SYMENGINE_SUBS_H
3 
4 #include <symengine/logic.h>
5 #include <symengine/visitor.h>
6 
7 namespace SymEngine
8 {
9 // xreplace replaces subtrees of a node in the expression tree
10 // with a new subtree
11 RCP<const Basic> xreplace(const RCP<const Basic> &x,
12  const map_basic_basic &subs_dict, bool cache = true);
13 // subs substitutes expressions similar to xreplace, but keeps
14 // the mathematical equivalence for derivatives and subs
15 RCP<const Basic> subs(const RCP<const Basic> &x,
16  const map_basic_basic &subs_dict, bool cache = true);
17 // port of sympy.physics.mechanics.msubs where f'(x) and f(x)
18 // are considered independent
19 RCP<const Basic> msubs(const RCP<const Basic> &x,
20  const map_basic_basic &subs_dict, bool cache = true);
21 // port of sympy's subs where subs inside derivatives are done
22 RCP<const Basic> ssubs(const RCP<const Basic> &x,
23  const map_basic_basic &subs_dict, bool cache = true);
24 
25 class XReplaceVisitor : public BaseVisitor<XReplaceVisitor>
26 {
27 
28 protected:
29  RCP<const Basic> result_;
30  const map_basic_basic &subs_dict_;
31  map_basic_basic visited;
32  bool cache;
33 
34 public:
35  XReplaceVisitor(const map_basic_basic &subs_dict, bool cache = true)
36  : subs_dict_(subs_dict), cache(cache)
37  {
38  if (cache) {
39  visited = subs_dict;
40  }
41  }
42  // TODO : Polynomials, Series, Sets
43  void bvisit(const Basic &x)
44  {
45  result_ = x.rcp_from_this();
46  }
47 
48  void bvisit(const Add &x)
49  {
51  RCP<const Number> coef;
52 
53  auto it = subs_dict_.find(x.get_coef());
54  if (it != subs_dict_.end()) {
55  coef = zero;
56  Add::coef_dict_add_term(outArg(coef), d, one, it->second);
57  } else {
58  coef = x.get_coef();
59  }
60 
61  for (const auto &p : x.get_dict()) {
62  auto it
63  = subs_dict_.find(Add::from_dict(zero, {{p.first, p.second}}));
64  if (it != subs_dict_.end()) {
65  Add::coef_dict_add_term(outArg(coef), d, one, it->second);
66  } else {
67  it = subs_dict_.find(p.second);
68  if (it != subs_dict_.end()) {
69  Add::coef_dict_add_term(outArg(coef), d, one,
70  mul(it->second, apply(p.first)));
71  } else {
72  Add::coef_dict_add_term(outArg(coef), d, p.second,
73  apply(p.first));
74  }
75  }
76  }
77  result_ = Add::from_dict(coef, std::move(d));
78  }
79 
80  void bvisit(const Mul &x)
81  {
82  RCP<const Number> coef = one;
84  for (const auto &p : x.get_dict()) {
85  RCP<const Basic> factor_old;
86  if (eq(*p.second, *one)) {
87  factor_old = p.first;
88  } else {
89  factor_old = make_rcp<Pow>(p.first, p.second);
90  }
91  RCP<const Basic> factor = apply(factor_old);
92  if (factor == factor_old) {
93  // TODO: Check if Mul::dict_add_term is enough
94  Mul::dict_add_term_new(outArg(coef), d, p.second, p.first);
95  } else if (is_a_Number(*factor)) {
96  if (down_cast<const Number &>(*factor).is_zero()) {
97  result_ = factor;
98  return;
99  }
100  imulnum(outArg(coef), rcp_static_cast<const Number>(factor));
101  } else if (is_a<Mul>(*factor)) {
102  RCP<const Mul> tmp = rcp_static_cast<const Mul>(factor);
103  imulnum(outArg(coef), tmp->get_coef());
104  for (const auto &q : tmp->get_dict()) {
105  Mul::dict_add_term_new(outArg(coef), d, q.second, q.first);
106  }
107  } else {
108  RCP<const Basic> exp, t;
109  Mul::as_base_exp(factor, outArg(exp), outArg(t));
110  Mul::dict_add_term_new(outArg(coef), d, exp, t);
111  }
112  }
113 
114  // Replace the coefficient
115  RCP<const Basic> factor = apply(x.get_coef());
116  RCP<const Basic> exp, t;
117  Mul::as_base_exp(factor, outArg(exp), outArg(t));
118  Mul::dict_add_term_new(outArg(coef), d, exp, t);
119 
120  result_ = Mul::from_dict(coef, std::move(d));
121  }
122 
123  void bvisit(const Pow &x)
124  {
125  RCP<const Basic> base_new = apply(x.get_base());
126  RCP<const Basic> exp_new = apply(x.get_exp());
127  if (base_new == x.get_base() and exp_new == x.get_exp()) {
128  result_ = x.rcp_from_this();
129  } else {
130  result_ = pow(base_new, exp_new);
131  }
132  }
133 
134  void bvisit(const OneArgFunction &x)
135  {
136  apply(x.get_arg());
137  if (result_ == x.get_arg()) {
138  result_ = x.rcp_from_this();
139  } else {
140  result_ = x.create(result_);
141  }
142  }
143 
144  template <class T>
145  void bvisit(const TwoArgBasic<T> &x)
146  {
147  RCP<const Basic> a = apply(x.get_arg1());
148  RCP<const Basic> b = apply(x.get_arg2());
149  if (a == x.get_arg1() and b == x.get_arg2())
150  result_ = x.rcp_from_this();
151  else
152  result_ = x.create(a, b);
153  }
154 
155  void bvisit(const MultiArgFunction &x)
156  {
157  vec_basic v = x.get_args();
158  for (auto &elem : v) {
159  elem = apply(elem);
160  }
161  result_ = x.create(v);
162  }
163 
164  void bvisit(const FunctionSymbol &x)
165  {
166  vec_basic v = x.get_args();
167  for (auto &elem : v) {
168  elem = apply(elem);
169  }
170  result_ = x.create(v);
171  }
172 
173  void bvisit(const Contains &x)
174  {
175  RCP<const Basic> a = apply(x.get_expr());
176  auto c = apply(x.get_set());
177  if (not is_a_Set(*c))
178  throw SymEngineException("expected an object of type Set");
179  RCP<const Set> b = rcp_static_cast<const Set>(c);
180  if (a == x.get_expr() and b == x.get_set())
181  result_ = x.rcp_from_this();
182  else
183  result_ = x.create(a, b);
184  }
185 
186  void bvisit(const And &x)
187  {
188  set_boolean v;
189  for (const auto &elem : x.get_container()) {
190  auto a = apply(elem);
191  if (not is_a_Boolean(*a))
192  throw SymEngineException("expected an object of type Boolean");
193  v.insert(rcp_static_cast<const Boolean>(a));
194  }
195  result_ = x.create(v);
196  }
197 
198  void bvisit(const FiniteSet &x)
199  {
200  set_basic v;
201  for (const auto &elem : x.get_container()) {
202  v.insert(apply(elem));
203  }
204  result_ = x.create(v);
205  }
206 
207  void bvisit(const ImageSet &x)
208  {
209  RCP<const Basic> s = apply(x.get_symbol());
210  RCP<const Basic> expr = apply(x.get_expr());
211  auto bs_ = apply(x.get_baseset());
212  if (not is_a_Set(*bs_))
213  throw SymEngineException("expected an object of type Set");
214  RCP<const Set> bs = rcp_static_cast<const Set>(bs_);
215  if (s == x.get_symbol() and expr == x.get_expr()
216  and bs == x.get_baseset()) {
217  result_ = x.rcp_from_this();
218  } else {
219  result_ = x.create(s, expr, bs);
220  }
221  }
222 
223  void bvisit(const Union &x)
224  {
225  set_set v;
226  for (const auto &elem : x.get_container()) {
227  auto a = apply(elem);
228  if (not is_a_Set(*a))
229  throw SymEngineException("expected an object of type Set");
230  v.insert(rcp_static_cast<const Set>(a));
231  }
232  result_ = x.create(v);
233  }
234 
235  void bvisit(const Piecewise &pw)
236  {
237  PiecewiseVec pwv;
238  pwv.reserve(pw.get_vec().size());
239  for (const auto &expr_pred : pw.get_vec()) {
240  const auto expr = apply(*expr_pred.first);
241  const auto pred = apply(*expr_pred.second);
242  pwv.emplace_back(
243  std::make_pair(expr, rcp_static_cast<const Boolean>(pred)));
244  }
245  result_ = piecewise(std::move(pwv));
246  }
247 
248  void bvisit(const Derivative &x)
249  {
250  auto expr = apply(x.get_arg());
251  for (const auto &sym : x.get_symbols()) {
252  auto s = apply(sym);
253  if (not is_a<Symbol>(*s)) {
254  throw SymEngineException("expected an object of type Symbol");
255  }
256  expr = expr->diff(rcp_static_cast<const Symbol>(s));
257  }
258  result_ = expr;
259  }
260 
261  void bvisit(const Subs &x)
262  {
263  auto expr = apply(x.get_arg());
264  map_basic_basic new_subs_dict;
265  for (const auto &sym : x.get_dict()) {
266  insert(new_subs_dict, apply(sym.first), apply(sym.second));
267  }
268  result_ = subs(expr, new_subs_dict);
269  }
270 
271  RCP<const Basic> apply(const Basic &x)
272  {
273  return apply(x.rcp_from_this());
274  }
275 
276  RCP<const Basic> apply(const RCP<const Basic> &x)
277  {
278  if (cache) {
279  auto it = visited.find(x);
280  if (it != visited.end()) {
281  result_ = it->second;
282  } else {
283  x->accept(*this);
284  insert(visited, x, result_);
285  }
286  } else {
287  auto it = subs_dict_.find(x);
288  if (it != subs_dict_.end()) {
289  result_ = it->second;
290  } else {
291  x->accept(*this);
292  }
293  }
294  return result_;
295  }
296 };
297 
299 inline RCP<const Basic> xreplace(const RCP<const Basic> &x,
300  const map_basic_basic &subs_dict, bool cache)
301 {
302  XReplaceVisitor s(subs_dict, cache);
303  return s.apply(x);
304 }
305 
306 class SubsVisitor : public BaseVisitor<SubsVisitor, XReplaceVisitor>
307 {
308 public:
309  using XReplaceVisitor::bvisit;
310 
311  SubsVisitor(const map_basic_basic &subs_dict_, bool cache = true)
312  : BaseVisitor<SubsVisitor, XReplaceVisitor>(subs_dict_, cache)
313  {
314  }
315 
316  void bvisit(const Pow &x)
317  {
318  RCP<const Basic> base_new = apply(x.get_base());
319  RCP<const Basic> exp_new = apply(x.get_exp());
320  if (subs_dict_.size() == 1 and is_a<Pow>(*((*subs_dict_.begin()).first))
321  and not is_a<Add>(
322  *down_cast<const Pow &>(*(*subs_dict_.begin()).first)
323  .get_exp())) {
324  auto &subs_first
325  = down_cast<const Pow &>(*(*subs_dict_.begin()).first);
326  if (eq(*subs_first.get_base(), *base_new)) {
327  auto newexpo = div(exp_new, subs_first.get_exp());
328  if (is_a_Number(*newexpo) or is_a<Constant>(*newexpo)) {
329  result_ = pow((*subs_dict_.begin()).second, newexpo);
330  return;
331  }
332  }
333  }
334  if (base_new == x.get_base() and exp_new == x.get_exp()) {
335  result_ = x.rcp_from_this();
336  } else {
337  result_ = pow(base_new, exp_new);
338  }
339  }
340 
341  void bvisit(const Derivative &x)
342  {
343  RCP<const Symbol> s;
344  map_basic_basic m, n;
345  bool subs;
346 
347  for (const auto &p : subs_dict_) {
348  // If the derivative arg is to be replaced in its entirety, allow
349  // it.
350  if (eq(*x.get_arg(), *p.first)) {
351  RCP<const Basic> t = p.second;
352  for (auto &sym : x.get_symbols()) {
353  if (not is_a<Symbol>(*sym)) {
354  throw SymEngineException("Error, expected a Symbol.");
355  }
356  t = t->diff(rcp_static_cast<const Symbol>(sym));
357  }
358  result_ = t;
359  return;
360  }
361  }
362  for (const auto &p : subs_dict_) {
363  subs = true;
364  if (eq(*x.get_arg()->subs({{p.first, p.second}}), *x.get_arg()))
365  continue;
366 
367  // If p.first and p.second are symbols and arg_ is
368  // independent of p.second, p.first can be replaced
369  if (is_a<Symbol>(*p.first) and is_a<Symbol>(*p.second)
370  and eq(
371  *x.get_arg()->diff(rcp_static_cast<const Symbol>(p.second)),
372  *zero)) {
373  insert(n, p.first, p.second);
374  continue;
375  }
376  for (const auto &d : x.get_symbols()) {
377  if (is_a<Symbol>(*d)) {
378  s = rcp_static_cast<const Symbol>(d);
379  // If p.first or p.second has non zero derivates wrt to s
380  // p.first cannot be replaced
381  if (neq(*zero, *(p.first->diff(s)))
382  || neq(*zero, *(p.second->diff(s)))) {
383  subs = false;
384  break;
385  }
386  } else {
387  result_
388  = make_rcp<const Subs>(x.rcp_from_this(), subs_dict_);
389  return;
390  }
391  }
392  if (subs) {
393  insert(n, p.first, p.second);
394  } else {
395  insert(m, p.first, p.second);
396  }
397  }
398  auto t = x.get_arg()->subs(n);
399  for (auto &p : x.get_symbols()) {
400  auto t2 = p->subs(n);
401  if (not is_a<Symbol>(*t2)) {
402  throw SymEngineException("Error, expected a Symbol.");
403  }
404  t = t->diff(rcp_static_cast<const Symbol>(t2));
405  }
406  if (m.empty()) {
407  result_ = t;
408  } else {
409  result_ = make_rcp<const Subs>(t, m);
410  }
411  }
412 
413  void bvisit(const Subs &x)
414  {
415  map_basic_basic m, n;
416  for (const auto &p : subs_dict_) {
417  bool found = false;
418  for (const auto &s : x.get_dict()) {
419  if (neq(*(s.first->subs({{p.first, p.second}})), *(s.first))) {
420  found = true;
421  break;
422  }
423  }
424  // If p.first is not replaced in arg_ by dict_,
425  // store p.first in n to replace in arg_
426  if (not found) {
427  insert(n, p.first, p.second);
428  }
429  }
430  for (const auto &s : x.get_dict()) {
431  insert(m, s.first, apply(s.second));
432  }
433  RCP<const Basic> presub = x.get_arg()->subs(n);
434  if (is_a<Subs>(*presub)) {
435  for (auto &q : down_cast<const Subs &>(*presub).get_dict()) {
436  insert(m, q.first, q.second);
437  }
438  result_ = down_cast<const Subs &>(*presub).get_arg()->subs(m);
439  } else {
440  result_ = presub->subs(m);
441  }
442  }
443 };
444 
445 class MSubsVisitor : public BaseVisitor<MSubsVisitor, XReplaceVisitor>
446 {
447 public:
448  using XReplaceVisitor::bvisit;
449 
450  MSubsVisitor(const map_basic_basic &d, bool cache = true)
452  {
453  }
454 
455  void bvisit(const Derivative &x)
456  {
457  result_ = x.rcp_from_this();
458  }
459 
460  void bvisit(const Subs &x)
461  {
462  map_basic_basic m = x.get_dict();
463  for (const auto &p : subs_dict_) {
464  m[p.first] = p.second;
465  }
466  result_ = msubs(x.get_arg(), m);
467  }
468 };
469 
470 class SSubsVisitor : public BaseVisitor<SSubsVisitor, SubsVisitor>
471 {
472 public:
473  using XReplaceVisitor::bvisit;
474 
475  SSubsVisitor(const map_basic_basic &d, bool cache = true)
477  {
478  }
479 
480  void bvisit(const Derivative &x)
481  {
482  apply(x.get_arg());
483  auto t = result_;
484  multiset_basic m;
485  for (auto &p : x.get_symbols()) {
486  apply(p);
487  m.insert(result_);
488  }
489  result_ = Derivative::create(t, m);
490  }
491 
492  void bvisit(const Subs &x)
493  {
494  map_basic_basic m = x.get_dict();
495  for (const auto &p : subs_dict_) {
496  m[p.first] = p.second;
497  }
498  result_ = ssubs(x.get_arg(), m);
499  }
500 };
501 
503 inline RCP<const Basic> msubs(const RCP<const Basic> &x,
504  const map_basic_basic &subs_dict, bool cache)
505 {
506  MSubsVisitor s(subs_dict, cache);
507  return s.apply(x);
508 }
509 
511 inline RCP<const Basic> ssubs(const RCP<const Basic> &x,
512  const map_basic_basic &subs_dict, bool cache)
513 {
514  SSubsVisitor s(subs_dict, cache);
515  return s.apply(x);
516 }
517 
518 inline RCP<const Basic> subs(const RCP<const Basic> &x,
519  const map_basic_basic &subs_dict, bool cache)
520 {
521  SubsVisitor b(subs_dict, cache);
522  return b.apply(x);
523 }
524 
525 } // namespace SymEngine
526 
527 #endif // SYMENGINE_SUBS_H
T begin(T... args)
The base class for representing addition in symbolic expressions.
Definition: add.h:27
static RCP< const Basic > from_dict(const RCP< const Number > &coef, umap_basic_num &&d)
Create an appropriate instance from dictionary quickly.
Definition: add.cpp:140
static void coef_dict_add_term(const Ptr< RCP< const Number >> &coef, umap_basic_num &d, const RCP< const Number > &c, const RCP< const Basic > &term)
Updates the numerical coefficient and the dictionary.
Definition: add.cpp:261
const RCP< const Number > & get_coef() const
Definition: add.h:142
The lowest unit of symbolic representation.
Definition: basic.h:95
RCP< const Basic > subs(const map_basic_basic &subs_dict) const
Substitutes 'subs_dict' into 'self'.
Definition: basic.cpp:27
RCP< T > rcp_from_this()
Get RCP<T> pointer to self (it will cast the pointer to T)
virtual RCP< const Basic > create(const vec_basic &x) const
Method to construct classes with canonicalization.
Definition: functions.cpp:1862
static void as_base_exp(const RCP< const Basic > &self, const Ptr< RCP< const Basic >> &exp, const Ptr< RCP< const Basic >> &base)
Convert to a base and exponent form.
Definition: mul.cpp:315
static RCP< const Basic > from_dict(const RCP< const Number > &coef, map_basic_basic &&d)
Create a Mul from a dict.
Definition: mul.cpp:116
virtual RCP< const Basic > create(const vec_basic &v) const =0
Method to construct classes with canonicalization.
virtual vec_basic get_args() const
Returns the list of arguments.
Definition: functions.h:159
virtual RCP< const Basic > create(const RCP< const Basic > &arg) const =0
Method to construct classes with canonicalization.
RCP< const Basic > get_arg() const
Definition: functions.h:36
RCP< const Basic > get_exp() const
Definition: pow.h:42
RCP< const Basic > get_base() const
Definition: pow.h:37
virtual RCP< const Basic > create(const RCP< const Basic > &a, const RCP< const Basic > &b) const =0
Method to construct classes with canonicalization.
RCP< const Basic > get_arg1() const
Definition: functions.h:91
RCP< const Basic > get_arg2() const
Definition: functions.h:96
T emplace_back(T... args)
T end(T... args)
T find(T... args)
T insert(T... args)
T make_pair(T... args)
T move(T... args)
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a_Number(const Basic &b)
Definition: number.h:130
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:426
bool eq(const Basic &a, const Basic &b)
Checks equality for a and b
Definition: basic-inl.h:21
RCP< const Basic > msubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
Subs which treat f(t) and Derivative(f(t), t) as separate variables.
Definition: subs.h:503
RCP< const Basic > exp(const RCP< const Basic > &x)
Returns the natural exponential function E**x = pow(E, x)
Definition: pow.cpp:270
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:347
void insert(T1 &m, const T2 &first, const T3 &second)
Definition: dict.h:83
int factor(const Ptr< RCP< const Integer >> &f, const Integer &n, double B1)
Definition: ntheory.cpp:367
RCP< const Basic > xreplace(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
Mappings in the subs_dict are applied to the expression tree of x
Definition: subs.h:299
bool neq(const Basic &a, const Basic &b)
Checks inequality for a and b
Definition: basic-inl.h:29
RCP< const Basic > ssubs(const RCP< const Basic > &x, const map_basic_basic &subs_dict, bool cache=true)
SymPy compatible subs.
Definition: subs.h:511
T reserve(T... args)
T size(T... args)